湯文武,李耀國,柳建新,劉春明

(1.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室,江西南昌330013;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,湖南長沙410083;3.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,湖南長沙410083;4.Center for Gravity,Electrical & Magnetic Studies,Department of Geophysics,Colorado School of Mines,Golden CO80401,USA)

基于二次電場(chǎng)的可控源電磁法三維矢量有限元正演模擬

湯文武1,2,3,李耀國4,柳建新2,3,劉春明2,3

(1.東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室,江西南昌330013;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,湖南長沙410083;3.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,湖南長沙410083;4.Center for Gravity,Electrical & Magnetic Studies,Department of Geophysics,Colorado School of Mines,Golden CO80401,USA)

可控源電磁法作為地震勘探的重要補(bǔ)充手段,可通過研究地下介質(zhì)電性的變化達(dá)到監(jiān)測(cè)油藏的目的。為了從理論上研究可控源電磁法在油藏監(jiān)測(cè)方面的應(yīng)用前景,在前人研究的基礎(chǔ)上,對(duì)三維可控源電磁法的正演模擬算法開展了進(jìn)一步的研究。為克服場(chǎng)源附近總場(chǎng)變化快、基于總場(chǎng)求解時(shí)難以精確模擬的問題,采用均勻半空間或水平層狀介質(zhì)模擬的一次場(chǎng)作為場(chǎng)源;針對(duì)傳統(tǒng)的節(jié)點(diǎn)型有限元模擬電場(chǎng)時(shí)存在散度條件不滿足的問題,采用矢量有限元對(duì)基于二次電場(chǎng)的偏微分方程進(jìn)行離散求解；多頻點(diǎn)正演時(shí)引入頻率適應(yīng)網(wǎng)格,可以在保持正演精度的同時(shí)提高計(jì)算速度。對(duì)比3層電導(dǎo)率模型的有限元數(shù)值解與解析解,驗(yàn)證了算法的正確性。通過一個(gè)簡單的油藏模型正演模擬,從理論上證明了可控源電磁法應(yīng)用于油藏監(jiān)測(cè)的可能性。

二次電場(chǎng);可控源電磁法;矢量有限元;正演模擬;油藏監(jiān)測(cè)

電磁法是以電性差異為基礎(chǔ)的一類勘探方法,作為地震勘探的一個(gè)重要補(bǔ)充手段已逐步應(yīng)用于油氣勘探及油藏監(jiān)測(cè)等領(lǐng)域[1-2]。電磁法數(shù)值模擬是電磁資料解釋的基礎(chǔ),受計(jì)算機(jī)資源及數(shù)值模擬方法的限制,目前實(shí)際生產(chǎn)中主要應(yīng)用一、二維正反演技術(shù)。近年來,隨著計(jì)算機(jī)可用內(nèi)存的擴(kuò)大及計(jì)算速度的不斷提高,三維電磁法數(shù)值模擬研究得到了極大的發(fā)展。如Newman等[3]發(fā)表了關(guān)于三維交錯(cuò)網(wǎng)格有限差分法的航空電磁數(shù)值模擬研究成果,通過計(jì)算二次場(chǎng)避免了源的奇異性。Druskin等[4]應(yīng)用頻譜Lanczos分解法進(jìn)行三維感應(yīng)電磁法的模擬計(jì)算,在避免偽解出現(xiàn)的同時(shí)有效地加快了線性方程組數(shù)值解的收斂。Haber等[5]、Badea等[6]、Puzyrev等[7]先后研究了基于電磁耦合勢(shì)的可控源三維電磁正演模擬,避免了直接基于電場(chǎng)求解低頻電磁響應(yīng)時(shí)可能出現(xiàn)的線性方程組數(shù)值解收斂慢的問題,但通過耦合勢(shì)求取電場(chǎng)需引入的數(shù)值導(dǎo)數(shù)運(yùn)算會(huì)損失部分精度。Streich[8]開展三維有限差分法可控源電磁模擬,利用MUMPS直接法求解器求解線性方程組,特別適用于多激發(fā)源的正演計(jì)算。Avdeev[9]對(duì)三維電磁數(shù)值模擬的正反演研究現(xiàn)狀進(jìn)行了較為系統(tǒng)的論述。國內(nèi)也有許多學(xué)者開展了三維電磁數(shù)值模擬研究,如張繼鋒等[10]基于電場(chǎng)的矢量波動(dòng)方程進(jìn)行了可控源電磁的有限元三維模擬研究,通過引入散度條件避免了偽解的出現(xiàn)。徐志鋒等[11]基于磁矢量勢(shì)及電標(biāo)量勢(shì)對(duì)三維可控源電磁進(jìn)行了模擬計(jì)算,通過引入罰項(xiàng)及穩(wěn)定化算法克服了數(shù)值不穩(wěn)定性。湯井田等[12]采用有限元-無限元耦合法對(duì)頻率域三維可控源電磁進(jìn)行了正演計(jì)算,較大程度地減少了網(wǎng)格計(jì)算區(qū)域,提高了正演速度。

