原源， 湯井田， 任政勇*， 周聰， 肖曉， 張林成

1 中南大學有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室， 長沙 4100832 中南大學地球科學與信息物理學院應(yīng)用地球物理系， 長沙 410083

?

基于非結(jié)構(gòu)化網(wǎng)格的任意復雜2D RMT有限元模擬

原源1,2， 湯井田1,2， 任政勇1,2*， 周聰1,2， 肖曉1,2， 張林成1,2

1 中南大學有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室， 長沙 4100832 中南大學地球科學與信息物理學院應(yīng)用地球物理系， 長沙 410083

射頻大地電磁法（RMT）是以潛艇天線發(fā)射的射線源等作為場源的一種地球物理勘探方法，目前被廣泛應(yīng)用于近地表工程和環(huán)境地球物理勘探.RMT數(shù)據(jù)解釋常采用基于靜態(tài)假設(shè)的大地電磁法（MT）程序，往往會反演出不真實的淺層目標體.為解決這一長期困擾RMT資料解釋的問題，本文實現(xiàn)了考慮位移電流的RMT有限元數(shù)值模擬.為了處理任意復雜模型，如起伏地形，非結(jié)構(gòu)的三角形網(wǎng)格被用于離散RMT模型.首先通過算例對比，驗證了程序的正確性和可靠性.通過Dike模型討論了空氣層厚度對RMT數(shù)值解的影響，結(jié)果表明當空氣層厚度大于1/4波長即可滿足精度要求.以山脊模型為例計算了位移電流對RMT響應(yīng)的影響，表明位移電流的影響會隨著山脊高程的增加而增大. 最后通過舒家店實際模型進一步驗證了位移電流在RMT中的重要性.

大地電磁； RMT； 有限元； 非結(jié)構(gòu)網(wǎng)格； 位移電流

1 引言

射頻大地電磁法（Radio-magnetotelluric）是近幾年才逐步發(fā)展的一種淺地表地球物理勘探方法.RMT法的測量頻段為10～250 kHz，勘探深度為數(shù)米到數(shù)十米，目前多應(yīng)用于工程勘探及環(huán)境監(jiān)測（Pedersen et al.,2005,2006; Tezkan et al.,2000,2005; Tezkan and Saraev,2008; Ismail et al.,2011a,2011b; Bastani et al.,2013）.

由于RMT法與MT法都是在相應(yīng)源的遠區(qū)觀測水平電磁場，因而目前在進行RMT資料處理與解釋時多是直接套用MT的正反演程序（Newman et al.,2003; Linde and Pedersen,2004; Candansayar and Tezkan,2006,2008），這會導致反演得到的地下電阻率值出現(xiàn)不切實際的極小極大值（Persson and Pedersen, 2002; Kalscheuer et al.,2008），從而得到不合理的解釋結(jié)果.由于MT法測量頻率在數(shù)百赫茲內(nèi)，在該頻段，電磁場以感應(yīng)擴散為主，因此在數(shù)值模擬時通常進行準靜態(tài)假設(shè)，忽略位移電流.RMT法的測量頻段為10～250 kHz，在該頻段內(nèi)，位移電流不是遠遠小于傳導電流，因而不可忽略位移電流，當?shù)叵聨r體為結(jié)晶類高阻巖石（如石英巖等）時尤其如此（Chave and Jones,2012）.因此，實現(xiàn)考慮位移電流的RMT正演程序具有重要意義.

Persson和Pedersen （2002）討論了一維RMT響應(yīng)中，位移電流對視電阻率和相位的影響，表明考慮位移電流得到的視電阻率和相位值均低于準靜態(tài)計算結(jié)果，同時文中對層狀模型進行反演，得出結(jié)論：不考慮位移電流的反演電阻率與實際偏差較大.繼均勻半空間和層狀模型后，Kalscheuer等（2008）首次研究了二維介質(zhì)中位移電流對RMT響應(yīng)的影響，說明了在RMT頻段內(nèi)，當?shù)叵码娮杪蚀笥?0000 Ωm后，位移電流會對計算精度產(chǎn)生較大影響，視電阻率和相位值均低于準靜態(tài)結(jié)果，同時準靜態(tài)條件下的反演結(jié)果會帶來虛假構(gòu)造，導致錯誤的解釋.然而，Kalscheuer等（2008）的工作均是基于平地形的有限差分正演程序，在野外實際工作中，通常面臨的是復雜起伏地形，因而開發(fā)帶地形的RMT程序更具實際意義.

