積分方程法模擬層狀介質(zhì)中薄板的瞬態(tài)響應(yīng)

劉愿愿，席振銖，蔣 歡，王 鶴

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083；2.中南大學(xué) 海洋礦產(chǎn)探測(cè)技術(shù)與裝備研究所，湖南 長(zhǎng)沙 410083

?

積分方程法模擬層狀介質(zhì)中薄板的瞬態(tài)響應(yīng)

劉愿愿1,2，席振銖1,2，蔣 歡1,2，王 鶴1,2

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083；2.中南大學(xué) 海洋礦產(chǎn)探測(cè)技術(shù)與裝備研究所，湖南 長(zhǎng)沙 410083

瞬變電磁法高維正反演技術(shù)發(fā)展緩慢，主要是由于正演方法不夠成熟，計(jì)算速度慢，精度低，達(dá)不到反演的要求。文闡明了薄板的定義，采用積分方程法模擬層狀介質(zhì)中薄板的瞬態(tài)響應(yīng)，簡(jiǎn)化了網(wǎng)格剖分，控制矩陣大小，采取直接解法求取系數(shù)矩陣，與均勻半空間的解析解對(duì)比驗(yàn)證了該算法和程序的可靠性，通過(guò)誤差分析體現(xiàn)了該算法的高精度。同時(shí)，研究了覆蓋層與基巖對(duì)薄板響應(yīng)結(jié)果的影響，發(fā)現(xiàn)高阻覆蓋層和低阻覆蓋層對(duì)薄板的響應(yīng)影響都很小，高阻基巖會(huì)擴(kuò)大薄板的響應(yīng)，低阻基巖會(huì)屏蔽薄板的響應(yīng)。該研究為瞬變電磁二維反演奠定了基礎(chǔ)，也為以后生產(chǎn)中薄板解釋提供了理論借鑒。

瞬變電磁法；積分方程法；層狀介質(zhì)；薄板

1 引 言

瞬變電磁法在一次脈沖磁場(chǎng)的間歇期測(cè)量二次場(chǎng)，簡(jiǎn)化了對(duì)異常響應(yīng)的研究，因此具有更高的分辨率和靈敏度，近年來(lái)已成為物探熱門方法之一[1,2]。但是和MT，CSAMT等相比，瞬變電磁法正反演大多都還停留在一維的水平，高維正反演尚不成熟[3-5]。

正演速度是制約反演速度的關(guān)鍵因素[6]。有限單元法需要進(jìn)行離散化，網(wǎng)格剖分?jǐn)?shù)目較多，計(jì)算速度慢，對(duì)內(nèi)存要求也很高[7,8]。電磁波速度變化范圍比較大，因此時(shí)間步長(zhǎng)一般比較小，有限差分法需要迭代很多次才能完成正演過(guò)程，速度同樣較慢[9-11]。

積分方程法是最早實(shí)現(xiàn)瞬變電磁高維模擬的數(shù)值算法，因?yàn)樗挥?jì)算小體積異常區(qū)的場(chǎng)，計(jì)算速度快。1985年，SanFilipo和Hohmann首先利用時(shí)間域積分方程法求出均勻半空間中的對(duì)稱異常體的三維響應(yīng)，源是矩形的發(fā)射回線[12]。1986年至1988年，Newman和Hohmann等學(xué)者利用積分方程求得頻率域時(shí)間響應(yīng)，然后通過(guò)傅式變換轉(zhuǎn)到時(shí)間域的響應(yīng)[13]。國(guó)內(nèi)的學(xué)者主要致力于研究電磁張量格林函數(shù)的數(shù)值計(jì)算。2007年，魏寶君等人提出了并矢格林函數(shù)的遞推算法，引發(fā)了人們對(duì)早先使用的計(jì)算公式的修正與校驗(yàn)[14]。2009年，陳桂波系統(tǒng)地推導(dǎo)了積分方程的計(jì)算公式，并應(yīng)用于各向異性地層中進(jìn)行電磁場(chǎng)三維模擬[15]。2014年，胡俊華利用積分方程法進(jìn)行三維正演[16]，賴劉保等人利用積分方程法模擬接地線源近區(qū)的三維薄板的響應(yīng)[17]，說(shuō)明積分方程法在瞬變電磁響應(yīng)計(jì)算中日益成熟。