可控源電磁法是繼大地電磁測(cè)深法后發(fā)展起來的一類人工源電磁勘探方法,因克服了天然場(chǎng)源信號(hào)的隨機(jī)性而得到了廣泛的應(yīng)用?？焖?、高精度的三維可控源電磁法正演模擬對(duì)于電磁法的靈敏度分析及反演至關(guān)重要。當(dāng)前三維可控源電磁法的正演模擬可直接基于電場(chǎng)求解或基于耦合勢(shì)求解。基于耦合勢(shì)的方式由于需要進(jìn)一步的數(shù)值求導(dǎo)(以得到電場(chǎng)),因而會(huì)損失部分精度?；陔妶?chǎng)的正演方式主要利用傳統(tǒng)的節(jié)點(diǎn)型有限元進(jìn)行求解,一般需要施加散度校正以獲得精確解。本文在前人研究工作的基礎(chǔ)上,對(duì)三維可控源電磁法的正演算法開展了進(jìn)一步的研究。首先從電磁場(chǎng)理論出發(fā),采用基于二次電場(chǎng)的偏微分方程進(jìn)行求解,應(yīng)用基于棱邊的矢量有限元法,采用頻率適應(yīng)網(wǎng)格進(jìn)行多頻點(diǎn)正演計(jì)算；然后設(shè)計(jì)了一個(gè)簡單的儲(chǔ)層模型,通過對(duì)電場(chǎng)分量的模擬計(jì)算,顯示了可控源電磁法在油氣勘探及油藏監(jiān)測(cè)中的應(yīng)用前景。

1 理論基礎(chǔ)

1.1 偏微分方程的建立

電磁類勘探方法均基于電磁場(chǎng)理論,在頻率域中，準(zhǔn)靜態(tài)條件下的電磁場(chǎng)滿足以下麥克斯韋方程組[13]:

式中:E指電場(chǎng)矢量，H指磁場(chǎng)矢量，Jc為外電流密度矢量；μ0為真空中的磁導(dǎo)率;σ為電導(dǎo)率。我們假設(shè)時(shí)諧因子為e-i ωt，對(duì)公式(1) 兩邊求取旋度并代入公式(2)，得到以下關(guān)于電場(chǎng)矢量的雙旋度方程：

(3)

對(duì)于電偶源、磁偶源等致密激發(fā)源,電場(chǎng)在源所處的位置是奇異的,并且在源的周圍變化極快。若直接計(jì)算總場(chǎng),需要在源的周圍加密網(wǎng)格,這將大大增加計(jì)算量。為了消除源的奇異性,可將總場(chǎng)分解為一次場(chǎng)和二次場(chǎng)(異常場(chǎng))。首先計(jì)算均勻半空間或水平層狀介質(zhì)模型的一次場(chǎng),再將一次場(chǎng)作為場(chǎng)源項(xiàng)來求解二次場(chǎng)，有:

(4)

式中:Ep為一次場(chǎng)，Es為二次場(chǎng)。與推導(dǎo)的總場(chǎng)雙旋度方程類似，對(duì)應(yīng)于一次場(chǎng)的雙旋度方程表示如下：

(5)

式中:σp對(duì)應(yīng)計(jì)算一次場(chǎng)時(shí)的電導(dǎo)率分布,一般采用均勻半空間或水平層狀介質(zhì)模型。將公式(3)減去公式(5)并進(jìn)行化簡后得到:

(6)

