起伏地形下可控源音頻大地電磁三維數(shù)值模擬

鐘蘇美，林昌洪，謝裕春

(1.上?？睖y設(shè)計研究院有限公司，上海 200434； 2.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院，北京 100083；3.北京城建勘測設(shè)計研究院有限責(zé)任公司，北京 100101)

0 引 言

加拿大學(xué)者Goldstein和Strangway[1]在20世紀(jì)70年代提出可控源音頻大地電磁法(CSAMT)。近40多年來，CSAMT得到了快速的發(fā)展，目前已廣泛地應(yīng)用于礦產(chǎn)、油氣和地?zé)岬荣Y源探測以及水文環(huán)境和地質(zhì)工程等領(lǐng)域[2-17]。然而，大多數(shù)的CSAMT工作是在存在地形起伏的區(qū)域開展的，很少在水平地表條件下采集CSAMT數(shù)據(jù)。前人研究已表明：地形對CSAMT數(shù)據(jù)存在強烈的影響，因此在CSAMT數(shù)據(jù)反演處理解釋時需要考慮地形。

近年來，隨著計算機技術(shù)的進(jìn)步，CSAMT二維正反演算法的研究逐漸趨于成熟，并且很多二維正反演算法也開始考慮地形的影響。例如，Mitsuhata[3]采用有限單元法實現(xiàn)了考慮地形的CSAMT二維正演。李予國等[2]采用有限單元法實現(xiàn)了考慮海底地形起伏的海洋可控源電磁法(CSEM)的二維數(shù)值模擬。雷達(dá)等[13]采用有限元法實現(xiàn)了起伏地形下CSAMT二維正反演研究。然而，使用二維正反演方法的前提是假設(shè)地下電阻率結(jié)構(gòu)是二維分布的，并且測線與二維電阻率結(jié)構(gòu)走向垂直。但是，實際地下地質(zhì)情況和地表的地形情況通常比較復(fù)雜，地質(zhì)結(jié)構(gòu)和地形大部分情況下都是三維的。前人對水平地表的大地電磁數(shù)據(jù)的研究表明，采用二維反演技術(shù)處理三維異常體產(chǎn)生的數(shù)據(jù)，很可能得到不可靠的地電模型，有時甚至得到畸變的結(jié)果。采用二維地形處理技術(shù)去逼近三維地形，也會導(dǎo)致類似的錯誤解釋結(jié)果。因而，需要有考慮三維地形起伏條件的三維CSAMT正、反演算法用于CSAMT數(shù)據(jù)的處理與解釋，以此能夠得到更加可靠的數(shù)據(jù)解釋結(jié)果。

開展帶地形三維CSAMT反演研究的前提是：先研究帶地形三維CSAMT正演模擬算法。Zhdanov等[5]2006年采用積分方程方法實現(xiàn)了考慮地形的三維海洋CSEM數(shù)值模擬。Sasaki[6]2011年采用有限差分法實現(xiàn)了考慮地形的三維海洋CSEM數(shù)值模擬。Schwarzbach等[10]2011年采用有限單元法實現(xiàn)了考慮地形的三維海洋CSEM數(shù)值模擬。然而，CSEM數(shù)值模擬中關(guān)于海底地形處理，與CSAMT數(shù)值模擬中關(guān)于地表地形的處理存在一些差別：海底地形分界面兩側(cè)的海水(其電阻率約為0.3 Ω·m)與海底介質(zhì)(約1 Ω·m)的電阻率差異不明顯，而地表地形分界面兩側(cè)的空氣(無窮大)與地下介質(zhì)(小于105Ω·m)的電阻率差異非常大。因此，CSAMT數(shù)值模擬中處理地形，需要更加謹(jǐn)慎。前人關(guān)于水平地表CSAMT三維數(shù)值模擬的研究較多[11-17],但目前還沒有考慮三維地表地形的CSAMT三維數(shù)值模擬報道。

