耦合勢有限體積法模擬海洋移動(dòng)導(dǎo)線源三維頻率域電磁響應(yīng)

周建美 李 貅 戚志鵬(長安大學(xué)地質(zhì)工程與測繪學(xué)院,陜西西安 710054)

1 引言

海洋可控源電磁法(MCSEM)是一種重要的海上油氣勘探方法[1,2]。典型的MCSEM勘探采用水平導(dǎo)線源(長度一般為幾百米)隨船拖拽移動(dòng)發(fā)射低頻信號(hào)，在海底布設(shè)接收機(jī)陣列測量水平電場和三分量磁場響應(yīng)[3]。在海洋電磁勘探設(shè)計(jì)及資料處理和解釋過程中,均需要進(jìn)行大量的數(shù)值模擬。三維正演算法研究一直都是電磁法的研究熱點(diǎn)。MCSEM的三維正演相比其他電磁方法，有其特殊性：由于海水的濾波作用，MCSEM的工作頻率非常低(0.1～10Hz)，導(dǎo)致控制方程呈病態(tài)，出現(xiàn)所謂的低感應(yīng)數(shù)問題[4]； MCSEM的數(shù)值模擬需要同時(shí)計(jì)算近區(qū)和遠(yuǎn)區(qū)的電磁場響應(yīng)(如發(fā)射源頻率為0.25Hz時(shí)一般需要求解發(fā)射源附近±15km范圍內(nèi)的電磁場分布)，由于在大區(qū)域范圍內(nèi)電磁場的近場和遠(yuǎn)場的振幅相差很大(10個(gè)數(shù)量級(jí)以上),同時(shí)由于空氣層的存在使得模型電導(dǎo)率差很大,導(dǎo)致離散化方程的條件數(shù)非常大[5]； MCSEM是一種典型的移動(dòng)源電磁方法，需要對(duì)多次移動(dòng)源進(jìn)行三維模擬[6]。

目前主要采用有限差分法(FDM)[7-10]、有限體積法(FVM)[11]和有限元法(FEM)[12]模擬復(fù)雜三維地形的MCSEM響應(yīng)。這三種方法的共同點(diǎn)是直接對(duì)Maxwell方程進(jìn)行離散處理，得到一個(gè)大型的稀疏代數(shù)方程，借助于迭代法[13]或直接法[14]得到方程的數(shù)值解。直接法由于其對(duì)于矩陣的條件數(shù)不敏感,矩陣求解穩(wěn)定,能有效處理移動(dòng)源的三維正演問題[15]，廣泛應(yīng)用于MCSEM電磁響應(yīng)的數(shù)值模擬[12,14]，但直接法的內(nèi)存需求和計(jì)算時(shí)間隨著未知數(shù)個(gè)數(shù)呈指數(shù)增長[16]，極大地制約了直接法用于復(fù)雜模型三維正演計(jì)算的效率。低感應(yīng)數(shù)問題和大條件數(shù)的系數(shù)矩陣使常規(guī)的迭代算法很難有效求解MCSEM的離散線性方程[5]。利用電場的矢勢和標(biāo)勢分解可有效克服低感應(yīng)數(shù)問題[17]，結(jié)合有效的預(yù)處理能夠?qū)崿F(xiàn)MCSEM離散的大條件數(shù)線性方程的求解[5,13,18]。

在MCSEM三維正演計(jì)算中，一般將導(dǎo)線源近似為電偶極源。對(duì)于電偶極源采用直接離散源項(xiàng)的方法，但電偶極源的奇異性將導(dǎo)致源附近區(qū)域存在較大的離散誤差。目前主要采用一次場/二次場方法利用等效源[13,14]或者差異場技術(shù)[15]克服源的奇異性問題。一次場/二次場方法能夠有效模擬復(fù)雜場源形態(tài)的MCSEM三維響應(yīng)[14]，但需要額外計(jì)算所有網(wǎng)格節(jié)點(diǎn)上的一次場分布；差異場技術(shù)雖然能夠消除源附近的離散誤差，但需要兩次三維正演計(jì)算。對(duì)于實(shí)際MCSEM勘探中采用導(dǎo)線源的源項(xiàng)直接離散技術(shù)，這不同于電偶極源，在保證離散網(wǎng)格較導(dǎo)線源長度足夠小時(shí)可得到源附近較為精確的電磁場分布。

