趙東東， 趙亞飛， 李 昆， 周印明，王旭龍， 戴世坤

(1. 中國地質(zhì)調(diào)查局 南京地質(zhì)調(diào)查中心，南京 210016；2. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 長沙 410083；3. 有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室, 長沙 410083;4. 西南石油大學(xué) 地球科學(xué)與技術(shù)學(xué)院，成都 610500;5. 湖南科技大學(xué) 環(huán)境與安全工程學(xué)院，湘潭 411201;6. 中國石油天然氣集團(tuán)公司 川慶鉆探工程有限公司 長慶井下技術(shù)作業(yè)公司，西安 710018)

0 引言

磁異常正演數(shù)值模擬，是磁法勘探數(shù)據(jù)定量解釋的基礎(chǔ)。有關(guān)二度體的磁異常數(shù)值模擬研究最早是在空間域開始的，Talwani等[1]首次給出了多邊形截面二度體磁異常解析表達(dá)式；在此基礎(chǔ)上，Talwani[2]在計算機(jī)上實現(xiàn)了近似二維棱柱狀多邊形數(shù)值模擬；Won等[3]根據(jù)泊松關(guān)系導(dǎo)出了多邊形截面二度體磁異常的解析式，并利用fortran77實現(xiàn)了相關(guān)算法，優(yōu)勢在于可以模擬地下任意觀測點的磁場；Shuey等[4]，Rasmussen等[5]及Cady等[6]在前人研究基礎(chǔ)上對原來的算法進(jìn)行校正，使其能模擬實際地質(zhì)體所產(chǎn)生的磁異常；賈真[7]系統(tǒng)地推導(dǎo)了均勻多邊形截面二度體磁異常及磁異常梯度正演計算公式，并重點探討了正演公式中所包含的奇異性；Mehdi等[8]采用自適應(yīng)有限單元法實現(xiàn)了基于拉普拉斯方程的二度體磁異常正演；Adim等[9]在Talwani[2]解析公式的基礎(chǔ)上編寫了一套基于Matlab的二度體磁異常模擬軟件，并提高了數(shù)值精度。在頻率域，Bhattacharyya[10]推導(dǎo)了任意二度體在頻率域的磁場表達(dá)式，并進(jìn)行了反演研究；Pedersen[11]給出了任意二度體、二度半體(以三角形組合為其表面的模型)重磁異常頻率域表達(dá)式，并指出該公式相對空間域解析公式更為簡潔；吳宣志[12]將任意二維物體磁場的波譜解析表達(dá)式推廣至可用于任意變磁化強(qiáng)度模型；柴玉璞[13]將偏移抽樣方法發(fā)展為一整套完整的理論體系,有效提升了傅里葉變換的數(shù)值精度,使得頻率域方法逐漸在處理復(fù)雜模型正演問題中廣泛應(yīng)用；吳樂園[14]在偏移抽樣理論的基礎(chǔ)上，提出利用高斯FFT提高傅里葉變換的數(shù)值精度,降低了由于快速傅里葉變換引起的強(qiáng)制周期化、邊界震蕩效應(yīng)等影響，并將其應(yīng)用于頻率域二度體磁異常正演模擬，取得了不錯的效果。

筆者從磁異常積分表達(dá)式出發(fā)，提出了一種基于積分解的空間-波數(shù)混合域磁異常數(shù)值模擬方法。該方法利用高斯FFT[14]將空間域磁位滿足的二維積分問題轉(zhuǎn)化為不同波數(shù)相互獨立的一維積分問題，將一個二維問題分解為波數(shù)之間相互獨立的一維問題，保留深度方向為空間域，便于淺層網(wǎng)格適當(dāng)加密，深層網(wǎng)格適當(dāng)稀疏，兼顧計算精度、計算效率及模擬復(fù)雜地形的需要，并采用基于二次插值的形函數(shù)實現(xiàn)垂向一維積分，提升積分精度和效率。模型算例中分別設(shè)計了突變介質(zhì)模型和起伏地形模型，首先利用突變介質(zhì)模型驗證了該方法的正確性；其次通過對比分析基于三角函數(shù)變化的起伏地形模型磁異?？焖偎惴?，與常規(guī)積分算法的計算精度與計算效率，驗證了起伏地形條件下基于形函數(shù)積分的磁異常快速算法的有效性和高效性，并研究了地形起伏最低點與最高點之間剖分節(jié)點個數(shù)對磁異常場計算精度的影響規(guī)律；最后設(shè)計非規(guī)則復(fù)雜模型，研究了不同觀測高度的磁異常響應(yīng)特征。

