粟　琪　戴世坤　趙東東

(①中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙 410083；②中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實驗室，湖南長沙 410083)

1　引言

隨著頻率域可控源電磁法和計算機(jī)技術(shù)的不斷發(fā)展，二維和三維電磁場的數(shù)值模擬及反演得到了很大的重視，而2.5維表現(xiàn)為地球物理模型的二維性和場源的三維性，相比于二維場源模擬更加接近于實際場源勘探的結(jié)果。早期關(guān)于2.5維可控源電磁法的研究大多基于各向同性介質(zhì)，Coggon[1]首先把有限單元法應(yīng)用于2.5維電磁法的模擬中； Stoyer等[2]利用有限差分法進(jìn)行頻率域2.5維磁偶極子源電磁響應(yīng)的正演模擬； Doherty[3]運(yùn)用積分方程法求解2.5維磁偶極子源的電磁響應(yīng)； Unsworth等[4]基于有限元詳細(xì)研究了2.5維電性源激發(fā)產(chǎn)生的電磁場響應(yīng)特征； Torres等[5]提出了一種2.5維數(shù)值模擬的快速算法； Mitsuhata[6]實現(xiàn)了基于總場的可控源電磁法2.5維正演模擬，利用網(wǎng)格單元場函數(shù)插值公式與等參單元有效模擬了地形；Li等[7-9]等基于非結(jié)構(gòu)化三角網(wǎng)格、采用預(yù)制條件的共軛梯度法和準(zhǔn)最小殘差法求解方程，實現(xiàn)了2.5維海洋電磁自適應(yīng)有限元正演模擬； 柳建新等[10]重點(diǎn)研究了基于雙二次插值的頻率域可控源電磁法有限元正演模擬，并對波數(shù)選取進(jìn)行研究，給出了一組兼顧計算速度和模擬精度的波數(shù)； Unsworth等[11]在有限元正演的基礎(chǔ)上提出了一種子空間海洋可控源電磁的二維反演算法，引入二維OCCAM算法并得以成功應(yīng)用；Abubakar等[12]研究了2.5維低頻電磁法正反演，實現(xiàn)了一種快速可靠的反演算法； Key等[13-15]基于自適應(yīng)有限元提出了面向目標(biāo)的2.5維可控源電磁法正演算法，并在此基礎(chǔ)上發(fā)展了Occam的反演算法； 陳光源等[16]實現(xiàn)了2.5維海洋可控源電磁法的共軛梯度反演，探討了2.5維海洋可控源電磁法的不同反演算法，并分析了不同參數(shù)對反演速度及精度的影響。

目前研究各向異性介質(zhì)可控源電磁法數(shù)值模擬和反演計算的文獻(xiàn)并不多。Li等[17,18]研究了均勻半空間各向異性介質(zhì)情況下的CSAMT響應(yīng)，并分析了視電阻率與地下真實電阻率的關(guān)系； Loseth等[19]提出了一種模擬層狀各向異性介質(zhì)中波數(shù)域水平電偶極源產(chǎn)生的電磁場的算法； 羅鳴等[20,21]探討了一維電阻率各向異性對海洋可控源電磁響應(yīng)的影響，并實現(xiàn)了一維垂直各向異性介質(zhì)頻率域海洋可控源電磁資料的反演； Kong等[22]運(yùn)用有限元法進(jìn)行水平各向異性條件下的2.5維海洋可控源電磁法數(shù)值模擬； Li等[23]采用自適應(yīng)有限元法實現(xiàn)了傾斜各向異性情形下基于二次場的二維海洋可控源電磁法的正演算法； 劉穎等[24]詳細(xì)研究了層狀各向異性介質(zhì)中任意取向電偶源的海洋電磁響應(yīng)以及不同發(fā)射姿態(tài)對勘探能力的影響； 趙東東[25]研究了任意各向異性介質(zhì)中2.5維可控源電磁法正演模擬算法及觀測系統(tǒng)的設(shè)計。

