李穎梅， 戴世坤， 李 昆， 陳輕蕊， 凌嘉宣

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083)

0 引言

重力場及重力張量場的數(shù)值模擬在地球物理勘探中有著重要作用，但是復(fù)雜二度體的重力數(shù)值模擬不易實現(xiàn)，因此研究高效、高精度的復(fù)雜形體重力場及重力張量場數(shù)值模擬具有重要意義。目前關(guān)于二度體重力場及重力張量場的正演計算方法主要可劃分為解析方法與數(shù)值方法兩大類。

1)解析式法。解析式法是利用數(shù)學(xué)計算方法計算不同模型的解析解，主要是對規(guī)則二度體和密度呈現(xiàn)不同函數(shù)特征分布的重力場與重力張量場解析解的求解。Kwok[1]用共軛復(fù)數(shù)法推導(dǎo)了截面為任意多邊形、密度隨水平和深度呈線性變化的二度半體的源內(nèi)外重力異常；Atienza[2]采用解析與數(shù)值相結(jié)合的方法推導(dǎo)了密度隨水平和垂直方向二次變化的二度體重力異常解析表達式；安玉林等[3-5]介紹了具有連續(xù)物性分布和任意邊界的二維和三維重磁場正反演方法，并推導(dǎo)了截面為規(guī)則形狀的二度體及其組合體復(fù)重磁場表達式；賈真[6]對Won和Bevis提出的新算法的數(shù)值不穩(wěn)定性作了適當修正；王彥國[7]推導(dǎo)了二度體背斜及三度體球冠兩種模型的重力公式；孫明國[8]推導(dǎo)了截面近似圓形和截面為任意形狀的二度體礦脈重力表達式，并給出了重力水平一階導(dǎo)數(shù)、垂向一階導(dǎo)數(shù)和垂向二階導(dǎo)數(shù)的解析公式；蘇和明等[9]重新推導(dǎo)了傾斜臺階重力異常公式，消除了原表達式間斷點的現(xiàn)象；求解析解時，存在奇異點計算和數(shù)值不穩(wěn)定所造成的誤差問題，而且當異常體模型形狀不規(guī)則、密度分布很復(fù)雜時，解析解的求解相當困難。

2)數(shù)值方法。數(shù)值方法主要包括微分法和積分法。微分法利用泊松方程和邊界條件求解整個區(qū)域節(jié)點場值的線性方程組，微分法的關(guān)鍵在于邊界條件的加載。朱自強等[10]在二維重力張量的正演計算中引入有限元法；Reeder等[11]提出一種二度體重力異常計算的有限元方法，采用滑動窗口法來減少邊界效應(yīng)，使用卷積算子來提高計算效率。積分法的優(yōu)勢在于只需對異常體進行剖分，自動滿足邊界條件。主要分為兩大類，一類是通過多個具有解析表達式的規(guī)則二度體逼近，以擬合復(fù)雜地質(zhì)構(gòu)造，或利用求積分的數(shù)值算法求解。例如：Ku[12]提出采用高斯數(shù)值積分法和三次樣條法計算二度體重力異常；徐世浙等[13]提出一種利用高斯數(shù)值積分計算三維變密度體重力異常的正演方法；張嶺[14]利用三角棱柱體組合的重力異常逼近截面為任意形狀的非均勻密度二度體重力場；另一類是譜方法，這種方法結(jié)合傅里葉變換等數(shù)值算法的優(yōu)點，將卷積問題變?yōu)槌朔e問題，大大提高了計算效率。吳樂園[15]提出一種用于重磁位場正演的頻率域高精度正演的Gauss-FFT新方法，Gauss-FFT法避免了零波數(shù)的計算，在同等計算精度的要求下，計算速度比標準FFT擴邊法要快得多。

針對復(fù)雜截面形狀、任意密度分布二度體重力場和重力張量場正演問題，筆者從積分方程出發(fā)，給出一種新的組合矩形模型波數(shù)域重力場和重力張量場表達式，并借助Gauss-FFT方法和三次樣條插值方法，克服了傳統(tǒng)FFT截斷邊界效應(yīng)問題，實現(xiàn)了起伏地形上重力場與重力張量場高效、高精度的二維正演，為實現(xiàn)快速反演與人機交互解釋提供借鑒。

1 正演方法

為模擬任意截面形狀、任意密度分布的復(fù)雜二度體，采用矩形剖分(圖1)。將研究區(qū)域剖分成Nx×Nz個小矩形，(ξi,ζj)表示小矩形中心坐標。每個剖分矩形內(nèi)密度為常值，剖分矩形個數(shù)越多，模型刻畫越精細，越逼近實際的復(fù)雜二度體。上述剖分方式簡單、靈活，適用于復(fù)雜二度體反演和人機交互解釋建模。依據(jù)重力異常的疊加原理，復(fù)雜二度體可以轉(zhuǎn)化為矩形二度體組合模型下重力場及重力張量場的正演問題。

