2.5維頻率域可控源電磁法水力壓裂監測正演模擬

余勝紅,唐新功,熊治濤

(1.“油氣資源與勘探技術”教育部重點實驗室(長江大學),湖北武漢430100;2.南方科技大學地球與空間科學系,廣東深圳518055)

水力壓裂是用于非常規油氣資源開采的一項重要技術手段之一,目前已被廣泛應用于煤層氣和頁巖氣以及干熱巖的開采中[1-2]。其實質是將壓裂液以高壓的形式注入到儲層中,致使儲層中原有的裂縫擴展并產生新的裂縫,從而增大儲層滲透率,提高油氣開采效率[3-4]。裂縫的幾何形狀、有效壓裂體積以及擴展范圍是評價水力壓裂效果的主要參數。目前,微地震監測因獲取信息多、操作簡便等優勢成為水力壓裂監測最常用的方法之一[5-6]。然而,當壓裂過程結束后,部分不含支撐劑的裂縫會發生閉合現象,而由于微地震監測與壓裂進程的同步性,無法準確獲取壓裂后的有效壓裂體積參數,導致對產能的預測與評價與實際情況產生偏差[7]。而水力壓裂時所注入壓裂液的電阻率通常比圍巖電阻率低很多,電磁法具有對于地下電阻率變化反應靈敏,尤其對低阻體反應更加靈敏的優勢,為電磁法用于監測水力壓裂裂縫的定位和后續擴展范圍評估提供了物理基礎[8-9]。

目前國內外已有一些學者開展了電磁法監測水力壓裂正反演的研究。YANG等[10]研究了鉆孔電阻率測量中裂縫橫截面積、形狀及傾斜角度等參數對電磁場響應的影響,研究結果表明電磁法對裂縫的面積、形狀及方位的評價具有較好的效果;李洋[11]采用有限元法實現了大型儲層壓裂裂縫三維空間接收響應的快速正演模擬,并模擬了甚低頻(very low frequency,VLF)電磁法非對稱裂縫的發育情況;KYUBO等[12]提出了一種粒子映射反演方法進行井地電磁水力壓裂監測數值模擬,利用磁化率對壓裂后支撐劑的分布進行了成像處理;LI等[13]將鋼套管作為長電極源進行了壓裂液成像的可檢測性和可恢復性研究;CURCIO[14]采用幾種不同的接收裝置對各向異性組合模型進行了水力壓裂監測試驗,結果表明,電磁法監測水力壓裂有助于確定裂縫的幾何形狀;范濤等[15]和陳理等[16]的研究結果表明,電磁數據的反演結果有助于實現對裂縫的識別和定位;由于反演多解性的問題,有學者提出將電磁法監測與微地震法聯合應用的思路,以彌補各自方法的不足,提高水力壓裂監測的精度[17-18]。這些關于電磁法監測水力壓裂的正反演研究結果為電磁法監測水力壓裂的實際應用提供了理論基礎。

近年來,可控源電磁法已被逐漸應用到實際的煤層氣和頁巖氣的水力壓裂監測中。張瑞林等[19]采用瞬變電磁法監測了水力壓裂時壓裂液的流動方向及擴展區域,證明了瞬變電磁法用于監測水力壓裂的可行性;REES等[20]采用大地電磁測深法(magnetotellurics method,MT)對頁巖氣儲層進行了水力壓裂監測,根據獲取的電磁場數據,推測出了裂縫的方向和連通性;YAN等[21]采用時間域可控源電磁法(time domain electromagnetic method,TDEM)對頁巖氣壓裂情況進行了連續監測,通過對電阻率作歸一化處理分析了裂縫的空間展布形態,并與實際的鉆井及地震資料相結合,驗證了TDEM監測頁巖水力壓裂的有效性;ZHANG等[22]分析了電磁法監測評價裂縫參數的原理,通過對磁場信號的分析實現了對裂縫的長度和高度的評估;電磁法除了應用于非常規油氣資源的水力壓裂監測之外,也已成功應用于地熱等資源的水力壓裂監測中[23-24]。電磁法監測水力壓裂的成功案例表明,該方法提高了對裂縫識別定位的精度,因而提高了油氣資源開發的效率和產量。

