内容摘要：在番外篇1的时候，我们曾经谈到了Tesseroid的单元定义。但是没有提到怎么生成该单元的网格。该类型网格单元可以直接采用经纬度坐标来剖分，对于解决大尺度的问题更方便，因此我们说说网格生成批处理方法，TesseroidMesh类的使用。
1、Tesseroid网格
在Geosit中的inversion模块mesh.py里的TesseroidMesh类，实现了Tesseroid单元的自动网格剖分。在球坐标系下，Tesseroid形式的单元，直接经纬度剖分网格即可，不用变换为直角坐标。这对大尺度模型解释时候常常是有用的，图1是该网格的示意图。
图1 Tesseroid网格示意图
2、网格数据生成
在Geoist中，inversion模块下geometry.py定义了各种几何单元类型，mesh.py是以不同单元剖分网格的批处理类实现。已经实现的几种类型见下表。
| 单元 | 网格 | 类型 |
|---|---|---|
| Prism | PrismMesh | 六面体 |
| Tesseroid | TesseroidMesh | 球面六面体 |
| Prism | PrismRelief | 起伏界面 |
| Square | SquareMesh | 2D正方形 |
| Sphere | PointGrid | 点源-球体 |
实例化后的网格对象，都可以通过addprop的方法添加属性。同时，如果需要得到这个网格对象中每个单元的几何属性，也可以使用for循环通过遍历实现（for cell in mesh: ...）。
下面是定义一个TesseroidMesh，并添加每个单元属性的过程。具体代码如下：
from geoist import gridder
from geoist.pfm import giutils
from geoist.inversion import mesh
w, e = -2, 2
s, n = -2, 2
bounds = (w, e, s, n, 500000, 0)
x, y = gridder.regular((w, e, s, n), (50, 50))
height = (250000 +
-100000 * giutils.gaussian2d(x, y, 1, 5, x0=-1, y0=-1, angle=-60) +
250000 * giutils.gaussian2d(x, y, 1, 1, x0=0.8, y0=1.7))
mesh = mesh.TesseroidMesh(bounds, (20, 50, 50))
mesh.addprop('density', (2670 for i in range(mesh.size)))
mesh.carvetopo(x, y, height)
一句话总结：今天以TesseroidMesh入手，简单总结了一下场源剖分常用的单元定义和离散化（网格）剖分方法。地球物理以场求源的过程中，经常需要将场源进行离散化处理，进一步计算格林函数矩阵，然后再求解，以上这些内容都非常基础，万丈高楼平地起，明白原理和逻辑后，很快你就可以写出非常优秀的代码啦。
网友评论