海洋可控源電磁法2.5維有限元正演模擬研究

韋艷剛， 郭榮文， 劉 穎， 柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083;2.中國(guó)海洋大學(xué) 海洋地球科學(xué)學(xué)院，青島 266003)

海洋可控源電磁法2.5維有限元正演模擬研究

韋艷剛1， 郭榮文1， 劉 穎2， 柳建新1

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083;2.中國(guó)海洋大學(xué) 海洋地球科學(xué)學(xué)院，青島 266003)

為了消除源點(diǎn)的奇異性問(wèn)題，在源點(diǎn)利用delta函數(shù)離散。設(shè)計(jì)一維的層狀模型，驗(yàn)證海洋可控源 2.5D 有限元數(shù)值模擬算法的正確性和精度；考慮到在海中洋流、海底地形起伏變化大以及海底介質(zhì)電性結(jié)構(gòu)復(fù)的影響，海洋可控源電磁的發(fā)射源會(huì)發(fā)生傾斜、旋轉(zhuǎn)、位置偏離測(cè)線，因此定義了任意取向的發(fā)射源。實(shí)際應(yīng)用證明，發(fā)射源姿態(tài)變化以及不同發(fā)射頻率和油氣層參數(shù),對(duì)海洋可控源電磁場(chǎng)響應(yīng)都有很大地影響。

海洋可控源電磁法; 2.5維; 有限元; 發(fā)射源

0 引言

海洋可控源電磁法(controlled-source electromagnetic methods,CSEM)，是勘探海底油氣資源和地質(zhì)構(gòu)造研究的一種海洋地球物理方法。在上個(gè)世紀(jì)80-90年代，海洋可控源電磁法在海洋地殼和上地幔等領(lǐng)域中得到廣泛地應(yīng)用[1-2]。在近十多年里，由于海洋可控源電磁法對(duì)高電阻率儲(chǔ)層具有很強(qiáng)的探測(cè)能力，其在海洋油氣資源和海底天然氣探測(cè)中也發(fā)揮著重要的作用[3-4]。目前已有不少關(guān)于海洋可控源電磁法的文章被發(fā)表，在海洋可控源電磁法正演模擬方法研究的初期，主要以一維模型正演為主，但是1D正演方法不適合探測(cè)海底較為復(fù)雜的結(jié)構(gòu)。對(duì)于三維模型,雖然Newman等[5]、Weiss等[6]使用交錯(cuò)網(wǎng)格有限差分法可以模擬三維不均勻結(jié)構(gòu)模型，但是這種方法采用結(jié)構(gòu)化網(wǎng)格，計(jì)算的精度和速度并不能滿足當(dāng)代的工作需求。

實(shí)際情況下，將三維問(wèn)題轉(zhuǎn)化為二維問(wèn)題，大大減少計(jì)算機(jī)內(nèi)存，提高計(jì)算速度，并且Ali Moradi Tenhrani等[5]分析了不同海底儲(chǔ)層大小、不同埋藏深度和不同發(fā)射頻率對(duì)2.5D和3D海洋電磁場(chǎng)的影響。結(jié)果分析表明，2.5D海洋可控源電磁響應(yīng)與3D電磁響應(yīng)的相對(duì)誤差非常小，應(yīng)用2.5D海洋可控源電磁法代替3D是具有合理性的，推動(dòng)了2.5D海洋可控源電磁法的快速發(fā)展。

通過(guò)有限元法離散得到一個(gè)大型、稀疏、復(fù)系數(shù)對(duì)稱線性方程組，求解大型稀疏線性方程組的方法有很多種，而很多學(xué)者經(jīng)過(guò)廣泛深入地研究[6-12]，Krylov子空間方法是最受歡迎的方法。Roland W.Freund等[13]提出了Quasi-Minnimal Residual(QMR)方法。這種方法基于Lanczos雙共軛過(guò)程，核心思想是極小化處在Krylov子空間的解的殘量范數(shù)的短遞歸方法。Weiss[14]首次將QMR法運(yùn)用到3D電磁感應(yīng)線性方程組求解,QMR法在迭代求解中數(shù)值能夠平滑穩(wěn)定收斂并且逼近解的準(zhǔn)確殘量范數(shù)，不需要顯示系數(shù)矩陣的存儲(chǔ)，因此減小了內(nèi)存空間，大大提高線性方程組的求解速度。

