正如“地震岩石物理驱动的裂缝预测技术研究现状与进展(Ⅰ)”的引言所述, 地震预测技术的探测深度大、范围广等特征是地质、岩心和测井技术等非地震手段所不具备的, 而且地震P波振幅^[1]、阻抗^[2-4]、衰减^[5-6]等属性随入射角和方位角变化(五维地震)技术较其它地震裂缝预测技术具有更高的分辨率和抗噪性, 可获得更多的储层信息且纵波数据易获得。所以在第(Ⅰ)部分^[7]综述了裂缝储层等效各向异性岩石物理理论研究现状与进展的基础上, 第二部分(Ⅱ)主要阐述地震P波属性随入射角和方位角变化(五维地震)的裂缝预测技术研究现状与进展。

五维地震数据较常规叠后数据增加了方位和偏移距维度, 蕴含着更加丰富的储层裂缝等各向异性信息^[8-9], 所以基于五维地震数据的裂缝信息解耦与提取得到了越来越多的关注^[10-14]。根据裂缝指示因子的差异, 目前基于五维地震数据的裂缝预测技术可大致分为基于椭圆分析的裂缝预测方法和基于岩石物理理论的裂缝参数反演方法, 前者利用拟合椭圆参数指示裂缝分布, 后者利用反演得到的岩石物理参数及其组合实现裂缝参数的空间描述。所以我们将从两方面阐述基于五维地震数据随方位角和入射角变化的裂缝预测理论与方法的研究现状。

1 基于椭圆分析的裂缝预测方法

基于椭圆分析裂缝预测方法的科学实质是裂缝介质诱导的地球物理响应在方位角域呈周期性变化, 且这种变化的极值点与裂缝对称轴构成映射关系。目前研究中, 可用于椭圆分析的参数主要包括运动学参数和动力学参数, 运动学参数包括动校正速度^[15-17]和各向异性剩余时差^[18-19]等, 动力学参数包括AVO梯度和杨氏模量等。前者主要利用地层裂缝诱导的地震波旅行时方位差异实现裂缝参数预测, 目前在HTI型^[17-19], OA(orthorhombic medium)型^[15]、TOA(titled orthorhombic medium)型^{[16, 20]}裂缝储层中均有研究和应用; 后者的理论基础是地震波振幅的方位差异, 目前主要应用于HTI型裂缝储层。所以我们将以各向同性背景中包含一组垂直平行裂缝介质(parallel vertical fractures in isotropic background, VFI)为例, 详细介绍基于AVO梯度椭圆拟合^[21]和方位杨氏模量椭圆拟合^[22]的裂缝预测技术。

1.1 基于AVO梯度椭圆拟合的裂缝预测方法

RÜGER^[23]给出了HTI介质反射系数近似公式, 可将其简写为:

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right)=R_{\mathrm{PP}}^{\mathrm{iso}}+R_{\mathrm{PP}}^{\mathrm{ani}}=P+ \\ G(\varphi) \sin ^{2} \theta+C(\varphi) \sin ^{2} \theta \tan ^{2} \theta \end{gathered} $

(1)

式中: φ=φ_obs－φ_sym; φ_sym和φ_obs分别为裂缝法向对称轴的方位角和观测方位角; R_PP(θ, φ_obs), R_PP^iso和R_PP^ani分别表示入射角为θ、观测方位角为φ_obs时的反射系数, 反射系数各向同性部分和各向异性部分; P, G(φ)和C(φ)分别为反射系数的截距, AVO梯度和曲率。当入射角小于30°时, 可舍掉sin²θtan²θ项, 得到估计AVO梯度的基础方程, 即:

$ R_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right)=R_{\mathrm{PP}}^{\text {iso }}+R_{\mathrm{PP}}^{\text {ani }}=P+G(\varphi) \sin ^{2} \theta $

(2)

根据公式(2), 可逐方位对五维地震数据进行AVO分析^[24], 将估计的多方位AVO梯度进行椭圆拟合^[25], 拟合椭圆的主轴指示裂缝对称轴方位(后文简称裂缝方位), 椭圆率指示裂缝密度分布特征。但逐方位AVO分析获得的自激自收反射系数P很难保证相等, 影响AVO梯度G(φ)的提取精度, 所以可联合五维地震数据构建反演方程, 以保证自激自收反射系数P在方位角域的常值特征^[26]。此外, 如图 1所示, 对于不同模型参数(表 1), 裂缝走向和裂缝对称轴方向的AVO梯度相对大小不同, 即AVO梯度拟合椭圆的长轴对应裂缝走向还是裂缝对称轴方向不明确, 通常将此现象称为90°模糊性。导致90°模糊性的根本原因为AVO梯度可表达为G(φ)=A+Bcos[2(φ－φ_sym)], 呈现周期为π的余弦变化规律, 式中: A, B均为介质弹性参数和各向异性参数的函数; φ_sym为裂缝对称轴。当参数B为正数时, 裂缝对称轴方向的AVO梯度取最大值, 即G(φ_sym)=A+|B|; 当参数B为负值时, 裂缝对称轴方向的AVO梯度取最小值, 即G(φ_sym)=A－|B|。由于B可为正数, 也可为负数, 取决于岩石弹性参数和各向异性参数, 而且周期为π的余弦曲线最大函数值和最小函数值的相位差为90°, 所以在岩石性质未知的情况下, 利用AVO梯度椭圆拟合的裂缝方位预测值存在90°模糊性, 需要在先验方位约束下开展裂缝方位预测。图 2是利用所述方法在A工区目的层裂缝参数预测结果, 与测井解释结果相符, 说明了方法的合理性。

1.2 基于方位杨氏模量椭圆拟合的裂缝预测

VFI介质中的纵、横波相速度v_P、v_S是地震波传播方位的函数, 同时杨氏模量与地震波速度存在函数关系, 即E=[v_S²(3v_P²－4v_S²)]/(v_P²－v_S²), 所以杨氏模量也是方位角的周期函数。SAYERS^[27]研究了通过动、静态测量两种方式获得的VFI介质对称轴方向杨氏模量E′和各向同性面上杨氏模量E之间的关系(图 3)。由图可见, 基本所有样点满足E>E′, 即VFI介质中裂缝走向的杨氏模量大于裂缝对称轴方向的杨氏模量。所以, 与AVO梯度一样, 方位杨氏模量也可以拟合为一个椭圆, 且椭圆率指示裂缝密度, 但此时不需要裂缝先验方位信息作为约束, 椭圆长轴直接指示裂缝走向。

宗兆云等^[28]推导了杨氏模量表示的反射系数方程和相应的弹性阻抗方程^[29]:

$ \mathrm{EI}=\mathrm{EI}_{0} \cdot\left(\frac{E}{E_{0}}\right)^{a_{E}(\theta)} \cdot\left(\frac{\sigma}{\sigma_{0}}\right)^{b_{\sigma}(\theta)} \cdot\left(\frac{\rho}{\rho_{0}}\right)^{c_{\rho}(\theta)} $

(3a)

其中,

$ a_{E}(\theta)=\frac{1}{2} \sec ^{2} \theta-4 \gamma^{-1} \sin ^{2} \theta $

(3b)

$ b_{\sigma}(\theta)=\frac{1}{2} \sec ^{2} \theta \frac{(2-3 \gamma)(2-\gamma)^{2}}{4-3 \gamma}+4 \gamma^{-1} \sin ^{2} \theta \frac{\gamma-2}{3 \gamma-4} $

(3c)

$ c_{\rho}(\theta)=1-\frac{1}{2} \sec ^{2} \theta $

(3d)

式中: γ=β²/α²; α和β分别为纵、横波速度; EI为弹性阻抗; E、σ和ρ分别代表杨氏模量、泊松比和密度; 下角标0代表归一化常数, 一般为各弹性参数的平均值。基于公式(3), 采用岩石物理驱动的阻尼最小二乘算法^[30]或贝叶斯反演算法^[31]等即可求得每个方位的杨氏模量, 然后利用方位杨氏模量拟合椭圆参数实现储层裂缝预测。图 4为基于方位杨氏模量椭圆拟合的A工区裂缝预测结果, 与测井解释结果基本吻合。但和图 2对比发现, 基于AVO梯度与基于方位杨氏模量椭圆拟合的裂缝预测结果不完全相同, 可能是由于两种方法的抗噪性不同和AVO梯度受子波影响较大等原因导致的。

2 基于等效各向异性岩石物理理论的裂缝预测方法

地下储层物性特征复杂多变, 根据储层空间的连通性, 可分为孤立孔缝型储层、连通孔缝弹性储层及连通孔缝衰减型储层。针对不同类型储层, 应采用不同的岩石物理模型进行表征(如前文^[7]所述), 即其模型参数化方法不一样。所以, 本节将以等效各向异性岩石物理理论为基础, 分别介绍这3种储层的模型参数化方法及裂缝参数预测的研究进展。

2.1 基于孤立孔缝岩石物理理论的裂缝预测方法 2.1.1 基于各向异性梯度的裂缝预测

