戴世坤 王旭龍* 趙東東 劉志偉 張錢江 孫金飛

(①有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室,湖南長(zhǎng)沙 410083； ②中南大學(xué)地球科學(xué)與信息物理學(xué)院,湖南長(zhǎng)沙 410083; ③中國能源建設(shè)集團(tuán)廣東省電力設(shè)計(jì)研究院有限公司，廣東廣州 510663；④桂林理工大學(xué)地球科學(xué)學(xué)院，廣西桂林 541004； ⑤東方地球物理公司大慶物探一公司，黑龍江大慶 163357)

0 引言

任意復(fù)雜條件下的重力異常的高效、高精度數(shù)值模擬，在重力異常反演成像及人機(jī)交互解釋、建模中至關(guān)重要。目前重力異常的正演方法按照求解域的不同一般分為空間域和波數(shù)域兩類。

在空間域正演中，Hubbert[1]給出了二度體重力異常線積分公式； Talwani等[2]推導(dǎo)了截面為多邊形的常密度二度體重力異常解析式，并指出對(duì)于任意二度體都可以用該模型逼近； Rao[3]推導(dǎo)了截面為矩形、密度隨深度呈二次變化的重力異常解析表達(dá)式。后來關(guān)于截面為多邊形、物性連續(xù)變化的二度體重力異常解析式的計(jì)算方法取得了較大進(jìn)展[4-8]。李明等[9]采用級(jí)數(shù)展開計(jì)算了重力垂直一次導(dǎo)數(shù)；徐世浙[10]在二維重力場(chǎng)垂直分量及梯度張量的正演計(jì)算中引入了有限單元法；姜效典等[11]利用樣條插值函數(shù)計(jì)算了變密度地質(zhì)體重力異常；金剛燮等[12]利用等參數(shù)有限元的高斯數(shù)值積分實(shí)現(xiàn)了復(fù)雜形體重力異常正演； Reeder等[13]提出一種采用卷積算子提高二度體重力異常計(jì)算效率的有限元方法； Ren等[14]提出了一種基于非結(jié)構(gòu)化網(wǎng)格的重磁正演快速多極子算法，很好地解決了復(fù)雜邊界和幾何形狀情況下重磁異常的正演問題；肖寶澤等[15]采用均勻多邊形截面二度體重力異常公式實(shí)現(xiàn)了復(fù)雜模型的重力異常數(shù)值模擬。

在波數(shù)域正演中，Sharma等[16]給出了任意斷層面傾角的斷層重力場(chǎng)的波數(shù)域表達(dá)式； Rao等[17]推導(dǎo)一個(gè)截面為等腰三角形的二度體模型的波數(shù)域表達(dá)式； Bhattacharyya等[18]推導(dǎo)了任意二度體的重力異常場(chǎng)波數(shù)域表達(dá)式； 在此基礎(chǔ)上，Pedersen[19]首次推導(dǎo)了以三角形組合作為其截面形狀的二度體重力異常波數(shù)域表達(dá)式，并分析了在波數(shù)域如何以簡(jiǎn)單方式建模，從而解決波數(shù)域反演中的一些問題； 吳宣志[20]將均質(zhì)物性的任意二度體模型推廣至密度隨深度呈線性或指數(shù)變化的模型； 柴玉璞[21]運(yùn)用偏移抽樣方法提高了反傅里葉變換數(shù)值精度； Tontini等[22]實(shí)現(xiàn)了重力異常波數(shù)域正演，并研究了快速傅里葉變換(FFT)擴(kuò)邊法與誤差的關(guān)系； Wu等[23]利用Gauss-FFT提高了反傅里葉變換的數(shù)值精度，降低了FFT引起的強(qiáng)制周期化、邊界震蕩效應(yīng)等影響；商宇航等[24]推導(dǎo)了雙曲線密度模型的波數(shù)域重力異常表達(dá)式。