MCSEM是一種典型的移動(dòng)源電磁方法。對(duì)于移動(dòng)源三維正演計(jì)算，常規(guī)迭代算法設(shè)置迭代初始解為零，總的計(jì)算時(shí)間隨著移動(dòng)源個(gè)數(shù)呈線性增長。設(shè)定移動(dòng)發(fā)射源位置為網(wǎng)格中心位置，發(fā)射源移動(dòng)，離散網(wǎng)格大小不變，但整體隨著發(fā)射源移動(dòng)。由于發(fā)射源移動(dòng)前后地下電性分布變化不大，因此發(fā)射源移動(dòng)前后的電磁場分布也非常接近。選擇發(fā)射源移動(dòng)前的電磁場分布作為發(fā)射源移動(dòng)后電磁場的初始解，將有效減少迭代次數(shù)從而減少總的計(jì)算時(shí)間。

本文從Maxwell方程出發(fā)，首先通過引入電場的矢勢標(biāo)勢分解，將Maxwell方程轉(zhuǎn)化為一個(gè)關(guān)于矢勢和標(biāo)勢的混合亥姆霍茲方程定解問題[17]；采用洛倫茲規(guī)范[16]，借助Yee氏交錯(cuò)網(wǎng)格和FVM對(duì)旋度和散度算子進(jìn)行離散，得到稀疏對(duì)稱系數(shù)矩陣的離散方程；采用直接加源技術(shù)實(shí)現(xiàn)導(dǎo)線源的離散；采用基于ILUT預(yù)處理的BICGSTAB迭代算法[17]實(shí)現(xiàn)離散方程的迭代求解。采用改進(jìn)的迭代算法實(shí)現(xiàn)移動(dòng)源的快速正演，即采用發(fā)射源移動(dòng)前的電磁場響應(yīng)作為發(fā)射源移動(dòng)后的電磁場分布的初始解，從而提高迭代收斂速度，減少計(jì)算時(shí)間。

2 耦合勢有限體積算法

2.1 控制方程

設(shè)電磁場隨時(shí)間變化關(guān)系為eiω t，MCSEM三維正演對(duì)應(yīng)的控制方程為

(1)

(2)

(3)

(4)

為了保證控制方程的完整性，對(duì)式(1)兩邊求散度，得到

(5)

將式(5)展開并利用洛倫茲規(guī)范進(jìn)行簡單的整理，得到

(6)

將式(4)與式(6)聯(lián)立，表示為矩陣形式

(7)

(8)

由式(8)可知，該控制方程的系數(shù)矩陣是稀疏對(duì)稱陣。采用x方向的線電流源，式(8)進(jìn)一步簡化為

D1Ax+Cxφ=-iωμ0Jx

(9a)

D1Ay+Cyφ=0

(9b)

D1Az+Czφ=0

(9c)

(9d)

邊界處采用簡單的截?cái)噙吔鐥l件[9]

(A,φ)|?Ω=0

(10)

式中Ω為求解區(qū)域。

2.2 有限體積離散

式(9a)對(duì)應(yīng)于x方向的稀疏矩陣方程的有限體積離散，引入x方向?qū)?yīng)的半整數(shù)控制體積單元

(11)

(12)

定義控制體積單元內(nèi)的場值是均一的，則有

(13)

(14)

式(9b)～式(9d)左邊項(xiàng)的有限體積離散形式與式(9a)類似。

2.3 源項(xiàng)直接離散

源項(xiàng)離散包括式(9a)和式(9d)右邊項(xiàng)的有限體積離散。

對(duì)于式(12)右邊項(xiàng)的離散，本文采用直接加源方式。為了簡單起見，設(shè)定網(wǎng)格與發(fā)射源位置重合，發(fā)射源的兩個(gè)端點(diǎn)分別為(xa,yc,zd)和(xb,yc,zd)，其中xa和xb分別對(duì)應(yīng)x方向相鄰的兩個(gè)網(wǎng)格節(jié)點(diǎn)坐標(biāo)，yc和zd分別為y方向和z方向某一網(wǎng)格節(jié)點(diǎn)坐標(biāo)，則沿x方向長度為xb-xa的線電流源可以定義為

