李林， 張廣智*

1 中國石油大學(xué)(華東)深層油氣重點實驗室， 青島 266580 2 中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580

0 引言

隨著能源需求的增加和常規(guī)油氣資源的減少，非常規(guī)油氣資源(頁巖氣、致密氣和致密油儲層)成為石油行業(yè)的關(guān)鍵勘探目標.非常規(guī)儲層具有低孔低滲的特點，需要通過大規(guī)模的水力壓裂以提高油氣產(chǎn)量.地下裂縫能夠有效改善儲層的滲透率，地應(yīng)力是控制水力壓裂的重要因素(陳志剛等，2020；王生奧等，2021)，因此利用高品質(zhì)“五維”地震數(shù)據(jù)開展裂縫參數(shù)及地應(yīng)力參數(shù)地震反演對于地下裂縫預(yù)測及水力壓裂具有重要意義(印興耀等，2018a，b).

早期的地應(yīng)力預(yù)測模型假設(shè)地層為各向同性，將各向同性地應(yīng)力預(yù)測模型應(yīng)用于各向異性地層會導(dǎo)致地應(yīng)力預(yù)測不準確，因此各向異性地應(yīng)力模型更加符合非常規(guī)儲層的實際情況(趙小龍等，2020).許多學(xué)者考慮頁巖具有各向異性的微觀結(jié)構(gòu)，提出了適用于垂直橫向各向同性(Vertical Transverse Isotropic, VTI)介質(zhì)的地應(yīng)力預(yù)測模型(Thiercelin and Plumb，1994；Higgins et al.，2008；鄧金根等，2013).張廣智等(2015)提出利用頁巖巖石物理等效模型進行地應(yīng)力預(yù)測的方法.考慮到頁巖儲層中高角度裂縫較為發(fā)育的特征，Gray等(2012)推導(dǎo)了水平橫向各向同性(Horizontal Transverse Isotropic, HTI)介質(zhì)中的水平應(yīng)力預(yù)測方程，并指出將水平應(yīng)力差異比(Differential Horizontal Stress Ratio，DHSR)作為衡量儲層是否能壓裂成網(wǎng)的指示因子.Far等(2016)研究了兩組正交的垂直裂縫嵌入VTI背景誘導(dǎo)的正交各向異性介質(zhì)中的水平應(yīng)力預(yù)測模型.馬妮等(2017，2018)推導(dǎo)了單組垂直裂縫嵌入VTI背景誘導(dǎo)的正交各向異性介質(zhì)中的水平應(yīng)力預(yù)測方程.

彈性阻抗(Elastic Impedance, EI)在預(yù)測彈性參數(shù)方面有著廣泛的應(yīng)用前景(Connolly，1999；Whitcombe，2002；Wang et al.，2006；印興耀等，2013；Zong et al.，2013).Martins(2006)將彈性阻抗擴展到各向異性介質(zhì)中，隨后許多學(xué)者研究了利用各向異性彈性阻抗提取各向異性參數(shù)的反演方法(陳懷震等，2014；Chen et al., 2018；Pan et al.，2020b；Li et al.，2020).Whitcombe等(2002)提出擴展彈性阻抗(Extended Elastic Impedance, EEI)的概念，通過擴展彈性阻抗可轉(zhuǎn)換為其他彈性參數(shù)、物性參數(shù)或地質(zhì)力學(xué)參數(shù)，進而用于巖性、流體識別及地應(yīng)力預(yù)測(Hafez and Castagna，2016；Aleardi，2018；Sharifi et al.，2019；Sharifi and Mirzakhanian，2019).研究表明深度、壓實趨勢、厚度和巖性等因素會影響EEI分析中的相關(guān)性趨勢(Ball et al.，2013；Thomas et al.，2013).Connolly(2019)指出彈性參數(shù)與擴展彈性阻抗的相關(guān)性取決于橫縱波速度平方比，且由于反演存在的誤差，最優(yōu)的旋轉(zhuǎn)角應(yīng)通過對反演結(jié)果而不是測井數(shù)據(jù)進行EEI分析得到.Russell和Hedlin(2019)基于Biot-Gassmann多孔彈性理論提出了擴展多孔彈性阻抗的概念，并指出擴展多孔彈性阻抗與彈性參數(shù)的相關(guān)性取決于干巖石的橫縱波速度平方比.

