陳博 汪宏年? 楊守文 殷長春

1) (吉林大學(xué)物理學(xué)院， 計算方法與軟件國際中心， 長春 130012)

2) (吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026)

1 引 言

海洋可控源電磁測量(marine controlled-source electromagnetic measurement， MCSEM)主要應(yīng)用于海底構(gòu)造勘探、油水特征識別和海上油氣儲層定量評價， 以便于降低海上盲鉆率和勘探成本.其典型測量方式是用隨船拖曳移動的低頻(0.1—10 Hz)水平電偶極子發(fā)射天線， 通過布設(shè)在海底的接收機(jī)測量電磁場強(qiáng)度， 并根據(jù)多頻、多源距接收信號變化特征判斷地下電導(dǎo)率空間分布情況[1，2].由于海底地形構(gòu)造復(fù)雜以及地層橫向電阻率分布不均勻， 在海洋可控源電磁采集方案設(shè)計以及海洋電磁資料處理和解釋過程中， 往往需要進(jìn)行大量的數(shù)值模擬、靈敏度計算以及反演成像等工作.因此，近幾十年來， 圍繞著海洋可控源電磁響應(yīng)的正演模擬、靈敏度矩陣計算以及反演成像技術(shù)均得到快速發(fā)展[3，4].

在正演模擬方面， 各種解析和數(shù)值方法在海洋可控源電磁理論中均得到廣泛研究與應(yīng)用， 例如一維水平層狀地層中的解析法[5-8]、三維有限元法[9，10]、三維有限體積法[11-14]、三維積分方程法[15]以及2.5 維混合法[16]等.在反演成像方面， 主要以高斯-牛頓算法或共軛梯度算法等迭代反演算法為主， 原理是通過逐步修改地層電導(dǎo)率參數(shù)實現(xiàn)理論合成數(shù)據(jù)與實際輸入資料的最佳擬合， 其中包括:一維Occam 反演[17]和一維多參數(shù)迭代反演[18-21];以積分方程為基礎(chǔ)的二維、三維迭代反演[22]以及以有限體積法或有限元法為基礎(chǔ)的二維、三維反演成像[23-27]和參數(shù)化反演技術(shù)[28-30]等.

在整個反演成像中均涉及到電磁響應(yīng)顯式靈敏度的計算， 以便確定目標(biāo)函數(shù)下降方向.通過研究不同接收、不同收發(fā)距電磁場靈敏度矩陣的變化特征， 也能夠為接收點位置以及收發(fā)距優(yōu)化提供理論依據(jù)， 并有助于提高數(shù)據(jù)采集效率以及減少反演工作量等.因此， 如何快速有效地計算海洋可控源電磁響應(yīng)的顯式靈敏度已發(fā)展成為一個專門的研究問題[31，32].在電磁反演成像的所有算法中， 除了一維層狀地層和含井眼的軸對稱電導(dǎo)率地層中的Fréchet 導(dǎo)數(shù)可以通過電磁場解析解[17-22]或半解析解[29，30]進(jìn)行有效計算外， 在其他的二維和三維電磁反演技術(shù)中， 電磁響應(yīng)的Fréchet 導(dǎo)數(shù)只能夠通過數(shù)值方法近似計算.目前主要有兩種不同的計算路線: 直接法和伴隨法[31].伴隨法主要利用電磁場伴隨方程的Green 函數(shù)與一階Born 近似， 將電導(dǎo)率攝動產(chǎn)生的電磁場變化表示成體積分的形式，從而得到電磁場微小變化與電導(dǎo)率攝動量間的線性關(guān)系[32-34].在二維和三維反演問題中， 由于伴隨方程的Green 函數(shù)通常沒有解析解， 需要通過數(shù)值計算技術(shù)確定各個接收點上的伴隨Green 函數(shù)，當(dāng)接收點很多時會大大增加Fréchet 導(dǎo)數(shù)的計算工作量.直接法則是通過各種數(shù)值方法求解攝動方程， 建立相應(yīng)的電磁場微小變化與電導(dǎo)率攝動量間的關(guān)系， 目前應(yīng)用最多的方法是對電場離散方程=S 進(jìn)行攝動處理， 并通過迭代法求解攝動方程確定電磁響應(yīng)的Fréchet 導(dǎo)數(shù)[23，24，32]， 顯然隨著攝動電導(dǎo)率單元數(shù)量的不斷增加， 靈敏度矩陣計算的工作量也會線性增加.為了提高三維電磁反演效率， 利用地層電導(dǎo)率空間分布特點， 可以分別采用像素反演或參數(shù)化反演兩種不同方法[28，35]， 因而迫切需要為像素反演或參數(shù)化反演建立一套統(tǒng)一高效的靈敏度矩陣計算技術(shù).

本文將借助電場耦合勢Helmholtz 方程的三維有限體積法， 結(jié)合直接法求解技術(shù)與嚴(yán)格的偏微分方程攝動理論， 為海洋可控源電磁三維像素反演或參數(shù)化反演問題研究建立一套顯式靈敏度矩陣的準(zhǔn)確高效算法.首先， 利用Yee 氏交錯網(wǎng)格節(jié)點和有限體積法對電場耦合勢Helmholtz 方程進(jìn)行離散[10，11]， 建立移動源耦合勢離散方程， 并將直接法求解器與三維線性插值技術(shù)結(jié)合， 給出同時計算多移動電磁場的一種新算法.在此基礎(chǔ)上， 進(jìn)一步針對像素電導(dǎo)率模型(假定電導(dǎo)率空間連續(xù)變化)和分塊常數(shù)電導(dǎo)率模型(參數(shù)化模型)， 將電導(dǎo)率函數(shù)和其攝動量表示成Yee 剖分網(wǎng)格 Vi，j，k上分片塊狀函數(shù)組合形式， 借助一階Born 近似與多移動源電場耦合勢正演結(jié)果， 將電導(dǎo)率攝動產(chǎn)生的一次空間散射電流定量表示為各個剖分單元上散射電流的疊加， 再利用有限體積法對所有散射電流進(jìn)行離散處理， 得到表征電導(dǎo)率攝動量與耦合勢微小變化之間關(guān)系的大型離散方程組.為快速求解攝動方程并建立顯式靈敏度矩陣表達(dá)式， 進(jìn)一步引入插值算子和投影算子確定離散方程的解， 這樣能夠同時計算所有發(fā)射天線在各個接收器上電場與磁場的顯式靈敏度.最后， 針對參數(shù)化電導(dǎo)率模型與像素電導(dǎo)率模型分別給出靈敏度矩陣的計算結(jié)果， 研究考察海洋可控源電磁測量的空間探測特征.