為了解決復雜起伏地形情況下RMT的數(shù)據(jù)解釋問題，本文開發(fā)了基于非結(jié)構(gòu)三角形網(wǎng)格的2D帶地形RMT有限元正演程序.首先通過一算例對比驗證了程序的正確性.由于RMT模擬中存在位移電流，那么TM模式的計算也必須考慮空氣層的影響，為此，文中討論了高精度RMT模擬中空氣層厚度的最佳選擇.接著，文中討論了起伏地形下位移電流對數(shù)值解的影響.最后通過一實際算例進一步說明在RMT數(shù)值模擬中必須考慮位移電流的影響.

2 含空氣層的二維RMT問題

2.1 邊值問題

圖1 求解域Ω=Ω1∪Ω2，上邊界Γ1，外邊界ΓD為無窮遠邊界，Γint為內(nèi)邊界，x為二維構(gòu)造走向方向

假設(shè)x方向為構(gòu)造走向方向（圖1），z方向垂直向下，電磁場分量只沿y、z方向有變化，因此二維構(gòu)造下，電磁場方程可寫為（取時間因子eiω t）：

（1）

（2）

（3）

（4）

（5）

（6）

其中ω為角頻率，σ為電導率，ε為介電常數(shù)，μ0為磁導率.（1）—（3）式定義了TE模式，（4）—（6）式定義了TM模式.

將（6）式兩邊求?/?y，（7）式兩邊求?/?z，并代入（5）式即得TE模式下電場Ex分量滿足的Helmholtz方程：

（8）

（7）式（8）式可統(tǒng)一表示為如下所示的偏微分方程：

（9）

在低頻情況下，由（4）、（5）式可知，磁場在空氣中幾乎不變化，因此對于低頻問題，如MT問題，TM模式一般不考慮空氣.對于RMT來說，由于位移電流在空氣中占主導，因此由（4）、（5）式可知，磁場在空氣中存在變化，因此RMT問題的TE和TM模式，均要考慮空氣層.

本文中計算區(qū)域如圖1所示，外邊界（ΓD）上電位通過層狀介質(zhì)解析解給出，上邊界（Γ1）取u=1，內(nèi)邊界（Γint）根據(jù)場切向分量連續(xù)給出.綜上，邊值問題中的邊界條件如下：

u=r（x,z）,onΓD

（10）

（11）

2.2 有限元近似

網(wǎng)格剖分可分為兩種類型，結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格.結(jié)構(gòu)化的網(wǎng)格在區(qū)域加密時會導致不必要的外部節(jié)點增加，從而增加額外的計算量，且結(jié)構(gòu)化網(wǎng)格在對復雜地形及地下構(gòu)造模擬中不夠靈活（Li and Spitzer,2002; Mitsuhata and Uchida,2004；譚捍東等，2003；胡祖志等，2006；胡祥云等，2012；董浩等，2014）.因而，非結(jié)構(gòu)網(wǎng)格以其靈活性近些年廣泛應(yīng)用于地球物理數(shù)值模擬領(lǐng)域（Franke et al.,2007; Ren and Tang, 2010；Ren et al.，2013；劉長生等，2010；吳小平等，2015）.

本文中，求解域Ω進行網(wǎng)格離散時采用非結(jié)構(gòu)的三角形網(wǎng)格（Shewchuk,1996），以保證靈活模擬復雜地質(zhì)構(gòu)造及地形.在離散過程中，通過設(shè)定每個三角形單元的最小角不低于30°來保證單元質(zhì)量，同時通過給定局部區(qū)域單元片面積最大值來實現(xiàn)網(wǎng)格的局部加密.進行網(wǎng)格剖分時，首先根據(jù)最低頻率確定最大趨膚深度δmax，剖分區(qū)域邊界距離不低于10δmax；在電阻率突變界面附近控制單元大小進行局部加密；與MT法不同，RMT法計算TM模式的磁場時，空氣中的磁場不能看做常數(shù)，因此TM模式也需要考慮空氣層，而空氣中的網(wǎng)格剖分也要仔細考慮.在高頻情況下，電磁場在空氣中以電磁波的形式向外傳播，空氣中電磁波的傳播速度為光速c=3×108m·s-1，根據(jù)c=λ·f，其中λ為波長，f為頻率，據(jù)此可計算出不同頻率下的波長，本文中網(wǎng)格剖分時最短波長內(nèi)保證有20～30個節(jié)點.