本文從各向異性的二維地電模型出發(fā)推導(dǎo)出基于總場的三軸各向異性的耦合微分方程組，避免了二次場計算時繁雜的地形計算問題。文中采用基于非結(jié)構(gòu)化網(wǎng)格的自適應(yīng)有限元法求解正演控制方程組，網(wǎng)格由僅局部加密而不改變插值函數(shù)階次的h型自適應(yīng)策略生成[26]，參考開源程序MARE2DEM的正演部分算法，使用KDTree管理數(shù)據(jù)以提高自適應(yīng)過程中對節(jié)點(diǎn)和單元訪問的效率。此外，對正演算法程序做了并行化處理以提高計算效率，并調(diào)用MKL庫中Pardiso_64位求解器求解大型稀疏復(fù)線性方程組，然后在理論模型算例計算中分別與一維解析解和MARE2DEM計算結(jié)果進(jìn)行對比，驗證了本文正演算法具有較高的計算精度和計算效率。反演時采用對初始模型要求較低的Occam反演算法，并設(shè)計了多個典型的各向異性地電模型進(jìn)行反演計算； 最后通過實測數(shù)據(jù)對本文反演算法進(jìn)行了驗證。

2　方法理論

2.1　正演理論

2.1.1各向異性介質(zhì)的2.5維問題

考慮具有一般性的各向異性的二維導(dǎo)電模型，設(shè)結(jié)構(gòu)走向為y，x方向軸水平向右，z軸垂直向上。假定時諧因子為eiω t，則帶源的頻率域麥克斯韋方程為

(1)

(2)

(3)

(4)

式中:E是電場強(qiáng)度；H是磁場強(qiáng)度；MS為外加源的等效磁化強(qiáng)度;JS為外加源的電流密度; 下標(biāo)“S”代表源；σ是任意各向異性的電導(dǎo)率張量；ε是任意各向異性的介電常數(shù)張量；=iωμ為阻抗率，μ為磁導(dǎo)率。式(1)和式(2)為三維方程，沿構(gòu)造走向 做傅里葉變換得到二維方程。聯(lián)立變換后的偏微分方程，可得到關(guān)于兩個平行于走向的電磁場分量和的耦合微分方程

(5)

(6)

(7)

(8)

其他系數(shù)是與波數(shù)和電導(dǎo)率張量有關(guān)的函數(shù)，具體表達(dá)式如下

(9)

(10)

(11)

(12)

(13)

(14)

2.1.2自適應(yīng)有限元算法

模型的網(wǎng)格剖分對正演模擬結(jié)果的精度有重要影響，而依靠經(jīng)驗進(jìn)行網(wǎng)格剖分僅適用于簡單的地電模型，對于實際勘探中不清楚的地下構(gòu)造復(fù)雜模型，通常難以進(jìn)行合理可靠的網(wǎng)格剖分。自適應(yīng)有限元法能夠自動細(xì)化網(wǎng)格并在不顯著增加計算時間的條件下提供可靠的計算結(jié)果，是一種能夠自動調(diào)整算法以改進(jìn)求解過程的數(shù)值方法。自適應(yīng)有限元法主要包括偏微分方程求解、誤差估計、網(wǎng)格標(biāo)記和網(wǎng)格優(yōu)化四個步驟，而其中研究重點(diǎn)是誤差估計和自適應(yīng)網(wǎng)格優(yōu)化。

(15)

式中：T和C為耦合微分方程式(13)和式(14)的系數(shù)矩陣；f為源項向量。為簡化表述，令

(16)

則式(15)的離散形式可表示為

B(un,v)=F(v)

(17)

式中：u為精確解；un為有限元數(shù)值解。令全局誤差估計算子ξn≈u-un，由于精確解u無法預(yù)知，所以全局誤差估計算子不能直接求解，可通過近似誤差函數(shù)計算

B(ξn,v)=B(u-un)=F(v)-B(un,v)

(18)

由于二維電磁對于沿走向方向的場及其空間梯度的精確度有較高的要求，參考Key等[13]定義的誤差函數(shù)

(19)

式中：e0和e1是兩個調(diào)節(jié)誤差函數(shù)的常數(shù)[13]； Δs表示包含接收點(diǎn)的可能不連續(xù)區(qū)域的三角單元。利用加權(quán)余量法計算得到耦合方程的有限元解un，然后通過式(18)求解全局誤差估計算子ξn， 代入式(19)計算誤差函數(shù)G，將G值代入B(v,wn)=

G(v)，求出B(un,v)的對偶解wn，通過B(v,δn)=G(v)-B(v,wn)得到對偶解的離散誤差δn，最后計算出三角形單元△的局部單元誤差指示因子uΔ=|F(δn)-B(un,δn)|。

