基于系統迭代法的線源瞬變電磁近區正演模擬

賴劉保，席振銖,*，王 鶴，侯海濤，朱偉國，劉愿愿，蔣 歡

(1. 中南大學 海洋礦產探測技術與裝備研究所，地球科學與信息物理學院，長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083; 3．長沙五維地科勘察技術有限責任公司，長沙 410205)

0 引言

瞬變電磁法(TEM)已廣泛應用于礦產資源調繪、水資源調查、地熱勘探[1]、巖溶探測以及油氣田勘查等領域[2]，國內、外幾十年的實踐工作說明了該方法在生產實踐中應用的有效性[3]。

根據TEM的場源性質，大體分為磁性源和電性源兩大類。一般來說磁性源的探測深度在1 km左右，不能滿足1 km～3 km的資源勘查、地熱探測、區域深部工程地質構造等方面的需求。目前，電性源瞬變電磁主要在遠區觀測，即長偏移距工作方式(LOTEM)，其采用1 km到3 km的接地長導線發射電磁場，在偏移距3 km～10 km主要觀測水平電場分量與垂直磁場分量[4]。前蘇聯及歐美國家廣泛將LOTEM應用于地熱勘探、地殼構造調查及地震勘探方法勘查時效果不好的地區，并形成了一套較完整的系統，但是在遠區觀測時存在信號強度弱、體積效應大等問題，而在近區線源觀測則有以下優點：①信噪比高、體積效應小、探測靈敏度高；②電壓衰減速度慢，發送磁矩大，探測深度相對較大，裝置輕便，適合山區作業等優點[5]。

當前線源近區的研究多局限于二維數值模擬[6]，一維反演解釋[7]。隨著計算機技術的發展，瞬變電磁的三維正演快速計算模擬可以實現。與需要全空間離散的有限元[8]、有限差分[9-10]所不同的是：積分方程法只需對異常體進行離散，具有離散數目少，求解精度高等優點[11-12]。與較成熟的LOTEM相比，國內、外學者對于線源近區三維正演研究相對較少。主要問題在于三維計算中消耗計算機大量內存和時間，尤其在求解大型稀疏矩陣時，系數矩陣的條件數很大，很難獲得高精度的結果[13]。同時在近場區，線源不能采用電偶極子近似[14]，場源剖分時也增加了計算量。為了同時解決三維計算的速度和精度問題，作者引入了系統迭代法求解積分方程來降低系數矩陣條件數，節省計算機的內存及計算時間。

作者首先介紹了層狀介質中三維電性異常體的積分方程，給出了第二類Fredholm方程，重點闡述了系統迭代法的求解過程，然后利用數字濾波法將頻率域響應轉換到時間域；最后在均勻半空間中將正演結果與解析解對比，檢驗算法的有效性，并結合簡單模型與復雜模型的計算，證明了該方法對異常體探測的可靠性。

1 基本理論

1.1 積分方程法

三維模型如圖1所示，水平層狀電導率的張量形式：σb=diag(σb,σb,σb)，其中σb= σb′ +iωεb為大地的復電導率。層狀大地與異常體的電導率之差σ=σaI-σb，其中I為3×3的單位矩陣，σa為異常體電導率，則三維異常體的散射電流為Js=σ·E，E為異常體中的電場強度。假設入射場在時諧因子eiwt，外加電流激勵下垂直向下傳播，且忽略位移電流。根據MAXWELL方程及相應的電磁格林張量理論，異常體V滿足第二類Fredholm積分方程：

(1)

其中θ=diag(1/σa-σb,1/σa-σb,1/σa-σb)；Ep為一次電場；GE=-iwμ0A+(1/σb)▽▽·A為電場格林函數[15]。

將異常體區域V剖分為M個小單元，得到如下矩陣方程：

[Γ][Js]=[Ep]

(2)

圖1 三維地電模型示意圖Fig．1 Sketch of 3D geoelectric model

1.2 線源的處理

由于觀測點常常比較靠近接地導線源，需把發射源看做多個偶極的疊加，計算電磁場分量時將偶極場表達式沿線源積分，文獻[14]給出了層狀大地表面的電磁場表達式。假設導線源的中心位于原點，沿水平軸向兩側延伸至-L和L。