基于線性三角形有限單元，通過求離散邊值問題，從而獲得關(guān)于方程（9）中未知數(shù)u的線性方程組：

AU=0,

（12）

式中，U為方程（9）中u的一維向量，剛度矩陣A為稀疏對稱復數(shù)矩陣.本文通過三個一維數(shù)組僅存儲矩陣的非零元素，從而節(jié)省計算內(nèi)存.采用Intel Math Kernel Library（Intel，2010）中的直接求解器PARDISO（Schenk and G?rtner，2004）來求解線性方程組.求出TE及TM模式的電場Ex和磁場Hx后，本文采用四點差分求得相應(yīng)的磁場Hy和電場Ey.最后，阻抗、視電阻率ρij和相位φij通過如下計算公式獲得：

（13）

（14）

φij=arctan（Im（Zij）/Re（Zij））,ij=xy,yx,

（15）

其中，Re、Im分別表示復數(shù)的實部與虛部.

3 數(shù)值算例

本節(jié)首先計算了Kalscheuer等（2008）文中的模型，并進行了結(jié)果對比，驗證了程序正確性；然后以Dike模型為例研究了空氣層厚度對數(shù)值解的影響；接著討論了起伏地形下，位移電流的影響規(guī)律；最后通過舒家店的實際地質(zhì)模型說明了位移電流對實際數(shù)據(jù)的影響.文中所有計算均在CPU為2.3GHz，RAM為2GB的個人計算機上完成.

3.1 均勻半空間中一矩形異常體

圖2 均勻半空間中賦存一矩形異常體，半空間電阻率為10000 Ωm，相對介電常數(shù)為5，矩形電阻率為1000 Ωm

圖4、圖5分別為TE、TM模式下計算得到的視電阻率及相位曲線，其中線條為本文有限元計算結(jié)果，符號為Kalscheuer等（2008）的有限差分程序計算結(jié)果，二者均考慮位移電流.從圖中可看出本文計算結(jié)果與Kalscheuer等（2008）的結(jié)果吻合，表明程序正確.

3.2 空氣層厚度對RMT數(shù)值解的影響

在RMT有限元模擬中，由于位移電流項的存在，即使TM模式也不可忽略空氣層，為此，本節(jié)討論空氣層厚度對數(shù)值解的影響.圖6為Dike模型示意圖，背景電阻率ρ2=10000 Ωm，相對介電常數(shù)為5，中間斷層電阻率ρ1=1000 Ωm.

圖3 局部網(wǎng)格剖分示意圖

圖4 圖2所示模型考慮位移電流的TE模式（a）視電阻率及（b）相位曲線，頻率分別為10 kHz、100 kHz、250 kHz

圖5 圖2所示模型考慮位移電流的TM模式（a）視電阻率及（b）相位曲線，圖例與圖4相同

圖6 Dike模型示意圖

表1為不同空氣層厚度下網(wǎng)格剖分后的單元節(jié)點信息、計算誤差及計算耗時.計算頻率為250 kHz，空氣層厚度從0.01λ～10λ，其中λ為波長，根據(jù)c=λ·f計算而來，c為光速，因此λ=1200 m.從表1中可看出，當空氣層厚度大于0.2λ， TE/TM模式視電阻率及相位的最大誤差均在1%以內(nèi)；當空氣層厚度小于0.2λ時，TE/TM模式視電阻率及相位的最大相對誤差均隨空氣層厚度的減小而增大.因此，在數(shù)值計算中，保證空氣層厚度大于0.2λ即可滿足精度要求.

表1 對不同空氣層厚度剖分后的網(wǎng)格信息、計算誤差及計算量