根據(jù)計算的誤差因子對網(wǎng)格進(jìn)行恰當(dāng)?shù)膬?yōu)化，加密局部網(wǎng)格以增加有限元解的精度。本文的網(wǎng)格優(yōu)化方法采用的是Delaunay三角約束剖分程序Triangle[27]，該程序使用Ruppert網(wǎng)格優(yōu)化策略，主要是通過增加節(jié)點(diǎn)數(shù)直到所有三角單元滿足權(quán)重約束從而更新和優(yōu)化網(wǎng)格。對于最初產(chǎn)生的粗網(wǎng)格，求解相應(yīng)的有限元解，計算每個單元上的局部誤差指示因子u△，標(biāo)記誤差較大的單元以便進(jìn)行優(yōu)化，產(chǎn)生新的細(xì)化網(wǎng)格并在該網(wǎng)格上重新求解偏微分方程，再次計算u△并標(biāo)記需要細(xì)化的網(wǎng)格單元，重復(fù)上述四個步驟直至誤差泛函滿足設(shè)定的要求，其算法流程如圖1所示。

圖1　自適應(yīng)有限元正演模擬算法流程圖

2.1.3靈敏度計算

靈敏度直接影響反演的精度和效率，目前常用的有Rodi算法、互易原理和伴隨方法[28]。伴隨方法通過設(shè)置合理的伴隨源，計算其產(chǎn)生的伴隨電磁場，對伴隨電場與真實場電場的點(diǎn)積進(jìn)行積分來計算電磁場分量U與模型電導(dǎo)率σj的偏導(dǎo)數(shù)，在此用通式U表示電場或磁場。McGillivray等[29]、Unsworth等[11]、Farquharson等[30]對該方法進(jìn)行了詳細(xì)論述，并將其應(yīng)用于二維電磁問題。劉穎[28]、Key[31]分別將伴隨方法應(yīng)用于各向同性和三軸各向異性介質(zhì)2.5維海洋可控源電磁法正反演計算中。在模型參數(shù)個數(shù)大于接收點(diǎn)數(shù)量時，利用伴隨方法計算效率相對最高。

本文利用上標(biāo)“a”表示引入伴隨源產(chǎn)生的相關(guān)量，對于有限區(qū)域D，通過添加伴隨場，靈敏度可由如下公式計算

(20)

通過設(shè)置恰當(dāng)?shù)陌殡S源即可求出電磁場分量U對電導(dǎo)率的偏導(dǎo)數(shù)

(21)

對于2.5維問題，沿走向y方向無限延伸，對于電導(dǎo)率為σj的單元，積分區(qū)域為Aj，故上式可寫成

(22)

其中

(23)

定義

(24)

將式(23)代入式(22)，并分別對y和ky求積分，將式(24)代入，可得靈敏度計算公式

(25)

式中，yr和ys分別表示沿走向方向的接收點(diǎn)和源的位置。

2.2　反演算法

反演算法采用基于高斯—牛頓法改進(jìn)的Occam算法[32]，基于正則化反演策略的目標(biāo)函數(shù)為

Φ=μ(λ‖Rm‖2+β‖m-m*‖2)+

(26)

對于給定的初始模型m(k)，為使目標(biāo)函數(shù)最小，新的迭代模型m(k+1)為

m(k+1)=[μRTR+(WJ(k))T(WJ(k))]-1×

(27)

(28)

模型粗糙度算子R通過提供對模型變化的衡量促使反演穩(wěn)定，并使目標(biāo)函數(shù)在最小化的過程中避免反演出虛假結(jié)構(gòu)。對于非結(jié)構(gòu)化三角網(wǎng)格，粗糙度范數(shù)取為

(29)

反演迭代過程通過計算模型修改量Δm確定新的模型m(k+1)，計算并判斷模型正演響應(yīng)與觀測數(shù)據(jù)之間的擬合差是否達(dá)到要求，實現(xiàn)對目標(biāo)函數(shù)Φ最小化的求解。

3　正演算法測試

3.1　算法驗證

為驗證正演算法的正確性，設(shè)計如圖2所示的4層VTI介質(zhì)模型，使用本文正演算法計算其正演響應(yīng)，并分別與解析解和MARE2DEM計算結(jié)果進(jìn)行對比。其中：第一層介質(zhì)為電阻率為0.33Ω·m的海水層，厚度為200m;第二層介質(zhì)厚度為800m，