Jx=Iδ(y-yc)δ(z-zd)×

[H(x-xb)-H(x-xa)]

(15)

(16)

對(duì)于式(9d)右邊源項(xiàng)的有限體積離散，對(duì)應(yīng)的控制體積單元為

(17)

對(duì)于x方向的發(fā)射源

[δ(x-xb)-δ(x-xa)]

(18)

則對(duì)應(yīng)的有限體積離散為

(19)

2.4 移動(dòng)源控制方程的快速求解

對(duì)離散結(jié)果進(jìn)行整理,得到關(guān)于矢勢A和標(biāo)勢φ的方程組

(20)

3 算法驗(yàn)證

3.1 精度驗(yàn)證

通過與一維模型的半解析計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證FVM的計(jì)算精度； 通過對(duì)比偶極源與線源的計(jì)算精度，驗(yàn)證源項(xiàng)直接離散方法的有效性。一維模型如圖1所示，空氣層電阻率設(shè)為1.0×106Ω·m，海水層電阻率為0.25Ω·m，海水層厚度為1000m，海底地層電阻率為1.0Ω·m。分別計(jì)算發(fā)射源為x方向的水平電偶極源(長度為1m)和線源(長度為400m)兩種情況，發(fā)射源極矩均為1A·m，發(fā)射源位于海底上方100m，發(fā)射頻率為0.25Hz。接收機(jī)位于海底，接收沿測線方向的電場Ex。剖分網(wǎng)格單元數(shù)為121×121×101，最小網(wǎng)格長度為100m，即線源為4個(gè)網(wǎng)格長度，偶極源為1個(gè)網(wǎng)格長度，總的網(wǎng)格剖分范圍如圖1所示。分別采用基于ILUT預(yù)處理的BICGSTAB算法和基于Jacobi預(yù)處理的QMR算法[14]求解該一維模型。圖2為采用電偶極源時(shí)兩種算法的收斂曲線對(duì)比，得到小于1.0×10-8的歸一化殘差，基于Jacobi預(yù)處理的QMR算法需要3316次迭代，計(jì)算時(shí)間為918s；而基于ILUT預(yù)處理的BICGSTAB算法只需要336次迭代，計(jì)算時(shí)間為416s， 可見后者計(jì)算速度更快，本文在后面的數(shù)值計(jì)算中都采用該算法。本文的FVM計(jì)算結(jié)果與開源軟件Dipole1D[19]的半解析計(jì)算結(jié)果對(duì)比如圖3所示,圖4為線源與偶極源響應(yīng)結(jié)果對(duì)比。

圖1 海洋一維模型

圖2 兩種迭代求解算法的收斂曲線對(duì)比

對(duì)比圖4a與圖4b可知，線源與偶極源的響應(yīng)只是在發(fā)射源附近區(qū)域(偏移距小于2500m)存在差異(圖3中的紅色曲線清晰地顯示了在發(fā)射源附近區(qū)域兩種發(fā)射源響應(yīng)的差別)，在偏移距大于2500m的區(qū)域，兩種場源的響應(yīng)基本相同，與文獻(xiàn)[20]的結(jié)論一致。對(duì)于MCSEM，異常一般都出現(xiàn)在大偏移距區(qū)域，因此MCSEM采用偶極源近似線源是合理的。由圖4c可知，在偏移距大于2500m的區(qū)域， FVM計(jì)算的線源與電偶極源的振幅的相對(duì)誤差基本一致，均小于3%；在發(fā)射源附近，F(xiàn)VM計(jì)算的電偶極源振幅的相對(duì)誤差較大，達(dá)到了23%，線源振幅的最大相對(duì)誤差小于7%。由圖4d可知，不論是電偶極源還是線源，在整個(gè)求解區(qū)域(包括源所在位置)的相位差均小于2°； 線源發(fā)射的情況下，源區(qū)域的相位差更小。綜合圖4c和圖4d可知：①對(duì)電偶極源采用直接加源技術(shù)，會(huì)導(dǎo)致源附近區(qū)域的振幅存在較大誤差；②對(duì)線源采用直接加源技術(shù)，并保證源附近區(qū)域網(wǎng)格長度相比發(fā)射源長度足夠小(本文中發(fā)射源長度為400m，源附近網(wǎng)格長度為100m)，能夠得到發(fā)射源附近區(qū)域較為精確的響應(yīng)。

