高 妍 馬 超*③ 張向宇

(①?gòu)V州海洋地質(zhì)調(diào)查局，廣東廣州 511458；②南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(廣州)，廣東廣州 511458；③天然氣水合物勘查開發(fā)國(guó)家工程研究中心，廣東廣州 511458)

0 引言

海洋可控源電磁(MCSEM)法是有效探測(cè)海底電阻率異常體的重要方法之一，在海洋地質(zhì)調(diào)查(海洋底構(gòu)造、洋脊擴(kuò)張、火山運(yùn)動(dòng)等)、海底資源(油氣、水合物、地下水、多金屬硫化物等)探測(cè)等諸多領(lǐng)域得到了成功應(yīng)用[1-7]。該方法的優(yōu)勢(shì)在于：①相較于海洋地震方法，勘探成本低；②電磁響應(yīng)對(duì)異常體的電阻率差異靈敏度高，可有效降低油藏勘探和開發(fā)的風(fēng)險(xiǎn)；③是地震勘探方法的有效補(bǔ)充，如在天然水合物調(diào)查中，對(duì)于地震方法無(wú)法識(shí)別的天然水合物頂界，電磁法可以有效識(shí)別。淺海的可燃冰儲(chǔ)層具有高電阻率特性，利用MCSEM勘探不僅可以估算可燃冰的飽和度分布，也可與地震勘探結(jié)合，優(yōu)勢(shì)互補(bǔ)，顯著降低鉆探風(fēng)險(xiǎn)、提高鉆探成功率[8]。中國(guó)的MCSEM勘探始于2016年在瓊東南海域利用自主研制的海底電磁接收機(jī)開展的可燃冰探測(cè)試驗(yàn)[9]。

MCSEM的正演是反演以及實(shí)際應(yīng)用研究的基礎(chǔ)。由于海底地層的電性結(jié)構(gòu)可能非常復(fù)雜[10]，因此開發(fā)一個(gè)適用于MCSEM復(fù)雜三維模型的穩(wěn)定且有效的正演算法至關(guān)重要。本文提出一種基于預(yù)條件迭代求解的頻率域矢量有限元電磁模擬方法，可高效地進(jìn)行復(fù)雜海底地形條件下的MCSM模型正演模擬。

MCSEM正演模擬的數(shù)值方法主要包括有限差分法、有限元方法、積分方程法、多網(wǎng)格法[11]等，這些方法各有優(yōu)、缺點(diǎn)。有限元方法構(gòu)造的源項(xiàng)采用單元積分方法，相比有限差分的加權(quán)平均方法更精確，因此有限元方法的求解精度理論上更高，并且更適用于復(fù)雜地質(zhì)模型的構(gòu)建。已證明有限元方法是一個(gè)適用于地球物理問(wèn)題中任何復(fù)雜的二維或三維結(jié)構(gòu)的數(shù)值分析的有效工具，其中矢量有限元方法因其算法穩(wěn)定、求解效率高且適合于復(fù)雜幾何形狀模型的模擬，一直是電磁模型領(lǐng)域的研究熱點(diǎn)[12-15]。矢量有限元(比如六面體單元)方法中，單元內(nèi)的電場(chǎng)之間具有更緊密的耦合關(guān)系，且插值函數(shù)具有嚴(yán)格的無(wú)散性。Mitsuhata等[16]提出了一種基于矢量有限元的大地電磁建模算法。Schwarzbach等[17]基于線性高階矢量有限元，利用四面體單元離散對(duì)MCSEM數(shù)據(jù)進(jìn)行建模，精確地再現(xiàn)了真實(shí)海洋深水模型。

三維模型模擬中，曲面的離散可使用普通網(wǎng)格、六面體單元或四面體單元。對(duì)于三維矢量有限元模擬，四面體單元具有結(jié)構(gòu)可任意離散化的特點(diǎn)[18-19]。比較而言，基于六面體單元的矢量有限元?jiǎng)t更易實(shí)現(xiàn)結(jié)構(gòu)化網(wǎng)格[20]、模型加密，具有正反演網(wǎng)格操作簡(jiǎn)單[21]等優(yōu)勢(shì)，集合了有限差分易于離散和矢量有限元對(duì)復(fù)雜模型適應(yīng)性強(qiáng)的優(yōu)勢(shì)，得到了廣泛的應(yīng)用[22-24]。Nam等[25]利用規(guī)則六面體單元的矢量有限元建立了大地電磁三維模型，研究了海洋火山地形與水深對(duì)不同頻段大地電磁數(shù)據(jù)的影響，并進(jìn)行了數(shù)值分析。Kordy等[26]將直接求解法應(yīng)用于不規(guī)則六面體單元的矢量有限元方法，并進(jìn)行了散度校正，結(jié)果證明該方法可用于準(zhǔn)確求解復(fù)雜地形的電磁響應(yīng)。

