吳 娟，席振銖，王 鶴

（中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

0 前言

有限元作為一種高效的數(shù)值模擬方法，由于不需要考慮物體的內(nèi)部邊界條件，而被廣泛應(yīng)用于地球物理的正演計算。作為大地電磁的數(shù)值模擬，最早是由Coggon［1］于1971年引入的。此后Wannamaker［2］、Mitsuhata［4］、徐世浙［5］、李予國［6］、阮百堯［7］、王緒本［8］等，相繼對有限元的二維三維模擬，研究了地形、插值函數(shù)以及介質(zhì)的電導(dǎo)率連續(xù)性等問題。隨著計算機的不斷發(fā)展，人們對正演模擬的精度要求越來越高，計算結(jié)果的誤差直接影響著后續(xù)反演的效果。因此，許多學(xué)者又從多個方面對有限元的算法加以改進。陳小斌［9］提出了一種新的有限元直接迭代算法，馬為［10］等又在此基礎(chǔ)上選用二次插值基函數(shù)計算輔助場的函數(shù)值；張秀敏［11］改進了電磁場計算中大型稀疏方程組的迭代解法；柳建新［12］等將帶預(yù)處理的雙共軛梯度算法運用到帶地形電磁的二維、三維正演中，使求解線性方程的精度進一步提高。

在眾多的研究中，都只是在有限元的插值函數(shù)或者解方程的精度上做了改進，而對于有限元前處理中的細節(jié)問題，討論得并不多。眾所周知，有限元研究區(qū)域的尺寸問題，直接影響到正演計算的精度和速度，由網(wǎng)格剖分的不合理性所導(dǎo)致的誤差，遠大于解方程所引起的誤差。如果網(wǎng)格剖分得不夠合理，那么再高精度的求解系統(tǒng)，也得不到正確的解。由于矩形網(wǎng)格線性插值在求解地球物理問題中具有典型的代表性，因此，作者在具體研究了矩形網(wǎng)格雙線性插值中網(wǎng)格尺寸以及邊界條件對計算結(jié)果的影響后，總結(jié)出相應(yīng)規(guī)律，以供后續(xù)計算。

1 變分問題

構(gòu)造如圖1所示（見下頁）坐標系。

有限元解大地電磁二維問題的邊值問題為：

圖1 大地電磁有限元模擬研究區(qū)域及坐標Fig.1 Study area and coordinates for MT modeling use finite element method

其中 u代表解線性方程組后得到的節(jié)點場值，對于 TE模式：u＝Ez；τ＝1／（iωμ）；λ＝σ。對于 TM模式：u＝ Hz；τ＝iωμ；λ＝1／σ；ω為角頻率；σ和μ分別代表介質(zhì)電導(dǎo)率和磁導(dǎo)率。

經(jīng)過變分之后，都得到了一個關(guān)于有限元系數(shù)矩陣的線性方程組，即：Ku＝0，其中K即有限元的剛度矩陣。在計算剛度矩陣時，先要利用插值函數(shù)來計算單元系數(shù)矩陣Ke，而等參單元的插值基函數(shù)，其系數(shù)不僅與單元的節(jié)點坐標相關(guān)，而且與單元的寬度和高度相關(guān)。

在求解線性方程組時，要帶入邊界條件，其中內(nèi)部邊界條件在泛函取極值的過程中自動滿足，屬于自然邊界條件，將不予考慮。剩下就是外邊界條件，即：

一般的處理是，對于上邊界AB，由于電磁波在空氣中幾乎不衰減，要使異常場在邊界上的偏導(dǎo)數(shù)為零，就必須把其放置在無窮遠，這樣就可以令u＝1。對于左右邊界，也采取同樣的方式。而對于下邊界，同樣也要使異常場在邊界上產(chǎn)生的值為零，由于下邊界以下為均勻半空間，因此，它所滿足的邊界條件為第三類邊界。