式中:σa為異常電導(dǎo)率分布，且σa=σ-σp。因此，右端項(xiàng)iωμ0σaEp即為計(jì)算二次場(chǎng)的場(chǎng)源項(xiàng)，可通過計(jì)算異常電導(dǎo)率區(qū)域的一次場(chǎng)得到。而均勻半空間及水平層狀介質(zhì)模型的電場(chǎng)可通過快速漢克爾變換[14]計(jì)算得到。相比總場(chǎng)計(jì)算方式，基于二次場(chǎng)的模擬計(jì)算方式的另一個(gè)優(yōu)勢(shì)是，通過一次場(chǎng)的計(jì)算來加載源項(xiàng)，統(tǒng)一了源項(xiàng)的處理，因而適用于不同激發(fā)源下的電磁場(chǎng)響應(yīng)的計(jì)算。

1.2 邊界條件

由于采用基于二次場(chǎng)的偏微分方程進(jìn)行正演計(jì)算,當(dāng)邊界距離目標(biāo)研究區(qū)域足夠遠(yuǎn)時(shí),二次場(chǎng)可近似為0。因此,可采用邊界值為0的狄利克雷邊界條件:

(7)

式中:n為邊界面上的法向單位矢量。值得注意的是,當(dāng)一次場(chǎng)是通過均勻半空間模型計(jì)算得到時(shí),該邊界條件相當(dāng)于基于總場(chǎng)方式時(shí)采用遠(yuǎn)區(qū)的均勻半空間值作為邊界值的人工截?cái)噙吔鐥l件[7-8]。當(dāng)采用更復(fù)雜的模型計(jì)算一次場(chǎng)時(shí),該邊界條件則變得更優(yōu)。

2 矢量有限元分析

為了精確表征電場(chǎng)法向分量在不同電導(dǎo)率界面上的不連續(xù)性,這里采用矢量有限元[15]對(duì)基于二次電場(chǎng)的偏微分方程進(jìn)行離散求解,采用伽遼金有限元方法推導(dǎo)其有限元方程。

2.1 有限元網(wǎng)格剖分

為了討論方便,采用立方體剖分單元對(duì)整個(gè)模擬計(jì)算區(qū)域進(jìn)行剖分。如圖1所示,中間的目標(biāo)研究區(qū)域采用較密的網(wǎng)格,由內(nèi)而外采用漸變稀疏的網(wǎng)格進(jìn)行剖分,注意網(wǎng)格大小遞增的比例不超過1.5。采用這種非均勻漸變稀疏網(wǎng)格,一方面能在目標(biāo)研究區(qū)域內(nèi)獲得精確的電磁場(chǎng)值,另一方面,在邊界條件得到滿足的情況下,網(wǎng)格數(shù)更少,從而可以減少正演計(jì)算量。對(duì)于包含多個(gè)頻率的正演計(jì)算,通常將邊界置于2～3倍最低頻率下的趨膚深度之外。

圖1 網(wǎng)格剖分

2.2 有限元方程推導(dǎo)

在基于棱邊的矢量有限元中,電場(chǎng)分量被置于各個(gè)剖分單元的棱邊上的中點(diǎn)處,與交錯(cuò)網(wǎng)格有限差分法的做法類似。如圖2所示,每個(gè)剖分單元內(nèi)電場(chǎng)分量的采樣數(shù)為12。假設(shè)對(duì)應(yīng)于各分量的插值函數(shù)(形函數(shù))為uj,則一次電場(chǎng)及二次電場(chǎng)可分別表示如下:

(8)

(9)

式中:Epj及Esj分別為各個(gè)棱邊中點(diǎn)的一次電場(chǎng)及二次電場(chǎng)值,N為剖分區(qū)域內(nèi)總的棱邊數(shù)。將公式(8)及公式(9)代入公式(6)中得到:

(10)

求取每個(gè)插值函數(shù)ui與公式(10)等號(hào)兩邊各項(xiàng)的內(nèi)積,并對(duì)其在整個(gè)剖分區(qū)域V上進(jìn)行積分,可得到如下方程:

圖2 電場(chǎng)分量在各個(gè)單元的布置

(11)

(12)

式中:S代表外邊界面,n為外邊界面上的外法向單位矢量。當(dāng)應(yīng)用邊界值為0的狄利克雷邊界條件時(shí),公式(12)中的面積分項(xiàng)為0,公式(11)可化為:

(13)

3 線性方程組的求解