在有限元求解框架下搭載穩(wěn)定、高效的矩陣可顯著提高求解效率。迭代方法是計(jì)算有限元稀疏矩陣的重要方法，其優(yōu)點(diǎn)是計(jì)算量小、對(duì)計(jì)算機(jī)的存儲(chǔ)要求低[13,27]。預(yù)條件方法在矩陣迭代求解中具有非常重要的作用，最重要的一種預(yù)條件方法是對(duì)矩陣的不完全分解，常見的是不完全LU(ILU)分解?？紤]矩陣對(duì)稱情形，當(dāng)矩陣正定性較強(qiáng)時(shí)，不完全Cholesky分解[28]也是非常有效的預(yù)條件手段。Um等[29]提出電磁場(chǎng)散射公式的有效ILU預(yù)條件的有限元迭代策略，并在多源、多頻、多模型等不同條件下全方位地對(duì)幾種預(yù)條件方法的效率進(jìn)行了評(píng)價(jià)，提出正確選擇預(yù)條件方法可以大大減少?gòu)?fù)雜海底模型的計(jì)算量和計(jì)算時(shí)間。

針對(duì)復(fù)雜電性結(jié)構(gòu)海底地層模型，本文提出一種基于預(yù)條件迭代求解的頻率域矢量有限元電磁正演模擬方法，基于ILU分解的預(yù)條件方法，提高了迭代求解效率。建立層狀海洋電磁模型，通過(guò)模擬電場(chǎng)響應(yīng)，并與解析解進(jìn)行對(duì)比，驗(yàn)證了本文算法的有效性?；诮Y(jié)構(gòu)化任意六面體單元建立了中國(guó)南海某海域真實(shí)地電框架模型，并針對(duì)不同構(gòu)造礦藏儲(chǔ)層特征，分別建立了深埋藏油氣儲(chǔ)層模型、可燃冰礦體模型及復(fù)雜基底地質(zhì)模型，應(yīng)用前述方法計(jì)算電場(chǎng)響應(yīng)，對(duì)不同激發(fā)、接收模式下的MCSEM響應(yīng)特征進(jìn)行分析，并進(jìn)一步總結(jié)出適宜的探測(cè)參數(shù)設(shè)置條件。

1 方法原理

1.1 控制方程

地球物理勘探領(lǐng)域中常用到低頻電磁場(chǎng)需滿足的宏觀電動(dòng)力方程，即Maxwell方程組。在穩(wěn)場(chǎng)近似下，可忽略位移電流。假設(shè)時(shí)間諧變因子為eiωt，頻率域Maxwell方程組可表示為

(1)

式中：E表示空間電場(chǎng)；H表示空間磁場(chǎng)；ω表示圓頻率；磁導(dǎo)率μ=μtμ0=4π×10-7H/m，其中μt和μ0分別表示相對(duì)磁導(dǎo)率和真空中的磁導(dǎo)率；J和M分別表示電流密度和磁流密度；q表示體電荷密度；σ為介質(zhì)電導(dǎo)率；B表示磁感應(yīng)強(qiáng)度；D表示電位移矢量。

由于背景場(chǎng)在源附近急劇變化，用有限元法直接求解總場(chǎng)時(shí)，在邊界處可采用稀疏網(wǎng)格，在場(chǎng)源附近則需采用精細(xì)網(wǎng)格。為了提高計(jì)算效率，可以采用求解異常場(chǎng)的方法求解總場(chǎng)。以電場(chǎng)E為例，首先將其分解為異常場(chǎng)Ea與背景場(chǎng)Eb之和，即E=Ea+Eb，Eb可以通過(guò)求解全空間或半空間的解析解得到，這樣可得到基于電場(chǎng)異常場(chǎng)Ea的二階偏微分方程

?×(μ-1?×Ea)+iωσEa=-iω(σ-σb)Eb

(2)

式中上標(biāo)“b”代表背景。對(duì)使用偶極子源進(jìn)行激發(fā)的情形，采用異常場(chǎng)方法構(gòu)造的源項(xiàng)更加平滑，有利于數(shù)值求解。背景場(chǎng)為均勻全空間或者層狀介質(zhì)的電場(chǎng)分布，可由快速漢克爾變換方法求得[30]。

為了表示的簡(jiǎn)潔性，后文中用E代表電場(chǎng)異常場(chǎng)Ea。假設(shè)求解目標(biāo)在三維長(zhǎng)方體域Ω內(nèi)滿足狄利克雷邊界條件的解，式(2)可以表示為算子的簡(jiǎn)單形式

AE+iωσE=b

(3)