圖1 復(fù)雜二度體正演示意圖Fig.1 Sketch of forward modeling for complex two-dimensional gravity anomaly

研究重點在于快速、準確計算矩形二度體組合模型產(chǎn)生的重力場。在圖1給出的剖分方式下，重力場計算“精確”表達式為式(1)～式(2)。

gz(x,z)=

(1)

gx(x,z)=

(2)

重力張量表達式為式(3)～式(4)。

gzx(x,z)=

(3)

gxx(x,z)=

(4)

式中：G為萬有引力常數(shù)；Nx、Nz分別表示研究區(qū)域x,z方向剖分矩形的個數(shù)；Δx、Δz表示矩形x、z方向采樣間隔，每個矩形水平方向采樣間隔Δx相同；可以根據(jù)密度垂向變化決定Δz大小，這樣可以進行靈活采樣，有助于減小計算量。矩形二度體的中心坐標記為(ξi,ζj)，密度為常值，記為ρ(ξi,ζj)。

本文提出的波數(shù)域重力高效、高精度二維正演方法，重要環(huán)節(jié)之一是對式(1)～式(4)兩邊施加一維傅里葉變換，可以推導(dǎo)得到新的組合模型正演波數(shù)域表達式為式(5)～式(8)。

(5)

(6)

(7)

(8)

式中：k表示x方向的波數(shù)，

(9)

式(9)在形式上類似于一維離散傅里葉變換，在數(shù)值計算實現(xiàn)時，可利用快速傅里葉算法(FFT)實現(xiàn)[16]。

式(5)中，sign(z-ζ)和sign(k)分別表示z和k的符號函數(shù)，表達式如下：

(10)

(11)

式(5)～式(8)為求解的四個空間-波數(shù)混合域重力場與重力張量場的解析表達式，進行一維傅里葉反變換，可得空間域重力場及重力張量場。以重力異常gz為例，反傅里葉變換的表達式為式(12)。

(12)

式(5)～式(8)和式(12)是計算二度體重力場與重力張量場的理論基礎(chǔ)。

在實現(xiàn)一維傅里葉反變換數(shù)值時，采用Gauss-FFT方法，既保證計算效率，又克服傳統(tǒng)FFT算法存在的截斷效應(yīng)及零波數(shù)奇異問題，可提高傅里葉反變換數(shù)值計算精度。

上述新算法能夠?qū)崿F(xiàn)起伏地形上場的快速計算，筆者采用的策略是：根據(jù)地形的分布，確定最高點和最低點之間等間距分布的多個不同高度水平網(wǎng)格上的場值；再采用三次樣條插值，獲得起伏地形上的場值。這種方法的優(yōu)勢在于可利用插值獲得地勢復(fù)雜難以觀測的區(qū)域內(nèi)的場值。

方法流程如下：

1)給定剖分區(qū)域矩形單元個數(shù)以及觀測平面參數(shù)，起伏地形包含在剖分區(qū)域內(nèi)。

2)根據(jù)地質(zhì)體密度分布給每個矩形的剩余密度進行賦值。

3)根據(jù)式(1)～式(4)，利用一維傅里葉變換計算獲得空間-波數(shù)混合域重力場/張量式(5)～式(8)。

4)對空間-波數(shù)混合域重力場/張量反傅里葉變換得到空間域下重力場/張量。

5)采用三次樣條插值，獲得起伏地形上的場值。

2 模型算例

為了驗證本文提出的二維正演快速算法的計算精度，以及對復(fù)雜條件的適應(yīng)性，設(shè)計矩形模型，分別計算了水平線和起伏地形上的場值。將數(shù)值解與解析解進行對比，驗證算法的正確性和高效性。算例中取x軸為東西方向，y軸為北南方向，z軸垂直向下。測試計算機配置為CPU-Intel(R) Core(TM) i7-4770，主頻為3.4 GHz，內(nèi)存為16 GB。傅里葉反變換采用4個高斯點的Gauss-FFT方法實現(xiàn)。

圖2 算例模型和起伏觀測線示意圖Fig.2 Synthetic model and fluctuation of the observation lines

圖3 水平測線重力場數(shù)值解與解析解的曲線圖及相對均方根誤差Fig.3 The numerical solution and the analytical solution of gravity fields when profile is horizontal as well as the relative error between the two

圖4 水平測線重力張量數(shù)值解與解析解的曲線圖及相對均方根誤差Fig.4 The numerical solution and the analytical solution of gradient tensors when profile is horizontal as well as the relative error between the two