考慮單組垂直裂縫誘導(dǎo)的方位各向異性的影響，筆者將擴展彈性阻抗推廣到HTI介質(zhì)中，提出擴展方位彈性阻抗的概念.首先推導(dǎo)HTI介質(zhì)中的擴展方位彈性阻抗方程及其傅里葉級數(shù)展開表達式；基于推導(dǎo)的傅里葉系數(shù)(Fourier Coefficient，F(xiàn)C)方程，分析傅里葉系數(shù)與裂縫參數(shù)及DHSR之間的關(guān)系，進而提出一種基于擴展方位彈性阻抗的裂縫參數(shù)及DHSR預(yù)測方法，最后通過模型測試和實際應(yīng)用驗證該方法在裂縫參數(shù)及DHSR預(yù)測方面的可靠性.

1 方法原理

1.1 HTI介質(zhì)擴展方位彈性阻抗方程及傅里葉系數(shù)

HTI介質(zhì)中的縱波反射系數(shù)方程為(潘新朋等，2018)：

+aδN(θ,φ)ΔδN+aδT(θ,φ)ΔδT，

(1)

式(1)可重寫為：

RPP(θ,φ)=A+B(φ)sin2θ+C(φ)sin2θtan2θ，

(2)

式中：

(3)

(4)

(5)

Whitcombe等(2002)推導(dǎo)了擴展彈性阻抗方程，通過改變單一的旋轉(zhuǎn)角即可由擴展彈性阻抗轉(zhuǎn)化為其他彈性參數(shù).類比各向同性介質(zhì)中的擴展彈性阻抗推導(dǎo)過程，將tan=sin2θ代入式(2)，左右兩邊同時乘以cos，得到：

Rcos=Acos+B(φ)sin

(6)

RPP(,φ)=Rcos(cos-sin)

=Ap()+B(φ)q()+C(φ)r()，

(7)

其中，p()=cos(cos-sin)，q()=sin(cos-sin)，r()=sin2.RPP(,φ)表示尺度化的方位反射系數(shù).

彈性阻抗的定義為(Connolly, 1999)：

(8)

結(jié)合式(7)、(8)，通過推導(dǎo)可得到HTI介質(zhì)中的擴展方位彈性阻抗方程為：

