CSAMT 2.5維有限元數(shù)值模擬

姚大為， 王大勇， 王　剛， 朱　威， 李永博， 張振宇

(中國地質(zhì)科學院　地球物理地球化學勘查研究所, 國家現(xiàn)代地質(zhì)勘查工程技術(shù)研究中心，廊坊　065000)

?

CSAMT 2.5維有限元數(shù)值模擬

姚大為， 王大勇， 王剛， 朱威， 李永博， 張振宇

(中國地質(zhì)科學院地球物理地球化學勘查研究所, 國家現(xiàn)代地質(zhì)勘查工程技術(shù)研究中心，廊坊065000)

摘要：對三維場源二維地電模型的正演計算稱為2.5維數(shù)值模擬。這里從麥克斯韋方程組出發(fā)，①分別求解電磁場的一次場和二次場，將三維場源降為二維；②利用傅氏變換將空間域方程轉(zhuǎn)化為波數(shù)域，應用有限元求解波數(shù)域電磁場方程，引入無限元解決無窮遠邊界的收斂問題；③根據(jù)電磁場值實虛部的曲線特征，按對數(shù)等間隔選取21個波數(shù)。編寫代碼計算均勻半空間與解析解結(jié)果對比，電阻率的均方相對誤差均小于0.5%，證明其有效性。計算了三種地電模型的電磁場響應，對單一高、低阻體和高低阻組合體的模擬效果真實，異常中心位置基本吻合。結(jié)果證明,這里的CSAMT2.5-D正演算法可以模擬較為復雜的地電模型，并取得良好的效果。

關(guān)鍵詞：CSAMT； 2.5-D； 有限元； 無限元

0引言

可控源音頻大地電磁法(CSAMT)，是頻率域電磁測深方法的之一。它是通過人工發(fā)射不同頻率可控的電磁信號，以達到探測地下地質(zhì)體的目的，提高信噪比，從而成為物探工作應用廣泛且高效的電磁測深法[1-4]。當前CSAMT二維正反演以小的計算量和優(yōu)越的計算速度而被廣泛研究[5-6]，三維正反演可以真實地模擬地質(zhì)構(gòu)造的電性分布，但由于當前商用計算機的限制，三維正反演算法計算速度和精度二者相互制約，提高其中之一，必以降低另一者為代價。目前，隨著計算機技術(shù)的提高，三維正反演算法處理解釋實際數(shù)據(jù)已經(jīng)開始應用。然而，實際應用中的地電結(jié)構(gòu)一般是只沿傾向發(fā)生變化，在走向方向上變化很小的二維模型，場源卻是三維的。這種地電結(jié)構(gòu)是二維的，場源是三維的CSAMT觀測，可用有限元2.5維方法進行數(shù)值模擬。

針對該問題，國內(nèi)、外學者做了大量工作，COGGON[7]提出2.5 維有限元電磁模擬方法，開啟了有限元法在電磁法中的應用;UNSWORTH 等[8]通過求解二次場和引入無限元，實現(xiàn)了CSAMT 2.5D 正演計算，并討論有限長度電性源激發(fā)時場的特征； LU 等[9]采用二次場算法進行CSAMT 2.5D 正演模擬；MITSUHATA 等[10-11]用總場算法進行了CSAMT2.5-D正演模擬，采用偽δ函數(shù)代替了源的影響，利用等參有限元技術(shù)解決了地形問題；LI 等[12]應用有限元進行了2.5D自適應非結(jié)構(gòu)化的三角形網(wǎng)格海洋電磁有限元正演計算，同時在解線性方程時運用了預置條件共軛梯度法，再用準最小殘差法求取其他分量；孟永良等[4]在實現(xiàn)有限元法CSAMT2.5D正演時，采用了總場算法；陳金窗等[13]研究了并行計算環(huán)境下2.5D CSAMT正演模擬；陳小斌等[14]為消除矩形單元中的節(jié)點，而采用了線源有限元直接迭代算法，從而減小了線性方程組的階數(shù)；底青云等[15-17]做了CSAMT2.5 D有限元正演，計算了不同模型的電磁響應；沈金松等[18-19]在海洋可控源電磁法模擬中研究2.5 維電磁響應與波數(shù)的關(guān)系。上述表明，電磁場2.5D 有限元正演可以平衡計算速度和精度的問題，對復雜的物性分布，有較好地電磁響應。作者在求解頻率域有源麥克斯韋方程組時，①將電磁場分為一次場和二次場分別求解，從而將三維場源降為二維；②研究區(qū)邊界處引入無限單元；③波數(shù)的選取則依據(jù)電磁場值實虛部的曲線特征，按對數(shù)等間隔選取，盡量減少波數(shù)個數(shù)的選取，以便減小相應的計算量；④模擬均勻半空間驗證代碼有效性；⑤設(shè)計三種地電模型進行了正演試算，取得較好的效果。

