任意四邊形剖分二維MT有限元正演模擬

畢 武，劉 云

（1．烏魯木齊金維圖文信息科技有限公司，烏魯木齊 830091；2．新疆地礦局物化探大隊，昌吉 831100；3．中國科學院地球化學研究所，貴陽 550002）

任意四邊形剖分二維MT有限元正演模擬

畢 武1，2，劉 云3

（1．烏魯木齊金維圖文信息科技有限公司，烏魯木齊 830091；2．新疆地礦局物化探大隊，昌吉 831100；3．中國科學院地球化學研究所，貴陽 550002）

為了模擬實際的起伏地形，在前人研究的基礎上，單元網格設計為任意四邊形網格，單元內的場值雙線性變化，采用了高斯數值積分法計算單元系數矩陣，并給出二維大地電場的輔助場表達式。通過模型算例表明，模擬結果與解析解的均方誤差在1%以內，地形模型與前人的模擬結果相吻合。分析對比了兩個不同坡度模型對視電阻率和相位的影響。采用任意四邊形網格剖分法，降低了編程難度，可以方便地適應野外地形的起伏變化。

大地電磁；起伏地形；任意四邊形網格；有限單元法

0 引言

大地電磁測深法（MT）的野外觀測，往往在地形條件復雜的山區進行，起伏地形對于電磁場的影響非常大，有必要研究帶地形的 MT數值模擬技術。

Wannamaker［1］、陳小斌［2］等采用自適應地形的三角網格剖分，雙線性插值有限元數值模擬起伏地形影響。劉云［3－4］等采用自適應地形的曲邊四邊形網格剖分雙二次插值方法，以及采用四邊形結合三角形剖分的電導率連續變化的二維有限元模擬技術，對起伏地形的大地電磁法觀測取得較好的效果。

作者在前人研究基礎上，導出了任意四邊形單元剖分、雙線性插值的大地電磁二維數值計算方法。實現了以高斯數值積分計算單元系數矩陣以及輔助場和視電阻率的計算。

1 變分問題

選擇如圖1所示的坐標系，假定地下電性結構是二維的，取走向為z軸，那么三維地電模型可以簡化為二維地電模型來進行處理。經過推導與 MT二維邊值問題相應的變分問題為［11］：

其中 ：σ為介質電導率；μ為磁導率取為常數μ≈μ0＝4π×10－7；ω是電磁波的角頻率；k是傳播系數，k

2 有限單元法

2．1 區域剖分

對于形如圖1所示的地形，采用四邊形網格剖分，可在地面各點加入高程信息形成二維起伏地形。由于網格單元不再是常規的矩形單元，將不滿足矩形剖分的公式，需導出任意四邊形單元線性插值場分量計算方法。公式（1）離散之后為

圖1 區域剖分示意圖Fig．1 Division of region with FEM grids

在圖1中，加入高程信息的計算公式為

其中：Yp是點p的高程數據（水平地形條件下Yp＝0，正地形時Yp＞0，負地形Yp＜0）；yi是水平地形下，y方向上第i個節點的y坐標；Y為水平條件下y方向的最大距離；N為y方向上剖分的網格節點個數。帶入此公式運算之后，Yi即為y方向上第i個節點添加了高程信息之后的y坐標。這樣高程便以一個較小的變化依次加入到y方向上的節點中［13］。對于TE模式，空氣層中高程信息的加入類比于TM模式中對于地下網格的處理。

圖2 任意四邊形單元Fig．2 Arbitrary quadrilateral element

2．2 單元分析

單元如圖2所示。

雙線性插值的形函數為

其中ξi、ηi是點i（i＝1，2，3，4）的坐標。圖2（a）所示的二維等參單元分別表示場值u和坐標x、y為：

將每個單元的矩陣按照單元和節點編號帶入總體矩陣中，采用變帶寬存儲總體矩陣。F（u）＝UTKU，令其變分為0，得到KU＝0，帶入上邊界條件之后，解方程組，就可以得到各個節點的U值。

當網格較多時，對所有單元系數均使用高斯數值積分計算是比較費時的。考慮到剖分區域時大量的單元是矩形的，只有涉及地形的部分是四邊形的，所以在單元分析的程序部分，加入一個條件判斷，如果該單元是矩形，那么直接使用常規系數矩陣，若不是矩形單元，則使用上述方法計算系數矩陣。下面給出矩形線性插值的單元系數矩陣：

3 計算視電阻率和阻抗相位

3．1 輔助場的求取

使用四邊形剖分，根據單元分析中的推導有：

3．2 視電阻率和阻抗相位的計算