在Schoenberg理论指导下, BAKULIN^[32-33]指出各向异性梯度是很好的裂缝密度指示因子, 为实现稳定的裂缝方位和各向异性梯度反演, 李春鹏等^[13]根据DOWNTON等^[34]的分析, 将公式(2)改写为:

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi_{\text {obs }}\right)=R_{\mathrm{PP}}^{\mathrm{iso}}+R_{\mathrm{PP}}^{\mathrm{ani}}=P+G\left(\varphi_{\text {obs }}-\varphi_{\mathrm{sym}}\right) \sin ^{2} \theta \\ =P+\left[G_{\mathrm{iso}}+G_{\text {ani }} \cos ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)\right] \sin ^{2} \theta \\ =P+\left(B+\frac{1}{2} C \cos 2 \varphi_{\text {obs }}+\frac{1}{2} D \sin 2 \varphi_{\text {obs }}\right) \sin ^{2} \theta \end{gathered} $

(4a)

其中,

$ B =G_{\text {iso }}+\frac{1}{2} G_{\text {ani }} $

(4b)

$ C =G_{\text {ani }} \cos 2 \varphi_{\text {sym }} $

(4c)

$ D =G_{\text {ani }} \sin 2 \varphi_{\mathrm{sym}} $

(4d)

式中: G_iso和G_ani分别称为各向同性梯度和各向异性梯度。根据公式(4), 在贝叶斯理论框架下, 利用多方位地震数据便可实现P, B, C和D的反演^[13]。然后, 可估计各向异性梯度绝对值和裂缝方位, 即:

$ \left|G_{\text {ani }}\right|=\sqrt{C^{2} \times D^{2}} $

(4e)

$ \varphi_{\mathrm{sym}}=\frac{1}{2} \arctan \frac{D}{C}+n \frac{{\rm{ \mathit{ π} }}}{2}, n=0, 1, \cdots $

(4f)

从公式(4f)可发现, 裂缝方位估计值存在90°模糊性, 这是由于公式(4e)中的各向异性梯度符号难以确定而导致的, 为此需要裂缝方位先验信息约束以消除90°模糊性。图 5为该方法的模型测试结果。模型信噪比为4∶1, 可以看出各向异性梯度分布能很好地预测裂缝密度分布, 裂缝对称轴方位预测结果与模型也具有较高的吻合度。XIE等^[17]又重新改写了公式(4a), 利用最小二乘算法实现了裂缝方位和各向异性梯度的估计, 但依旧存在90°模糊性问题。陈怀震等^[35]基于弹性阻抗反演思想对该方法作了进一步发展, 得到了更加稳定的结果。

由于AVAZ反演(P-wave amplitude versus angle and azimuth inversion)的抗噪性比弹性阻抗反演低^[36-37], 所以基于公式(4)的裂缝参数直接反演很不稳定。曲寿利等^[2-3]提出了阻抗随方位角变化(impedance versus azimuth, IPVA)的裂缝预测技术以提高稳定性。他们将阻抗写为观测方位角的余弦函数, 即:

$ \operatorname{IP}(\varphi)=A_{\mathrm{IP}}+B_{\mathrm{IP}} \cos \left[2\left(\varphi-\varphi_{\mathrm{sym}}\right)\right] $

(5)

式中: IP(φ)为观测方位为φ的波阻抗, A_IP和B_IP为介质弹性参数和各向异性参数的函数, 且B_IP是裂缝密度的指示因子。在获得不同方位阻抗数据的基础上, 利用公式(5)进行拟合, 便可获得A_IP, B_IP和φ_sym的估计值, 从而实现裂缝特征的刻画。季玉新等^[4]将构造正反演裂缝预测技术与IPVA技术进行对比, 发现二者预测结果基本吻合, 进一步说明IPVA技术的有效性。但该技术依旧存在90°模糊性, 而且裂缝诱导的各向异性比较弱, 裂缝信息在阻抗数据中的比重较小, 降低了裂缝参数预测的稳定性。为了弱化90°模糊性, 进一步提高预测合理性, MA等^[38-39]基于傅里叶级数展开, 将裂缝方位先验信息作为约束条件, 构建一个新的VFI介质方位弹性阻抗方程:

$ \begin{gathered} \ln \frac{\mathrm{EI}\left(\varphi_{\mathrm{obs}}, \theta\right)}{\mathrm{EI}_{0}}=r_{\mathrm{EI}_{0}}(\theta)+\operatorname{sign}\left\{r _ { \mathrm { EI } _ { 2 } } ( \theta ) \operatorname { cos } \left[2 \left(\varphi_{\text {well }}-\right.\right.\right. \\ \left.\left.\left.\varphi_{\mathrm{sym}}\right)\right]\right\} \times\left|r_{\mathrm{EI}_{2}}(\theta)\right| \cos \left[2\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)\right] \end{gathered} $

(6)

式中: sign[·]是取符号算子; |·|代表取绝对值; φ_well代表来自地质、测井数据的裂缝对称轴方位先验信息; r_EIn是方位弹性阻抗的n阶傅里叶系数, 可用离散傅里叶变换进行计算; EI₀为平均波阻抗, 用于方位弹性阻抗的标准化; EI(φ_obs, θ)代表观测方位为φ_obs、入射角为θ的方位弹性阻抗, 可以利用约束稀疏脉冲反演^[40]、贝叶斯反演^[36]等方法从地震数据中获得。方程(6)的系数矩阵存在待反演参数, 所以MA等^[38-39]采用迭代理论完成求解, 获得无90°模糊性的裂缝方位估计。利用二阶傅里叶系数r_EI2对方位弹性阻抗的各向异性梯度进行估计G_ani^AEI≈r_EI2/sin²(θ)。需要注意的是, 方位弹性阻抗各向异性梯度G_ani^AEI=g[δ_T－(1－2γ)δ_N]是层属性, 也是裂缝密度指示因子。

图 6为B工区基于傅里叶级数的VFI介质裂缝参数预测结果。图 6a为方位弹性阻抗各向异性梯度层位切片, 黑线代表裂缝方位, 预测结果与该地区的地质构造相契合。图 6b为裂缝方位预测结果玫瑰统计, 与图 7所示的测井解释结果基本一致。MA等^[41]基于傅里叶级数展开, 也实现了VTI(vertically transverse isotropy)背景包含一组垂直定向裂缝诱导的OA介质(orthorhombic anisotropy)的裂缝参数预测, 并取得较好的应用效果。

通过地震属性提取也可以直接获得各向异性梯度和裂缝对称轴方位^[42]。以VFI介质为例, 陈怀震^[42]首先根据公式(4a)对3个观测方位φ₁, φ₂和φ₃的地震数据进行AVO分析, 提取截距P₁, P₂, P₃和梯度属性G₁, G₂, G₃, 然后估计裂缝方位φ_sym, 则有:

$ \begin{aligned} \varphi_{\mathrm{sym}}=& \frac{1}{2} \arctan \left\{\left[\cos \left(2 \varphi_{3}\right)-\cos \left(2 \varphi_{2}\right)-\right.\right.\\ &\left.D \cos \left(2 \varphi_{2}\right)+D \cos \left(2 \varphi_{1}\right)\right] /\left[D \sin \left(2 \varphi_{2}\right)-\right.\\ &\left.\left.D \sin \left(2 \varphi_{1}\right)-\sin \left(2 \varphi_{3}\right)+\sin \left(2 \varphi_{2}\right)\right]\right\} \end{aligned} $

(7a)

或

$ \begin{aligned} \varphi_{\mathrm{sym}}&=\frac{1}{2} \arctan \left\{\left[\cos \left(2 \varphi_{3}\right)-\cos \left(2 \varphi_{2}\right)-\right.\right. \\ &\left.D \cos \left(2 \varphi_{2}\right)+D \cos \left(2 \varphi_{1}\right)\right] /\left[D \sin \left(2 \varphi_{2}\right)-\right. \\ &\left.\left.D \sin \left(2 \varphi_{1}\right)-\sin \left(2 \varphi_{3}\right)+\sin \left(2 \varphi_{2}\right)\right]\right\}+90^{\circ} \end{aligned} $

(7b)

其中,

$ D=\frac{G_{3}-G_{2}}{G_{2}-G_{1}}=\frac{\cos ^{2}\left(\varphi_{3}-\varphi_{\mathrm{sym}}\right)-\cos ^{2}\left(\varphi_{2}-\varphi_{\mathrm{sym}}\right)}{\cos ^{2}\left(\varphi_{2}-\varphi_{\mathrm{sym}}\right)-\cos ^{2}\left(\varphi_{1}-\varphi_{\mathrm{sym}}\right)} $

(7c)

将估计值φ_sym带入G_ani=(G₃－G₂)/[cos²(φ₃－φ_sym)－cos²(φ₂－φ_sym)]完成各向异性梯度G_ani的估计。从公式(7)中发现, 裂缝方位估计也存在90°模糊性^[43], 需要裂缝方位先验信息作为约束。MESDAG^[44-45]指出, 在多数情况下, 各向异性梯度为正值, 将其作为约束, 也可弱化90°模糊性(注: 文献^[44]原文有误, 后与作者讨论, 确认大多数情况下各向异性梯度为正值)。对于OA型储层, 也可采用相似的方法实现各向异性梯度和裂缝方位的估计^[42]。但该方法提取的各向异性梯度属性分辨率较低, 且抗噪性较差, 只能作为勘探前期裂缝发育区的圈定手段, 或为后期裂缝参数的稳定预测提供低频约束。

