有限單元法與有限差分法在MT一維正演模擬中的對(duì)比分析

王 濤，柳建新，童孝忠＊，曹創(chuàng)華，譚神湘

（1.中國(guó)海洋大學(xué) 海洋地球科學(xué)學(xué)院，青島 266100；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083）

0 前言

大地電磁測(cè)深法（Magnetotelluric Sounding，簡(jiǎn)稱MT）是一種重要的地球物理勘探方法，由蘇聯(lián)學(xué)者 Tikhonov［1］和Cagniard［2］于上世紀(jì)五十年代初期提出。與其它地球物理勘探方法一樣，正演問題是大地電磁測(cè)深法的理論基礎(chǔ)，同時(shí)也是我們認(rèn)識(shí)各種地電條件下大地電磁場(chǎng)響應(yīng)特征的良好途徑，其研究一直受到廣泛關(guān)注。通過對(duì)不同地質(zhì)模型的正演研究，我們總結(jié)出在不同地質(zhì)條件下大地電磁場(chǎng)的分布規(guī)律；同時(shí)通過正演，也能準(zhǔn)確了解在不同的地形起伏情況下，大地電磁場(chǎng)的分布特點(diǎn)。

大地電磁測(cè)深正演模擬的數(shù)值方法主要分為三種：①有限差分法；②積分方程法；③有限單元法。作者主要針對(duì)有限單元法和有限差分法進(jìn)行討論。有限單元法（Finite Element Method或簡(jiǎn)稱FEM）是將要分析的連續(xù)場(chǎng)分割為很多較小的區(qū)域。Coggon［3］首先將有限單元法應(yīng)用在電磁法正演模擬中；七十年代末，朱伯芳［4］將有限單元法引入國(guó)內(nèi)。有限差分法（Finite Difference Method或簡(jiǎn)稱FDM）是一種經(jīng)典的數(shù)值模擬方法。地球物理工作者從八十年代開始研究有限差分法正演計(jì)算問題，為有限差分法的發(fā)展做了相當(dāng)大的貢獻(xiàn)（周熙襄等［5］；羅延鐘等［6－7］。根據(jù)長(zhǎng)期的研究，有限單元法和有限差分法在實(shí)際應(yīng)用中均有各自的優(yōu)勢(shì)和不足，作者在本文中主要對(duì)有限單元法和有限差分法進(jìn)行對(duì)比分析。

1 數(shù)值模擬

1.1 大地電磁測(cè)深正演的基本理論

麥克斯韋方程組是電磁場(chǎng)必須遵從的微分方程組，利用傅里葉變換可將任意隨時(shí)間變化的電磁場(chǎng)分解為一系列諧變場(chǎng)的組合，通常我們以e－iωt表示諧變場(chǎng)的時(shí)間因子（即以負(fù)諧時(shí)表示），電場(chǎng)強(qiáng)度和磁場(chǎng)強(qiáng)度可表示為式（1）與式（2）。

大地電磁測(cè)深所討論的電磁場(chǎng)頻率是極低的，故在大地介質(zhì)中可忽略位移電流對(duì)場(chǎng)分布的影響，即在大地電磁測(cè)深正演中研究的是似穩(wěn)電磁場(chǎng)問題。于是，導(dǎo)電介質(zhì)低頻諧變場(chǎng)的麥克斯韋方程組為式（3）～式（6）。

式中 ▽·E＝0是因?yàn)閷?dǎo)電介質(zhì)內(nèi)部體電荷密度實(shí)際上為零，公式中時(shí)間因子都隱含在電場(chǎng)E和磁場(chǎng)H 中，方程組（3）～方程組（6）是大地電磁測(cè)深理論研究的出發(fā)點(diǎn)。

1.2 一維有限單元正演算法

電場(chǎng)滿足的微分方程為公式（7）。

根據(jù)邊值問題所對(duì)應(yīng)的變分問題為

采用有限單元法進(jìn)行計(jì)算，首先應(yīng)將區(qū)域剖分為若干單元。在單元內(nèi)電導(dǎo)率（或電阻率）必須是連續(xù)的，也就是說，電導(dǎo)率的間斷點(diǎn)不能在單元內(nèi)。對(duì)于高頻，電磁波衰減非常迅速，如果采用線性插值，必須取很小的單元，從而增加了計(jì)算工作量。