1電磁場方程

頻率域麥克斯韋方程：

(1)

(2)

(3)

(4)

對式(2)兩邊取旋度，結(jié)合狀態(tài)方程和邊界條件，可得到電場矢量的非齊次亥姆霍茲方程表示如式(5)所示。

(5)

忽略位移電流，麥克斯韋方程電場和磁場的空間和時間對偶關(guān)系表述如式(6)、式(7)所示。

(6)

(7)

取式(6)的旋度，并代入式(7)，電場矢量的擴散方程可寫成式(8)。

(8)

同理，磁場矢量擴散方程表示為式(9)。

(9)

利用矢量恒等式：

(10)

(11)

在均勻大地不存在異常體的介質(zhì)中，一次場在全區(qū)滿足式(12)和式(13)。

(12)

(13)

σ1為電導率得背景值，一次場由解析法直接求解。從式(10)中減去式(12)，從式(11)中減去式(13)，得到二次場方程，如式(14)和式(15)所示。

(14)

(15)

(16)

同法，消去電偶源項，得到二次磁場微分方程如式(17)所示。

(17)

取時間因子為eiωt，對式(17)進行傅立葉變換，得到頻率域總場、二次場的微分方程[19]：

(18)

(19)

(20)

(21)

(22)

(23)

(24)

2波數(shù)域有限單元方程

設(shè)有限單元法的平面算域為Ω，其外圍邊界為?Ω，利用Galerkin加權(quán)余量導出控制式(23)、式(24)的有限元方程，則波數(shù)域電場和磁場的殘差為：

(25)

(26)

Galerkin有限單元法的加權(quán)余量應滿足如下積分方程：

(27)

式中Ni表示第i個節(jié)點的插值函數(shù)。利用格林定理：

(28)

將式(27)轉(zhuǎn)化為如下形式：

(29)

同法可將磁場分量表述如式(30)所示。

(30)

式中nx和nz分別表示法線與坐標軸之間夾角α的余弦函數(shù)和正弦函數(shù)。運用格林式(28)對式(29)、式(30)逐項積分得：

(31)

(32)

3有限單元積分

研究區(qū)剖分成矩形單元[22-23]，矩形頂點定義為節(jié)點，采用雙線性插值進行等參單元分析。圖1(a)中的正方形單元為母單元，形函數(shù)如式(33)所示。

圖1　等參單元示意圖Fig．1　Isoparametric element schematic(a)母單元；(b)子單元

(33)

式中：ξi、ηi是點i的坐標(i=1,2,3,4)。形函數(shù)需滿足式(34)。

(34)

將其展開，得到式(35)。

Ni=a1ξη+a2ξ+a3η+a4

(35)

式中：a1、a2、a3、a4為常數(shù)。

圖1(b)中子單元頂點的坐標分別為：(x1,y1)、(x2,y2)、(x3,y3)、(x4,y4)，函數(shù)值為：u1、u2、u3、u4。子單元上的u、x 和z關(guān)系如式(36)。

當源為水平電偶極子時，式(31)和式(32)可簡化處理為如下形式的離散化方程：

(37)

(38)

求解此方程，要進行坐標轉(zhuǎn)換，坐標轉(zhuǎn)換過程參考文獻[21]。

每個單元內(nèi)4個節(jié)點的合成的電磁場方程如式(39)所示。

(39)

矩陣形式為：

(40)

4無限元的應用

無限元是一種在物理上進行界定的、幾何意義無限大的有限單元[24-26]。它解決了無窮遠邊界問題，改善了有限元的局限性。二者的結(jié)合使用，首先真正意義上處理了外邊界條件，保證了無窮遠處電場場值衰減為零；其次它縮減了有限元的計算區(qū)域，可提高計算速度。