式中：A=?×μ-1?×；b=-iω(σ-σb)Eb。數(shù)值求解得到電場(chǎng)E后，可根據(jù)法拉第定律直接求解磁場(chǎng)

(4)

1.2 頻率域有限元離散

本文使用矢量有限元方法求解式(3)，將矢量形函數(shù)定義在單元的邊上，然后通過(guò)插值表示單元中各點(diǎn)的場(chǎng)值。這里采用六面體單元，六面體單元中節(jié)點(diǎn)和電/磁場(chǎng)的位置如圖1所示。

圖1 六面長(zhǎng)方體單元中節(jié)點(diǎn)和電/磁場(chǎng)位置示意圖

定義單元e沿x、y和z方向的邊長(zhǎng)分別為Δxe、Δye和Δze，單元內(nèi)任一點(diǎn)的電場(chǎng)分量可以表示為

(5)

矢量基函數(shù)在單元邊界也是連續(xù)的，因此基于矢量有限元的公式也適用于無(wú)源條件下的電磁場(chǎng)模擬。采用Galerkin方法[30]，設(shè)殘差函數(shù)(式(2))的弱形式為零

f(E)=?ΩNj·[?×(μ-1?×E)+

iωσE+iω(σ-σb)Eb]dV=0

(6)

式中Nj表示第j個(gè)插值函數(shù)?；谑噶亢愕仁胶偷谝皇噶扛窳侄ɡ?，對(duì)于單元e，式(6)可寫成

iω(σ-σb)MeEb]=0

(7)

式中：ne表示與單元e插值相關(guān)的單元個(gè)數(shù)；Ge和Me分別為單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣，其表達(dá)式為

(8)

將單元?jiǎng)偠染仃嘒e整理到總剛度矩陣G即可得到最后的線性方程組，借助相應(yīng)的矩陣求解器，即可直接求解三維電場(chǎng)E的分布。

通常情況下，僅僅使用長(zhǎng)方體網(wǎng)格無(wú)法精確模擬復(fù)雜地球物理模型[31]。在海洋電磁模擬中，電磁響應(yīng)受地形起伏的影響很明顯，粗糙的網(wǎng)格剖分會(huì)嚴(yán)重影響模擬結(jié)果對(duì)目標(biāo)的分辨率。因此，基于不規(guī)則六面體單元的有限元方法在結(jié)果驗(yàn)證和數(shù)據(jù)補(bǔ)充方面具有重要的意義?？刹捎米鴺?biāo)變換將基于xyz坐標(biāo)系的任意六面體單元(圖2a)映射為基于ξηζ坐標(biāo)系的長(zhǎng)方體單元(圖2b)

圖2 六面體單元從xyz坐標(biāo)系(a)映射到ξηζ坐標(biāo)系(b)

(9)

其中

(10)

為了建立平行于ξ軸的四個(gè)邊的形函數(shù)，首先定義ξ為常數(shù)的面。由于該面垂直于ξ軸，在面上?ξ只有法向分量，這樣通過(guò)?ξ定義的形函數(shù)只在平行于ξ軸的邊上才有非零切向分量。根據(jù)式(10)建立六面體單元的矢量形函數(shù)

(11)

矢量形函數(shù)滿足插值電場(chǎng)在切向分量上的連續(xù)性，因而可有效壓制偽解的產(chǎn)生。對(duì)六面體單元e，與式(8)對(duì)應(yīng)的形式為

(12)

式中J為Jacobian矩陣。

將式(9)和式(10)代入式(12)可得

(13)

六面體單元的矢量形函數(shù)式(11)的旋度為

(14)

確定了形函數(shù)及形函數(shù)的旋度后，式(13)可用三維數(shù)值積分的方法(如八點(diǎn)高斯積分)進(jìn)行求解[20]，得到

(15)

式中：ξi、ηj和ζk分別為高斯積分點(diǎn)；Wi、Wj和Wk分別為加權(quán)系數(shù)。

1.3 基于ILU分解的預(yù)條件迭代方法

ILU分解預(yù)條件方法給大型稀疏方程組的求解帶來(lái)很多新的思路和靈活選擇。比如多重網(wǎng)格方法，在規(guī)則網(wǎng)格的假設(shè)下可進(jìn)行完美的多級(jí)求解[30]，但對(duì)非規(guī)則網(wǎng)格問(wèn)題往往不能達(dá)到理想效果。利用ILU以及改進(jìn)ILU(MILU)可實(shí)現(xiàn)有效的多級(jí)求解。對(duì)系數(shù)矩陣C，定義最簡(jiǎn)單的ILU分解為

C=LU+δ

(16)

