王新宇 嚴(yán)良俊* 毛玉蓉

(①油氣資源與勘探技術(shù)教育部重點實驗室(長江大學(xué))，湖北武漢 430100； ②非常規(guī)油氣省部共建協(xié)同創(chuàng)新中心，湖北武漢 430100)

0 引言

近年來，隨著油氣資源需求的急劇增加，剩余油氣勘探及油氣藏動態(tài)監(jiān)測新技術(shù)、新方法成為業(yè)內(nèi)研究熱點。時移地球物理方法作為油氣藏動態(tài)監(jiān)測的有效途徑之一，逐漸應(yīng)用于石油開發(fā)領(lǐng)域，取得了良好的效果[1-3]。時移地震技術(shù)是目前發(fā)展最為成熟的時移地球物理方法[4-6]，但對于多數(shù)油氣儲層，油氣藏被流體驅(qū)替引起的聲波阻抗差異較小，致使時移地震監(jiān)測資料解釋困難，且該方法經(jīng)濟成本高，對儲層條件、注采方式等要求嚴(yán)苛。然而，在油氣藏開采過程中，油氣藏被流體驅(qū)替后引起的儲層電性變化明顯，這為時移電磁方法應(yīng)用于油氣藏動態(tài)監(jiān)測提供了地球物理基礎(chǔ)。

隨著地球物理勘探方法的快速發(fā)展，電磁勘探方法逐漸向高維發(fā)展，并成為地震勘探方法的重要補充。電性源瞬變電磁法作為電磁法的重要分支，相比于可控源電磁(CSEM)法和大地電磁(MT)法，具有探測效率高、勘探深度大、信噪比高等優(yōu)點，廣泛應(yīng)用于環(huán)境調(diào)查、油氣、礦產(chǎn)、地?zé)豳Y源勘探及深部地殼研究等領(lǐng)域[7-14]。電磁資料的合理、精細(xì)解釋離不開可靠的正、反演技術(shù)，而正演模擬作為反演的關(guān)鍵步驟，一直是地球物理工作者的研究重點。Edwards等[15]采用頻時轉(zhuǎn)換方法研究了電偶源在海水、基底雙層介質(zhì)中產(chǎn)生的階躍、脈沖響應(yīng)特征。Commer等[16]基于直流電法，通過求解泊松方程得到初始電場，進(jìn)而實現(xiàn)長偏移距瞬變電磁三維正演，并研究了層狀介質(zhì)中油氣儲層的電磁場響應(yīng)特征。Avdeeav等[17]基于有限差分法對比了頻域和時間域海洋可控源電磁法對淺海油氣儲層的探測能力，發(fā)現(xiàn)時間域方法具有更好的勘探效果。Um等[18]采用隱式時間步長，基于一階后退歐拉算法實現(xiàn)了非結(jié)構(gòu)化網(wǎng)格的電性源瞬變電磁三維正演，并提出了一種自適應(yīng)時間步長算法提高計算精度，極大地促進(jìn)了隱式時間步長方法在瞬變電磁非結(jié)構(gòu)有限元三維模擬中的應(yīng)用。此后，隨著地電模型復(fù)雜程度的增加，基于隱式時間步長的非結(jié)構(gòu)有限元瞬變電磁法正、反演得到快速發(fā)展，并廣泛應(yīng)用于航空、半航空、地面、井中、海洋時域電磁法的數(shù)值計算[19-26]。

時移電磁法作為一種有效的地下介質(zhì)電性變化監(jiān)測手段，在油氣藏、地?zé)?、地下水動態(tài)監(jiān)測等方面具有良好的應(yīng)用前景，相關(guān)研究取得了一定的進(jìn)展，但整體上其理論方法、儀器研發(fā)、采集方式、數(shù)據(jù)解釋等仍處于初始階段。