應用Astley 無限單元，在有限元網(wǎng)格的邊界至無限遠方向，通過無限元坐標映射將整體坐標映射到局部坐標上，實現(xiàn)無限遠處的積分，無需考慮邊界條件。位場的插值函數(shù)稱為形函數(shù)，坐標映射函數(shù)叫做映射函數(shù)，二者有所不同。

圖2　一維無限元映射Fig．2　1-D infinite element mapping

圖3　二維無限元映射Fig．3　2-D infinite element mapping

x=G1(ξ)x1+G2(ξ)x2

(41)

方程中：G1(ξ)和G2(ξ)分別為：

(42)

這種映射適用于多維且具唯一性。

圖3為基于矩形的二維無限元幾何映射，原點O是無限單元映射的出發(fā)點。幾何形式的無限單元是以O(shè)為起點，沿同一方向的兩條射線，與有限矩形單元的交點為:①、②。在射線上再選取兩個點③和④，使其滿足①為O和③的中點，②為O和④的中點。無限遠處兩個節(jié)點的場值為零，不參與運算。圖3(b)中的矩形即是無限單元映射到有限元后的效果，從而可進行積分運算。在二維映射中，ξ表示無窮遠方向，映射方式和一維無限元映射方式相同。2.5維正演，假定電磁場在y方向上不變，不做任何處理，只對z方向上做無限單元映射和函數(shù)插值，在ξ-η平面內(nèi)，采用的映射形式和形函數(shù)與有限單元完全相同。

坐標映射函數(shù)如式(43)所示。

(43)

式中，Ci為長度的坐標。

從圖3中可以看出：

(44)

四節(jié)點的無限單元節(jié)點坐標映射函數(shù)為式(45)。

(45)

其中：

(46)

合并式(45)、式(46)得式(47)。

(47)

無限單元的上述坐標映射，使得電磁場延伸至無窮遠處且衰減為零，因此不需要再進行邊界積分，從而不再考慮邊界條件。

(48)

此形函數(shù)使無限元系數(shù)矩陣具有對稱性，目的是使其與有限元剛度矩陣一致。在系數(shù)矩陣中添加單元編號，則稀疏對稱復系數(shù)矩陣和方程組表述如下：

(49)

(50)

(51)

(52)

5波數(shù)選取

(53)

圖4　均勻大地條件下不同ζ的隨τ的變化曲線Fig．4　Different curve along with the change of τ ζ in uniform earth conditions(a)實部變化曲線;(b)虛部變化曲線

6模型算例

6.1驗證算法有效性

對均勻半空間模型進行了2.5維正演和解析解結(jié)果對比。通過對比異常區(qū)的四條曲線及數(shù)據(jù)誤差，驗證了CSAM2.5維有限元正演算法的正確性和精度。2.5維正演剖分的網(wǎng)格節(jié)點如圖5所示，X方向上節(jié)點個數(shù)64個，Z方向節(jié)點個數(shù)55個，另外空氣層層數(shù)為11層，采用x方向電偶極子。接收位置沿x方向，范圍從200 m～5 000 m。對比了收發(fā)距4 000 m、4 400 m、4 600 m和5 000 m四個測點的電阻率結(jié)果，圖6為四個測點不同頻率電阻率的對比結(jié)果，紅色實線為一維正演計算得到的電阻率曲線，綠色圓點為2.5維正演計算出來的電阻率曲線。由圖6可以看出,兩種方法計算結(jié)果在過渡區(qū)和近區(qū)電阻率相對誤差較大，遠區(qū)相對誤差較小。經(jīng)過計算，電阻率的均方相對誤差均小于0.5 %。通過對比驗證，說明這里研究的2.5維正演算法的計算精度是有保障的,可以用于設(shè)計模型的正演計算。

圖5　網(wǎng)格剖分示意圖Fig．5　Grid section

圖6　不同收發(fā)距電阻率對比結(jié)果Fig．6　Contrast results of different receiving and transmitting distance(a)4 000 m處對比結(jié)果；(b)4 400 m處對比結(jié)果；(c)4 600 m處對比結(jié)果；(d)5 000 m處對比結(jié)果

6.2地電模型數(shù)值模擬