在TE模式中，由于磁場的連續性，起伏地形條件下Ez和Hx均可測得，所以起伏地形TE模式的視電阻率和阻抗相位的定義與水平地形條件相同。

在TM模式中，磁場分量Hz可以測得，但是電場在傾斜的地形線上不連續，電場水平分量Ex無法測得，只能測得沿地形線的El。根據矢量分解的定義，El分解為Ex和Ey分量［2］，即

此時，視電阻率和阻抗相位的定義為：

易知當水平地形時，由于Ey為零，同樣可以應用該表達式。

4 模型算例

在以下模型算例中使用的網格剖分相同，研究區域兩側分別取18個稀疏網格，TE模式下，空氣層取14個網格，地下介質剖分為56個網格（y方向），求解區域x方向剖分的網格數取為測點數。

4．1 層狀模型

模型如圖3所示，水平層狀KH型地電模型（4層），模型參數為：第一層ρ1＝20Ω·m，h1＝200 m；第二層ρ2＝200 Ω·m，h2＝1 000 m；第三層ρ3＝20Ω·m，h3＝2 500 m；基底層ρ4＝200Ω·m。地面測線長4 km，測點間距100 m，一共41個測點，頻率范圍為1 000 Hz～0．01 Hz，對數等間隔采樣，一共41個頻點。

使用有限單元法求解，在CPU為T6600（2．2 GHz），內存2 GB的計算機上計算，TE模式每個頻點的計算時間為0．58 s，TM模式每個頻點的計算時間為0．31 s。

圖3 KH型地電模型Fig．3 Schematic of KH－type geoelectric model

取21號測點（測區中間）的數值解與該模型的解析解進行對比，如圖4所示。

圖4 KH型模型解析解數值解對比Fig．4 Comparison between analytical and numerical solutions by FEM

從圖4中可以看到，數值解與解析解是基本重合的。通過計算其他模型與解析解對比發現，視電阻率和阻抗相位值的均方誤差均小于1%。

4．2 梯形山模型

模型具體參數如圖5所示，該模型引自文獻［1］，測線水平寬度為4 km，測點距50 m，梯形山位于測區中心，背景電阻率ρ＝100Ω·m，求解的頻率為2 Hz。

模擬結果如圖6所示。

對比文獻［1］中的結果，可以看出響應曲線的形態和值的大小是一致的。

4．3 斜坡地形1

圖5 梯形山模型Fig．5 Trapezoidal hill

如圖7所示，測線水平寬度為6 km，測點距50 m，共121個測點，斜坡高度為1．2 km，斜坡跨度為3．0 km，斜坡坡度為23．6°，背景電阻率ρ＝100Ω·m，求解的頻率為1 Hz。模擬結果的視電阻率以及相位曲線，如圖8所示。

圖6 梯形山模型2 Hz TE，TM模式視電阻率以及相位響應曲線Fig．6 TE／TM response curve at 2Hz of trapezoidal hill

對比文獻［1，12］中的模擬結果，視電阻率以及相位的響應曲線基本一致。

從梯形山模型以及斜坡模型的數值模擬結果可以看出，TM模式的視電阻率以及相位受地形的影響較為嚴重，而TE模式受到影響相對較小。

4．4 斜坡地形2

下面通過對兩個不同坡度模型的模擬結果，對比不同頻率以及坡度對大地電磁場測量的影響。

如圖9所示，斜坡高度為2 km，斜坡跨度為3．0 km，斜坡坡度為41．8°，背景電阻率ρ＝100Ω·m，測量參數與4．3中一致，求解頻率為0．1 Hz、1 Hz以及10 Hz。

兩個坡度的視電阻率模擬結果如圖10所示。

圖8 斜坡模型10 Hz的TE，TM響應曲線Fig．8 TE／TM response curve at 10Hz of the ramp terrain model

圖9 斜坡地形2Fig．9 Ramp terrain model 2

在圖10（a）、（b）中，為了突出對視電阻率的影響，縱軸為線性變化，從模擬結果可以看出，低頻時，TE模式受到的影響較小，越往高頻，受到地形的影響越大，同時對比圖（a）、（b）兩圖可以看出，斜坡2視電阻率受到的影響相對于斜坡1更大，說明坡度越大，對視電阻率的影響越大。圖10（c）、（d）中縱軸為對數值，可以看出TM模式下，視電阻率在低頻和高頻均受到很大的影響，并且坡度越大，其視電阻率畸變越大。在遠離坡腳1 km處，高頻視電阻率趨于背景電阻率，而低頻視電阻率依然受到斜坡的影響。