求解方程組，得到各個節(jié)點的場值之后，還必須求得相應(yīng)的視電阻率值：

在計算偏導(dǎo)數(shù)時，可以采取二點差分公式，也可以采用三點差分公式或者更高次的差分公式。為了提高視電阻率精度，作者在本文中采用了四點差分公式。取近地表的連續(xù)的四個節(jié)點，為了計算方便，一律把位于近地表一層連續(xù)的四個節(jié)點采用等距剖分，ui、ui＋1、ui＋2、ui＋3分別代表從地表連續(xù)向下分布的四個節(jié)點，則：

其中 l是四個節(jié)點之間的距離，帶入公式（3）、公式（4）中，就能得到TE和TM模式下的視電阻率了。

從公式（5）中，我們可以得出以下結(jié)論：

（1）采用四點差分格式計算微分時，必須是等步長，也即地表以下連續(xù)三個網(wǎng)格必須為相同的間距。

（2）偏導(dǎo)數(shù)的值與網(wǎng)格步長息息相關(guān)，即網(wǎng)格剖分的疏密程度將直接影響著計算結(jié)果的精度。

根據(jù)有限元滿足的邊界，又必須將網(wǎng)格邊界設(shè)置在無限遠處，這樣必將增加計算量。那么，究竟如何選擇合適網(wǎng)格參數(shù)來既保證計算精度又不太消耗計算機內(nèi)存呢？網(wǎng)格的邊界取值的遠近是否會影響計算結(jié)果呢？下面將具體從邊界影響和網(wǎng)格橫、縱網(wǎng)格步長三個方面來討論。

圖2 一維和二維介質(zhì)模型圖Fig.2 One－dimensional and two－dimensional models

2 模型計算

2.1 左右邊界的影響

為了研究左右邊界對計算結(jié)果的影響，作者分別以均勻半空間、水平層狀介質(zhì)和二維的低阻嵌入體模型為例，模型如圖2所示。

其中，均勻半空間電阻率為100Ω·m。采用30（縱）×32（橫）的網(wǎng)格，橫縱單元均為等距剖分，網(wǎng)格步長均為10m。在此基礎(chǔ)上，令橫向第一個和最后一個網(wǎng)格的左右邊界，分別為±160m、±500m、±1km、±10km、±100km，頻率范圍從10－3Hz～103Hz，對數(shù)等間距分布21個頻點，所得視電阻率測深曲線如下頁圖3所示。

在圖3（a）中，視電阻率的幅值在100Ω·m附近，最大誤差不超過0.5%，由此驗證了程序的正確性。結(jié)合圖3（b），在左右邊界不斷延伸的過程中，其數(shù)值模擬的結(jié)果也與解析解十分吻合，由此可以看出，左右邊界對均勻和一維介質(zhì)幾乎沒有影響，因此，可以不用將左右邊界設(shè)置在無窮遠處。

圖3（c）中的二維介質(zhì)，由于沒有解析解，特給出一個參考解。根據(jù)有限元基本理論，當研究區(qū)域無限大，網(wǎng)格足夠細時，任何一個模型都能得到一個近似的真解。以此為依據(jù)，對于尺寸為200m×100m，埋深為100m的低阻嵌入體，將左右邊界設(shè)置在±100km處，下邊界設(shè)置在－10km處，網(wǎng)格60（縱）×62（橫），以橫向±150m和縱向0m～－250m范圍為目標區(qū)，采用5m的步長等距剖分。由于頻率最高為103Hz，趨膚深度為160m，因此5m的步長可認為足夠細。在目標區(qū)域以外，步長逐漸放大，直至網(wǎng)格的邊界，由此解出的數(shù)值解定為二維介質(zhì)的參考解。

