基于雙二次插值的2.5維FCSEM有限元正演模擬

柳建新 ，湯文武 ，童孝忠

(1. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙，410083；2. 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙，410083)

2.5 維數(shù)值模擬指的是源為三維的、正演模型為二維的地球物理數(shù)值模擬問題[1?3]。眾所周知，2.5維數(shù)值模擬問題有著重要的實(shí)用價(jià)值，這是因?yàn)閷⑷S問題簡(jiǎn)化為2.5維模型后只需要對(duì)剖面而不是整個(gè)體積空間進(jìn)行離散處理，這樣就大大減小了矩陣的尺寸[4]。

2.5維頻率域可控源電磁法(frequency domain controlled-source electromagnetic, FCSEM)的有限元正演模擬，國(guó)內(nèi)外都進(jìn)行了研究。Unsworth等[5]將一次場(chǎng)與二次場(chǎng)分離，利用有限元對(duì)二維地電模型下水平電偶源的電磁場(chǎng)響應(yīng)進(jìn)行了模擬計(jì)算；Mitsuhata[6]借鑒地震模擬中的偽delta函數(shù)模擬電偶源，只需一次性計(jì)算總場(chǎng)，并利用等參有限元實(shí)現(xiàn)了帶地形的2.5維頻率域電磁正演模擬；底青云等[7?8]對(duì)2.5維頻率域電磁進(jìn)行了有限元模擬并給出了一組波數(shù)，對(duì)層狀模型的模擬計(jì)算效果較好；薛東川等[9?10]對(duì)2.5維頻率域電磁的吸收邊界條件進(jìn)行了討論，發(fā)現(xiàn)吸收邊界條件能夠有效壓制邊界處的反射。在此，本文作者在前人研究基礎(chǔ)上，重點(diǎn)研究基于雙二次插值的FCSEM有限元正演模擬，并對(duì)波數(shù)選取進(jìn)行研究，給出一組兼顧計(jì)算速度及模擬精度的波數(shù)。

1 有限元方程的推導(dǎo)

圖1所示為2.5維FCSEM正演模擬示意圖：其y方向?yàn)樽呦蚍较?，水平電偶源Idy沿走向方向布置。電磁法勘探是基于電磁場(chǎng)理論的一類地球物理勘探方法，2.5維頻率域可控源電磁(FCSEM)也應(yīng)滿足麥克斯韋方程。假設(shè)時(shí)諧因子為eiωt，則麥克斯韋方程表示如下[3]：

式中：E為電場(chǎng)強(qiáng)度；H為磁場(chǎng)強(qiáng)度；Js為外加電偶源；阻抗率=iμω；導(dǎo)納率=σ+εω。對(duì)式(1)和(2)中的所有電場(chǎng)分量及磁場(chǎng)分量都進(jìn)行傅里葉變換得：

可以得到如下的微分控制方程[6]：

圖1 2.5-D FCSEM模擬示意圖Fig. 1 Simulation schematic diagram of 2.5-D FCSEM

首先利用有限元方法求解式(4)和(5)可得到傅氏域中的電場(chǎng)分量y及磁場(chǎng)分量y，再利用式(6)～(9)求得其余4個(gè)分量，最后通過傅氏逆變換可求得真實(shí)電磁場(chǎng)。這里直接給出由伽遼金有限元方法推導(dǎo)得到的有限元方程如下[6]：

2 有限元求解過程

先對(duì)模擬區(qū)域進(jìn)行剖分。由于同等計(jì)算量下雙二次插值一般較雙線性插值精度高，因此，本文采用雙二次插值以提高計(jì)算精度；接著對(duì)各個(gè)單元進(jìn)行分析、計(jì)算，得到各個(gè)小型剛度矩陣；最后，對(duì)所有單元總體集成得到總體剛度矩陣[11?15]。

2.1 網(wǎng)格剖分及插值

首先在研究區(qū)域內(nèi)進(jìn)行矩形單元剖分。為了提高計(jì)算精度及減小利用式(6)～(9)求解時(shí)帶來的誤差，本文采用精度較高的雙二次插值進(jìn)行模擬。圖2所示為坐標(biāo)變換關(guān)系及局部節(jié)點(diǎn)編號(hào)，其中圖 2(a)所示為母單元，圖2(b)所示為子單元；a和b分別表示子單元的寬度及高度；(x0，z0)為矩形單元中心點(diǎn)的坐標(biāo)。文中采用的插值函數(shù)表示如下：

圖2 矩形單元剖分示意圖Fig. 2 Schematic diagram of rectangle meshing element

子單元與母單元之間的坐標(biāo)變換關(guān)系式表示如下：

2.2 單元分析

首先對(duì)式(10)進(jìn)行單元分析。