考慮到時域電磁法勘探能力優(yōu)于頻域電磁法，本文基于非結(jié)構(gòu)化網(wǎng)格矢量有限元法進(jìn)行電性源瞬變電磁三維正演，并采用一種二階后退歐拉算法(BDF2)變步長差分格式，避免了迭代過程中需向前查找上一個時間點。該方法實現(xiàn)簡單，理論上可更準(zhǔn)確地計算電場對時間的導(dǎo)數(shù)項，計算精度較高。本文對二階后退歐拉法定步長、變步長差分格式的電性源瞬變電磁數(shù)值結(jié)果與解析解結(jié)果進(jìn)行比對，對三維復(fù)雜地電模型數(shù)值結(jié)果與有限體積法計算結(jié)果進(jìn)行對比分析，以驗證本文算法的有效性與精度。最后，通過計算復(fù)雜油藏模型及實際頁巖氣壓裂模型的動態(tài)監(jiān)測平行電場分量，分析相對異常的變化，驗證了電性源瞬變電磁法對陸地油氣藏監(jiān)測的有效性。

1 正演算法

1.1 非結(jié)構(gòu)網(wǎng)格時域矢量有限元法

對于設(shè)定有限區(qū)域的地電模型，時域電磁法電場擴散方程為[18]

(1)

式中：E(r,t)和js(r,t)分別為時刻t、任意位置r上的電場和發(fā)射源電流密度，下標(biāo)“s”表示外部電流源；μ為介質(zhì)磁導(dǎo)率；σ為電導(dǎo)率。

采用非結(jié)構(gòu)四面體網(wǎng)格對地電模型進(jìn)行剖分(圖1)，并引入矢量插值基函數(shù)，將自由度賦予棱邊上，則任意四面體單元內(nèi)的電場可展開為

(2)

(3)

圖1 四面體單元e的矢量電場分布圖

利用Galerkin方法離散式(1)，得到單元e的殘差矢量

(4)

并確保該單元的加權(quán)余量為零

(5)

將所有四面體單元加權(quán)余量組合，可得

(6)

式中：V表示四面體單元體積；M是四面體單元總數(shù)；A表示單元質(zhì)量矩陣；B表示單元剛度矩陣；S表示電流源項。對于四面體單元e，Ae、Be和Se的積分表達(dá)式分別為

(7)

(8)

(9)

式中i、j均取1～6，為單元棱邊索引。

利用非結(jié)構(gòu)四面體網(wǎng)格離散靈活的特點，將長導(dǎo)線源分解為多段，每段導(dǎo)線可近似為電偶極子，每個電偶極子的電流密度可表示為[28]

js(r,t)=δ(r-rs)I(t)dl

(10)

式中：δ是脈沖函數(shù)；rs為電偶源位置；I是電流矢量； dl是電偶極子長度。

1.2 二階后退歐拉變步長差分格式

求解式(6)需對時間進(jìn)行離散，本文采用精度較高的二階后退歐拉法[18]，經(jīng)Taylor展開得到

(11)

式中：Δt為時間步長；k為迭代次數(shù)。將式(11)代入式(6)可得

(3A+2ΔtB)E(t)(k)=A[4E(t)(k-1)-E(t)(k-2)]-2ΔtS(k)

(12)

上式左端項系數(shù)矩陣與時間步長Δt相關(guān)。當(dāng)使用直接求解器(Pardiso)時，若Δt保持不變，只需更新方程右端項，重復(fù)將儲存的矩陣分解結(jié)果帶回右端項求解(圖2a)； 若Δt是變化的，需向前查找上一個時刻nΔt1的電場值，重新分解系數(shù)矩陣并帶入右端項求解。本文采用一種二階后退歐拉變步長差分格式(圖2b)[29]，此方法僅需利用第k個時間道的前兩個時間步長(Δt1、Δt2)的結(jié)果進(jìn)行計算(為更直觀地表示變步長差分格式，用Δt2替代定步長差分格式的Δt)，具有較高的精度和良好的穩(wěn)定性，比使用定步長差分格式更有效，其具體表達(dá)式為

(13)

將式(13)代入式(6)，得

(14)

圖2 時間步長離散示意圖(a)定步長差分格式； (b)變步長差分格式

當(dāng)Δt1=Δt2時，式(13)和式(14)退化為與式(11)和式(12)相同的結(jié)果，后文對該差分格式的效果進(jìn)行了驗證。

