基于四叉樹網(wǎng)格的MT二維正演

王培杰 胡 華 徐 菲 郭 或 陳連木

(①中國地震局地質(zhì)研究所，北京 100029； ②長江大學(xué)地球科學(xué)學(xué)院,湖北武漢 430100； ③油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué)),湖北武漢 430100； ④長江大學(xué)地球物理與石油資源學(xué)院,湖北武漢 430100)

0 引言

在大地電磁測深(MT)二維正演數(shù)值模擬中，有限差分法(FDM)和有限元法(FEM)是主要的計算方法[1]。FEM的核心基礎(chǔ)是變分原理和加權(quán)余量法，Coggon[2]最先實現(xiàn)FEM的電磁正演模擬，隨后William等[3]、Rijo[4]、Wannamaker等[5]和Franke等[6]分別采用矩形元、三角元、混合元和非結(jié)構(gòu)化三角剖分等實現(xiàn)了FEM的MT二維正演模擬。陳樂壽等[7-8]、胡建德等[9-10]先期實現(xiàn)了有限元MT二維正演模擬； 徐世浙等[11-13]深入研究了FEM及其網(wǎng)格剖分； 陳小斌等[14-15]提出了有限元直接迭代算法，取得了很好的效果。FDM的原理是將連續(xù)微分方程及邊界條件用含有限個未知數(shù)的離散差分方程近似代替，差分方程的解就是微分方程邊值問題的數(shù)值解。Madden等[16]率先采用交錯網(wǎng)格FDM進(jìn)行三維大地電磁正演； 隨后Mackie等[17-20]圍繞交錯采樣的FDM進(jìn)行了大量研究； 譚捍東等[21]在此基礎(chǔ)上引入了雙共軛穩(wěn)定梯度法，提高了正演的精度和效率； 肖騎彬等[22]對不同網(wǎng)格精度、系統(tǒng)方程、邊界條件以及預(yù)條件線性算子等進(jìn)行對比，提出了FDM的MT正演的最優(yōu)方案。

FDM和FEM各有優(yōu)缺點。前者運算簡便易行、迭代收斂穩(wěn)定、計算速度快，但僅適用于地下介質(zhì)形狀規(guī)則的情況；后者不要求規(guī)則的網(wǎng)格剖分，適用于復(fù)雜結(jié)構(gòu)的問題，但其存儲量大，求解速度慢，算法的實用性差。為了解決這一難題，不少學(xué)者研究了多種新型的正演方法，如無網(wǎng)格法[23-24]、多重網(wǎng)格法[25-26]等，取得了一定的效果。四叉樹(Quadtree)和八叉樹(Octree)方法[27]是一類適應(yīng)能力強且高效的網(wǎng)格剖分方法。四叉樹是一種樹狀數(shù)據(jù)結(jié)構(gòu)，將平面模型的外接矩形等分成四個小的矩形,然后遞歸劃分，直至達(dá)到精度要求。四叉樹數(shù)據(jù)結(jié)構(gòu)在計算機領(lǐng)域用于研究圖像分割、遙感圖像處理和地理信息系統(tǒng)，Yerry等[28-29]首先將四叉樹數(shù)據(jù)結(jié)構(gòu)應(yīng)用于網(wǎng)格生成，隨后四叉樹網(wǎng)格的數(shù)值模擬被廣泛應(yīng)用于流體力學(xué)領(lǐng)域[30]。Haber等[31]、Lior 等[32]將八叉樹數(shù)據(jù)結(jié)構(gòu)應(yīng)用于可控源電磁法的數(shù)值模擬，但是用于MT正演的研究很少。

本文將四叉樹網(wǎng)格剖分應(yīng)用于MT二維FDM正演，推導(dǎo)出四叉樹網(wǎng)格的差分公式，并對其計算結(jié)果進(jìn)行了驗證。

1 地電模型的四叉樹網(wǎng)格剖分

大地電磁場是從高空垂直入射到地面的平面電磁波，在地球空間中廣泛存在，其外邊界為無窮遠(yuǎn)，實際正演模擬時需取足夠大的矩形區(qū)域。這個矩形區(qū)域?qū)?yīng)四叉樹的根節(jié)點。四叉樹網(wǎng)格將大地電磁場的計算區(qū)域按如圖1所示的四個象限(NE,NW,SW,SE)進(jìn)行遞歸分割，逐步分解成一系列矩形區(qū)域，直至整個區(qū)域的各個部分達(dá)到預(yù)定的最高分辨率，所分解的最小矩形區(qū)域即是一個具有單一屬性(電性參數(shù))的矩形單元，對應(yīng)四叉樹數(shù)據(jù)結(jié)構(gòu)中的葉節(jié)點，每一個葉節(jié)點與實際模型中的某一矩形區(qū)域?qū)?yīng)。