應用上述研究的CSAMT2.5-D有限元算法，對三種典型的二維地電模型進行數(shù)值模擬。繪制典型測點的測深曲線，選取異常體附近1 km范圍，繪制視電阻率擬斷面圖，分析異常體的位置和響應異常特征，說明算法對不同地電模型的響應規(guī)律。

6.2.1高阻異常體數(shù)值模擬結(jié)果

場源在原點(0，0)處，異常體大小為200 m×200 m，埋深為160 m，四個角點的坐標為(4 600，160)、(4 600，360)、(4 400，360)和(4 400，160)。背景電阻率為100 Ω·m，異常體電阻率為1 000 Ω·m(圖7(a))。

圖7　體模型示意圖Fig．7　Model diagram(a)高阻體; (b)低阻體; (c)高低阻組合體

圖8　高阻體模擬結(jié)果圖Fig．8　High resistivity simulation results(a)不同測點曲線; (b)不同頻率曲線; (c)模擬結(jié)果擬斷面圖

高阻異常體不同測點視電阻率測深曲線,如圖8(a)所示： 4 000 m和5 000 m 兩個測點在高頻段視電阻率均在100 Ω·m左右，4 000 m和5 000 m曲線分別在50 Hz和32 Hz進入過渡區(qū)。4 520 m響應異常值最高，4400 m和4 600 m的響應異常值相當。說明正演模擬對高阻異常體兩側(cè)邊界和中心位置反應靈敏，能夠較好控制異常體的位置。

高阻異常體不同頻率視電阻率曲線如圖8(b)所示：在選取的5個頻率中，隨著頻率的降低，高阻異常體的響應增強。高頻率值是異常體上方背景場的響應；隨頻率的降低異常體的響應占據(jù)主導，在異常體中心幅值最大，向兩側(cè)又趨于背景場的值。

高阻異常體數(shù)值模擬電阻率擬斷面圖如圖8(c)所示：在4 400 m～4 600 m之間有一個明顯的高阻異常，中心位置在4 520 m。高阻異常體上下方，等值線收斂較快，邊界確定清晰；左右兩側(cè)，等值線相對稀疏，左右邊界不容易清晰劃定。從圖8的分析可知，CSAMT2.5-D正演算法對高阻體的模擬效果很好，異常中心清晰，范圍基本能夠圈定，在異常體外圍不大的范圍內(nèi)，可以較快恢復到背景場。

6.2.2低阻異常體數(shù)值模擬結(jié)果

低阻異常體電阻率為10 Ω·m，其余各參與高阻異常體相同(圖7(b))。

低阻異常體不同測點視電阻率測深曲線圖如9(a)所示： 4 000 m和5 000 m 兩個測點在高頻段的視電阻率均在100 Ω·m左右，低頻段視電阻率成下降趨勢。4 520 m 響應異常值最低，4 400 m和4 600 m的響應異常值相對較高，二者的值相當。這說明了正演模擬對低阻異常體兩側(cè)邊界和中心位置反應比較靈敏，能夠較好控制異常體的位置；50 Hz以后各頻率相繼進入過渡區(qū)，低阻范圍向下延伸。低阻異常體不同頻率視電阻率曲線如9(b)所示：在選取的5個頻率中，隨著頻率的降低，低阻異常體的響應增強。高頻率值是異常體上方背景場的響應；隨頻率的降低異常體的響應占據(jù)主導，在異常體中心幅值最大，向兩側(cè)又趨于背景場的值。

低阻異常體數(shù)值模擬結(jié)果電阻率擬斷面圖如9(c)所示：在4 400 m～4 600 m之間有一個相對明顯的低阻異常，中心位置在4 520 m處。低阻異常體上方和兩側(cè)，等值線收斂較快，低阻異常體邊界相對清晰；低阻異常體下方，等值線沒有閉合，異常向下延伸，沒有恢復到背景場，這是50 Hz以后進入過渡區(qū)所致。從圖9中的分析可知，CSAMT2.5-D正演算法對低阻異常體的模擬效果稍差，中心位置能夠吻合，低頻段進入過渡區(qū)很難恢復背景場的響應異常特征。

6.2.3高低阻組合體數(shù)值模擬結(jié)果