空間域方法數(shù)值模擬精度高，網(wǎng)格剖分靈活，但計(jì)算大量位置點(diǎn)的異常場(chǎng)數(shù)據(jù)時(shí)，速度較慢。相對(duì)空間域方法，波數(shù)域方法速度較快，但精度相對(duì)較低，且垂向網(wǎng)格等間隔采樣不利于模擬復(fù)雜地形或復(fù)雜模型。為了兼顧數(shù)值模擬的精度與效率，充分利用空間域方法和波數(shù)域方法的優(yōu)勢(shì)，本文針對(duì)任意形狀、任意密度分布的二度體重力異常計(jì)算效率低的問題，提出一種空間—波數(shù)混合域二度體重力異常正演方法。該方法對(duì)空間域引力位滿足的二維偏微分方程沿水平方向進(jìn)行一維傅里葉變換，把空間域引力位滿足的二維偏微分方程轉(zhuǎn)化為不同波數(shù)相互獨(dú)立的一維常微分方程，將一個(gè)復(fù)雜問題分解為多個(gè)小問題。該方法實(shí)現(xiàn)了不同波數(shù)之間常微分方程相互獨(dú)立，具有高度并行性。該方法在垂向上為空間域，便于淺層網(wǎng)格剖分適當(dāng)加密，深層網(wǎng)格剖分適當(dāng)稀疏，有利于適應(yīng)復(fù)雜地形的模擬，兼顧了計(jì)算精度與計(jì)算效率；采用有限單元法求解不同波數(shù)滿足的常微分方程，并充分利用追趕法求解定帶寬線性方程組的高效性，進(jìn)一步提高了計(jì)算效率。

1 理論方法

1.1 控制方程

引力位滿足二維泊松方程[25]

(1)

式中：U為引力位；G為萬有引力常數(shù)； Δρ為異常體剩余密度； (x,z)表示空間任意一點(diǎn)坐標(biāo)。

(2)

特別地，當(dāng)k=0時(shí)，空間—波數(shù)混合域引力位滿足常微分方程

(3)

1.2 邊界條件

在無源區(qū)域，式(2)的通解為

(4)

式中A和B為任意常數(shù)。在笛卡爾坐標(biāo)系中，令坐標(biāo)軸z向下為正，模型的上邊界z=zmin，下邊界z=zmax，如圖1所示。則上、下邊界條件為

(5)

圖1 二度體模型上、下邊界示意圖

1.3 重力場(chǎng)和重力梯度張量的計(jì)算

在無源區(qū)域，可以通過頻率域求導(dǎo)得到重力場(chǎng)和重力梯度張量。空間域重力場(chǎng)g、重力梯度張量T與引力位分別為

(6)

(7)

由式(6)可得空間—波數(shù)混合域引力位與重力場(chǎng)的關(guān)系式

(8)

由式(7)可得空間—波數(shù)混合域引力位與重力梯度張量的關(guān)系式

(9)

2 常微分方程求解

式(2)為引力位在空間—波數(shù)混合域滿足的常微分方程，本文采用基于二次插值的一維有限單元法對(duì)其進(jìn)行數(shù)值模擬。該方法一方面能實(shí)現(xiàn)復(fù)雜模型和復(fù)雜地形的高精度模擬，另一方面在一維條件下利用追趕法可實(shí)現(xiàn)對(duì)角線性方程組(五對(duì)角)的高效、高精度求解。

聯(lián)立式(2)與式(5)可得空間—波數(shù)混合域引力位的邊值問題

(10)

基于變分原理[25]構(gòu)造泛函，可得到與邊值問題(式(10))等價(jià)的變分問題

(11)

(12)

(13)

由于δu≠0，故

Ku=P

(14)

對(duì)該五對(duì)角線性方程組(式(14))的求解參見文獻(xiàn)[26]。通過重力場(chǎng)、重力梯度張量與引力位之間的關(guān)系式(式(6)和式(7))，得到空間—波數(shù)混合域的重力場(chǎng)和重力梯度張量，采用Gauss-FFT對(duì)空間波數(shù)混合域引力位、場(chǎng)及梯度張量進(jìn)行反變換，可得到空間域的引力位、重力場(chǎng)和重力梯度張量。