為了提高正演速度，本文在闡明薄板概念的基礎(chǔ)上，采用積分方程法模擬層狀介質(zhì)下任意位置任意形態(tài)薄板的瞬態(tài)響應(yīng)，并采用LU分解法求解矩陣。通過(guò)與均勻半空間的對(duì)比驗(yàn)證了該程序的正確性和精確度，研究了層狀介質(zhì)中覆蓋層和基巖對(duì)于薄板響應(yīng)值的影響。

2 理論基礎(chǔ)

薄板，從理論上來(lái)說(shuō)是厚度為0的導(dǎo)電板，這意味著感應(yīng)電流將會(huì)集中在薄板所在的那個(gè)面上。實(shí)際上，薄板的定義必須根據(jù)實(shí)際情況來(lái)定。5 m厚的薄板在早期或者高頻階段不可以被當(dāng)做薄板計(jì)算，但是在晚期或者低頻階段，40 m厚的薄板也可以當(dāng)做薄板來(lái)模擬，因?yàn)閷?dǎo)電覆蓋層會(huì)濾掉高頻，從而擴(kuò)大薄板的定義范圍。

由電磁理論可知，任意點(diǎn)r的電磁場(chǎng)E或者H可表示為源在圍巖中的背景場(chǎng)Eb(Hb)和薄板引起的散射場(chǎng)Es(Hs)的疊加：

E(r)=Eb(r)+Es(r)

(1-1)

H(r)=Hb(r)+Hs(r)

(1-2)

根據(jù)并矢格林函數(shù)理論，散射磁場(chǎng)可表示為：

(2)

(3)

Js(r)=Δσ(r) ·E(r)

(4)

其中，Δσ是圍巖與薄板電導(dǎo)率之差。把式(2)和式(4)代入式(1)，可得：

(5)

將異常體所在區(qū)域進(jìn)行網(wǎng)格剖分，并假定每個(gè)單元內(nèi)電場(chǎng)與電導(dǎo)率等于單元中心的值，則式(5)可寫為：

(6)

薄板是二維體，因此系數(shù)矩陣一般不會(huì)很大，采取LU分解法求解式(6)，可求出空間中任意位置的磁場(chǎng)。

3 正確性驗(yàn)證

電阻率為100歐姆的均勻半空間中磁通密度B隨時(shí)間t變化的解析表達(dá)式為：

(7)

同樣，在電阻率為100Ω的均勻半空間指定一塊100m×100m的區(qū)域看作薄板，分別是X方向(-50m，50m)，Z方向是(-50m，-150m)，在地面布置測(cè)線一條，(-200m,0m)到(200m,0m)，一共11個(gè)點(diǎn)，發(fā)射線圈邊長(zhǎng)100m×100m，電流為5A，使用該算法進(jìn)行正演模擬，模擬結(jié)果同解析解進(jìn)行對(duì)比，對(duì)比結(jié)果見(jiàn)圖1，具體誤差見(jiàn)表1。

由圖1可見(jiàn)，積分方程法計(jì)算的響應(yīng)值和解析解的響應(yīng)值得出的曲線形態(tài)完全一致，驗(yàn)證了該算法的正確性。由表1可見(jiàn)，除開(kāi)早期的幾個(gè)節(jié)點(diǎn)，該算法的誤差基本控制在3%以內(nèi)，確保了該算法的高精度，可以計(jì)算更為復(fù)雜的模型。

圖1 均勻半空間解析解與積分方程法對(duì)比圖Fig.1 Comparison between the analytical solution and the integral equation method for the homogeneous half space