圖2　一維VTI介質(zhì)模型示意圖

電阻率分別為ρx=ρy=1Ω·m,ρz=100Ω·m； 第三層介質(zhì)電阻率為50Ω·m，厚度為200m；最下層介質(zhì)電阻率與第二層一致。發(fā)射源置于150m深度，發(fā)射頻率為0.5Hz， 120個接收機(jī)在水平方向-10～10km范圍內(nèi)均勻分布，位于海底上方0.1m處。將本文計算的電磁響應(yīng)與解析解進(jìn)行對比，結(jié)果如圖3所示。水平各向異性介質(zhì)正演數(shù)值解與解析解吻合得很好，相對誤差均在1%以下。本文算法與MARE2DEM在運(yùn)用自適應(yīng)有限元求解2.5維電磁問題時都是基于總場離散，與MARE2DEM計算結(jié)果進(jìn)行對比，結(jié)果如圖4所示，可見相對誤差相對較小。

圖3　一維VTI介質(zhì)模型響應(yīng)對比圖

圖4　本文算法與MARE2DEM計算結(jié)果對比

3.2　計算效率

參考Li等[23]設(shè)計的各向異性海洋油氣藏模型(圖5)，空氣電阻率取為1.0×1012Ω·m，空氣下方為電阻率0.3Ω·m、深度為1km的海水層。采用電偶極源，發(fā)射頻率為0.25Hz，發(fā)射電流為1A，發(fā)射源位于海底上方50m的海水中，接收裝置均勻放置在海底上方0.1m的海水中。

假設(shè)海底地層圍巖為各向異性介質(zhì)ρw， 其中各主軸電阻率分別取ρwx=ρwy=1.0Ω·m，ρwz=8.0Ω·m，異常體為各向同性的高阻油氣藏位于海底以下1km處，電阻率為ρa(bǔ)=50Ω·m。 針對上述模型進(jìn)行正演，網(wǎng)格采用算法中的自適應(yīng)網(wǎng)格剖分，由于沿海底地形加密測點(diǎn)，模型正演初次優(yōu)化網(wǎng)格剖分為由3328個節(jié)點(diǎn)組成的6552個三角形單元，

圖5　二維各向異性介質(zhì)模型示意圖

剖分網(wǎng)格的中間部分如圖6a所示； 自適應(yīng)網(wǎng)格優(yōu)化共經(jīng)過7次細(xì)化剖分，最終產(chǎn)生了17189個節(jié)點(diǎn)和34267個三角形單元，如圖6b所示?？梢钥闯觯捎跍y點(diǎn)附近空間剖分對正演響應(yīng)的計算精度影響較大，所以在創(chuàng)建正演模型時對測點(diǎn)附近的網(wǎng)格進(jìn)行了加密，使得剖分的初始網(wǎng)格在測點(diǎn)附近已經(jīng)占據(jù)了較多單元。從自適應(yīng)細(xì)化的過程可以看出，經(jīng)過7次網(wǎng)格優(yōu)化，深部的粗網(wǎng)格沒有進(jìn)行太多的加密，而源和測點(diǎn)附近及異常體邊界的網(wǎng)格得到了加密。

快速、高精度的正演是反演的基礎(chǔ)，而線性方程的求解占據(jù)了正演計算的大部分時間，自適應(yīng)有限元法在對網(wǎng)格進(jìn)行優(yōu)化之后需要計算單元誤差因子，相對于傳統(tǒng)有限元其不僅增加了線性方程組的求解次數(shù)，還可能增加線性方程組的維數(shù)。因此，含有大型稀疏矩陣的線性方程組的穩(wěn)定、快速求解是正演的關(guān)鍵。網(wǎng)格的每一次自適應(yīng)優(yōu)化，都會增加線性方程組的行數(shù)。圖7a～圖7c展示的三個稀疏矩陣取自于自適應(yīng)正演過程中線性方程Ax=b中的左端項矩陣A，隨著網(wǎng)格優(yōu)化其矩陣維數(shù)明顯增加。為驗證本文程序在求解大型稀疏矩陣線性方程組時的性能，選取每次網(wǎng)格優(yōu)化后的線性系統(tǒng)，在設(shè)置OpenMP線程為1的情況下，分別調(diào)用并記錄本文程序AFEM2D所用的Pardiso求解部分算法與MARE2DEM的Umfpack求解部分算法的耗時，結(jié)果如圖7d所示?？梢钥闯?，Pardiso在求解大型稀疏矩陣線性方程組時求解速度相對較快，尤其是稀疏矩陣維數(shù)比較大的時候相對于MARE2DEM使用的Umfpack具有一定優(yōu)勢。