高低阻組合體各參數(shù)和低阻高阻異常相同，只有異常體位置有所變化，組合體中：低阻異常體的四個角點坐標為(4 120，160)、(4 320，360)、(4 120，360)和(4 320，160)；高阻異常體的四個角點坐標為(4 520，160)、(4 720，360)、(4 520，360)和(4 720，160)，其余各參數(shù)見圖7(c)。

高低阻組合異常體不同測點視電阻率測深曲線見圖10(a)， 4 000 m和5 000 m 兩個測點的曲線形狀和圖8(a)和圖9(a)相似， 4 000 m比5 000 m較早進入過渡區(qū)，電阻率逐漸下降。4 240 m曲線和圖9(a)中4 520 m曲線相似，4 640 m曲線和圖8(a)中4 520 m曲線相似，說明對高低阻組合體的反應和單個異常體是相似的。

高低阻組合異常體不同頻率視電阻率曲線見圖10(b)：各頻率曲線在異常體上方的響應均和單個異常體的響應相似，具體分析見圖8(b)和圖9(b)。

高低阻組合異常體數(shù)值模擬結(jié)果電阻率擬斷面圖見圖10(c)：在4 120 m～4 320 m之間有明顯的低阻異常，中心位置在4 240 m處；在4 520 m～4 720 m之間有明顯的高阻異常，中心位置在4 640 m處。低阻異常體和高阻異常體的響應異常特征和單個異常體時，基本上具有相同的形態(tài)和變化規(guī)律。圖10(a)、圖10(b)、圖10(c)并結(jié)合單個異常體正演分析可知：CSAMT2.5-D正演算法對高低阻異常體組合模型和單個異常體的數(shù)值模擬具有相同響應異常特征和變化規(guī)律，可以模擬較為復雜的地電模型，并取得良好的效果。

圖9　低阻體模擬結(jié)果圖Fig．9　Low resistivity simulation results(a)不同測點曲線; (b)不同頻率曲線; (c)模擬結(jié)果擬斷面圖

圖10　高低阻組合體模擬結(jié)果圖Fig．10　High and low resistance group simulation results(a)不同測點曲線; (b)不同頻率曲線; (c)模擬結(jié)果擬斷面圖

7結(jié)論

1)從有源MAXWELL方程組出發(fā)，采用Gelerkin加權(quán)余值法建立CSAMT2.5-D有限元方程?；诰匦尉W(wǎng)格剖分，得到有限元系數(shù)矩陣，最終形成剛度矩陣。通過傅氏變換將三維源轉(zhuǎn)換為二維，根據(jù)電磁場值實虛部的曲線特征，按對數(shù)等間隔選取21個波數(shù)，計算電磁場響應，最終得到視電阻率。

2)推導了4節(jié)點無限單元的坐標映射函數(shù)和形函數(shù)，并將其與有限元網(wǎng)格結(jié)合，取代外邊界的Dirichlet 邊界條件，有效減少了有限元網(wǎng)格的剖分范圍。通過數(shù)值模擬均勻半空間，與一維模擬解析解結(jié)果進行對比，電阻率的均方相對誤差均小于1%，滿足精度要求。

3)應用CSAMT2.5-D正演算法三種地電模型進行了數(shù)值模擬，對單一的高、低阻體的模擬效果很好，中心位置基本能夠吻合；對高低阻組合異常體的數(shù)值模擬和單個異常體的數(shù)值模擬具有類似的響應異常特征和變化規(guī)律，證實這里研究的CSAMT2.5-D正演算法可以模擬較為復雜的地電模型，并取得良好的效果。

參考文獻：[1]何繼善.廣域電磁測深法研究[J].中南大學學報:自然科學版,2010,41(3):1065-1072.

HE J S.Wide field electromagnetic sounding methods[J].Journal of Central South University:Science and Technology,2010,41(3):1065-1072.(In Chinese)

[2]湯井田,周聰,肖曉.復雜介質(zhì)條件下CSAMT最小發(fā)收距的選擇[J].中國有色金屬學報,2013,23(6):1681-1693.

TANG J T,ZHOU C,XIAO X.Selection of minimum transmit-receive distance of CSAMT on complicated media[J].The Chinese Journal of Nonferrous Metals,2013,23(6):1681-1693.(In Chinese)

