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

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

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

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

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

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

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

0 引言

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

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

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

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

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

1 方法理論

1.1 異常場方程

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

E=EP+ES,H=HP+HS

(1)

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

▽×ES=iωμ0HS

▽×HS-σES=σSEP

(2)

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

通過傅立葉變換，場可以轉(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個三角單元，單元編號記為e=1、2、…、ne。方程(13)和方程(14)的積分則分解為各個單元積分之和：

(15)

(16)

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

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

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

(17)

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

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

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

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

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

4 數(shù)值計算

4.1 算法正確性檢驗—一維層狀模型的計算

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

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

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

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

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

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

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

模擬一個二維海洋地電模型。如圖5所示，空氣電阻率設(shè)為1012Ω·m，空氣層下方有厚度為1 000 m的海水層，電阻率為0.3 Ω·m，海底埋深1 000 m處有一個電阻率為100 Ω·m的油氣層，其長為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個接收站，分布范圍在y軸上從-6 000 m到6 000 m,間隔為100 m。

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

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

圖6 在不同發(fā)射頻率、不同埋藏深度和不同海底儲層大小下Ey振幅相對誤差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)不同油氣層的水平長度

5 結(jié)論與討論

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

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

發(fā)射頻率對海洋可控源電磁響應(yīng)影響較大，針對海洋可控源電磁法的工作頻率范圍0.001 Hz ～10 Hz，設(shè)計0.25 Hz、0.1 Hz、1 Hz三個發(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é)報,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

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

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

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

1001-1749(2017)05-0579-08