圖7、圖8分別為不同空氣層厚度下計算的所有測點的TE/TM模式視電阻率、相位及其相應(yīng)誤差.圖中展示的是幾個典型空氣層厚度的結(jié)果.所有的誤差計算均與10λ的結(jié)果進行對比，無論TE還是TM模式，當空氣層厚度大于0.2λ時，TE/TM模式的視電阻率及相位所有測點的誤差均在1%以內(nèi)，而當空氣層厚度為0.1λ時，TE模式視電阻率最大誤差為1.66%，相位最大誤差為2.08%；TM模式視電阻率最大誤差為1.65%，相位最大誤差為1.50%.對比TE模式所有測點的誤差曲線可發(fā)現(xiàn)，除了局部個別測點外，幾乎所有測點的誤差均隨著空氣層厚度的減小而增大，而TM模式的誤差曲線對這一規(guī)律的反映則不夠明顯，這是因為TM模式受橫向異常體影響較大所致.

3.3 起伏地形模型

本節(jié)采用圖9所示的模型（Wannamaker et al.,1986）討論不同高程下位移電流的影響.圖9中共給出了3個測試模型，余弦型山脊的高程分別為100 m、300 m、600 m，電阻率均為10000 Ωm，相對介電常數(shù)為5，計算的三個頻率為10 kHz、100 kHz、250 kHz.

為保證結(jié)果可靠，本節(jié)三個模型均采用較密的網(wǎng)格剖分，以降低計算誤差.網(wǎng)格剖分參數(shù)和計算耗時如表2所示.

表2 起伏地形模型網(wǎng)格剖分參數(shù)表

圖7 不同空氣層厚度下得到的TE模式（a）視電阻率、（c）相位及其相應(yīng)誤差（b）（d）.計算頻率f=250 kHz，空氣層厚度分別取10λ、1λ、0.25λ、0.1λ；計算誤差均與10λ的結(jié)果進行對比

圖8 同圖7，但為Dike模型的TM模式結(jié)果對比

圖9 三個模型均為余弦型山脊，山脊高程分別為100 m、300 m、600 m，三個模型的電阻率均為10000 Ωm，相對介電常數(shù)為5

圖10為考慮位移電流和準靜態(tài)情況下計算的視電阻率及相位曲線，計算頻率為10 kHz、100 kHz、250 kHz.從不同頻點的結(jié)果對比來看，10 kHz下考慮位移電流和不考慮位移電流的曲線都基本重合，說明本算例中位移電流在10 kHz的頻段處影響不大，而計算頻率為100 kHz、250 kHz時，考慮位移電流的視電阻率曲線和相位曲線均出現(xiàn)下掉現(xiàn)象，其中相位曲線更為嚴重，表明對于本算例中的模型參數(shù)而言，當頻率大于100 kHz后，位移電流對視電阻率的影響不可忽略，而對相位曲線的影響從頻率為20 kHz后就必須考慮.

圖11是頻率為250 kHz下考慮位移電流和準靜態(tài)條件下計算得到的視電阻率及相位絕對誤差，將三個模型計算結(jié)果對比，發(fā)現(xiàn)當山脊高程越大，坡度越陡，位移電流的影響就越大.從圖11a,11c可看出，TE模式的視電阻率及相位的最大誤差發(fā)生在山脊兩側(cè)，而山脊頂部則出現(xiàn)局部極大值；而圖11b，11d顯示TM模式視電阻率及相位的最大誤差也發(fā)生在山脊兩側(cè)，但山脊頂部無較明顯的極值.

3.4 實際模型

實際模型來自安徽銅陵礦集區(qū)舒家店地區(qū)某剖面（圖12），該剖面長3.8 km，具備以下兩個特點：該區(qū)域存在高阻閃長巖，位移電流可能會造成顯著影響；地形起伏明顯，有兩個山脊，地形對位移電流的影響不可忽略.