[3]湯井田,何繼善.可控源音頻大地電磁法及其應用[M].長沙:中南大學出版社,2005.

TANG J T,HE J S.Controlled source audio magnetotelluric method and its application[M].Changsha:Central South University Press,2005.(In Chinese)

[4]孟永良,羅延鐘.電偶源CSAMT法二維正演的有限單元算法[C]//勘查地球物理勘查地球化學文集,北京:地質(zhì)出版社,1996:103-114.

MENG Y L,LUO Y Z.Two-dimensional finite element method for electric dipole source CSAMT[C]//Essays on geophysics and geochemistry exploration.Beijing:Geological Publishing House,1996:103-114.(In Chinese)

[5]王若,王妙月,底青云.頻率域線源大地電磁法有限元正演模擬[J].地球物理學報,2006,49(6):1858-1866.

WANG R,WANG M Y,DI Q Y.Electromagnetic modeling due to line source in frequency domain using finite element method[J].Chinese Journal of Geophysics,2006,49(6):1858-1866.(In Chinese)

[6]張繼鋒,湯井田,喻言,等.基于二次插值的線源可控源有限元數(shù)值模擬[J].吉林大學學報:地球科學版,2009,39(5):929-935.

ZHANG J F,TANG J T,YU Y,et al.Finite element numerical simulation on line controlled source based on quadratic interpolation[J].Journal of Jilin University:Earth Science Edition,2009,39(5):929-935.(In Chinese)

[7]COGGON J H.Electromagnetic and electrical modeling by the finite element method[J].Geophysics,1971,36(1):132-155.

[8]UNSWORTH M J,TRAVIS B J,CHAVE A D.Electromagnetic induction by a finite electric dipole source over a 2-D earth[J].Geophysics,1993,58(2):198-214.

[9]LU X,UNSWORTH M,BOOKER J.Rapid relaxation inversion of CSAMT data[J].Geophys J Int,1999,138:381-392.

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

[11]MITSUHATA Y,UCHIDA T,AMANO H.2.5D inversion of frequency-domain electromagnetic data generated by a grounded-wire source[J].Geophysics,2002,67(6):1753-1768.

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

[13]陳金窗，戴光明.微機網(wǎng)絡(luò)并行計算及2.5維CSAMT正演的并行實現(xiàn)[J].物探化探計算技術(shù),1997,19(2):103-107.

CHEN J C,DAI G M.Micro-computer networked computing and 2.5-D CSAMT forward parallel computing[J].Computing Techniques for Geophysical and Geochemical Exploration,1997,19(2):103-107.(In Chinese)

[14]陳小斌，胡文寶.有限元直接迭代算法及其在線源頻率域電磁響應計算中的應用[J].地球物理學報，2002,45(1):119-130.

CHEN X B,HU W B.Direct iterative finite element algorithm and its application to electromagnetic response modeling of line current source[J].Chinese Journal of Geophysics,2002,45(1):119-130.(In Chinese)

[15]底青云,UNSWORTH M,王妙月.有限元法2.5 維CSAMT數(shù)值模擬[J].地球物理學進展,2004,19(2):317-324.

DI Q Y,UNSWORTH M,WANG M Y.2.5-D CSAMT modeling with the finite element method[J].Progress in Geophysics,2004,19(2):317-324.(In Chinese)

[16]底青云,UNSWORTH M,王妙月.復雜介質(zhì)有限元法2.5維CSAM T數(shù)值模擬[J].地球物理學報,2004,47(4):723-730.

DI Q Y,UNSWORTH M,WANG M Y.2.5-D CSAMT modeling with the finite element method over 2-D complex earth media[J].Chinese Journal of Geophysics,2004,47(4):723-730.(In Chinese)

[17]底青云,UNSWORTH M,王妙月.2.5維有限元法CSAM T數(shù)值反演[J].石油地球物理勘探,2006,41(1):100-106.

DI Q Y,UNSWORTH M,WANG M Y.2.5-D finite element method CSAMT numerical inversion[J].Oil Geophysical Prospecting,2006,41(1):100-106.(In Chinese)

[18]沈金松,孫文博.2.5維電磁響應的有限元模擬與波數(shù)取值研究[J].物探化探計算技術(shù),2008,30(2):135-144.

