湯井田，薛 帥

中南大學地球科學與信息物理學院，長沙 410083

0 引言

國內(nèi)外學者很早就進行了大地電磁場的數(shù)值模擬研究。1971年Coggon［1］首先提出將有限單元法應(yīng)用在電磁法數(shù)值模擬中，并詳細闡述了模擬大地電磁場的有限元格式。有限單元法作為一種便捷通用的數(shù)值分析方法，相較于有限差分法和積分方程法，具有易于擬合復雜結(jié)構(gòu)、易于對復雜介質(zhì)目標建模、程序通用性強等優(yōu)勢，從而被迅速廣泛地應(yīng)用于電磁場的數(shù)值模擬中。國外學者Rijo［2］、Wannamaker等［3－4］分別對電磁場的有限元模擬進行了發(fā)展和完善，Wannamaker等［3－4］將有限元中的矩形網(wǎng)格繼續(xù)剖分為4個三角形網(wǎng)格，可以更好地模擬地形和地下構(gòu)造，提高了數(shù)值計算精度。在國內(nèi)，1981年，陳樂壽［4］首先詳細介紹了有限元在大地電磁場正演計算中的應(yīng)用，并利用矩形－三角形網(wǎng)格剖分對有限元模擬大地電磁進行了改進。之后胡建德等［6］、徐世浙等［7］、楊長福［8］、陳小斌［9］分別對地球物理場中的有限單元法進行了研究和發(fā)展，至今，有限單元法已經(jīng)成為二維和三維大地電磁數(shù)值模擬的主要手段。

在利用有限元方法模擬二維和三維大地電磁場時，電磁場正演計算時的一些邊界條件需要在“無窮遠處”才能夠滿足。而有限元的模擬區(qū)域是有限的，所以網(wǎng)格邊界處存在截斷邊界影響。在有限元模擬大地電磁場時，為了避免截斷邊界的影響，一般將網(wǎng)格邊界放在距離研究區(qū)域外很遠處，如 Rijo［2］、Wannamaker［3－4］、陳樂壽等［10］、徐世浙［11］、趙廣茂等［12］、劉長生等［13］在進行大地電磁有限元數(shù)值模擬時網(wǎng)格邊界都放在足夠遠處來滿足邊界條件。但是“足夠遠距離”一直沒有適當?shù)膮⒖贾?，網(wǎng)格邊界選擇小了會影響計算精度，而選擇過大則需要較高的內(nèi)存和付出較長的計算時間。Sharma等［14］曾利用趨膚深度δ來進行網(wǎng)格剖分和劃分網(wǎng)格邊界，認為將網(wǎng)格邊界設(shè)置在5δ處時即可忽略截斷邊界的影響，但是沒有進行相應(yīng)的研究說明。作者利用有限元二維大地電磁正演程序，計算總結(jié)出適合有限元大地電磁模擬的參考網(wǎng)格邊界。

1 邊界條件和截斷邊界

假定地下電性結(jié)構(gòu)為二維結(jié)構(gòu)，選取右旋直角坐標系的原點在地面上，z軸正方向垂直向下，x軸的方向平行于走向方向，即介質(zhì)的電性參數(shù)僅與y、z方向有關(guān)。根據(jù)麥克斯韋方程，二維地電結(jié)構(gòu)下的大地電磁滿足以下二式，稱之為亥娒霍茲方程式［10］：

其中：Hx和Ex是x 方向磁場和電場；k2＝iωμσ；ω是角頻率；μ是導磁率；σ是導電率。

圖1 E型波（a）和 H型波（b）的研究區(qū)域［11］Fig.1 Research area of E type wave（a）and H type wave（b）［11］

邊界如圖1所示，設(shè)u為需要求解的電磁場，在上邊界AB處，電磁場滿足第一類邊界條件；值得注意的是，TE（transverse electric）模式的上邊界要離開地面一段距離。在左右邊界AD和BC處，要求由橫向不均勻體引起的二次場為0。為此，有2種方式的邊界條件加載［10］：一種是強加邊界條件，即計算出水平層狀介質(zhì)中的電磁場分布來給出邊界上的值；另一種是采用齊次的第二類邊界條件，即在側(cè)邊界上的場分量u的法向?qū)?shù)為零。一般采用后一種方式作為側(cè)邊界的邊界條件。而下邊界CD處以下被認為是均質(zhì)巖石，要求局部異常體在CD上的異常場為0。所以u滿足的邊界條件為

