附加列插值法在起伏地形條件下二維大地電磁正演中的研究與應用

吳 頔，辛慧翠，卜傳新

(1.湖南省有色地質勘查局二四七隊，湖南 長沙 410129；2.承德石油高等專科學校 計算機與信息工程系，河北 承德 067000)

起伏地形對大地電磁場影響很大，研究表明二維帶地形條件下TM 和TE模式電阻率曲線具有不同的畸變特征，這種畸變既可能是電阻率數值發生畸變也可能是發生靜態位移, 或者是二者的綜合影響，畸變一般與相對波長的地形幾何尺度有關，大地電磁測深地形影響的復雜性使得對它的有效校正較為困難, 本文利用雙二次插值有限元法分別推導出了TM和TE模式在復雜地形條件下的二維正演公式，采用整體法和附加列法兩種方法進行計算，并對計算結果進行了分析和評價。

1 變分問題

1.1 基本理論

根據麥克斯韋方程，假定二維地下電性結構，各坐標軸關系如圖1所示，在均勻半空間里當平面電磁波以任何角度入射地面時，地下介質中的電磁波總以平面波形式近似垂直地向下傳播，將(1)、(2)式按分量展開，并考慮到?/?z=0，得兩個獨立的方程組，以垂直分量為準，分別命名為TE和TH極化模式[1]。

TE極化模式：

(1)

TM極化模式：

(2)

式(1)、(2)可以統一寫成下式：

(3)

其中：μ是介質的磁導率，σ是電導率，ε是介電常數，ω是角頻率， e-iωt是時間因子，是二維哈密頓算符。

1.2 邊值問題和變分問題

求解微分方程必須給出邊界條件，建立以AB、CD為邊界的模型(見圖1和圖2)，邊值問題可以歸納為:

(4)

與邊值問題等價的變分問題可以歸納為：

(5)

2 有限元法

采用不均勻矩形單元對區域進行剖分，每個矩形網格單元上取4個頂點和4條邊的中點, 每個矩形網格的編號如下圖所示[2-3]，圖3為母單元, 圖4為子單元, 兩者的坐標關系為:

(6)

其中x0、y0為子單元的中點,a、b為子單元的兩個邊長。

將(5)式分解為單元積分之和，當單元中的τ和λ是常數時，各項單元積分可分別寫為：

(7)

(8)

對(8)式求變分，由于δu具有任意性，因此可得：

Ku=0

(9)

加載式(4)中AB的邊界條件，形成右端項進而可得到大型稀疏方程組，利用迭代法求解線性方程組。視電阻率計算公式如下：

TE型極化：

(10)

TM型極化：

(11)

在求得地表各節點的u后還要計算其對y方向的導數，取近地表的4個等距節點的u及節點1與節點4間的距離來計算導數，將計算結果代入公式(10)、(11)即可計算電阻率[4-6]。

3 程序實現及驗證

3.1 程序實現

一般利用有限元模擬計算大地電磁2D正演分為以下幾個步驟：1)網格剖分，單元格內電導率連續變化；2)單元內插值并進行單元分析，一般以整個模型為研究對象，對所有網格單元進行編號；3)總體合成；4)代入邊界條件，解線性方程組；5)解出各節點u值，求導后計算視電阻率，計算過程中對模型整體進行插值因此這種方法也被稱為整體法。在此基礎上作者提出了一種改進的插值方法，以第N列(N=1…L，L為網格列數)網格為對象，增加一列與其完全一致的網格，以這兩列為對象進行編號和計算，然后取這兩列網格地表三個節點中的第二個節點及其下面的三個節點的U值來計算視電阻率值，所得結果即為整個模型中第N列網格表層中間點的視電阻率值，以此類推，從第一列循環計算到最后一列，從本質上講該方法是對若干個層狀模型進行組合，可有效改善一般算法中因地形起伏引起的分界面上的電阻率畸變，該方法稱之為附加列法，計算思路見圖5。

3.2 程序驗證

為驗證程序可靠性，設計三層一維模型，第一層厚度為300 m，電阻率值500 Ω·m，第二層厚度為500 m，電阻率值50 Ω·m，第三層厚度為900 m，電阻率值1 500 Ω·m，經正演計算出不同頻率的電阻率計算(圖6)。

一維模型電阻率只在深度方向發生變化，理論上無論TE模擬或者TM模式，電阻率計算結果與解析解應近似，由上圖可以看出整體法和列分解法在TE和TM模式下都能較好的模擬電阻率變化，這也驗證了方法的正確性。

4 二維模型計算

設計二維模型長度4 km，第一層中間厚度250 m，最大厚度為450 m，最小厚度為100 m，分別模擬山脊和溝谷；第二層厚度為300 m；第三層厚度為450 m(圖7)。

圖8～圖11是應用整體法和列分解法計算的視電阻率擬斷面圖。其中圖8是TE模式下整體法電阻率計算結果，該結果能較清晰地反映地形引起的電阻率變化且與地形構造大致吻合，但在地形變化分界面上存在由地下電性結構變化引起的干擾，表層高阻電性層未能較好體現。圖9是TM模式整體法電阻率計算結果，正演結果在模擬山脊和溝谷處都出現了貫穿了整個頻率帶的低阻假異常，且兩個低阻假異常間還出現了高阻假異常。圖10和圖11是TE和TM模式下附加列法電阻率計算結果，由圖可知兩種極化模式下的正演效果都得到了改善，其中TE模式下表層高阻得到了連續且清晰地反映，也很好的表現出了其隨地形變化的特征，在TM模式下正演結果取得了很大的提高，消除了淺部(高頻部分)高阻假異常以及模擬山脊和溝谷處所對應的低阻假異常，清楚地反映了模型的地形特征。

5 結語

本文利用有限元數值模擬技術在整體法基礎上提出了更為精細的附加列插值法，對起伏地形進行了正演模擬，通過不同方法在TM和TE模式下計算結果的對比分析得出結論：1)整體法模擬中，TE模式正演結果比TM模式正演結果更加接近理論模型，TE模式正演結果基本反映出了地形模型的地電結構特征，模型表層存在些許干擾形成的高阻假異常，但影響較輕微。2)整體法TM模式下電阻率計算結果發生嚴重的畸變和位移，在地形發生變化的分界面上出現了嚴重干擾，假異常貫穿了整個頻帶且幾乎掩蓋了地形。3)附加列法模擬中，TE和TM模式下的正演結果基本相同，能夠很好地反映表層地電特征，基本消除了由地形引起的干擾及假異常。4)通過對模型正演計算結果的分析和研究, 總結了起伏地形電阻率畸變規律，為正確識別和選取極化模式及大地電磁數據校正與處理提供了參考。