內(nèi)容摘要：在番外篇1的時(shí)候，我們?cè)?jīng)談到了Tesseroid的單元定義。但是沒有提到怎么生成該單元的網(wǎng)格。該類型網(wǎng)格單元可以直接采用經(jīng)緯度坐標(biāo)來剖分，對(duì)于解決大尺度的問題更方便，因此我們說說網(wǎng)格生成批處理方法，TesseroidMesh類的使用。
1、Tesseroid網(wǎng)格
在Geosit中的inversion模塊mesh.py里的TesseroidMesh類，實(shí)現(xiàn)了Tesseroid單元的自動(dòng)網(wǎng)格剖分。在球坐標(biāo)系下，Tesseroid形式的單元，直接經(jīng)緯度剖分網(wǎng)格即可，不用變換為直角坐標(biāo)。這對(duì)大尺度模型解釋時(shí)候常常是有用的，圖1是該網(wǎng)格的示意圖。
2、網(wǎng)格數(shù)據(jù)生成
在Geoist中，inversion模塊下geometry.py定義了各種幾何單元類型，mesh.py是以不同單元剖分網(wǎng)格的批處理類實(shí)現(xiàn)。已經(jīng)實(shí)現(xiàn)的幾種類型見下表。
| 單元 | 網(wǎng)格 | 類型 |
|---|---|---|
| Prism | PrismMesh | 六面體 |
| Tesseroid | TesseroidMesh | 球面六面體 |
| Prism | PrismRelief | 起伏界面 |
| Square | SquareMesh | 2D正方形 |
| Sphere | PointGrid | 點(diǎn)源-球體 |
實(shí)例化后的網(wǎng)格對(duì)象，都可以通過addprop的方法添加屬性。同時(shí)，如果需要得到這個(gè)網(wǎng)格對(duì)象中每個(gè)單元的幾何屬性，也可以使用for循環(huán)通過遍歷實(shí)現(xiàn)（for cell in mesh: ...）。
下面是定義一個(gè)TesseroidMesh，并添加每個(gè)單元屬性的過程。具體代碼如下：
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)
一句話總結(jié)：今天以TesseroidMesh入手，簡單總結(jié)了一下場源剖分常用的單元定義和離散化（網(wǎng)格）剖分方法。地球物理以場求源的過程中，經(jīng)常需要將場源進(jìn)行離散化處理，進(jìn)一步計(jì)算格林函數(shù)矩陣，然后再求解，以上這些內(nèi)容都非?；A(chǔ)，萬丈高樓平地起，明白原理和邏輯后，很快你就可以寫出非常優(yōu)秀的代碼啦。