四叉樹編碼分為顯式四叉樹和線性四叉樹。顯式四叉樹通過父節(jié)點和子節(jié)點之間設(shè)立指針的方式建立整個樹結(jié)構(gòu)，又稱為有指針?biāo)牟鏄?。顯式四叉樹節(jié)點不僅要記錄節(jié)點所對應(yīng)的矩形區(qū)域的信息，還要記錄一個父節(jié)點和四個子節(jié)點的指針，整個四叉樹不僅要記錄葉節(jié)點還要記錄非葉節(jié)點，所以占用儲存空間較多。線性四叉樹只記錄從根節(jié)點到葉節(jié)點的路徑以及葉節(jié)點的屬性，占用儲存空間小。但是顯式四叉樹由于使用指針連接各節(jié)點，所以其鄰域查詢更快捷[33]。在正演計算建立差分方程時要頻繁地查詢各矩形單元的鄰域，所以本文選擇使用顯式四叉樹法。

顯式四叉樹的鄰域查詢采用最近公共祖先方法，這種方法不考慮節(jié)點的位置坐標(biāo)，只需要四叉樹的樹結(jié)構(gòu)。下面以圖2中節(jié)點D鄰域的查詢?yōu)槔榻B顯式四叉樹的鄰域查詢方法。對矩形單元D的四個節(jié)點建立差分方程時，需要獲取與單元D相鄰的六個單元A、B、C、E、F、G的信息，就要通過四叉樹查詢這六個葉節(jié)點。對于節(jié)點D，沿四叉樹向上查詢的到其父節(jié)點，該父節(jié)點含A、B、C、D四個子節(jié)點，且A、B、C均為葉節(jié)點，故A、B、C都為D節(jié)點的鄰域，且他們的層級(距離根節(jié)點的層數(shù))相同。在節(jié)點D的父節(jié)點再查詢上一級父節(jié)點，該節(jié)點包含E、F、G三個葉節(jié)點和節(jié)點D的父節(jié)點，故節(jié)點E、F、G即為節(jié)點D的鄰域，且層級比D低1。若節(jié)點E為非葉節(jié)點，則需要向下查詢節(jié)點E的子節(jié)點，直至找到與節(jié)點D相鄰的某個葉節(jié)點。

在查詢鄰域時，要檢查節(jié)點的鄰域與節(jié)點的層級差。如圖3所示，左圖中存在相鄰單元層級差大于1，要對層級較小的單元進(jìn)行二次劃分，從而得到右圖所示更光滑的網(wǎng)格，以限制離散誤差[34]。

圖1 四叉樹網(wǎng)格剖分示意圖 圖2 四叉樹網(wǎng)格鄰域查找示意圖

圖3 網(wǎng)格的光滑處理

2 四叉樹網(wǎng)格的MT差分公式

在四叉樹網(wǎng)格剖分中，網(wǎng)格節(jié)點除邊界點外，可以分為如圖4所示的五種節(jié)點，其中類型a網(wǎng)格節(jié)點與常規(guī)的有限差分相同，采用五點中心差分即可；其他四種節(jié)點呈T形，需要另行推導(dǎo)差分公式。下面用待定系數(shù)法[35]推導(dǎo)各種節(jié)點的差分公式。

取二維構(gòu)造的走向方向為x，由Helmholtz方程可以推導(dǎo)出二維情況下MT響應(yīng)所滿足的偏微分方程為

TE模式

(1)

TM模式

(2)

圖4 五種節(jié)點類型

2.1 TE模式

對圖4中類型a,可以離散為五點差分公式

E=a0E0+a1E1+a2E2+a3E3+a4E4=0

(3)

分別對圖4a中的電場分量E1～E4在E0處泰勒展開

(4)

把式(4)代入式(3)，整理得到

a2O[(dy1)2]+a3O[(dy2)2]+a4O[(dz2)2]=0

(5)

式(5)與式(1)等價。略去高階余項，則

(6)

解之即可得5點差分格式(圖4中類型a)的系數(shù)

(7)

同樣對于T型節(jié)點(類型b～e)也可以同樣推導(dǎo)出差分公式，例如類型d可以離散為

d0E0+d1E1+d2E2+d3E3+d4E4+d5E5=0