式中：δ表示殘差；L和U分別代表下三角和上三角矩陣，且與矩陣C的上三角和下三角部分具有完全相同的稀疏性，即非零元素的結(jié)構(gòu)一致，這種情況被稱為ILU0分解。ILU0分解可通過(guò)按行進(jìn)行的高斯消元法實(shí)現(xiàn)，并舍棄所有零元素對(duì)應(yīng)位置的元素?；谶@樣的定義，ILU0分解對(duì)矩陣的分解精度不高，在一定程度上甚至?xí)茐膶?duì)角線元素的占優(yōu)性。

本文提出一種Jacobian+ILU0的預(yù)條件方法，具體分為兩步：首先，將系數(shù)矩陣C分成上三角和下三角部分，并把上三角部分賦給U；然后，用ILU0分解計(jì)算L，L的對(duì)角線元素均為1。

2 數(shù)值算例

2.1 正演模擬算法有效性驗(yàn)證

2.1.1 精度驗(yàn)證

建立兩個(gè)MCSEM層狀模型，這兩個(gè)模型的區(qū)別在于異常層的厚度和深度不同，通過(guò)模擬這兩個(gè)模型的電磁響應(yīng)分析本文求解方法的精度。

模型(圖3)包括海水層、海底地層、高阻異常層以及基底，兩個(gè)模型的區(qū)別在于高阻異常層的厚度(分別為100、650m)和深度(分別為2000、1200m)。假設(shè)水平電偶極子(HED)源沿x方向置于海底以上50m，源中心點(diǎn)的x坐標(biāo)為0，以100A·m的電偶極矩發(fā)射方波電信號(hào)。接收器置于海底，測(cè)量x方向的電場(chǎng)Ex。下文針對(duì)非異常體和異常體區(qū)域分別對(duì)采用長(zhǎng)方體網(wǎng)格和不規(guī)則網(wǎng)格情況下的電磁模擬精度進(jìn)行對(duì)比。

(1)非異常體區(qū)域

對(duì)圖3a模型的z∈[-1km,1km]分別用長(zhǎng)方體和不規(guī)則矩形網(wǎng)格進(jìn)行劃分，結(jié)果見圖4a和圖4b。設(shè)背景電場(chǎng)響應(yīng)振幅為Eb，解析解為Et，數(shù)值解為Es。根據(jù)下式計(jì)算電場(chǎng)模擬的相對(duì)誤差

圖3 高阻層不同厚度和深度的MCSEM模型

圖4 圖3a模型非異常體區(qū)域網(wǎng)格劃分方案

(17)

0.25Hz的電場(chǎng)分量Ex模擬結(jié)果見圖5。迭代過(guò)程中迭代殘差范數(shù)設(shè)定為1×10-6?？梢钥闯觯邢拊治鼋Y(jié)果與半解析解吻合較好，兩種網(wǎng)格的模擬結(jié)果相對(duì)誤差均小于4.0%，規(guī)則網(wǎng)格剖分和非規(guī)則網(wǎng)格剖分的模擬結(jié)果最大相對(duì)誤差分別為3.8%和2.9%?？梢姡且?guī)則矩形網(wǎng)格劃分情況下的模擬精度略高，尤其是在收發(fā)距較大的情況下。

圖5 圖4兩種網(wǎng)格剖分的MCSEM模型電場(chǎng)響應(yīng)Ex振幅(a)及相對(duì)誤差(b)

(2)異常體區(qū)域

對(duì)圖3b模型中的高阻異常層分別進(jìn)行長(zhǎng)方形網(wǎng)格和不規(guī)則矩形網(wǎng)格剖分(圖6)。0.25Hz時(shí)的Ex振幅模擬結(jié)果見圖7?？梢钥闯?，使用長(zhǎng)方體網(wǎng)格和不規(guī)則矩形網(wǎng)格對(duì)高阻異常層進(jìn)行剖分時(shí)，對(duì)應(yīng)的模擬最大振幅相對(duì)誤差分別為3.7%和2.3%，可見整體上非規(guī)則矩形網(wǎng)格剖分情況下的模擬誤差更小。

圖6 圖3b模型中高阻異常層的長(zhǎng)方體網(wǎng)格剖分(a)和不規(guī)則矩形網(wǎng)格剖分(b)

圖7 圖6兩種網(wǎng)格剖分的模型電場(chǎng)響應(yīng)Ex振幅(a)及相對(duì)誤差(b)

2.1.2 預(yù)條件算法效率分析

利用前文提到的三種預(yù)條件方法(Jacobian預(yù)條件方法、ILU0預(yù)條件方法及ILU0 + Jacobian預(yù)條件方法)進(jìn)行數(shù)值算例測(cè)試，以評(píng)估不同預(yù)條件方法的有效性和可行性。