對(duì)公式(13)離散并積分后,可得到一個(gè)大型的對(duì)稱復(fù)系數(shù)線性方程組。對(duì)于三維電磁法的數(shù)值模擬計(jì)算而言,該線性方程組的求解尤為重要。求解方法通常要兼顧內(nèi)存需求及計(jì)算速度，本文考慮三維電磁模擬程序在個(gè)人計(jì)算機(jī)上運(yùn)行的需要,采用內(nèi)存需求方面更經(jīng)濟(jì)的迭代法。對(duì)剛度矩陣采用只存儲(chǔ)非零元素的稀疏存儲(chǔ)方式,同時(shí),由于線性方程組條件數(shù)極大,采用不完全喬列斯基分解法對(duì)方程組進(jìn)行預(yù)條件處理[16]。綜合對(duì)比了各種常用的Krylov迭代法后,選取了最適合本文研究課題的QMR(Quasi-Minimal Residual)方法[17]對(duì)方程組進(jìn)行迭代求解。

4 頻率適應(yīng)網(wǎng)格及其應(yīng)用效果

進(jìn)行頻率域可控源電磁法正演模擬及反演成像時(shí),通常需要同時(shí)獲得多個(gè)頻點(diǎn)的電磁場(chǎng)數(shù)據(jù)。而電磁法勘探的探測(cè)深度與頻率有關(guān),為了保證計(jì)算精度,設(shè)計(jì)網(wǎng)格時(shí)要考慮特定的頻率大小，通常采用兩種方法來實(shí)現(xiàn):一是針對(duì)每個(gè)頻率設(shè)計(jì)相對(duì)應(yīng)的計(jì)算網(wǎng)格,這樣可使得網(wǎng)格最小,正演計(jì)算時(shí)間最短,但是需要花費(fèi)額外的時(shí)間來設(shè)計(jì)多套網(wǎng)格,而且不易應(yīng)用到多頻點(diǎn)數(shù)據(jù)的反演中;二是設(shè)計(jì)一個(gè)同時(shí)兼顧所有頻點(diǎn)下電磁場(chǎng)數(shù)據(jù)計(jì)算精度的網(wǎng)格,這將使得網(wǎng)格最大化,正演時(shí)間會(huì)大大增加。

考慮到文中采用邊界值為0的第一類邊界條件,在中心網(wǎng)格之外需要額外的非均勻網(wǎng)格擴(kuò)展到2～3倍最大趨膚深度(與最低頻率對(duì)應(yīng))以盡量消除邊界的影響。為了論述方便,我們將這些額外用來消除邊界影響的網(wǎng)格稱為附加網(wǎng)格,中間區(qū)域的網(wǎng)格稱為目標(biāo)網(wǎng)格。實(shí)際上,隨著頻率升高,趨膚深度減小,參與計(jì)算的附加網(wǎng)格數(shù)也可減少。因此,我們首先根據(jù)所有頻率設(shè)計(jì)好包含目標(biāo)網(wǎng)格及附加網(wǎng)格的固定網(wǎng)格,當(dāng)從低頻至高頻進(jìn)行正演計(jì)算時(shí),逐步減少附加網(wǎng)格的數(shù)目(保持目標(biāo)網(wǎng)格不變),從而逐步減少高頻時(shí)正演的計(jì)算量,最終達(dá)到減少總的正演計(jì)算時(shí)間的目的。同時(shí),由于目標(biāo)網(wǎng)格及各個(gè)頻點(diǎn)利用到的附加網(wǎng)格的網(wǎng)格間距都始終保持不變,因此可以直接將其應(yīng)用到反演過程中,減少網(wǎng)格設(shè)計(jì)的工作量。

圖3為一個(gè)3層水平電阻率模型。第1層的電阻率為100Ω·m,厚度為1075m;第2層的電阻率為10Ω·m,厚度為250m;第3層的電阻率為100Ω·m。單位電偶源置于坐標(biāo)原點(diǎn)O處,分別采用非頻率適應(yīng)網(wǎng)格(一套同時(shí)兼顧各個(gè)頻率的固定網(wǎng)格)及頻率適應(yīng)網(wǎng)格對(duì)頻率范圍為1～128Hz的8個(gè)頻點(diǎn)進(jìn)行正演計(jì)算。