在圖3（c）中，隨著邊界的不斷擴大，其曲線中段的幅值總體向上抬升，但曲線的起點和終點卻沒有移動，曲線的基本形態(tài)沒有變化。我們可以這樣解釋，對于同樣的200m×100m的低阻嵌入體，在左右邊界不斷擴大的過程中，其圍巖的范圍也就擴大了。相對而言，模型尺寸與周圍的圍巖比例就縮小了，而由于視電阻率測深曲線，是地下介質(zhì)體的綜合反映，因此，曲線的幅值也就有所偏向圍巖。但當邊界為500m以上時，曲線形態(tài)基本和參考解重合，由此可見，當邊界為目標體尺寸五倍以上時，就可以顯著削弱邊界條件的影響了。

縱觀圖3地表視電阻率的計算，其左右邊界并不需要取到無窮遠處。雖然對于二維體，不同邊界的曲線并沒與參考解完全重合，但它所反映的淺層和深部的視電阻率信息與實際是沒有區(qū)別的，而且對于曲線形態(tài)的改變并沒有質(zhì)的影響。這是由于在泛函取極值的過程中，要使得異常場在左右邊界上的法向偏導(dǎo)數(shù)的值為零的真正目的，是為了消去變分過程中左右邊界的積分項，使得計算更為便捷。而實現(xiàn)這一目的最簡單的方式，就是讓邊界取到無窮遠。然而，對于圖3所示的介質(zhì)和模型，都是關(guān)于左右對稱的，由于其特殊的對稱性，使得其在邊界上的影響總是大小相等方向相反的，這樣一來，使得函數(shù)在邊界上的積分和為零，自然滿足條件。因此，可以得到這樣的結(jié)論，即對于左右對稱良好的介質(zhì)，其左右邊界并不需取到無窮遠。

圖3 左右邊界對MT計算精度的影響Fig.3 Effects of left and right boundary on accuracy of MT

2.2 邊界的影響

為了研究下邊界對MT計算精度的影響，作者仍然以模型1和模型2為例。網(wǎng)格仍然采用30（縱）×32（橫）的剖分方式，橫縱網(wǎng)格分別為均勻剖分，步長10m，此時最后一個網(wǎng)格所示下邊界為地表以下260m。在此基礎(chǔ)上，分別讓最后一層的網(wǎng)格邊界改為500m、1km、5km，所得視電阻率曲線如圖4所示。

在圖4中，層狀介質(zhì)四條不同底邊界的曲線與解析解吻合良好。低阻嵌入體同樣也引入和2.2節(jié)中相同的參考解，不同底邊界的響應(yīng)和參考解也十分吻合。這也就表明，下邊界的遠近，并不會影響大地電磁的計算精度。

2.3 橫向網(wǎng)格尺寸的影響

為了研究橫向網(wǎng)格尺寸對計算結(jié)果的影響，作者分別以均勻半空間100Ω·m、模型1和模型2為例，采用30（縱）×32（橫）的網(wǎng)格，橫縱單元均為等距剖分，縱向網(wǎng)格步長為10m，橫向網(wǎng)格分別取單元寬度為100m、50m、20m、10m、2m，相對于單元高度為10m，其橫縱比分別為：10／1、5／1、2／1、1、1／5?？紤]到左右邊界的因素，將水平方向第一個和最后一個網(wǎng)格坐標設(shè)為±1km。其相應(yīng)的視電阻率曲線如圖5所示。

圖5（a）和圖5（b）中的數(shù)值模擬解和解析解都對應(yīng)的十分良好，由此可見，橫向網(wǎng)格寬度以及橫縱比對均勻介質(zhì)和一維介質(zhì)都沒有影響。對于圖5（c）所對應(yīng)的模型2，為低阻嵌入體，當左右邊界相同的情況下，網(wǎng)格橫縱比小于10／1時，都與參考解對應(yīng)良好。

由此可見，大地電磁的有限元數(shù)值模擬，其橫向網(wǎng)格的剖分可相對隨意些。

2.4 縱向網(wǎng)格尺寸的影響

