吳 頔，辛慧翠，卜傳新

(1.湖南省有色地質(zhì)勘查局二四七隊(duì)，湖南 長(zhǎng)沙 410129；2.承德石油高等?？茖W(xué)校 計(jì)算機(jī)與信息工程系，河北 承德 067000)

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

1 變分問(wèn)題

1.1 基本理論

根據(jù)麥克斯韋方程，假定二維地下電性結(jié)構(gòu)，各坐標(biāo)軸關(guān)系如圖1所示，在均勻半空間里當(dāng)平面電磁波以任何角度入射地面時(shí)，地下介質(zhì)中的電磁波總以平面波形式近似垂直地向下傳播，將(1)、(2)式按分量展開(kāi)，并考慮到?/?z=0，得兩個(gè)獨(dú)立的方程組，以垂直分量為準(zhǔn)，分別命名為T(mén)E和TH極化模式[1]。

TE極化模式：

(1)

TM極化模式：

(2)

式(1)、(2)可以統(tǒng)一寫(xiě)成下式：

(3)

其中：μ是介質(zhì)的磁導(dǎo)率，σ是電導(dǎo)率，ε是介電常數(shù)，ω是角頻率， e-iωt是時(shí)間因子，是二維哈密頓算符。

1.2 邊值問(wèn)題和變分問(wèn)題

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

(4)

與邊值問(wèn)題等價(jià)的變分問(wèn)題可以歸納為：

(5)

2 有限元法

采用不均勻矩形單元對(duì)區(qū)域進(jìn)行剖分，每個(gè)矩形網(wǎng)格單元上取4個(gè)頂點(diǎn)和4條邊的中點(diǎn), 每個(gè)矩形網(wǎng)格的編號(hào)如下圖所示[2-3]，圖3為母單元, 圖4為子單元, 兩者的坐標(biāo)關(guān)系為:

(6)

其中x0、y0為子單元的中點(diǎn),a、b為子單元的兩個(gè)邊長(zhǎng)。

將(5)式分解為單元積分之和，當(dāng)單元中的τ和λ是常數(shù)時(shí)，各項(xiàng)單元積分可分別寫(xiě)為：

(7)

(8)

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

Ku=0

(9)

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

TE型極化：

(10)

TM型極化：

(11)

在求得地表各節(jié)點(diǎn)的u后還要計(jì)算其對(duì)y方向的導(dǎo)數(shù)，取近地表的4個(gè)等距節(jié)點(diǎn)的u及節(jié)點(diǎn)1與節(jié)點(diǎn)4間的距離來(lái)計(jì)算導(dǎo)數(shù)，將計(jì)算結(jié)果代入公式(10)、(11)即可計(jì)算電阻率[4-6]。

3 程序?qū)崿F(xiàn)及驗(yàn)證

3.1 程序?qū)崿F(xiàn)

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

3.2 程序驗(yàn)證

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

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

4 二維模型計(jì)算

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

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

5 結(jié)語(yǔ)

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