圖4為兩種不同網(wǎng)格的正演計(jì)算精度對(duì)比,所計(jì)算的頻率為16Hz。其中圖4a為二次電場(chǎng)x分量,圖4b為二次電場(chǎng)y分量。由圖4可見,兩種不同網(wǎng)格正演計(jì)算得到的二次水平電場(chǎng)分量基本一致,頻率適應(yīng)網(wǎng)格與非頻率適應(yīng)網(wǎng)格的相對(duì)誤差小于1%,說明了頻率適應(yīng)網(wǎng)格的可行性。表1為不同頻率下兩種網(wǎng)格的正演計(jì)算時(shí)間對(duì)比情況。由表1 可知,在多頻點(diǎn)正演計(jì)算時(shí),采用頻率適應(yīng)網(wǎng)格可以使得高頻點(diǎn)的正演速度逐步提高,單頻點(diǎn)的最大加速比約為2,從而可減少總體正演計(jì)算的時(shí)間。這是因?yàn)殡S著頻率的逐步提高,實(shí)際被用于計(jì)算的網(wǎng)格逐步減少。由此可知,頻率適應(yīng)網(wǎng)格在多頻點(diǎn)正演計(jì)算時(shí),既能加快正演計(jì)算的速度,又能保持正演計(jì)算的精度。

圖3 3層電阻率模型

圖4 非頻率適應(yīng)及頻率適應(yīng)網(wǎng)格正演計(jì)算精度對(duì)比

表1 非頻率適應(yīng)及頻率適應(yīng)網(wǎng)格正演計(jì)算時(shí)間對(duì)比

頻率/Hz1248163264128非頻率適應(yīng)網(wǎng)格/s180155148137125110106102頻率適應(yīng)網(wǎng)格/s1811511239984675952加速比1．001．031．201．381．491．641．801．96

5 算法精度驗(yàn)證及模型算例

5.1 算法精度驗(yàn)證

為了驗(yàn)證本文算法的正確性及精度,我們以圖3 中的3層水平電導(dǎo)率模型為例進(jìn)行有限元正演計(jì)算。由于水平層狀介質(zhì)模型有解析積分形式,通過快速漢克爾變換[14]可以計(jì)算出相應(yīng)的解析解,對(duì)比驗(yàn)證本文算法的正確性及精度。

將沿x方向的單位電偶源置于坐標(biāo)原點(diǎn),在y=2.5km，長為2km的關(guān)于y軸對(duì)稱的測(cè)線上觀測(cè)二次電場(chǎng)的水平分量,測(cè)點(diǎn)間距為100m。基于電導(dǎo)率為100Ω·m的均勻半空間模型計(jì)算得到一次場(chǎng),將該一次場(chǎng)作為場(chǎng)源計(jì)算對(duì)應(yīng)的二次場(chǎng)。采用的剖分網(wǎng)格大小為71×61×50,即x,y及z方向分別剖分為71,61及50個(gè)單元,其中z方向包括20個(gè)空氣層及30個(gè)地下介質(zhì)層。中間的目標(biāo)研究區(qū)域在水平方向上采用較小的網(wǎng)格進(jìn)行均勻剖分,外圍層逐漸擴(kuò)大到約3倍最低頻率時(shí)的趨膚深度。在垂直方向上靠近地面處采用較小的網(wǎng)格剖分,同時(shí)在異常電導(dǎo)率區(qū)域剖分較細(xì),遠(yuǎn)離界面時(shí)采用逐漸增大的網(wǎng)格擴(kuò)展到約3倍趨膚深度處。

圖5a是二次電場(chǎng)x分量的有限元數(shù)值解及解析解對(duì)比,計(jì)算的頻率為1Hz,圖5b為數(shù)值解與解析解之間的相對(duì)誤差。分析圖5a可知,數(shù)值解的實(shí)部及虛部分別與解析解的實(shí)部及虛部吻合。由圖5b可以看出,數(shù)值解實(shí)部的相對(duì)誤差均在1%以下,而虛部的相對(duì)誤差在2%～3%,滿足精度要求。觀察圖5a發(fā)現(xiàn)，實(shí)部數(shù)值遠(yuǎn)大于虛部,這可能是導(dǎo)致虛部誤差較實(shí)部誤差大的原因之一。實(shí)部及虛部的誤差分布呈鋸齒狀則與剖分網(wǎng)格情況有關(guān),因?yàn)樵谒椒较蛏系淖钚【W(wǎng)格間距為200m,而測(cè)點(diǎn)間距為100m,分布在兩個(gè)單元之間的測(cè)點(diǎn)的誤差相對(duì)偏大。