本文介紹一種采用有限差分法模擬三維起伏地形的三維CSAMT響應(yīng)算法。與Sasaki[6]采用的有限差分法模擬海底地形的三維海洋CSEM響應(yīng)算法不同的是：(1)本文模擬三維地表地形；(2)本文采用不同的方式處理場源奇異性。由于場源附近的電場和磁場變化非常迅速，如何處理場源奇異性的策略通常有兩種：(1)采取將電磁場總場分離成背景場和二次場的策略；(2)采用偽δ函數(shù)來代替場源直接計算總場的策略。正如Mitsuhata[3]在考慮地形的CSAMT二維數(shù)值模擬的研究中所述：當(dāng)?shù)叵陆橘|(zhì)復(fù)雜的情況下，很難找到一個簡單的背景場構(gòu)造，而采用偽δ函數(shù)來代替場源直接計算總場的策略不存在這個問題，因此更適合于模擬復(fù)雜地質(zhì)結(jié)構(gòu)。本文中，首先將簡單介紹采用有限差分法如何模擬考慮地形的三維CSAMT響應(yīng)， 然后給出理論模型的模擬結(jié)果，最后對三維數(shù)值模擬算法的正確性進(jìn)行檢驗。

1 有限差分法模擬起伏地形的三維CSAMT響應(yīng)

本算法是在林昌洪等[11]提出的有限差分法求解水平地表條件下三維CSAMT 響應(yīng)算法的基礎(chǔ)上修改代碼實現(xiàn)的，因此下面將簡單介紹有限差分算法，重點介紹與水平地表三維CSAMT 數(shù)值模擬算法[11]的不同之處。

1.1 電磁場耦合方程

可控源音頻大地電磁法工作頻率范圍為0.1～10 000 Hz，在此頻段可忽略位移電流的作用，地下介質(zhì)和空氣的磁導(dǎo)率相接近。在電偶極子或者水平線圈激發(fā)情況下,總電場強度矢量E和總磁場強度矢量H相互作用的關(guān)系滿足麥克斯韋方程：

∮Hdl=?(σE+Js)ds

(1)

∮Edl=?iωμ0Hds

(2)

式中：ω表示角頻率；μ0為真空磁導(dǎo)率；σ為介質(zhì)的電導(dǎo)率；E表示電場強度；H表示磁場強度；Js為源電流密度。

1.2 場源奇異性處理

在水平地表三維CSAMT 數(shù)值模擬算法[11]中，也是采用將總場分離成背景場和二次場，最后合成總場的策略處理場源奇異性。由于直接模擬總場的策略更適合于模擬復(fù)雜地質(zhì)結(jié)構(gòu)，參考Mitsuhata[3]提出的CSAMT二維數(shù)值模擬的策略，在本文的考慮地形三維CSAMT數(shù)值模擬中，使用偽δ函數(shù)代替場源直接模擬總場以解決場源的奇異性，偽δ函數(shù)表達(dá)式如下：

δs(x-x0)=

(3)

其中：x0是源的橫坐標(biāo)；τ是常數(shù)，控制等效場源作用的范圍和幅值。

當(dāng)源和研究區(qū)域都是三維時，偽δ函數(shù)的表達(dá)式為:

δ(z-z0)

(4)

式中：r0為場源的矢量坐標(biāo)；x0、y0、z0分別表示場源X軸坐標(biāo)、Y軸坐標(biāo)和Z軸坐標(biāo)。

根據(jù)偽δ函數(shù)的表達(dá)式，式(1)中的源電流密度Js可寫成：

(5)

式中：I為場源電流強度；L為場源長度。

1.3 三維正演方程

在笛卡兒右手坐標(biāo)系中，用交錯采樣剖分網(wǎng)格[18]在研究區(qū)域內(nèi)對電阻率和電磁場總場滿足的麥克斯韋積分方程(1)和(2)進(jìn)行離散化，可獲得關(guān)于地下各網(wǎng)格單元采樣點處總磁場的正演方程：

AU=b

(6)