在海洋可控源電磁法探測(cè)過(guò)程中，由于海底水流動(dòng)性和儀器本身等問(wèn)題，發(fā)射源容易發(fā)生旋轉(zhuǎn)和傾斜，場(chǎng)源形態(tài)的變化均會(huì)對(duì)海洋CSEM數(shù)據(jù)產(chǎn)生影響。Streich等[15]認(rèn)為發(fā)射源方位和位置偏差等對(duì)電磁場(chǎng)響應(yīng)產(chǎn)生影響，因此我們引入了水平旋轉(zhuǎn)角θ和傾角Φ來(lái)定義任意取向的發(fā)射源，分別討論發(fā)射源發(fā)生旋轉(zhuǎn)以及傾斜與無(wú)旋轉(zhuǎn)傾斜時(shí)對(duì)探測(cè)油氣層的影響。

筆者基于麥克斯韋方程組出發(fā)，推導(dǎo)出2.5D海洋可控源電磁法波數(shù)域電磁場(chǎng)的耦合微分方程；通過(guò)有限元法進(jìn)行數(shù)值離散得到一個(gè)大型線性方程組，使用QMR法迭代求解。先設(shè)計(jì)一維的層狀模型，驗(yàn)證海洋可控源 2.5D 有限元數(shù)值模擬算法的正確性和精度，并考慮電偶極子發(fā)生偏差與無(wú)偏差時(shí)誤差。最后采用最優(yōu)觀測(cè)方式探索不同海底儲(chǔ)層大小、不同埋藏深度和不同發(fā)射頻率，對(duì)海底探測(cè)油氣層的影響。

1 方法理論

1.1 異常場(chǎng)方程

為了消除源點(diǎn)的奇異性問(wèn)題，在有限元中應(yīng)用疊加原理，將電磁場(chǎng)分解由電偶源在一維水平層狀(電導(dǎo)率σp)中產(chǎn)生的一次場(chǎng)(EP和HP)和由二維異常體(電導(dǎo)率為σs=σ-σp)產(chǎn)生的二次場(chǎng)(ES和HS)：

E=EP+ES,H=HP+HS

(1)

設(shè)計(jì)一個(gè)二維海洋地電模型，其走向方向假定為沿x軸,垂直向下的方向?yàn)閦軸,與x軸相垂直的是y軸，水平電偶源位于海底上方海水中。假定時(shí)間因子為 e-iωt,似穩(wěn)態(tài)情形下，二次場(chǎng)滿足偏微分方程：

▽×ES=iωμ0HS

▽×HS-σES=σSEP

(2)

式中：E為電場(chǎng)強(qiáng)度(V/m);H為磁場(chǎng)強(qiáng)度(A/m);μ0=4π×10-7(H/m)為介質(zhì)磁導(dǎo)率；ω為角頻率；σ為介質(zhì)電導(dǎo)率張量(S/m)

通過(guò)傅立葉變換，場(chǎng)可以轉(zhuǎn)換為：

(3)

(4)

(5)

1.2 有限元方程

(6)

(7)

式中:n表示邊界上的外法線單位向量；ny和nz分別表示沿著y軸和z軸的單位向量。

(8)

利用矢量公式

a▽·b=▽·(ab)-▽a·b

(9)

和散度公式

(10)

以及格林公式

(11)

(12)

方程(8)則變成為

(13)

(14)

將求解區(qū)域分解成ne個(gè)三角單元，單元編號(hào)記為e=1、2、…、ne。方程(13)和方程(14)的積分則分解為各個(gè)單元積分之和：

(15)

(16)

采用規(guī)則的矩形單元對(duì)模擬區(qū)域進(jìn)行有限元離散，在單元積分內(nèi)采用雙線性插值,得到一個(gè)大型、稀疏、復(fù)系數(shù)對(duì)稱線性方程組Ku=p，對(duì)其優(yōu)化處理后采用QMR(Quasi-Residual-Method，擬最小殘差法)迭代求解,分別得到總場(chǎng)和二次場(chǎng)各個(gè)分量。

