復(fù)雜條件二維重力場及重力張量場空間波數(shù)域正演方法

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

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

0 引言

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

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

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

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

1 正演方法

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

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

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

gz(x,z)=

(1)

gx(x,z)=

(2)

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

gzx(x,z)=

(3)

gxx(x,z)=

(4)

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

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

(5)

(6)

(7)

(8)

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

(9)

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

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

(10)

(11)

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

(12)

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

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

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

方法流程如下：

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

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

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

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

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

2 模型算例

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

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

圖3 水平測線重力場數(shù)值解與解析解的曲線圖及相對(duì)均方根誤差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ù)值解與解析解的曲線圖及相對(duì)均方根誤差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

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

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

(13)

起伏地形上最高點(diǎn)與最低點(diǎn)的z方向坐標(biāo)為-300 m和-100 m。

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

(14)

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

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

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

圖6 起伏觀測網(wǎng)格節(jié)點(diǎn)重力場與重力張量的數(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

利用文獻(xiàn)[17]中格爾木—花海子剖面巖石圈二維密度結(jié)構(gòu)模型，擬合密度分布的模型如圖7所示，計(jì)算地面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é)論

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

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

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

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

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