2.1.2 基于岩石物理参数的裂缝预测

1) 各向同性背景中包含一组垂直定向裂缝(VFI)假设。从公式(6)可以看出, 当孤立裂缝充填流体时, VFI介质裂缝切向弱度δ_T仅与裂缝密度有关, 与裂缝充填流体类型无关, 所以可利用切向弱度直接计算裂缝密度^[32]:

$ e=\frac{3\left(3-2 \gamma_{\mathrm{b}}\right)}{16} \delta_{\mathrm{T}} $

(8)

印兴耀等^{[37, 46-56]}已经研究了多种流体识别方法, 针对孤立定向裂缝介质, SCHOENBERG等^[57]提出了一种有效识别裂缝流体类型的指示因子, 即:

$ \frac{Z_{\mathrm{N}}}{Z_{\mathrm{T}}}=\gamma_{\mathrm{b}} \frac{\delta_{\mathrm{N}}\left(1-\delta_{\mathrm{T}}\right)}{\delta_{\mathrm{T}}\left(1-\delta_{\mathrm{N}}\right)} $

(9)

式中: Z_N和Z_T分别为裂缝法向和切向柔度参数。可以看出, 裂缝密度e和流体指示因子Z_N/Z_T与裂缝弱度参数δ_N, δ_T和背景横、纵波速度比的平方γ_b有关, 所以只要反演出裂缝弱度参数即可实现裂缝密度预测及流体识别。

CHEN等^{[30, 58-59]}将VFI介质反射系数方程舍去密度项后, 改写为:

$ \begin{gathered} R_{\mathrm{PP}}(\theta, \varphi)=\sec ^{2} \theta R_{\mathrm{P}}-8 \gamma_{\mathrm{b}} \sin ^{2} \theta R_{\mathrm{S}}-\left(\gamma_{\mathrm{b}} \cos ^{2} \varphi\right. \\ \left.\sin ^{2} \theta\right)\left(1-2 \gamma_{\mathrm{b}}\right) \Delta_{\delta \mathrm{N}}+\left(\gamma_{\mathrm{b}} \cos ^{2} \varphi \sin ^{2} \theta\right) \Delta_{\delta_{\mathrm{T}}} \end{gathered} $

(10)

式中: φ=φ_obs－φ_sym; R_P=(I_P2－I_P1)/(I_P2+I_P1)=1/2(ΔI_P)/(I_P), R_S=(I_S2－I_S1)/(I_S2+I_S1)=1/2(ΔI_S/I_S), I_P, I_S为背景纵、横波阻抗, Δ代表界面上、下层介质属性差异, 参数上方的横线代表上下界面属性平均; Δ_δT=δ_T2－δ_T1, Δ_δN=δ_N2－δ_N1, δ_T和δ_N分别代表裂缝的切向和法向弱度参数。该式中只包含4个待求参数, 一定程度提高了反演稳定性。在获得裂缝方位φ_sym的基础上, 通过岩石物理约束的阻尼最小二乘算法执行反演, 最后利用道积分思想完成纵、横波阻抗及弱度参数的叠前地震预测^[30]:

$ \left[\begin{array}{c} \log \left(I_{\mathrm{P}}\right) \\ \log \left(I_{\mathrm{S}}\right) \\ \delta_{\mathrm{N}} \\ \delta_{\mathrm{T}} \end{array}\right]=\int\left(\begin{array}{c} \frac{\Delta I_{\mathrm{P}}}{I_{\mathrm{P}}} \\ \frac{\Delta I_{\mathrm{S}}}{I_{\mathrm{S}}} \\ \Delta_{\delta \mathrm{N}} \\ \Delta_{\delta_{\mathrm{T}}} \end{array}\right) \mathrm{d} t $

(11)

公式(11)中的黑体变量代表与时间相关的函数, 即地震道记录。在得到γ_b=(β/α)²=(I_S/I_P)²后, 便可实现裂缝密度预测和流体识别。图 8为C工区裂缝流体因子预测剖面, 低值区代表饱含水层, 高值区代表高含气层, 中间值从小到大依次代表油水层和气水层。预测剖面很好地勾勒出含气区域的边界, 与测井解释结果一致。陈怀震等^[60-61]基于弹性阻抗反演思想对该方法作了进一步发展, 以获得更加稳定的结果。

VFI介质裂缝诱导的弱各向异性信息在地震数据中的比重远远小于各向同性信息, 很容易被各向同性信息掩盖, 所以CHEN等^[62]根据SHAW等^[63-64]提出的散射系数与反射系数的关系推导了VFI介质反射系数方程, 即

$ \begin{gathered} R_{\mathrm{PP}}(\theta, \varphi)=\frac{1}{2 \cos ^{2} \theta} \frac{\Delta \alpha}{\alpha}-4 \gamma_{b} \sin ^{2} \theta \frac{\Delta \beta}{\beta}+ \\ \left(\frac{1}{2}-2 \gamma_{b} \sin ^{2} \theta\right) \frac{\Delta \rho}{\rho}+a(\theta, \varphi) \Delta_{\delta_{\mathrm{N}}}+b(\theta, \varphi) \Delta_{\delta_{\mathrm{T}}} \end{gathered} $

(12)

根据振幅方位差异思想, 得到不同观测方位反射系数差异表达:

$ \begin{gathered} \Delta R(\theta, \varphi)=\left[a\left(\theta, \varphi_{2}\right)-a\left(\theta, \varphi_{1}\right)\right] \Delta_{\delta \mathrm{N}}+ \\ {\left[b\left(\theta, \varphi_{2}\right)-b\left(\theta, \varphi_{1}\right)\right] \Delta_{\delta_{\mathrm{T}}}} \end{gathered} $

(13a)

其中,

$ a(\theta, \varphi)=-\frac{1}{4 \cos ^{2} \theta}\left[2\left(\sin ^{2} \theta \sin ^{2} \varphi+\cos ^{2} \theta\right) \gamma_{\mathrm{b}}-1\right]^{2} $

(13b)

$ b(\theta, \varphi)=-\gamma_{\mathrm{b}}\left(\sin ^{2} \theta \tan ^{2} \theta \sin ^{2} \varphi \cos ^{2} \varphi-\sin ^{2} \theta \cos ^{2} \varphi\right) $

(13c)

公式(13)实现了裂缝各向异性信息与各向同性背景信息的解耦, 可提高裂缝弱度参数预测稳定性。同时, 法向弱度与切向弱度参数的相关性可能影响二者反演的稳定性, CHEN等^[62]通过岩石物理分析, 获得二者的线性关系δ_N=ξδ_T+Γ, 代入公式(13), 消除待反演参数的相关性, 最后利用贝叶斯反演方法获得裂缝弱度参数, 实现裂缝储层描述。PAN等^[65]基于类似的思想, 采用马尔科夫-蒙特卡洛反演算法(MCMC)进行裂缝储层描述。为了消除待反演参数的空间相关性, PAN等^[66]也采用多尺度频率域反演方法, 将地震数据变换到傅里叶频率域, 再反演裂缝弱度参数。

LI等^[67]将方位弹性阻抗方程写为傅里叶级数形式, 即

$ \begin{gathered} L_{\mathrm{EI}}\left(\varphi_{\mathrm{obs}}, \theta\right)=a_{0}(\theta)+a_{2}(\theta) \cos \left(2 \varphi_{\mathrm{obs}}\right)+b_{2}(\theta) \cdot \\ \sin \left(2 \varphi_{\mathrm{obs}}\right)+a_{4}(\theta) \cos \left(4 \varphi_{\mathrm{obs}}\right)+b_{4}(\theta) \sin \left(4 \varphi_{\mathrm{obs}}\right) \end{gathered} $

(14)

其中,

$ \begin{gathered} a_{0}(\theta)=A(\theta) \ln \frac{\alpha}{\alpha_{0}}+B(\theta) \ln \frac{\beta}{\beta_{0}}+C(\theta) \ln \frac{\rho}{\rho_{0}}+ \\ F(\theta) \delta_{\mathrm{N}}+G(\theta) \delta_{\mathrm{T}} \end{gathered} $

$ a_{2}(\theta)=H(\theta) \cos \left(2 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{N}}+I(\theta) \cos \left(2 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{T}} $

$ b_{2}(\theta)=H(\theta) \sin \left(2 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{N}}+I(\theta) \sin \left(2 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{T}} $

$ a_{4}(\theta)=J(\theta) \cos \left(4 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{N}}+K(\theta) \cos \left(4 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{T}} $

$ b_{4}(\theta)=J(\theta) \sin \left(4 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{N}}+K(\theta) \sin \left(4 \varphi_{\mathrm{sym}}\right) \delta_{\mathrm{T}} $

A(θ), B(θ), C(θ), F(θ), G(θ), H(θ), I(θ), J(θ)和K(θ)均是横、纵波速度比的平方γ_b和入射角θ的函数。公式(14)同样解耦了裂缝信息与背景信息, 提高了裂缝弱度参数反演稳定性。此时, 裂缝方位被估计为:

$ \varphi_{\mathrm{sym}}=\frac{1}{2} \arctan \left[\frac{b_{2}(\theta)}{a_{2}(\theta)}\right] $

(15)