時間域電磁法是目前采用電磁法進行水力壓裂監測最常用的方法[25-26]。該方法采用井下激發和接收的方式獲得較為精確的結果,但是方法成本高,數據處理與解釋難度大,工作效率較低?？煽卦匆纛l大地電磁法(controlled source audio-frequency magnetotellurics,CSAMT)是一種頻率電磁測深的方法,采用接地導線或不接地線圈作為發射場源,在波區測量正交的電磁場切向分量,并計算其卡尼亞視電阻率[27]。該方法具有操作簡便、工作效率高、信噪比高、探測深度深等優點,目前已廣泛應用于金屬礦產等資源的勘探中[28-31],而在非常規油氣資源勘探開發中的應用相對較少[32-34]。實際的地質結構通常具有二維性,而CSAMT的場源具有三維性,將三維場源二維模型稱作2.5維。相比于二維場源的模擬結果,2.5維的模擬結果更加符合實際地質情況[35]。CSAMT數值模擬的方法主要有積分方程法[36-37]、有限差分法[38-39]和有限元法[40-42]。相比于積分方程法和有限差分法,有限元法在模擬幾何形狀復雜模型方面更具有靈活性,可通過采用非結構網格精確剖分復雜的地電模型,因此在復雜結構的電磁模擬中得到了廣泛應用[43-46]。本文基于2.5維CSAMT有限元正演算法開展水力壓裂正演數值模擬研究,通過設計包含鋼套管和裂縫的二維模型來模擬水力壓裂的過程。對于數據的處理,本文采用對壓裂前、后裂縫的電場響應值進行殘差處理的方法,以獲得壓裂區域的純異常響應,通過對純異常響應特征的分析,闡明了采用CSAMT進行水力壓裂監測的可行性。本文的研究將為CSAMT法在水力壓裂監測中的應用提供理論基礎,并將CSAMT法拓展到水力壓裂監測的應用中。

1 2.5維CSAMT有限元法正演理論

1.1 2.5維CSAMT偏微分方程

準靜態條件下電性源頻率域的Maxwell方程組為(取正時諧eiωt):

(1)

假設異常體的走向沿著y軸方向無限延伸,地下為各向同性介質,即沿著y軸方向電導率σ、磁導率μ、介電常數ε均保持不變,只在x-z平面內發生變化,電偶極子置于y軸與異常體走向平行,且偶極子中心位于坐標原點,將(1)式沿y方向作傅里葉變換后,再聯立方程組即可得到波數域中其它方向電磁場分量的求解表達式:

(2)

(3)

(4)

(5)

將(2)～(5)式分別代入(1)式,可得:

(6)

(7)

1.2 2.5維CSAMT有限元求解過程

1.2.1 有限元計算方程

本文采用伽遼金(Galerkin)法,以(6)式和(7)式為基礎,通過推導相應的有限元計算方程,得到節點電磁場的線性方程組:

(8)

(9)

1.2.2 單元分析

本文采用規則的矩形單元進行網格剖分,考慮到計算速度的問題,在源和觀測區域網格剖分較密,在遠離源和觀測區域網格逐漸稀疏。本文采用4節點雙線性插值單元,第e個單元的節點局部編碼見圖1,雙線性插值單元的插值函數表達式如(10)式所示。

圖1 雙線性差值剖分單元示意

(10)

(8)式和(9)式經坐標域轉換后的表達式為:

(11)

(12)

(11)式和(12)式可簡化為以下形式:

(13)

其中,

(14)

(15)

(16)

(17)

1.2.3 剛度矩陣的合成

(13)式為單元e的電磁場線性方程組,可以寫成如下形式:

(18)

1.2.4 場源的加載

為解決場源的奇異性問題,本文采用偽δ函數來等效場源。當偽δ函數中心位于坐標原點處時,其表達式如下[47]:

(19)

式中:r表示當電偶極子源中心位于坐標原點時,空間任意一點相對于電偶極子源點的位置;參數τ控制著發射源的規模和寬度,其值應設為場源附近網格步長的整數倍,本文取τ=10。由于本文采用4節點矩形網格,因此場源被分配到了緊鄰源點附近的16個網格單元上。將(19)式的源項代入(11)式和(12)式,在等效源加載的節點處作相應的場源加載即可。