2 基本原理

本節(jié)主要研究電場耦合勢海洋可控源電磁響應(yīng)三維有限體積法離散與直接法求解、插值算子和投影算子、攝動方程離散與求解技術(shù)以及像素電導(dǎo)率模型和分塊常數(shù)電導(dǎo)率模型空間中顯式靈敏度快速算法.

2.1 三維有限體積正演與投影算子

引入滿足庫侖規(guī)范條件 ? ·A=0 的矢勢A 和標(biāo)勢 φ 后， 可以將非均勻各向同性地層中海洋可控源電磁響應(yīng)數(shù)值模擬表示為求解如下耦合勢的Helmholtz 方程(時間變化關(guān)系假定為 eiωt)[11，12]:

空間任意位置 r 處的電場與磁場強(qiáng)度可以應(yīng)用如下公式通過方程(1)中的耦合勢得到:

為用三維有限體積法求解方程(1)， 選擇足夠大的計算區(qū)域Ω 并在區(qū)域外邊界 ? Ω 上選擇簡單的截斷邊界條件×E(r，rS)|?Ω=0 ， 這時相應(yīng)的矢勢與標(biāo)勢截斷邊界條件為

基 于Yee 氏 交 錯 網(wǎng) 格 Vi+1/2，j，k， Vi，j+1/2，k，Vi，j，k+1/2及Vi，j，k(i=1，··· ，Nx; j =1，2，3，··· ，Ny;k =1，2，··· ，Nz)對計算區(qū)域Ω 進(jìn)行剖分， 同時在各個交錯網(wǎng)格的中心分別定義離散矢勢分量和以 及 標(biāo) 勢 φi，j，k，并借助三維有限體積法對方程(1)進(jìn)行離散， 得到如下離散耦合勢 A (r，rS) 和 φ (r，rS) 的線性代數(shù)方程組[12，14]:

海洋可控源電磁勘探往往采用移動源測量方式， 在同一計算區(qū)域Ω 中包含多個不同位置的發(fā)射源 JSδ(r-rS，k)，(k =1， 2， ··· ，K) ， 為提高多發(fā)射源電磁場數(shù)值模擬效率， 將計算區(qū)域Ω 內(nèi)所有不同位置 rS，k發(fā)射源產(chǎn)生的右端項b(JS，rS，k)，(k =1， 2， ··· ，K) 進(jìn) 行合并， 形成一個 N ×K 階的稀疏帶狀矩陣Bˉ =(b(JS，1，rS，1)，b(JS，2，rS，2)，··· ，b(JS，K，rS，K))， 同時將所有發(fā)射源產(chǎn)生的離散耦合勢A 和φ也組合起來， 即 X ˉ =(x(rS，1)，x(rS，2)，··· ，x(rS，K)).因為在同一個計算區(qū)域中， 每個發(fā)射源對應(yīng)的方程(1)左邊的離散矩陣完全相同， 利用方程(4)的離散結(jié)果可以得到所有移動源耦合勢滿足的方程為

利用直接法求解器PARDISO[36]計算出離散矩陣的逆并代入方程(5)， 可以同時計算出所有發(fā)射源離散耦合勢Xˉ =(x(rS，1)， x(rS，2)，··· ，x(rS，K))的數(shù)值解:

將方程(6)代入方程(7)中并整理， 得到由方程(5)右端項直接計算接收點 rR處的電場與磁場強(qiáng)度的轉(zhuǎn)換公式:

其中，

分別稱為電場和磁場投影算子.

在海洋可控源電磁勘探中 L ?K (即接收器個數(shù)遠(yuǎn)少于發(fā)射源個數(shù))， 由方程(9)計算投影算子所花費的時間遠(yuǎn)遠(yuǎn)小于直接應(yīng)用方程(7)計算離散耦合勢的時間， 因此， 通過引入投影算子(9)能夠有效提高整個的正演模擬效率.

2.2 攝動方程求解與顯式靈敏度

為計算海洋可控源電磁響應(yīng)顯式靈敏度矩陣，應(yīng)用攝動原理可以由方程(1)得到電導(dǎo)率微小攝動量 δ σ(r) 引起的耦合勢變化 δ A(r，rS) 和δφ(r，rS)滿足的方程:

圖1 地層模型、像素和分塊電導(dǎo)率模型示意圖 (a)非均質(zhì)地層模型與剖分網(wǎng)格; (b)像素電導(dǎo)率模型; (c)分塊常數(shù)電導(dǎo)率模型Fig.1.Formation model and two different model spaces: (a) Inhomogeneous formation model and grids; (b) pixel model; (c) block model.

以及截斷邊界條件:

方程(10)的右端項δJ(r，rS)=δσ(r)E(r，rS)是電導(dǎo)率微小攝動產(chǎn)生的散射電流密度.對比耦合勢攝動方程(10)與耦合勢正演方程(1)可以看出，除右端項不同外， 兩者其他部分完全相同， 因此，在散射電流密度已知的情況下， 可采用與方程(1)完全相同的方法求解方程(10).