2 源點(diǎn)奇異性處理

海洋可控源電磁法與大地電磁法最大的區(qū)別，在于引入人工源，源的處理一直是一個(gè)難點(diǎn)。在不需要源點(diǎn)附近增加網(wǎng)格數(shù)量的前提下，避免源點(diǎn)在遠(yuǎn)處產(chǎn)生奇異性，而獲得高精度的有限元解，需要采用 HERRMANN[16-17]提出的偽δ函數(shù)來(lái)等效場(chǎng)源的作用，其中偽δ的函數(shù)表達(dá)式為式(17)。

(17)

其中:x為離散點(diǎn)坐標(biāo)；x0為中心點(diǎn)坐標(biāo)；參數(shù)τ是控制源分布寬度和幅值的參數(shù)。顯然，使用偽δ函數(shù)描述源項(xiàng)，可使總體方程右端項(xiàng)不為“0”的節(jié)點(diǎn)增多，改善了總體方程組的不穩(wěn)定性，可以精準(zhǔn)地模擬場(chǎng)源的分布，提高計(jì)算精度。

3 任意取向發(fā)射源的定義

海洋可控源電磁法采用船載可移動(dòng)水平電偶極子源在海水上方一定高度上勻速前進(jìn)，由于海底洋流和電纜本身重量等問(wèn)題,不可能達(dá)到精確的水平沿測(cè)線方向,會(huì)存在一定的偏差,因此發(fā)射源取向往往是任意的。

如圖1所示,假設(shè)定義發(fā)射源水平旋轉(zhuǎn)角θ為電偶極子源沿順時(shí)針?lè)较虺痽軸旋轉(zhuǎn)的角度,發(fā)射源傾角Φ為電偶極子源朝z軸正向偏離xoy面的角度。發(fā)射源沿x方向布置時(shí)(HED-x),此時(shí)為赤道裝置(broadside 觀測(cè)),激發(fā)的電場(chǎng)主要為x分量，y分量和z分量十分微弱，當(dāng)發(fā)射源沿y方向布置時(shí)(HED-y),此時(shí)為軸向裝置(inline 觀測(cè)),激發(fā)的電場(chǎng)主要為y分量，x分量和z分量極微弱，發(fā)射源沿z方向布置時(shí),傾角φ=90°,為垂直偶極源(VED)。

圖1 電偶極子方位示意圖Fig.1 The orientation of an electric dipole

4 數(shù)值計(jì)算

4.1 算法正確性檢驗(yàn)—一維層狀模型的計(jì)算

為了驗(yàn)證算法和程序的正確性，計(jì)算一個(gè)高阻儲(chǔ)層的一維經(jīng)典模型(Constable[21]),如圖2所示,海水電阻率為0.3 Ω·m,海底深度為 1 000 m ;海底下方為 1 000 m 有一個(gè)厚度100 m電阻率為100 Ω·m的高阻油氣薄層；油氣層上方與下方均為巖層，其電阻率為1 Ω·m;電偶源布置于海底上方50 m處,長(zhǎng)度為1 m,發(fā)射電流為1 A ,頻率為0.25 Hz ;測(cè)線沿著y方向布置，測(cè)線長(zhǎng)度為12 km，共有121個(gè)接收點(diǎn),間隔為100 m均勻布置于海底。

圖2 1D層狀模型Fig.2 1D layered model

