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

王 濤，柳建新，童孝忠＊，曹創華，譚神湘

（1.中國海洋大學 海洋地球科學學院，青島 266100；2.中南大學 地球科學與信息物理學院，湖南 長沙 410083）

0 前言

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

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

1 數值模擬

1.1 大地電磁測深正演的基本理論

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

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

式中 ▽·E＝0是因為導電介質內部體電荷密度實際上為零，公式中時間因子都隱含在電場E和磁場H 中，方程組（3）～方程組（6）是大地電磁測深理論研究的出發點。

1.2 一維有限單元正演算法

電場滿足的微分方程為公式（7）。

根據邊值問題所對應的變分問題為

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

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

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

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

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

求解上述線性方程組式（12）即可得到節點處的電場值，從而也可以進一步計算模型響應的視電阻率和相位。

1.3 一維有限差分正演算法

根據電場所滿足的微分方程，將一維地電模型離散化，如圖1所示。

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

然后對電場進行離散處理，采用有限差分近似計算一階、二階偏導數。寫成矩陣形式如式（13）。

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

求解方程組式（14）可得節點處電場值，從而可以進一步計算模型響應的視電阻率和相位。

2 模型算例

2.1 均勻半空間模型

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

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

2.2 層狀介質模型

2.2.1 二層層狀介質

設計一個二層層狀介質模型，ρ1＝10Ω·m，ρ2＝100Ω·m，h1＝1 000m。

（1）解析法正演結果。正演計算的曲線如圖4所示。

（2）有限單元法正演結果。有限單元法的網格剖分如表1所示，有限單元法計算結果與解析法計算結果比較如圖5所示。

公式（13）如下：

公式（14）如下：

（3）有限差分法正演結果。有限差分法的網格剖分如表2所示，有限差分法計算結果與解析法計算結果比較如圖6所示。

2.2.2 三層層狀介質

（1）解析法正演結果。正演計算的曲線如圖7所示。

（2）有限單元法正演結果。有限單元法的網格剖分如表3所示，有限單元法計算結果與解析法計算結果比較如圖8所示。

（3）有限差分法正演結果。有限差分法的網格剖分如表4所示，有限差分法計算結果與解析法計算結果比較如圖9所示。

圖4 二層層狀介質大地電磁響應曲線Fig.4 Two－tier media magnetotelluric response curve

表1 有限單元法網格剖分參數Tab.1 Finite element method meshing parameters

表2 有限差分法網格剖分參數Tab.2 Finite difference method meshing parameters

圖7 三層層狀介質大地電磁響應曲線Fig.7 Three－tier media magnetotelluric response curve

表3 有限單元法網格剖分參數Tab.3 Finite element method meshing parameters

表4 有限差分法網格剖分參數Tab.4 Finite difference method meshing parameters

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

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

3 結論

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

［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］朱伯芳.有限單元法原理與應用［M］.北京：水利出版社，1979.

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

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

［7］羅延鐘，萬樂.用數值模擬方法構組保角變換坐標網［J］.物探化探計算技術，1986，8（1）：23－31.

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

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

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

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