為便于研究散射電流密度對電磁場的影響、有效計算顯示靈敏度矩陣， 假定地層模型由已知電導(dǎo)率為 σS(r) 層狀地層形成的圍巖和多個電導(dǎo)率不同的異常體組成.圖1 是非均質(zhì)地層模型、Yee 剖分網(wǎng)格上像素電導(dǎo)率模型與分塊常數(shù)電導(dǎo)率模型示意圖.其中， 圖1(a)表示電導(dǎo)率為 σS(r) 的圍巖和多個電導(dǎo)率分別為 σBlock，γ(r)，(γ =1， 2，··· ，Γ) 、空間分布范圍為 Vγ，(γ =1，2，··· ，Γ) 的三維塊狀異常體.圖1(b)和圖1(c)分別是像素電導(dǎo)率模型和分塊電導(dǎo)率模型在Yee 氏交錯網(wǎng)格上電導(dǎo)率的空間分布情況， 在像素電導(dǎo)率模型中(圖1(b))， 每個像素單元 Vi，j，k均有一個獨立的電導(dǎo)率， 而在分塊電導(dǎo)率模型中(圖1(c))， 每個塊狀體內(nèi)的所有像素單元的電導(dǎo)率完全相同.下面針對圖1(b)和圖1(c)兩種不同的電導(dǎo)率模型， 分別研究建立方程(10)右端項的離散方法以及相應(yīng)的顯式靈敏度快速算法.

2.3 像素電導(dǎo)率模型中的顯式空間靈敏度

對于圖1(b)中的像素電導(dǎo)率模型， 為簡單起見， 將Yee 氏交錯網(wǎng)格中的整數(shù)網(wǎng)格作為基本像素單元， 并假定圍巖電導(dǎo)率保持不變， 而電導(dǎo)率異常體中的所有基本像素單元由M 個整數(shù)剖分網(wǎng)格組成， 即 Viα，jα，kα，(α=1， 2，··· ，M) ， 其像素電導(dǎo)率為 σiα，jα，kα，(α=1， 2，··· ，M) ， 這時在像素電導(dǎo)率模型中， 各個異常體上電導(dǎo)率攝動量的空間分布可以表示為如下分片常數(shù)函數(shù)形式:

這時， 因電導(dǎo)率攝動產(chǎn)生的散射電流可以看作是如下一系列電偶極子元的疊加:

將方程(14)代入方程(10a)中， 并分別在半整數(shù)單元 Vi+1/2，j，k， Vi，j+1/2，k和 Vi，j，k+1/2上計算其體平均，則得到右端項離散向量的各個分量[12]:

對(15)式的離散結(jié)果進(jìn)行適當(dāng)整理并將其右端項的離散向量表示成矩陣的乘積形式:， 其 中， h 是 對 應(yīng) 單 元 的 寬 度，δν=(δνi1，j1，k1， δνi2，j2，k2， ··· ， δνiM，jM，kM)T是M 個像素電導(dǎo)率相對攝動量組成的M 維列向量，是 (15)式中的離散系數(shù)經(jīng)過適當(dāng)組合后得到的N × M 階矩陣.這時， 方程(10)在(11)式的邊界條件下的離散結(jié)果可表示為如下線性方程:

這里， δ X(rS) 是像素電導(dǎo)率相對攝動向量 δ ν 引起的離散耦合勢微小變化， 為N × M 階未知矩陣， 而是方程(10)左端的離散矩陣， 并且與方程(5)中的離散矩陣完全相同.利用方程(9)中的電場和磁場投影算子， 從方程(16)可以得到電磁場強(qiáng)度微小變化與像素電導(dǎo)率相對攝動向量間的線性關(guān)系:

這里

是 3 ×M 階復(fù)矩陣， 表示 rS處的發(fā)射天線在接收點rR處電場強(qiáng)度與磁場強(qiáng)度像素顯式靈敏度矩陣， 靈敏度矩陣中每個元素反映出單元 Viα，jα，kα中電導(dǎo)率相對變化量 δ νiα，jα，kα，(α=1， 2， ··· ，M) 對電磁場各個分量的定量影響.

這里特別需要指出的是， 在進(jìn)行像素靈敏度矩陣計算或進(jìn)行像素反演時， 像素單元總數(shù)M 可能達(dá)到十萬次以上的量級， 而網(wǎng)格節(jié)點上未知離散矢勢和標(biāo)勢總數(shù)N 達(dá)到百萬次的量級.例如若像素總數(shù) M =40×40×56=89600 個， 離散矢勢和標(biāo)勢總數(shù) N =1857844.如果直接從線性化方程(16)出發(fā)， 確定 δ X(rS) 和像素電導(dǎo)率相對攝動向量 δν 間的顯式線性關(guān)系再通過插值方法計算各個接收器上電場和磁場顯式靈敏度， 則必須先計算出N × M 階的大型復(fù)矩陣和N × M 階復(fù)矩陣的 乘積計算工作量是十分巨大的， 消耗的CPU 時間往往也是難以承受的.而通過方程(9)的投影算子以及(18)式計算電場和磁場的顯式靈敏度， 只需要進(jìn)行3 × M 階復(fù)矩陣與離散右端項N × M 階復(fù)矩陣的乘積， 大大減少了中間的計算過程.此外， 需要說明的是投影算子只與接收點的位置有關(guān)， 因此， 可以根據(jù)接收點位置事先計算出并反復(fù)使用， 而復(fù)矩陣只需要利用各個像素單元上電流密度值通過方程(15)就能夠快速確定， 所以在正演過程中只需增加一項復(fù)矩陣的計算， 然后通過方程(8)和方程(18)就能夠同時得到各個測點上的電磁場與電磁場的顯示靈敏度， 這種顯式靈敏度的快速算法為三維非線性反演建立了非常有利的研究基礎(chǔ).

