內(nèi)容摘要：在重力資料的處理分析解釋中，地形、盆地、莫霍面等通常都具有三維曲面特性。前文我們討論過(guò)正六面體這種離散單元形式，對(duì)曲面網(wǎng)格數(shù)據(jù)進(jìn)行六面體近似，可以一定程度上解決曲面類(lèi)的異常正反演問(wèn)題，包括：完全布格重力異常改正和Moho面的反演等。
1、曲面的六面體近似
一般曲面的近似采用三角網(wǎng)格或四面體網(wǎng)格較多，如常用表示的數(shù)字高程模型的不規(guī)則三角網(wǎng)(簡(jiǎn)稱(chēng)TIN,即Triangulated Irregular Network)算法。
但是，不規(guī)則四面體的重磁位場(chǎng)異常計(jì)算比較復(fù)雜，不利于離散化后的模型快速計(jì)算異常。而正六面體單元的位場(chǎng)異常計(jì)算簡(jiǎn)單，在一定精度水平上，更適合逼近曲面異常體。如圖1所示，采用六面體單元離散化的帶地形起伏的場(chǎng)源模型。
今天，我們就Geoist軟件包中，提供的正六面體曲面擬合功能進(jìn)行講解。
2、認(rèn)識(shí)PrismRelief類(lèi)
在Geosit中的inversion模塊mesh.py里的PrismRelief類(lèi)。在PrismRelief類(lèi)的實(shí)例化函數(shù)里面，有三個(gè)參數(shù)分別為：reference, shape, nodes，其中，reference為參考平面（起算面高度）；shape網(wǎng)格點(diǎn)數(shù)元組；nodes規(guī)則網(wǎng)格上的坐標(biāo)和高度，元組類(lèi)型。
一個(gè)典型的實(shí)現(xiàn)代碼如下：
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對(duì)象的實(shí)例relief具有遍歷特性，for函數(shù)可以實(shí)現(xiàn)單元的遍歷，可以對(duì)其進(jìn)行屬性參數(shù)設(shè)置和提取，進(jìn)而計(jì)算觀(guān)測(cè)點(diǎn)或網(wǎng)格位置的重磁異常。
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()
一句話(huà)總結(jié)：基于正六面體單元可以實(shí)現(xiàn)地形、盆地基底、Moho面等起伏曲面的異常模擬。在布格重力異常的地形改正，Moho面深度反演等方面，都具有重要的應(yīng)用前景。另外，在重磁位場(chǎng)正反演計(jì)算中，還可以采用頻率域算法來(lái)實(shí)現(xiàn)曲面類(lèi)異常計(jì)算（Parker算法），在后續(xù)的文中我們會(huì)專(zhuān)門(mén)介紹，感興趣的小盆友們歡迎關(guān)注哦！
