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

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

0 引言

瞬變電磁法(TEM)已廣泛應(yīng)用于礦產(chǎn)資源調(diào)繪、水資源調(diào)查、地?zé)峥碧絒1]、巖溶探測以及油氣田勘查等領(lǐng)域[2]，國內(nèi)、外幾十年的實(shí)踐工作說明了該方法在生產(chǎn)實(shí)踐中應(yīng)用的有效性[3]。

根據(jù)TEM的場源性質(zhì)，大體分為磁性源和電性源兩大類。一般來說磁性源的探測深度在1 km左右，不能滿足1 km～3 km的資源勘查、地?zé)崽綔y、區(qū)域深部工程地質(zhì)構(gòu)造等方面的需求。目前，電性源瞬變電磁主要在遠(yuǎn)區(qū)觀測，即長偏移距工作方式(LOTEM)，其采用1 km到3 km的接地長導(dǎo)線發(fā)射電磁場，在偏移距3 km～10 km主要觀測水平電場分量與垂直磁場分量[4]。前蘇聯(lián)及歐美國家廣泛將LOTEM應(yīng)用于地?zé)峥碧健⒌貧?gòu)造調(diào)查及地震勘探方法勘查時(shí)效果不好的地區(qū)，并形成了一套較完整的系統(tǒng)，但是在遠(yuǎn)區(qū)觀測時(shí)存在信號(hào)強(qiáng)度弱、體積效應(yīng)大等問題，而在近區(qū)線源觀測則有以下優(yōu)點(diǎn)：①信噪比高、體積效應(yīng)小、探測靈敏度高；②電壓衰減速度慢，發(fā)送磁矩大，探測深度相對(duì)較大，裝置輕便，適合山區(qū)作業(yè)等優(yōu)點(diǎn)[5]。

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

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

1 基本理論

1.1 積分方程法

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

(1)

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

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

[Γ][Js]=[Ep]

(2)

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

1.2 線源的處理

由于觀測點(diǎn)常常比較靠近接地導(dǎo)線源，需把發(fā)射源看做多個(gè)偶極的疊加，計(jì)算電磁場分量時(shí)將偶極場表達(dá)式沿線源積分，文獻(xiàn)[14]給出了層狀大地表面的電磁場表達(dá)式。假設(shè)導(dǎo)線源的中心位于原點(diǎn)，沿水平軸向兩側(cè)延伸至-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 系統(tǒng)迭代法解矩陣方程

如果把三維異常體V分成m塊，V1、…、Vi、…、Vm，然后再將第i塊異常體Vi剖分為Mi個(gè)小單元，式(1)變?yōu)椋?/p>

(3)

式(3)表明，Vi區(qū)域外的異常體對(duì)Vi的影響相當(dāng)于外加激發(fā)源的作用。

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

根據(jù)矩陣?yán)碚摚?3)相當(dāng)于將散射矩陣Γ離散成塊矩陣,式(2)變?yōu)椋?/p>

(4)

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

(5)

為提高收斂速度采用逐次超松弛(SOR)迭代法，式(5)變?yōu)椋?/p>

(6)

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

選擇合適的松弛因子可以保證SOR的收斂。經(jīng)過多次驗(yàn)算，本文w=1.25。

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

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

1.4 時(shí)頻轉(zhuǎn)換

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

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

2 算法檢驗(yàn)

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

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

假設(shè)均勻半空間的電阻率為100 Ω·m,線源長1 km，激發(fā)電流1 A,偏移距1 000 m處接收。積分方程計(jì)算中，在頻率域范圍0.01 Hz～1 000 000 Hz中計(jì)算了49個(gè)頻點(diǎn)，每個(gè)數(shù)量級(jí)取6個(gè)頻點(diǎn)，向時(shí)間域0.01 ms～1 000 ms轉(zhuǎn)換中，共計(jì)算了31個(gè)三維體的瞬態(tài)響應(yīng)，從而保證了計(jì)算精度。