圖6　正演自適應(yīng)網(wǎng)格剖分示意圖

4　反演算法測試

4.1　合成數(shù)據(jù)測試

為驗證算法的準(zhǔn)確性和可靠性，設(shè)計如圖8所示的復(fù)雜各向異性介質(zhì)模型，模擬真實構(gòu)造。設(shè)海底為起伏地形，海水電阻率為0.3Ω·m，海底表層為一與海底地形同樣起伏、厚度約為500m的覆蓋層，電阻率為0.7Ω·m； 覆蓋層下方為一正斷層，上盤、下盤介質(zhì)電阻率分別為2.0Ω·m和4.0Ω·m，厚度為1.5～2.0km。上盤中賦存有各向異性的異常體D，其三個主軸方向的電阻率分別為：ρx=80Ω·m、ρy=ρz=2Ω·m，即異常體D僅在x方向表現(xiàn)出電阻率異常。斷層下方存在A、B和C三個電阻率為1000Ω·m的高阻異常體； 基巖電阻率為100Ω·m； 其上方沉積層電阻率為1Ω·m，厚度約為2km。

圖8　二維各向異性海洋模型

反演所用的觀測數(shù)據(jù)采用自適應(yīng)有限元正演響應(yīng)加上4%高斯隨機(jī)噪聲的合成數(shù)據(jù)；正演采用電偶極發(fā)射源，發(fā)射頻率為0.1、0.25、1.0和3.0Hz，發(fā)射源和接收裝置皆均勻分布在-11～11km范圍內(nèi)，發(fā)射源與接收裝置交錯排布。發(fā)射源共65個，位于海底起伏地形上方50m的海水中，間距為346m； 65個接收裝置放置在海底上方0.1m的海水中，間隔也為346m。觀測數(shù)據(jù)包含電磁場分量Ex的幅值和相位，加入噪聲和去除發(fā)射源附近的部分?jǐn)?shù)據(jù)后共包含18288個數(shù)據(jù)點(diǎn)。

反演初始模型包含空氣、海水和海底地層，其中空氣和海水電阻率不參與反演，給定地層的初始電阻率為1.0Ω·m。計算區(qū)域所劃分的單元越多，所需要的計算量越大，因此只對包含異常體的部分采用較細(xì)的剖分網(wǎng)格，其余部分則采用粗網(wǎng)格以減少不必要的計算量。如圖9所示，選取x范圍為-10～10km、海底起伏面到8.5km深的區(qū)域進(jìn)行精細(xì)剖分，共生成25766個三角形單元，其余部分共剖分為2108個三角形單元。反演過程中針對特定的發(fā)射源、測點(diǎn)和頻率，利用自適應(yīng)有限元正演法對正演網(wǎng)格進(jìn)行自適應(yīng)網(wǎng)格優(yōu)化，靠近源和測點(diǎn)位置附近的單元通常做更細(xì)的剖分，以滿足正演計算的精度需求。對反演初始模型進(jìn)行剖分時采用粗細(xì)網(wǎng)格組合的方式，對異常體可能存在的區(qū)域進(jìn)行較精細(xì)的網(wǎng)格剖分，有利于對異常體邊界的界定，而粗網(wǎng)格的使用可以減少不必要的計算量。

圖9　反演初始模型示意圖