（1）單元積分為式（9）。

（2）單元積分為式（10）。

然后將各單元的擴(kuò)展矩陣相加，則得K＝∑K1e－∑K2e＋K3。再根據(jù)δF（Ex）＝0，可得

將電場(chǎng)所滿足的上邊界條件和下邊界條件代入上面的方程組，則有

求解上述線性方程組式（12）即可得到節(jié)點(diǎn)處的電場(chǎng)值，從而也可以進(jìn)一步計(jì)算模型響應(yīng)的視電阻率和相位。

1.3 一維有限差分正演算法

根據(jù)電場(chǎng)所滿足的微分方程，將一維地電模型離散化，如圖1所示。

圖1 地電模型離散化Fig.1 Geoelectric model discretization

然后對(duì)電場(chǎng)進(jìn)行離散處理，采用有限差分近似計(jì)算一階、二階偏導(dǎo)數(shù)。寫成矩陣形式如式（13）。

根據(jù)下邊界條件，電場(chǎng)衰減為“0”，當(dāng)上邊界條件取Hy＝1，則（E1＋1／2－E1－1／2）＝iωμ，得到公式（14）。

求解方程組式（14）可得節(jié)點(diǎn)處電場(chǎng)值，從而可以進(jìn)一步計(jì)算模型響應(yīng)的視電阻率和相位。

2 模型算例

2.1 均勻半空間模型

設(shè)計(jì)一個(gè)均勻半空間模型，ρ＝10Ω·m，分別用有限單元法和有限差分法進(jìn)行正演模擬。圖2和圖3分別顯示有限單元法和有限差分法模擬的視電阻率的相對(duì)誤差。

由圖2和圖3可知，兩幅曲線圖的變化趨勢(shì)基本一致，在低頻段，這兩種方法的計(jì)算結(jié)果與實(shí)際值基本一致；而在高頻段，則出現(xiàn)隨著頻率升高視電阻率的相對(duì)誤差增大的現(xiàn)象。從整體上看，這兩種方法的計(jì)算結(jié)果都基本符合正演的要求，但是有限差分法正演結(jié)果的相對(duì)誤差比有限單元法正演結(jié)果的相對(duì)誤差要小，說明有限差分法在均勻半空間模型的正演效果比有限單元法略好。通過以上模擬的計(jì)算結(jié)果對(duì)比，證明了兩種方法正演算法的正確性，以及所編寫程序的正確性。

2.2 層狀介質(zhì)模型

2.2.1 二層層狀介質(zhì)

設(shè)計(jì)一個(gè)二層層狀介質(zhì)模型，ρ1＝10Ω·m，ρ2＝100Ω·m，h1＝1 000m。

（1）解析法正演結(jié)果。正演計(jì)算的曲線如圖4所示。

（2）有限單元法正演結(jié)果。有限單元法的網(wǎng)格剖分如表1所示，有限單元法計(jì)算結(jié)果與解析法計(jì)算結(jié)果比較如圖5所示。

公式（13）如下：

公式（14）如下：

（3）有限差分法正演結(jié)果。有限差分法的網(wǎng)格剖分如表2所示，有限差分法計(jì)算結(jié)果與解析法計(jì)算結(jié)果比較如圖6所示。

2.2.2 三層層狀介質(zhì)

（1）解析法正演結(jié)果。正演計(jì)算的曲線如圖7所示。

（2）有限單元法正演結(jié)果。有限單元法的網(wǎng)格剖分如表3所示，有限單元法計(jì)算結(jié)果與解析法計(jì)算結(jié)果比較如圖8所示。

（3）有限差分法正演結(jié)果。有限差分法的網(wǎng)格剖分如表4所示，有限差分法計(jì)算結(jié)果與解析法計(jì)算結(jié)果比較如圖9所示。

圖4 二層層狀介質(zhì)大地電磁響應(yīng)曲線Fig.4 Two－tier media magnetotelluric response curve

表1 有限單元法網(wǎng)格剖分參數(shù)Tab.1 Finite element method meshing parameters