當?shù)叵麓嬖诓痪鶆蝮w時，網(wǎng)格邊界需要在足夠遠處才能夠滿足電磁場的邊界條件。

2 截斷邊界研究

由于在利用大地電磁（magnetotellurics，MT）有限元正演程序進行大地電磁數(shù)值模擬時，網(wǎng)格剖分會對計算結(jié)果產(chǎn)生很大的影響，所以在進行截斷邊界研究時，要先避免因網(wǎng)格剖分而造成的計算結(jié)果差異。筆者為確保計算結(jié)果的正確與可靠，在橫向上（y方向）－1 000～1 000m、縱向上（z方向）0～2 000m，按照10m的網(wǎng)格間距進行剖分；而在該研究區(qū)域以外，無論橫向還是縱向的網(wǎng)格間距，皆按1.1倍率進行增長剖分。對于TE模式參照Wannamaker等［3－4］在PW2D中的作法，以總網(wǎng)格橫向距離的一半作為上網(wǎng)格邊界的距離。

數(shù)值模擬時，在不改變網(wǎng)格剖分的情況下，只對外網(wǎng)格邊界進行擴大和減小，從而避免網(wǎng)格剖分的影響。根據(jù)電磁場理論，某一頻率的電磁場在該頻點的一個趨膚深度內(nèi)迅速衰減至1／e，并在傳播中繼續(xù)衰減至0；因此，筆者以某一頻點的趨膚深度作為該頻點網(wǎng)格邊界選擇的量度：

式中：δ為趨膚深度；ρ為電阻率；f為頻率。可見，δ與f相關(guān)。為了體現(xiàn)普遍性，在MT儀器常用的4個頻段中分別選擇一個頻率作為研究頻率，即0.25、2.50、25.00、250.00Hz。

在模型選擇上，筆者分別對一維模型（均勻半空間模型和層狀介質(zhì)模型）和二維地電模型進行模擬計算。在進行數(shù)值模擬前，將2個方向的網(wǎng)格邊界分別放在20倍的趨膚深度距離處，計算出一組數(shù)據(jù)作為參考解，認為該參考解是在“足夠遠”的網(wǎng)格邊界條件下計算的較精確的數(shù)值解。

2.1 一維模型模擬

首先對電阻率為100Ω·m的均勻半空間模型進行模擬計算，分別計算左右網(wǎng)格邊界和下網(wǎng)格邊界同時為100m和同時為500m的情況，計算結(jié)果如表1。

再建立一個層狀介質(zhì)模型（圖2），分別計算左右網(wǎng)格邊界和下網(wǎng)格邊界同時為500m和同時為1 000m情況下的視電阻率，結(jié)果如表2所示。

表1 均勻半空間視電阻率計算結(jié)果Table1 Apparent resistivity results of homogenous half space

表2 H型層狀介質(zhì)視電阻率計算結(jié)果Table2 Apparent resistivity results of layered medium model

在一維模型中不存在局部異常體，自然滿足有限元模擬的外邊界條件。通過對以上2個模型的計算，進一步說明了在一維有限元數(shù)值模擬時邊界條件自然滿足，不存在截斷邊界影響，從而無需進行網(wǎng)格邊界的擴展。

圖2 層狀介質(zhì)模型Fig.2 Layered medium model

2.2 二維地電模型

有限元模擬地下存在二度體的二維地電模型時，由于局部不均勻體的存在，不合適的網(wǎng)格邊界會產(chǎn)生較大的誤差。對模型2D－1（圖3）進行模擬計算，設(shè)置左右網(wǎng)格邊界和下網(wǎng)格邊界均為500m，分別在O點和B點進行觀測，表3是O點觀測結(jié)果。

從表3中可以看出，不合適的網(wǎng)格邊界對計算結(jié)果影響較大。因此，對二維地電模型進行有限元模擬時，需要選擇合適的網(wǎng)格邊界。

圖3 二維地電模型Fig.3 Two－dimensional geoeletric model