建立兩種MCSEM模型進(jìn)行預(yù)條件算法算例測(cè)試，即圖8模型(平坦海底地形)和圖9模型(崎嶇海底地形)。將HED源中點(diǎn)置于x=-6.0km、z=1.6km處進(jìn)行激發(fā)，并用有限元方法正演計(jì)算電磁響應(yīng)。對(duì)這兩個(gè)模型分別剖分為長(zhǎng)方體網(wǎng)格和不規(guī)則四邊形網(wǎng)格。

本文采用擬最小殘差法(QMR)迭代求解，不同預(yù)條件方法下的迭代收斂曲線見圖10?？梢钥闯觯翰捎肑acobian預(yù)條件方法的迭代收斂效果最差，平坦海底模型(圖8)經(jīng)過(guò)2000次迭代后的殘差范數(shù)為5.290×10-3，而崎嶇海底地形模型(圖9)經(jīng)3000次迭代后的殘差范數(shù)為3.489×10-4；ILU0預(yù)條件方法較Jacobian預(yù)條件方法的迭代收斂效果略好，圖8模型經(jīng)758次迭代后的殘差范數(shù)為9.931×10-6，而圖9模型經(jīng)796次迭代后的殘差范數(shù)為9.996×10-6；采用Jacobian+ILU0預(yù)條件方法的迭代效果最優(yōu)，圖8模型經(jīng)498次迭代后的殘差范數(shù)為9.884×10-6，圖9模型經(jīng)388次迭代后的殘差范數(shù)為9.761×10-6。

圖10 圖8模型(a)和圖9模型(b)不同預(yù)條件方法下的迭代收斂曲線

與其他兩種預(yù)條件方法相比，本文提出的Jacobian+ILU0預(yù)條件方法對(duì)圖8和圖9模型所需迭代次數(shù)更少，收斂更平滑，尤其對(duì)崎嶇海底地形模型的迭代表現(xiàn)出了較強(qiáng)的適用性。因此，Jacobian+ILU0預(yù)條件方法在MCSEM模型有限元正演模擬中表現(xiàn)出了較好的收斂性能，即便對(duì)崎嶇海底地形模型也表現(xiàn)出較好的模擬收斂效果。

圖8 平坦海底層狀MCSEM模型及網(wǎng)格剖分

圖9 崎嶇海底層狀MCSEM模型及網(wǎng)格剖分

2.2 MCSEM深水模型仿真模擬

受淺水域空氣波的影響，深水模型的異常體較淺水模型的異常體響應(yīng)更明顯[32]，但前提是海底地形是平坦的。當(dāng)海底地形崎嶇時(shí)，深水模型異常體的響應(yīng)特征有待研究，為此，建立深水崎嶇海底模型進(jìn)行仿真模擬?；诠_的中國(guó)南海某海域海底地形數(shù)據(jù)(圖11)，建立不同的三維儲(chǔ)層模型，進(jìn)行MCSEM有限元仿真模擬。

圖11 中國(guó)南海某海域三維海底地形

2.2.1 考慮海底地形的油氣儲(chǔ)層模型

圖12所示三維儲(chǔ)層模型包含海底地形信息，油氣儲(chǔ)層為4km(x)×4km(y)×100m(z)的長(zhǎng)方體，油藏頂面位于海平面以下3km處。HED源(Tx)沿x方向在-10～10km范圍內(nèi)拖曳前行，距離海底的高度始終保持為50m。接收器(Rx)放置于海底，沿x方向均勻布設(shè)，測(cè)量-10～10km范圍內(nèi)的水平電場(chǎng)Ex，源Tx與測(cè)線位于同一xOz平面[33]。

圖12 考慮海底地形的三維MCSEM儲(chǔ)層模型及網(wǎng)格剖分

圖13a為水平電場(chǎng)Ex振幅隨Tx-Rx偏移距(這里指Tx中點(diǎn)在海底的投影點(diǎn)與Rx的x坐標(biāo)之差，后文簡(jiǎn)稱收發(fā)距)的變化曲線，參考背景幅度為直線?？梢钥闯觯谑瞻l(fā)距較小(1～3km)的情況下，異常體的電場(chǎng)分量Ex振幅較大，即異常體所在的平面位置附近Ex振幅出現(xiàn)正異常，但異常幅值總體不大，局部被地形影響抵消?？梢娛瞻l(fā)距在3km以內(nèi)時(shí)，地形對(duì)水平電場(chǎng)的影響較為明顯，即使在儲(chǔ)層埋藏較淺的情形下也難以明確判斷儲(chǔ)層范圍。隨著收發(fā)距的增加(4～6km)，異常體響應(yīng)逐漸變得明顯，由于地形對(duì)電場(chǎng)振幅的影響較小，曲線較為光滑。因此，對(duì)于該模型，在地形變化不太劇烈情況下進(jìn)行儲(chǔ)層識(shí)別的適宜收發(fā)距是4～6km。當(dāng)收發(fā)距進(jìn)一步增大(7～8km)時(shí)，異常體響應(yīng)減弱，曲線劇烈震蕩，這是地形與異常體共同作用的結(jié)果，在這個(gè)收發(fā)距范圍內(nèi)，僅通過(guò)Ex振幅曲線已經(jīng)無(wú)法有效識(shí)別儲(chǔ)層。