計(jì)算結(jié)果如圖3所示，將本文計(jì)算的電磁場(chǎng)分量的振幅、相位與水平層狀模型的一次場(chǎng)解析解做對(duì)比，其中圖3(a)、圖3(b)、圖3(e)、圖3(f)分別為軸向裝置(發(fā)射源沿y方向布置)所對(duì)應(yīng)的Hx、Ey、Ez分量振幅相位及其誤差，圖3(c)、圖3(d)、圖3(g)、圖3(h)分別為赤道裝置(發(fā)射源沿x方向布置)所對(duì)應(yīng)的Ex、Hy、Hz分量振幅相位及其誤差。由于在海洋CSEM勘探中,發(fā)射天線通常由較長(zhǎng)的導(dǎo)線構(gòu)成，此時(shí)若仍將發(fā)射源看作電偶極子源，將導(dǎo)致在源點(diǎn)附近的電磁場(chǎng)計(jì)算結(jié)果存在明顯誤差。從圖3(a)～圖3(h)可以發(fā)現(xiàn)，與1D解析解算法相比，在收發(fā)距較小的情況下，電磁場(chǎng)分量的振幅及相位誤差都比較大，但除在小收發(fā)距處外，各個(gè)分量振幅的相對(duì)誤差都小于4%，相位絕對(duì)誤差都小于5°，有限元解與解析解吻合得比較好，說(shuō)明本文的算法和程序是正確的和有效的。通過(guò)計(jì)算和對(duì)比inline 觀測(cè)與broadside 觀測(cè)的電磁場(chǎng)響應(yīng)，我們得到一個(gè)結(jié)論，inline 觀測(cè)得到振幅相對(duì)誤差和相位絕對(duì)誤差比broadside 觀測(cè)小，并且在偏移距較遠(yuǎn)時(shí)，都趨向于0，inline 觀測(cè)為電磁場(chǎng)響應(yīng)較優(yōu)的觀測(cè)方式。

4.2 發(fā)射源位置偏離對(duì)海洋CSEM響應(yīng)的影響

由于發(fā)射源方向的任意性，為了研究發(fā)射源位置偏離對(duì)海洋CSEM響應(yīng)的影響，采用圖2 的一維層狀模型進(jìn)行模擬，分析不同水平旋轉(zhuǎn)角和傾角海洋CSEM響應(yīng)的影響。圖4左邊為傾角φ=0°，水平旋轉(zhuǎn)角θ=90°、87°、85°、83°、81°得到的電磁場(chǎng)振幅值，右邊以θ=90°，φ=0°為基準(zhǔn)，水平旋轉(zhuǎn)角偏差為3°、5°、7°、9°與水平旋轉(zhuǎn)角無(wú)偏差的誤差圖。當(dāng)發(fā)射源由y軸逐漸向x軸旋轉(zhuǎn)時(shí)，點(diǎn)偶極矩在y軸上的投影越來(lái)越小，使得y方向上的電場(chǎng)分量越來(lái)越小。從圖4右邊發(fā)現(xiàn)，隨著偏移距的增大，振幅的相對(duì)誤差越來(lái)越小，但是水平旋轉(zhuǎn)角偏離的角度越大，振幅的相對(duì)誤差就越大，當(dāng)水平旋轉(zhuǎn)角偏離9°時(shí)，在收發(fā)距小于2 500 m時(shí)，y方向的電場(chǎng)和z方向的磁場(chǎng)的振幅相對(duì)誤差都接近于10%，相比較而言，水平方向的磁場(chǎng)受角度的影響較小，因此在實(shí)際海洋CSEM勘探中，在小偏移距下水平旋轉(zhuǎn)角偏離測(cè)線大于9°時(shí)，有必要將發(fā)射源進(jìn)行調(diào)整或者計(jì)算時(shí)考慮誤差的存在。

圖4 1D模型不同水平旋轉(zhuǎn)角振幅及振幅相對(duì)誤差Fig.4 Different horizontal rotation angle amplitude and it’s relative errof of 1D model

4.3頻率及油氣層模型參數(shù)對(duì)2.5DMCSEM電磁響應(yīng)的影響

模擬一個(gè)二維海洋地電模型。如圖5所示，空氣電阻率設(shè)為1012Ω·m，空氣層下方有厚度為1 000 m的海水層，電阻率為0.3 Ω·m，海底埋深1 000 m處有一個(gè)電阻率為100 Ω·m的油氣層，其長(zhǎng)為4 000 m,高為100 m,在y軸上范圍從-2 000 m到2 000 m,在z軸上的范圍從1 000 m到1 100 m，圍巖電阻率為1 Ω·m，發(fā)射源位置為海底上方50 m，y= 0 m處,在海底共有 121個(gè)接收站，分布范圍在y軸上從-6 000 m到6 000 m,間隔為100 m。