時(shí)間/ms誤差百分比/%0.30631.110.40422.920.50217.110.60012.810.6989.700.8455.191.0412.641.2371.091.4330.121.6290.612.3152.133.0991.773.4911.764.8631.785.6471.376.4311.147.2150.738.3911.539.9590.8911.5270.6014.6630.2917.0150.2520.1510.68

4 實(shí)例模擬

4.1 兩層介質(zhì)下直立薄板模擬

為了清晰直觀地反應(yīng)層狀介質(zhì)對(duì)薄板的響應(yīng)，首先定義直立薄板的電阻率為1 Ω，沿X方向厚度為0，位于零點(diǎn)處，Y方向是(-200 m，200 m)，Z方向是(-25 m，-75 m)，在地表Z=0的平面布置測(cè)線一條，沿(-400 m，0 m)至(400 m，0 m)，一共41個(gè)測(cè)點(diǎn)，發(fā)射線圈邊長(zhǎng)為100 m×100 m，電流為5 A，然后設(shè)定三種基巖模型：①是電阻率為100 Ω的均勻半空間；②是電阻率為100 Ω的基巖上覆20 m厚電阻率為1 000 Ω的覆蓋層；③是電阻率為1 000 Ω的基巖上覆80 m厚的電阻率為100 Ω的覆蓋層。模型示意圖見(jiàn)圖2，Z方向和X方向的響應(yīng)分別見(jiàn)圖3～圖6。

由圖3至圖6可見(jiàn)，模型1和模型2的Z方向和X方向的響應(yīng)在0.5 ms和2 ms曲線形態(tài)完全一致，數(shù)值相差也很小，可得高阻覆蓋層對(duì)于薄板響應(yīng)影響很小，幾乎可以忽略。但是模型3靠近薄板時(shí)響應(yīng)值變大，曲線光滑性明顯變差，極值點(diǎn)處的響應(yīng)更加尖銳，且越靠近晚期，區(qū)別更加明顯，可見(jiàn)低阻基巖會(huì)增大上方介質(zhì)中的薄板響應(yīng)值，并且增幅會(huì)隨著關(guān)斷時(shí)間的增加而增加。

圖2 直立薄板模型示意圖Fig.2 Sketch map of vertical thin plate

4.2 兩層狀介質(zhì)下水平薄板模擬

定義水平薄板的電阻率為1 Ω，沿Z方向厚度為0，位于Z=-50 m處，Y方向是(-200 m，200 m)，X方向是(-25 m，25 m)，在地表Z=0的平面布置測(cè)線一條，沿(-400 m，0 m)至(400 m，0 m)，一共41個(gè)測(cè)點(diǎn)，發(fā)射線圈邊長(zhǎng)為100 m×100 m，電流為5 A，然后設(shè)定三種基巖模型：①是電阻率為1 000 Ω的均勻半空間；②是電阻率為1 000 Ω的基巖上覆20 m厚電阻率為100 Ω的覆蓋層；③是電阻率為100 Ω的基巖上覆80 m厚的電阻率為1 000 Ω的覆蓋層。模型示意圖見(jiàn)圖7，Z方向和X方向的響應(yīng)分別見(jiàn)圖8～圖11。

圖3 t=0.5 ms直立薄板Z分量響應(yīng)圖Fig.3 Z component response diagram of vertical thin plate at t=0.5 ms

圖4 t=2 ms直立薄板Z分量響應(yīng)圖Fig.4 Z component response diagram of vertical thin plate at t=2 ms

圖5 t=0.5 ms直立薄板X分量響應(yīng)圖Fig.5 X component response diagram of vertical thin plate at t=0.5 ms

圖6 t=2 ms直立薄板X分量響應(yīng)圖Fig.6 X component response diagram of vertical thin plate at t=2 ms