圖14為TE/TM模式下考慮位移電流和準靜態(tài)條件下的視電阻率及相位誤差.從圖中可看出，誤差整體隨著頻率的增加而增大；在x=-1000 m和x=500 m的高阻石英閃長巖體區(qū)域，考慮位移電流的結(jié)果和準靜態(tài)結(jié)果存在明顯的誤差，而在x=-1000 m的山脊處這種誤差更大.值得一提的是，相位受位移電流影響比視電阻率更為顯著，即使x=1300 m處，粉砂巖電阻率為420 Ωm的情況下，在高頻和山脊的共同作用下，相位也出現(xiàn)了明顯的誤差.ρTE的最大相對誤差為14.5%，位于x=-950 m處，而在該處附近的誤差也都在10%左右，此外，在x=-1600 m、x=500 m以及x=1500 m附近的誤差也都在5%以上；ρTM的最大相對誤差高達78%，這是因為TM模式橫向分辨率高，在電阻率分界面處二者計算結(jié)果有較大誤差，ρTM也主要分布在兩個山脊處和x=400 m的分界面處.

圖10 （a）（c）為model-1的TE/TM模式的視電阻率曲線，（b）（d）為model-1的TE/TM模式的相位曲線，計算頻率為10 kHz、100 kHz、250 kHz

圖11 考慮位移電流和不考慮位移電流的視電阻率及相位誤差對比，頻率為250 kHz

圖12 舒家店地質(zhì)綜合剖面圖（安徽省地質(zhì)調(diào)查院提供）

圖13 舒家店地球物理模型及網(wǎng)格剖分示意圖

圖14 TE/TM模式下考慮位移電流和準靜態(tài)條件下的視電阻率及相位誤差曲線

圖15 考慮位移電流條件下計算的位移電流密度和傳導電流密度，計算頻率f=250 kHz

圖15為考慮位移電流時計算的位移電流密度和傳導電流密度，計算頻率為250 kHz.圖15a，15b為電流密度的模.圖15c為位移電流密度占總電流密度的百分比.從圖15c中可看出，空氣中位移電流密度占100%，也就是說空氣中僅存在位移電流，無傳導電流.在地下砂巖、粉砂巖對應(yīng)的低阻區(qū)，位移電流密度在總電流密度中所占的比例為1%.花崗閃長巖區(qū)域的位移電流密度占3%左右.對于電阻率很高的石英閃長巖，位移電流密度所占比例達到10%，因此，在RMT數(shù)值模擬中不可忽略位移電流.

4 討論與結(jié)論

大多數(shù)RMT的資料解釋多是直接應(yīng)用MT程序，忽略了位移電流的影響，這可能導致錯誤的解釋結(jié)果；并且，在起伏地形條件下，位移電流的影響更為顯著.針對這一問題，本文編寫了起伏地形下的全電流2D RMT有限元程序.通過非結(jié)構(gòu)的三角形單元實現(xiàn)了任意復雜模型的離散化，同時，程序采用局部加密策略獲得了高精度的數(shù)值解.

由于位移電流的存在，空氣層中磁場的偏導數(shù)不能近似為0，因而在對TM模式進行有限元離散時也必須考慮空氣層.通過Dike模型討論了空氣層厚度對數(shù)值解的影響，當空氣層厚度大于1/4波長即可保證TE/TM模式的有限元解的誤差均在1%以內(nèi).

根據(jù)電流公式，位移電流隨頻率的升高而增大，傳導電流隨介質(zhì)電阻率的增大而減小，因此在高頻高阻情況下，位移電流在總電流中所占的比重就越大.文中結(jié)果表明，考慮位移電流情況下得到的視電阻率和相位曲線相比于準靜態(tài)條件下的計算結(jié)果會出現(xiàn)下掉現(xiàn)象，且頻率越高下掉越明顯，這是因為位移電流項的引入增大了總電流密度，而總電流密度的增大會使得地下介質(zhì)的電導率特性放大，從而導致視電阻率曲線下掉.數(shù)值算例表明，當背景電阻率為10000 Ωm時，頻率大于100 kHz后，位移電流對視電阻率的影響不可忽略，而對相位曲線的影響從頻率為20 kHz后就必須考慮.在背景電阻率與頻率不變的情況下，地形起伏也會影響位移電流，地形高程越大，考慮位移電流計算的視電阻率與準靜態(tài)條件下得到的視電阻率差異也越大，同時，數(shù)值結(jié)果表明差異最大的地方位于山脊的兩側(cè).