其中：A為n×n階對稱大型稀疏系數(shù)矩陣；U為待求解總磁場分量構(gòu)成的n×1維列向量；b為一個與邊界場值和場源有關(guān)的n×1維列向量；n為待求解的未知磁場分量的個數(shù)。采用預(yù)條件的雙共軛梯度穩(wěn)定算法[18]，求解該線性方程組，從而獲得總磁場值列向量。

1.4 邊界條件

在水平地表三維CSAMT 數(shù)值模擬算法[10]中，由于求解的是二次磁場，采取研究區(qū)域的空中頂邊界、地下底邊界和4個側(cè)邊界處的二次磁場值為零的方法處理邊界條件。當(dāng)選擇合適的背景電阻率求解一次場時，二次場基本由異常體的剩余電阻率(電阻率-背景電阻率)產(chǎn)生，而且研究區(qū)域的空中頂邊界、地下底邊界和4個側(cè)邊界處離異常體較遠(yuǎn)，因此求解二次磁場時可以采用邊界二次磁場值為零的邊界條件。但是，在本文中是直接計算總場，總場在上述邊界處的磁場值不滿足為零的條件，因此，不能采用相同的策略。

本文中，筆者仍然采用第一類邊界條件。在有限差分?jǐn)?shù)值模擬中，常常把區(qū)域剖分得足夠大，讓異常體和場源距離邊界足夠遠(yuǎn)，邊界上的總電磁場基本滿足遠(yuǎn)區(qū)條件。因此，可以認(rèn)為在邊界上的總磁場是均勻半空間中電偶極子在遠(yuǎn)區(qū)產(chǎn)生的一次磁場響應(yīng)。參考Ward和Hohmann[9]的分析，地表水平電偶極子在均勻半空間中產(chǎn)生一次磁場表達(dá)式，加上準(zhǔn)靜態(tài)和遠(yuǎn)區(qū)場條件，可以將原一次磁場表達(dá)式進(jìn)一步簡化，得到X方向電偶極子在地表、接收在空中的遠(yuǎn)區(qū)一次磁場響應(yīng)的計算表達(dá)式如下：

(7)

(8)

(9)

X方向電偶極子在地表、接收在地下的遠(yuǎn)區(qū)一次磁場響應(yīng)的計算表達(dá)式如下：

(10)

(11)

(12)

Y方向電偶極子在地表、接收在空中的遠(yuǎn)區(qū)一次磁場響應(yīng)的計算表達(dá)式如下：

(13)

(14)

(15)

Y方向電偶極子在地表、接收在地下的遠(yuǎn)區(qū)一次磁場響應(yīng)的計算表達(dá)式如下：

(16)

(17)

(18)

圖1 地表處用于計算地表總電場的差分關(guān)系Fig.1 View of the difference relationship used to calculate the total electric fields on the surface

1.5 地形處理

求解方程(6)可得到地下交錯采樣點上的總磁場值U。得到總磁場值U后，需要計算地表(水平或是起伏地形)處的總電場和總磁場值，進(jìn)而利用林昌洪等[11]提出的地表總電場和總磁場計算表視電阻率和相位的表達(dá)式計算得到地表處的視電阻率和相位響應(yīng)。

不考慮磁導(dǎo)率參數(shù)，磁場在空氣-地下介質(zhì)分界面處的法向磁場和切向磁場都是連續(xù)的，因此，無論水平地表的磁場還是起伏地形分界面處的磁場，其求解表達(dá)式是一樣的?？梢岳貌蓸狱c處的總磁場值，通過和水平地表同樣的線性插值公式計算得到地表網(wǎng)格頂面中心處的總磁場的三個分量。

對于水平地表的電場，通過麥克斯韋方程▽×E=-iμ0ωH，地表網(wǎng)格頂面中心處的總電場的x分量(圖1)可以通過下列表達(dá)式計算：

(19)

ρ(i,j,k)

(20)

(21)