圖7 水平薄板模型示意圖Fig.7 Sketch map of horizontal thin plate model

圖8 t=0.5 ms水平薄板Z分量響應(yīng)圖Fig.8 Z component response diagram of horizontal thin plate at t=0.5 ms

圖9 t=2 ms水平薄板Z分量響應(yīng)圖 Fig.9 Z component response diagram of horizontal thin plate at t=2 ms

圖10 t=0.5 ms水平薄板X分量響應(yīng)圖Fig.10 X component response diagram of horizontal thin plate at t=0.5 ms

圖11 t=2 ms水平薄板X分量響應(yīng)圖Fig.11 X component response diagram of horizontal thin plate at t=2 ms

由圖8至圖11可見(jiàn)，模型4和模型5的Z方向和X方向在0.5 ms和2 ms時(shí)曲線基本重合，響應(yīng)值相差很小，說(shuō)明低阻覆蓋層對(duì)于薄板的影響可以忽略。但是模型6基巖是低阻，所以Z方向背景場(chǎng)較大，但是Z方向和X方向薄板的響應(yīng)值很小，和模型4和模型5的響應(yīng)值相差約1個(gè)數(shù)據(jù)級(jí)，并且隨著關(guān)斷時(shí)間的增加，背景場(chǎng)逐步趨于重合，但是薄板響應(yīng)值的差距不會(huì)減少，依然十分明顯?？梢?jiàn)低阻基巖會(huì)部分屏蔽掉上方介質(zhì)中薄板的響應(yīng)，并且這種屏蔽效果不會(huì)隨著關(guān)斷時(shí)間的增加而減少。

5 結(jié) 論

1)積分方程法適合模擬層狀介質(zhì)中薄板的響應(yīng)，計(jì)算速度快，計(jì)算結(jié)果準(zhǔn)確，為反演奠定了扎實(shí)的基礎(chǔ)。

2)無(wú)論是高阻覆蓋層或者低阻覆蓋層對(duì)于下方基巖中薄板的響應(yīng)值影響很小，幾乎可以忽略。

3)基巖對(duì)于上方覆蓋層中薄板的響應(yīng)值有著顯著的影響。高阻基巖會(huì)增大薄板的響應(yīng)值，且增幅隨著關(guān)斷時(shí)間的增加而增加；低阻基巖會(huì)抑制薄板的響應(yīng)值，這種影響不會(huì)隨著關(guān)斷時(shí)間的變化而變化。

致 謝 感謝湖南五維地質(zhì)科技有限公司提供的幫助。

[1]牛之鏈.時(shí)間域電磁法原理[M].長(zhǎng)沙:中南大學(xué)出版社,2007.

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

[3]薛國(guó)強(qiáng),李貅,底青云.瞬變電磁法正反演問(wèn)題研究進(jìn)展[J].地球物理學(xué)進(jìn)展,2008,23(4):1 165-1 172.

[4]高陽(yáng),陳志軍,熊華山.矩形框瞬變電磁法一維正演研究[J].工程地球物理學(xué)報(bào),2015,12(6):726-731.

[5]李建慧,朱自強(qiáng),曾思紅,等.瞬變電磁法正演計(jì)算進(jìn)展[J].地球物理學(xué)進(jìn)展,2012,27(4):1 393-1 400.

[6]米薩克 N.納比吉安.勘探地球物理電磁法[M].趙經(jīng)祥,王艷君譯.北京:地質(zhì)出版社,1992.

[7]M.I.Epov, E.P.Shurina, O.V.Nechaev. 3D forward modeling of vector field for induction logging problems[J]. Russian Geology and Geophysics,2007,48(9):770-774.

[8]殷長(zhǎng)春,張博,劉云鶴,等.2.5維起伏地表?xiàng)l件下時(shí)間域航空電磁正演模擬[J].地球物理學(xué)報(bào),2015,58(4):1 411-1 424.