但仍存在90°模糊性。LI等^[67]采用各向异性岩石物理建模方法减弱该效应。类比于方位弹性阻抗, MESDAG等^[44]通过重组反射系数方程, 提出等效弹性参数, 并将其表达为傅里叶级数形式:

$ \begin{gathered} A^{\prime}=b_{0}+b_{1} \cos \left[2\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)\right]+b_{2} \cos \\ {\left[4\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)\right]} \end{gathered} $

(16)

式中: A′代表等效弹性参数; b_i代表等效弹性参数第i阶傅里叶系数。对于不同弹性参数, b_i具体表达不同。MESDAG等^[45]分析发现, 公式(16)中二阶傅里叶系数一般为正值, 所以基于公式(16)可实现无90°模糊性的裂缝方位估计及裂缝密度预测。

以上方法通常要求方位分扇叠加数据, JENNER^[68]指出, 方位分扇叠加数据会造成: ①丢失了部分方位信息; ②即使某个方位扇具有较低的覆盖次数, 即该方位扇叠加数据具有较低的信噪比, 但反演时依旧认为每个方位扇具有相同的权重, ③方位扇分析结果经常受到采集系统的影响。据此, JENNER改写Ruger反射系数方程为:

$ \begin{gathered} R\left(\theta, \varphi_{\mathrm{obs}}\right)=P+\left[W_{11} \cos ^{2} \varphi_{\mathrm{obs}}+2 W_{12} \cos \varphi_{\mathrm{obs}}\right. \\ \left.\sin \varphi_{\mathrm{obs}}+W_{22} \sin ^{2} \varphi_{\mathrm{obs}}\right] \sin ^{2} \theta \end{gathered} $

(17)

式中: W_ij是裂缝弱度和方位的函数。基于不分扇的方位数据, 采用最小二乘算法, 估计参数W_ij和截距P, 各向异性梯度和裂缝方位依次可估计为:

$ G_{\mathrm{ani}}=\sqrt{\left(W_{11}-W_{22}\right)^{2}+4 W_{12}^{2}} $

$ \varphi_{\mathrm{sym}}=\arctan \left[\frac{W_{11}-W_{22}+\sqrt{\left(W_{11}-W_{22}\right)^{2}+4 W_{12}^{2}}}{2 W_{12}}\right] $

2) 各向同性背景中包含一组倾斜定向裂缝(TTI)假设。CHEN等^[69]利用Bond变换研究了倾斜定向裂缝分布在各向同性基岩情况(图 9)的刚度矩阵, 然后利用逆散射理论, 在小入射角和裂缝倾角大于70°的假设下, 推导了反射系数方程, 即:

$ \begin{aligned} R&_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right) \approx R_{\mathrm{PP}}^{\mathrm{iso}}(\theta)+\left[2 \gamma^{2} \sin ^{2} \theta \tan ^{2} \theta \sin ^{2} \varphi_{\mathrm{obs}}\right. \\ &\left.\cos ^{2} \varphi_{\mathrm{obs}}-\gamma \tan ^{2} \theta \cos ^{2} \varphi_{\mathrm{obs}}-2 \gamma\left(\gamma-\frac{1}{2}\right)\right] \beta_{\mathrm{N}}+ \\ &{\left[\gamma^{2} \sin ^{2} \theta \tan ^{2} \theta \cos ^{4} \varphi_{\mathrm{obs}}+2 \gamma^{2} \sin ^{2} \theta \cos ^{2} \varphi_{\mathrm{obs}}+\right.} \\ &\left.2 \gamma\left(\gamma-\frac{1}{2}\right)-\left(\gamma-\frac{1}{2}\right)^{2} \sec ^{2} \theta\right] \delta_{\mathrm{N}}+ \\ &\gamma\left(\sin ^{2} \theta \tan ^{2} \theta \sin ^{2} \varphi_{\mathrm{obs}} \cos 2 \varphi_{\mathrm{obs}}+\sin ^{2} \theta \cos ^{2} \varphi_{\mathrm{obs}}+\right. \\ &\left.\sin ^{2} \theta \cos 2 \varphi_{\mathrm{obs}}+\cos ^{2} \varphi_{\mathrm{obs}}\right) \beta_{\mathrm{T}}+ \\ &\gamma\left(\sin ^{2} \theta \sin ^{4} \varphi_{\mathrm{obs}}-\tan ^{2} \theta \cos ^{4} \varphi_{\mathrm{obs}}-\cos ^{2} \theta\right) \delta_{\mathrm{T}} \end{aligned} $

(18)

式中: R_PP^iso(θ)代表反射系数各向同性部分; β_N=δ_Nsin²ν, β_T=δ_Tsin²ν; ν代表裂缝倾角。基于振幅方位差异反演策略, 构建三步法预测弱度参数和裂缝倾角, 具体步骤为: ①将弱度参数设为零, 反演β_N和β_T; ②用原始地震数据与β_N和β_T产生的新数据差异反演弱度参数, 且在反演弱度参数时, 将预测的β_N和β_T作为初始模型; ③估计的裂缝倾角为$ \nu = {\rm{arcsin}}(\sqrt {{\beta _{\rm{N}}}/{\delta _{\rm{N}}}} )$。

3) 裂缝诱导的正交各向异性介质(OA)假设。裂缝诱导正交各向异性包括: ①VTI背景包含一组垂直定向裂缝; ②各向同性背景包含两组垂直定向旋转不变裂缝; ③VTI背景包含两组垂直定向旋转不变裂缝; ④各向同性背景包含一组水平旋转不变裂缝和一组垂直定向旋转不变裂缝等情况, 如图 10所示。

针对情况①, CHEN等^[71]采用逆散射理论可以获得反射系数方程, 即

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right)=a(\theta)\left(\frac{1}{2} \cdot \frac{\Delta I_{\mathrm{P}}}{I_{\mathrm{P}}}\right)+b(\theta)\left(\frac{1}{2} \cdot \frac{\Delta I_{\mathrm{S}}}{I_{\mathrm{S}}}\right)+ \\ c(\theta)\left(\frac{1}{2} \cdot \frac{\Delta \rho}{\rho}\right)+e(\theta)\left(\frac{1}{2} \Delta \varepsilon_{\mathrm{b}}\right)+f(\theta)\left(\frac{1}{2} \Delta \delta_{\mathrm{b}}\right)+ \\ l\left(\theta, \varphi_{\mathrm{obs}}\right)\left(\frac{1}{2} \Delta_{\delta \mathrm{N}}\right)+m\left(\theta, \varphi_{\mathrm{obs}}\right)\left(\frac{1}{2} \Delta_{\delta \mathrm{V}}\right)+ \\ n\left(\theta, \varphi_{\mathrm{obs}}\right)\left(\frac{1}{2} \Delta_{\delta H}\right) \end{gathered} $

(19)

其中,

$ a(\theta)=\sec ^{2} \theta, b(\theta)=-8 \gamma_\text{b} \sin ^{2} \theta $

$ c(\theta)=4 \gamma_\text{b} \sin ^{2} \theta-\tan ^{2} \theta $

$ e(\theta)=\sin ^{2} \theta \tan ^{2} \theta, f(\theta)=\sin ^{2} \theta $

$ \begin{gathered} l\left(\theta, \varphi_{\text {obs }}\right)=-\frac{1}{2 \cos ^{2} \theta}\left[1-2 \gamma_{\mathrm{b}}\left(\sin ^{2} \theta \sin ^{2} \varphi_{\mathrm{obs}}+\right.\right. \\ \left.\left.\cos ^{2} \theta\right)\right]^{2} \end{gathered} $

$ m\left(\theta, \varphi_{\mathrm{obs}}\right)=2 \gamma_{\mathrm{b}} \sin ^{2} \theta \cos ^{2} \varphi_{\mathrm{obs}} $

$ n\left(\theta, \varphi_{\mathrm{obs}}\right)=-2 \gamma_{\mathrm{b}} \sin ^{2} \theta \tan ^{2} \theta \sin ^{2} \varphi_{\mathrm{obs}} \cos ^{2} \varphi_{\mathrm{obs}} $

CHEN等^[71]进一步推导了弹性阻抗方程, 并基于方位地震数据, 采用马尔科夫-蒙特卡洛反演算法实现了各向同性背景模量、背景各向异性参数和裂缝弱度参数的反演。PAN等^[72]引入准弱度参数和相对各向异性参数来线性化方位弹性阻抗方程, 提高参数反演稳定性。对于情况②和③, PAN等^[11]利用散射理论分别推导了线性化反射系数方程。

对于情况④, ZHANG等^[73]推导了线性化反射系数方程, 并将其表达为傅里叶级数形式, 即:

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right)=a_{0}(\theta)+a_{2}(\theta) \cos \left(2 \varphi_{\mathrm{obs}}\right)+b_{2}(\theta) \cdot \\ \sin \left(2 \varphi_{\mathrm{obs}}\right)+a_{4}(\theta) \cos \left(4 \varphi_{\mathrm{obs}}\right)+b_{4}(\theta) \sin \left(4 \varphi_{\mathrm{obs}}\right) \\ =r_{0}(\theta)+r_{2}(\theta) \cos \left\{2\left[\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right]\right\}+ \\ r_{4}(\theta) \cos \left\{4\left[\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right]\right\} \end{gathered} $