設(shè)定目標(biāo)均方根擬合差為1.0，經(jīng)過13次反演迭代，最終反演模型均方根擬合差達(dá)到1.0031。在CPU為Intel i5、內(nèi)存為16G的臺式機(jī)上，反演總耗時為1.688×105s。反演結(jié)果如圖10和圖11所示，其中圖10是各向異性x軸方向的電阻率反演結(jié)果，圖11為各向異性y方向和z方向的電阻率反演結(jié)果。整體來看，淺部反演效果較好，海底下方地層的電阻率和厚度與模型基本一致，斷層兩側(cè)地層也能較好地進(jìn)行區(qū)分，對于基巖上方相對較深的沉積層，其電阻率反演結(jié)果與模型比較接近，但是由于上方高阻異常體的影響，導(dǎo)致基巖起伏界面不太明確，高阻異常體下方呈現(xiàn)相對低阻致使基巖起伏面下延。反演結(jié)果中除了各向異性異常體存在一定的誤差，總體來看兩圖在非各向異性區(qū)域的反演結(jié)果吻合很好，很接近真實模型的地質(zhì)特征。各向同性異常體A、B和C電阻率皆為1000Ω·m，反演結(jié)果中其形態(tài)范圍和中心電阻率與模型都基本接近，但是還存在一定誤差，其中異常體C的電阻率反演結(jié)果相對更好。對于各向異性異常體D來說，僅在x軸方向表現(xiàn)出電阻率異常，而y軸和z軸電阻率相對于圍巖相同，從x軸方向反演結(jié)果看，盡管其電阻率和輪廓與真實模型還存在差距，但都得到了較好的還原，相對于其他主軸電阻率反演結(jié)果能夠很好地識別該異常體。

圖10　電阻率ρx反演結(jié)果

4.2　實測數(shù)據(jù)測試

數(shù)據(jù)源自斯克里普斯海洋學(xué)研究所Weitemey-er等[33,34]的項目，采集于美國俄勒岡州的水合物脊，該區(qū)屬于距俄勒岡州西海岸約80km的近海，是卡斯卡底古陸板塊匯聚邊緣，它以海底蘊(yùn)藏大量甲烷水合物著稱，該地區(qū)流體、氣體噴發(fā)活躍，有巨大的硫氧化物細(xì)菌和蛤蚌組成的群落，存在大量的自生碳酸鹽巖。測線處海水深度約800～1200m。采用一個拖曳式水平電偶極子源進(jìn)行59次發(fā)射，發(fā)射源位于海底上方約100m的海水中，發(fā)射頻率固定為5Hz，接收裝置由線性排布的25個接收機(jī)組成，測點(diǎn)間距約為600m，測點(diǎn)位置如圖12所示。

圖11　電阻率ρy/ρz反演結(jié)果

圖12　測點(diǎn)位置示意圖

根據(jù)測點(diǎn)和發(fā)射源的位置大致建立海底地形，對實測數(shù)據(jù)進(jìn)行帶地形的各向異性反演。由于建立地形過程中的手動誤差，初始模型中測點(diǎn)與海底地層之間距離與實際情況存在差距，因此也可能給反演帶來一定的誤差。本文反演的自適應(yīng)性體現(xiàn)在其所采用的自適應(yīng)有限元正演引擎。反演初始模型網(wǎng)格的粗細(xì)只影響反演結(jié)果中異常體邊界的識別，理論上精細(xì)網(wǎng)格能夠更加精確地限定異常體，但是也會使得初始模型的劃分單元過多而增加計算時間，且反演過程本身是對初始模型各單元電阻率的不斷估計，因此要求反演迭代過程中模型單元數(shù)保持一致。對于實測數(shù)據(jù)，初始模型應(yīng)根據(jù)測區(qū)地質(zhì)情況和勘探任務(wù)建立，本例中初始模型和網(wǎng)格剖分如圖13所示，由于天然氣水合物通常賦存于水深幾百米的沉積層中，所以初始模型僅對淺部進(jìn)行細(xì)化剖分，其余部分采用粗網(wǎng)格。設(shè)定初始模型為三軸各向同性介質(zhì)(ρx=ρy=ρz=1.0Ω·m)，網(wǎng)格共剖分為4732個節(jié)點(diǎn)、27447個三角單元。

反演經(jīng)8次迭代，耗時7.35×104s，反演模型均方根擬合差為6.095，反演結(jié)果如圖14所示。可以看出，三個軸的電阻率反演結(jié)果非常接近，說明該測線處地下水合物應(yīng)該是近似于電性各向同性的。反演結(jié)果中A、B、C和D四處電阻率在對數(shù)化處理后依然存在相對高阻異常，根據(jù)此處地質(zhì)環(huán)境的情況，推斷上述四處應(yīng)該是天然氣水合物或游離氣，與Weitemeyer等[34]基于有限差分法和有限元正演反演的研究結(jié)果基本一致。

圖13　實測數(shù)據(jù)反演初始模型剖分示意圖

圖14　實測數(shù)據(jù)電阻率反演剖面圖

5　結(jié)束語