根據(jù)所選頻率和模型電阻率，將下網(wǎng)格邊界設(shè)置在“足夠遠”的距離，放在200 000m的地方，計算趨膚深度作為變化量度來改變左右網(wǎng)格邊界的邊界距離，將計算結(jié)果與參考解進行比較。同時約定下網(wǎng)格邊界以地下異常體的下邊界為基準進行擴展，左右網(wǎng)格邊界以異常體的左右邊界為基準進行擴展。約定計算趨膚深度時的電阻率為地下電性結(jié)構(gòu)中電阻率最高者。

對背景介質(zhì)中存在低阻異常體和高阻異常體2種模型進行模擬計算，并在二度體中心上方O點和邊界上方B點同時進行TE模式和TM（transverse magnetic）模式觀測。

2.2.1 左右網(wǎng)格邊界變化

表3 模型2D－1在O點觀測結(jié)果Table3 Observed results of model 2D－1on point O

表4 左右網(wǎng)格邊界變化下視電阻率觀測結(jié)果Table4 Apparent resistivity results when changing left and right mesh boundary

從表4可看出：對于模型2D－1，TE模式在1δ網(wǎng)格邊界時，O點和B點觀測結(jié)果的相對誤差除250Hz外大部降至1%以下；在2δ的網(wǎng)格邊界時，O點觀測結(jié)果的相對誤差已全降至0.1%以下，B點雖相對于1δ網(wǎng)格邊界條件時的相對誤差有所減小，但部分頻率的誤差仍在0.1%數(shù)量級上；繼續(xù)增大網(wǎng)格邊界至3δ時，O點和B點觀測的相對誤差已全部降至0.1%以下；繼續(xù)加大網(wǎng)格邊界至5δ和10δ，觀測結(jié)果與參考結(jié)果相等，誤差為0。

對于TM模式：在1δ網(wǎng)格邊界時，計算的結(jié)果相對誤差全部在1%以下；增大到2δ網(wǎng)格邊界時，相對誤差大部為0；繼續(xù)加大網(wǎng)格邊界，相對誤差保持為0。對于模型2D－2，在1δ網(wǎng)格邊界時，大部分計算結(jié)果與參考解相等，相對誤差為0。

值得注意的是，對模型2D－2進行模擬時，網(wǎng)格邊界變化的量度采用的趨膚深度已相當于模型2D－1的3倍。

2.2.2 下網(wǎng)格邊界變化

在對下邊界研究時，則把左右網(wǎng)格邊界設(shè)置在“足夠遠”的距離，以趨膚深度作為量度來改變下網(wǎng)格邊界的邊界距離，計算結(jié)果并與參考解進行比較。仍然以模型2D－1和模型2D－2為計算模型。模型2D－1的部分計算結(jié)果如表5。

從表5中看出：根據(jù)TE模式計算結(jié)果，在1δ網(wǎng)格邊界時，O點和B點觀測結(jié)果的相對誤差已全部降至1%以下；在2δ網(wǎng)格邊界時，相對誤差繼續(xù)減小，但大部仍在0.1%的數(shù)量級上；繼續(xù)增大網(wǎng)格邊界至3δ時，O點和B點觀測的相對誤差已全部降至0.1%以下；加大邊界至5δ和10δ，觀測結(jié)果與參考結(jié)果相等，誤差為0。

對于TM模式：在1δ網(wǎng)格邊界時，計算的結(jié)果相對誤差大部在1%以下；增大到2δ網(wǎng)格邊界時，除B點觀測的25Hz和250Hz外，大部分結(jié)果的相對誤差降至0.1%以下；繼續(xù)增大網(wǎng)格邊界，誤差繼續(xù)減少直至為0。

對于下網(wǎng)格邊界變化下的高阻模型2D－2，計算的結(jié)果與對應(yīng)的左右網(wǎng)格邊界的變化結(jié)果相似。即在1δ時，大部分計算結(jié)果已與參考解相等，相對誤差為0。

3 截斷邊界選擇

3.1 綜合網(wǎng)格邊界計算

以上的計算結(jié)論都是在某一方向的網(wǎng)格邊界位于“足夠遠處”得到的，現(xiàn)在把以上得出的最佳網(wǎng)格邊界距離綜合在一起進行計算比較。