2.4 分塊常數(shù)電導(dǎo)率模型中的顯式靈敏度

對于圖1(c)中分塊常數(shù)電導(dǎo)率模型， 各個塊狀異常體 Vγ，(γ =1，2，··· ，Γ) 內(nèi)部的電導(dǎo)率假定為常 數(shù) σBlock，γ， 其電導(dǎo)率相對攝動量 為δνBlock，γ=δσBlock，γ/σBlock，γ， 這時由于各個異常體電導(dǎo)率攝動導(dǎo)致的地層電導(dǎo)率相對攝動量的空間分布也可表示分片常數(shù)函數(shù)形式:

這里 Ξ (r，Vγ) 是塊狀異常體 Vγ上的特征函數(shù).此外， 為了便于計算， 因為塊狀異常體電導(dǎo)率攝動對電磁場的影響， 假定異常體 Vγ中包含的像素單元為這時方程(10)右端的散射源可近似表示為

采用與方程(14)類似的離散方法， 可以得到塊狀體模型中右端項的離散公式:

與方程(14)右端項類似， 方程(21)的離散結(jié)果也可表示成矩陣形式:其中， δ νBlock=(δνBlock，1， δνBlock，2， ··· ， δνBlock，Γ)T是Γ 個塊狀異常體上電導(dǎo)率相對攝動量組成的Γ 階列向量，是 (21)式中的離散系數(shù)組合形成的 N ×Γ 階已知矩陣.因此， 分塊常數(shù)電導(dǎo)率模型中各個異常體電導(dǎo)率相對攝動對離散耦合勢微小變化影響也可以表示為線性代數(shù)方程:

其中， δ XBlock(rS) 是 N ×Γ 階未知數(shù)矩陣， 其每個列向量對應(yīng)于單個塊狀體電導(dǎo)率相對攝動量引起的離散耦合勢的微小變化.采用與方程(17)和方程(18)類似的方法， 引入分塊常數(shù)電導(dǎo)率模型中顯式靈敏度矩陣( 3 ×Γ 階復(fù)矩陣):

則得到電磁場強(qiáng)度微小變化與塊狀體模型中各個異常電導(dǎo)率的微小相對攝動量間的線性關(guān)系:

2.5 電磁場振幅與相位顯式靈敏度

在海洋電磁測量中， 往往以電磁場的振幅與相位作為輸出結(jié)果， 如x 方向的電場強(qiáng)度常表示為其中AEx(rR，rS)和 θEx(rR，rS) 分 別是 Ex(rR，rS) 的振幅與相位， 可證明電場強(qiáng)度微小變化與其振幅和相位微小變化間的關(guān)系為進(jìn)一步結(jié)合方程(17)， 可以得到像素電導(dǎo)率模型空間中電場強(qiáng)度 Ex(rR，rS) 的振幅與相位顯式靈敏度矩陣:

以及像素電導(dǎo)率攝動模型空間中電場強(qiáng)度Ex(rR，rS)振幅與相位的線性化響應(yīng):

對于磁場強(qiáng)度 Hy(rR，rS) 的振幅與相位的顯式靈敏度矩陣與線性化響應(yīng)， 以及分塊模型中電磁場振幅與相位的顯式靈敏度矩陣與線性化響應(yīng)也可以用類似方法得到.

3 數(shù)值結(jié)果

本節(jié)首先將前面的高效直接求解方法(efficient direct method， EDM)得到的顯式靈敏度與有限差分方法(finite different method， FDM)[31]得到的靈敏度進(jìn)行對比， 驗證靈敏度快速算法的有效性.在此基礎(chǔ)上利用數(shù)值模擬結(jié)果進(jìn)一步研究考察不同收發(fā)距情況下分塊常數(shù)電導(dǎo)率模型空間與像素電導(dǎo)率模型空間中電磁場水平分量 Ex和 Hy的振幅和相位靈敏度， 同時也將進(jìn)行主測線與旁測線上 Ex振幅相位靈敏度的對比.在下面的數(shù)值結(jié)果中， 工作頻率為0.25 Hz， 發(fā)射天線離海底的距離為50 m， 發(fā)射源間距為100 m.模型電導(dǎo)率是各向同性的， 包含空氣層、海水、沉積層以及異常體， 空氣層電阻率為 1 06Ω·m， 海水層厚度為1 km、電阻率為0.3 Ω·m， 沉積層電阻率為1 Ω·m.整個計算區(qū)域大小為(100， 100， 100) km， 并且x， y 和z 三個正交方向的離散網(wǎng)格數(shù)分別為74， 74 和84 個，離散矢勢和標(biāo)勢總數(shù) N =1857844 ， 在整個考察區(qū)域內(nèi)采用均勻網(wǎng)格， 而考察區(qū)域外到計算區(qū)域邊界間的網(wǎng)格間距按對數(shù)增加.3 個接收器R1， R2 和R3 分別位于x = —6， —4， —2 km 的海床面上， 發(fā)射源在x 軸上的變化范圍是—8—8 km， 間距250 m，每條測線有65 個不同位置的發(fā)射源， 天線的電偶極矩假定等于1 A·m， 圖2 是y = 0 的垂直截面上等間距整數(shù)網(wǎng)格、接收點和發(fā)射源位置示意圖.其中， 異常體在x， y 和z 方向的長度分別為3， 3和0.3 km， 其中心位置為 (0， 0， 1) km， 電阻率為100 Ω·m， 3 個接收點位置固定不動.下面的所有數(shù)值結(jié)果均是在惠普Z820 工作站上運行的， 內(nèi)存配置為256 Gb.

3.1 靈敏度快速算法驗證