1 方法原理

1.1 空間-波數(shù)混合域磁位

弱磁性體空間域磁位積分公式為式(1)。

(1)

式中：(x,z)和(x′,z′)分別表示觀測點和場源任意一點的坐標(biāo)；U為異常體在空間任意一點產(chǎn)生的磁位；Mx和Mz分別為磁化強(qiáng)度矢量M在x和z方向上的分量；π為圓周率常數(shù)；μ0為真空磁導(dǎo)率，μ0= 4π×10-7H/m；如果場源均勻，磁化強(qiáng)度M為常矢量。

對式(1)沿x方向做一維傅里葉變換，得到空間-波數(shù)混合域磁異常積分公式為式(2)。

(2)

(3)

(4)

式(2)即空間-波數(shù)混合域磁異常積分公式。對式(1)沿x方向進(jìn)行傅里葉變換，把空間域滿足的二維積分問題轉(zhuǎn)化為不同波數(shù)相互獨立的一維積分問題，并行性高；垂向一維積分采用基于二次插值的形函數(shù)積分實現(xiàn)，將z方向離散成M個單元，每個單元積分均可得到解析表達(dá)式，計算精度高；垂向保持空間域，便于淺層網(wǎng)格適當(dāng)加密，深層網(wǎng)格適當(dāng)稀疏，兼顧計算精度、計算效率及模擬復(fù)雜地形的需求。

1.2 磁場和磁梯度張量

根據(jù)空間域磁場、磁梯度張量與磁位的關(guān)系

(5)

式中，B為空間域磁感應(yīng)強(qiáng)度。

(6)

式中，T為空間域磁梯度張量。由式(5)得空間波數(shù)混合域磁位與磁場關(guān)系

(7)

(8)

1.3 形函數(shù)一維積分

垂向空間域一維積分，將一維積分沿深度方向剖分成M個單元，在每個單元內(nèi)磁化強(qiáng)度滿足二次變化，筆者采用二次插值的形函數(shù)進(jìn)行積分，在磁場變化快的地方加密網(wǎng)格，磁場變化慢的地方剖分稀疏，既能保證計算精度，又能減少計算量。

對空間-波數(shù)混合域磁位積分公式(2)離散

(9)

sign (k))e-|k||z-z'|dz'

(10)