(20)

其中, 加权振幅表示为:

$ \begin{gathered} a_{0}(\theta)=a^{\prime}(\theta) \frac{\Delta M_{\mathrm{b}}}{M_{\mathrm{b}}}+b^{\prime}(\theta) \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}+c^{\prime}(\theta) \frac{\Delta \rho}{\rho}+ \\ d^{\prime}(\theta) \Delta \delta_{\mathrm{N}}^{\mathrm{HOR}}+e^{\prime}(\theta) \Delta \delta_{\mathrm{T}}^{\mathrm{HOR}}+h^{\prime}(\theta) \Delta \delta_{\mathrm{N}}^{\mathrm{VER}}+ \\ i^{\prime}(\theta) \Delta \delta_{\mathrm{T}}^{\mathrm{VER}} \end{gathered} $

(21a)

$ \begin{gathered} a_{2}(\theta)=j^{\prime}(\theta) \cos \left(2 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{N}}^{\mathrm{VER}}+k^{\prime}(\theta) \cdot \\ \cos \left(2 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{T}}^{\mathrm{VER}} \end{gathered} $

(21b)

$ \begin{gathered} b_{2}(\theta)=j^{\prime}(\theta) \sin \left(2 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{N}}^{\mathrm{VER}}+k^{\prime}(\theta) \cdot \\ \sin \left(2 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{T}}^{\mathrm{VER}} \end{gathered} $

(21c)

$ \begin{gathered} a_{4}(\theta)=l^{\prime}(\theta) \cos \left(4 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{N}}^{\mathrm{VER}}+m^{\prime}(\theta) \cdot\\ \cos \left(4 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{T}}^{\mathrm{VER}} \end{gathered} $

(21d)

$ \begin{gathered} b_{4}(\theta)=l^{\prime}(\theta) \sin \left(4 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{N}}^{\mathrm{VER}}+m^{\prime}(\theta) \cdot \\ \sin \left(4 \varphi_{\mathrm{sym}}\right) \Delta \delta_{\mathrm{T}}^{\mathrm{VER}} \end{gathered} $

(21e)

式中: δ_N^VER, δ_T^VER, δ_N^HOR和δ_T^HOR分别为垂直和水平裂缝的法向和切向弱度参数。公式(20)中解耦了各向同性背景、水平裂缝、垂直裂缝信息, 加权振幅a_i和b_i也可通过离散傅里叶变换估计。ZHANG等^[73]经过分析发现, 二阶傅里叶系数r₂对垂直裂缝弱度参数敏感性大于零阶傅里叶系数r₀和四阶傅里叶系数r₄。所以采用二阶傅里叶系数反演垂直裂缝弱度参数, 利用零阶傅里叶系数反演各向同性背景参数和水平裂缝弱度参数, 最后裂缝方位可估计为: φ_sym=1/2arctan[b₂(θ)/a₂(θ)], 且需要裂缝方位先验信息以消除90°模糊性。BARONE等^[74]也将傅里叶级数形式的反射系数方程应用到Haynesville页岩工区(该工区可看作VTI背景中包含一组垂直定向裂缝, 即情况①), 并取得了较好的效果。对于复杂裂缝介质描述, 构建岩石物理模型, 为未知参数提供低频约束, 也是提高反演稳定性的有效途径^[12]。

2.2 基于连通孔缝低频岩石物理的裂缝预测方法

本节将从裂缝诱导VFI介质(一种HTI型介质)、TTI介质、OA介质三个角度分别阐述没有衰减特征或衰减特征不明显、孔缝连通储层的裂缝参数五维地震预测方法研究进展。

2.2.1 各向同性背景中包含一组垂直定向裂缝(VFI)假设

GUREVICH^[75]和HUANG等^[76]基于广义Gassmann流体替换方程, 研究了饱和VFI介质刚度矩阵。假设含气致密砂岩或页岩储层满足: ①流体体积模量远小于矿物模量; ②小裂缝密度; ③小孔隙度; ④等应变介质假设, 即K_b^dry=K_g(1－ϕ), 其中ϕ代表孔隙度。则饱和VFI介质刚度系数可简化为^[77]:

$ C_{11}^{\mathrm{sat}} \approx M_{\mathrm{b}}^{\text {dry }}\left(1-\delta_{\mathrm{N}}^{\text {dry }}\right)+\phi K_{f}+2 \delta_{\mathrm{N}}^{\text {dry }} K_{f} $

(22a)

$ C_{12}^{\mathrm{sat}} \approx \lambda_{\mathrm{b}}^{\text {dry }}\left(1-\delta_{\mathrm{N}}^{\text {dry }}\right)+\phi K_{f}+\left(\chi_{\mathrm{b}}^{\text {dry }}+1\right) \delta_{\mathrm{N}}^{\text {dry }} K_{f} $

(22b)

$ C_{23}^{\mathrm{sat}} \approx \lambda_{\mathrm{b}}^{\text {dry }}\left(1-\chi_{\mathrm{b}}^{\text {dry }} \delta_{\mathrm{N}}^{\text {dry }}\right)+\phi K_{f}+2 \chi_{\mathrm{b}}^{\text {dry }} \delta_{\mathrm{N}}^{\text {dry }} K_{f} $

(22c)

$ C_{33}^{\mathrm{sat}} \approx M_{\mathrm{b}}^{\text {dry }}\left[1-\left(\chi_{\mathrm{b}}^{\text {dry }}\right)^{2} \delta_{\mathrm{N}}^{\text {dry }}\right]+\phi K_{f}+2 \chi_{\mathrm{b}}^{\text {dry }} \delta_{\mathrm{N}}^{\text {dry }} K_{f} $

(22d)

$ C_{44}^{\text {sat }}=\mu_{\mathrm{b}}^{\text {dry }} $

(22e)

$ C_{55}^{\mathrm{sat}}=\mu_{\mathrm{b}}^{\text {dry }}\left(1-\delta_{\mathrm{T}}^{\text {dry }}\right) $

(22f)

然后利用逆散射理论可推导流体体积模量表征的反射系数方程^[77]:

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi_{\mathrm{obs}}\right)=a_{M_{\mathrm{b}}^{\mathrm{dry}}}(\theta) \frac{\Delta M_{\mathrm{b}}^{\mathrm{dry}}}{M_{\mathrm{b}}^{\mathrm{dry}}}+a_{\mu_{\mathrm{b}}^{\mathrm{dry}}} \frac{\Delta \mu_{\mathrm{b}}^{\mathrm{dry}}}{\mu_{\mathrm{b}}^{\mathrm{dry}}}+ \\ a_{\rho}(\theta) \frac{\Delta \rho}{\rho}+a_{f}(\theta)\left(\frac{\Delta K_{f}}{K_{f}}+\frac{\Delta \phi}{\phi}\right)+a_{\delta_{\mathrm{N}}^{\mathrm{dry}}}(\theta \\ \left.\varphi_{\mathrm{obs}}\right) \Delta_{\delta \mathrm{N}}+a_{\delta_{\mathrm{T}}^{\text {dry }}}\left(\theta, \varphi_{\mathrm{obs}}\right) \Delta_{\delta_{\mathrm{T}}} \end{gathered} $

(23)

其中,

$ a_{M_{\mathrm{b}}^{\mathrm{dry}}}(\theta)=\frac{1}{4 \cos ^{2} \theta} \quad a_{\mu_{\mathrm{b}}^{\mathrm{dry}}}(\theta)=-2 \gamma_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \theta $

$ a_{\rho}(\theta)=\frac{\cos 2 \theta}{4 \cos ^{2} \theta} \quad a_{f}(\theta)=\frac{1}{4 \cos ^{2} \theta}\left(1-\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}}\right) $

$ \begin{gathered} a_{\delta_{\mathrm{N}}^{\mathrm{dry}}}\left(\theta, \varphi_{\mathrm{obs}}\right)=-\frac{1}{4 \cos ^{2} \theta} \frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}}\left\{1-2 \gamma_{\mathrm{b}}^{\mathrm{dry}}\left[\sin ^{2} \theta\right.\right. \\ \left.\left.\sin ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)+\cos ^{2} \theta\right]\right\}^{2} \end{gathered} $

$ \begin{gathered} a_{\delta_{\mathrm{T}}^{\mathrm{dry}}}\left(\theta, \varphi_{\mathrm{obs}}\right)=-\gamma_{\mathrm{b}}^{\mathrm{sat}} \tan ^{2} \theta \cos ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right) \\ {\left[\sin ^{2} \theta \sin ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)-\cos ^{2} \theta\right]^{2}} \end{gathered} $

δ_T^dry和δ_N^dry是干岩石裂缝弱度参数, 可直接作为裂缝密度指示因子; K_f代表流体体积模量, 可用于区分流体类型。公式(23)实现了孔裂隙流体和孔裂隙度及岩石骨架的解耦, 可以直接反演K_f和裂缝弱度参数。其实介质等应变假设会导致比较大的误差, 所以在只满足前两条假设时, 饱和VFI介质刚度系数可简化为^[78]

$ C_{11}^{\mathrm{sat}} \approx\left(\lambda_{\mathrm{b}}^{\mathrm{dry}}+2 \mu_{\mathrm{b}}\right)\left(1-\delta_{\mathrm{N}}^{\mathrm{dry}}\right)+f $