3 數(shù)值試驗(yàn)

為驗(yàn)證本文算法的正確性及可靠性，分別設(shè)計(jì)了截面為矩形的常密度和變密度二度體模型及任意密度分布的Teplow 模型[13]。最后，對(duì)本文算法的計(jì)算效率隨網(wǎng)格剖分大小的變化規(guī)律進(jìn)行了統(tǒng)計(jì)分析。以下算法均利用Fortran95語言編程串行計(jì)算實(shí)現(xiàn)，筆記本電腦配置為：CPU-Inter Core i4-4790，主頻為4.0GHz，內(nèi)存為32GB。

3.1 常密度二度體模型

設(shè)計(jì)如圖2所示的截面為矩形的常密度二度體模型。研究區(qū)域?yàn)椋簒方向[-500m，500m]，z方向[0，500m]。網(wǎng)格數(shù)為200×100，水平和垂向采樣間隔均為5m。異常體范圍沿x方向?yàn)閇-100m，100m]，z方向?yàn)閇200m，300m]，異常體剩余密度為100kg/m3。重力異常及梯度張量的解析解[27]和數(shù)值解及與地面位置各觀測(cè)點(diǎn)的相對(duì)誤差分別如圖3、圖4所示。

圖2 二度體常密度模型示意圖

圖3 常密度模型重力場(chǎng)解析解和數(shù)值解以及相對(duì)誤差(a)重力異常x分量解析解; (b)重力異常x分量數(shù)值解; (c)地表處相對(duì)誤差； (d)重力異常解析解; (e)重力異常數(shù)值解; (f)地表處相對(duì)誤差

圖4 常密度模型重力梯度張量解析解與數(shù)值解以及地面位置的相對(duì)誤差(a)重力異常垂向梯度解析解; (b)重力異常垂向梯度數(shù)值解; (c)圖a與圖b相對(duì)誤差； (d)重力異常水平梯度解析解; (e)重力異常水平梯度數(shù)值解; (f)圖d與圖e相對(duì)誤差

從圖3可以看出，重力場(chǎng)x分量的數(shù)值解與解析解吻合度較高，地面位置的最大相對(duì)誤差小于1.0×10-4。從圖4可以看出，重力梯度張量的數(shù)值解與解析解吻合較好，地面各觀測(cè)點(diǎn)的重力梯度z分量在零值點(diǎn)附近相對(duì)誤差最大，約為2.0×10-3，這是零值點(diǎn)附近數(shù)值計(jì)算誤差大引起的，在其他位置相對(duì)誤差均小于1.0×10-3。綜合圖3和圖4可以看出，重力場(chǎng)和重力梯度張量的相對(duì)誤差較小，驗(yàn)證了本文算法的正確性和高精度。

3.2 變密度二度體模型

設(shè)計(jì)一個(gè)截面為矩形的變密度二度體模型，研究區(qū)域?yàn)椋簒方向[-500m，500m]，z方向[0,500m]。網(wǎng)格數(shù)為200×100，水平和垂向采樣間隔均為5m。異常體沿x方向范圍為[-100m，100m]，z方向?yàn)閇150m，300m]。異常體的密度隨深度的關(guān)系為

ρ(z)=1.54+0.24z-0.035z2150≤z≤300

當(dāng)z=0時(shí)(即地面)，重力異常和梯度張量解析解[27]及相對(duì)誤差如圖5所示?？梢钥闯觯瑪?shù)值解與解析解形態(tài)十分吻合，重力異常和重力梯度張量的相對(duì)誤差均小于2.0×10-4，與常密度模型計(jì)算結(jié)果相比，變密度模型的計(jì)算結(jié)果精度更高。說明本文算法更適用于變密度二度體模型的重力異常數(shù)值模擬。

圖5 變密度模型重力異常及重力梯度張量解析解與數(shù)值解以及z=0處的相對(duì)誤差(a)重力異常解析解和數(shù)值解; (b)z=0處重力異常相對(duì)誤差; (c)重力異常 水平梯度解析解和數(shù)值解；(d)z=0處重力異常水平梯度相對(duì)誤差