對于模型2D－1：計算TE模式時，左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處；計算TM模式時，則將左右網(wǎng)格邊界和下網(wǎng)格邊界放在2δ處。對于模型2D－2，把網(wǎng)格邊界全部設(shè)置在1δ處即可。綜合網(wǎng)格邊界條件下部分觀測結(jié)果如表6所示。

由表6可見，以上利用綜合網(wǎng)格邊界計算的結(jié)果相對誤差基本在0.1%以下，說明該網(wǎng)格邊界條件下計算結(jié)果的可靠性。

3.2 綜合模型計算

為了進一步驗證綜合的最佳網(wǎng)格邊界計算結(jié)果的可靠性，建立1個較復雜的二維模型，并采用新的頻率進行計算。模型2D－3（圖4）是在背景介質(zhì)中賦存3種異常體的較復雜模型，分別在O、O1、O2和B點觀測，計算的頻率為0.01、0.10、1.00、10.00、100.00、1 000.00Hz，剖分情況與前面相同。

在模型2D－3中存在多個異常體，在確定網(wǎng)格邊界變化量度時，需要幾個邊界變化量度進行比較才能確定。在TE模式計算時，對于低阻體ρ2和ρ4，當它們賦存在較高電阻率的背景介質(zhì)中時，趨膚深度計算采用背景電阻率，邊界變化量度為3δ＝503；而當背景介質(zhì)中存在高阻體ρ3時，趨膚深度計算采用高阻體電阻率ρ3，其邊界變化量度為δ ＝503-。為了確保計算結(jié)果正確和保證精度，網(wǎng)格邊界變化量度應(yīng)為TM模式計算時網(wǎng)格邊界變化量度的確定與之類似。表7為模型2D－3的部分計算結(jié)果。

表5 模型2D－1下網(wǎng)格邊界變化下O點視電阻率觀測結(jié)果Table5 Apparent resistivity results of model 2D－1on point Owhen changing bottom boundary

表6 綜合網(wǎng)格邊界O點觀測結(jié)果Table6 Observed results on point Owith compositive boundaries

表7 模型2D－3在O點的觀測結(jié)果Table7 Observed results of model 2D－3on point O

圖4 模型2D－3Fig.4 Model 2D－3

由表7可見，利用綜合的網(wǎng)格邊界計算的結(jié)果與參考解的相對誤差全部在0.1%以下，進一步驗證了綜合邊界的正確性。

值得注意的是，以上的結(jié)果是在背景電阻率和異常體電阻率最大相差一個數(shù)量級的情況下得出的。作者又利用綜合邊界計算了背景電阻率為100 Ω·m、異常體電阻率分別為1、0.1、0.01、0.001 Ω·m的情況以及背景電阻率為100Ω·m、異常體電阻率為200Ω·m的情況，計算結(jié)果表明TE模式中誤差最大為0.40%，TM模式誤差基本為0。另外作者還計算了一些極端的情況，如背景電阻率為1Ω·m、異常體電阻率為2Ω·m的模型和背景電阻率為2Ω·m、異常體電阻率為1Ω·m的模型，發(fā)現(xiàn)參考邊界在此也適用。

4 結(jié)論與建議

利用有限元方法對大地電磁進行模擬計算時，由于存在截斷邊界影響，從而影響計算結(jié)果和精度，或占用大量內(nèi)存花費較長的時間。筆者通過大量的數(shù)值模擬計算與對比，以趨膚深度作為網(wǎng)格邊界變化量度，考慮到大地電磁的精度要求，總結(jié)得出適合二維大地電磁有限元正演模擬計算的參考網(wǎng)格邊界：

1）利用有限元對一維地電模型進行模擬計算時，不存在截斷邊界的影響，邊界條件自然滿足。

2）對二維低阻模型進行模擬計算時，以背景電阻率計算趨膚深度作為網(wǎng)格邊界變化量度。對于TE模式，當左右網(wǎng)格邊界和下網(wǎng)格邊界都放在3δ處時，截斷邊界的影響已可忽略。對于TM模式，左右網(wǎng)格邊界和下網(wǎng)格邊界在2δ處，截斷邊界影響可以忽略。