記U1=(u1，u2，…，u8)T，V1=(v1，v3，…，v8)T，則

式(10)中左邊第1項(xiàng)：

式(10)中左邊第2項(xiàng)：

式(10)中左邊第3項(xiàng)：

式(10)中左邊第4項(xiàng)：

式(10)中左邊第5項(xiàng)：

再對(duì)式(11)進(jìn)行單元分析。

式(11)中左邊第1項(xiàng)：

式(11)中左邊第2項(xiàng)：

式(11)中左邊第3項(xiàng)：

式(11)中左邊第4項(xiàng)：

式(11)中左邊第5項(xiàng)：

觀察式(15)～(24)，借助符號(hào)計(jì)算工具maple可求得：

2.3 總體合成

將式(15)～(19)中各項(xiàng)累加，從而可將式(10)化為

其中：J代表電流源項(xiàng)。將式(20)～(24)中各項(xiàng)累加，從而可將式(11)化為

再對(duì)式(30)及(31)化簡(jiǎn)，令

則式(30)和(31)可記為：

將式(36)等價(jià)為

將式(37)中左邊的剛度矩陣記作

經(jīng)分析可知：M矩陣為復(fù)對(duì)稱矩陣。求解線性方程組(38)即可得到傅氏域的 2個(gè)電磁場(chǎng)分量，再利用式(6)～(9)求得傅氏域中的其余4個(gè)電磁場(chǎng)分量；最后，利用傅里葉逆變換求得各個(gè)電磁場(chǎng)分量。

3 波數(shù)的選取

2.5 維數(shù)值模擬是先在波數(shù)域里計(jì)算電磁場(chǎng)，再通過傅里葉逆變換求得空間域的電磁場(chǎng)。由于在數(shù)值模擬時(shí)不可能取無窮多個(gè)波數(shù)來還原真實(shí)的電磁場(chǎng)，因此，需要選取若干個(gè)合適的波數(shù)以盡可能、快速地求得真實(shí)的電磁場(chǎng)。如何選取波數(shù)及確定波數(shù)的選取范圍便是一個(gè)關(guān)鍵性問題。

本文通過實(shí)際計(jì)算及分析，選取一組兼顧計(jì)算速度與精度的波數(shù)。在收發(fā)距為4 km，頻率f為2 Hz和1 024 Hz時(shí)，分別進(jìn)行正演模擬計(jì)算并繪制了及隨波數(shù)的變化曲線圖，如圖3所示。其中，圖3(a)所示為隨波數(shù)的變化曲線圖，圖3(b)所示為隨波數(shù)的變化曲線圖。從圖3(a)可知：當(dāng)波數(shù)較小時(shí)，基本保持不變；隨著波數(shù)繼續(xù)增大時(shí)，開始下降；當(dāng)波數(shù)約為0.01時(shí)到達(dá)波谷，而后又稍有回升。由圖 3(b)可知：當(dāng)波數(shù)較小時(shí)，很小；隨著波數(shù)增大，開始先上升，到達(dá)波峰后以較陡傾角下降；同樣地，當(dāng)波數(shù)約為0.01時(shí)到達(dá)低谷，而后又稍有回升，但變化曲線沒有規(guī)律。據(jù)文獻(xiàn)[6]，波谷右側(cè)的曲線是不可靠的。事實(shí)上，曲線到達(dá)波谷后繼續(xù)往下遞減，且由于此時(shí)取值已經(jīng)很小，對(duì)傅里葉逆變換中的積分結(jié)果的影響可以忽略，因此，只需要在波谷左邊部分選取足夠的能夠精確計(jì)算積分的一組波數(shù)即可。本文通過不斷進(jìn)行正演計(jì)算及對(duì)比分析，選取由20個(gè)波數(shù)構(gòu)成的一組波數(shù)，如表1所示。

圖3 傅氏域中電場(chǎng)和磁場(chǎng)隨波數(shù)變化示意圖Fig. 3 Schematic diagram of electric field and magnetic field in Fourier domain of changing with wave number

表1 波數(shù)值選取Table 1 Choice of wavenumbers

事實(shí)驗(yàn)證，這一組波數(shù)能夠有效且快速、精確計(jì)算2.5維FCSEM的電磁場(chǎng)響應(yīng)。本文后續(xù)的工作都是基于這一組波數(shù)進(jìn)行的。

4 數(shù)值算例

4.1 算法驗(yàn)證