水平電場公式為：

R=[(x-x′)2+y2]1/2;

R1=[(x+L)2+y2]1/2;

R2=[(x-L)2+y2]1/2。

垂直磁場公式為：

1.3 系統迭代法解矩陣方程

如果把三維異常體V分成m塊，V1、…、Vi、…、Vm，然后再將第i塊異常體Vi剖分為Mi個小單元，式(1)變為：

(3)

式(3)表明，Vi區域外的異常體對Vi的影響相當于外加激發源的作用。

若把整個異常體區域看成一個系統，則每個塊狀異常體在各自區域內形成子系統。我們可以單獨的處理各個異常體，獲得子系統間的相互影響，從而得到整個系統的值，這就是系統迭代法。該方法只需要保存處于計算中的單個異常體的離散單元，而不需要存儲整個異常體的全部離散，可以節省時間，節約內存，而且隨著剖分越來越細，系數矩陣的條件數越來越大，線性方程組的解越來越不穩定。系統迭代法只跟剖分后異常體的個數有關，而與未知數無關，這樣大大降低了矩陣的條件數，提高了精度。

根據矩陣理論，式(3)相當于將散射矩陣Γ離散成塊矩陣,式(2)變為：

(4)

用高斯-賽德爾(GS)迭代法求解式(4)：

(5)

為提高收斂速度采用逐次超松弛(SOR)迭代法，式(5)變為：

(6)

其中ω為松弛因子且1≤ω<2。當w=1時，SOR法為GS法，即式(6)等于式(5)。

選擇合適的松弛因子可以保證SOR的收斂。經過多次驗算，本文w=1.25。

Γii可以通過LU或者Cholesky分解并將離散矩陣用于迭代，Ep可以用快速漢克爾變換計算，異常體部分可以用文獻[15]所提到的格林函數法求解。如果迭代的相對誤差小于ε，則矩陣方程收斂，迭代終止。

求得Js后，根據并矢格林函數，二次電磁場分別為：

1.4 時頻轉換

在頻率域中求得電磁響應后，利用安德森的滯后褶積的快速數字濾波法[16]計算如下正弦與余弦變換。

Im[E(ω)]和Im[H(ω)]分別是頻率域中電、磁場分量的虛部。求得磁場強度對時間的導數?H(t)/?t后，再乘真空磁導率μ0=4π×10-7Wb/(A·m)， 可得磁感應強度對時間的導數 ?B(t)/?t。

2 算法檢驗

為了檢驗本算法的正確性，本研究采用線源赤道軸線上的解析式與積分方程法得到的數值解對比。線源赤道平面如圖2。

圖2 線源赤道平面圖Fig．2 Plan view of the line source

假設均勻半空間的電阻率為100 Ω·m,線源長1 km，激發電流1 A,偏移距1 000 m處接收。積分方程計算中，在頻率域范圍0.01 Hz～1 000 000 Hz中計算了49個頻點，每個數量級取6個頻點，向時間域0.01 ms～1 000 ms轉換中，共計算了31個三維體的瞬態響應，從而保證了計算精度。

從圖3中可以看出，不管垂直感應磁場時間倒數還是水平電場強度的衰減曲線結果，都與解析解十分吻合。這也驗證本文的數值模擬精度高，計算可靠。

圖3 半空間中的解析解與數值解的對比Fig．3 Contrast between numeric solution and analytic solution in a half-space(a) dBz/dt衰減曲線; (b)水平電場Ey衰減曲線

3 模型計算

為了說明線源對異常體的反應效果，分別設計簡單模型和復雜模型對瞬態電磁響應加以分析。計算區域示意圖4顯示：場源(箭頭)中心位于原點，長1 000 m；異常體位于5 000 m×4 000 m的地面下，埋深750 m；發送電流1 A，三角形為測點。兩個模型的以上參數相同，分別計算垂直感應時間導數與平行于場源方向的水平電場Ey，然后作出等值線圖與隨時間的衰減曲線圖。