從圖11中發現，斜坡對于TE模式相位的影響相對較小，對TM模式的影響較大，并且坡度越大，阻抗相位的率畸變越大。在遠離坡腳1 km處，高頻的阻抗相位趨于45°，而低頻的阻抗相位依然受到斜坡的影響。這些影響特點，與斜坡對視電阻率的影響是一致的。

5 結論

使用任意四邊形剖分方式可以方便的適應地形變化，同時因為總體網格的剖分是采用矩形剖分，所以構建網格的部分是比較簡單的。通過對層狀模型以及地形模型的模擬，驗證了方法的正確性。對斜坡地形，TE模式視電阻率和相位受地形影響較小，TM模式受到的影響較大。坡腳越大，視電阻率與阻抗相位畸變越大。遠離斜坡測點的視電阻率和阻抗相位依然受到斜坡的影響，并且低頻受到的影響更為嚴重。

圖10 視電阻率曲線圖Fig．10 Curve of apparent resistivity

圖11 阻抗相位曲線圖Fig．11 Curve of impedance phase

［1］ WANNAMAKER P E，STODT J A，RIJOFI L．Two－dimensional topographic responses in magnetotellurics modeled using finite elements［J］．Geophysics，1986，51（11）：2131－2144．

［2］ 陳小斌．MT二維正演計算中地形影響的研究［J］．石油地球物理勘探，2000，39（3）：112－120．

［3］ 劉云，王緒本．大地電磁二維自適應地形有限元正演模擬［J］．地震地質，2010，32（3）：382－391．

［4］ 劉云，王緒本．電性參數分塊連續變化二維MT有限元數值模擬［J］．地球物理學報，2012，55（6）：2079－2086．

［5］ 樸化榮．電磁測深法原理［M］．北京：地質出版社，1990．

［6］ 王緒本，李永年，高永才．大地電磁測深二維地形影響及其校正方法研究［J］．物探化探計算技術，1999，21（4）：327－332．

［7］ 史明娟，徐世浙，劉斌．大地電磁二次函數插值的有限元法正演模擬［J］．地球物理學報，1997，40（3）：421－430．

［8］ 徐世浙，王慶乙，王軍．用邊界單元法模擬二維地形對大地電磁場的影響［J］．地球物理學報，1992，35（3）：380－388．

［9］ 徐世浙，于濤，李予國，等．電導率分塊連續變化的二維MT有限元模擬（Ⅰ）［J］．高校地質學報，1995，1（2）：65－73．

［10］李予國，徐世浙，劉斌，等．電導率分塊連續變化的二維MT有限元模擬（Ⅱ）［J］．高校地質學報，1996，2（4）：448－452．

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

［12］尤淼．電性參數分塊連續的大地電磁二維有限元數值模擬［D］．成都：成都理工大學，2011．

［13］呂玉增，阮百堯．復雜地形條件下四面體剖分電阻率三維有限元數值模擬［J］．地球物理學進展，2006，21（4）：1302－1308．

［14］彭國倫．Fortran 95程序設計（第一版）［M］．北京：中國電力出版社，2002．

［15］徐士良．計算機常用算法（第二版）［M］．北京：清華大學出版社，1995．

Two－dimensional magnetotelluric forward modeling by arbitrary quadrilateral finite element method

BI Wu1，2，LIU Yun3
（1．Urumqi Jinwei Map－character Information And Science＆Technology Co．Ltd，Urumqi 830091，China；2．Geophysical and Geochemical Party，Xinjiang Bureau of Geology and Mineral Resource Exploration and Development，Changji 831100，China；3．Institute of Geochemistry Chinese Academy of Sciences，Guiyang 550002，China）

In order to easily model topography in field work，on the basis of predecessors＇achievements，the arbitrary quadrilateral element grid was used in the finite element method（FEM），and the electromagnetic field of models were designed to bilinear variation within each quadrilateral element in our numeric modeling．We deduced the specific formula to calculate the unit coefficient matrix by Gaussian numerical integral and the expression of the auxiliary field when using arbitrary quadratic element．Model calculation shows that the mean square error between the simulation results with the analytical solution is less than 1 percent，and the results of terrain model is consistent with the one of other scholar＇s．Through modeling of two ramp terrain modes with different pitch，comparing and analyzing the effect of the apparent resistivity and the phase．Due to the use of arbitrary quadrilateral element grid，the difficulty of programing was reduced，forward modeling can be more easily adapt to the topography variation．

MT；topography；arbitrary quadratic element；FEM

P 631．3＋25

A

10．3969／j．issn．1001-1749．2014．05．02

1001-1749（2014）05-0521-08

2014-01-17 改回日期：2014-07-17

畢武（1970－），男，高級工程師，主要從事物化探數據處理解釋等方面的研究工作，E-mail：382280581＠qq．com。