首先對圖2 中層狀背景電導(dǎo)率中單一異常體模型， 利用EDM 和FDM 兩種方法分別計算電場振幅和相位相對于整個異常體的靈敏度.圖3(a)和圖3(b)是3 個不同位置接收點上電場水平分量Ex振幅與偏移距(magnitude versus offset， MVO)正演曲線和相位與偏移距(phase versus offset，PVO)正演曲線.可以看出， 在3 個接收器位置， 因為源的奇異性， 振幅響應(yīng)強(qiáng)度最大， 由于電磁場的傳播效應(yīng)， 電場振幅隨源距的增加快速衰減， 相位隨源距增加也呈現(xiàn)出快速變化.然而， 在高阻異常體的上方及右邊附近， MVO 衰減曲線產(chǎn)生輕微的波動效應(yīng)， 距離異常體更近的接收線圈(R3)上的波動效應(yīng)更加明顯， 此外PVO 衰減曲線也顯示出類似的波動效應(yīng)， 這些都是因為高阻異常體的影響.

圖2 三維MCSEM 檢驗?zāi)Ｐ团c網(wǎng)格剖分示意圖(y = 0 截面)Fig.2.Verification model and mesh grid of three-dimensional MCSEM (cross-section of y = 0).

圖3 3 個不同位置接收器上電場 E x 分量的MVO 和PVO 曲線以及由EDM 和FDM 計算得到的振幅與相位靈敏度對比圖 (a)Ex的MVO 正演結(jié)果; (b) E x 的PVO 正演結(jié)果; (c) E x 振幅靈敏度對比; (d) E x 相位靈敏度對比Fig.3.The E x magnitudes (MVO) and phases (PVO) of 3 receivers at different positions and comparison of Fréchet sensitivity obtained by two different algorithms of EDM and FDM: (a) MVO of E x ; (b) PVO of E x ; (c) comparison of E x amplitude sensitivity; (d) comparison of E x phase sensitivity.

為了定量考察高阻異常體微小攝動引起的水平電場振幅以及相位變化， 同時用EDM和FDM 計算了塊狀體模型中水平電場振幅顯式靈敏度與相位差顯式靈敏度這 里 的表 示 異 常體電導(dǎo)率的相對攝動量.圖3(c)和圖3(d)分別是EDM 和FDM 計算得到的MVO 和PVO 靈敏度變化曲線.其中， EDM 得到的靈敏度結(jié)果由3 種不同類型的線表示， FDM 得到的靈敏度結(jié)果由3 種不同類型的離散符號表示.從圖3(c)和圖3(d)可以看出， 兩種不同方法得到的靈敏度曲線符合得非常好， 兩種方法最大相對誤差不超過2%， 從而證明了EDM 在計算顯式靈敏度方面的有效性.此外， 從顯式靈敏度曲線可以清楚看出接收點位置與發(fā)射天線的偏移距對水平電場靈敏度的影響非常明顯， 例如R1 接收器由于離異常體的水平距離最遠(yuǎn)， 其靈敏度值明顯小于另兩個離異常體較近的接收器R2 和R3 上水平電場靈敏度.此外， 從靈敏度曲線的大小可以看出， 發(fā)射天線位于異常體中心右上方0—5 km 范圍內(nèi)， 靈敏度最大， 應(yīng)用這個區(qū)段的測量結(jié)果進(jìn)行異常體反演或?qū)Ξ惓ｓw電導(dǎo)率變化進(jìn)行檢測， 效果會更好.

圖4 3 個不同位置接收器上磁場 H y 分量MVO 和PVO 曲線以及由EDM 和FDM 計算得到的振幅與相位靈敏度對比圖 (a)Hy的MVO 正演結(jié)果; (b) H y 的PVO 正演結(jié)果; (c) H y 振幅靈敏度對比; (d) H y 相位靈敏度對比Fig.4.The H y magnitudes (MVO) and phases (PVO) of 3 receivers at different position and comparison of Fréchet sensitivity obtained by two different algorithms of EDM and FDM: (a) MVO of H y ; (b) PVO of H y ; (c) comparison of H y amplitude sensitivity; (d) comparison of H y phase sensitivity.

圖4 (a)和圖4(b)是與圖3 相同位置接收點上磁場水平分量的 MVO 和PVO 正演曲線.與圖3 中電場振幅與相位的變化特征類似， 隨著源距增加磁場振幅快速衰減， 同時， 在異常體上方以及右邊附近， 可以看到高阻異常體對MVO 與PVO的影響.為了定量考察高阻異常體微小攝動引起的水平磁場振幅與相位變化， 同時用EDM 和FDM 計算了塊狀體模型中水平磁場振幅顯式靈敏度與相位的顯式靈敏度圖4(c)和圖4(d)是兩種不同方法計算得到的MVO 和PVO 顯式靈敏度正演曲線對比圖.其中， EDM 得到的靈敏度結(jié)果由3 種不同類型的線表示， 而FDM 得到的靈敏度結(jié)果由3 種不同類型的離散符號表示.從圖4(c)和圖4(d)可以看出， 兩種不同方法得到的靈敏度曲線符合得非常好， 兩者最大相對誤差也不超過2%， 從而進(jìn)一步證明了EDM 在計算顯式靈敏度方面的有效性.此外， 從顯式靈敏度曲線可以清楚看出， 接收點位置與發(fā)射天線的偏移距對水平磁場靈敏度的影響也非常明顯， 例如R1 接收器由于離異常體的水平距離最遠(yuǎn)， 其靈敏度值明顯小于另兩個離異常體較近的接收器R2 和R3 上水平磁場靈敏度.從靈敏度曲線的大小可以看出， 發(fā)射天線位于異常體中心右上方0—5 km 范圍內(nèi)， 靈敏度最大.對比電場的靈敏度曲線可以看出， 在有效的偏移距范圍內(nèi)(6—10 km)， 磁場的振幅和相位靈敏度高于電場，電場和磁場的相位包含的靈敏度信息均稍多于振幅， 綜合應(yīng)用磁場和相位測量結(jié)果進(jìn)行異常體反演或?qū)Ξ惓ｓw電導(dǎo)率變化進(jìn)行檢測， 效果會更好.