以1500m水深的均勻海底地層模型(含海底地形)的電場(chǎng)響應(yīng)作為背景場(chǎng)，計(jì)算Ex振幅與背景場(chǎng)的比值，得到歸一化曲線(圖13b)?？梢姎w一化曲線與圖13a振幅曲線的變化趨勢(shì)一致，但根據(jù)歸一化曲線可較直觀地判斷儲(chǔ)層的平面分布范圍。整體而言，收發(fā)距較小時(shí)的Ex振幅歸一化曲線，尤其是2km時(shí)，曲線異常范圍與高阻油藏的水平范圍基本吻合。隨著收發(fā)距逐步增大，歸一化振幅隨之增大，歸一化曲線異常邊界逐漸模糊。當(dāng)收發(fā)距大于6km時(shí)，歸一化振幅開始變小，同時(shí)地形的影響開始增大，油藏的響應(yīng)逐漸被淹沒(méi)。

與圖13對(duì)應(yīng)的Ex相位及相位差見圖14。由圖可見，相位曲線(圖14a)與振幅曲線(圖13a)特征基本一致。需指出的是，相位對(duì)地形起伏不敏感，尤其是在儲(chǔ)層響應(yīng)最明顯的4～6km收發(fā)距范圍內(nèi)，相位曲線非常平滑。當(dāng)收發(fā)距增大到7～8km時(shí)，在相位曲線上可以同時(shí)觀察到明顯的儲(chǔ)層(低頻信號(hào))和地形(高頻信號(hào))響應(yīng)疊加在一起形成的鋸齒狀響應(yīng)。與歸一化Ex振幅(圖13b)相比，相位差曲線(圖14b)上出現(xiàn)相同特征所對(duì)應(yīng)的收發(fā)距大于1km，且曲線更平滑。

圖13 圖12模型模擬Ex振幅(a)及歸一化振幅(b)曲線

圖14 圖12模型模擬Ex相位(a)及相位差(b)

2.2.2 考慮海底地形和近海底可燃冰分布的油氣儲(chǔ)層模型

中國(guó)南海海域具有豐富的可燃冰資源?？扇急纬捎诘蜏馗邏涵h(huán)境，因此通常分布于近海底區(qū)域，在可燃冰穩(wěn)定區(qū)域之下常常會(huì)有高飽和的伴生天然氣，造成近海底地層電阻率明顯增大?？扇急诤５滓话闶欠蔷鶆蚍植嫉??；趫D12所示模型，在海底增加四個(gè)不連續(xù)分布的可燃冰礦體，礦體沿y方向的延伸范圍與儲(chǔ)層范圍一致(4km)，建立考慮地形和近海底可燃冰分布的油氣儲(chǔ)層模型(圖15)，網(wǎng)格剖分方案同圖12。

圖15 考慮海底地形和近海底可燃冰分布的三維油氣儲(chǔ)層模型及網(wǎng)格剖分

應(yīng)用本文正演方法得到不同收發(fā)距下縱測(cè)線水平電場(chǎng)Ex振幅及歸一化振幅曲線(圖16)。與僅考慮海底地形的情況(圖12)相比，含可燃冰礦體模型的歸一化振幅的最大值增大約2倍，曲線扭曲嚴(yán)重，僅收發(fā)距較小時(shí)所對(duì)應(yīng)的歸一化振幅曲線可較清楚地分辨可燃冰礦體的分布范圍(圖16b中虛線橢圓所示)，但儲(chǔ)層分布范圍難以判斷。據(jù)圖16還可以看出，含海底地形和可燃冰礦體的油藏模型電場(chǎng)響應(yīng)非常復(fù)雜，可燃冰礦體的電場(chǎng)響應(yīng)與海底地形耦合所產(chǎn)生的響應(yīng)疊合在一起，使儲(chǔ)層響應(yīng)難以識(shí)別?？扇急V體分布的非均質(zhì)性越強(qiáng)，響應(yīng)越復(fù)雜。特別地，當(dāng)可燃冰礦體與油氣儲(chǔ)層的水平分布范圍大致相同時(shí)，二者的響應(yīng)會(huì)疊加在一起，對(duì)目標(biāo)儲(chǔ)層的識(shí)別則愈加困難。進(jìn)一步地，若目標(biāo)儲(chǔ)層水平分布范圍較小，要對(duì)其進(jìn)行識(shí)別就更難以實(shí)現(xiàn)。