為了驗(yàn)證本文算法的正確性，對(duì)1個(gè)電阻率為100?·m的均勻半空間模型進(jìn)行正演模擬計(jì)算及與解析解的對(duì)比分析。圖4所示為解析解與數(shù)值解對(duì)比結(jié)果圖(其中，x軸代表的是收發(fā)距)。由圖4可知：電磁場(chǎng)分量數(shù)值解與解析解相吻合，且頻率低時(shí)精度更高。計(jì)算結(jié)果表明：對(duì)于電場(chǎng)分量，當(dāng)f=2 Hz時(shí)，最小相對(duì)誤差為 0.106%，最大相對(duì)誤差為 0.639%，平均相對(duì)誤差為0.460%；當(dāng)f=1 024 Hz時(shí)，最小相對(duì)誤差為0.451%，最大相對(duì)誤差為 1.763%，平均相對(duì)誤差為0.920%。對(duì)于磁場(chǎng)分量，當(dāng)f=2 Hz時(shí)，最小相對(duì)誤差為 0.016%，最大相對(duì)誤差為 0.513%，平均相對(duì)誤差為0.220%；當(dāng)f=1 024 Hz時(shí)，最小相對(duì)誤差為0.545%，最大相對(duì)誤差為 1.921%，平均相對(duì)誤差為 1.028%。可見，完全滿足正演計(jì)算的要求。

4.2 模型計(jì)算

4.2.1 低阻體模型

圖5所示為 1個(gè)低阻體模型的示意圖，其中，x方向標(biāo)出的刻度代表垂直收發(fā)距。圖6和圖7所示為正演模擬結(jié)果。其中，圖6(a)和圖6(b)所示分別為TE模式下的視電阻率及相位等值線圖，這2個(gè)圖中都反映出了異常，且異常的中心位置與理論模型一致，異常都有往右下方拉長(zhǎng)的趨勢(shì)，這正是場(chǎng)源陰影效應(yīng)的影響；視電阻率等值線圖異常最低為30 ?·m，接近異常體的真實(shí)值。圖7(a)和圖7(b)所示分別為TM模式下的視電阻率及相位等值線圖，這2個(gè)圖中也都反映出了異常，且相比 TE模式異常垂向拉長(zhǎng)明顯。這是由于TM模式下測(cè)定的是電場(chǎng)的法向分量，界面積累電荷效應(yīng)使得受靜態(tài)效應(yīng)明顯，造成垂向拉長(zhǎng)，這提高了方法的靈敏度，有利于發(fā)現(xiàn)小規(guī)模低阻體，故TM模式發(fā)現(xiàn)異常的能力比TE模式的強(qiáng)。

圖4 均勻半空間模型解析解與數(shù)值解對(duì)比結(jié)果Fig. 4 Comparison between analytic solution and numerical solution of homogenous half space

圖5 低阻體模型示意圖Fig. 5 Schematic diagram of low resistivity body

圖6 TE模式下低阻體模擬結(jié)果Fig. 6 Schematic simulation result of low resistivity body of TE mode

4.2.2 高阻體模型

圖7 TM模式下低阻體模擬結(jié)果Fig. 7 Simulation results of low resistivity body of TM mode

圖8 高阻體模型示意圖Fig. 8 Schematic diagram of high resistivity body

圖8 所示為1個(gè)高阻體模型的示意圖，圖9和圖10所示為正演模擬結(jié)果。其中，圖9(a)和(b)所示分別為TE模式下的視電阻率及相位等值線圖，這2個(gè)圖中都反映出了異常，視電阻率等值線圖中異常最高為130 ?·m，而真實(shí)異常為500 ?·m，說明頻率域可控源電磁對(duì)高阻體反映不靈敏，這是電磁法類勘探方法普遍存在的一個(gè)問題。圖 10(a)和 10(b)所示分別為 TM模式下的視電阻率及相位等值線圖，這2個(gè)圖中都反映出異常，視電阻率等值線圖中異常最高為200 ?·m，說明TM模式較TE模式對(duì)高阻體反映較靈敏。

圖9 TE模式下高阻體模擬結(jié)果Fig. 9 Simulation results of high resistivity body of TE mode

圖10 TM模式下高阻體模擬結(jié)果Fig. 10 Simulation results of high resistivity body of TM mode

5 結(jié)論

(1) 通過均勻半空間模型下低頻及高頻時(shí)的有限元數(shù)值解與解析解對(duì)比，結(jié)果顯示有限元數(shù)值解與解析解基本吻合，從而從一定程度上驗(yàn)證了基于雙二次插值的2.5維FCSEM正演模擬算法的正確性。

(2) 本文2.5維FCSEM有限元數(shù)值模擬結(jié)果能夠正確反映二維模型的淺部異常及構(gòu)造，進(jìn)一步驗(yàn)證了該算法的可靠性，可以有效地模擬2.5維地電模型。

(3) TM模式下低阻及高阻異常體響應(yīng)均有往下拉長(zhǎng)的趨勢(shì)，對(duì)低阻及高阻異常體反應(yīng)較TE模式明顯，且橫向分辨能力更強(qiáng)。