圖3 海洋一維模型正演響應(yīng)(a)電偶極源的Ex振幅； (b)電偶極源的Ex相位； (c)線源的Ex振幅; (d)線源的Ex相位

圖4 FVM計(jì)算的電偶極源與線源計(jì)算結(jié)果對(duì)比(a)線源與偶極源Ex振幅比值； (b)線源與偶極源Ex相位比值; (c)電偶極源和線源相

3.2 單異常體模型

本文采用迭代法計(jì)算移動(dòng)源模型的三維正演，利用解的初值重置技術(shù)減少迭代次數(shù)從而減少計(jì)算時(shí)間。為了分析算法的有效性，設(shè)計(jì)如圖5所示的三維模型。在海底下方1000m處存在一個(gè)直徑4km、厚度100m、電阻率為100Ω·m的圓形高阻油氣藏。發(fā)射源為沿x方向、長度為400m的線電流源，置于海底上方100m，發(fā)射頻率為1.0Hz。發(fā)射源沿x方向移動(dòng)，設(shè)定初始發(fā)射源與油氣藏中心在x方向的距離為2500m，每隔500m移動(dòng)一次，共移動(dòng)11次。接收機(jī)位于海底，接收沿測線方向的電場Ex。采用FVM算法計(jì)算該移動(dòng)源模型，剖分網(wǎng)格單元數(shù)為141×121×112，最小網(wǎng)格長度為100m(圖6)。

圖5 單異常體海洋三維模型

圖6 移動(dòng)源海洋3D模擬網(wǎng)格剖分示意圖(a)發(fā)射源移動(dòng)前； (b)發(fā)射源移動(dòng)后

圖7 FVM計(jì)算的海洋移動(dòng)源三維模型正演響應(yīng)(a)Ex振幅； (b)Ex相位

圖8 移動(dòng)源正演計(jì)算迭代收斂曲線

圖9 移動(dòng)源正演響應(yīng)計(jì)算時(shí)間對(duì)比

3.3 雙異常體模型

雙異常體三維模型如圖10所示。在海底下方1000m深度處存在兩個(gè)直徑為2km、厚度為100m、電阻率為100Ω·m的高阻油氣藏。發(fā)射源為沿x方向、長度為400m的線電流源，位于海底上方100m，發(fā)射頻率為1.0Hz, 發(fā)射源沿x方向移動(dòng)。設(shè)初始發(fā)射源距離第一個(gè)目標(biāo)體邊界在x方向上的距離為500m，每隔1000m移動(dòng)一次，共移動(dòng)7次。接收機(jī)位于海底，接收沿測線方向的電場Ex。剖分網(wǎng)格單元數(shù)為201×121×112，最小網(wǎng)格長度為100m，網(wǎng)格以發(fā)射源為中心。

圖10 雙異常體模型

圖11 雙異常體模型正演響應(yīng)(a)Ex振幅； (b)Ex相位

圖12 雙異常體模型移動(dòng)源正演計(jì)算迭代收斂曲線

圖13 雙異常體模型移動(dòng)源正演響應(yīng)計(jì)算時(shí)間對(duì)比

3.4 非均勻覆蓋層模型

為了進(jìn)一步研究本文算法的適用性，設(shè)計(jì)如圖14所示的非均勻覆蓋層的三維模型。 在海底下方1000m深度處存在一個(gè)直徑為4km、厚度為100m、電阻率為100Ω·m的圓形高阻油氣藏。海水下方覆蓋層中有兩個(gè)電阻率分別為10Ω·m和5Ω·m的三維異常體。發(fā)射源為沿x方向、長度為400m的線電流源，發(fā)射源始終保持在海底上方100m處，發(fā)射頻率為1.0Hz，發(fā)射源沿x方向移動(dòng)。設(shè)定初始發(fā)射源距離高阻油氣藏異常體邊界在x方向上的距離為500m，每隔500m移動(dòng)一次，共移動(dòng)11次。接收機(jī)位于海底沿地形布設(shè)，接收沿測線方向的電場Ex。網(wǎng)格剖分單元為141×121×112，最小網(wǎng)格長度為100m，網(wǎng)格以發(fā)射源為中心。

