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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

TE模式

(1)

TM模式

(2)

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

2.1 TE模式

對(duì)圖4中類型a,可以離散為五點(diǎn)差分公式

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

(3)

分別對(duì)圖4a中的電場(chǎng)分量E1～E4在E0處泰勒展開(kāi)

(4)

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

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

(5)

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

(6)

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

(7)

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

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

(8)

對(duì)E1～E5泰勒展開(kāi)后，代入式(8)解得

(9)

式中

(10)

2.2 TM模式

將式(2)展開(kāi)為

(11)

(12)

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

(13)

各點(diǎn)處的泰勒展開(kāi)為

(14)

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

(15)

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

(16)

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

(17)

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

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

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

(18)

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

(19)

(20)

3 算例分析

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

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

圖5 四叉樹(shù)有限差分法MT正演流程圖

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

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

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

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

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

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

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

4 結(jié)論與建議

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

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

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

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

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