[1]底青云, 王若. 可控源音頻大地電磁數(shù)據(jù)正反演及方法應(yīng)用[M]. 北京: 科學(xué)出版社, 2008: 74?84.DI Qingyun, WANG Ruo. Forward modeling and inversion of CSAMT data and its application in earth[M]. Beijing: Science Press, 2008: 74?84.

[2]湯井田, 何繼善. 可控源音頻大地電磁法及其應(yīng)用[M]. 長(zhǎng)沙:中南大學(xué)出版社, 2005: 9?33.TANG Jingtian, HE Jishan. Audio-frequency of controlled source electromagnetic method and its application in earth[M].Changsha: Central South University Press, 2005: 9?33.

[3]米薩克 N 納比吉安. 勘查地球物理電磁法. 1卷: 理論)[M].趙經(jīng)祥, 王艷君, 譯. 北京: 地質(zhì)出版社, 1992: 330?333.Nabighian M. Electromagnetic methods in applied geophysics:Ⅰ. Theory[M]. ZHAO Jingxiang, WANG Yanjun, trans. Beijing:Geological Publishing House, 1992: 330?333.

[4]熊彬. 關(guān)于瞬變電磁法2.5維正演中的幾個(gè)問題[J]. 物探化探計(jì)算技術(shù), 2005, 28(2): 124?128.XIONG Bin. Several problems related to 2.5 dimensional forward modeling of transient electromagnetic method[J].Computing Techniques for Geophysical and Geochemical Exploration, 2005, 28(2): 124?128.

[5]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.

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

[7]底青云, Unsworth M, 王妙月. 復(fù)雜介質(zhì)有限元法2.5維可控源音頻大地電磁法數(shù)值模擬[J]. 地球物理學(xué)報(bào), 2004, 47(4):723?730.DI Qingyun, Unsworth M, WANG Miaoyue. 2.5 D CSAMT modeling with the finite element method over 2D complex earth media[J]. Journal of Geophysics, 2004, 47(4): 723?730.

[8]底青云, Unsworth M, 王妙月. 有限元法2.5維CSAMT數(shù)值模擬[J]. 地球物理學(xué)進(jìn)展, 2004, 19(2): 317?324.DI Qingyun, Unsworth M, WANG Miaoyue. 2.5 D CSAMT modeling with finite element method[J]. Progress in Geophysics,2004, 19(2): 317?324.

[9]薛東川, 戴世坤. 頻率域2.5維電磁測(cè)深有限元模擬中的吸收邊界條件[J]. 中國(guó)石油大學(xué)學(xué)報(bào)(自然科學(xué)版), 2008, 32(6):57?61.XUE Dongchuan, DAI Shikun. Absorbing boundary condition for simulating 2.5-D electromagnetic sounding in frequency domain by finite element method[J]. Journal of China University o f Petroleum (Natural Science Edition), 2008, 32(6): 57?61.

[10]王華軍, 羅延鐘. 中心回線瞬變電磁法2.5維有限單元算法[J].地球物理學(xué)報(bào), 2003, 46(6): 855?862.WANG Huajun, LUO Yanzhong. Algorithm of a 2.5 dimensional finite element method for transient electromagnetic with a central loop[J]. Journal of Geophysics, 2003, 46(6):855?862.

[11]TONG Xiaozhong, LIU Jianxin, XIE Wei, et al.Three-dimensional forward modeling for magnetotelluric sounding by finite element method[J]. Journal of Central South University of Technology, 2009, 16(1): 136?142.

[12]柳建新, 蔣鵬飛, 童孝忠, 等. 不完全 LU 分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2009, 40(2): 484?491.LIU Jianxin, JIANG Pengfei, TONG Xiaozhong, et al.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling[J]. Journal of Central South University (Science and Technology), 2009, 40(2): 484?491.

[13]TANG Jingtian, WANG Feiyan, REN Zheng-yong, et al. 3-D direct current resistivity forward modeling by adaptive multigrid finite element method[J]. Journal of Central South University of Technology, 2010, 17(3): 587?592.

[14]劉長(zhǎng)生, 湯井田, 任政勇, 等. 基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2010, 41(5): 1855?1859.LIU Changsheng, TANG Jingtian, REN Zhengyong, et al.Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes[J]. Journal of Central South University (Science and Technology), 2010, 41(5):1855?1859.

[15]張繼鋒, 湯井田, 王燁, 等. 基于預(yù)處理共軛梯度的大地電磁快速正演[J]. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2010, 41(5):1877?1882.ZHANG Jifeng, TANG Jingtian, WANG Ye, et al.Magnetotellurics fast forward based on preconditioning conjugate gradient[J]. Journal of Central South University(Science and Technology), 2010, 41(5): 1855?1859.