粟　琪　戴世坤　趙東東

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

1　引言

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

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

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

2　方法理論

2.1　正演理論

2.1.1各向異性介質的2.5維問題

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

(1)

(2)

(3)

(4)

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

(5)

(6)

(7)

(8)

其他系數是與波數和電導率張量有關的函數，具體表達式如下

(9)

(10)

(11)

(12)

(13)

(14)

2.1.2自適應有限元算法

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

(15)

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

(16)

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

B(un,v)=F(v)

(17)

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

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

(18)

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

(19)

式中：e0和e1是兩個調節(jié)誤差函數的常數[13]； Δs表示包含接收點的可能不連續(xù)區(qū)域的三角單元。利用加權余量法計算得到耦合方程的有限元解un，然后通過式(18)求解全局誤差估計算子ξn， 代入式(19)計算誤差函數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)|。

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

圖1　自適應有限元正演模擬算法流程圖

2.1.3靈敏度計算

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

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

(20)

通過設置恰當的伴隨源即可求出電磁場分量U對電導率的偏導數

(21)

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

(22)

其中

(23)

定義

(24)

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

(25)

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

2.2　反演算法

反演算法采用基于高斯—牛頓法改進的Occam算法[32]，基于正則化反演策略的目標函數為

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

(26)

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

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

(27)

(28)

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

(29)

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

3　正演算法測試

3.1　算法驗證

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

圖2　一維VTI介質模型示意圖

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

圖3　一維VTI介質模型響應對比圖

圖4　本文算法與MARE2DEM計算結果對比

3.2　計算效率

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

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

圖5　二維各向異性介質模型示意圖

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

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

圖6　正演自適應網格剖分示意圖

4　反演算法測試

4.1　合成數據測試

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

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

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

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

圖9　反演初始模型示意圖

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

圖10　電阻率ρx反演結果

4.2　實測數據測試

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

圖11　電阻率ρy/ρz反演結果

圖12　測點位置示意圖

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

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

圖13　實測數據反演初始模型剖分示意圖

圖14　實測數據電阻率反演剖面圖

5　結束語

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

[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有限元正演模擬.中南大學學報(自然科學版),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維海洋可控源電磁反演算法及影響參數研究.地球物理學進展,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]羅鳴,李予國.一維電阻率各向異性對海洋可控源電磁響應的影響研究.地球物理學報,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]羅鳴,李予國,李剛.一維垂直各向異性介質頻率域海洋可控源電磁資料反演方法.地球物理學報,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]劉穎,李予國.層狀各向異性介質中任意取向電偶源的海洋電磁響應.石油地球物理勘探,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]趙東東.各向異性介質可控源電磁法2.5維數值模擬及觀測系統(tǒng)研究[學位論文].湖南長沙：中南大學，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]劉穎.海洋可控源電磁法二維有限元正演及反演[學位論文].山東青島：中國海洋大學,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.