通過銅陵舒家店礦床的實際地質(zhì)模型，分別計算了傳導電流和位移電流，結(jié)果表明在高阻花崗閃長巖地區(qū)位移電流不可忽略，會對模擬、觀測結(jié)果造成顯著的影響.如忽略位移電流，采用現(xiàn)有的MT程序進行RMT反演解釋，勢必帶來結(jié)果的偏差，進而影響資料解釋的準確度.而采用本文提供的全電流起伏地形模擬程序，可以獲得更為準確的模擬結(jié)果.

基于本文的研究成果，作者將進一步開發(fā)針對RMT的反演程序.

致謝 感謝Shewchuk提供的非結(jié)構(gòu)三角形網(wǎng)格剖分開源代碼Triangle，感謝Schenk提供的方程組求解器PARDISO，感謝網(wǎng)格剖分可視化軟件PARAVIEW的作者，感謝安徽省地質(zhì)調(diào)查院提供舒家店礦床地質(zhì)模型.同時感謝舒家店野外物性測量人員王顯瑩、張弛等.

Bastani M, Persson L, Beiki M, et al. 2013. A radio magnetotelluric study to evaluate the extents of a limestone quarry in Estonia.GeophysicalProspecting, 61（3）: 678-687.

Candansayar M E, Tezkan B. 2006. A comparison of different radiomagnetotelluric data inversion methods for buried waste sites.JournalofAppliedGeophysics, 58（3）: 218-231. Candansayar M E, Tezkan B. 2008. Two-dimensional joint inversion of radiomagnetotelluric and direct current resistivity data.GeophysicalProspecting, 56（5）: 737-749.

Chave A D, Jones A G. 2012. The Magnetotelluric Method: Theory and Practice. Cambridge: Cambridge University Press. Dong H, Wei W B, Ye G F, et al. 2014. Study of Three-dimensional magnetotelluric inversion including surface topography based on finite-difference method.ChineseJ.Geophys. （in Chinese）, 57（3）: 939-952, doi: 10.6038/cjg20140323.

Franke A, B?rner R U, Spitzer K. 2007. Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography.GeophysicalJournalInternational, 171（1）: 71-86.

Hu X Y, Li Y, Yang W C, et al. 2012. Three-dimensional magnetotelluric parallel inversion algorithm using data space method.ChineseJ.Geophys. （in Chinese）, 55（12）: 3969-3978, doi: 10.6038/j.issn.0001-5733.2012.12.009.

Hu Z Z, Hu X Y, He Z X. 2006. Pseudo-Three-Dimensional magnetotelluric inversion using nonlinear conjugate gradients.ChineseJ.Geophys. （in Chinese）, 49（4）: 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

Intel. 2010. Intel Math Kernel Library. Linear Solver Basics.Ismail N, Schwarz G, Pedersen L B. 2011a. Investigation of groundwater resources using controlled-source radio magnetotellurics （CSRMT） in glacial deposits in Heby, Sweden.JournalofAppliedGeophysics, 73（1）: 74-83.

Ismail N, Pedersen L. 2011b. The electrical conductivity distribution of the Hallands?s Horst, Sweden: a controlled source radiomagnetotelluric study.NearSurfaceGeophysics, 9（1）: 45-54.

Kalscheuer T, Pedersen L B, Siripunvaraporn W. 2008. Radiomagnetotelluric two-dimensional forward and inverse modelling accounting for displacement currents.GeophysicalJournalInternational, 175（2）: 486-514.

Li Y G, Spitzer K. 2002. Three-dimensional DC resistivity forward modelling using finite elements in comparison with finite-difference solutions.GeophysicalJournalInternational, 151（3）: 924-934.

Linde N, Pedersen L B. 2004. Characterization of a fractured granite using radio magnetotelluric （RMT） data.Geophysics, 69（5）: 1155-1165.

Liu C S, Tang J T, Ren Z Y, et al. 2010. Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes.JournalofCentralSouthUniversity（ScienceandTechnology） （in Chinese）, 41（5）: 1855-1859.

Mitsuhata Y, Uchida T. 2004. 3D magnetotelluric modeling using the T-finite-element method.Geophysics, 69（1）: 108-119.