圖6為二次電場(chǎng)y分量的有限元數(shù)值解及解析解的對(duì)比及其相對(duì)誤差。同樣可見有限元數(shù)值解的實(shí)部及虛部分別與解析解的實(shí)部及虛部吻合(圖6a),實(shí)部及虛部的誤差均在1%以下(圖6b),滿足精度要求。

5.2 模型算例

圖7為一個(gè)簡化的儲(chǔ)層模型示意圖,用以研究油氣開采過程中從油氣底部注入鹽水時(shí)地面上的電磁響應(yīng)變化規(guī)律。模型為3層電導(dǎo)率介質(zhì),有一個(gè)代表油氣的高阻梯形區(qū)域埋藏其中。第1層的電阻率為100Ω·m,厚度為0.5km;第2層的電阻率為10Ω·m,厚度為2km;底層(基巖)的電阻率為200Ω·m。發(fā)射源為處于原點(diǎn)并沿x方向的單位電偶極。水平x方向長為6km,y方向頂部寬為3km,底部寬為4km的梯形油藏關(guān)于y軸對(duì)稱,其中心距離發(fā)射源的水平距離為2.8km。油藏頂部距離地面1.5km,油氣及注入鹽水在深度方向的總厚度為0.5km。其中,油氣的電阻率設(shè)為500Ω·m,鹽水的電阻率設(shè)為1Ω·m。當(dāng)發(fā)射源發(fā)射時(shí)諧電磁信號(hào)時(shí),在地面上觀測(cè)電場(chǎng)信號(hào)。

圖5 二次電場(chǎng)x分量精度驗(yàn)證

圖6 二次電場(chǎng)y分量精度驗(yàn)證

圖7 簡單儲(chǔ)層模型

為了研究模型中油藏的電磁響應(yīng)情況,我們?cè)O(shè)計(jì)如下模擬計(jì)算方案:當(dāng)油氣層的厚度分別為400,300及200m,水層厚度分別為100,200及300m時(shí),計(jì)算相應(yīng)3個(gè)時(shí)刻下的二次電場(chǎng)(異常場(chǎng)),并計(jì)算二次場(chǎng)與總場(chǎng)的百分比,采用的頻率為1Hz。這里只展示二次電場(chǎng)是因?yàn)橐淮螆?chǎng)在激發(fā)源附近占主導(dǎo)地位,在總場(chǎng)的分布變化圖中很難區(qū)分不同時(shí)刻的異常場(chǎng)變化。圖8為電場(chǎng)x分量異常場(chǎng)的實(shí)部變化情況,其中圖8a,圖8b,圖8c分別對(duì)應(yīng)油氣厚度400,300及200m時(shí)刻的二次場(chǎng);圖8d,圖8e,圖8f 分別為這3個(gè)不同時(shí)刻二次場(chǎng)占總場(chǎng)的百分比。觀測(cè)圖8可知,異常場(chǎng)的幅度最大達(dá)10-11V·m-1,若實(shí)際儀器的可探測(cè)精度為10-9V·m-1,則將激發(fā)源的電偶極距增大為100A·m時(shí),異常場(chǎng)便處于可探測(cè)的范圍。從圖8a 至圖8c,異常場(chǎng)的分布也發(fā)生了變化,總體上場(chǎng)值呈增大的趨勢(shì),這是由于水層厚度逐漸變大并逐漸靠近激發(fā)源。觀察每一時(shí)刻的二次場(chǎng)分布圖,發(fā)現(xiàn)油氣區(qū)域邊界附近的異常值較大,為勘探油氣的有利觀測(cè)區(qū)域。同時(shí),在實(shí)際勘探中,電磁場(chǎng)數(shù)據(jù)不可避免會(huì)夾雜噪聲,若假設(shè)噪聲水平為5%,則當(dāng)異常場(chǎng)與總場(chǎng)的比值大于5%時(shí)即可發(fā)現(xiàn)異常的存在。觀察圖8d至圖8f 可知,在油氣區(qū)域邊界附近,相對(duì)異常較大(部分區(qū)域遠(yuǎn)大于5%),有利于目標(biāo)儲(chǔ)層的勘探,并且隨著水層厚度增大,有利探測(cè)區(qū)域也增大。圖9為油氣厚度變化時(shí)的異常差值(時(shí)移異常),其中圖9a 為油藏厚度從400m減少到300m時(shí)的異常差值,圖9b為油藏厚度從300m減少到200m時(shí)的異常差值。由圖9 可知,當(dāng)從油氣底部注入水時(shí),油氣在地面上的投影位置區(qū)域可以觀測(cè)到較大的時(shí)移異常,從而為油藏的監(jiān)測(cè)提供指示。