sign (z-z'))e-|k||z-z'|dz'

(11)

式中:zj、zm分別為單元積分的下限和上限。

將磁化強(qiáng)度表示成單元節(jié)點磁化強(qiáng)度的二次形函數(shù)

(12)

(13)

將式(12)～式(13)代入式(10)～式(11)中，整理可得

(14)

(15)

1.4 起伏地形條件下的快速算法

對于起伏地形條件下的二度體磁異常數(shù)值模擬，將地下研究區(qū)域海拔最低點與最高點之間剖分為n個深度節(jié)點(起伏地形模型見圖1)。常規(guī)算法需要計算海拔最高點與最低點之間n個深度節(jié)點所在水平測線的磁異常，然后通過插值求得起伏地形上的磁異常，計算量是水平地形條件下的n倍，計算量大，效率低。針對這一問題，筆者提出一種針對起伏地形條件下二度體磁異常數(shù)值模擬的快速算法。

以x分量積分為例，形函數(shù)單元積分公式為

e-|k||zp-z'|dz'

(16)

其中，zp為起伏地形上不同深度的節(jié)點。

垂向一維積分Ix(k,zp)積分區(qū)域由三部分組成(圖1所示中的(1)、(2)、(3))，故式(16)可以表示為

圖1 地形起伏模型示意圖Fig.1 Topographical model

Ix(k,zp)=Ix1(k,zp)e-|k|zp+Ix2(k,zp)e|k|zp+

Ix3(k,zp)e|k|zp

(17)

其中

(18)

其中：m1、m2、m3分別為三部分積分區(qū)域的單元個數(shù)。

由式(16)可知，①計算不同深度節(jié)點zp所在測線磁異常時，Ix3(k,zp)是相同的，即對于不同的深度節(jié)點zp，地形最低點以下的積分只需計算一次；②隨著深度節(jié)點zp從地形最低點z1向上移動，深度節(jié)點zp的積分Ix1(k,zp)、Ix2(k,zp)與深度節(jié)點zp-1的積分Ix1(k,zp-1)、Ix2(k,zp-1)存在積分區(qū)域重疊，重疊區(qū)域可只計算一次；③當(dāng)計算深度節(jié)點z1時，積分區(qū)域變?yōu)镮x2(k,z1)和Ix3(k,z1)，當(dāng)計算深度節(jié)點zn時，積分區(qū)域變?yōu)镮x1(k,zn)和Ix3(k,zn)兩部分。

根據(jù)式(16)積分特點，當(dāng)計算深度節(jié)點zp所在測線磁異常場時，可利用上一個深度節(jié)點zp-1的積分計算結(jié)果，避免了重復(fù)計算，提高了計算效率。快速算法主要流程如下：

1)計算起伏面最低點z1節(jié)點所在線磁異常場。由式(18)計算z1節(jié)點上下兩部分積分Ix2(k,z1)、Ix3(k,z1)，并對Ix3(k,z1)及Ix2(k,z1)每個單元積分結(jié)果進(jìn)行存儲。

2)計算下一個深度節(jié)點z2。由積分公式特點可知，z2節(jié)點的Ix3(k,z2)與z1節(jié)點的Ix3(k,z1)相同，可直接調(diào)用已存儲結(jié)果；z2節(jié)點的Ix2(k,z2)積分區(qū)域包含在z1節(jié)點Ix2(k,z1)積分區(qū)域以內(nèi)，每個單元計算結(jié)果已經(jīng)計算存儲，可直接調(diào)用；計算Ix1(k,z2)并存儲。

3)其他深度節(jié)點依次類推，直到計算出最高點zn節(jié)點所在平面磁異常場，然后通過插值得到起伏地形上磁異常場。

2 模型算例

在模型算例中，首先設(shè)計突變介質(zhì)模型用于驗證本文算法的正確性和可靠性，即采用基于4個高斯點的高斯FFT，計算了空間-波數(shù)混合域二度體磁場及其梯度異常，并與空間域解析解[15]對比，驗證本文算法的正確性和可靠性。其次通過對比復(fù)雜地形條件下的磁異常場快速算法與常規(guī)算法的計算時間與計算精度，驗證基于形函數(shù)積分快速算法的有效性和高效性，并分析了地形起伏最低點與最高點之間垂向剖分節(jié)點個數(shù)對磁異常場計算精度的影響規(guī)律。最后設(shè)計非規(guī)則復(fù)雜模型用于研究不同觀測高度的磁場分量響應(yīng)特征。

2.1 正確性驗證

設(shè)計如圖2所示的異常體模型。模擬區(qū)域大小為1 000 m×500 m，網(wǎng)格剖分為201×101，水平和垂向網(wǎng)格間距均為5 m；異常體模型位于研究區(qū)域中心位置，異常體寬為400 m，厚度為100 m，頂面埋深為200 m，底面埋深為300 m，異常體的磁化率為 0.01，地球正常磁場強(qiáng)度為45 000 nT，磁傾角為45°，磁偏角為0°。

圖2 異常體模型示意圖Fig.1 Anomalous body model