圖4 三維異常體的平面圖和斷面圖位置A、B、C、D分別對應于偏移距r=-250 m、250 m、1 500 m、3 000 mFig．4 Plan view and cross-section for 3D conductive prism(a)X-Y平面圖; (b)X-Z斷面圖

3.1 低阻模型

圖5是一個低阻異常體存在于均勻半空間中的模型。半空間的電阻率為50 Ω·m，異常體電阻率為5 Ω·m，尺寸1 000 m×1 000 m×500 m。將異常體剖分為10×10×10個網格單元。

圖5 半空間中的單個三維異常體模型Fig．5 Model of a single anomalous 3-D body in a half-space

圖7表示垂直感應磁場的時間導數在不同觀測點隨時間的衰減曲線。在x=-250 m 處，前支出現反向電流，瞬態值為負，在120 ms方向發生逆轉變為正值；x=250 m 剛好相反，只是在80 ms 時刻方向發生偏轉，因為B點離異常體更近。這個符號變化特征可以為縱向上確定異常體提供參考[19]，經計算，在80 ms磁場電流擴散到異常體。觀測點在橫縱向上離異常體750 m左右，而偏移距只有250 m。這表明近區也能探測較遠的異常體，隨著偏移距的增加，曲線形態逐漸平緩，這點和文獻[20]給出的一維計算結果相同。在x=1 500 m、3 000 m時，如同磁性源及LOTEM，感應磁場時間導數隨時間按t-5/2衰減，這說明了本文方法的可靠性。

圖6 不同時間的水平電場Ey(V/m)等值線Fig．6 Contour maps of horizontal electrical field Ey at different times(a)t=1 ms；(b)t=3 ms；(c)t=10 ms；(d)t=30 ms；(e)t=100 ms；(f)t=300 ms；箭頭表示場源，長方形表示異常體的平面位置

圖7 不同偏移距的dBz/dt曲線Fig．7 Curves with different offset

3.2 復雜模型

假設電阻率為50 Ω·m 的均勻半空間中分別賦存高、低阻兩個異常體：高阻電阻值為500 Ω·m，尺寸1 000 m×1 000 m×500 m,中心點位于(1 500,-1 000,750)；低阻值為5 Ω·m，尺寸1 000 m×1 000 m×500 m，中心位于(1 500,1 000,750)。在計算中將兩個異常體各剖分為10×10×10個網格單元，共 2 000個網格(圖8)。

圖8 高阻和低阻兩個異常體的三維模型Fig．8 Schematic drawing of 3D model with two anomalies in homogeneous half space

水平電場切片圖9表明：剛開始，電流主要集中在發射源附近，這點和低阻模型相同，也符合電磁波的傳播規律；在100 ms時，在異常體位置正上方出現高、低電場值，如圖中矩形所示；當t=300 ms時，異常體的位置反應得更加明顯，而且不隨時間變化。與圖6(f)類似，圖9(f)的兩個極值對應異常體邊界，而且高阻異常體對應高電場值，低阻異常體對應低電場值。因此電場值不僅能反應低阻的存在，而且對高阻體也有較高的分辨率，這是其他裝置所不具備的優勢在圖10中，在低阻模型的同一深度的加入高阻體后，感應磁場時間導數的衰減曲線和單個低阻體的曲線基本相同，隨偏移距增加，曲線仍然變緩，只是在x=250 m、t=100 ms瞬態方向發生逆轉，比后者晚了20 ms。

以上數值計算是在雙核2.4GHz，內存2.0 G的個人筆記本電腦上實現的，在精度達到0.000 1的情況下，低阻模型只需8 min，高低阻模型只需25 min。如果剖分網格單元變小，需要的時間更少，說明系統迭代法能節省時間。

4 結論

本次研究采用基于系統迭代的積分方程法計算了線源近區三維瞬變場值。通過與解析解的對比及兩個算例的模擬，得出如下結論：

(1)系統迭代法求解積分方程，不僅節省了內存和計算時間，而且降低了離散矩陣的條件數，提高了計算精度。