圖17所示為電場(chǎng)分量Ex相位及相位差曲線。從圖17a所示相位曲線可以看出，可燃冰礦體引起的異常范圍隨著收發(fā)距的增大逐漸減小，當(dāng)收發(fā)距達(dá)到6km時(shí)已無(wú)法觀察到可燃冰礦體引起的明顯擾動(dòng)。收發(fā)距逐漸增大，可燃冰礦體的電磁響應(yīng)與儲(chǔ)層響應(yīng)逐步發(fā)生疊合，因此異常極大值點(diǎn)明顯漸漸向右偏移。從圖17b所示相位差曲線可以看出，收發(fā)距為1km時(shí)，可清楚識(shí)別可燃冰礦體的分布范圍，但可燃冰礦體和儲(chǔ)層響應(yīng)難以區(qū)分；收發(fā)距增加到2km時(shí)，較1km收發(fā)距時(shí)相位差增大，但可燃冰礦體和儲(chǔ)層響應(yīng)仍難以區(qū)分；當(dāng)收發(fā)距逐漸增大至6km時(shí)，可燃冰礦體的響應(yīng)逐步減弱，相位差曲線的中心點(diǎn)逐漸向儲(chǔ)層中心移動(dòng)，可較準(zhǔn)確地判斷儲(chǔ)層的平面位置，這與圖16b揭示的特征一致。因此，若存在可燃冰礦體時(shí)，在保證信噪比的前提下，為更好地識(shí)別油氣藏，應(yīng)盡量采用較大的收發(fā)距。

圖16 圖15模型模擬Ex振幅(a)及歸一化振幅(b)曲線

圖17 圖15模型模擬Ex相位(a)及相位差(b)

2.2.3 考慮海底地形和基底地層起伏的復(fù)雜地質(zhì)模型

利用電磁數(shù)據(jù)可較準(zhǔn)確地解釋電性異常體的橫向(平面)分布范圍，其缺點(diǎn)是縱向分辨率較低，難以精確重建地層模型，而地震勘探的優(yōu)點(diǎn)在于縱向分辨率高，與電磁方法形成互補(bǔ)關(guān)系。因而，開展高精度電磁勘探時(shí)，結(jié)合地震解釋成果對(duì)電性層位進(jìn)行深度標(biāo)定，是提高電磁勘探效果的重要手段。利用地震資料的處理和解釋成果能較精確地確定海底標(biāo)志層，結(jié)合電磁法對(duì)電性異常體的平面分布解釋成果，可建立區(qū)域地球物理框架，建立相應(yīng)的地電模型，為高精度電磁勘探提供基礎(chǔ)數(shù)據(jù)。

根據(jù)中國(guó)海南某海域地震解釋成果和電磁資料，建立一個(gè)三層的地球物理框架，框架包括海底地形及海底地層。基于此框架，分別建立四個(gè)模型M0、M1、M2和M3(圖18)：M0只考慮海底地形和基底；M1考慮海底地形、次高阻層和基底；M2綜合考慮海底地形、次高阻層、基底及高阻油氣儲(chǔ)層；M3考慮海底地形、基底和高阻油氣儲(chǔ)層。

圖18 基于海南某海域地球物理框架建立的電阻率模型

采用2.2.1節(jié)的激發(fā)、接收參數(shù)，這里僅分析源位于-5km時(shí)的情形。圖19是這四個(gè)模型的模擬Ex振幅剖面，頻率為0.25Hz。從圖19b和圖19c可以看出，雖然模型中存在厚度較大的次高阻層，但由于與背景電阻率差異(2倍)較小，Ex振幅曲線并無(wú)明顯異常。從圖19d可以看出，雖然模型M3中油氣儲(chǔ)層的厚度和分布范圍都比模型M1和M2中的次高阻體小得多，但卻引起了較明顯的電場(chǎng)變化，這是由于儲(chǔ)層與圍巖的電阻率差較大(50倍)。

圖19 CSEM模型模擬Ex振幅剖面

以模型M0的模擬數(shù)據(jù)為背景場(chǎng)，對(duì)模型M1、M2、M3的Ex振幅(圖20a)進(jìn)行歸一化，結(jié)果見圖20b。從圖20b所示歸一化曲線上可以看出，收發(fā)距為2km曲線上次高阻體的特征比較明顯(圖20b中紅色虛線框所示)，而油氣儲(chǔ)層的響應(yīng)(圖20b中的綠色虛線框所示)則在收發(fā)距大于5km才可看出明顯異常。這說(shuō)明不同埋藏深度的電阻率異常體所對(duì)應(yīng)的電磁響應(yīng)會(huì)體現(xiàn)在不同收發(fā)距的電場(chǎng)數(shù)據(jù)中。同時(shí)還可以看出，在收發(fā)距大于9km(即對(duì)應(yīng)坐標(biāo)4km位置)時(shí)，油氣儲(chǔ)層的響應(yīng)依然很明顯，但振幅較小，這樣的有用信號(hào)在實(shí)際勘探中容易被噪聲淹沒(méi)。