圖3為磁位和磁場的解析解、數(shù)值解以及地面z=0各節(jié)點的相對誤差。從圖3中可以看出：磁位和磁場的數(shù)值解與解析解吻合較好，地面節(jié)點磁位的最大相對誤差為0.48%，地面節(jié)點磁場Bx分量和Bz分量最大相對誤差均為0.52%。圖4為磁梯度張量的解析解和數(shù)值解，以及地面節(jié)點的相對誤差。從圖4可以看出，磁梯度張量的數(shù)值解與解析解吻合較好，地面節(jié)點磁梯度張量的最大相對誤差為0.22%。綜合圖3和圖4可以得出如下結(jié)論：磁位、磁場和磁梯度張量的數(shù)值解和解析解對相對誤差均很小，僅在模型突變點附近誤差略有增大，故驗證了本文算法的正確性，且計算結(jié)果表明其精度較高。

圖3 磁位、磁場數(shù)值解、解析解及地面節(jié)點相對誤差Fig.3 Numerical and analytical solutions of magnetic potential, fields and relative errors of ground nodes(a)磁位U解析解;(b)磁位U數(shù)值解;(c)磁位U相對誤差;(d)Bx解析解;(e)Bx數(shù)值解; (f)Bx相對誤差;(g)Bz解析解;(h)Bz數(shù)值解;(i)Bz相對誤差

圖4 磁梯度張量數(shù)值解、解析解及地面節(jié)點相對誤差Fig.4 Numerical and analytical solutions of magnetic gradient tensor and relative errors of ground nodes(a)Txx解析解;(b)Txx數(shù)值解;(c)Txx相對誤差;(d)Txz解析解;(e)Txz數(shù)值解; (f)Txz相對誤差;(g)Tzz解析解;(h)Tzz數(shù)值解;(i)Tzz相對誤差

2.2 計算精度與計算效率

設(shè)計如圖5所示起伏地形組合模型，驗證本文提出的快速算法對復(fù)雜地形的適應(yīng)性。地形起伏用式(19)所示的函數(shù)描述，模擬區(qū)域地形起伏最低點與最高點之間高差為2 km。

圖5 組合模型地形起伏示意圖Fig.5 Combined model with topography

x∈[-25000,25000]

(19)

起伏地形最低點以下計算區(qū)域為50 km×15 km，網(wǎng)格剖分為1001×501，水平網(wǎng)格間距為50 m，垂向網(wǎng)格間距為30 m，起伏地形最低點與最高點之間垂向剖分11、21、31、41、51節(jié)點。在計算區(qū)域分布著三個大小相同的異常體，異常體寬為6 km，厚度為3 km，頂面埋深距地形起伏最低點為1.5 km。左右兩個異常體的磁化率均為 0.01，中間異常體的磁化率為0.03；地球正常場的磁化強(qiáng)度為45 000 nT，磁偏角為0°，磁傾角為45°。

針對該模型，分析對比了起伏地形最低點與最高點之間垂向剖分節(jié)點個數(shù)對快速算法與常規(guī)算法計算效率和計算精度的影響，計算效率和計算精度統(tǒng)計結(jié)果如圖6所示。本文算法均采用fortran 95語言編寫,模型算例均在配置為CPU-Inter Core i5,主頻3.30 GHz,物理內(nèi)存32.00 GB的臺式機(jī)上串行測試得到。

圖6 快速算法的計算精度與計算效率圖Fig.6 Calculation accuracy and efficiency of the fast algorithm(a) 計算效率；(b)計算精度

從圖6(a)可以看出常規(guī)算法的計算時間隨著N(起伏地形最低點與最高點之間剖分節(jié)點數(shù))的增大近似線性增長，快速算法的計算時間隨著N的增大時間幾乎不變；快速算法計算起伏地形最低點與最高點深度方向剖分11個節(jié)點的磁異常時間為0.72 s，約為常規(guī)算法的1/5；剖分51個節(jié)點的磁異常時間為0.88 s，約為常規(guī)算法的1/17；隨著深度方向網(wǎng)格剖分節(jié)點個數(shù)N的增大，快速算法效率優(yōu)勢表現(xiàn)更為明顯。