3）對二維高阻模型進行模擬計算時，以高阻體電阻率計算趨膚深度作為網(wǎng)格邊界變化量度。不管TE模式還是TM模式，當左右網(wǎng)格邊界和下網(wǎng)格邊界位于1δ時，截斷邊界影響已可忽略。

4）當二維地電模型中，同時存在低阻異常體和高阻異常體時，分別以背景介質(zhì)電阻率和高阻體中電阻率最高者計算趨膚深度，然后根據(jù)參考網(wǎng)格邊界選擇最佳的網(wǎng)格邊界。

（References）：

［1］Coggon J H.Electromagnetic and Electrical Modeling by Finite Element Method［J］.Geophysics，1971，36（1）：132－155.

［2］Rodi W L.A Technique for Improving the Accuracy of Finite Element Solution for Magnetotelluric Data［J］.Geophysics，1976，44（2）：483－506.

［3］Wannamaker P E，Stodt J A，Rijo W L.Two－Dimensional Topographic Responses in Magnetotelluric Modeling Using Finite Element［J］.Geophysics，1986，51（11）：2131－2144.

［4］Wannamaker P E，Stodt J A，Rijo W L.A Stable Finite Element Solution for Two－Dimensional Magnetotelluric Modelling［J］.Geophysics，1987，88（1）：277－296.

［5］陳樂壽.有限元法在大地電磁場正演計算中的應(yīng)用與改進［J］.石油物探，1981，20（3）：84－103.Chen Leshou.The Improvement and Application of Finite Element Method to 2DForward Solution of Magnetotelluric Sounding［J］.Geophysical Prospecting for Petroleum，1981，20（3）：84－103.

［6］胡建德，蔡剛.用三角形二次插值有限元法計算二維大地電磁測深曲線［J］.石油地球物理勘探，1984，19（4）：358－367.Hu Jiande，Cai Gang.Calculating 2DMagnetotelluric Sounding Curve with the Twice Interpolation Triangle Finite Element Method ［J］. Oil Geophysical Prospecting，1984，19（4）：358－367.

［7］史明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學報，1997，40（3）：421－430.Shi Mingjuan，Xu Shizhe，Liu Bin.Finite Element Method Using Quadratic Element in MT Forward Modeling［J］.Acta Geophysica Sinica，1997，40（3）：421－430.

［8］楊長福.大地電磁二維對稱各向異性介質(zhì)的有限元數(shù)值模擬［J］.西北地震學報，1997，19（2）：27－33.Yang Changfu.MT Numerical Simulation of Symmetrically 2DAnisotropic Media Based on the Finite Element Method ［J］. Northwestern Seismological Journal，1997，19（2）：27－33.

［9］陳小斌.有限元直接迭代算法［J］.物探化探計算技術(shù)，1999，21（2）：165－171.Chen Xiaobin.Direct Iterative Algorithm of Finite Element［J］.Computing Techniques for Geophysical and Geochemical Exploration，1999，21（2）：165－171.

［10］陳樂壽，劉任，王天生.大地電磁測深資料處理與解釋［M］.北京：石油工業(yè)出版社，1989：1－160.Chen Leshou， Liu Ren， Wang Tiansheng.Magnetotelluric Sounding Data Processing and Interpretation［M］.Beijing：Petroleum Industry Press，1989：1－160.

［11］徐世浙.地球物理中的有限單元法［M］.北京：科學出版社，1994：220－241.Xu Shizhe.The Finite Element in Geophysics［M］.Beijing：Science Press，1994：220－241.

［12］趙廣茂，李桐林，王大勇，等.基于二次場二維起伏地形MT有限元數(shù)值模擬［J］.吉林大學學報：地球科學版，2008，38（6）：1055－1059.Zhao Guangmao，Li Tonglin，Wang Dayong，et al.Secondary Field－Based Two－Dimensional Topographic Numerical Simulation in Magnetoteilurics by Finite Element Method［J］.Journal of Jilin University：Earth Science Edition，2008，38（6）：1055－1059.

［13］劉長生，湯井田，任政勇，等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬［J］.中南大學學報：自然科學版，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.

［14］Sharma P S，Kaikkonen P.An Automated Finite Element Mesh Generation and Element Coding in 2－D Electromagnetic Inversion［J］.Geophysics，1998，34（3）：93－114.