圖5 2D油氣模型Fig.5 The 2D reservoir model

在圖5二維油氣模型基礎(chǔ)上，圖6為inline 觀測(cè)下y方向的電場(chǎng)在不同發(fā)射源頻率,不同油氣層的厚度,不同油氣層的埋藏深度,不同油氣層的水平長(zhǎng)度振幅的相對(duì)誤差，設(shè)計(jì)0.25 Hz、0.1 Hz、1 Hz三個(gè)發(fā)射頻率，在圖6(a)中發(fā)現(xiàn)，不同的發(fā)射頻率電場(chǎng)振幅的相對(duì)誤差不同，隨著發(fā)射頻率的增大，Ey振幅的相對(duì)誤差也增大，在偏移距較小尤其明顯，f=0.25 Hz所得到振幅的相對(duì)誤差比較小，因此選擇適當(dāng)?shù)念l率對(duì)海洋可控源電磁法的觀測(cè)有著重要的意義。在圖6(b)～圖6(d)中，在發(fā)射頻率為0.25 Hz下，不同油氣層的厚度，不同水平長(zhǎng)度以及不同油氣層埋藏深度的得到的y方向的電場(chǎng)的振幅相對(duì)誤差接近吻合，因此，不同油氣層的厚度，不同水平長(zhǎng)度以及不同油氣層埋藏深度，對(duì)2.5DMCSEM響應(yīng)的影響不大。

圖6 在不同發(fā)射頻率、不同埋藏深度和不同海底儲(chǔ)層大小下Ey振幅相對(duì)誤差Fig.6 Eymplitude relative error for different frequencies of the emitted field. different depths of the reservoir and different size of the reservoir(a)不同發(fā)射源頻率；(b)不同油氣層的厚度；(c)不同油氣層的埋藏深度;(d)不同油氣層的水平長(zhǎng)度

5 結(jié)論與討論

筆者詳細(xì)推導(dǎo)了有限元2.5D海洋CSEM電磁場(chǎng)響應(yīng)，采用疊加原理，并驗(yàn)證算法的正確性和有效性。在一維層狀模型探測(cè)中，inline 觀測(cè)較broadside 觀測(cè)反應(yīng)靈敏。

發(fā)射源位置的偏離對(duì)海洋可控源電磁響應(yīng)影響較大，尤其是在收發(fā)距較小的情況，水平方向的電場(chǎng)振幅隨著偏離角度的增大而增大，當(dāng)偏離角度達(dá)到9° 時(shí)，振幅相對(duì)誤差接近10%，在實(shí)際海洋CSEM勘探中，偏離角度較大時(shí)有必要考慮對(duì)海洋可控源電磁響應(yīng)影響。

發(fā)射頻率對(duì)海洋可控源電磁響應(yīng)影響較大，針對(duì)海洋可控源電磁法的工作頻率范圍0.001 Hz ～10 Hz，設(shè)計(jì)0.25 Hz、0.1 Hz、1 Hz三個(gè)發(fā)射頻率，在上述模型中，f=0.25 Hz是比較適當(dāng)?shù)陌l(fā)射頻率，在海洋可控源電磁響應(yīng)中選擇適當(dāng)?shù)念l率很重要。

[1] YOUNG P D,COX C S.Electromagnetic active source sounding near the East Pacific Rise[J].Geophysical Research Letters,1981(8):1043-1046.

[2] MACGREGOR L M, CONSTABLE S, SINHA M C. The RAMESSES experiment III: Controlled-source electromagnetic sounding of the Reykjanes Ridge at 57°45’57N [J]. Geophysical Journal International, 1998,135:773-789.

[3] EIDESMO,T.S.ELLINGSRUD,L.MACGREGOR,et al. Sea bed logging(SBL), a new method for remote and direct identification of hydrocarbon filled layers in deep water areas[J]. First Break,2002(20):144-152.

[4] CONSTABLE S, WEISS C J. Mapping thin resistors and hydrocarbons with marine em methods: Insights from 1D modeling[J]. Geophysical, 2006,(71): G43- G51.

