內(nèi)容摘要：在重力資料的處理分析解釋中攻人，地形距芬、盆地、莫霍面等通常都具有三維曲面特性穗酥。前文我們討論過正六面體這種離散單元形式护赊，對曲面網(wǎng)格數(shù)據(jù)進行六面體近似，可以一定程度上解決曲面類的異常正反演問題砾跃，包括：完全布格重力異常改正和Moho面的反演等骏啰。
1、曲面的六面體近似
一般曲面的近似采用三角網(wǎng)格或四面體網(wǎng)格較多抽高，如常用表示的數(shù)字高程模型的不規(guī)則三角網(wǎng)(簡稱TIN,即Triangulated Irregular Network)算法判耕。
但是，不規(guī)則四面體的重磁位場異常計算比較復(fù)雜翘骂，不利于離散化后的模型快速計算異常壁熄。而正六面體單元的位場異常計算簡單，在一定精度水平上雏胃，更適合逼近曲面異常體请毛。如圖1所示，采用六面體單元離散化的帶地形起伏的場源模型瞭亮。
今天，我們就Geoist軟件包中固棚，提供的正六面體曲面擬合功能進行講解统翩。
2仙蚜、認識PrismRelief類
在Geosit中的inversion模塊mesh.py里的PrismRelief類。在PrismRelief類的實例化函數(shù)里面厂汗，有三個參數(shù)分別為：reference, shape, nodes委粉，其中，reference為參考平面（起算面高度）娶桦；shape網(wǎng)格點數(shù)元組贾节；nodes規(guī)則網(wǎng)格上的坐標和高度，元組類型衷畦。
一個典型的實現(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對象的實例relief具有遍歷特性，for函數(shù)可以實現(xiàn)單元的遍歷祈争，可以對其進行屬性參數(shù)設(shè)置和提取斤程，進而計算觀測點或網(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()
一句話總結(jié)：基于正六面體單元可以實現(xiàn)地形菩混、盆地基底忿墅、Moho面等起伏曲面的異常模擬。在布格重力異常的地形改正沮峡，Moho面深度反演等方面疚脐，都具有重要的應(yīng)用前景。另外邢疙，在重磁位場正反演計算中亮曹，還可以采用頻率域算法來實現(xiàn)曲面類異常計算（Parker算法），在后續(xù)的文中我們會專門介紹秘症，感興趣的小盆友們歡迎關(guān)注哦照卦！