圖14 海洋非均勻覆蓋層三維模型(a)xoz平面圖； (b)xoy平面圖

圖15 非均勻覆蓋層三維模型響應(yīng)曲線(a)Ex振幅; (b)Ex相位

圖16 非均勻?qū)痈采w層模型的正演計(jì)算迭代收斂曲線對(duì)比

圖中實(shí)線為移動(dòng)源11次正演常規(guī)算法(即迭代初始解X0=0)的迭代收斂曲線，虛線為本文改進(jìn)算法(即迭代初始解X0=X1)計(jì)算的后10次正演(不包括第1次正演)的迭代收斂曲線

圖17 非均勻覆蓋層三維模型正演響應(yīng)計(jì)算時(shí)間對(duì)比

4 結(jié)束語

本文采用耦合勢有限體積法建立了一套快速計(jì)算海洋移動(dòng)導(dǎo)線源電磁響應(yīng)的三維數(shù)值模擬技術(shù)：基于電場矢勢和標(biāo)勢分解將Maxwell方程轉(zhuǎn)換為關(guān)于矢勢與標(biāo)勢的亥姆霍茲方程；采用洛倫茲規(guī)范整理得到對(duì)稱形式的離散系數(shù)矩陣；通過基于ILUT預(yù)處理的BICGSTAB算法迭代求解離散線性方程；采用直接加源技術(shù)實(shí)現(xiàn)導(dǎo)線源的精確快速模擬；針對(duì)移動(dòng)源正演響應(yīng)計(jì)算，提出采用解的初值重置技術(shù)加快正演速度，即采用發(fā)射源移動(dòng)前的電磁場響應(yīng)作為發(fā)射源移動(dòng)后電磁場分布的初始解，提高迭代收斂速度，減少迭代次數(shù)。

數(shù)值實(shí)例計(jì)算結(jié)果表明，采用直接加源技術(shù)，通過將導(dǎo)線源表示為多個(gè)電偶極源的疊加，能夠顯著減少源附近區(qū)域的計(jì)算誤差。對(duì)于移動(dòng)源的正演計(jì)算，采用解的初值重置技術(shù)，選擇發(fā)射源移動(dòng)前的電磁場分布作為移動(dòng)后控制方程迭代計(jì)算的初始解，當(dāng)發(fā)射源移動(dòng)前后空間中電磁場分布變化較小時(shí)，能顯著地提高正演速度；當(dāng)發(fā)射源移動(dòng)前后空間中電磁場分布變化較大時(shí)，正演速度提高有限。

參考文獻(xiàn)

[1] Constable S.Ten years of marine CSEM for hydrocarbon exploration.Geophysics,2010，75(5):A67-A81.

[2] 何展翔,孫衛(wèi)斌,孔繁恕等.海洋電磁法.石油地球物理勘探,2006,41(4):451-457.

He Zhanxiang,Sun Weibin,Kong Fanshu et al.Marine electromagnetic approach.OGP,2006,41(4):451-457.

[3] Constable S and Srnka L.An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration.Geophysics,2007,72(2):WA3-WA12.

[4] Newman G A and Alumbaugh D L.Three-dimensio-nal induction logging problems,Part 2:A finite-diffe-rence solution.Geophysics,2002,67(2):484-491.

[5] Puzyrev V.A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling.Geophysical Journal International,2003,193(2):678-693.

[6] Plessix R E,Darnet M and Mulder W A.An approach for 3D multisource,multi frequency CSEM modeling.Geophysics,2007,72(5):SM177-SM184.

[7] 殷長春,賁放,劉云鶴等.三維任意各向異性介質(zhì)中海洋可控源電磁法正演研究.地球物理學(xué)報(bào),2014,57(12):4110-4122.