[5] ALI MORADI TENHRANI, EVERT SLOD. Applicability of 1D and 2.5D marine controlled source electromagnetic modelling[J]. Geophysical Prospecting, 2013(61):602-613.

[6] LI Y, DAI S. Finite element modelling of marine controlled-source electromagnetic responses in two-dimensional dipping anisotropic conductivity structures[J]. Geophysical Journal International, 2011,185:622-636.

[7] SAAD Y,SCHULTZ MH.GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems[J].SIAMJ. Sci.Stat.Comput,1986(7):856-869.

[8] ROL,W.FREUND,NOEL M.NACHTIGAL,QMR: a quasiminimal residual method for non-Hermitian linear systems[J].Numer,Math,1991,60:315-339.

[9] FREUND R W,HOCHBRUCK M.On the use of two QMR algorithms for solving singular systems and applications in Markov chain modeling[J].Numerical Linear Algebra with Applications,2010,1(4):403-420.

[10] A.ESSAI,WEIGHTED FOM.Gmres for solving nonsymmetric linear systems[J].Numer, Algorithms,1998(18):277-292.

[11] Z.H.CAO, X.Y.YU. A note on weighted FOM and GMRES for solving nonsymmetric linear systems[J]. Appl.Math.Comput,2004,151:719-727.

[12] SAAD Y. Iterative Methods for Sparse Linear Systems[M]. second edition, SIAM, Philadelphia,2003.

[13] FREUND R W. Nachtigal N M. QMR: a quasi-minimal residual method for non-Hermitian linear systems[J]. Numerische Mathematik, 1991, 60(1):315-339.

[14] WEISS C J. A Matrix: free Approach to Solving the Fully 3D Electromagnetic Induction Problem[J]. Seg Technical Program Expanded Abstracts, 2001, 20(1):1451.

[15] STREICH R,MICHAEL B.Electromagnetic fileds generted by finite-length wire sources:comparison with ponit dipole solutions[J].Geophysical Prospecting,2001,59(2):361-374.

[16] HERRMANN R B. SH-wave generation by dislocation source-A numerical study[J]. Bulletin of the Seismological Society of America,1979, 69(1):1-15.

[17] 張繼鋒,湯井田,喻言,等.基于二次插值的線源可控源有限元數(shù)值模擬[J].吉林大學(xué)學(xué)報(bào),2009,39(5):929-935.

ZHANG J F, TANG J T, YU Y, et al. Finite element numerical simulation of controlled-source based on the quadratic interpolation line source[J].Jilin university journals,2009,39(5):929-935.(In Chinese)

2.5Dfiniteelementmodelingofmarinecontrolled-sourceelectromagneticfields

WEI Yangang1, GUO Rongwen1, LIU Ying2, LIU Jianxin1

(1.Central South University, School of Geosciences and Info-Physics，Changsha 410083，China;2.Ocean University of China, College of Marinca Geosciences, Qingdao 266003，China)

To eliminate singularity of source point, the delta function is used in discrete source point. We have designed a one-dimensional layered model to test the correctness and accuracy of marine controlled source's 2.5 D ( finite element method) algorithm. In the midst of the ocean currents, considering the influence of ocean current, fluctuation of seabed topography and electrical structure of seabed, the emission source of marine controlled electromagnetism will incline, rotate and deviate from the line. Therefore, we define the emission source with arbitrary orientation. The practical application shows that the variation of emission source attitude and the different emission frequency and reservoir parameters have great influence on the electromagnetic response of marine controlled source.

marine controlled-source electromagnetic method; 2.5 D; finite element; emission source

P 631.2

A

10.3969/j.issn.1001-1749.2017.05.01

2016-10-18 改回日期： 2016-11-07

國(guó)家自然科學(xué)基金(4167407965,4167408070)

韋艷剛(1991-)，女，碩士，主要從事海洋可控源電磁法研究，E-mail：1312877493@qq.com。

郭榮文(1980-)，男，副教授，主要從事電磁法勘探研究，E-mail：573377360@qq.com。

1001-1749(2017)05-0579-08