表2 有限差分法網(wǎng)格剖分參數(shù)Tab.2 Finite difference method meshing parameters

圖7 三層層狀介質(zhì)大地電磁響應(yīng)曲線Fig.7 Three－tier media magnetotelluric response curve

表3 有限單元法網(wǎng)格剖分參數(shù)Tab.3 Finite element method meshing parameters

表4 有限差分法網(wǎng)格剖分參數(shù)Tab.4 Finite difference method meshing parameters

（1）由圖5和圖6可得，在二層層狀介質(zhì)模型中，當(dāng)頻率較高時(shí)，有限單元法與有限差分法計(jì)算的視電阻率曲線與理論值曲線吻合非常好；而當(dāng)頻率較低時(shí)，兩種方法計(jì)算的視電阻率曲線逐漸地脫離了理論的視電阻率曲線，在尾端出現(xiàn)了上升的趨勢(shì)。

（2）由圖8和圖9可得，在三層層狀介質(zhì)模型中，有限單元法與有限差分法計(jì)算的視電阻率曲線與理論值曲線基本上一致，但二者在曲線的極大值處附近與理論值有所偏離。①有限單元法計(jì)算的視電阻率曲線與理論值曲線對(duì)比的趨勢(shì)為：開始偏離后漸漸吻合，再過極大值點(diǎn)偏離，后漸漸吻合；②有限差分法計(jì)算的視電阻率曲線與理論值曲線對(duì)比的趨勢(shì)為：開始吻合，過極小值點(diǎn)后偏離，再過極大值后吻合。

3 結(jié)論

作者在本文以大地電磁測(cè)深的基本理論為基礎(chǔ)，從麥克斯韋方程組出發(fā)，以有限單元法和有限差分法為計(jì)算工具，研究了大地電磁測(cè)深的一維正演問題。通過對(duì)一些特定的模型進(jìn)行正演模擬，驗(yàn)證了有限單元法和有限差分法在大地電磁測(cè)深正演計(jì)算中的有效性。從以上的二層層狀介質(zhì)模型和三層層狀模型對(duì)比分析中，我們可以看出，有限單元法和有限差分法的計(jì)算結(jié)果與解析法的計(jì)算結(jié)果（即理論值）基本一致，兩種方法的正演結(jié)果均真實(shí)地反映了模型的地電參數(shù)，說明了有限單元正演算法和有限差分正演算法的正確性，同時(shí)也說明了編寫程序的正確性。

［1］TIKNONOV A N.On determining electrical characteristics of the deep layers of the earth’s crust［J］.Deki Akud Nuck，1950（73）：295－297.

［2］CAGNIARD L.Basic theory of the magnetotelluric methods of geophysical prospecting［J］.Geophysics，1953（18）：605－635.

［3］COGGON J H.Electromagnetic and electrical modeling by the finite element method［J］.Geophysics，1971，36（2）：132－151.

［4］朱伯芳.有限單元法原理與應(yīng)用［M］.北京：水利出版社，1979.

［5］周熙襄，鐘本善，江東玉.點(diǎn)源二維電阻率法有限差分正演計(jì)算［J］.物探化探電子計(jì)算技術(shù)，1983，5（3）：34－38.

［6］羅延鐘，萬樂.二維地形不平條件下外電場(chǎng)的有限差分模擬［J］.物探化探計(jì)算技術(shù)，1984，6（4）：15－26.

［7］羅延鐘，萬樂.用數(shù)值模擬方法構(gòu)組保角變換坐標(biāo)網(wǎng)［J］.物探化探計(jì)算技術(shù)，1986，8（1）：23－31.

［8］陳小斌.有限元直接迭代算法［J］.物探化探計(jì)算技術(shù)，1999，21（2）：165－171.

［9］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994.

［10］陳樂壽.有限元法在大地電磁測(cè)深正演計(jì)算中的應(yīng)用與改進(jìn)［J］.石油物探，1981，20（3）：84－103.

［11］陳樂壽，孫必俊.有限元法在大地電磁測(cè)深中正演計(jì)算的改進(jìn)［J］.石油地球物理勘探，1982，17（3）：69－72.