基于擴(kuò)展方位彈性阻抗的裂縫參數(shù)及地應(yīng)力參數(shù)預(yù)測(cè)方法

李林， 張廣智*

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

0 引言

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

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

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

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

1 方法原理

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

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

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

(1)

式(1)可重寫(xiě)為：

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

(2)

式中：

(3)

(4)

(5)

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

EEI(,

(9)

(10)

(11)

(12)

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

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

(13)

其中，LEEI(,

對(duì)式(13)進(jìn)行傅里葉級(jí)數(shù)展開(kāi)，得到：

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ù)可通過(guò)式(18)計(jì)算：

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

(18)

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

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

(19)

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

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

圖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)可得到各向異性梯度表達(dá)式：

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

(22)

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

(23)

(24)

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

(25)

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

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

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

(26)

其中，上標(biāo)T表示矩陣的轉(zhuǎn)置，符號(hào)diag表示對(duì)角矩陣.

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

(27)

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

式(27)可進(jìn)一步簡(jiǎn)化為矩陣表達(dá)式：

d=Gm，

(28)

其中：

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

(29)

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

(30)

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

(31)

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

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

2 模型試算

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

3 實(shí)際資料處理

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

圖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é)果與測(cè)井曲線對(duì)比Fig.11 Comparison between nearby well inversion results and logging curves

4 結(jié)論

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