SHEN J S,SUN W B.A study on the 2.5-D electromagnetic modeling by the finite element method and the wave number selection[J].Computing Techniques forGeophysical and Geochemical Exploration,2008,30(2):135-144.(In Chinese)

[19]沈金松,孫文博.二維海底地層可控源海洋電磁響應的數(shù)值模擬[J].石油物探,2009,48(2):187-194.

SHEN J S,SUN W B. Numerical simulation of marine controlled-source electromagnetic (MCSEM) response to 2-Dseabed formation[J].Geophysics Prospecting for Petroleum,2009,48(2):187-194.(In Chinese)

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

NABIGIAN N M.Electromagnetic methods in appliedgeophysics (Volume 1)[M].ZHAO J X,WANG Y J,transl. Beijing:Geological Publishing House,1992.(In Chinese)

[21]葛偉男.2.5-D CSAMT有限元正演模擬[D].長沙：中南大學.2012.

Ge W N.Finite element numerical simulate of 2.5-D CSAMT[D].Changsha:Central South University,2012. (In Chinese)

[22]徐世浙.地球物理中的有限單元法[M].北京:科學出版社,1994.

XU S Z.The finite element method in geophysics[M].Beijing:Science Press,1994.(In Chinese)

[23]姚姚.地球物理反演基本理論與應用方法[M].武漢:中國地質(zhì)大學出版社,2005.

YAO Y.Basic theory on inversion and its application method[M]. Wuhan:China University of Geosciences,2005.(In Chinese)

[24]吳小平,徐果明.電阻率三維反演中偏導數(shù)矩陣的求取與分析[J].石油地球物理勘探,1999,34(4):363-372.

WU X P,XU G M. Derivation and analysis of partial derivative matrix in resistivity 3-D inversion[J].Oil Geophysical Prospecting,1999,34(4):363-372.(In Chinese)

[25]公勁喆.有限元·無限元耦合法在三維直流電和電磁數(shù)值模擬中的應用[D].長沙:中南大學信息物理工程學院,2009.

GONG J Z.Application of finite-infinite element coupling method on 3D electric and electromagnetic numerical modeling[D].Changsha:Central South University.School of Information Physical Engineering,2009.(In Chinese)

[26]ASTLEY R J,COYETTE J P.Conditioning of infinite element schemes for wave problems[J].Communications in Numerical Methods in Engineering,2001,17:31-41.

收稿日期：2015-03-31改回日期：2015-04-12

基金項目:中國地質(zhì)科學院地球物理地球化學勘查研究所基本科研業(yè)務(wù)專項資金項目(AS2013J01)

作者簡介：姚大為(1984-)，男，碩士，主要從事地球物理電磁法方法技術(shù)的工作和研究， E-mail: yaodawei@igge.cn。

文章編號：1001-1749(2016)03-0297-011

中圖分類號：P 631.3

文獻標志碼：A

DOI：10.3969/j.issn.1001-1749.2016.03.02

2.5-D finite element numerical simulation of CSAMT

YAO Da-wei, WANG Da-yong, WANG Gang, ZHU Wei, LI Yong-bo, ZHANG Zhen-yu

(Institute of Geophysical and Geochemical Exploration, CAGS The National Center of Geological Exploration Technology, Langfang065000, China)

Abstract:The forward calculation of 3-D field source 2-D geoelectric model was called 2.5-D numerical simulation. Starting from Maxwell's equations, respectively to solve the electromagnetic field of once field and secondary field in order to put the 3-d field source reduced to 2-D. By using Fourier transform wavenumber domain and spatial domain equations could be converted to finite element solution of wave number domain electromagnetic fields, introducing infinite element to solve the convergence problem of infinite boundary. According to characteristic curves real and imaginary parts of the electromagnetic field value, selected 21 wavenumbers at equal logarithmic. Comparison between the results of the code for the calculation of the homogeneous half space and the analytical solution, it was confirming effectiveness of the code that the resistivity of the mean squared relative error was less than 0.5%. Then to calculate the responses of three geoelectric model of electromagnetic field, simulation results of single high and low resistivity body and low resistance group fit true anomaly center position are consistent. This study confirmed that the CSAMT2.5-D forward algorithm can simulate complex geoelectric model, and achieved good results.

Key words:CSAMT; 2.5-D; finite element; infinite element