1.2.5 反傅里葉變換

(20)

式中:F(x,ky,z,ω)為被積核函數,表示波數域的電磁場分量。根據(2)式至(5)式即可求得波數域的其它電磁場結果,進而可以求得空間域的Ex和Hx。

1.2.6 波數的選取

前文對2.5維CSAMT進行有限元分析及求解過程中均在波數域內進行,要得到空間域的電磁場值還需進行傅里葉逆變換。根據傅里葉逆變換的公式可知,理論上傅里葉逆變換的積分范圍為無窮大,即需要進行無窮多次的正演計算,這是無法實現的?？紤]到計算精度和速度的問題,本文參照了MITSUHATA[47]選取波數的方式,令波數ky取值上限為:

(21)

式中:Δ為區域離散的最小網格步長。波數的選取個數及分布決定著反傅里葉變換的精度。在給定的取值范圍內(ky≤1/Δ),波數選取的個數越多,計算精度越高,同時計算量也越大。為了兼顧計算精度和計算量,一般在每個對數等間隔內取5～10個波數,便能夠保證足夠的計算精度。在本文計算中,波數的取值范圍為10-6～1/Δ,在每個對數間隔內取7個波數,共36個波數。

2 正演數值模擬

2.1 算法有效性驗證

為了驗證2.5維CSAMT正演算法的正確性和精度,首先設計了均勻半空間和一維層狀介質模型進行算法驗證。設置的均勻半空間模型的電阻率為100Ω·m,發射源沿構造走向y方向布設,源中心點的坐標位于原點,長度為100m,電流強度為1A,測點位于x=7044m的位置。設計的3層H型模型的第1層電阻率為100Ω·m,厚度為528m;第2層電阻率為10Ω·m,厚度為200m;第3層電阻率為100Ω·m,其它參數均與均勻半空間模型保持一致,頻率范圍取為10-1～103Hz。圖2和圖3分別為采用本文算法計算的均勻半空間模型和3層H型模型的視電阻率解析解[27]與數值解結果的對比及其相對誤差曲線。從圖2和圖3可以看出,均勻半空間模型和層狀模型的視電阻率在高頻段(遠區)時數值解與解析解的相對誤差較小,二者基本擬合;在中低頻段(近區)時,均勻半空間模型的視電阻率相對誤差在1.8%以內,層狀模型的視電阻率相對誤差在1.5%以內,誤差在可接受范圍內,驗證了本文算法的準確性并保證了足夠的精度。

圖2 均勻半空間模型視電阻率(ρs)解析解與數值解結果對比(a)及相對誤差曲線(b)

圖3 3層H型模型視電阻率解析解與數值解結果對比(a)及相對誤差曲線(b)

為了進一步驗證本文算法的準確性,計算了100Ω·m均勻半空間下距離場源50m的測點的電場強度數值解與解析解結果,如圖4所示。從圖4可以看出,隨著頻率的降低,電場強度迅速下降;當頻率降低到約25Hz后,電場幅值變化緩慢。從圖4還可以看出,數值解與解析解在高頻段誤差較小,而在低頻段二者誤差略微增大,但整體的誤差均小于1.6%,表明了本文算法在近區同樣具有足夠的精度。

圖4 均勻半空間下距離場源50m測點電場幅值數值解與解析解結果對比(a)及相對誤差曲線(b)

在實際應用中,場源下方的地電結構對電磁場的影響規律至關重要。因此,首先設計了在均勻半空間條件下,場源下方存在一個電阻率為10Ω·m、埋深為428m、大小為200m×300m的二維低阻異常體模型,發射源沿構造走向y方向布設,源中心點的坐標位于原點,長度為100m,電流強度為1A。分別計算了場源下存在異常體和不存在異常體時電場和磁場的空間分布,計算結果如圖5所示。從圖5可以看出,當場源下方存在低阻異常體時,在場源附近及異常體所對應的深度(lgf約為1.5Hz)范圍內,電場和磁場會有略微的增強,但由于場源的影響,其變化并不明顯。在遠區,無論場源下方是否存在異常體,二者電場和磁場的響應基本相同。計算結果表明,采用可控源電磁法探測地下低阻裂縫時,與場源下方是否存在異常體的關系不明顯。