Newman G A, Recher S, Tezkan B, et al. 2003. 3D inversion of a scalar radio magnetotelluric field data set.Geophysics, 68（3）: 791-802. Pedersen L B, Bastani M, Dynesius L. 2005. Groundwater exploration using combined controlled-source and radiomagnetotelluric techniques.Geophysics, 70（1）: G8-G15. Pedersen L B, Bastani M, Dynesius L. 2006. Some characteristics of the electromagnetic field from radio transmitters in Europe.Geophysics, 71（6）: G279-G284.

Persson L, Pedersen L B. 2002. The importance of displacement currents in RMT measurements in high resistivity environments.JournalofAppliedGeophysics, 51（1）: 11-20.

Ren Z Y, Tang J T. 2010. 3D direct current resistivity modeling

with unstructured mesh by adaptive finite-element method.Geophysics, 75（1）: H7-H17.

Ren Z Y, Kalscheuer T, Greenhalgh S, et al. 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling.GeophysicalJournalInternational, 194（2）: 700-718.

Schenk O, G?rtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO.FutureGenerationComputerSystems, 20（3）: 475-487.

Shewchuk J R. 1996. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator. ∥Lin M C, Manocha D eds. Applied Computational Geometry Towards Geometric Engineering. Berlin Heidelberg: Springer, 1148: 203-222.

Tan H D, Yu Q F, Booker J, et al. 2003. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method.ChineseJ.Geophys. （in Chinese）, 46（5）: 705-711, doi: 10.3321/j.issn:0001-5733.2003.05.019. Tezkan B, H?rdt A, Gobashy M. 2000. Two-dimensional radiomagnetotelluric investigation of industrial and domestic waste sites in Germany.JournalofAppliedGeophysics, 44（2-3）: 237-256.

Tezkan B, Georgescu P, Fauzi U. 2005. A radiomagnetotelluric survey on an oil-contaminated area near the Brazi Refinery, Romania.GeophysicalProspecting, 53（3）: 311-323.

Tezkan B, Saraev A. 2008. A new broadband radiomagnetotelluric instrument: applications to near surface investigations.NearSurfaceGeophysics, 6（4）: 245-252. Wannamaker P E, Stodt J A, Rijo L. 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements.Geophysics, 51（11）: 2131-2144.

Wu X P, Liu Y, Wang W. 2015. 3D resistivity inversion incorporating topography based on unstructured meshes.ChineseJ.Geophys. （in Chinese）, 58（8）: 2706-2717, doi: 10.6038/cjg20150808.

附中文參考文獻

董浩, 魏文博, 葉高峰等. 2014. 基于有限差分正演的帶地形三維大地電磁反演方法. 地球物理學報, 57（3）: 939-952, doi: 10.6038/cjg20140323.

胡祥云, 李焱, 楊文采等. 2012. 大地電磁三維數(shù)據(jù)空間反演并行算法研究. 地球物理學報, 55（12）: 3969-3978, doi: 10.6038/j.issn.0001-5733.2012.12.009.

胡祖志, 胡祥云, 何展翔. 2006. 大地電磁非線性共軛梯度擬三維反演. 地球物理學報, 49（4）: 1226-1234, doi: 10.3321/j.issn:0001-5733.2006.04.039.

劉長生, 湯井田, 任政勇等. 2010. 基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬. 中南大學學報（自然科學版）, 41（5）: 1855-1859.

譚捍東, 余欽范, Booker J等. 2003. 大地電磁法三維交錯采樣有限差分數(shù)值模擬. 地球物理學報, 46（5）: 705-711, doi: 10.3321/j.issn:0001-5733.2003.05.019. 吳小平, 劉洋, 王威. 2015. 基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演. 地球物理學報, 58（8）: 2706-2717, doi: 10.6038/cjg20150808.

（本文編輯 何燕）

Two-dimensional complicated radio-magnetotelluric finite-element modeling using unstructured grids

YUAN Yuan1,2， TANG Jing-Tian1,2， REN Zheng-Yong1,2*， ZHOU Cong1,2， XIAO Xiao1,2， ZHANG Lin-Cheng1,2

1KeyLaboratoryofMetallogenicPredictionofNonferrousMetalsandGeologicalEnvironmentMonitoring,MinistryofEducation，CentralSouthUniversity，Changsha410083，China2InstituteofAppliedGeophysics,SchoolofGeosciencesandInfo-Physics，CentralSouthUniversity，Changsha410083，China