(8)

對E1～E5泰勒展開后，代入式(8)解得

(9)

式中

(10)

2.2 TM模式

將式(2)展開為

(11)

(12)

差分公式推導(dǎo)以類型d為例說明，設(shè)其差分公式為

(13)

各點處的泰勒展開為

(14)

將式(14)代入式(13)，整理得

(15)

式(15)與式(13)等價，可得

(16)

由式(12)～式(16)可解得

(17)

式(17)即為TM模式下節(jié)點類型d的差分方程的系數(shù)。

用上述方法可以推導(dǎo)出TE模式、TM模式各類節(jié)點處的差分方程，在此不贅述。由推導(dǎo)過程可知差分方程具有二階精度，且節(jié)點類型a處TE、TM模式的差分方程與常規(guī)的有限差分方程一致。

通過上述差分公式可以求得各節(jié)點處TE模式下的Ex和TM模式下的Hx，此外還需根據(jù)對應(yīng)的正交輔助場Hy和Ey計算阻抗ZTE、ZTM、視電阻率ρTE、ρTM和相位φTE、φTM。同樣用待定系數(shù)法可以求得Hy和Ey，如對節(jié)點類型a有

(18)

式(18)同樣具有二階精度。為提高地表阻抗的計算精度，可考慮加密阻抗計算點附近的網(wǎng)格。阻抗、視電阻率和相位的計算公式分別為

(19)

(20)

3 算例分析

根據(jù)上述差分公式(式(9)和式(17))，按圖5所示的流程利用Fortran語言編寫了四叉樹有限差分法(Quadtree-FDM，簡稱QFDM)二維MT正演程序，并與FDM和FEM的計算結(jié)果進(jìn)行對比。

首先用一維模型進(jìn)行驗證。設(shè)計一個三層模型，從上至下各層電阻率分別為：ρ1=50Ω·m，ρ2=1000Ω·m，ρ3=10Ω·m，對應(yīng)各層厚度分別為：h1=2600m，h2=400m。對該模型分別用解析法、FEM、FDM和QFDM進(jìn)行正演計算，其中FEM和FDM采用相同的網(wǎng)格剖分(共100行100列，中間部分dy和dz均為100m，共30×30個網(wǎng)格，網(wǎng)格大小按1.2倍率向四周擴展)，QFDM對FDM的網(wǎng)格地表部分局部加密。計算結(jié)果如圖6所示，可以看出在10-3～101Hz頻率范圍內(nèi)，幾種方法計算的視電阻率值和阻抗相位值與一維解析解吻合很好； 101～103Hz的高頻部分，幾種方法均有不同程度的誤差(表1)，且頻率越高誤差越大，其中FEM的誤差最小，QFDM其次，F(xiàn)DM誤差最大。

圖5 四叉樹有限差分法MT正演流程圖

圖6 一維模型不同方法正演結(jié)果

方法極化模式視電阻率相對誤差/%相位相對誤差/%最大最小平均最大最小平均FEMTE0.3580.0010.0870.2680.0010.059TM0.3350.0070.0890.3210.0270.072QFDMTE1.8470.0100.5790.9130.0020.201TM2.4260.0340.8461.8030.0880.349FDMTE2.8610.0110.5951.2680.0010.308TM7.9110.0131.4254.9880.1320.909

圖7 二維理論模型示意圖

用二維理論模型分析不同剖分方法的效果。設(shè)計如圖7的理論模型： 上覆地層的電阻率ρ1=50Ω·m，厚度為30km；下伏低阻層電阻率ρ2=10Ω·m，厚度為70km； 在2km深處有一個二維矩形低阻異常體，沿x方向(長度)無限延伸，寬度為2km,厚度為0.4km，電阻率為ρ3=10Ω·m。其四叉樹網(wǎng)格剖分如圖8所示，可以看出，正演核心區(qū)域網(wǎng)格剖分較密，向四周逐漸變稀，在近地表和異常體邊界處再次加密，符合大地電磁場分布的基本特征。圖9為用QFDM正演得到的視電阻率和相位擬斷面圖(繪圖

范圍為[-4000m，4000m])。同樣用FEM和FDM對相同理論模型進(jìn)行正演計算。FDM采用常規(guī)的網(wǎng)格剖分，其加密區(qū)域與四叉樹網(wǎng)格基本一致，最小網(wǎng)格大小也與QFDM一樣，F(xiàn)EM的網(wǎng)格在FDM的基礎(chǔ)上再次加密，用FEM計算結(jié)果為標(biāo)準(zhǔn)解，在y=-2000～2000m區(qū)間以100m為間隔共設(shè)41個觀測點，頻率范圍從10-3～103Hz按對數(shù)等間隔取61個頻點，正演得到共2501個數(shù)據(jù)，統(tǒng)計FDM和QFDM的計算效率和效果如表2所示。從表2可以看出，QFDM的計算量和計算時間僅約為FDM的60%，但計算精度優(yōu)于FDM。

