基于扩展方位弹性阻抗的裂缝参数及地应力参数预测方法

李林， 张广智*

1 中国石油大学(华东)深层油气重点实验室， 青岛 266580 2 中国石油大学(华东)地球科学与技术学院， 青岛 266580

0 引言

随着能源需求的增加和常规油气资源的减少，非常规油气资源(页岩气、致密气和致密油储层)成为石油行业的关键勘探目标.非常规储层具有低孔低渗的特点，需要通过大规模的水力压裂以提高油气产量.地下裂缝能够有效改善储层的渗透率，地应力是控制水力压裂的重要因素(陈志刚等，2020；王生奥等，2021)，因此利用高品质“五维”地震数据开展裂缝参数及地应力参数地震反演对于地下裂缝预测及水力压裂具有重要意义(印兴耀等，2018a，b).

早期的地应力预测模型假设地层为各向同性，将各向同性地应力预测模型应用于各向异性地层会导致地应力预测不准确，因此各向异性地应力模型更加符合非常规储层的实际情况(赵小龙等，2020).许多学者考虑页岩具有各向异性的微观结构，提出了适用于垂直横向各向同性(Vertical Transverse Isotropic, VTI)介质的地应力预测模型(Thiercelin and Plumb，1994；Higgins et al.，2008；邓金根等，2013).张广智等(2015)提出利用页岩岩石物理等效模型进行地应力预测的方法.考虑到页岩储层中高角度裂缝较为发育的特征，Gray等(2012)推导了水平横向各向同性(Horizontal Transverse Isotropic, HTI)介质中的水平应力预测方程，并指出将水平应力差异比(Differential Horizontal Stress Ratio，DHSR)作为衡量储层是否能压裂成网的指示因子.Far等(2016)研究了两组正交的垂直裂缝嵌入VTI背景诱导的正交各向异性介质中的水平应力预测模型.马妮等(2017，2018)推导了单组垂直裂缝嵌入VTI背景诱导的正交各向异性介质中的水平应力预测方程.

弹性阻抗(Elastic Impedance, EI)在预测弹性参数方面有着广泛的应用前景(Connolly，1999；Whitcombe，2002；Wang et al.，2006；印兴耀等，2013；Zong et al.，2013).Martins(2006)将弹性阻抗扩展到各向异性介质中，随后许多学者研究了利用各向异性弹性阻抗提取各向异性参数的反演方法(陈怀震等，2014；Chen et al., 2018；Pan et al.，2020b；Li et al.，2020).Whitcombe等(2002)提出扩展弹性阻抗(Extended Elastic Impedance, EEI)的概念，通过扩展弹性阻抗可转换为其他弹性参数、物性参数或地质力学参数，进而用于岩性、流体识别及地应力预测(Hafez and Castagna，2016；Aleardi，2018；Sharifi et al.，2019；Sharifi and Mirzakhanian，2019).研究表明深度、压实趋势、厚度和岩性等因素会影响EEI分析中的相关性趋势(Ball et al.，2013；Thomas et al.，2013).Connolly(2019)指出弹性参数与扩展弹性阻抗的相关性取决于横纵波速度平方比，且由于反演存在的误差，最优的旋转角应通过对反演结果而不是测井数据进行EEI分析得到.Russell和Hedlin(2019)基于Biot-Gassmann多孔弹性理论提出了扩展多孔弹性阻抗的概念，并指出扩展多孔弹性阻抗与弹性参数的相关性取决于干岩石的横纵波速度平方比.

考虑单组垂直裂缝诱导的方位各向异性的影响，笔者将扩展弹性阻抗推广到HTI介质中，提出扩展方位弹性阻抗的概念.首先推导HTI介质中的扩展方位弹性阻抗方程及其傅里叶级数展开表达式；基于推导的傅里叶系数(Fourier Coefficient，FC)方程，分析傅里叶系数与裂缝参数及DHSR之间的关系，进而提出一种基于扩展方位弹性阻抗的裂缝参数及DHSR预测方法，最后通过模型测试和实际应用验证该方法在裂缝参数及DHSR预测方面的可靠性.

1 方法原理

1.1 HTI介质扩展方位弹性阻抗方程及傅里叶系数

HTI介质中的纵波反射系数方程为(潘新朋等，2018)：

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

(1)

式(1)可重写为：

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

(2)

式中：

(3)

(4)

(5)

Whitcombe等(2002)推导了扩展弹性阻抗方程，通过改变单一的旋转角即可由扩展弹性阻抗转化为其他弹性参数.类比各向同性介质中的扩展弹性阻抗推导过程，将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(,φ)表示尺度化的方位反射系数.

弹性阻抗的定义为(Connolly, 1999)：

(8)

结合式(7)、(8)，通过推导可得到HTI介质中的扩展方位弹性阻抗方程为：