(24a)

$ C_{12}^{\mathrm{sat}} \approx \lambda_{\mathrm{b}}^{\mathrm{dry}}\left(1-\delta_{\mathrm{N}}^{\mathrm{dry}}\right)+f $

(24b)

$ C_{23}^{\mathrm{sat}} \approx \lambda_{\mathrm{b}}^{\mathrm{dry}}\left(1-\chi_{\mathrm{b}}^{\mathrm{dry}} \delta_{\mathrm{N}}^{\mathrm{dry}}\right)+f $

(24c)

$ C_{33}^{\mathrm{sat}} \approx\left(\lambda_{\mathrm{b}}^{\mathrm{dry}}+2 \mu_{\mathrm{b}}\right)\left[1-\left(\chi_{\mathrm{b}}^{\mathrm{dry}}\right)^{2} \delta_{\mathrm{N}}^{\mathrm{dry}}\right]+f $

(24d)

式中: f≈K_f(1－K_b^dry/K_g)²/ϕ。由于f受孔隙度的影响较大, 所以CHEN等^[78]提出一个新的VFI介质流体指示因子F=f/ϕ, 该指示因子受孔隙度的影响非常小, 将其带入公式(24), 便可推导F和干裂缝弱度参数表示的线性反射系数方程以实现裂缝与流体检测。PAN等^[79]也仅在前两条假设下, 使用临界孔隙度模型研究了关于等效流体体积模量K_f和干裂缝弱度参数的线性反射系数方程及裂缝和流体类型预测方法。另外, YIN等^[80]根据RUSSELL等^[81]流体因子的定义方式, 构建了一个融合裂缝影响的各向异性Russell流体因子, 实现裂缝储层流体识别。

当仅考虑VFI介质中背景孔隙之间的流体流动时, 可采用各向同性Gassman方程, 将基于Schoenberg理论表达的VFI介质刚度矩阵中的各向同性背景岩石模量写为Gassmann流体因子与干岩石骨架模量的和, 即M_b=f+γ_b^dryμ, γ_b^dry为背景干岩石骨架横纵波速度比的平方。然后根据推导的刚度矩阵表达式, 便可推导相应的反射系数方程^[82], 采用合理的反演算法, 实现流体指示因子f及裂缝弱度参数的反演。

2.2.2 各向同性背景中包含一组倾斜定向裂缝(TTI)假设

对于TTI型裂缝储层, PAN等^[70]仅考虑背景孔隙之间的流体流动, 采用各向同性Gassman方程, 将VFI介质的各向同性背景岩石模量写为Gassmann流体因子与干岩石骨架模量的和, 然后利用Bond变换得到TTI介质刚度系数, 即:

$ \begin{aligned} &C_{11}^{\mathrm{sat}}=\left(f+\gamma_{\mathrm{b}}^{\mathrm{sat}} \mu_{\mathrm{b}}\right)\left(\sin ^{4} \nu+\cos ^{4} \nu\right)-M_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}\left[\left(\chi_{\mathrm{b}}^{\mathrm{sat}}\right)^{2} \cdot\right. \\ &\left.\sin ^{4} \nu+\cos ^{4} \nu\right]+\sin ^{2} 2 \nu\left[\frac{1}{2} M_{\mathrm{b}}^{\mathrm{sat}}+\frac{1}{2} \lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}-\mu_{\mathrm{b}} \delta_{\mathrm{T}}\right] \end{aligned} $

(25a)

$ C_{12}^{\mathrm{sat}}=f+\left(\gamma_{\mathrm{b}}^{\mathrm{sat}}-2\right) \mu-\lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}\left(\chi_{\mathrm{b}}^{\mathrm{sat}} \cos ^{2} \nu+\sin ^{2} \nu\right) $

(25b)

$ \begin{gathered} C_{13}^{\mathrm{sat}}=f+\left(\gamma_{\mathrm{b}}^{\mathrm{sat}}-2\right) \mu_{\mathrm{b}}-\lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}+ \\ \frac{1}{4} \sin ^{2} 2 \nu\left\{2 \lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}-\left[1+\left(\chi_{\mathrm{b}}^{\mathrm{sat}}\right)^{2}\right] M_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}+4 \mu_{\mathrm{b}} \delta_{\mathrm{T}}\right\} \end{gathered} $

(25c)

$ \begin{aligned} C_{15}^{\mathrm{sat}}=& \frac{1}{4} \sin 4 \nu\left(\lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}+2 \mu_{\mathrm{b}} \delta_{\mathrm{T}}\right)+\frac{1}{2} \sin ^{2} 2{\nu} \\ & {\left[\delta_{\mathrm{N}} M_{\mathrm{b}}^{\mathrm{sat}}\left(\sin ^{2} \nu-\chi_{\mathrm{b}}^{\mathrm{sat}} \cos ^{2} \nu\right)\right] } \end{aligned} $

(25d)

$ C_{22}^{\mathrm{sat}}=f+\gamma_{\mathrm{b}}^{\mathrm{sat}} \mu_{\mathrm{b}}-\left(\chi_{\mathrm{b}}^{\mathrm{sat}}\right)^{2} M_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}} $

(25e)

$ C_{23}^{\mathrm{sat}}=f+\left(\gamma_{\mathrm{b}}^{\mathrm{sat}}-2\right) \mu_{\mathrm{b}}-\lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}\left(\chi_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \nu+\cos ^{2} \nu\right) $

(25f)

$ C_{25}^{\mathrm{sat}}=\frac{1}{2}\left[f+\left(\gamma_{\mathrm{b}}^{\mathrm{sat}}-2\right) \mu_{\mathrm{b}}-\chi_{\mathrm{b}}^{\mathrm{sat}} \lambda_{\mathrm{b}}^{\mathrm{sat}}\right] \delta_{\mathrm{N}} \sin 2 \nu $

(25g)

$ C_{44}^{\mathrm{sat}}=\mu_{\mathrm{b}}-\mu_{\mathrm{b}} \delta_{\mathrm{T}} \cos ^{2} \nu $

(25h)

$ C_{46}^{\mathrm{sat}}=\frac{1}{2} \mu_{\mathrm{b}} \delta_{\mathrm{T}} \sin 2 \nu $

(25i)

$ \begin{gathered} C_{35}^{\mathrm{sat}}=-\frac{1}{2}\left(f+\gamma_{\mathrm{b}}^{\mathrm{sat}} \mu_{\mathrm{b}}\right) \delta_{\mathrm{N}} \sin 2 \nu\left\{\left(\chi_{\mathrm{b}}^{\mathrm{sat}}\right)^{2} \sin ^{2} \nu+\right. \\ \left.\cos ^{2} \nu\right]+\frac{1}{4} \sin 4 \nu\left(2 \mu_{\mathrm{b}} \delta_{\mathrm{T}}-\lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}\right) \end{gathered} $

(25j)

$ \begin{gathered} C_{55}^{\mathrm{sat}}=\frac{1}{4} \sin ^{2} 2 \nu\left\{2 \lambda_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}-M_{\mathrm{b}}^{\mathrm{sat}} \delta_{\mathrm{N}}\left[\left(\chi_{\mathrm{b}}^{\mathrm{sat}}\right)^{2}+1\right]-\right. \\ \left.4 \mu_{\mathrm{b}}\right\}+\cos ^{2} 2 \nu\left[\mu_{\mathrm{b}}\left(1-\delta_{\mathrm{T}}\right)\right\} \end{gathered} $

(25k)

$ C_{66}^{\mathrm{sat}}=\mu_{\mathrm{b}}-\mu_{\mathrm{b}} \delta_{\mathrm{T}} \sin ^{2} \nu $

(25l)

然后推导出Gassmann流体因子f和裂缝弱度参数线性表征的反射系数方程和弹性阻抗方程, 在获得裂缝倾角ν先验信息的前提下, 实现TTI型裂缝储层流体分布和裂缝特征的线性化反演^[70]。

2.2.3 裂缝诱导的正交各向异性介质(OA)假设

对于各向同性岩石中包含两组垂直正交定向裂缝系统的情况, 可以假设每组定向裂缝系统诱导的各向异性比较微弱, 则岩石整体的等效刚度系数可表达为等效流体体积模量和裂缝弱度参数的函数^[79]:

$ \begin{gathered} C_{11}^{\text {sat }}=M_{\mathrm{b}}^{\mathrm{dry}}\left[1-\delta_{\mathrm{N} 1}-\left(\chi_{\mathrm{b}}^{\mathrm{dry}}\right)^{2} \delta_{N 2}\right]+G_{n}(\phi) K_{f}+\\ 2 F_{n}(\phi) K_{f}\left(\delta_{\mathrm{N} 1}+\chi_{\mathrm{b}}^{\mathrm{dry}} \delta_{\mathrm{N} 2}\right) \end{gathered} $

(26a)

$ \begin{gathered} C_{12}^{\mathrm{sat}}=\lambda_{\mathrm{b}}^{\mathrm{dry}}\left(1-\delta_{\mathrm{N} 1}-\delta_{\mathrm{N} 2}\right)+G_{n}(\phi) K_{f}+F_{n}(\phi) \cdot\\ K_{f}\left(1+\chi_{\mathrm{b}}^{\text {dry }}\right)\left(\delta_{\mathrm{N} 1}+\delta_{\mathrm{N} 2}\right) \end{gathered} $