圖8 電場(chǎng)x分量異常場(chǎng)實(shí)部變化

圖9 電場(chǎng)x分量實(shí)部時(shí)移異常

類似地,圖10為相同條件下電場(chǎng)x分量異常場(chǎng)虛部變化圖。觀察可知,異常場(chǎng)虛部的最大幅度達(dá)10-11V·m-1,因此當(dāng)激發(fā)源增大到100A·m時(shí),虛部分量也處于可探測(cè)范圍內(nèi),且在油氣區(qū)域邊界處有較大異常值。比較分析虛部異常場(chǎng)與總場(chǎng)的比值可知,在油氣區(qū)域的邊界及以外區(qū)域,相對(duì)異常較大,處于可探測(cè)范圍；同時(shí),當(dāng)水位上升時(shí),可探測(cè)范圍增大。圖11為油氣厚度變化時(shí)的時(shí)移異常,與實(shí)部情況類似,虛部的時(shí)移異常也可為油藏的變化提供指示信息。

圖10 電場(chǎng)x分量異常場(chǎng)虛部變化

圖11 電場(chǎng)x分量虛部時(shí)移異常

6 結(jié)束語

本文利用矢量有限元方法實(shí)現(xiàn)了三維可控源電磁法的正演模擬。該方法能夠精確模擬電場(chǎng)法向分量在不同電導(dǎo)率界面的不連續(xù)性,且在剖分單元內(nèi)自動(dòng)滿足散度條件，采用的頻率適應(yīng)網(wǎng)格能夠在保持計(jì)算精度的同時(shí),加快多頻點(diǎn)正演計(jì)算速度。模擬計(jì)算結(jié)果表明,可控源電磁法能有效探測(cè)到高阻油藏,顯示了可控源電磁法在油藏監(jiān)測(cè)方面的應(yīng)用前景。

下一步研究工作是將該方法與傳統(tǒng)的節(jié)點(diǎn)型有限元可控源電磁法進(jìn)行對(duì)比，同時(shí)開展三維反演研究,進(jìn)一步檢驗(yàn)頻率適應(yīng)網(wǎng)格在反演中的應(yīng)用效果。

[1] 沈金松,孫文博.二維海底地層可控源海洋電磁響應(yīng)的數(shù)值模擬[J].石油物探,2009,48(2):187-194 Shen J S,Sun W B.The numerical simulation of controlled source marine electromagnetic response for 2-D seabed stratum[J].Geophysical Prospecting for Petroleum,2009,48(2):187-194

[2] 唐新功,胡文寶,嚴(yán)良俊,等.瞬變電磁法油藏動(dòng)態(tài)監(jiān)測(cè)模擬[J].石油物探,2004,43(2):192-195 Tang X G,Hu W B,Yan L J,et al.The application of transient electromagnetic in reservoir monitoring[J].Geophysical Prospecting for Petroleum,2004,43(2):192-195

[3] Newman G A,Alumbaugh D L.Frequency-domain modeling of airborne electromagnetic responses using staggered finite differences[J].Geophysical Prospecting,1995,43(8):1021-1042

[4] Druskin V L,Knizhnerman L A,Lee P.New spectral Lanczos decomposition method for induction modeling in arbitrary 3-D geometry[J].Geophysics,1999,64(3):701-706

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

[6] Badea E A,Everett M E,Newman G A,et al.Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials[J].Geophysics,2001,66(3):786-799

[7] Puzyrev V,Koldan J,de la Puente J,et al.A parallel finite-element method for three-dimensional controlled source electromagnetic forward modelling[J].Geophysical Journal International,2013,193(2):678-693

[8] Streich R.3D finite-difference frequency-domain modeling of controlled-source electromagnetic data:direct solution and optimization for high accuracy[J].Geophysics,2009,74(5):F95-F105

[9] Avdeev D B.Three-dimensional electromagnetic modelling and inversion from theory to application[J].Surveys in Geophysics,2005,26(6):767-799