EEI(,

(9)

(10)

(11)

(12)

对式(9)两边同时取对数，得到归一化的扩展方位弹性阻抗方程：

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

(13)

其中，LEEI(,

对式(13)进行傅里叶级数展开，得到：

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)表示与旋转角相关的傅里叶系数，Bani=g(δT-kδN)表示各向异性梯度.

傅里叶系数可通过式(18)计算：

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

(18)

1.2 傅里叶系数与裂缝参数及DHSR的关系

Gray等(2012)推导了杨氏模量、泊松比及法向柔度表征的DHSR表达式，通过推导可得到由纵、横波阻抗(IP和IS)及法向弱度表征的DHSR，其表达式为：

(19)

式(19)表明，DHSR仅与地层系数g和法向弱度δN有关.根据式(15)、(16)、(17)可知，零阶傅里叶系数与各向同性参数及裂缝参数均有关，二阶和四阶傅里叶系数主要与裂缝参数有关.通过变化角，四阶傅里叶系数只能转化为而不能转化为裂缝参数和DHSR，因此下面以某工区A井数据为例，分析零阶和二阶傅里叶系数与各向同性参数、裂缝参数及DHSR之间的关系.

图1为A井中的各向同性参数及裂缝参数，包括纵、横波速度、密度、法向及切向裂缝弱度.结合式(15)，利用井数据合成不同角情况下的零阶傅里叶系数如图2a所示，可以看出零阶傅里叶系数随着角的变化而变化.由于零阶傅里叶系数与各向同性参数和裂缝参数均相关，因此分别计算零阶傅里叶系数与各向同性参数、裂缝参数及DHSR的相关系数，如图2b所示，可以看出纵、横波速度、密度、法向弱度、切向弱度及DHSR与零阶傅里叶系数的相关系数随着角的变化而变化，纵、横波速度、密度、法向弱度、切向弱度及DHSR与零阶傅里叶系数的最大相关系数绝对值分别为0.96、0.99、0.83、0.72、0.74和0.74.可以看出法向弱度、切向弱度及DHSR与零阶傅里叶系数的相关系数远低于各向同性参数与零阶傅里叶系数的相关系数.李林等(2020)研究表明，裂缝参数对零阶傅里叶系数的贡献远小于各向同性参数的贡献，导致零阶傅里叶系数对裂缝参数的敏感性较低，因此零阶傅里叶系数仅可用于估测各向同性参数，而不能用于裂缝参数和DHSR估测.

图1 A井的各向同性参数及裂缝参数Fig.1 Isotropic parameters and fracture parameters of well A

图2 合成的(a)零阶傅里叶系数，(b)零阶傅里叶系数与各向同性参数、裂缝参数及DHSR的相关系数Fig.2 Synthetic (a) zeroth FC, and (b) correlation coefficients of the zeroth FC with isotropic parameters, fracture parameters and DHSR

图3 合成的(a)二阶傅里叶系数，(b)二阶傅里叶系数与裂缝参数和DHSR的相关系数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)

结合式(20)、(21)可得到各向异性梯度表达式：

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

(22)

在叠前地震反演中，通常假设一段时窗内的地层系数g为一常数，例如假设g=0.25(Whitcombe et al.，2002；Russell and Hedlin，2019).结合各向异性梯度与裂缝参数的关系表达式，通过推导可得到：

(23)

(24)

结合式(19)、(23)可得到HTI介质中的DHSR简化表达式：

(25)

因此，当假设地层系数g为0.25时，利用式(23)、(24)、(25)即可由二阶傅里叶系数分别估测HTI介质中的裂缝参数及DHSR.Connolly(2019)指出弹性参数与扩展弹性阻抗的相关性取决于地层系数g.本文研究表明，二阶傅里叶系数与裂缝参数及DHSR的相关性也与地层系数g相关.由于地层系数g是随着时间或深度变化的，因此在实际应用中应通过寻找二阶傅里叶系数与裂缝参数及DHSR的最大相关系数，分别确定最优的角进行裂缝参数及DHSR预测.

1.3 分方位叠前地震贝叶斯反演

由于计算傅里叶系数要先反演得到截距阻抗、梯度阻抗及曲率阻抗，因此提出利用叠前地震三参数贝叶斯反演方法，通过分方位地震数据分别进行截距阻抗、梯度阻抗及曲率阻抗的同时反演.对于一个方位，m个入射角和n个时间采样点情况下，式(2)可重写为矩阵表达式：

(26)

其中，上标T表示矩阵的转置，符号diag表示对角矩阵.

考虑地震子波的影响，式(26)变为：

(27)

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

式(27)可进一步简化为矩阵表达式：

d=Gm，

(28)

其中：

假设地震噪声服从高斯分布，待反演模型参数的先验信息服从柯西分布，则后验概率密度函数可表示为：

(29)

通过推导式(29)得到初始目标函数为：

(30)

其中，λQ表示柯西权重因子.由于地震数据缺乏低频信息，引入低频模型约束项能够补充待反演参数的低频信息，并提高反演的稳定性及横向连续性(印兴耀等，2013；Zong et al.，2013；潘新朋等，2018)，因此最终的目标函数为：

(31)

式中，λA、λB和λC分别为截距项、梯度项及曲率项的低频权重因子，AI0、BI0(φ)和CI0(φ)分别为与截距阻抗、梯度阻抗及曲率阻抗相关的低频模型，P为积分矩阵.

