基于積分解的空間-波數混合域二度體磁異常數值模擬

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

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

0 引言

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

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

1 方法原理

1.1 空間-波數混合域磁位

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

(1)

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

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

(2)

(3)

(4)

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

1.2 磁場和磁梯度張量

根據空間域磁場、磁梯度張量與磁位的關系

(5)

式中，B為空間域磁感應強度。

(6)

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

(7)

(8)

1.3 形函數一維積分

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

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

(9)

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

(10)

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

(11)

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

將磁化強度表示成單元節點磁化強度的二次形函數

(12)

(13)

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

(14)

(15)

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

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

以x分量積分為例，形函數單元積分公式為

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

(16)

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

垂向一維積分Ix(k,zp)積分區域由三部分組成(圖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分別為三部分積分區域的單元個數。

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

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

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

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

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

2 模型算例

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

2.1 正確性驗證

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

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

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

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

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

2.2 計算精度與計算效率

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

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

x∈[-25000,25000]

(19)

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

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

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

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

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

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

圖7 起伏地形上磁場、磁梯度張量數值解與解析解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 復雜模型

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

圖8 復雜地形模型Fig.8 The complex topographic model

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

3 結論

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