EEI(,

(9)

(10)

(11)

(12)

對式(9)兩邊同時取對數(shù)，得到歸一化的擴展方位彈性阻抗方程：

LEEI(,φ)=p()LAI+q()LBI(φ)+r()LCI(φ)，

(13)

其中，LEEI(,

對式(13)進行傅里葉級數(shù)展開，得到：

LE E I(,φ)=R0()+R2()cos(2φ)+R4()cos(4φ)，

(14)

R0()

(15)

R2()=Baniq()+[g(g-1)δN]r()，

(16)

R4()，

(17)

式中，Rn(),(n=0,2,4)表示與旋轉(zhuǎn)角相關(guān)的傅里葉系數(shù)，Bani=g(δT-kδN)表示各向異性梯度.

傅里葉系數(shù)可通過式(18)計算：

Rn(,φi)cos(nφ)dφ/cos(nφsym)

(18)

1.2 傅里葉系數(shù)與裂縫參數(shù)及DHSR的關(guān)系

Gray等(2012)推導(dǎo)了楊氏模量、泊松比及法向柔度表征的DHSR表達式，通過推導(dǎo)可得到由縱、橫波阻抗(IP和IS)及法向弱度表征的DHSR，其表達式為：

(19)

式(19)表明，DHSR僅與地層系數(shù)g和法向弱度δN有關(guān).根據(jù)式(15)、(16)、(17)可知，零階傅里葉系數(shù)與各向同性參數(shù)及裂縫參數(shù)均有關(guān)，二階和四階傅里葉系數(shù)主要與裂縫參數(shù)有關(guān).通過變化角，四階傅里葉系數(shù)只能轉(zhuǎn)化為而不能轉(zhuǎn)化為裂縫參數(shù)和DHSR，因此下面以某工區(qū)A井數(shù)據(jù)為例，分析零階和二階傅里葉系數(shù)與各向同性參數(shù)、裂縫參數(shù)及DHSR之間的關(guān)系.

圖1為A井中的各向同性參數(shù)及裂縫參數(shù)，包括縱、橫波速度、密度、法向及切向裂縫弱度.結(jié)合式(15)，利用井數(shù)據(jù)合成不同角情況下的零階傅里葉系數(shù)如圖2a所示，可以看出零階傅里葉系數(shù)隨著角的變化而變化.由于零階傅里葉系數(shù)與各向同性參數(shù)和裂縫參數(shù)均相關(guān)，因此分別計算零階傅里葉系數(shù)與各向同性參數(shù)、裂縫參數(shù)及DHSR的相關(guān)系數(shù)，如圖2b所示，可以看出縱、橫波速度、密度、法向弱度、切向弱度及DHSR與零階傅里葉系數(shù)的相關(guān)系數(shù)隨著角的變化而變化，縱、橫波速度、密度、法向弱度、切向弱度及DHSR與零階傅里葉系數(shù)的最大相關(guān)系數(shù)絕對值分別為0.96、0.99、0.83、0.72、0.74和0.74.可以看出法向弱度、切向弱度及DHSR與零階傅里葉系數(shù)的相關(guān)系數(shù)遠低于各向同性參數(shù)與零階傅里葉系數(shù)的相關(guān)系數(shù).李林等(2020)研究表明，裂縫參數(shù)對零階傅里葉系數(shù)的貢獻遠小于各向同性參數(shù)的貢獻，導(dǎo)致零階傅里葉系數(shù)對裂縫參數(shù)的敏感性較低，因此零階傅里葉系數(shù)僅可用于估測各向同性參數(shù)，而不能用于裂縫參數(shù)和DHSR估測.

圖1 A井的各向同性參數(shù)及裂縫參數(shù)Fig.1 Isotropic parameters and fracture parameters of well A

圖2 合成的(a)零階傅里葉系數(shù)，(b)零階傅里葉系數(shù)與各向同性參數(shù)、裂縫參數(shù)及DHSR的相關(guān)系數(shù)Fig.2 Synthetic (a) zeroth FC, and (b) correlation coefficients of the zeroth FC with isotropic parameters, fracture parameters and DHSR

圖3 合成的(a)二階傅里葉系數(shù)，(b)二階傅里葉系數(shù)與裂縫參數(shù)和DHSR的相關(guān)系數(shù)Fig.3 Synthetic (a) second FC, and (b) correlation coefficients of the second FC with fracture parameters and DHSR

R2(

(20)

R2(=90°)=-Bani+g(g-1)δN.

(21)

結(jié)合式(20)、(21)可得到各向異性梯度表達式：

Bani=2R2(=45°)-R2(=90°).

(22)

在疊前地震反演中，通常假設(shè)一段時窗內(nèi)的地層系數(shù)g為一常數(shù)，例如假設(shè)g=0.25(Whitcombe et al.，2002；Russell and Hedlin，2019).結(jié)合各向異性梯度與裂縫參數(shù)的關(guān)系表達式，通過推導(dǎo)可得到：

(23)

(24)

結(jié)合式(19)、(23)可得到HTI介質(zhì)中的DHSR簡化表達式：

(25)

因此，當假設(shè)地層系數(shù)g為0.25時，利用式(23)、(24)、(25)即可由二階傅里葉系數(shù)分別估測HTI介質(zhì)中的裂縫參數(shù)及DHSR.Connolly(2019)指出彈性參數(shù)與擴展彈性阻抗的相關(guān)性取決于地層系數(shù)g.本文研究表明，二階傅里葉系數(shù)與裂縫參數(shù)及DHSR的相關(guān)性也與地層系數(shù)g相關(guān).由于地層系數(shù)g是隨著時間或深度變化的，因此在實際應(yīng)用中應(yīng)通過尋找二階傅里葉系數(shù)與裂縫參數(shù)及DHSR的最大相關(guān)系數(shù)，分別確定最優(yōu)的角進行裂縫參數(shù)及DHSR預(yù)測.

1.3 分方位疊前地震貝葉斯反演

由于計算傅里葉系數(shù)要先反演得到截距阻抗、梯度阻抗及曲率阻抗，因此提出利用疊前地震三參數(shù)貝葉斯反演方法，通過分方位地震數(shù)據(jù)分別進行截距阻抗、梯度阻抗及曲率阻抗的同時反演.對于一個方位，m個入射角和n個時間采樣點情況下，式(2)可重寫為矩陣表達式：

(26)

其中，上標T表示矩陣的轉(zhuǎn)置，符號diag表示對角矩陣.

考慮地震子波的影響，式(26)變?yōu)椋?/p>

(27)

式中，W(θi)表示角度子波矩陣.

式(27)可進一步簡化為矩陣表達式：

d=Gm，

(28)

其中：

假設(shè)地震噪聲服從高斯分布，待反演模型參數(shù)的先驗信息服從柯西分布，則后驗概率密度函數(shù)可表示為：

(29)

通過推導(dǎo)式(29)得到初始目標函數(shù)為：

(30)

其中，λQ表示柯西權(quán)重因子.由于地震數(shù)據(jù)缺乏低頻信息，引入低頻模型約束項能夠補充待反演參數(shù)的低頻信息，并提高反演的穩(wěn)定性及橫向連續(xù)性(印興耀等，2013；Zong et al.，2013；潘新朋等，2018)，因此最終的目標函數(shù)為：

(31)

式中，λA、λB和λC分別為截距項、梯度項及曲率項的低頻權(quán)重因子，AI0、BI0(φ)和CI0(φ)分別為與截距阻抗、梯度阻抗及曲率阻抗相關(guān)的低頻模型，P為積分矩陣.

在求得A、B(φ)和C(φ)之后，利用道積分即可求得截距阻抗、梯度阻抗及曲率阻抗.通過調(diào)試不同的角合成擴展方位彈性阻抗，對其進行傅里葉級數(shù)展開得到傅里葉系數(shù)，計算二階傅里葉系數(shù)與井中裂縫參數(shù)及DHSR的最大相關(guān)系數(shù)，分別確定最優(yōu)的角，進而實現(xiàn)裂縫參數(shù)及DHSR預(yù)測.

2 模型試算

利用某工區(qū)井數(shù)據(jù)開展模型測試，以驗證提出方法的有效性.將圖1中的井曲線和主頻為35 Hz的雷克子波褶積合成方位角道集，入射角范圍為5°～40°，以5°間隔，方位角分別為0°、45°、90°和135°.向合成的方位角道集中分別添加信噪比(Signal to noise ratio，S/N)為5和2的高斯噪聲以模擬觀測地震記錄.圖4展示了不同信噪比的合成方位角道集.圖5展示了不同信噪比情況下反演的歸一化的截距阻抗、梯度阻抗及曲率阻抗，圖6展示了不同信噪比情況下反演的歸一化的截距阻抗、梯度阻抗及曲率阻抗與真實值的相對誤差.可以看出截距阻抗反演結(jié)果的相對誤差最小，梯度阻抗和曲率阻抗反演結(jié)果的相對誤差相對較大.截距阻抗反演結(jié)果的總體相對誤差小于5%，而梯度阻抗和曲率阻抗反演結(jié)果的總體相對誤差小于10%，在局部區(qū)域有時接近20%，但總體而言，在信噪比為5和2時均能得到較好的截距阻抗、梯度阻抗及曲率阻抗反演結(jié)果.圖7展示了不同信噪比情況下反演的法向弱度、切向弱度及DHSR.圖8展示了不同信噪比情況下反演的法向弱度、切向弱度及DHSR與真實值的相對誤差.可以看出，隨著信噪比降低，法向弱度、切向弱度及DHSR反演結(jié)果的相對誤差均有所增加，但總體相對誤差均小于10%.分別計算不同信噪比情況下的反演結(jié)果與真實值的相關(guān)系數(shù)，信噪比為5情況下，反演的法向弱度、切向弱度及DHSR與真實值的相關(guān)系數(shù)分別為0.90、0.93和0.91；信噪比為2情況下，反演的法向弱度、切向弱度及DHSR與真實值的相關(guān)系數(shù)分別為0.83、0.84和0.84.因此，即使在信噪比為2情況下，利用提出的方法也能夠得到合理可靠的法向弱度、切向弱度及DHSR預(yù)測結(jié)果，驗證了該方法在預(yù)測裂縫參數(shù)及DHSR方面的有效性.

圖4 合成方位角道集(a) 信噪比為5； (b) 信噪比為2.Fig.4 Synthetic azimuthal angle gather(a) S/N=5； (b) S/N=2.

圖5 不同信噪比情況下的模型參數(shù)反演結(jié)果(a) 歸一化的截距阻抗； (b) 歸一化的梯度阻抗； (c) 歸一化的曲率阻抗.Fig.5 Inversion result of model parameters for different noise levels(a) Normalized intercept impedance； (b) Normalized gradient impedance； (c) Normalized curvature impedance.

圖6 不同信噪比情況下模型參數(shù)反演結(jié)果與真實值的相對誤差(a) 歸一化的截距阻抗； (b) 歸一化的梯度阻抗； (c) 歸一化的曲率阻抗.Fig.6 Relative error between inversion result of model parameters and the true value for different noise levels(a) Normalized intercept impedance； (b) Normalized gradient impedance； (c) Normalized curvature. impedance.

圖7 不同信噪比情況下反演的法向弱度、切向弱度及DHSRFig.7 Inverted normal weakness, tangential weakness, and DHSR for different noise levels

圖8 不同信噪比情況下反演的法向弱度、切向弱度及DHSR與真實值的相對誤差Fig.8 Relative error between inversion result of normal weakness, tangential weakness, DHSR and the true value for different noise levels

3 實際資料處理

研究區(qū)域位于四川盆地新場工區(qū)，該工區(qū)儲層屬超低孔、超低滲型致密砂巖儲層，構(gòu)造縫較為發(fā)育，多為高角度縫和立縫，構(gòu)造縫多數(shù)未充填，對產(chǎn)能貢獻極大，因此可將儲層等效為HTI介質(zhì).經(jīng)過處理后的疊前方位角道集首先被劃分為四個方位，具體方位部分角度疊加方案為22.5°(0°～45°)、67.5°(45°～90°)、112.5°(90°～135°)和157.5°(135°～180°).對每一個方位角道集又分別進行入射角部分角度疊加，得到三個平均入射角，分別為15°(10°～19°)、22°(17°～26°)和29°(24°～33°).利用提出的分方位疊前地震貝葉斯反演方法得到的歸一化截距阻抗、梯度阻抗及曲率阻抗反演結(jié)果如圖9所示.最終預(yù)測的法向弱度、切向弱度及DHSR結(jié)果如圖10所示.提取井旁道地震反演結(jié)果與測井曲線進行對比，如圖11所示，可以看出在井旁道的反演結(jié)果與井曲線之間盡管存在一定誤差，但總體趨勢與井曲線基本一致.分別計算反演的裂縫參數(shù)及DHSR與井曲線的相關(guān)系數(shù)，反演的法向弱度、切向弱度及DHSR與井中相應(yīng)曲線的相關(guān)系數(shù)分別為0.85、0.84、0.85.反演的法向弱度和切向弱度剖面在儲層位置處(黑色橢圓)均表現(xiàn)為相對高值，可靠地指示了儲層的裂縫發(fā)育情況；反演的DHSR剖面在儲層位置處也表現(xiàn)為相對高值，在壓裂時應(yīng)選取DHSR數(shù)值相對較小的區(qū)域.利用提出的方法預(yù)測得到的裂縫參數(shù)及DHSR剖面具有較好的橫向連續(xù)性，有助于裂縫發(fā)育區(qū)域及有利壓裂區(qū)域的橫向識別.

圖9 反演的(a)歸一化截距阻抗，(b)歸一化梯度阻抗及(c)歸一化曲率阻抗(黑線表示井的位置)Fig.9 Inverted (a) normalized intercept impedance, (b) normalized gradient impedance, and (c) normalized curvature impedance (the black line indicates the location of the well)

圖10 反演的(a)法向弱度，(b)切向弱度及(c)DHSRFig.10 Inverted (a) normal weakness, (b) tangential weakness, and (c) DHSR

圖11 井旁道反演結(jié)果與測井曲線對比Fig.11 Comparison between nearby well inversion results and logging curves

4 結(jié)論

裂縫參數(shù)及DHSR的準確估測對于地下裂縫預(yù)測及水力壓裂具有重要意義，本文提出了一種基于擴展方位彈性阻抗的裂縫參數(shù)及DHSR預(yù)測方法.首先推導(dǎo)了HTI介質(zhì)中的擴展方位彈性阻抗方程及傅里葉系數(shù)方程，傅里葉系數(shù)相關(guān)性分析表明，零階傅里葉系數(shù)與各向同性參數(shù)具有較好的相關(guān)性，但與裂縫參數(shù)及DHSR相關(guān)性較差；二階傅里葉系數(shù)與裂縫參數(shù)及DHSR具有較好的相關(guān)性.通過結(jié)合不同的角和分方位疊前地震貝葉斯反演得到截距阻抗、梯度阻抗及曲率阻抗，分析二階傅里葉系數(shù)與井中裂縫參數(shù)及DHSR的相關(guān)性以確定最優(yōu)的角，進而利用二階傅里葉系數(shù)實現(xiàn)裂縫參數(shù)及DHSR預(yù)測.模型測試表明，即使在信噪比為2情況下，利用提出的方法也能夠得到合理可靠的裂縫參數(shù)及DHSR預(yù)測結(jié)果.實際應(yīng)用表明，利用提出的方法預(yù)測的裂縫參數(shù)及DHSR與井中曲線較為符合，且反演剖面具有較好的橫向連續(xù)性，有助于裂縫發(fā)育區(qū)域及有利壓裂區(qū)域的橫向識別.