(26b)

$ \begin{gathered} C_{13}^{\mathrm{sat}}=\lambda_{\mathrm{b}}^{\mathrm{dry}}\left(1-\delta_{\mathrm{N} 1}-\chi_{\mathrm{b}}^{\mathrm{dry}} \delta_{\mathrm{N} 2}\right)+G_{n}(\phi) K_{f}+\\ F_{n}(\phi) K_{f}\left(1+\chi_{\mathrm{b}}^{\mathrm{dry}}\right) \delta_{\mathrm{N} 1}+2 F_{n}(\phi) K_{f} \chi_{\mathrm{b}}^{\mathrm{dry}} \delta_{\mathrm{N} 2} \end{gathered} $

(26c)

$ \begin{gathered} C_{22}^{\text {sat }}=M_{\mathrm{b}}^{\text {dry }}\left[1-\left(\chi_{\mathrm{b}}^{\text {dry }}\right)^{2} \delta_{\mathrm{N} 1}-\delta_{\mathrm{N} 2}\right]+G_{n}(\phi) K_{f}+\\ 2 F_{n}(\phi) K_{f}\left(\chi_{\mathrm{b}}^{\text {dry }} \delta_{\mathrm{N} 1}+\delta_{\mathrm{N} 2}\right) \end{gathered} $

(26d)

$ \begin{gathered} C_{23}^{\text {sat }}=\lambda_{\mathrm{b}}^{\text {dry }}\left(1-\chi_{\mathrm{b}}^{\text {dry }} \delta_{\mathrm{N} 1}-\delta_{\mathrm{N} 2}\right)+G_{n}(\phi) K_{f}+\\ 2 F_{n}(\phi) K_{f} \chi_{\mathrm{b}}^{\mathrm{dry}} \delta_{\mathrm{N} 1}+F_{n}(\phi) K_{f}\left(1+\chi_{\mathrm{b}}^{\mathrm{dry}}\right) \delta_{\mathrm{N} 2} \end{gathered} $

(26e)

$ \begin{gathered} C_{33}^{\mathrm{sat}}=M_{\mathrm{b}}^{\mathrm{dry}}\left[1-\left(\chi_{\mathrm{b}}^{\mathrm{dry}}\right)^{2}\left(\delta_{\mathrm{N} 1}+\delta_{\mathrm{N} 2}\right)\right]+G_{n}(\phi) K_{f}+\\ 2 F_{n}(\phi) K_{f} \chi_{\mathrm{b}}^{\mathrm{dry}}\left(\delta_{\mathrm{N} 1}+\delta_{\mathrm{N} 2}\right) \end{gathered} $

(26f)

$ C_{44}^{\mathrm{sat}}=\mu_{\mathrm{b}}\left(1-\delta_{\mathrm{T} 2}\right) $

(26g)

$ C_{55}^{\mathrm{sat}}=\mu_{\mathrm{b}}\left(1-\delta_{\mathrm{T} 2}\right) $

(26h)

$ C_{66}^{\mathrm{sat}}=\mu_{\mathrm{b}}\left(1-\delta_{\mathrm{T} 1}-\delta_{\mathrm{T} 2}\right) $

(26i)

式中: G_n(ϕ)=(1－K_b^dry/K_g)²/ϕ, F_n(ϕ)=(K_b^dry/K_g)(1－K_b^dry/K_g)/ϕ; δ_Ni和δ_Ti分别为第i组裂缝的法向和切向弱度参数。PAN等^[79]利用散射理论获得了等效流体体积模量和两组裂缝弱度参数线性表征的反射系数方程:

$ \begin{aligned} &R_{\mathrm{PP}}(\theta, \varphi)=\left[\left(1-\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}}\right) \frac{\sec ^{2} \theta}{4}\right] \frac{\Delta K_{f}}{K_{f}}+\left(\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}} \cdot\right. \\ &\left.\frac{\sec ^{2} \theta}{4}-2 \gamma_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \theta\right) \frac{\Delta f_{m}}{f_{m}}+\left(\frac{\sec ^{2} \theta}{4}-\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}} \frac{\sec ^{2} \theta}{2}+\right. \\ &\left.2 \gamma_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \theta\right) \frac{\Delta \varphi}{\varphi}+\left(\frac{1}{2}-\frac{\sec ^{2} \theta}{4}\right) \frac{\Delta \rho}{\rho}-\\ &\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}} \frac{\sec ^{2} \theta}{4}\left[2 \gamma_{\mathrm{b}}^{\mathrm{sat}}\left(\sin ^{2} \theta \sin ^{2} \varphi+\cos ^{2} \theta\right)-1\right]^{2} \Delta \delta_{\mathrm{N} 1}+ \\ &\gamma_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \theta \cos ^{2} \varphi\left(1-\tan ^{2} \theta \sin ^{2} \varphi\right) \Delta \delta_{\mathrm{T} 1}- \\ &\frac{\gamma_{\mathrm{b}}^{\mathrm{sat}}}{\gamma_{\mathrm{b}}^{\mathrm{dry}}} \frac{\sec ^{2} \theta}{4}\left[2 \gamma_{\mathrm{b}}^{\mathrm{sat}}\left(\sin ^{2} \theta \cos ^{2} \varphi+\cos ^{2} \theta\right)-1\right]^{2} \Delta \delta_{\mathrm{N} 2}+ \\ &\gamma_{\mathrm{b}}^{\mathrm{sat}} \sin ^{2} \theta \sin ^{2} \varphi\left(1-\tan ^{2} \theta \cos ^{2} \varphi\right) \Delta \delta_{\mathrm{T} 2} \end{aligned} $

(27)

式中: f_m=φμ_b。基于公式(27), 采用合理的反演算法便可实现裂缝储层定量表征。类似于TTI介质和VFI介质, 对于各向同性背景包含两组垂直正交裂缝系统或VTI背景包含一组垂直裂缝系统诱导的OA介质, 我们也可将背景各向同性部分的饱和岩石模量分解为Gassman流体因子与干岩石骨架部分, 再推导Gassman流体因子与裂缝弱度参数表征的刚度矩阵和线性化反射系数方程^[82], 最后采用合理的反演算法实现裂缝储层表征和流体预测。

2.3 基于衰减岩石物理理论的裂缝预测方法

MAULTZSCH等^{[5, 83]}在研究Clair油田VSP资料时发现, 沿裂缝法向和走向传播的地震波能量衰减值不一样, 如图 11所示。EKANEM等^[84]和尹志恒等^[85-86]通过物理模型模拟方法也发现裂缝走向的品质因子大于裂缝法向。所以在包含定向裂缝的岩石里品质因子将随观测方位变化而变化, 将这种现象称为衰减各向异性(Q-factor versus offset and azimuth, QVOA)^[87-91]。VASCONCELOS等^[92]也证实, 衰减各向异性特征要强于速度方位各向异性, 因此, 利用QVOA反演有利于改善裂缝预测精度^[93]。

对于VFI介质, CHICHININA等^[94-98]利用复值弱度、ZHU等^[99-100]利用衰减各向异性参数均将衰减因子α_Q近似表达为:

$ \alpha_{Q} \approx A+B\left(\varphi_{\mathrm{obs}}\right) \sin ^{2} \theta+C\left(\varphi_{\mathrm{obs}}\right) \sin ^{4} \theta $

(28)

其中,

$ B\left(\varphi_{\mathrm{obs}}\right)=B_{0}+B_{1} \cos ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right) $

$ \begin{gathered} C\left(\varphi_{\mathrm{obs}}\right)=C_{0}+C_{1} \cos ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)+ \\ C_{2} \cos ^{4}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right) \end{gathered} $

并且很多数值实验和分析表明, 在大多数情况下, 裂缝法向的衰减系数大于裂缝走向的衰减系数。所以结合不同方位的衰减系数预测值和公式(28), 利用最小二乘算法, 可实现裂缝方位φ_sym的预测^{[90, 101]}。其中, 不同方位衰减系数可以采用谱比法提取^[98], 具体地讲, 计算地层上下界面同相轴频谱振幅比的自然对数, 即:

$ \ln S(f)=\ln \frac{A_{\mathrm{b}}(f)}{A_{t}(f)}=\eta f+R_{0} $

(29)

式中: A_t(f)和A_b(f)分别为顶底界面同相轴的频谱振幅; f为频率, η=τd, τ为层间旅行时, 可采用射线追踪算法进行估计。如图 12所示, 选用红线标识的代表频带中lnS-f的关系计算η, 然后即可估计衰减系数α_Q=$1 - {{\rm{e}}^{ - 2\frac{\eta }{\tau }}} $。但是BARBOSA等^[102]研究表明, 对于多孔饱和裂缝介质中流体流动诱导地震波衰减的情况, 在某些条件下, 其衰减值在裂缝法向和走向是相等的, 此时可能会导致基于QVOA预测裂缝方位的90°模糊性, 所以需要裂缝方位先验信息作为约束。