表1 列出了利用EDM 和FDM 兩種算法計算靈敏度矩陣的計算耗時， 表中第1 和第2 行為3 個接收器時的兩種算法的計算耗時， 第3 和第4 行為5 個接收器情況下的計算耗時對比.結(jié)果顯示，隨著發(fā)射器接收器的增加， FDM 的耗時明顯增加，而EDM 得益于投影算子的采用， 計算耗時幾乎沒有增加， 所以可以快速計算多源多接收器的結(jié)果，體現(xiàn)了算法的高效性.

表1 兩種算法的計算耗時Table 1.Computational time cost of EDM and FDM.

3.2 多塊狀體模型的顯式靈敏度

為進(jìn)一步考察多個塊狀體情況下， 塊狀體水平位置以及埋藏深度等對顯式靈敏度的影響， 假定地層中存在3 個幾何尺寸與電阻率與圖2 完全相同的異常體， 中心位置分別為(—3.5， 0， 0.7)， (0， 0， 1.0)和(3.5， 0， 1.3) km.圖5 是含3 個塊狀異常體模型在y = 0 垂直截面上的等效電導(dǎo)率分布、網(wǎng)格剖分以及發(fā)射與接收點位置分布圖.為簡潔起見， 這里僅給出位于R2 位置處的接收器上水平電場 Ex和水平磁場 Hy的振幅、相位分別相對于3 個塊狀體電導(dǎo)率靈敏度的數(shù)值結(jié)果.

圖6(a)和圖6(b)是位于R2 處的接收器電場水平分量 Ex的MVO 和PVO 正演曲線.可以看出， 電場振幅隨源距的增加而快速衰減， 然而， 在異常體上方以及異常體右端， 由于3 個高阻異常體的作用， MVO 衰減速度變慢， 從PVO 曲線也可以看到3 個高阻異常體的影響.圖6(c)和圖6(d)是用EDM 計算得到的MVO 和PVO 分別相對于3 個塊狀體電導(dǎo)率(Block 1， Block 2 和Block 3)的顯式靈敏度曲線.可以清楚地看出， 接收點位置與發(fā)射天線的偏移距以及異常體中心深度對水平電場靈敏度影響均非常明顯， 異常體Block 1 離接收點最近且中心深度最淺， 相應(yīng)的靈敏度值最大，而異常體Block 3 離接收點最遠(yuǎn)且中心深度最深，相應(yīng)的靈敏度值最小.此外， 從靈敏度曲線還可以看出， 對于每個異常體的靈敏度曲線均存在靈敏度絕對值的最大值， 說明在海洋可控源探測中存在著最佳測量位置， 對應(yīng)的靈敏度最強(qiáng).

圖5 3 個塊狀異常體模型在y = 0 垂直截面上的網(wǎng)格剖分與等效電導(dǎo)率分布圖Fig.5.Conductivity distribution and mesh grid of computation model including three different blocks at cross-section of y = 0.

圖6 位于R2 處的接收器上 E x 的MVO 和PVO 曲線以及振幅和相位相對于3 個塊狀體電導(dǎo)率靈敏度 (a) E x 振幅; (b) E x 相位; (c) E x 振幅分別對Block 1， Block 2， Block 3 的靈敏度; (d) E x 相位分別對Block 1， Block 2， Block 3 的靈敏度響應(yīng)Fig.6.Magnitudes (MVO) and phases (PVO) versus offset of E x at receiver R2 and Fréchet sensitivity about three different Blocks: (a) MVO of E x ; (b) PVO of E x ; (c) sensitivity of E x amplitude about Block 1， Block 2， Block 3; (d) sensitivity ofEx phase about Block 1， Block 2， Block 3.

圖7 (a)和圖7(b)是位于R2 處的接收器磁場水平分量 Hy的MVO 和PVO 正演曲線.可以看出， 磁場振幅隨源距的增加而快速衰減， 然而， 在異常體上方以及異常體右端， 由于3 個高阻異常體的作用， MVO 衰減速度變慢， 從PVO 曲線也可以看到3 個高阻異常體的影響.圖7(c)和圖7(d)是用EDM 計算得到的MVO 和PVO 分別相對于3 個塊狀體電導(dǎo)率(Block 1， Block 2 和Block 3)的顯式靈敏度曲線.同樣可以看出， 接收點位置與發(fā)射天線的偏移距以及異常體中心深度對水平磁場靈敏度影響均非常明顯， 但有別于電場靈敏度結(jié)果， 異常體Block 1 雖然離接收點最近且中心深度最淺， 相應(yīng)的磁場靈敏度也是最早出現(xiàn)， 但其靈敏度響應(yīng)峰值小于異常體Block 2， 而異常體Block 3雖然離接收點最遠(yuǎn)且中心深度最深， 其磁場響應(yīng)靈敏度值最小， 但峰值非常接近Block 1， 明顯好于對應(yīng)的電場靈敏度結(jié)果， 說明磁場的最佳探測深度和探測范圍要大于電場.同時， 綜合對比電場和磁場的數(shù)值模擬結(jié)果可以發(fā)現(xiàn)， 電場和磁場靈敏度的異常響應(yīng)范圍會隨著異常體埋深、收發(fā)距的增加而逐漸減小， 磁場靈敏度的值對異常體的響應(yīng)要大于電場靈敏度的值， 磁場靈敏度曲線的峰值明顯大于電場靈敏度峰值， 但靈敏度曲線大于零的持續(xù)范圍比電場靈敏度的范圍要小.由于磁場比電場的靈敏度更大， 更快速達(dá)到最大峰值， 因此選用磁場進(jìn)行電導(dǎo)率動態(tài)監(jiān)測效果可能會更好.此外， 從靈敏度曲線還可以看出， 對于每個異常體的靈敏度曲線均存在靈敏度絕對值最大值， 說明在海洋可控源探測中存在著最佳的測量位置， 對應(yīng)的靈敏度最強(qiáng).