圖20 模型M0、M1、M2和M3的模擬Ex振幅(a)及歸一化振幅(b)曲線

圖21是模型M0、M1、M2和M3正演電場(chǎng)分量Ex的相位及相位差曲線，其中背景響應(yīng)采用模型M0的相位數(shù)據(jù)。

由圖21可見，相位參數(shù)可更直觀地顯示不同模型電場(chǎng)響應(yīng)的差異。根據(jù)圖21a，四條曲線整體上可分成兩簇，一簇是模型M0和模型M3響應(yīng)曲線，另一簇是模型M1和模型M2響應(yīng)曲線。這兩簇曲線在收發(fā)距為-2.5km處出現(xiàn)分歧，這與淺部次高阻層的位置有關(guān)；隨著收發(fā)距逐漸增至約5km時(shí)，兩簇曲線中的兩支曲線均發(fā)生分離，相位差逐漸增大，這與油氣儲(chǔ)層的平面分布范圍有關(guān)；隨著收發(fā)距進(jìn)一步增大，四條曲線最終在收發(fā)距為10km處趨于一致，表明此處已無(wú)法觀測(cè)到淺部次高阻體或深部高阻油藏引起的電磁異常。圖21b相位差曲線也顯示出與圖21a類似的特征，其中模型M2和M3中油藏產(chǎn)生的相位差分別約為60°和30°，而淺部異常體產(chǎn)生的相位差則很小。

圖21 模型M0、M1、M2和M3模擬Ex相位(a)及相位差(b)曲線

3 結(jié)論

基于預(yù)條件迭代求解的頻率域矢量有限元模擬，本文實(shí)現(xiàn)了結(jié)構(gòu)化任意六面體單元的三維海洋可控源電磁響應(yīng)模擬，通過(guò)與層狀模型的解析解對(duì)比驗(yàn)證了算法的正確性和有效性。本文提出了簡(jiǎn)化的ILU0預(yù)條件方法，有效提高了迭代求解的效率，收斂曲線更平滑，在崎嶇海底地形模型的模擬中表現(xiàn)出較好的適應(yīng)性。

基于中國(guó)南海某海域真實(shí)深水(1600～2000m)地形數(shù)據(jù)，分別建立了常規(guī)油氣儲(chǔ)層模型、考慮可燃冰分布的油氣儲(chǔ)層模型及考慮基底起伏的復(fù)雜地質(zhì)模型。通過(guò)對(duì)這些模型的模擬及分析得出以下認(rèn)識(shí)。

(1)在近收發(fā)距區(qū)域，海底地形對(duì)水平電場(chǎng)的影響較明顯，即使在儲(chǔ)層埋藏較淺的情形下也難以分辨其水平分布范圍；隨著收發(fā)距的增加，儲(chǔ)層響應(yīng)逐漸增大，同時(shí)海底地形對(duì)電場(chǎng)振幅的影響變小，在地形起伏不太劇烈的情況下可識(shí)別儲(chǔ)層的平面分布范圍；當(dāng)收發(fā)距進(jìn)一步增大，異常體響應(yīng)隨之減弱，海底地形與油藏響應(yīng)共同作用，造成電場(chǎng)曲線劇烈震蕩，無(wú)法對(duì)高阻油藏進(jìn)行有效識(shí)別。

(2)若可燃冰礦體埋藏淺，其電場(chǎng)響應(yīng)明顯，但會(huì)與海底地形耦合，因而只在近收發(fā)距區(qū)域能比較清楚地識(shí)別可燃冰礦體的分布范圍。特別地，可燃冰礦體與油氣儲(chǔ)層的水平分布范圍接近，二者的響應(yīng)會(huì)產(chǎn)生重疊和耦合，則需基于更大收發(fā)距的電場(chǎng)數(shù)據(jù)才能識(shí)別出深部高阻油藏。

(3)埋深不同的電阻率異常體，其響應(yīng)特征會(huì)體現(xiàn)在不同收發(fā)距范圍的電場(chǎng)數(shù)據(jù)中。收發(fā)距稍大(這是一個(gè)相對(duì)數(shù)值，不同條件對(duì)應(yīng)的適宜收發(fā)距會(huì)有所不同)條件下的電場(chǎng)數(shù)據(jù)中，深水油氣儲(chǔ)層的響應(yīng)依然明顯，而海底地形的影響則可忽略。