1.3 初始電場計算

(15)

直流電法的電位u滿足三維Poisson方程

?·(σ?u)=-Iδ(r-r0)

(16)

式中：r0為發(fā)射源位置；I為電流強度。

為保證空氣電場切向分量不為零，采取總場方法求解式(16)，并采用與時間域電磁場同一套非結(jié)構(gòu)化網(wǎng)格及對應(yīng)的Dirichlet邊界條件

u|Γ=0

(17)

2 數(shù)值模擬

2.1 解析解驗證

本算例開展正演計算的工作站配置為：處理器AMD Ryzen 9-5950，CPU主頻3.4GHz，內(nèi)存96GB。為驗證本文算法的正確性，對比后退歐拉法定步長與變步長差分格式的計算精度。設(shè)計一個均勻半空間模型進(jìn)行數(shù)值計算，均勻半空間的介質(zhì)電阻率為100Ω·m，空氣電阻率為1×108Ω·m，電性源長度為200m，沿y方向布設(shè)，發(fā)射電流為1A。測點也沿y方向布設(shè)，與源的偏移距為1000m。將長導(dǎo)線等間距分割為100段電偶極子，為保證計算精度，對發(fā)射源與接收點處的網(wǎng)格進(jìn)行局部加密，最終生成201906個四面體、32983個節(jié)點、235656條棱邊(圖3)。計算時間為2×10-7～4s，離散為1122個時間道，分別計算二階后退歐拉法定步長、變步長差分格式在測點處平行電場分量Ey的響應(yīng)曲線，并與解析解對比。利用下式計算相對誤差

(18)

結(jié)果見圖4。式中：EN表示有限元數(shù)值解；EA表示解析解或有限體積解。

由圖4可見，兩種差分格式的計算結(jié)果均與解析解擬合較好，變步長法在早期的計算精度高于定步長法，晚期基本一致，說明變步長差分格式具有較高的精度和穩(wěn)定性。整體而言，兩種方法均可以高精度地計算電場響應(yīng)，表明本文算法的可行性，適用于電性源瞬變電磁的三維正演模擬。

圖3 發(fā)射源區(qū)域(左)與測點區(qū)域(右)網(wǎng)格局部加密示意圖

圖4 不同步長差分格式Ey計算結(jié)果(a)及其與數(shù)值解的相對誤差(b)

2.2 三維模型驗證

為驗證二階后退歐拉法變步長差分格式對復(fù)雜地電模型的計算精度，設(shè)計一個三維模型[20]進(jìn)行計算，模型切面如圖5a所示。電性源長度為1000m，沿y方向布設(shè)，中心點位于地表(0,0,0)處，發(fā)射電流為1A。將長導(dǎo)線等間距分割為400段電偶極子。沿x負(fù)方向布置12個測點，對發(fā)射源和測點區(qū)域進(jìn)行局部網(wǎng)格加密，最終生成1255209個網(wǎng)格、202819個節(jié)點、1458795條棱邊(圖5b)。

計算該模型占用內(nèi)存10.5G，總耗時為3296s。圖6a為測點A(-500m,0,0)、測點B(-1050m,0,0)、測點C(-2000m,0,0)的電場分量Ey數(shù)值計算結(jié)果與Liu等[20]有限體積計算結(jié)果對比。由圖可見，二階后退歐拉法變步長差分格式的矢量有限元解與有限體積解吻合良好。圖6b為誤差曲線，可見這三個測點的計算結(jié)果相對誤差均低于5%，證明了本文算法的正確性，可用于電性源瞬變電磁法油氣藏動態(tài)監(jiān)測數(shù)值研究。

2.3 油藏動態(tài)監(jiān)測模型響應(yīng)分析