3.3 Teplow密度分布模型

為了進(jìn)一步檢驗(yàn)本文算法對(duì)任意形狀截面、任意密度分布情況下二度體重力異常計(jì)算的適應(yīng)性，引入Reeder等[13]的Teplow密度分布模型(圖6)。研究區(qū)域范圍為：x方向[0，4268.3m]，z方向?yàn)閇0，1829.3m]，剖分網(wǎng)格數(shù)為994×743。

采用本文算法計(jì)算地面位置的重力異常，結(jié)果如圖7所示。由圖可知，傳統(tǒng)空間域累加算法[27]與本文算法的計(jì)算結(jié)果吻合很好，表明本文算法能夠計(jì)算復(fù)雜密度分布情況下的重力異常，且精度較高。

圖6 Teplow密度分布模型[13]

圖7 Teplow密度模型不同方法計(jì)算的重力異常曲線

3.4 計(jì)算效率統(tǒng)計(jì)

計(jì)算效率是評(píng)價(jià)數(shù)值模擬方法好壞的一個(gè)重要指標(biāo)。圖8給出了本文方法的計(jì)算時(shí)間隨網(wǎng)格剖分?jǐn)?shù)目變化的曲線?？梢钥闯?，計(jì)算時(shí)間隨著網(wǎng)格剖分?jǐn)?shù)目的大小呈近似線性增長(zhǎng)。當(dāng)網(wǎng)格剖分?jǐn)?shù)目為500×500時(shí)，即251001個(gè)節(jié)點(diǎn)，采用本文算法計(jì)算整個(gè)剖面耗時(shí)約0.24s，采用傳統(tǒng)空間域累加算法計(jì)算地面501個(gè)節(jié)點(diǎn)需耗時(shí)38.73s[27]，表明本文算法計(jì)算效率更高。

圖8 本文方法計(jì)算時(shí)間隨網(wǎng)格剖分節(jié)點(diǎn)數(shù)的變化曲線

4 結(jié)論

本文提出了一種高效、高精度的空間—波數(shù)混合域二度體重力異常正演方法。設(shè)計(jì)了二度體模型開展數(shù)值試驗(yàn)，通過對(duì)比數(shù)值解與解析解驗(yàn)證了本文方法的正確性和可靠性。通過計(jì)算Teplow密度分布模型重力異常，表明了該算法對(duì)任意密度分布的二度體模型具有很好的適應(yīng)性。數(shù)值算例結(jié)果表明，本文算法具有較高的計(jì)算精度和計(jì)算效率，為計(jì)算任意復(fù)雜形體的重力異常提供了一種新途徑，對(duì)提高二度體重力異常反演成像的效率具有現(xiàn)實(shí)意義。

附錄A

求解文中變分問題(式(11))時(shí)，需將整個(gè)區(qū)域的單元積分分解為各個(gè)單元的積分之和，則式(11)變?yōu)?/p>

(A-1)

其中

(A-2)

令

(A-3)

式中：j、m分別為單元兩端節(jié)點(diǎn)坐標(biāo)；p為單元中點(diǎn)坐標(biāo)。有

(A-4)

式(A-1)中第一項(xiàng)單元積分為

(A-5)

其中

(A-6)

式中l(wèi)為單元長(zhǎng)度。

式(A-1)中第二項(xiàng)單元積分為

(A-7)

其中

(A-8)

式(A-1)中第三項(xiàng)單元積分為

(A-9)

利用形函數(shù)將jaz表示為

(A-10)

其中

szq=(jazj,jazp,jazm)T

(A-11)

則式(A-9)中的單元積分為

(A-12)

其中

(A-13)

式(A-1)中z=zmin、z=zmax分別為第一個(gè)和最后一個(gè)節(jié)點(diǎn)，可將其分別擴(kuò)展成

(A-14)

其中

(A-15)

綜上，可將式(A-1)表示為

F(u)=uTKu-2uTP

(A-16)