從上兩節(jié)得出的結(jié)論可知，橫向網(wǎng)格對視電阻率結(jié)果的影響并不大，對稱介質(zhì)邊界的遠近對視電阻率的計算影響也很小。那么縱向網(wǎng)格對精度的影響又如何呢？在歐東新［13］等人的研究中提到，對于不同頻率的大地電磁測深，不能采用同一條網(wǎng)格體系，在高頻時要采用小網(wǎng)格，低頻時采用大網(wǎng)格。這主要針對的也是縱向網(wǎng)格的尺寸，卻沒有給出具體的類比參數(shù)。而且對于頻帶相對廣泛的MT正演模擬而言，采用二套網(wǎng)格也較不現(xiàn)實。

眾所周知，電磁波在地下傳播時，能量是衰減的，衡量其探測能力的一個重要參數(shù)，就是趨膚深度：

趨膚深度與電阻率成正比，和頻率成反比。顯然，當頻率越高或電阻率越低時，其探測深度越淺。

在正演模擬的過程中，作者根據(jù)公式（5），利用差分來代替微分，以求得地面節(jié)點值的偏導(dǎo)數(shù)。當頻率過高時，如果網(wǎng)格寬度大于一個趨膚深度，顯然其精度是不夠的。由此可見，一套網(wǎng)格最終決定它的縱向網(wǎng)格剖分的，是其最小的趨膚深度。為分析網(wǎng)格剖分與趨膚深度的相互關(guān)系，作者仍然以均勻半空間100Ω·m和模型1為例，分別分析影響規(guī)律。

首先，在均勻半空間，采用30（縱）×32（橫）的網(wǎng)格。由于計算視電阻率時，要求地表以下四個節(jié)點是等步長的，為簡單起見，橫縱單元均為等距剖分，橫向網(wǎng)格步長為10m，左右邊界為±1km，頻帶范圍10Hz～105Hz，重點分析高頻淺部的情況。由2.2小節(jié)計算得出，下邊界對計算精度影響不大，由此可見，由縱向網(wǎng)格變化引起的下邊界的遠近影響可忽略。根據(jù)趨膚深度的計算公式，電阻率100Ω·m的介質(zhì)，當頻率達到100kHz時，其趨膚深度約為16m。若分別令剖分步長為2m、4m、8m、16m，相應(yīng)的縱向網(wǎng)格與趨膚深度之比分別為1／8、1／4、1／2、1／1，可以做出相應(yīng)視電阻率曲線如圖6所示。

在圖6中，網(wǎng)格相對于趨膚深度越大，其誤差越大，并且其誤差隨著頻率的增高而加大。這可以解釋為同樣的網(wǎng)格，在低頻時遠遠小于趨膚深度，而在高頻，就不小于甚至?xí)笥谮吥w深度。在圖6中的反映就是10kHz以下的頻率，正演結(jié)果基本接近真值。對于超高頻的100kHz，其趨膚深度為16m，此時用大于1／4個趨膚深度的波長，其視電阻率的最大誤差就已超過3%，可見此時已不滿足精度要求了。

同樣的，從圖6中還能發(fā)現(xiàn)，尺寸越小精度越高。但一味的追求更小的尺寸，同時也增加了模型的計算量。用1／4個和1／8個趨膚深度步長達到的最大誤差，都不會超過1%。由此可見，滿足精度要求1%的條件下，一個趨膚深度范圍內(nèi)，用四個網(wǎng)格就可以達到很好的計算要求了，這樣既保證了精度，又節(jié)省了計算量。當然，此結(jié)果只針對于線性插值，當插值函數(shù)為二次或者更高次時，則只需要更少的網(wǎng)格數(shù)。

對模型1所示的層狀一維介質(zhì)，我們采用頻率范圍為10－3Hz～103Hz，對應(yīng)于地表一層為100Ω·m的介質(zhì)，其趨膚深度最小約為160m。為此我們采用和均與半空間相似的做法，保持橫向剖分10m的步長，縱向網(wǎng)格等距剖分，分別為10m、20m、50m和100m，那么其相對于趨膚深度的比值為1／16、2／16、5／16、10／16，并與解析解對比，如圖7所示。