但是，對于起伏地形分界面處總電場的x分量不能采用上述公式計算。例如，網(wǎng)格塊 (i,j,k)在地下，而(i-1,j,k)在空氣中時，靠近分界面處的法向電流密度變化梯度大，公式(20)采用界面處法向電流密度進(jìn)行插值，計算結(jié)果不夠精確。為了得到更為精確的結(jié)果，可以采用分界面右側(cè)2個網(wǎng)格處的電流密度(圖1)進(jìn)行線性插值，從而公式(20)變?yōu)椋?/p>

Jx2-Js]×ρ(i,j,k)

(22)

如Wannamaker等[4]所述，為了避免和空氣-地下分界面有關(guān)的場值偏導(dǎo)數(shù)的不連續(xù)性，應(yīng)該選擇計算偏導(dǎo)數(shù)的2個場值的節(jié)點都位于同一介質(zhì)中。網(wǎng)格塊 (i,j,k)在地下而(i-1,j,k)在空氣中時，計算公式(21)時也應(yīng)采用位于同一介質(zhì)內(nèi)2個節(jié)點的垂直電場進(jìn)行計算場值的偏導(dǎo)數(shù)，避免計算結(jié)果的不精確。

當(dāng)網(wǎng)格塊 (i,j,k)在地下，而(i+1,j,k)在空氣中時，公式(20)和(21)可以做類似的修改。對于網(wǎng)格塊 (i,j,k)在地下，而(i+1,j,k)和(i-1,j,k)都在空氣中的情況，在網(wǎng)格剖分的時候應(yīng)避免類似的剖分。

2 算法驗證

為了檢驗算法的有效性和正確性，筆者設(shè)計一個水平地表下的三維異常體模型，將本文算法的計算結(jié)果和水平地表CSAMT三維正演算法[11]的計算結(jié)果進(jìn)行對比。同時為了檢驗本文開發(fā)的帶地形CSAMT三維正演算法在計算三維起伏地形界面處響應(yīng)的正確性，筆者設(shè)計一個三維山峰純地形模型，并且把場源放置在距離三維地形較遠(yuǎn)處，以保證測區(qū)處于遠(yuǎn)區(qū)，將本文算法的計算結(jié)果和三維有限元大地電磁數(shù)值模擬算法[7]的計算結(jié)果進(jìn)行對比。

2.1 水平地表三維異常體模型

X軸方向場源，場源中心坐標(biāo)(-2.4 km,0 km,0 km)；測線在X軸上，在電阻率值為100 Ω·m的均勻半空間中埋藏一個長、寬、高為240 m×120 m×120 m，頂面埋深為100 m，電阻率為10 Ω·m低阻長方體。長方體的中心距離發(fā)射場源2.4 km，計算頻率為1 000 Hz。模型如圖2所示。

分別采用本文算法和水平地表CSAMT三維正演算法[11]利用32位 Windows XP系統(tǒng)臺式計算機對圖2模型進(jìn)行數(shù)值模擬計算，臺式計算機的配置為Intel TM 2.83 GHz 四核處理器、4G內(nèi)存。2種算法都采用56×58×38(空氣層：10)的網(wǎng)格，計算頻率為1 000 Hz。圖3為頻率1 000 Hz時，水平地表三維異常體模型(圖2)的本算法(黑色圓圈)和林昌洪等[11](黑色實線)模擬結(jié)果的對比圖。其中，視電阻率響應(yīng)的最大誤差0.23%，相位響應(yīng)的最大誤差0.28°。

圖2 水平地表三維異常體模型(a)異常體模型三維視圖；(b) 異常體模型X-Z平面視圖 Fig.2 Sketches of the horizontal surface three-dimensional anomalous body model

2.2 三維山峰模型