為從理論上研究電性源瞬變電磁法對油藏動態(tài)監(jiān)測的響應(yīng)能力，設(shè)計圖7所示的復(fù)雜油藏模型?？諝怆娮杪试O(shè)為1×108Ω·m。發(fā)射源沿y向布設(shè)，長度為1000m，發(fā)射電流為1A。將長導(dǎo)線源等間距分割為400段電偶極子。在地表沿y方向從-1500m到1500m共布置31個測點，點距為100m，測線與源的偏移距為2000m。在近地表存在兩個異常體作為對油氣藏動態(tài)監(jiān)測過程的干擾區(qū)域：一個是長方體異常體，大小為1000m×1000m×200m，傾角為20°，異常體中心埋深為400m，電阻率為2000Ω·m； 另一個是橢球狀異常體，x、y、z方向的半軸長度分別為200、200、100m，中心埋深為200m，電阻率為10Ω·m。假設(shè)一個金字塔狀的油藏，基底邊長為2000m，頂部邊長為1000m，高為900m，油藏頂面距地表2000m，電阻率為500Ω·m。模型背景為三層復(fù)雜起伏地層，從上至下電阻率分別為1000、200、100Ω·m。對發(fā)射源、測點、異常體周圍的網(wǎng)格進(jìn)行局部加密，最終生成1547756個網(wǎng)格、249739個節(jié)點、1798262條棱邊(圖7b)。按照圖中序號，分三個階段進(jìn)行電性源瞬變電磁法監(jiān)測，每個階段鹽水驅(qū)油后儲層電阻率變?yōu)?0Ω·m，由下至上每次驅(qū)油深度300m。該模型計算占用內(nèi)存13.1G，總求解時間為3823s。

圖5 三維復(fù)雜模型切面圖(a)和測點區(qū)域網(wǎng)格剖分示意圖(b)

圖6 電性源瞬變電磁法三維模型Ey分量數(shù)值計算結(jié)果(a)及三個測點的相對誤差曲線(b)

圖7 復(fù)雜地質(zhì)油藏模型3D示意圖(左)和局部網(wǎng)格剖分示意圖(右)右圖中序號①～③為鹽水驅(qū)替油層順序，長方體與橢球體為淺地表異常干擾區(qū)域

圖8為利用以下公式計算的三段水驅(qū)油動態(tài)監(jiān)測相對異常

式中：Ebd、Eal分別表示水驅(qū)油或壓裂前、后的電場響應(yīng)。由圖可見，隨著水驅(qū)油過程中地下介質(zhì)電阻率的不斷降低，電場響應(yīng)的相對異常Ra不斷增大，最大可達(dá)24.00%(絕對值)，且相對異?；静皇艿叵缕鸱缑婕捌渌惓ｓw的干擾，可以清晰地反映油藏動態(tài)監(jiān)測過程。還可以看出，相對異常區(qū)域的頂邊界隨鹽水驅(qū)油過程在時間道上向上移動。根據(jù)電磁理論，電性源瞬變電磁數(shù)據(jù)的早期時間道勘探深度小、晚期時間道勘探深度大，時間道與深度由晚期到早期、由深至淺互相對應(yīng)，很好地反映了鹽水驅(qū)油由下至上的動態(tài)過程。由于近地表異常體產(chǎn)生的電性干擾沒有發(fā)生變化，油藏監(jiān)測過程中干擾區(qū)域的電磁場可作為背景場(Ebd)從總場中減去。因此，近地表異常干擾對油藏動態(tài)監(jiān)測基本沒有影響(圖8)。此例計算結(jié)果表明，利用時移電磁法對油藏進(jìn)行動態(tài)監(jiān)測具有良好的應(yīng)用前景。

圖8 水驅(qū)油動態(tài)監(jiān)測過程階段①～③(左，中，右)相對異常Ra剖面

2.4 實際應(yīng)用

基于涪陵頁巖氣田區(qū)焦頁30井測井資料，Liu等[30]設(shè)計了圖9所示的頁巖氣儲層模型，并采用可控源電磁法在理論上對該地區(qū)頁巖氣壓裂監(jiān)測進(jìn)行了可行性論證。本文基于Liu等建立的模型，采用電性源瞬變電磁法對頁巖氣壓裂模型的動態(tài)監(jiān)測效果進(jìn)行分析。