图 13为衰减系数α_Q与sin²θ的关系曲线。由图可见, 衰减系数α_Q对sin²θ的依赖性并非是二次的, 所以基于公式(28)的裂缝方位预测结果可能存在较大的误差。因此, SABININ^[101]和CHICHININA等^[98]分别采用更为精确的衰减系数方程, 即:

$ \begin{aligned} \alpha_{Q}&=A+\sin ^{2} \theta\left[B_{0}+B_{1} \cos ^{2}\left(\varphi_{\text {obs }}-\varphi_{\text {sym }}\right)\right]+\\ &\sin ^{4} \theta\left[C_{0}+C_{1} \cos ^{2}\left(\varphi_{\text {obs }}-\varphi_{\text {sym }}\right)+C_{2} \cos ^{4}\left(\varphi_{\text {obs }}-\right.\right.\\ &\left.\left.\varphi_{\mathrm{sym}}\right)\right] /\left\{1+\sin ^{2} \theta\left[D_{0}+D_{1} \cos ^{2}\left(\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}\right)\right]+\right.\\ &\sin ^{4} \theta\left[E_{0}+E_{1} \cos ^{2}\left(\varphi_{\text {obs }}-\varphi_{\mathrm{sym}}\right)+\right.\\ &\left.\left.E_{2} \cos ^{4}\left(\varphi_{\text {obs }}-\varphi_{\text {sym }}\right)\right]\right\} \end{aligned} $

(30)

结合多方位衰减因子估计值, 利用最小二乘算法实现裂缝方位估计。公式(30)中, A, B₀, B₁, C₀, C₁, C₂, D₀, D₁, E₀, E₁和E₂分别为储层属性的函数, 为未知变量。

不同于QVOA裂缝参数预测, CHEN等^[6]采用修正的Schoenberg线性滑动理论研究部分饱和定向裂缝介质。在气体粘度为零、体积模量极小但气体体积分数不小的假设下, 简化修正Schoenberg线性滑动理论。在含油储层裂缝参数地震预测的现实需求下, 对简化后的复值$ {{\tilde U}_{11}}$和$ {{\tilde U}_{33}}$进行分析, 发现二者的虚部均可忽略, 可得到部分饱和定向裂缝诱导的VFI介质刚度系数为:

$ C_{11}=M_{\mathrm{b}}-M_{\mathrm{b}} \delta_{\mathrm{N}}^{\text {dry }}-\omega M_{\mathrm{b}} \varPsi_{n} $

(31a)

$ C_{12}=\lambda_{\mathrm{b}}-\lambda_{\mathrm{b}} \delta_{\mathrm{N}}^{\text {dry }}-\omega^{2} \lambda_{\mathrm{b}} \varPsi_{n} $

(31b)

$ C_{23}=\lambda_{\mathrm{b}}-\lambda_{\mathrm{b}}\left(1-2 \gamma_{\mathrm{b}}\right) \delta_{\mathrm{N}}^{\mathrm{dry}}-\omega^{2} \lambda_{\mathrm{b}}\left(1-2 \gamma_{\mathrm{b}}\right) \varPsi_{n} $

(31c)

$ C_{33}=M_{\mathrm{b}}-\lambda_{\mathrm{b}}\left(1-2 \gamma_{\mathrm{b}}\right) \delta_{\mathrm{N}}^{\text {dry }}-\omega^{2} \lambda_{\mathrm{b}}\left(1-2 \gamma_{\mathrm{b}}\right) \varPsi_{n} $

(31d)

$ C_{44}=\mu_{\mathrm{b}} $

(31e)

$ C_{55}=\mu_{\mathrm{b}}-\mu_{\mathrm{b}} \delta_{\mathrm{T}}^{\mathrm{dry}} $

(31f)

式中: Ψ_n=δ_N^dry(η_l/πμ_b)(1/α_ratio³){0.053[1+cos(π－qπ)]/[(1－q)(1－γ)]}, 可作为部分饱和VFI介质的流体指示因子; δ_N^dry和δ_T^dry分别代表干裂缝法向和切向弱度, η_l和q分别代表裂缝中液体的粘性因子和体积分数。然后可推导得到反射系数方程^[6]:

$ \begin{gathered} R_{\mathrm{PP}}\left(\theta, \varphi=\varphi_{\mathrm{obs}}-\varphi_{\mathrm{sym}}, \omega\right)=\frac{\cos 2 \theta}{4 \cos ^{2} \theta} \frac{\Delta \rho}{\rho}+\frac{1}{4 \cos ^{2} \theta} \\ \frac{\Delta M_{\mathrm{b}}}{M_{\mathrm{b}}}-2 \gamma_{\mathrm{b}} \sin ^{2} \theta \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}-\frac{1}{4 \cos ^{2} \theta}\left[1-2 \gamma_{\mathrm{b}}\left(\sin ^{2} \theta\right.\right. \\ \left.\left.\sin ^{2} \varphi+\cos ^{2} \theta\right)\right]^{2} \Delta \delta_{\mathrm{N}}^{\mathrm{dry}}-\frac{\omega^{2}}{4 \cos ^{2} \theta}[1- \\ \left.2 \gamma_{\mathrm{b}}\left(\sin ^{2} \theta \sin ^{2} \varphi+\cos ^{2} \theta\right)\right]^{2} \Delta \psi_{n}+ \\ \gamma_{\mathrm{b}}\left(\sin ^{2} \theta \cos ^{2} \varphi-\sin ^{2} \theta \tan ^{2} \theta \sin ^{2} \varphi \cos ^{2} \varphi\right) \Delta \delta_{\mathrm{T}}^{\mathrm{dry}} \end{gathered} $

(32)

基于公式(32), 采用分频地震反演方法^[103]即可实现流体因子Ψ_n和干裂缝弱度参数的预测。

3 应用实例

某工区内A井钻遇优质页岩储层, 该层段表现出早期水体较深、后期持续海退的沉积序列, 水平层理、水平微断层以及高角度裂缝较为发育, 可以将其等效为正交各向异性介质。所以我们采用MA等^[41]提出的基于傅里叶级数展开的正交各向异性介质裂缝参数预测方法进行裂缝特征空间刻画。

首先, 需要进行叠前方位地震数据处理。主要包括精细的排线扩散处理、震检组合效应的校正、反Q滤波、地表一致性处理、叠前剩余振幅补偿、精细的初至切除处理, 分方位角道集集成处理、宽方位速度分析、分方位各向异性叠前偏移处理, 方位叠前道集转换以及部分角度叠加等。最后得到6个方位地震叠加道集(0°~30°, 30°~60°, 60°~90°, 90°~120°, 120°~150°, 150°~180°)以及每个方位包含一个入射角信息(16°)。采用约束稀疏脉冲反演将得到的五维地震数据转换为方位弹性阻抗体, 利用电成像测井数据及地质先验信息模拟得到工区裂缝方位先验数据体, 然后开展相关的裂缝密度和方位预测。

图 14为该工区过水平井的地震剖面, 白色箭头所指的地震同相轴为目的层反射轴。从地震剖面中很难看到裂缝密度的横向变化。图 15为对应于图 14的裂缝密度预测结果, 图中红色曲线代表水平井, 曲线上的编号代表压裂段。图 16为水平井各压裂段的破裂压力统计分析, 显示第11段~第21段的破裂压力值较低, 由于高裂缝密度将导致低的破裂压力, 所以第11段~第21段的裂缝较为发育, 与图 15显示的裂缝密度预测结果比较吻合, 说明裂缝密度预测结果较为合理。图 17为A井目的层段裂缝方位玫瑰统计分析, 图 17a为成像测井解释结果, 图 17b为地震预测结果, 可以发现预测结果和测井解释结果较为吻合, 说明预测结果的正确性。

4 挑战和机遇

目前, 随着缝洞发育的碳酸盐岩储层和页岩储层等非常规油气藏勘探开发的不断深入, 裂缝参数五维地震预测已是国内外研究的热点和难点之一。同时, 地震勘探目标从浅层迈向深层、从勘探迈向开发, 从均匀弹性介质迈向非均匀复杂介质, 由此对裂缝参数地震预测技术的精度和可靠性提出了新的要求。所以, 根据目前油气勘探开发需求和现状以及五维地震裂缝预测技术的瓶颈问题, 未来我们应该在以下3个方面开展进一步的研究。

1) 五维地震数据中的裂缝信息比较微弱, 在地震噪声的干扰下, 实现复杂储层裂缝信息的准确提取是裂缝空间描述的难题之一, 所以构建合理稳定的地震反演策略^{[38, 104-106]}是未来裂缝参数稳定预测的有效途径之一;

2) 五维地震数据中包含丰富的方位、偏移距、振幅、频率和相位等信息, 是地下储层裂缝、流体、孔隙度和矿物等物性参数的综合体现。探索这些信息与储层物性参数之间的关系, 研究地震数据中地层物性参数的解耦方法, 将有助于充分挖掘五维地震数据, 有效提高裂缝检测和流体识别的精度;

3) 裂缝储层反演技术目前基本假定介质弱各向异性和小弹性差异, 从近似各向异性反射特征方程出发实现裂缝参数预测, 然而地下介质复杂多变, 各向异性程度范围较广, 所以基于精确各向异性拟Zoeppritz方程^[107-108]的裂缝储层描述将是未来油气勘探领域发展方向之一, 以适应页岩等各向异性程度较大的目标。