最后，利用二維半圓柱山谷模型和山脊模型驗證本文算法對地形的適應(yīng)性。Ward[36]給出了二維半圓柱山谷模型的解析響應(yīng)。在電阻率為100Ω·m的均勻半空間中，有一個直徑為100m的半圓柱狀山谷(圖10a)。用QFDM計算0.01Hz、TM極化模式的視電阻率響應(yīng)(圖10b)。從圖中可以看出QFDM計算結(jié)果與解析解基本一致。

方法網(wǎng)格節(jié)點數(shù)線性方程組非零元個數(shù)計算時間s極化模式視電阻率相對誤差/%相位相對誤差/%最大平均最大平均QFDM6948351436.67FDM125026084111.32TE0.594 0.088 0.423 0.055TM0.586 0.088 0.412 0.055TE0.393 0.103 0.170 0.037TM1.294 0.146 1.214 0.144

圖9 二維模型TE(左)、TM(右)模式下的視電阻率(上)和相位(下)擬斷面圖

圖10 半圓柱山谷模型(a)及TM視電阻率曲線(b)

二維山脊模型如圖11上所示，電阻率ρ=100Ω·m，求解區(qū)域(y方向)為-2000～2000m，點距為50m，共81個測點，起伏落差為100m，山脊底部寬2.4km，頻率為10Hz。Wannamaker等[37]采用三角網(wǎng)格、線性插值FEM，徐世浙等[38]采用邊界單元法、劉云等[39]采用自適應(yīng)地形的四邊形網(wǎng)格、雙二次插值FEM，均在10Hz頻點處對此模型做過模擬。本文用QFDM和FEM分別對該模型進(jìn)行計算,兩種方法計算的網(wǎng)格核心區(qū)域如圖11中、圖11下所示，計算結(jié)果如圖12所示。由圖可見，在10Hz頻點處，QFDM與FEM的計算結(jié)果基本一致，但QFDM算法由于采用矩形單元計算，并不能完全貼合地形起伏，只能以階梯狀模擬帶地形的模型，在地表曲率較大的點(如山脊兩邊的山腳和山脊的脊背)，粗糙網(wǎng)格的計算結(jié)果并不平滑，TM極化對地形起伏的響應(yīng)更強烈。關(guān)于這種臺階型地形，陳小斌[40]已有詳細(xì)論述。若要改善這些測點的計算效果，可以加密這些測點附近的網(wǎng)格，再調(diào)整模型，使模型在地表與空氣的接觸面更加平滑。以上兩個模型的正演計算表明QFDM對起伏地形的計算結(jié)果是正確的。

圖11 山脊模型(上)及QFDM(中)和FEM(下)網(wǎng)格剖分方案示意圖

圖12 山脊模型10Hz正演響應(yīng)視電阻率(a)和相位(b)曲線

4 結(jié)論與建議

QFDM相較于FDM有很多優(yōu)勢。主要體現(xiàn)在以下幾個方面。

(1)QFDM由于網(wǎng)格剖分自由度高，能根據(jù)實際地質(zhì)情況加密或抽稀網(wǎng)格，對復(fù)雜地質(zhì)情況適應(yīng)能力更強，能計算的頻帶也更寬。

(2)QFDM在電性突變界面加密網(wǎng)格可以有效減小離散誤差，在地表觀測點附近適當(dāng)加密網(wǎng)格，可以提升主場和輔助場的精度，使正演結(jié)果更加精確；在精度要求相同的情況下，網(wǎng)格剖分合理的QFDM計算量遠(yuǎn)小于FDM。

(3)QFDM相較于FEM，在相同網(wǎng)格剖分的情況下計算精度略有不足，但計算效率更高。若兼顧計算效率與計算精度，采用網(wǎng)格略密集的QFDM是最優(yōu)選擇。

本文采用QFDM算法在保證計算精度的前提下，計算效率比常規(guī)FDM提高約一倍，所以該算法可望推廣到三維，為三維反演提供快速算法。三維模型的剖分需采用八叉樹剖分方法，將是進(jìn)一步研究的內(nèi)容。