Yin Changchun,Ben Fang,Liu Yunhe et al.MCSEM 3D modeling for arbitrarily anisotropic media.Chinese Journal of Geophysics,2014,57(12):4110-4122.

[8] 付長民,底青云,王妙月.海洋可控源電磁法三維數(shù)值模擬.石油地球物理勘探,2009,44(3):358-363.

Fu Changmin,Di Qingyun,Wang Miaoyue.3D nume-ric simulation of marine controlled source electromagnetics (MCSEM).OGP,2009,44(3):358-363.

[9] 嚴(yán)波,李予國,韓波等.任意方位電偶源的MCSEM電磁場三維正演.石油地球物理勘探,2017,52(4):859-868.

Yan Bo,Li Yuguo,Han Bo et al.3D marine controlled-source electromagnetic forward modeling with arbitrarily orientated dipole source.OGP,2017,52(4):859-868.

[10] 朱成,李桐林,楊海斌等.帶地形頻率域可控源電磁法三維反演研究.石油地球物理勘探,2016,51(5):1031-1039.

Zhu Cheng,Li Tonglin,Yang Haibin et al.3D controlled source electromagnetic inversion with topography in the frequency domain.OGP,2016,51(5):1031-1039.

[11] 楊波,徐義賢,何展翔等.考慮海底地形的三維頻率域可控源電磁響應(yīng)有限體積法模擬.地球物理學(xué)報(bào),2012,55(4):1390-1399.

Yang Bo,Xu Yixian,He Zhanxiang et al.3D frequency-domain modeling of marine controlled source electromagnetic responses with topography using finite volume method.Chinese Journal of Geophysics,2012,55(4):1390-1399.

[12] 楊軍,劉穎,吳小平.海洋可控源電磁三維非結(jié)構(gòu)矢量有限元數(shù)值模擬.地球物理學(xué)報(bào),2015,58(8):2827-2838.

Yang Jun,Liu Ying,Wu Xiaoping.3D simulation of marine CSEM using vector finite element method on unstructured grids.Chinese Journal of Geophysics,2015,58(8):2827-2838.

[13] 蔡紅柱,熊彬,Michael Zhdanov.電導(dǎo)率各向異性的海洋電磁三維有限單元法正演.地球物理學(xué)報(bào),2015,58(8):2839-2850.

Cai Hongzhu,Xiong Bin,Zhdanov M.Three-dimensional marine controlled-source electromagnetic mode-lling in anisotropic medium using finite element method.Chinese Journal of Geophysics,2015,58(8):2839-2850.

[14] 韓波,胡祥云,黃一凡等.基于并行化直接解法的頻率域可控源電磁三維正演.地球物理學(xué)報(bào),2015,58(8):2812-2826.

Han Bo,Hu Xiangyun,Huang Yifan et al.3-D frequency-domain CSEM modeling using a parallel direct solver.Chinese Journal of Geophysics,2015,58(8):2812-2826.

[15] 周建美,張燁,汪宏年等.耦合勢有限體積法高效模擬各向異性地層中海洋可控源的三維電磁響應(yīng).物理學(xué)報(bào),2014,63(15):159101-1.

Zhou Jianmei,Zhang Ye,Wang Hongnian et al.Efficient simulation of three-dimensional marine controlled-source electromagnetic response in anisotropic formation by means of coupled potential finite volume method.Acta Physica Sinica,2014,63(15):159101-1.

[16] Streich R.3D finite-difference frequency-domain mode-ling of controlled-source electromagnetic data:Direct solution and optimization for high accuracy.Geophy-sics,2009,74(5):F95.

[17] Haber E,Ascher U,Aruliah D et al.Fast simulation of 3D electromagnetic problems using potentials.Journal of Computational Physics,2000,163(1):150-171.

[18] Weiss C J. Project APhiD:a Lorenz-gauged A-decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth.Computers & Geosciences,2013,58(1):40-52.

[19] Key K.1D inversion of multi component,multi frequency marine CSEM data:Methodology and synthetic studies for resolving thin resistive layers.Geophysics,2009，74(2):F9-F20.

[20] Streich R and Becken M.Electromagnetic fields ge-nerated by finite-length wire sources:comparison with point dipole solutions.Geophysical Prospecting,2011，59(2):361-374.