圖6(b)為起伏地形上磁異常場及磁梯度張量相對均方根誤差[14](rrms)隨起伏地形最低點與最高點之間垂向剖分節(jié)點個數(shù)變化的規(guī)律[14]。從圖6(b)中可以看出：起伏地形上磁異常場及磁梯度張量的計算精度隨著剖分節(jié)點個數(shù)的增加而提高。起伏地形最低點與最高點之間剖分9個節(jié)點時，磁異常場及磁梯度張量的相對均方根誤差小于0.5%。當(dāng)剖分節(jié)點個數(shù)大于9時，磁異常場及梯度張量的相對均方根誤差幾乎不再變化。

為了兼顧計算精度與計算效率，圖7給出了起伏地形條件下最低點與最高點之間剖分15個節(jié)點的起伏地形上磁場和磁梯度張量數(shù)值解與解析解。由圖7可知：磁場與磁梯度張量不同分量的解析解與本文算法結(jié)果基本重合，地面磁異常場的最大相對誤差為0.50%，地面磁梯度張量的最大相對誤差為0.15%，表明本文提出的起伏地形條件下磁異?？焖偎惴ň雀?，速度快，適用性強(qiáng)。

圖7 起伏地形上磁場、磁梯度張量數(shù)值解與解析解Fig.7 Numerical and analytical solutions of magnetic fields and gradient tensor on undulating terrai(a)Bx;(b)Bz;(c)Txx;(d)Txz

2.3 復(fù)雜模型

為了進(jìn)一步檢驗本文算法對非規(guī)則復(fù)雜二度體磁異常計算的有效性，設(shè)計了如圖8所示的復(fù)雜地形模型。研究區(qū)域范圍均為：x方向從-200 m到200 m，z方向從-200 m到200 m。z軸鉛垂向上為正，網(wǎng)格剖分大小為401×201。圖8起伏地形異常體部分由山谷和山脊構(gòu)成。非規(guī)則異常體的磁化率κ為0.02，背景介質(zhì)的磁化率為“0”，地球磁場強(qiáng)度為45 000 nT，磁傾角為45°，磁偏角為0°。 這里分別在起伏地形和觀測面為40 m和100 m的高度上觀測，計算結(jié)果如圖9所示。從圖9中可以看出，在起伏地形上觀測的磁異常分量在地形突變點處有明顯的異常拐點，在山脊上存在明顯的正異常響應(yīng)，山谷處存在明顯的負(fù)異常響應(yīng)，且Bz分量的最大異常幅度明顯大于Bx分量最大異常幅度；隨著觀測高度的增加，磁異常場響應(yīng)逐漸變緩；山脊磁異常體正上方磁場Bx分量的正異常幅值較大，Bz分量的負(fù)異常幅值較大，山谷磁異常體正上方Bx和Bz分量的正負(fù)異常幅值均相對較小，且比較平滑，這是由于左邊山脊磁異常體的體積明顯大于右邊山谷磁異常體。通過研究非規(guī)則復(fù)雜模型起伏地形上和不同觀測高度的磁異常響應(yīng)特征，為磁異常反演成像和地質(zhì)解釋提供借鑒。

圖8 復(fù)雜地形模型Fig.8 The complex topographic model

圖9 不同高度磁場分量Fig.9 Results of magnetic field components at different heights(a)Bx;(b)Bz

3 結(jié)論

這里提出一種基于積分解的空間-波數(shù)混合域二度體磁異常數(shù)值模擬方法。該方法充分利用傅里葉變換方法和形函數(shù)積分的特點，將一個二維積分問題轉(zhuǎn)化為一維準(zhǔn)解析積分問題，有效提升了計算效率。數(shù)值試驗結(jié)果表明：該算法數(shù)值精度高、計算速度快、適應(yīng)性強(qiáng)，尤其適用于任意起伏地形條件下的磁異常及其梯度場的高效、高精度模擬，為二度體或三度體位場高效數(shù)值模擬提供了一種新的研究思路。