圖5 均勻半空間條件下場源下方無異常體和有異常體時電場和磁場的空間分布a 場源下無異常體時電場空間分布; b 場源下存在異常體時電場空間分布; c 場源下無異常體時磁場空間分布; d 場源下存在異常體時磁場空間分布

根據頁巖氣水力壓裂的實際情形和電導率參數,設計了二維模型來模擬水力壓裂的過程,正演計算了不同模型的電場響應結果,分析了CSAMT對水力壓裂裂縫位置的探測情況。在進行數值模擬時,充分考慮了鋼套管的影響。

2.2 發射源靠近井

圖6為二維模型示意圖。在背景電導率σb=0.01S/m的均勻半空間內放置一個直徑為0.1m、由550m的垂直井和一個2000m的水平井組合而成的鋼套管井,以模擬水平井的壓裂過程。沿著y方向布設的電偶極子長度為100m,電流大小為1A,其中心與坐標原點重合,偶極子中心與井口的距離為50m。在進行數值模擬時,鋼套管的電導率設置為106S/m[48],壓裂液為泥質流體,電導率取為3S/m[8]。由于鋼套管和壓裂液的電導率均與背景電導率差異較大,本文將鋼套管和壓裂液作為一個整體,并令其電導率σc=106S/m,空氣電導率σair=10-8S/m。共模擬了3個水平壓裂段,自右至左分別為壓裂段1(stage1)、壓裂段2(stage2)和壓裂段3(stage3)。壓裂的順序從水平井的右端按照stage1、stage2、stage3依次向左進行。其中,stage1裂縫的水平位置為x=2000m;stage2裂縫的水平位置為x=1800m;stage3裂縫的水平位置為x=1600m。將每個壓裂段生成的裂縫等效為一個寬度為0.1m、長度為200m的垂直矩形薄板(圖6),裂縫內支撐劑和壓裂液混合物的有效電導率σf=2500S/m[48-49]。

圖6 發射源靠近井口時的水力壓裂模型示意

測線以50m為測點間距、沿著x方向從井口向右布設,共60個測點,覆蓋長度為2950m。研究區域部分網格剖分如圖7所示。x、z方向的網格數分別為Nx=180、Nz=72,空氣層數量為18層,地層數量為54層。其中場源、鋼套管及裂縫處區域網格剖分較密,其它區域網格剖分較稀疏。場源處網格最小步長為10m;鋼套管處網格最小步長為0.025m;裂縫處網格最小步長為0.02m。挑選了1,10,30,50,80,100Hz共6個頻點的電場響應結果進行展示。

圖7 發射源靠近井口時的水力壓裂模型研究區域網格剖分示意

為了突出每個壓裂段視電阻率的變化,采用(22)式對計算結果進行了相對差異處理:

(22)

式中:ρxi為第i個壓裂段的視電阻率;ρxc為僅含套管(無裂縫)時的視電阻率;R1為二者取對數后的相對差異。

壓裂前、后的計算結果如圖8所示。其中,s1-c表示stage1與僅含套管時視電阻率的相對差異結果;s2-c表示stage2與僅含套管時視電阻率的相對差異結果;s3-c表示stage3與僅含套管時視電阻率的相對差異結果。從圖8可以看出,當頻率從1Hz增加到100Hz時,在對應于1200～2000m的三段裂縫的范圍內,視電阻率相對差異值從負值逐漸過渡到正負相間再過渡到正值,這是因為頻率越低,CSAMT方法的探測深度越深,反之越淺。根據本文設置的模型參數,當頻率分別為1Hz、10Hz和30Hz時,視電阻率的相對差異值較大,此時主要反映了裂縫的視電阻率;當頻率上升為80Hz和100Hz時,視電阻率的相對差異值較小,曲線主要反映了淺部結構的特征,即背景視電阻率響應;當頻率大約為50Hz時,曲線主要反映了二者過渡區域的視電阻率響應特征。從圖8還可以看出,曲線左側壓裂前、后視電阻率的相對差異值均為0,表明在該區域尚未生成裂縫,不同頻率下3條曲線的極值均在1600～2000m,即分別對應著三段壓裂的地面投影位置。當頻率固定時,隨著壓裂進程從最右側的stage1向最左側的stage3移動時,每條曲線對應的極值也具有向左偏移的趨勢,這是因為隨著壓裂過程的進行,左邊不斷產生新的裂縫,曲線的峰值不斷地向新生裂縫的位置偏移,因此從曲線極值點可以定性和動態地反映出新生裂縫的位置。