圖7 位于R2 處接收器上磁場 H y 分量的MVO 和PVO 曲線以及振幅和相位相對于3 個塊狀體電導(dǎo)率靈敏度 (a) H y 的MVO曲 線; (b) H y 的PVO 曲 線; (c) H y 振 幅對Block 1， Block 2， Block 3 的靈敏度; (d) H y 相位對Block 1， Block 2， Block 3 的靈敏 度Fig.7.The magnitudes (MVO) and phases (PVO) versus offset of H y at receiver R2 and Fréchet sensitivity about three different Blocks: (a) MVO of H y ; (b) PVO of H y ; (c) sensitivity of H y amplitude about Block 1， Block 2， Block 3; (d) sensitivity ofHy phase about Block 1， Block 2， Block 3.

3.3 像素電導(dǎo)率模型空間中的靈敏度

下面進(jìn)一步考察單一異常體模型上像素靈敏度的空間分布， 了解不同位置的像素電導(dǎo)率相對攝動對電磁響應(yīng)的影響.限于篇幅， 這里僅給出水平分量 Ex的相關(guān)結(jié)果.圖8 給出了海洋地電模型和像素變化范圍， 其中圖8(a)是主測線(y = 0)垂直截面上， 接收器、發(fā)射源、異常體以及x 和z 方向上的像素分布， 圖8(b)是像素空間、異常體以及3 條測線在xy 平面上的投影.x 和y 方向上像素數(shù)均等于40 且尺寸均為250 m， 而z 方向的像素個數(shù)是56， 像素高度為50 m， 整個像素在x， y 和z三個方向上空間變化范圍是(-5， 5)×(-5， 5)×(0.3， 3) km， 像素總數(shù)M =40×40×56=89600個.其中主測線位置為y = 0， 兩條旁測線的位置分別為y = 1 km 和y = 2 km.3 條測線上接收點的位置固定不變， 坐標(biāo)分別為(—4， 0， 0)， (—4， 1， 0)和(—4， 2， 0) km.計算該模型顯式靈敏度與正演結(jié)果總共消耗CPU 時間是1 小時28 分， 占用的最大內(nèi)存為139.61 Gb.

在計算像素靈敏度之前， 首先給出3 條測線上水平電場 Ex的MVO 和PVO 曲線(圖9(a)和圖9(b)).可以看出， 由于3 條測線的位置差異， 這3 條不同測線產(chǎn)生的MVO 和PVO 曲線在異常體右側(cè)附近(0—4 km)產(chǎn)生明顯偏差， 在遠(yuǎn)離異常體區(qū)域， 不同測線產(chǎn)生的差異逐漸減小.圖9(c)和圖9(d)是3 條測線上水平電場 Ex的振幅與相位的塊狀體顯式靈敏度變化曲線， 在收發(fā)距位于2—10 km 的范圍內(nèi)， 主測線和旁測線上靈敏度均較為明顯， 但仔細(xì)比較可以看出， 相位靈敏度受測線位置影響更明顯， 這與不同測線上PVO 的差異加大完全符合.與塊狀體靈敏度類似， 像素靈敏度空間變化特征受收發(fā)距的影響也較大.

圖8 三維MCSEM 模型像素劃分與多測線位置示意圖 (a) y = 0 截面; (b) 俯視示意圖Fig.8.Pixel mesh and survey lines of three dimensional MCSEM model: (a) Cross-section of y = 0; (b) top view.

圖9 3 條不同測線時 E x 振幅和相位以及相對于塊狀異常體的靈敏度對比 (a)振幅; (b)相位; (c)不同測線上振幅靈敏度;(d)不同測線上相位靈敏度Fig.9.Comparison of magnitudes (MVO)， phases (PVO) and Fréchet sensitivity versus offset of E x obtained by three different survey lines: (a) MVO of E x ; (b) PVO of E x ; (c) comparison of E x amplitude sensitivity; (d) comparison of E x phase sensitivity.

根據(jù)圖9(c)和圖9(d)中塊狀體靈敏度曲線的變化特征， 在3 條測線上均選擇x = 4 km 作為發(fā)射源的位置， 利用每條測線上固定的發(fā)射與接收位置， 計算出 Ex振幅和相位的像素靈敏度， 圖10 是3 條測線上電場振幅和相位的像素靈敏度在垂直截面上的灰度圖， 可以看出， 發(fā)射與接收點周圍的像素靈敏度最強(qiáng)， 產(chǎn)生這一現(xiàn)象的主要原因是發(fā)射源周圍的電場較強(qiáng)， 其周圍各個像素單元上電導(dǎo)率的微小攝動會產(chǎn)生較強(qiáng)的散射電流， 引起空間電磁場的更大變化， 而接收點附近的電磁場雖然并不十分強(qiáng)， 但由于其周圍像素單元離接收點距離較近，像素單元中的散射電流對接收點上電磁場的影響也會較大.

在遠(yuǎn)離發(fā)射與接收點位置的其他像素單元上，在異常體邊界附近以及高阻異常體內(nèi)部的單元網(wǎng)格上的像素顯式靈敏度明顯比異常體外面低阻圍巖中的靈敏度的值大， 反映出海洋可控源電磁勘探對高阻異常體具有較強(qiáng)的探測能力.此外， 對比3 條不同測線上像素靈敏度的空間分布可以看出，因為y = 0 和y = 1 km 的垂直截面均穿過高阻異常體， 兩個不同垂直截面上空間靈敏度分布特征非常相似， 而y = 2 km 的垂直截面已經(jīng)在高阻異常體的外面， 該垂直截面上靈敏度空間分布與另兩個穿過異常體的截面上靈敏度存在著非常大的差異.