模型算例的計算區(qū)域如圖2所示，矩形的x方向和z方向均勻剖分，整個區(qū)域剖分的網(wǎng)格個數(shù)100×100，圍巖密度0 kg/m3。利用本文算法計算水平線上z=-100 m上的重力場與重力張量場的值，并與解析解進行對比。圖3所示為水平線上重力場gx,gz數(shù)值解與解析解的數(shù)值對比及相對誤差；圖4所示為水平線上重力張量場gxx,gzx數(shù)值解與解析解的數(shù)值對比及相對誤差；從圖3與圖4可以看出，重力場與重力張量場的數(shù)值解與解析解具有較高吻合度，相對誤差均較小，由于截斷效應(yīng)，重力場在邊界處相對誤差稍大，最大約為0.004%；在邊界處與零值附近，重力張量數(shù)值計算誤差稍大，最大約為0.4%。由圖3和圖4可知，張量的變化趨勢較場更加劇烈，場的計算精度高于張量。

為了研究起伏地形上場值的計算精度隨著插值曲線條數(shù)的變化，設(shè)計一個正弦函數(shù)變化的起伏地形，在觀測區(qū)域的函數(shù)表達式為式(13)。

(13)

起伏地形上最高點與最低點的z方向坐標為-300 m和-100 m。

引入相對均方根誤差，統(tǒng)計整條起伏觀測線上場值的誤差。該誤差統(tǒng)計方式能突出大異常值所占的比重，計算公式為式(14)。

(14)

其中：rrms表示相對均方根誤差；N表示x方向的節(jié)點數(shù)；gi表示重力異常(重力張量)數(shù)值解；Ti表示重力異常(重力張量)的解析解。

圖5為起伏地形上重力場及重力張量場相對均方根誤差隨水平插值個數(shù)的變化趨勢圖。從圖5中可以看出，起伏地形上重力場與重力張量場的計算精度隨著插值個數(shù)的增加而提高。當用3個不同高度進行插值時，起伏地形上重力場的rrms已小于0.5%，而重力張量場rrms較大，最大約為2.5%。當選取的高度網(wǎng)格節(jié)點個數(shù)大于11時，起伏地形上重力場與重力張量場的rrms幾乎不再變化，計算精度達到最優(yōu)。圖6所示插值個數(shù)11時，起伏地形上重力場與張量的曲線圖。從圖6中可以看出，起伏地形上數(shù)值解與解析解的吻合度很高；選取15個不同高度網(wǎng)格節(jié)點上的場值進行插值時，算法耗時僅為0.020 s，這說明本文提出的插值策略是可行的，同時也證實了算法在二度體重力場及重力張量場正演計算時，具有很高的計算精度和計算效率。

圖5 相對均方根誤差隨插值個數(shù)變化Fig.5 Rrms with a varying number of interpolation dots

圖6 起伏觀測網(wǎng)格節(jié)點重力場與重力張量的數(shù)值解與解析解Fig.6 The numerical solution and the analytical solution of gravitational fields when profile is horizontal as well as the relative error between the two

利用文獻[17]中格爾木—花海子剖面巖石圈二維密度結(jié)構(gòu)模型，擬合密度分布的模型如圖7所示，計算地面z=0 m水平線上重力場和重力梯度如圖8所示。

圖7 格爾木-花海子剖面巖石圈二維密度結(jié)構(gòu)模型Fig.7 2-D lithospheric density structure along the profile from Golmud to Huahaizi

圖8 水平測線重力場與重力張量數(shù)值解Fig.8 The numerical solution of gravity fields and gradient tensors when profile is horizontal

3 結(jié)論

筆者針對任意密度分布情況下、起伏觀測線二度體重力場及重力張量高效、高精度正演問題，提出了一種基于Gauss-FFT和三次樣條插值的空間波數(shù)混合域正演算法。結(jié)論如下：

1)算法采用矩形二度體組合模型的波數(shù)域計算公式，將二維空間域卷積拆分成波數(shù)域多個一維積分累加，大大減少了計算量。同時，不同波數(shù)場的計算相互獨立，采用并行算法可進一步提高計算速度。

2)每個波數(shù)對應(yīng)的空間域一維積分可得到解析解，計算精準；利用Gauss-FFT反傅里葉變換，避免零波數(shù)計算，能克服傳統(tǒng)FFT法的截斷效應(yīng)，進一步提高了計算精度。采用矩形二度體組合，能模擬地下任意復(fù)雜地質(zhì)條件，從而實現(xiàn)重力勘探復(fù)雜條件二度體高效、高精度正演計算。

3)網(wǎng)格剖分100×100時，水平網(wǎng)格節(jié)點重力場相對誤差最大僅為0.004%，重力張量場為0.4%；計算起伏地形15個不同高度水平網(wǎng)格上的場值時僅用0.02 s；網(wǎng)格數(shù)1 000×1 000時，計算時間1.3 s?？梢?，隨著網(wǎng)格數(shù)量的增大，算法耗時近似線性增長，算法具有良好的伸縮性。

4)筆者提出的算法，還可應(yīng)用于三維重磁問題的求解，為重磁數(shù)據(jù)的高效、高精度反演、人機交互正反演解釋提供了新途徑。