圖8 發射源接近井口時壓裂前后不同頻率下各壓裂段與僅含套管時視電阻率的相對差異曲線a～f 頻率分別為1Hz、10Hz、30Hz、50Hz、80Hz和100Hz

為了更清晰地闡明CSAMT方法對每個壓裂段裂縫位置的探測情況,對相鄰兩個壓裂段的視電阻率也進行相對差異處理,根據其視電阻率差異曲線進一步分析每個壓裂段裂縫的位置情況。相對差異計算公式為:

(23)

式中:ρxi為第i個壓裂段視電阻率的值;ρx(i-1)為第i-1個壓裂段視電阻率的值;R2為二者取對數后的相對差異。當i=1時,R2=R1。

圖9展示了對裂縫位置反映最明顯的1Hz和10Hz時相鄰兩個壓裂段視電阻率的相對差異曲線。其中,s1-c表示stage1與僅含套管時視電阻率的相對差異結果;s2-s1表示stage2與stage1視電阻率的相對差異結果;s3-s2表示stage3與stage2視電阻率的相對差異結果。從圖9可以看出,隨著壓裂過程的持續,曲線的極小值也明顯具有向左偏移的趨勢,即曲線的極小值與每個壓裂段所產生的裂縫的位置具有良好的對應關系。

圖9 發射源接近井口時1Hz(a)和10Hz(b)頻率下相鄰兩個壓裂段視電阻率的相對差異曲線

2.3 發射源遠離井口

由于CSAMT法通常是在遠區或波區進行觀測,因此設計了如圖10所示的發射源遠離井口的正演模型。與發射源接近井口模型不同的是,此處垂直鋼套管(井口)距發射源的水平距離增大為7600m。本次數值模擬也分為3個壓裂段,壓裂的順序改為從水平井左端依次向右分別按照壓裂段1(stage1)、壓裂段2(stage2)、壓裂段3(stage3)進行壓裂(圖10)。

圖10 發射源遠離井口時的水力壓裂模型示意

其中,stage1裂縫的水平位置為x=6000m;stage2裂縫的水平位置為x=6200m;stage3裂縫的水平位置為x=6400m。128個測點以點距50m沿著x方向從3000m到9350m,長度為6350m,其余參數均與前面的模型一致。網格剖分方式除改變鋼套管和裂縫的位置外,其它網格剖分方式均與圖7一致(圖11)。計算了1Hz,10Hz,30Hz,50Hz,80Hz和100Hz共6個頻點的視電阻率響應結果,采用(22)式對數據進行相對差異處理,計算結果如圖12所示。

圖11 發射源遠離井口時的水力壓裂模型研究區域網格剖分示意

圖12 發射源遠離井口時壓裂前后不同頻率下各壓裂段與僅含套管時視電阻率的相對差異曲線a～f 頻率分別為1Hz,10Hz,30Hz,50Hz,80Hz和100Hz

從圖12可以看出,當頻率從1Hz增加到100Hz時,對應于6000～6400m的裂縫位置范圍內,視電阻率的相對差異值從負值逐漸過渡到正負相間再過渡到正值,隨著頻率的增加,視電阻率的相對差異值也逐漸減小,其原因與圖8相同,不再贅述。從圖12還可以看出,在對應于6000～6400m裂縫上方的收發距位置,不論低頻還是高頻,視電阻率的相對差異均發生了顯著變化。在同一頻率下,隨著壓裂過程的進行,曲線的極值逐漸向右移動,即向裂縫新生成的方向發生偏移,曲線峰值不僅體現了先前壓裂的舊裂縫位置,也體現出了新生裂縫的位置,表明了根據頻率域可控源電磁法視電阻率的相對差異曲線對裂縫形成的位置進行初步定位是完全可行的。