為進(jìn)一步研究水平截面上空間靈敏度的分布情況， 圖11 給出了z = 0.3， 0.6 和1 km 三個不同深度的水平截面上像素靈敏度的空間分布， 這里的發(fā)射和接收點均位于主測線且與圖10 中主測線上發(fā)射和接收點位置相同.在z = 0.3 和0.6 km 的水平截面上， 由于發(fā)射與接收點位置離水平截面比較近， 發(fā)射與接收點周圍的靈敏度非常大， 此外由于z = 0.6 km 的水平截面更接近高阻異常體上邊界(0.85 km)， 可以看到高阻異常與邊界附近單元上的靈敏度有明顯顯示.而在z = 1.0 km 的水平截面上， 由于發(fā)射與接收點位置離水平截面較遠(yuǎn)， 且水平截面穿過了高阻異常體， 靈敏度的空間分布充分地反映出高阻異常體的位置.可以看到， 像素靈敏度能夠更好地展示靈敏度的空間分布規(guī)律和特征原因， 在靈敏度空間分布中， 發(fā)射機(jī)和接收器對靈敏度的影響等同， 靈敏度最佳區(qū)域分布在收發(fā)器之間， 并隨著深度的增加逐漸衰減， 異常體由于其高阻特性， 在異常體內(nèi)部區(qū)域， 靈敏度有一定程度的降低， 而在異常體邊界區(qū)域， 由于感應(yīng)電流的影響， 靈敏度有明顯的提高， 同樣地， 空間靈敏度分布也顯示出相位包含的信息稍多于振幅.

圖10 發(fā)射和接收天線固定在不同截面情況下， 主測線和兩條旁測線在垂直截面上 E x 振幅和相位的像素靈敏度空間分布( r S =(4000， 0， -50) m 和 r R =(-4000， 0， -50) m) (a) E x 振幅靈敏度; (b) E x 相位靈 敏度Fig.10.The xz cross-sections along inline and two sidelines of E x pixel sensitivity distribution for a given source-receiver combination( r S =(4000， 0， -50) m and r R =(-4000， 0， -50) m): (a) Sensitivity of E x amplitude; (b) sensitivity of E x phase.

圖11 發(fā)射和接收天線固定在不同截面情況下， 3 個不同深度的水平截面上 E x 振幅和相位的像素靈敏度空間分布(rS =(4000， 0， -50) m 和 r R =(-4000， 0， -50) m) (a) E x 振幅靈敏 度; (b) E x 相位靈敏度Fig.11.The xy cross-sections at different depth of E x pixel sensitivity distribution for a given source-receiver combination( r S =(4000， 0， -50) m and r R =(-4000， 0， -50) m): (a) Sensitivity of E x amplitude; (b) sensitivity of E x phase.

4 結(jié) 論

本文基于電場耦合勢Helmholtz 方程的三維有限體積法與直接法求解技術(shù)， 并結(jié)合攝動理論研究建立了一套三維海洋可控源電磁響應(yīng)顯式靈敏度矩陣的快速算法， 能夠針對像素和分塊兩種不同的電導(dǎo)率模型有效地確定電導(dǎo)率微小攝動與電磁響應(yīng)微小變化間的線性關(guān)系.

在正演過程中， 利用直接求解器得到離散矩陣的逆并結(jié)合接收器位置事先計算出插值算子和投影算子， 由于插值算子和投影算子與發(fā)射源位置無關(guān)， 而僅僅與離散矩陣的逆矩陣和接收器位置有關(guān)， 通過投影算子與發(fā)射源離散向量的乘積就可以快速確定每個接收點上的電磁響應(yīng)， 有效提高了多發(fā)射源電磁場數(shù)值模擬的效率.

不論是在像素電導(dǎo)率還是分塊電導(dǎo)率模型中，利用Yee 氏交錯網(wǎng)格可以將電導(dǎo)率攝動量表示為剖分網(wǎng)格上的分塊常數(shù)函數(shù)， 每個剖分網(wǎng)格上電場強(qiáng)度與電導(dǎo)率相對攝動量的乘積可以作為散射電流元， 在一階Born 近似情況下， 散射電磁場實質(zhì)上可以看作各個剖分網(wǎng)格上所有散射電流元的疊加， 因此， 整個散射電磁場被轉(zhuǎn)換為多發(fā)射源電磁場的疊加， 利用投影算子和每個散射電流元離散向量的乘積可以快速獲得電磁場微小變化與電導(dǎo)率相對攝動間的線性關(guān)系， 得到顯式靈敏度計算結(jié)果， 從而大大提高了散射電磁場的計算效率.

數(shù)值結(jié)果顯示， 塊狀體靈敏度能夠更好地評估接收器發(fā)射源的偏移距與異常體探測能力的變化規(guī)律， 對于單個異常體模型， 在0.25 Hz 工作頻率下， 電場和磁場有效靈敏度的最大偏移距可以達(dá)到10 km 左右， 最小偏離距與接收器到異常體邊界的距離有關(guān)， 變化范圍為2—6 km.對于多個異常體模型， 利用塊狀體靈敏度可以快速分析考察不同收發(fā)距情況下電磁場響應(yīng)特征、確定最佳的接收器發(fā)射機(jī)組合.此外， 多個數(shù)值模型的計算結(jié)果均顯示， 在有效的偏移距范圍內(nèi)， 磁場的振幅相位靈敏度高于電場， 相位包含的靈敏度信息稍多于振幅.像素靈敏度能夠展示靈敏度空間分布規(guī)律和特征， 從靈敏度空間分布中確定出最佳探測區(qū)域.