番外篇6_Geoist之曲面异常拟合

内容摘要：在重力资料的处理分析解释中，地形、盆地、莫霍面等通常都具有三维曲面特性。前文我们讨论过正六面体这种离散单元形式，对曲面网格数据进行六面体近似，可以一定程度上解决曲面类的异常正反演问题，包括：完全布格重力异常改正和Moho面的反演等。

1、曲面的六面体近似

一般曲面的近似采用三角网格或四面体网格较多，如常用表示的数字高程模型的不规则三角网(简称TIN,即Triangulated Irregular Network)算法。

但是，不规则四面体的重磁位场异常计算比较复杂，不利于离散化后的模型快速计算异常。而正六面体单元的位场异常计算简单，在一定精度水平上，更适合逼近曲面异常体。如图1所示，采用六面体单元离散化的带地形起伏的场源模型。

图1 六面体近似效果示意图

今天，我们就Geoist软件包中，提供的正六面体曲面拟合功能进行讲解。

2、认识PrismRelief类

在Geosit中的inversion模块mesh.py里的PrismRelief类。在PrismRelief类的实例化函数里面，有三个参数分别为：reference, shape, nodes，其中，reference为参考平面（起算面高度）；shape网格点数元组；nodes规则网格上的坐标和高度，元组类型。

一个典型的实现代码如下：

from geoist import gridder
from geoist.pfm import giutils
from geoist.inversion import mesh

area = (-150, 150, -300, 300)
shape = (100, 50)
x, y = gridder.regular(area, shape)
height = (-80 * giutils.gaussian2d(x, y, 100, 200, x0=-50, y0=-100, angle=-60) +
          100 * giutils.gaussian2d(x, y, 50, 100, x0=80, y0=170))

nodes = (x, y, -1 * height)  # -1 is to convert height to z coordinate
reference = 0  # z coordinate of the reference surface
relief = mesh.PrismRelief(reference, shape, nodes)
relief.addprop('density', (2670 for i in range(relief.size)))

上述代码中，PrismRelief对象的实例relief具有遍历特性，for函数可以实现单元的遍历，可以对其进行属性参数设置和提取，进而计算观测点或网格位置的重磁异常。

prop = 'density'
for prism in relief:
    if prism is None or (prop is not None and prop not in prism.props):
        continue
    x1, x2, y1, y2, z1, z2 = prism.get_bounds()

一句话总结：基于正六面体单元可以实现地形、盆地基底、Moho面等起伏曲面的异常模拟。在布格重力异常的地形改正，Moho面深度反演等方面，都具有重要的应用前景。另外，在重磁位场正反演计算中，还可以采用频率域算法来实现曲面类异常计算（Parker算法），在后续的文中我们会专门介绍，感兴趣的小盆友们欢迎关注哦！

最后编辑于：2020.05.05 18:11:51