本文針對各向異性介質(zhì)中的2.5維可控源電磁場正反演問題，從各向異性二維導(dǎo)電模型出發(fā)到三軸各向異性的耦合微分方程的推導(dǎo)，使用基于非結(jié)構(gòu)網(wǎng)格的自適應(yīng)有限元法求解正演控制方程，反演算法采用對初始模型要求較低的基于高斯—牛頓法改進(jìn)的Occam反演算法。對一維VTI介質(zhì)模型和經(jīng)典高阻油氣藏模型進(jìn)行正演計算，利用自適應(yīng)有限元算法進(jìn)行網(wǎng)格優(yōu)化加密，將計算結(jié)果與一維解析解和MARE2DEM結(jié)果進(jìn)行了對比，驗證了計算精度和計算效率；對于賦存有各向異性異常體的復(fù)雜海洋模型，由于非結(jié)構(gòu)化三角單元網(wǎng)格具有模擬地形和復(fù)雜構(gòu)造的優(yōu)勢，除了反演出各向同性介質(zhì)的構(gòu)造，還比較真實地還原了各向異性異常體的形態(tài)、大小和埋深；對實測數(shù)據(jù)給定各向異性初始模型進(jìn)行了反演，與已有文獻(xiàn)的研究結(jié)果對比，證明了反演算法能比較成功地還原模型的主要電性特征。目前本研究還存在一定局限性，僅推導(dǎo)了任意傾斜各向異性介質(zhì)的情況，今后將會把主軸與構(gòu)造走向夾角加入各向異性正、反演算法中。

[1]Coggon J H.Electromagnetic and electrical modeling by the finite element method.Geophysics,1971,36(2):132-155.

[2]Stoyer C H and Greenfield R J.Numerical solutions of the response of a two-dimensional earth to an oscillating magnetic dipole source.Geophysics,1976,41(3):519-530.

[3]Doherty J.EM modeling using surface integral equations.Geophysics,1988,53(6):644-668.

[4]Unsworth J M,Bryan J T,Alan D C.Electromagnetic induction by a finite electric dipole source over a 2-D earth.Geophysics,1993,58(2):198-214.

[5]Torres V C,Habashy T M.Rapid 2.5D forward mode-ling and inversion via a new nonlinear scattering approximation.Radio Science,1994,29(3):1051-1079.

[6]Mistsuhata Y.2-D electromagnetic modeling by finite-element method with a dipole source and topography.Geophysics,2000,65(2):465-475.

[7]Li Y.A finite-element algorithm for electromagnetic induction in two dimensional anisotropic conductivity structures.Geophysical Journal International,2002,148(3):389-401.

[8]Li Y,Key K.2D marine controlled-source electromag-netic modeling (Part 1):An adaptive finite-element algorithm.Geophysics,2007,72(2):WA51-WA62.

[9]Li Y,Constable S.2D marine controlled-source electro-magnetic modeling (Part 2):The effect of bathymetry.Geophysics,2007,72(2):WA63-WA71.

[10]柳建新，湯文武，童孝忠.基于雙二次插值的2.5維FCSEM有限元正演模擬.中南大學(xué)學(xué)報(自然科學(xué)版),2014,45(2):474-482.

Liu Jianxin,Tang Wenwu,Tong Xiaozhong.Finite element forward modeling of 2.5-D FCSEM based on biquadratic interpolation.Journal of Central South University (Science and Technology),2014(2):474-482.

[11]Unsworth M,Oldenburg D.Subspace inversion of electromagnetic data:Application to mid-ocean-ridge exploration.Geophysical Journal International,1995,123(1):161-168.

[12]Abubakar A,Habashy T M,Druskin V L et al.2.5D forward and inverse modeling for interpreting low-frequency electromagneticmeasurements.Geophysics,2008,73(4):F165-F177.

[13]Key K,Ovall J.A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling.Geophysical Journal International,2011,186(1):137-154.

[14]Key K.Marine EM inversion using unstructured grids:a 2D parallel adaptive finite element algorithm.SEG Technical Program Expanded Abstracts,2012，31:1-5.

[15]Myer D,Key K,Constable S.Marine CSEM of the Scarborough gas field,Part 2:2D inversion.Geophy-sics,2015,80(3):E187-E196.

[16]陳光源,杜立彬,景建恩等.2.5維海洋可控源電磁反演算法及影響參數(shù)研究.地球物理學(xué)進(jìn)展,2016,31(4):1796-1802.