[9]Michael Commer, Gregory Newman. A parallel finite-difference approach for 3D transient electromagnetic modeling with galvanicsources[J]. Geophysics,2004,69(5):1 192-1 202.

[10]孫懷鳳,李貅,李術(shù)才,等.考慮關(guān)斷時(shí)間的回線源激發(fā)TEM三維時(shí)域有限差分正演[J].地球物理學(xué)報(bào),2013,56(3):1 049-1 064.

[11]徐正玉,楊海燕,鄧居智,等.回線源三維地—井瞬變電磁法FDTD數(shù)值模擬[J].工程地球物理學(xué)報(bào),2015,12(3):327-332.

[12]Wiliam A Sanfilipo, Gerald W Hohmann. Integral equation solution for the transient electromagnetic response of a three dimensional body in a conductive half-space[J]. Geophysics, 1985,50(5):798-809.

[13] Gregory A Newman, Gerald W Hohmann. Transient electromagnetic responses of high-contrast prisms in a layered earth[J]. Geophysics,1988,53(5):691-706.

[14]魏寶君，王甜甜，王穎.用磁流源并矢Green函數(shù)的遞推矩陣方法計(jì)算層狀各向異性地層中多分量感應(yīng)測(cè)井響應(yīng).地球物理學(xué)報(bào)[J],2009,52(11):2 920-2 928.

[15]陳桂波.各向異性地層中電磁場(chǎng)三維數(shù)值模擬的積分方程算法及其應(yīng)用[D].長(zhǎng)春:吉林大學(xué),2009.

[16]胡俊華.瞬變電磁積分方程法正演模擬研究[D].武漢:中國(guó)地質(zhì)大學(xué),2014.

[17]賴劉保.接地線源近區(qū)三維薄板瞬變電磁響應(yīng)積分方程法數(shù)值模擬[D].長(zhǎng)沙:中南大學(xué),2014.

The Simulation of Transient Response of Thin Plate in Layered Media by Integral Equation Method Modeling

Liu Yuanyuan1,2,Xi Zhenzhu1,2,Jiang Huan1,2,Wang He1,2

(1.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,ChangshaHunan410083,China;2.MarineMarineExplorationandEquationresearchInstitute,CentralSouthUniversity,ChangshaHunan410083,China)

Transient electromagnetic method multidimensional forward and inversion technology development is slow, because its forward modeling method is not mature, the calculated speed is slow and the accuracy is low, which can not be up to the requirements of inversion. This paper clarifies the definition of sheet, uses integral equation method to simulate the transient response of sheets in layered medium, simplifies meshing and controls the size of matrix to take direct method to solve coefficient matrix. The reliability of this algorithm and this program is verified by the comparison of analytical solution of the homogeneous half space, and the accuracy of the algorithm is demonstrated by the error analysis. At the same time, the effect of the overburden and bedrock on the response of the thin plate is studied. It is found that the response of the high resistance and low resistance to thin plate is very small, and that high resistance bedrock would expand the response of thin plate, and that low resistivity bedrock would shield the response of thin plate. This study lays the foundation for 2-dimensional inversion of transient electromagnetic, and also provides a theoretical reference for the production of thin plate in the future.

transient electromagnetic method; integral equation; layered media; thin plate

1672—7940(2016)02—0143—06

10.3969/j.issn.1672-7940.2016.02.001

國(guó)家自然科學(xué)基金項(xiàng)目(編號(hào):41304090)；深圳市戰(zhàn)略新興產(chǎn)業(yè)發(fā)展專項(xiàng)資金項(xiàng)目(編號(hào)：CXZZ20120618165608947)

劉愿愿(1991-)，男，碩士研究生，主要從事于電磁法正反演研究。E-mail:csu_lyy@163.com

席振銖(1966-)，男，教授，博士生導(dǎo)師，主要從事瞬變電磁方法與技術(shù)研究。E-mail:xizhenzhu@163.com

P631.3

A

2015-10-09