本文采用Nam等[7]文中的三維山峰純地形模型，該模型與Wannamaker 等[4]提出的二維經(jīng)典地形模型類似。模型如圖4所示，山峰高度設(shè)置為0.45 km(Z軸：-0.45 km和0 km之間)，山頂寬度為0.45 km (X、Y軸：-0.225 km 和 0.225 km之間),山的底部寬度為2 km(X、Y軸：-1 km 和 1 km 之間)，地下介質(zhì)為電阻率為100 Ω·m的均勻半空間。將坐標(biāo)原點設(shè)于山頂中心，為了讓山峰地形處于遠(yuǎn)區(qū)，將一個X軸電偶極子源放置在離山頂中心65 km處(-65 km，0 km，0 km)。

圖3 頻率為1 000 Hz時水平地表三維異常體模型的本文算法(黑色圓圈)和林昌洪等[11]算法的模擬結(jié)果的對比(a)視電阻率；(b)相位響應(yīng)；(c)視電阻率數(shù)值模擬結(jié)果誤差；(d)相位響應(yīng)數(shù)值模擬結(jié)果誤差Fig.3 Results of 3D responses at frequency of 1,000 Hz calculated by the forward scheme (black circle) described in this paper(the horizontal surface three-dimensional anomalous body with a X-direction source located at -2.4 km,0 km,0 km) compared with responses calculated by the method given by Lin C H et al[11](black solid line,along the line x=0 m)

圖4 三維山峰地形模型(a) 山峰地形模型X-Y平面圖；(b) 山峰地形模型X-Z平面圖 Fig.4 Sketches of the three-dimensional trapezoidal-hill models

圖5為頻率2 Hz 時，本文算法計算的三維山峰模型(圖4)的XY和YX模式三維CSAMT視電阻率(圖5(a)和(b))和相位(圖5 (c)和(d))響應(yīng)(黑色虛線)與 Nam等[7]算法計算的大地電磁響應(yīng)(黑色圓圈)的對比圖。由于測區(qū)位于遠(yuǎn)區(qū)，場源影響小，三維CSAMT結(jié)果與大地電磁響應(yīng)基本一致。圖5的結(jié)果顯示本文有限差分的計算結(jié)果和有限元的計算結(jié)果吻合得很好，其中，XY和YX模式視電阻率響應(yīng)的最大誤差分別為4.03%、1.335%；相位響應(yīng)的最大誤差分別為0.9°、0.45°，驗證了本文提出的帶地形CSAMT三維數(shù)值模擬算法的準(zhǔn)確性。

圖5 頻率2 Hz 時XY和YX模式三維山峰模型(圖4)的本算法響應(yīng)(黑色圓圈)和 Nam等[7]算法計算的大地電磁響應(yīng)(黑色實線)模擬結(jié)果的對比(a)、(b) XY、YX模式視電阻率結(jié)果對比；(c)、(d) XY、YX模式兩種數(shù)值模擬視電阻率結(jié)果誤差；(e)、(f) XY、YX模式相位結(jié)果對比；(g)、(h)XY、YX模式兩種數(shù)值模擬相位結(jié)果誤差Fig.5 Comparison of 3D responses (at frequency of 2 Hz) of the method given by this paper (black circle) and 3D FEM MT responses of the method given by Nam et al[7] (black solid line) for models of XY and YX

3 結(jié) 論

本文在水平地表三維有限差分CSAMT數(shù)值模擬算法的基礎(chǔ)上，實現(xiàn)了模擬三維起伏地形CSAMT響應(yīng)的正演算法。在本算法中，提出了新的由地下交錯采樣點處的總磁場計算起伏地形下空氣-地下介質(zhì)分界面處電場的表達(dá)式；同時，為了避免了原有將總場分離成背景場和二次場的策略在復(fù)雜地質(zhì)條件下難以選擇合適背景電阻率構(gòu)造的問題，引入偽δ函數(shù)來代替麥克斯韋方程中的場源項，在三維數(shù)值模擬中直接計算總電磁場，從而解決場源奇異性問題；為了三維正演方程可以直接求解總電磁場，將原有的研究邊界處二次場值為零的邊界條件修改為均勻半空間中電偶極子在遠(yuǎn)區(qū)的一次場值的邊界條件，并通過公式推導(dǎo)，獲得計算均勻半空間中電偶極子在遠(yuǎn)區(qū)的一次場值的解析表達(dá)式。

