基于四叉樹網格的MT二維正演

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

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

0 引言

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

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

本文將四叉樹網格剖分應用于MT二維FDM正演，推導出四叉樹網格的差分公式，并對其計算結果進行了驗證。

1 地電模型的四叉樹網格剖分

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

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

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

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

圖1 四叉樹網格剖分示意圖 圖2 四叉樹網格鄰域查找示意圖

圖3 網格的光滑處理

2 四叉樹網格的MT差分公式

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

取二維構造的走向方向為x，由Helmholtz方程可以推導出二維情況下MT響應所滿足的偏微分方程為

TE模式

(1)

TM模式

(2)

圖4 五種節點類型

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)的系數

(7)

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

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

(8)

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

(9)

式中

(10)

2.2 TM模式

將式(2)展開為

(11)

(12)

差分公式推導以類型d為例說明，設其差分公式為

(13)

各點處的泰勒展開為

(14)

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

(15)

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

(16)

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

(17)

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

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

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

(18)

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

(19)

(20)

3 算例分析

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

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

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

圖6 一維模型不同方法正演結果

方法極化模式視電阻率相對誤差/%相位相對誤差/%最大最小平均最大最小平均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 二維理論模型示意圖

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

范圍為[-4000m，4000m])。同樣用FEM和FDM對相同理論模型進行正演計算。FDM采用常規的網格剖分，其加密區域與四叉樹網格基本一致，最小網格大小也與QFDM一樣，FEM的網格在FDM的基礎上再次加密，用FEM計算結果為標準解，在y=-2000～2000m區間以100m為間隔共設41個觀測點，頻率范圍從10-3～103Hz按對數等間隔取61個頻點，正演得到共2501個數據，統計FDM和QFDM的計算效率和效果如表2所示。從表2可以看出，QFDM的計算量和計算時間僅約為FDM的60%，但計算精度優于FDM。

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

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

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

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

4 結論與建議

QFDM相較于FDM有很多優勢。主要體現在以下幾個方面。

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

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

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

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