為了更清晰地闡明CSAMT方法對每個壓裂段裂縫位置的探測情況,同樣采用了(23)式對視電阻率進行相對差異處理,即計算了相鄰兩個壓裂段視電阻率的相對差異,根據設置的模型參數及對圖12的分析,此處只展示了對裂縫位置反映最明顯的1Hz和10Hz時視電阻率的相對差異結果,如圖13所示。從圖13可以看出,在6000～7000m范圍內,s2-s1和s3-s2視電阻率的相對差異曲線呈現先增大再減小最后再逐漸增大的趨勢,最終趨近于0。其中,s2-s1視電阻率相對差異曲線的極大值所對應的位置與stage2所產生的裂縫位置(x=6200m)相對應;s3-s2視電阻率的相對差異曲線的極大值所對應的位置則與stage3所產生的裂縫位置(x=6400m)相對應,即相鄰兩個壓裂段視電阻率的相對差異極大值所對應的位置與新生成裂縫的位置一致。

圖13 發射源遠離井口時頻率為1Hz(a)和10Hz(b)時相鄰兩個壓裂段視電阻率的相對差異曲線

為了更加直觀地顯示相鄰兩個壓裂段視電阻率特征的變化,對圖10所示模型相鄰兩個壓裂段的視電阻率直接進行了殘差處理,繪制出相應的擬斷面圖。圖14為發射場源遠離井口時相鄰兩個壓裂段的視電阻率殘差擬斷面圖。圖14涵蓋了3段水平井壓裂的區域。

圖14 發射場源遠離井口時相鄰兩個壓裂段的視電阻率殘差擬斷面a stage1與僅含套管時的視電阻率殘差擬斷面; b stage2與stage1的視電阻率殘差擬斷面; c stage3與stage2的視電阻率殘差擬斷面

從圖14可以看出,隨著壓裂過程的進行,相鄰兩個壓裂段的視電阻率發生了明顯的變化。在x=6000～6400m、f=1～30Hz范圍內,視電阻率的殘差值較大,體現出stage1裂縫所在的大致位置。由于場源的影響會使得異常體沿著遠離場源的一側稍稍被拉伸[50],因此異常區域的中心位置與stage1所產生裂縫的中心位置并不完全一致,而是稍微偏向裂縫的右側(圖14a)。在x=6200～6600m、f=1～30Hz范圍內,視電阻率的殘差值也產生了相當明顯的異常區域,異常區域的中心也對應于stage2所產生的裂縫中心的稍微右側位置(圖14b)。在x=6400～6800m,f=1～30Hz范圍內,視電阻率的殘差值異常也非常明顯,異常區域的中心也對應于stage3所產生的裂縫中心稍微右側的位置(圖14c)。雖然3個壓裂段的視電阻率殘差的峰值與裂縫位置并不完全一致,但差異很小,因此根據視電阻率的殘差響應,可以對每個壓裂段所產生的裂縫進行大致定位。

3 結論

1) 隨著壓裂過程的進行,壓裂前、后的視電阻率會發生變化,視電阻率的相對差異曲線的極值會朝著裂縫新生成的方向偏移,其移動的方向和范圍與新生裂縫的位置具有良好的對應關系。根據視電阻率的相對差異曲線極值的位置可以對裂縫進行大致定位。

2) 在裂縫所對應的位置處,視電阻率產生了足夠大的殘差值,但由于受到場源的影響,異常區域的中心與裂縫的中心并不一一對應,而是稍微偏向裂縫位置的右側。根據相鄰兩個壓裂段視電阻率的殘差結果可以對裂縫的位置進行初步定位。

3) 從計算結果還可以看出,無論是在近區還是遠區觀測,頻率域可控源電磁法的視電阻率對裂縫的響應特征是一致的,即視電阻率的相對差異極值偏移的方向均與裂縫新生成的方向一致。因此,可以根據曲線極值的偏移對新生裂縫的位置進行初步定位。計算結果表明了頻率域可控源電磁法對于水力壓裂裂縫的定位是可行的。本文的研究對于指導頁巖氣與干熱巖等致密儲層的水力壓裂監測具有參考價值。