在求得A、B(φ)和C(φ)之后，利用道积分即可求得截距阻抗、梯度阻抗及曲率阻抗.通过调试不同的角合成扩展方位弹性阻抗，对其进行傅里叶级数展开得到傅里叶系数，计算二阶傅里叶系数与井中裂缝参数及DHSR的最大相关系数，分别确定最优的角，进而实现裂缝参数及DHSR预测.

2 模型试算

利用某工区井数据开展模型测试，以验证提出方法的有效性.将图1中的井曲线和主频为35 Hz的雷克子波褶积合成方位角道集，入射角范围为5°～40°，以5°间隔，方位角分别为0°、45°、90°和135°.向合成的方位角道集中分别添加信噪比(Signal to noise ratio，S/N)为5和2的高斯噪声以模拟观测地震记录.图4展示了不同信噪比的合成方位角道集.图5展示了不同信噪比情况下反演的归一化的截距阻抗、梯度阻抗及曲率阻抗，图6展示了不同信噪比情况下反演的归一化的截距阻抗、梯度阻抗及曲率阻抗与真实值的相对误差.可以看出截距阻抗反演结果的相对误差最小，梯度阻抗和曲率阻抗反演结果的相对误差相对较大.截距阻抗反演结果的总体相对误差小于5%，而梯度阻抗和曲率阻抗反演结果的总体相对误差小于10%，在局部区域有时接近20%，但总体而言，在信噪比为5和2时均能得到较好的截距阻抗、梯度阻抗及曲率阻抗反演结果.图7展示了不同信噪比情况下反演的法向弱度、切向弱度及DHSR.图8展示了不同信噪比情况下反演的法向弱度、切向弱度及DHSR与真实值的相对误差.可以看出，随着信噪比降低，法向弱度、切向弱度及DHSR反演结果的相对误差均有所增加，但总体相对误差均小于10%.分别计算不同信噪比情况下的反演结果与真实值的相关系数，信噪比为5情况下，反演的法向弱度、切向弱度及DHSR与真实值的相关系数分别为0.90、0.93和0.91；信噪比为2情况下，反演的法向弱度、切向弱度及DHSR与真实值的相关系数分别为0.83、0.84和0.84.因此，即使在信噪比为2情况下，利用提出的方法也能够得到合理可靠的法向弱度、切向弱度及DHSR预测结果，验证了该方法在预测裂缝参数及DHSR方面的有效性.

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

图5 不同信噪比情况下的模型参数反演结果(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 不同信噪比情况下模型参数反演结果与真实值的相对误差(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 实际资料处理

研究区域位于四川盆地新场工区，该工区储层属超低孔、超低渗型致密砂岩储层，构造缝较为发育，多为高角度缝和立缝，构造缝多数未充填，对产能贡献极大，因此可将储层等效为HTI介质.经过处理后的叠前方位角道集首先被划分为四个方位，具体方位部分角度叠加方案为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°).利用提出的分方位叠前地震贝叶斯反演方法得到的归一化截距阻抗、梯度阻抗及曲率阻抗反演结果如图9所示.最终预测的法向弱度、切向弱度及DHSR结果如图10所示.提取井旁道地震反演结果与测井曲线进行对比，如图11所示，可以看出在井旁道的反演结果与井曲线之间尽管存在一定误差，但总体趋势与井曲线基本一致.分别计算反演的裂缝参数及DHSR与井曲线的相关系数，反演的法向弱度、切向弱度及DHSR与井中相应曲线的相关系数分别为0.85、0.84、0.85.反演的法向弱度和切向弱度剖面在储层位置处(黑色椭圆)均表现为相对高值，可靠地指示了储层的裂缝发育情况；反演的DHSR剖面在储层位置处也表现为相对高值，在压裂时应选取DHSR数值相对较小的区域.利用提出的方法预测得到的裂缝参数及DHSR剖面具有较好的横向连续性，有助于裂缝发育区域及有利压裂区域的横向识别.

图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 井旁道反演结果与测井曲线对比Fig.11 Comparison between nearby well inversion results and logging curves

4 结论

裂缝参数及DHSR的准确估测对于地下裂缝预测及水力压裂具有重要意义，本文提出了一种基于扩展方位弹性阻抗的裂缝参数及DHSR预测方法.首先推导了HTI介质中的扩展方位弹性阻抗方程及傅里叶系数方程，傅里叶系数相关性分析表明，零阶傅里叶系数与各向同性参数具有较好的相关性，但与裂缝参数及DHSR相关性较差；二阶傅里叶系数与裂缝参数及DHSR具有较好的相关性.通过结合不同的角和分方位叠前地震贝叶斯反演得到截距阻抗、梯度阻抗及曲率阻抗，分析二阶傅里叶系数与井中裂缝参数及DHSR的相关性以确定最优的角，进而利用二阶傅里叶系数实现裂缝参数及DHSR预测.模型测试表明，即使在信噪比为2情况下，利用提出的方法也能够得到合理可靠的裂缝参数及DHSR预测结果.实际应用表明，利用提出的方法预测的裂缝参数及DHSR与井中曲线较为符合，且反演剖面具有较好的横向连续性，有助于裂缝发育区域及有利压裂区域的横向识别.