發(fā)射源沿y方向布設(shè)，中心點位于地表(0,0,0)，長度為1000m，發(fā)射電流為100A。將長導(dǎo)線源等間距分割為400段電偶極子。在地表沿y方向-800～800m區(qū)域共布置41個測點，間距為40m，測線偏移距為5000m。圖9中頁巖氣規(guī)模為2300m×840m×300m，電阻率為42Ω·m，壓裂后電阻率變?yōu)?Ω·m。對發(fā)射源、測點、頁巖氣藏區(qū)域進(jìn)行局部網(wǎng)格加密，最終生成1914330個網(wǎng)格、305778個節(jié)點和2221195條棱邊。

該模型計算占用內(nèi)存14.6G，總求解時間為4075s。圖10a為測線中點(5000，0，0)壓裂前電場Ey

圖9 涪陵焦頁30井頁巖氣儲層電阻率模型

圖10 頁巖氣層狀儲層模型測點(5000m，0，0)壓裂前、后的Ey曲線及相對誤差(a)壓裂前有限元解和解析解及壓裂后的響應(yīng)曲線； (b)壓裂前有限元解與解析解的相對誤差； (c)壓裂后響應(yīng)與壓裂前有限元解的相對異常

的解析解、有限元數(shù)值解(均勻?qū)訝顑幽Ｐ?及頁巖氣壓裂后的電場Ey； 圖10b為壓裂前均勻?qū)訝顑幽Ｐ陀邢拊獢?shù)值解相對于解析解的相對誤差曲線，可見誤差不大于3.50%，具有較高的計算精度； 圖10c為壓裂前、后有限元數(shù)值解的相對異常，可以看出相對異常最大達(dá)27.39%，可以清楚地反映出頁巖氣藏壓裂引起的電性變化。

圖11a為頁巖氣儲層壓裂前、后測線上不同時間道的電場Ey分量相對異常等值線，可以明顯看出壓裂引起的電性異常，但同時也伴隨著假異常(圖中虛線框區(qū)域)的出現(xiàn)，這是由于壓裂過程中，隨著壓裂液的侵入，壓裂區(qū)的電阻率會降低，電流在低阻體中形成明顯的電流通道效應(yīng)，致使壓裂區(qū)域電場出現(xiàn)較大變化，最終表現(xiàn)為正、負(fù)異常特征。圖11b為壓裂前、后總場殘差(Eal-Ebd)等值線圖，同樣可以看出頁巖氣藏壓裂導(dǎo)致的電場殘差出現(xiàn)正、負(fù)特征。總之，頁巖氣藏壓裂前、后電場異常響應(yīng)明顯，證實了利用時移電磁法對油氣藏動態(tài)監(jiān)測的有效性和可行性。

圖11 頁巖氣層狀儲層壓裂前、后Ey相對異常Ra(a)及總場殘差(Eal-Ebd)(b)

3 結(jié)論

本文將非結(jié)構(gòu)化網(wǎng)格矢量有限元法應(yīng)用于電性源瞬變電磁法油氣藏動態(tài)監(jiān)測的正演模擬，得到以下結(jié)論。

(1)采用二階后退歐拉法變步長差分格式離散時間步長，可有效提高數(shù)值計算精度。對比復(fù)雜地電模型的有限體積解，有限元解可有效保證數(shù)值結(jié)果的穩(wěn)定性，適于復(fù)雜油氣藏模型動態(tài)監(jiān)測數(shù)值模擬。

(2)對理論油藏模型與實際頁巖氣模型的數(shù)值分析表明，電性源時移電磁法的相對異常響應(yīng)基本不受地層中其它異常體的影響，可有效刻畫油(氣)藏驅(qū)替(壓裂)前、后的電性差異，實現(xiàn)高精度油氣藏動態(tài)監(jiān)測。

電性源時移電磁法對油氣藏、地下水、金屬等資源動態(tài)監(jiān)測是有效的。作為時移地震勘探的重要補充，時移電磁法目前還處于起步階段，借助本文算法可分析電性源瞬變電磁法對礦產(chǎn)資源動態(tài)監(jiān)測的規(guī)律，為實際礦產(chǎn)資源動態(tài)監(jiān)測提供理論指導(dǎo)。