Chen Guangyuan,Du Libin,Jing Jian’en et al.2.5D inversion algorithm and parameters of marine controlled-source electromagnetic method.Progress in Geophysics,2016,31(4):1796-1802.

[17]Li X and Pedersen L B.The electromagnetic response of an azimuthally anisotropic half-space.Geophysics,1991,56(9):1462-1473.

[18]Li X and Pedersen L B.Controlled-source tensor magnetotelluric responses of a layered earth with azimu-thal anisotropy.Geophysical Journal International,1992,111(1):91-103.

[19]Loseth L O,Ursin B.Electromagnetic fields in planarly layered anisotropic media.Geophysical Journal International,2007,170(1):44-80.

[20]羅鳴,李予國.一維電阻率各向異性對海洋可控源電磁響應(yīng)的影響研究.地球物理學(xué)報,2015,58(8):2851-2861.

Luo Ming,Li Yuguo.Effects of the electric anisotropy on marine controlled-source electromagnetic responses.Chinese Journal of Geophysics,2015,58(8):2851-2861.

[21]羅鳴,李予國,李剛.一維垂直各向異性介質(zhì)頻率域海洋可控源電磁資料反演方法.地球物理學(xué)報,2016,59(11):4349-4359.

Luo Ming,Li Yuguo,Li Gang.Frequency domain inversion of marine CSEM data in one dimensional vertically anisotropic structures.Chinese Journal of Geophysics,2016,59(11):4349-4359.

[22]Kong F N,Johnstad S E,Rosten T et al.A 2.5D finite element-modeling difference method for marine CSEM modeling in stratified anisotropic media.Geophysics,2008,73(1):F9-F19.

[23]Li Y G,Dai S K.Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures.Geophysical Journal International,2011,185(2):622-636.

[24]劉穎,李予國.層狀各向異性介質(zhì)中任意取向電偶源的海洋電磁響應(yīng).石油地球物理勘探,2015,50(4):755-765.

Liu Ying,Li Yuguo.Marine controlled-source electromagnetic fields of an arbitrary electric dipole over a layered anisotropic medium.OGP,2015,50(4):755-765.

[25]趙東東.各向異性介質(zhì)可控源電磁法2.5維數(shù)值模擬及觀測系統(tǒng)研究[學(xué)位論文].湖南長沙：中南大學(xué)，2016.

Zhao Dongdong.Research on 2.5D Modeling and Observation System for Controlled-Source Electromagnetic Method in Anisotropic Medium [D].Central South University,Changsha,Hu’nan,2016.

[26]Kittur M G,Huston R L,Oswald F B.Mesh refinement in finite element analysis by minimization of the stiffness matrix trace.Computers & Structures,1989,31(6):891-896.

[27]Shewchuk J R.Reprint of delaunay refinement algorithms for triangular mesh generation.Computational Geometry,2014,47(7):741-778.

[28]劉穎.海洋可控源電磁法二維有限元正演及反演[學(xué)位論文].山東青島：中國海洋大學(xué),2014.

Liu Ying.2D Finite Element Modeling and Inversion for Marine Controlled-Source Electromagnetic Fields[D].Chinese Marine University,Qingdao,Shandong,2014.

[29]McGillivray P R,Oldenburg D W,Ellis R G et al.Calculation of sensitivities for the frequency-domain electromagnetic problem.Geophysical Journal International,1994,116(1):1-4.

[30]Farquharson C G,Oldenburg D W.Approximate sensi-

tivities for the electromagnetic inverse problem.Geophysical Journal International,1996,126(1):235-252.

[31]Key K.MARE2DEM:a 2-D inversion code for con-trolled-source electromagnetic and magnetotelluric data.Geophysical Journal International,2016,207(1):571-588.

[32]Constable S C,Parker R L,Constable C G.Occam’sinversion — A practical algorithm for generating smooth models from electromagnetic sounding data.Geophysics,1987,52(1):289-300.

[33]Weitemeyer K,Constable S,Key K et al.First results from a marine controlled-source electromagnetic survey to detect gas hydrates offshore Oregon.Geophysical Research Letters,2006,33(3):L03304.

[34]Weitemeyer K,Gao G,Constable S et al.The practical application of 2D inversion to marine controlled-source electromagnetic data.Geophysics,2010,75(6):F199-F211.