内容摘要：在了解基本的位场理论，异常计算方法，插值、平滑等数据处理方后，我们今天进一步介绍重磁异常的解析延拓，化极，求导，欧拉反演等内容。由于本节方法多需要在频率域内进行，需要FFT运算，因此对数据格式有一定要求。
在开始介绍今天内容之前，我们看一下重力位场数据处理的基本流程，如图1所示。图中将一个重力项目分成了项目设计、数据采集、数据处理和解释四个阶段，简述了每一个阶段需要考虑的问题和要点，理解图1内容，是更好应用每项技术的基础。
本专辑介绍的内容，都偏重怎么实现一个具体的技术内容，如果将每个技术看成招式，那么图1就是心法，两者结合才能更好地运行这些工具。
图1 重力位场数据处理和解释流程
首先，是加载数据，在上一节中，我们将网格化好的文件都保存成了surfer的grd格式，代码如下：
from geoist.others.utils import grid2Grid,map2DGrid
from geoist.pfm.grdio import grddata
rBGAg2d = grid2Grid(grid2.resBGA.easting, grid2.resBGA.northing, grid2.resBGA.values)
g1out = grddata()
g1out.cols = rBGAg2d.getGeoDict().nx
g1out.rows = rBGAg2d.getGeoDict().ny
g1out.xmin = rBGAg2d.getGeoDict().xmin
g1out.xmax = rBGAg2d.getGeoDict().xmax
g1out.ymin = rBGAg2d.getGeoDict().ymin
g1out.ymax = rBGAg2d.getGeoDict().ymax
g1out.data0 = grid2.resBGA.values
g1out.export_surfer(Path(DATA_PATH, 'ynyx_bgar.grd'), False, 'ascii')
然后，加载到网格文件，代码如下：
import matplotlib.pyplot as plt
from pathlib import Path
from geoist import DATA_PATH
from geoist.others.gdal import GDALGrid
from geoist.others.utils import map2DGrid
# 1.读取数据
filename1 = Path(DATA_PATH, 'ynyx_elev.grd')
gd1, ff = GDALGrid.getFileGeoDict(filename1)
grid1 = GDALGrid.load(filename1,gd1)
filename2 = Path(DATA_PATH, 'ynyx_fga.grd')
gd2, ff = GDALGrid.getFileGeoDict(filename2)
grid2 = GDALGrid.load(filename2,gd2)
filename3 = Path(DATA_PATH, 'ynyx_bgas.grd')
gd3, ff = GDALGrid.getFileGeoDict(filename3)
grid3 = GDALGrid.load(filename3,gd3)
filename4 = Path(DATA_PATH, 'ynyx_bgar.grd')
gd4, ff = GDALGrid.getFileGeoDict(filename4)
grid4 = GDALGrid.load(filename4,gd4)
结果如图2所示，由于区域比较小，坐标轴单位为km。
图2 加载的Surfer网格文件
1、位场的解析延拓
位场的延拓需要使用pftrans模块下的upcontinue函数，其中有两个参数，一个是shape，另一个是height需要设定。height单位为m，向上为正。代码如下：
from geoist.pfm import pftrans
from geoist.vis import giplt
from geoist.others.utils import Grid2Xyz
x, y ,bgas = Grid2Xyz(grid3)
shape = (grid3.getGeoDict().nx, grid3.getGeoDict().ny)
height = 1000 # How much higher to go
bgas_contf = pftrans.upcontinue(x, y, bgas, shape, height)
args = dict(shape=shape, levels=20, cmap=plt.cm.RdBu_r)
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes = axes.ravel()
plt.sca(axes[0])
plt.title("Original")
plt.axis('scaled')
giplt.contourf(x, y, bgas, **args)
plt.colorbar(pad=0).set_label('mGal')
giplt.m2km()
plt.sca(axes[1])
plt.title('Upward continuation 1000m')
plt.axis('scaled')
giplt.contourf(x, y, bgas_contf , **args)
plt.colorbar(pad=0).set_label('mGal')
giplt.m2km()
fig.tight_layout()
这些需要注意的是，upcontinue函数需要输入的位场数据格式为x,y,z形式，这三个都必须是1-d array，与读入的grid格式在排列上有一定区别，我们先用Grid2Xyz函数转换了一下。向上延拓1000m后的布格重力异常如图3所示。
图3 重力异常向上延拓1000m
向上延拓1000m后，异常是不是变光滑了！其实延拓也起到了一种滤波效果。
2、空间导数计算
空间导数的计算与延拓类似，都在pftrans模块下，函数名称为：derivx，derivy，derivz。分别对应三个方向。参数多了一个order，可以设置求取不同阶数。另外，还有一个函数tga，这个可以用于求位场的解析信号，
我们将刚才的数据分别求导数和解析信号，代码形式如下：
bgas_dx = pftrans.derivx(x, y, bgas_contf, shape)
bgas_dy = pftrans.derivy(x, y, bgas_contf, shape)
bgas_dz = pftrans.derivz(x, y, bgas_contf, shape)
bgas_tga = pftrans.tga(x, y, bgas_contf, shape)
args = dict(shape=shape, levels=20, cmap=plt.cm.RdBu_r)
fig, axes = plt.subplots(2, 2, figsize=(10, 10))
plt.sca(axes[0, 0])
plt.title("Derivative of BGA in X")
plt.axis('scaled')
giplt.contourf(x, y, bgas_dx, **args)
plt.colorbar(pad=0).set_label('mGal/m')
giplt.m2km()
plt.sca(axes[0, 1])
plt.title("Derivative of BGA in Y")
plt.axis('scaled')
giplt.contourf(x, y, bgas_dy, **args)
plt.colorbar(pad=0).set_label('mGal/m')
giplt.m2km()
plt.sca(axes[1, 0])
plt.title("Derivative of BGA in Z")
plt.axis('scaled')
giplt.contourf(x, y, bgas_dz, **args)
plt.colorbar(pad=0).set_label('mGal/m')
giplt.m2km()
plt.sca(axes[1, 1])
plt.title("Total gradient amplitude of BGA")
plt.axis('scaled')
giplt.contourf(x, y, bgas_tga, **args)
plt.colorbar(pad=0).set_label('mGal/m')
giplt.m2km()
效果如图4所示：
图4 空间导数计算
小结：今天讲的内容结合上次讲的算是本套教程的重磁资料处理专业部分开始。本想把重磁方法分开来讲，但是由于时间关系，没有来得及准备。后续我们会继续打磨，更新稿件。本套教程共设计了30回，期望在五一假期之前写完，而下一节课将开始三维位场反演方面的内容讲解，这些内容是Geoist软件包的精华，因为，很多技术都是我们团队的原创成果，甚至一些还是没有公开发表过的，如果您也对这方面工作感兴趣，欢迎继续关注吧！
网友评论