從圖3中可以看出，不管垂直感應(yīng)磁場時(shí)間倒數(shù)還是水平電場強(qiáng)度的衰減曲線結(jié)果，都與解析解十分吻合。這也驗(yàn)證本文的數(shù)值模擬精度高，計(jì)算可靠。

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

3 模型計(jì)算

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

圖4 三維異常體的平面圖和斷面圖位置A、B、C、D分別對(duì)應(yīng)于偏移距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是一個(gè)低阻異常體存在于均勻半空間中的模型。半空間的電阻率為50 Ω·m，異常體電阻率為5 Ω·m，尺寸1 000 m×1 000 m×500 m。將異常體剖分為10×10×10個(gè)網(wǎng)格單元。

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

圖7表示垂直感應(yīng)磁場的時(shí)間導(dǎo)數(shù)在不同觀測點(diǎn)隨時(shí)間的衰減曲線。在x=-250 m 處，前支出現(xiàn)反向電流，瞬態(tài)值為負(fù)，在120 ms方向發(fā)生逆轉(zhuǎn)變?yōu)檎?；x=250 m 剛好相反，只是在80 ms 時(shí)刻方向發(fā)生偏轉(zhuǎn)，因?yàn)锽點(diǎn)離異常體更近。這個(gè)符號(hào)變化特征可以為縱向上確定異常體提供參考[19]，經(jīng)計(jì)算，在80 ms磁場電流擴(kuò)散到異常體。觀測點(diǎn)在橫縱向上離異常體750 m左右，而偏移距只有250 m。這表明近區(qū)也能探測較遠(yuǎn)的異常體，隨著偏移距的增加，曲線形態(tài)逐漸平緩，這點(diǎn)和文獻(xiàn)[20]給出的一維計(jì)算結(jié)果相同。在x=1 500 m、3 000 m時(shí)，如同磁性源及LOTEM，感應(yīng)磁場時(shí)間導(dǎo)數(shù)隨時(shí)間按t-5/2衰減，這說明了本文方法的可靠性。

圖6 不同時(shí)間的水平電場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 復(fù)雜模型

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

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

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

以上數(shù)值計(jì)算是在雙核2.4GHz，內(nèi)存2.0 G的個(gè)人筆記本電腦上實(shí)現(xiàn)的，在精度達(dá)到0.000 1的情況下，低阻模型只需8 min，高低阻模型只需25 min。如果剖分網(wǎng)格單元變小，需要的時(shí)間更少，說明系統(tǒng)迭代法能節(jié)省時(shí)間。

4 結(jié)論

本次研究采用基于系統(tǒng)迭代的積分方程法計(jì)算了線源近區(qū)三維瞬變場值。通過與解析解的對(duì)比及兩個(gè)算例的模擬，得出如下結(jié)論：

(1)系統(tǒng)迭代法求解積分方程，不僅節(jié)省了內(nèi)存和計(jì)算時(shí)間，而且降低了離散矩陣的條件數(shù)，提高了計(jì)算精度。

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

(3)異常體的感應(yīng)磁場對(duì)時(shí)間導(dǎo)數(shù)的方向發(fā)生逆轉(zhuǎn)，而且隨偏移距的增加，曲線形態(tài)變緩。

圖9 高低阻模型不同時(shí)間水平電場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箭頭表示場源，長方形表示異常體的平面位置

參考文獻(xiàn)：

[1] 李文堯，廖忠.瞬變電磁法在騰沖尋找地?zé)嶂械膽?yīng)用[J].物探與化探，2002，26(5):368-371.

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

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

[3] 李貅.瞬變電磁測深的理論與應(yīng)用[M]．西安：陜西科學(xué)技術(shù)出版社，2002.

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

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

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

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

[8] 徐志鋒，吳小平.可控源電磁三維頻率域有限元模擬[J].地球物理學(xué)報(bào)，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] .趙經(jīng)祥，王艷君譯.北京：地質(zhì)出版社，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] 薛國強(qiáng)，陳衛(wèi)營，周楠楠，等.接地源瞬變電磁短偏移距深部探測技術(shù)[J].地球物理學(xué)報(bào)，2013,56(1)：255-261.