(2)水平電場不僅對低阻體有良好的分辨力，而且對高阻也反應靈敏。低阻異常體對應低場值，高阻體對應高場值，同時其極值對應異常體的左右邊界，這為今后更好利用電場分量提供了理論基礎。

(3)異常體的感應磁場對時間導數的方向發生逆轉，而且隨偏移距的增加，曲線形態變緩。

圖9 高低阻模型不同時間水平電場Ey(V/m)等值線圖Fig．9 Contour maps of horizontal electrical field Ey at different times(a)t=1 ms；(b)t=3 ms；(c)t=10 ms；(d)t=30 ms；(e)t=100 ms；(f)t=300 ms箭頭表示場源，長方形表示異常體的平面位置

參考文獻：

[1] 李文堯，廖忠.瞬變電磁法在騰沖尋找地熱中的應用[J].物探與化探，2002，26(5):368-371.

圖10 不同偏移距的dBz/dt曲線Fig．10 Curves with different offset

[2] 牛之璉.時間域電磁法原理[M].長沙：中南大學出版社，2007.

[3] 李貅.瞬變電磁測深的理論與應用[M]．西安：陜西科學技術出版社，2002.

[4] 陳明生.電偶源瞬變電磁測深研究(二)[J].煤田地質與勘探，1999，27(2)：54-57.

[5] 趙福元，嚴良俊，何展翔，等.線源時間域電測深全期全區視電阻率求解[J].地震地質，2010,32(3):473-481.

[6] 嚴良俊，徐世浙，陳小斌，等.線源二維瞬變電磁場的正演計算新方法[J].煤田地質與勘探，2004,32(5)：58-61.

[7] 劉瑞德，梅巖輝，黃力軍.電性源瞬變電磁快速反演解釋方法[J].工程地球物理學報，2005,2(1):18-21.

[8] 徐志鋒，吳小平.可控源電磁三維頻率域有限元模擬[J].地球物理學報，2010,53(8)：1931-1939.

[9] YUTAKA SASAKI, SEONG-JUN CHO.3D finite-difference modeling of time-domain electromagnetic data for mineral exploration[C].Proceedings of the 10th SEGJ Interational Symposium,2011:235-238.

[10] TSILI WANG, ALAN C TRIPP,Gerald W HOHMANN. Studying the TEM response of a 3-D conductor at a geologucal contact using the FDTD method[J].Geophysics,1995,60(4):1265-1269.

[11] GERALD W,HOHMANN.Three dimensio-nal em modeling[J].Geophysical Surveys,1983,6(1983):27-53.

[12] PHILIPE,WANNAMAKER,GERALDW.Hohmannand Willian A.SanFilipo.Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations[J].Geophysics,1984,49(1):60-74.

[13] ZONGHOU XIONG.Electromagnetic modeling of 3-D structures by the method of system iteration using integral equations[J].Geophysics,1992,57(12):1556-1561.

[14] 米薩克 N.納比吉安.勘察地球物理電磁法 第一卷:理論[M] .趙經祥，王艷君譯.北京：地質出版社，1992.

[15] ZONGHOU XIONG. Electromagnetic fields of electric dipoles embedded in a stratified anisotropic earth[J]. Geophysics,1989,54(12):1643-1646.

[16] WALTER L,ANDERSON. Fast Hankel transforms using realted and lagged convolutions[C]. ACM Transactions on Mathematical Software,1982,8(4):344-368.

[17] ORISTA GLIO ML..Diffusion of electromagnetic fields into the earth from a line source of current[J].Geophysics,1982,47(11):1585-1592.

[18] NABIGHIAN M N, ORISTAGLIO M L On the approximation of finite loop sources by two-dimensional line source[J].Geophysics,1984,49(7):1027-1029.

[19] BRIANM.GUNDERSION,GREGORY A.NEWMAN,Gerald W.HOHMANN. Three-dimensional transient electromagnetic response for a grounded source[J].Geophysics,1986,51(11):2117-2130.

[20] 薛國強，陳衛營，周楠楠，等.接地源瞬變電磁短偏移距深部探測技術[J].地球物理學報，2013,56(1)：255-261.