[10] 張繼鋒,湯井田,喻言,等.基于電場(chǎng)矢量波動(dòng)方程的三維可控源電磁法有限元法數(shù)值模擬[J].地球物理學(xué)報(bào),2009,52(12):3132-3141 Zhang J F,Tang J T,Yu Y,et al.Three-dimensional controlled source electromagnetic numerical simulation based on electromagnetic vector wave equation using finite element method[J].Chinese Journal of Geophysics(In Chinese),2009,52(12):3132-3141

[11] 徐志鋒,吳小平.可控源電磁三維頻率域有限元模擬[J].地球物理學(xué)報(bào),2010,53(8):1931-1939 Xu Z F,Wu X P.Controlled source electromagnetic 3-D in frequency-domain by finite element method[J].Chinese Journal of Geophysics(In Chinese),2010,53(8):1931-1939

[12] 湯井田,張林成,公勁喆,等.三維頻率域可控源電磁法有限元-無限元結(jié)合數(shù)值模擬[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,45(4):1251-1260 Tang J T,Zhang L C,Gong J Z,et al.3D frequency domain controlled source electromagnetic numerical modeling with coupled finite-infinite element method[J].Journal of Central South University(Science and Technology),2014,45(4):1251-1260

[13] Nabighian M N.Electromagnetic methods in applied geophysics:theory[M].Tulsa:Society of Exploration Geophysics,1987:131-138

[14] Guptasarma D,Singh B.New digital linear filters for Hankel J0 and J1 transforms[J].Geophysical Prospecting,1997,45(5):745-762

[15] 金建銘.電磁場(chǎng)有限元方法[M].王建國,譯.第1版.西安:西安電子科技大學(xué)出版社,2001:292-300 Jin J M.Electromagnetic FEM[M].Wang Jianguo,translator.1st ed.Xi’an:Xidian University Press,2001:292-300

[16] Smith J T.Conservative modeling of 3-D electromagnetic fields,part II:biconjugate gradient solution and an accelerator[J].Geophysics,1996,61(5):1319-1324

[17] Freund R W,Nachtigal N M.An implementation of the QMR method based on coupled two-term recurrences[J].SIAM Journal on Scientific Computing,1994,15(2):313-337

(編輯：戴春秋)

Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element using secondary electric field

Tang Wenwu1,2,3,Li Yaoguo4,Liu Jianxin2,3,Liu Chunming2,3

(1.FundamentalScienceonRadioactiveGeologyandExplorationTechnologyLaboratory,EastChinaInstituteofTechnology,Nanchang330013,China; 2.SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China; 3.KeyLaboratoryofNonferrousResourcesandGeologicalHazardExploration,Changsha410083,China; 4.CenterforGravity,Electrical&MagneticStudies,DepartmentofGeophysics,ColoradoSchoolofMines,GoldenCO80401,USA)

As an important complementary tool for seismic survey,the controlled-source electromagnetic (EM) method can be used to monitor reservoirs by detecting the subsurface electrical properties.In order to investigate theoretically feasibility of monitoring reservoirs by the controlled-source EM,a requisite component is to calculate the EM responses of three-dimensional conductivity model excited by the EM sources.Total fields varies rapidly and is difficult to simulate accurately for the region adjacent to the source locations,to this end,the primary field is calculated for a half-space or a layered-earth model as the source term.To bypass the problem that the divergence condition is unsatisfied when using traditional nodal element to model the electrical field,the partial deferential equation (PDE) for the controlled-source EM,based on the secondary electric field,is solved discretely by the edge-based finite element method.To model the responses at multiple frequencies,the frequency-adaptive mesh can be used to increase the speed of forward calculation while maintaining the accuracy.The algorithm is verified by comparing the analytic solution with the numerical solution for a three-layered conductivity model.The modeling results of a simple reservoir model demonstrate the possibility of the controlled-source electromagnetic method in reservoir monitoring.

three-dimensional,controlled-source electromagnetic method,edge-based finite element,forward modeling,reservoir monitoring

2015-04-07;改回日期：2015-08-19。

湯文武(1987—),男,博士,講師,主要研究方向?yàn)殡姶欧〝?shù)值模擬及反演成像。

劉春明(1974—),男,博士,講師,主要研究方向?yàn)殡姶欧碧郊皵?shù)值模擬。

國家自然科學(xué)基金項(xiàng)目(41174103)、國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2014AA06A615)、湖南省自然科學(xué)基金項(xiàng)目(2015JJ2151)聯(lián)合資助。

P631

A

1000-1441(2015)06-0665-09

10.3969/j.issn.1000-1441.2015.06.004