美文网首页
番外篇6_Geoist之曲面异常拟合

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

作者: 地学小哥 | 来源:发表于2020-05-05 18:08 被阅读0次

内容摘要:在重力资料的处理分析解释中,地形、盆地、莫霍面等通常都具有三维曲面特性。前文我们讨论过正六面体这种离散单元形式,对曲面网格数据进行六面体近似,可以一定程度上解决曲面类的异常正反演问题,包括:完全布格重力异常改正和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算法),在后续的文中我们会专门介绍,感兴趣的小盆友们欢迎关注哦!

相关文章

网友评论

      本文标题:番外篇6_Geoist之曲面异常拟合

      本文链接:https://www.haomeiwen.com/subject/omvnmhtx.html