為了驗證本文算法的可靠性，對水平地表三維異常體模型和三維山峰純地形模型進(jìn)行了三維數(shù)值模擬。采用本文算法模擬得到兩個三維模型的視電阻率響應(yīng)和相位響應(yīng)，計算結(jié)果與前人算法的結(jié)果吻合得很好，視電阻率響應(yīng)的最大偏差小于4.03%，相位響應(yīng)的最大偏差小于0.9°，驗證了本文所實現(xiàn)的起伏地形下三維CSAMT數(shù)值模擬算法的正確性。

參考文獻(xiàn)：

[1] GOLDSTEIN M A,STRANGWAY D W.Audio-frequency magetotellurics with a grounded electric dipole source[J].Geophysics,1975,40(4)：669-683.

[2] LI Y,KEY K.2D marine controlled-source electromagnetic modeling: Part 1—An adaptive finite-element algorithm[J].Geophysics,2007,72(2)：51-62.

[3] MITSUHATA Y.2-D electromagnetic modeling by finite-element method with a dipole source and topography[J].Geophysics,2000,65(2)：465-475.

[4] WANNAMAKER P E,STODT J A,RIJO L.Two-dimensional topography responses in magnetotelluric modeled using finite elements[J].Geophysics,1986,51(11)：2131-2144.

[5] ZHDANOV M S,LEE S K,YOSHIOKA K.Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity[J].Geophysics,2006,71(6)：G333-G345.

[6] SASAKI Y.Bathymetric effects corrections in marine CSEM data[J].Geophysics,2011,76(3)： F139-F146.

[7] NAM M J,KIM H J,SONG Y,et al.3D magnetotelluric modelling including surface topography[J].Geophysical Prospecting,2007,55：277-287.

[8] SCHWALENBERG K，EDWARDS R N.The effect of seafloor topography on magnetotelluric fields: an analytical formulation confirmed with numerical results[J].Geophysical Journal International,2004,159：607-621.

[9] WARD S H，HOHMANN G W.Geophysical Electromagnetic Theory,Electromagnetic Method in Applied Geophysics[M].Katris: Society of Exploration Geophysicists,1988:1.

[10] SCHWARZBACH C,BORNER R U,SPITZER K.Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics: a marine CSEM example[J].Geophysical Journal International,2011,187：63-74.

[11] 林昌洪，譚捍東，舒晴，等.可控源音頻大地電磁法三維共軛梯度反演研究[J].地球物理學(xué)報，2012，55(11): 3829-3838.

[12] 鄧居智，譚捍東，陳輝，等.CSAMT三維交錯采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)進(jìn)展，2011，26(6): 2026-2032.

[13] 雷達(dá).起伏地形下CSAMT正反演研究與應(yīng)用[J].地球物理學(xué)報，2010，53(4):982-993.

[14] 張繼鋒，湯井田，喻言，等.基于電場矢量波動方程的三維可控源電磁法有限單元數(shù)值模擬[J].地球物理學(xué)報，2009，52(12): 3132-3141.

[15] 李帝銓,王光杰，底青云，等.基于遺傳算法的CSAMT最小構(gòu)造反演[J].地球物理學(xué)報，2008，51(4)：1234-1245.

[16] 陳桂波，汪宏年，姚敬金，等.用積分方程法模擬各向異性中三維電性異常體的電磁響應(yīng)[J].地球物理學(xué)報，2009，52(8)：2174-2181.

[17] 王若，王妙月，盧元林.三維三分量CSAMT法有限元正演模擬研究初探[J].地球物理學(xué)進(jìn)展，2007，22(2)：579-585.

[18] 譚捍東，余欽范，BOOKER John，等.大地電磁法三維交錯采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)進(jìn)展，2003，46(5)：705-711.