在圖7中，可以發(fā)現(xiàn)仍然和圖6一樣，當尺寸小時精度高，在小于1／8個趨膚深度范圍內(nèi)計算的結(jié)果和解析解對應(yīng)良好，而當網(wǎng)格尺寸大于5／16時，其計算就開始偏離真值了，而且頻率越高，其畸變越明顯。分析其原因有二：①因為隨著頻率的升高，其趨膚深度不斷減小，因此對網(wǎng)格的尺寸要求更小；②模型2所示的層狀介質(zhì)，其層厚度為100m，當縱向網(wǎng)格尺寸大于50m時，三個網(wǎng)格的深度就已超過了100m，對于公式（3）、公式（4）所計算的地表節(jié)點場值的偏導(dǎo)數(shù)，就跨越了物性分界面，因此進一步加大了計算誤差。所以，在縱向網(wǎng)格剖分時，還應(yīng)當注意不但與最小趨膚深度有關(guān)，而且地表以下連續(xù)三個網(wǎng)格不能跨越物性分界面。

3 結(jié)論

（1）針對矩形網(wǎng)格雙線性插值的有限單元模擬MT問題，下邊界對正演精度影響可忽略；而左右對稱的介質(zhì)，左右邊界對正演精度影響不大，故可不用將邊界放置在無窮遠處。

（2）橫向網(wǎng)格寬度以及單元橫縱比，對計算結(jié)果的影響不大，剖分可相對隨意些。

（3）縱向網(wǎng)格的剖分與趨膚深度有密切的關(guān)系，對于線性插值，網(wǎng)格步長應(yīng)小于最小趨膚深度的1／4。

［1］ COGGON G H.Electromagnetic and electrical modeling by the finite element［J］.Geophysics，1971（36）：132.

［2］ WANNAMAKER P E，STODT J A，ROJO L.PW2D－Finite element program for solution of magnetotelluric responses of two－dimensional earth resistivity structure［D］.Univ.of Utah Res.Inst.Rep，ESL，1985.

［3］ WANNAMAKER P E，STODT J A，RIJO L.Twodimensional topographic responses in magnetotellurics modeled using finite elements［J］.Geophysics，1986（51）：2131.

［4］ MITSUHATA Y，UCHIDA T.3Dmagnetotelluric modeling using the T－Ωfinite－element method［J］.Geophysics，2004（69）：108.

［5］ 石明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學(xué)報，1997，40（3）：421.

［6］ 李予國，徐世浙，劉斌，等.電導(dǎo)率分塊連續(xù)變化的二維 MT有限元模擬（Ⅱ）［J］.高校地質(zhì)學(xué)報，1996，2（4）：448.

［7］ 阮百堯，徐世浙.電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬［J］.中國地質(zhì)大學(xué)學(xué)報，1998，23（3）：303.

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

［9］ 陳小斌，張翔，胡文寶.有限元直接迭代算法在 MT二維正演計算中的應(yīng)用［J］.石油地球物理勘探，2000，35（4）：487.

［10］馬為，陳小斌，趙國澤.大地電磁測深二維正演中輔助場的新算法 ［J］.地震地質(zhì)，2008，30（2）：525.

［11］張秀敏，苑津莎，徐永生.電磁場分析中大型稀疏方程組迭代解法的改進 ［J］.華北電力大學(xué)學(xué)報，2003，30（3）：34.

［12］柳建新，蔣鵬飛，童孝忠，等.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用 ［J］.中南大學(xué)學(xué)報，2009，40（2）：484.

［13］歐東新.計算機精度和網(wǎng)格大小對大地電磁有限單元法正演的影響 ［J］.桂林工學(xué)院學(xué)報，2007，27（3）：329.