As a newly developed geophysical exploration method, the radio-magnetotelluric （RMT） method is widely used in near-surface engineering and environment geophysical investigation. Since the interpretation of RMT data usually adopts the magnetotelluric （MT） forward modeling routine, in which displacement currents are generally neglected, the inverted conductivity model may not correctly reflect the true conductivity structure in the shallow subsurface. To solve this issue, we developed a finite-element forward modeling code for RMT data, in which displacement currents are considered. With this code, we analyze the dielectric effect of displacement currents on RMT responses over resistive subsurface models.First, we derived the boundary value problem （BVP） about the EM field components by the Galerkin method, in which the displacement currents were considered. Then we used the finite element method and PARDISO solver to calculate the electric and magnetic field components. To deal with complicated structure and surface topography, unstructured triangle grids were adopted for mesh generation. To improve the computation accuracy, the local refinement was used. At last, we developed a forward modeling code for RMT data with the consideration of displacement currents and analyzed the effect of displacement currents on 2D TM-mode, TE-mode data, which measured with the RMT method at frequencies between 10 and 300 kHz.First, a synthetic model was used to verify the correction of our new code. The result shows that the response calculated by our code agrees well with other results. Utilizing the Dike model, the effect of the thickness of the air layer on accuracy of numerical solutions was investigated. The result shows that when the thickness of the air layer is greater than 1/4 wavelength, highly accurate solutions can be guaranteed. Then the impact of displacement currents on RMT data with ridge terrain was studied on a hill model with complicated topography. From this model, the following results can be demonstrated: （1） The effect of displacement currents would increase with increasing height of the hill and the corner of hill was more easier to be affected. （2） The phase curves are more likely to be affected than apparent resistivity curves at high frequency. （3） The effect of displacement currents on apparent resistivity cannot be neglected when the frequency is larger than 100 kHz and the effect on phase must be considered when the frequency is larger than 20 kHz. Finally, a field model was studied to further demonstrate the importance of displacement currents in the RMT method. The results show that: （1） The error caused by displacement currents increases with frequency. （2） The apparent resistivity of TM-mode is more easily to be affected by displacement currents than TE-mode apparent resistivity. （3） For the area of quartz diorite with high resistivity, the percentage of displacement current density in all current density can be 10%. It is clear that displacement currents must be considered in RMT forward modeling.With numerical calculations, we demonstrated the effect of displacement currents on 2D RMT data. Forward modeling confirms that responses computed in the quasi-static approximation become increasingly inaccurate with rising frequency and strongly affected by terrain. However, RMT data processing and interpretation are mostly based on the MT program in recent years, which may result in incorrect structure. Based on the work in this paper, the author will develop RMT inversion code with consideration of displacement current.

Magnetotelluric; RMT; Finite element method; Unstructured grids; Displacement currents

國家深部探測專項第3項目（SinoProbe-03），“十二五”國家科技支撐計劃課題（2011BAB04B01），國家自然科學基金（41574120），國家高技術(shù)研究發(fā)展計劃（863計劃）（2014AA06A602），國家重點基礎(chǔ)研究發(fā)展計劃（973計劃）（2015CB060201），中南大學創(chuàng)新驅(qū)動計劃（2015CX008）資助.

原源，女，1989年生，博士研究生，從事電磁法有限元正反演研究. E-mail: 13657400521@163.com

*通信作者 任政勇，男，1983年生，副教授，主要從事電磁場數(shù)值模擬及反演研究. E-mail: renzhengyong@csu.edu.cn

10.6038/cjg20151229.

10.6038/cjg20151229

P631

2014-10-27，2015-06-18收修定稿

原源， 湯井田， 任政勇等. 2015. 基于非結(jié)構(gòu)化網(wǎng)格的任意復雜2D RMT有限元模擬.地球物理學報，58（12）:4685-4695,

Yuan Y, Tang J T, Ren Z Y, et al. 2015. Two-dimensional complicated radio-magnetotelluric finite-element modeling using unstructured grids.ChineseJ.Geophys. （in Chinese）,58（12）:4685-4695,doi:10.6038/cjg20151229.