電性源瞬變電磁法油氣藏動態監測模擬分析

王新宇 嚴良俊* 毛玉蓉

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

0 引言

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

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

時移電磁法作為一種有效的地下介質電性變化監測手段，在油氣藏、地熱、地下水動態監測等方面具有良好的應用前景，相關研究取得了一定的進展，但整體上其理論方法、儀器研發、采集方式、數據解釋等仍處于初始階段。

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

1 正演算法

1.1 非結構網格時域矢量有限元法

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

(1)

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

采用非結構四面體網格對地電模型進行剖分(圖1)，并引入矢量插值基函數，將自由度賦予棱邊上，則任意四面體單元內的電場可展開為

(2)

(3)

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

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

(4)

并確保該單元的加權余量為零

(5)

將所有四面體單元加權余量組合，可得

(6)

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

(7)

(8)

(9)

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

利用非結構四面體網格離散靈活的特點，將長導線源分解為多段，每段導線可近似為電偶極子，每個電偶極子的電流密度可表示為[28]

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

(10)

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

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

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

(11)

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

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

(12)

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

(13)

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

(14)

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

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

1.3 初始電場計算

(15)

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

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

(16)

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

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

u|Γ=0

(17)

2 數值模擬

2.1 解析解驗證

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

(18)

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

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

圖3 發射源區域(左)與測點區域(右)網格局部加密示意圖

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

2.2 三維模型驗證

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

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

2.3 油藏動態監測模型響應分析

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

圖5 三維復雜模型切面圖(a)和測點區域網格剖分示意圖(b)

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

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

圖8為利用以下公式計算的三段水驅油動態監測相對異常

式中：Ebd、Eal分別表示水驅油或壓裂前、后的電場響應。由圖可見，隨著水驅油過程中地下介質電阻率的不斷降低，電場響應的相對異常Ra不斷增大，最大可達24.00%(絕對值)，且相對異常基本不受地下起伏界面及其他異常體的干擾，可以清晰地反映油藏動態監測過程。還可以看出，相對異常區域的頂邊界隨鹽水驅油過程在時間道上向上移動。根據電磁理論，電性源瞬變電磁數據的早期時間道勘探深度小、晚期時間道勘探深度大，時間道與深度由晚期到早期、由深至淺互相對應，很好地反映了鹽水驅油由下至上的動態過程。由于近地表異常體產生的電性干擾沒有發生變化，油藏監測過程中干擾區域的電磁場可作為背景場(Ebd)從總場中減去。因此，近地表異常干擾對油藏動態監測基本沒有影響(圖8)。此例計算結果表明，利用時移電磁法對油藏進行動態監測具有良好的應用前景。

圖8 水驅油動態監測過程階段①～③(左，中，右)相對異常Ra剖面

2.4 實際應用

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

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

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

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

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

的解析解、有限元數值解(均勻層狀儲層模型)及頁巖氣壓裂后的電場Ey； 圖10b為壓裂前均勻層狀儲層模型有限元數值解相對于解析解的相對誤差曲線，可見誤差不大于3.50%，具有較高的計算精度； 圖10c為壓裂前、后有限元數值解的相對異常，可以看出相對異常最大達27.39%，可以清楚地反映出頁巖氣藏壓裂引起的電性變化。

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

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

3 結論

本文將非結構化網格矢量有限元法應用于電性源瞬變電磁法油氣藏動態監測的正演模擬，得到以下結論。

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

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

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