0 引言

能源是国民经济发展的重要支撑，能源安全直接影响国家安全、可持续发展以及社会稳定。近年来，油气勘探、开发不断向深水、深地、非常规等复杂油气领域推进，油气藏特征更复杂，其中最典型的裂缝型储层油气识别成为研究前沿和世界级难题。五维数据反演是以宽方位地震五维数据为基础的反演，五维数据是通过炮检距向量片(Offset Vector Tile，OVT)技术处理宽方位地震(横向与纵向排列尺度的比值大于0.5的地震勘探系统)资料得到的道集资料，其显著特点是同时包含高品质炮检距和方位角信息(即空间三维坐标+炮检距+方位角)^[1]。五维数据中丰富的方位角及炮检距等信息为裂缝型储层地震预测及油气识别奠定了资料基础，但在实际裂缝型储层五维数据解释方面仍存在挑战^[1]，瓶颈在于现有裂缝型储层流体识别方法在构建流体敏感参数过程中没有考虑岩石骨架、油气水及裂缝、孔隙(“固液缝”)耦合影响，对五维数据中丰富的方位信息蕴含的流体信息挖掘不够，并且针对裂缝型储层油气识别以间接计算参数为主^[2-4]。业界主要针对裂缝型储层地震岩石物理流体因子构建及油气识别方法开展了大量基础研究与应用^[5-9]。

裂缝型储层地震流体识别的关键是流体因子构建。基于地震岩石物理模型及储层流体地震响应分析，可将孔隙及裂缝对流体的异常特性表征为各向异性流体因子，进而判识非常规储层含流体类型。人们基于岩石物理理论或裂缝各向异性理论等从不同角度构建了不同的流体敏感参数^[10-11]，并研究地震信息与流体类型之间的关系，推导地震反射系数与流体敏感参数之间的关系进而由反演预测流体^[12-17]。目前利用五维数据研究裂缝型储层含油气性大多以常规油气识别方法在裂缝型储层中的拓展应用为主，包括利用各向异性Gassmann孔隙弹性理论流体体积模量识别含油气性^[18-19]以及基于各向异性Gassmann孔隙弹性理论和线性滑动理论构建各向异性流体因子^[20]。上述方法从岩石物理及各向异性机理出发构建综合考虑“固液缝”解耦的流体因子，为裂缝型储层各向异性流体识别提供了新途径，本文基于流体因子研究直接反演储层参数的方法。

由于裂缝型储层的特殊性及地下埋藏条件的复杂性，给裂缝型储层流体识别带来挑战。宽方位地震五维数据反演问题本身存在“病态解”及“不确定性”等问题，多裂缝参数同时反演更易导致储层参数反演精度下降。因此，基于间接代数组合的流体因子在计算过程中不可避免地累积计算误差，影响油气识别效果^{[6-7, 20]}。为了有效地结合叠前地震反演与流体识别，人们广泛研究了流体因子直接反演问题，主要包括推导包含流体因子的反射系数线性近似公式及研究直接反演方法^{[18-19, 21-24]}等。针对裂缝型储层，弹性阻抗反演和基于贝叶斯框架的AVAZ反演^[25-30]是常用的直接反演方法，除单纯利用地震时间域纵波信号反演外，还包括多波联合反演及非线性反演方法。基于贝叶斯框架的地震反演算法，采取不同的反演策略，充分利用五维数据中多波的特征信息(振幅、频率及相位)直接获取表征储层流体信息异常特性的流体敏感参数，是裂缝型储层流体识别的主要目的。

本文根据裂缝型储层岩石物理及各向异性理论，推导由“固液缝”解耦流体因子直接表征的裂缝型储层方位地震反射系数方程，分析模型的适用条件，从理论上阐明采用该方程获取“固液缝”解耦流体因子的可行性。提出了“固液缝”解耦油气识别方法，采用五维数据时频联合反演方法，基于振幅—频率联合反演策略，在贝叶斯理论框架下构建五维数据振幅—频率联合反演目标泛函，综合利用五维数据的频率和振幅信息及测井、地质资料直接获取“固液缝”解耦流体因子，消除因参数间接运算的累计误差引起的流体识别假象。模型测试及实际资料研究表明，该方法可提高流体识别精度，利用五维数据可直接识别裂缝型储层油气。

1 方法原理 1.1 裂缝型储层“固液缝”解耦流体因子刚度矩阵构建及分析

Gurevich^[31]在Schoenberg线性滑动模型的基础上，考虑裂缝与背景基质孔隙间的流体流动效应，将各向异性Gassmann方程改写为

$ C_{i j}^{\text {sat }}=C_{i j}^{\mathrm{dry}}+\beta_i \beta_j K^* \quad i, j=1, \cdots, 6 $

(1)

式中：C_ij^sat为饱和流体裂缝岩石的等效刚度弹性矩阵元素；C_ij^dry为干裂缝岩石的等效刚度弹性矩阵元素；K^*为各向异性Gassmann孔隙空间模量；β_i和β_j均为各向异性Biot系数，统一用β_m表示为

$ \left\{\begin{array}{l} \beta_m=\left(1-\frac{\sum\limits_{n=1}^3 C_{m n}^{\mathrm{dry}}}{3 K_{\mathrm{g}}}\right) \vartheta_m \\ \vartheta_m= \begin{cases}1 & m=1, 2, 3 \\ 0 & m=4, 5, 6\end{cases} \end{array}\right. $

(2)

式中K_g为固体矿物颗粒等效体积模量。

式(1)中K^*与干裂缝岩石刚度矩阵有关，结合线性滑动理论，有

$ K^*=\frac{K_{\mathrm{g}}}{1-\phi-\left(1-\alpha_0\right)\left(1-\frac{K_{\mathrm{d}}}{M_{\mathrm{d}}} \delta_{\mathrm{N}}\right)+\phi \frac{K_{\mathrm{g}}}{K_{\mathrm{f}}}} $

(3)

式中：ϕ为总孔隙度；K_f为流体体积模量；K_d为背景干岩石有效体积模量；M_d为干岩石纵波模量；α₀为各向同性Biot系数；δ_N为法向裂隙弱度。

将式(2)、式(3)代入式(1)，可得含孔隙裂缝介质饱和流体时的刚度矩阵解析表达式

$ \boldsymbol{C}_{\mathrm{HTI}}^{\text {sat }}=\left(\begin{array}{cccccc} C_{11}^{\text {sat }} & C_{13}^{\text {sat }} & C_{13}^{\text {sat }} & 0 & 0 & 0 \\ C_{13}^{\text {sat }} & C_{33}^{\text {sat }} & C_{23}^{\text {sat }} & 0 & 0 & 0 \\ C_{13}^{\text {sat }} & C_{23}^{\text {sat }} & C_{33}^{\text {sat }} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44}^{\text {sat }} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55}^{\text {sat }} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{55}^{\text {sat }} \end{array}\right) $

(4)

其中

$ \left\{\begin{array}{l} C_{11}^{\mathrm{sat}}=M_{\mathrm{d}}\left(1-\delta_{\mathrm{N}}\right)+\frac{\alpha_0^2}{\phi} K_{\mathrm{f}}+\frac{2 \alpha_0\left(1-\alpha_0\right)}{\phi} K_{\mathrm{f}} \delta_{\mathrm{N}} \\ C_{13}^{\mathrm{sat}}=\lambda_{\mathrm{d}}\left(1-\delta_{\mathrm{N}}\right)+\frac{\alpha_0^2}{\phi} K_{\mathrm{f}}+\frac{\alpha_0\left(1-\alpha_0\right)}{\phi} K_{\mathrm{f}}(1+\chi) \delta_{\mathrm{N}} \\ C_{23}^{\mathrm{sat}}=\lambda_{\mathrm{d}}\left(1-\chi \delta_{\mathrm{N}}\right)+\frac{\alpha_0^2}{\phi} K_{\mathrm{f}}+\frac{2 \alpha_0\left(1-\alpha_0\right)}{\phi} K_{\mathrm{f}} \chi \delta_{\mathrm{N}} \\ C_{33}^{\mathrm{sat}}=M_{\mathrm{d}}\left(1-\chi^2 \delta_{\mathrm{N}}\right)+\frac{\alpha_0^2}{\phi} K_{\mathrm{f}}+\frac{2 \alpha_0\left(1-\alpha_0\right)}{\phi} K_{\mathrm{f}} \chi \delta_{\mathrm{N}} \\ C_{44}^{\mathrm{sat}}=\mu \\ C_{55}^{\mathrm{sat}}=\mu\left(1-\delta_{\mathrm{T}}\right) \end{array}\right. $

(5)

式中：χ=λ_d/M_d，M_d和λ_d分别为干燥岩石纵波模量和拉梅参数；μ为剪切模量；δ_T为切向裂隙弱度。

令$G(\phi)=\alpha_0^2 / \phi$和$F(\phi)=\left[\alpha_0\left(1-\alpha_0\right)\right] / \phi$为与孔隙度相关的两个简化的颗粒函数，在弱各向异性条件下，即裂缝的参数较小、流体体积模量远小于固体颗粒的假设下，进而将“固液缝”流体因子^[29]f_ani改写成另外一种形式

$ f_{\text {ani }}=K_{\mathrm{f}} G(\phi)+2 F(\phi) K_{\mathrm{f}} \chi \delta_{\mathrm{N}} $

(6)

将式(6)代入式(4)，有

$ \left\{\begin{aligned} C_{11}^{\mathrm{sat}} \approx & M_{\mathrm{d}}\left(1-\delta_{\mathrm{N}}\right)+f_{\mathrm{ani}}+2 F(\phi) K_{\mathrm{f}}(1-\chi) \delta_{\mathrm{N}} \\ C_{12}^{\mathrm{sat}}= & C_{13}^{\mathrm{sat}} \approx \lambda_{\mathrm{d}}\left(1-\delta_{\mathrm{N}}\right)+f_{\mathrm{ani}}+ \\ & F(\phi) K_{\mathrm{f}}(1-\chi) \delta_{\mathrm{N}} \\ C_{22}^{\mathrm{sat}}= & C_{33}^{\mathrm{sat}} \approx M_{\mathrm{d}}\left(1-\chi^2 \delta_{\mathrm{N}}\right)+f_{\mathrm{ani}} \\ C_{23}^{\mathrm{sat}} \approx & \lambda_{\mathrm{d}}\left(1-\chi^2 \delta_{\mathrm{N}}\right)+f_{\mathrm{ani}} \\ C_{44}^{\mathrm{sat}}= & \mu \\ C_{55}^{\mathrm{sat}}= & C_{66}^{\mathrm{sat}}=\mu\left(1-\delta_{\mathrm{T}}\right) \end{aligned}\right. $

(7)

为了验证弱各向异性近似刚度参数(式(7))的精度，对比式(7)与精确刚度参数(式(5))的计算结果。假定裂缝型储层模型(岩石由石英和黏土组成，比例为2∶1，含水饱和度分别为30%和70%)裂隙饱含气和水，使用Wood公式计算混合流体模量，并使用VRH方程计算矿物的平均模量。

图 1为不同含气饱和度条件下裂缝型储层模型的刚度参数计算结果。由图可见：①在含气饱和度S_g相同情况下(图 1a与图 1b、图 1c与图 1d)，δ_N越小，则式(7)计算结果越准确，故式(7)更适用于弱各向异性介质；②在δ_N相同情况下(图 1a与图 1c、图 1b与图 1d)，S_g越大，则式(7)计算结果越准确，故式(7)更适用于含气(低体积模量)裂缝型储层。

1.2 裂缝储层“固液缝”解耦地震反射方程推导

式(7)是各向同性介质中含一组定向排列的垂直裂缝形成的等效HTI介质饱和流体弹性刚度参数。在此基础上，考虑反射界面两侧弹性参数弱差值近似及弱各向异性近似(即$\delta_{\mathrm{N}}、\delta_{\mathrm{T}} \ll 1$)，忽略含ΔMδ_N、Δλδ_N、Δf_ani(ϕ)δ_N、ΔK_fδ_N的项(符号Δ表示反射界面两侧背景介质参数的差值)，则HTI等效介质的饱和流体扰动刚度参数可表示为

$ \Delta \boldsymbol{C}_{\mathrm{HTI}}^{\mathrm{sat}}=\left(\begin{array}{cccccc} \Delta C_{11}^{\mathrm{sat}} & \Delta C_{12}^{\mathrm{sat}} & \Delta C_{13}^{\mathrm{sat}} & 0 & 0 & 0 \\ \Delta C_{12}^{\mathrm{sat}} & \Delta C_{22}^{\mathrm{sat}} & \Delta C_{23}^{\mathrm{sat}} & 0 & 0 & 0 \\ \Delta C_{13}^{\mathrm{sat}} & \Delta C_{23}^{\mathrm{sat}} & \Delta C_{33}^{\mathrm{sat}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \Delta C_{44}^{\mathrm{sat}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \Delta C_{55}^{\mathrm{sat}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \Delta C_{66}^{\text {sat }} \end{array}\right) $

(8)

其中

$ \left\{\begin{aligned} \Delta C_{11}^{\mathrm{sat}} \approx & \Delta M_{\mathrm{d}}-M_{\mathrm{d}} \Delta \delta_{\mathrm{N}}+\Delta f_{\text {ani }}+2 F(\phi) K_{\mathrm{f}}(1-\chi) \Delta \delta_{\mathrm{N}} \\ \Delta C_{12}^{\text {sat }}= & \Delta C_{13}^{\text {sat }} \approx \Delta \lambda_{\mathrm{d}}-\lambda_{\mathrm{d}} \Delta \delta_{\mathrm{N}}+\Delta f_{\text {ani }}+ \\ & F(\phi) K_{\mathrm{f}}(1-\chi) \Delta \delta_{\mathrm{N}} \\ \Delta C_{22}^{\text {sat }}= & \Delta C_{33}^{\text {sat }} \approx \Delta M_{\mathrm{d}}-M_{\mathrm{d}} \chi^2 \Delta \delta_{\mathrm{N}}+\Delta f_{\text {ani }} \\ \Delta C_{23}^{\text {sat }} \approx & \lambda_{\mathrm{d}}-M_{\mathrm{d}} \chi^2 \Delta \delta_{\mathrm{N}}+\Delta f_{\text {ani }} \\ \Delta C_{44}^{\text {sat }}= & \Delta \mu \\ \Delta C_{55}^{\text {sat }}= & \Delta C_{66}^{\text {sat }}=\Delta \mu-\mu \Delta \delta_{\mathrm{T}} \end{aligned}\right. $

基于Born积分和稳相方法^[21]，反射系数与散射函数之间的关系可简单表示为

$ R_{\mathrm{PP}}(\theta)=\frac{1}{4 \rho \cos ^2 \theta}\left(\Delta \rho \xi^{\mathrm{PP}}+\sum\limits_{i=1}^6 \sum\limits_{j=1}^6 \Delta C_{i j}^{\mathrm{Cat}} \eta_{\mathrm{pp}}^{\mathrm{pp}}\right) $

(9)

式中：θ为地震波入射角；ρ为介质密度；ΔC_ij^sat为扰动刚度参数；ξ^PP和η_ij^pp可由入射PP波、散射PP波的慢度和极化矢量表示。

将式(8)代入式(9)可以得到(附录A)

$ \begin{aligned} & R_{\mathrm{PP}}(\theta, \varphi)=a(\theta) \frac{\Delta M_{\mathrm{d}}}{M_{\mathrm{d}}}+b(\theta) \frac{\Delta \mu}{\mu}+c(\theta) \frac{\Delta \rho}{\rho}+ \\ & d(\theta) \frac{\Delta f_{\mathrm{ani}}}{f_{\mathrm{ani}}}+e(\theta, \varphi) \Delta \delta_{\mathrm{N}}+f(\theta, \varphi) \Delta \delta_{\mathrm{T}} \end{aligned} $

(10)

其中

$ \begin{aligned} & a(\theta)=\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4} ; b(\theta)=-\frac{2}{\gamma_{\mathrm{sat}}^2} \sin ^2 \theta \\ & c(\theta)=\frac{1}{2}-\frac{\sec ^2 \theta}{4} ; d(\theta)=\left(1-\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2}\right) \frac{\sec ^2 \theta}{4} \\ & e(\theta, \varphi)=-\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4}\left[\frac{2}{\gamma_{\mathrm{sat}}^2}\left(1-\sin ^2 \theta \cos ^2 \varphi\right)-1\right]^2 \\ & f(\theta, \varphi)=\frac{1}{\gamma_{\mathrm{sat}}^2} \cos ^2 \varphi \sin ^2 \theta\left(1-\tan ^2 \theta \sin ^2 \varphi\right) \end{aligned} $

式中：φ为方位角；ΔM_d、Δμ、Δρ、Δf_ani、Δδ_N和Δδ_T分别为介质扰动参数；γ_sat²=M_s/μ，γ_dry²=M_d/μ，M_s为含水岩石纵波模量，通过计算井上数据得到。

图 2为M_d、μ、ρ、f_ani、δ_N及δ_T扰动变化对R_PP的影响，由图可见：f_ani扰动(图 2d)在θ≤40^°范围内对R_PP均有贡献，因此可很好地反演f_ani；M_d(图 2a)和μ(图 2b)扰动几乎在所有θ对R_PP均有贡献；相比而言，ρ扰动(图 2c)在小θ才对R_PP有贡献，而在大θ几乎没有任何贡献；不同θ、φ的δ_N(图 2e)与δ_T(图 2f)扰动对R_PP贡献较小，但对大θ数据较敏感。因此大入射角或远炮检距数据更利于反演裂缝弱度。

1.3 裂缝储层“固液缝”直接油气反演方法

对于五维数据，时间域和频率域地震反演分别考虑不同域中的地震响应与模型合成地震记录的匹配程度，进而搜索拟合数据空间的最优模型参数。时间域地震反演的抗噪性较频率域强，但是时间域地震反演的分辨率不及频率域。为保持频率域地震反演的高分辨率优良特性和解决频率域地震信号对噪声敏感的问题，充分利用带限地震信号的时间域和频率域信息，利用傅里叶变换将地震道由时间域变换到频率域，然后综合时间域和频率域响应特征发展基于贝叶斯框架的时间—频率域联合反演方法。

将式(10)的背景介质参数的相对差值取自然对数近似替换相对差值，即ΔM_d/M_d≈Δ(lnM_d)、Δμ/μ≈Δ(lnμ)、Δρ/ρ≈Δ(lnρ)和Δf_ani/f_ani≈Δ(lnf_ani)，假设背景介质参数与裂缝参数变化连续且扰动范围较小，则Δ(lnM_d)≈d(lnM_d)、Δ(lnμ)≈d(lnμ)和Δ(lnρ)≈d(lnρ)，且Δδ_N≈dδ_N、Δδ_T≈dδ_T。

将式(10)分别与时间域和频率域地震子波褶积，得到时间域正演方程

$ \begin{aligned} & S(t)=W(t) a_{M_{\mathrm{d}}}(\theta) D \ln M_{\mathrm{d}}+W(t) a_\mu(\theta) D \ln \mu+ \\ & W(t) a_\rho(\theta) D \ln \rho+W(t) a_{f_{\mathrm{ani}}}(\theta) D \ln f_{\mathrm{ani}}+ \\ & W(t) a_{\delta_{\mathrm{N}}}(\theta, \varphi) D \delta_{\mathrm{N}}+W(t) a_{\delta_{\mathrm{T}}}(\theta, \varphi) D \delta_{\mathrm{T}} \end{aligned} $

(11)

和频率域正演方程

$ \begin{aligned} & S^{\prime}(\omega)=W^{\prime}(\omega) E a_{M_{\mathrm{d}}}(\theta) D \ln M_{\mathrm{d}}+W^{\prime}(\omega) E a_\mu(\theta) D \ln \mu+ \\ & W^{\prime}(\omega) E a_\rho(\theta) D \ln \rho+W^{\prime}(\omega) E a_{f_{\mathrm{ani}}}(\theta) D \ln f_{\mathrm{ani}}+ \\ & W^{\prime}(\omega) E a_{\delta_{\mathrm{N}}}(\theta, \varphi) D \delta_{\mathrm{N}}+W^{\prime}(\omega) E a_{\delta_{\mathrm{T}}}(\theta, \varphi) D \delta_{\mathrm{T}} \end{aligned} $

(12)

式中：S(t)和S^′(ω)分别为时间域、频率域地震记录，t为时间，ω为角频率；W(t)和W^′(ω)分别为时间域、频率域子波；a_Md(θ)、a_μ(θ)、a_ρ(θ)、a_fani(θ)、a_δN(θ, φ)以及a_δT(θ, φ)分别对应式(10)中的系数项；E为傅里叶算子；D为差分算子。

将式(11)和式(12)写成矩阵形式

$ \left(\begin{array}{c} \boldsymbol{S}(t) \\ \boldsymbol{S}^{\prime}(\omega) \end{array}\right)=\left(\begin{array}{cc} \boldsymbol{G}_{\text {iso }}(t) & \boldsymbol{G}_{\text {ani }}(t) \\ \boldsymbol{G}_{\text {iso }}^{\prime}(\omega) & \boldsymbol{G}_{\text {ani }}^{\prime}(\omega) \end{array}\right)\left(\begin{array}{c} \boldsymbol{m}_{\text {iso }} \\ \boldsymbol{m}_{\text {ani }} \end{array}\right) $

(13)

式中：G_iso (t)、G_iso^′ (ω)分别为时间域、频率域各向同性正演算子；G_ani (t)、G_ani^′ (ω)分别为时间域、频率域各向异性正演算子；m_iso、m_ani分别为各向同性、各向异性模型参数矩阵。

在贝叶斯框架下，基于五维数据时频域分步反演算法，借助目标参数协方差矩阵C_cov的特征值和特征向量U处理G_iso (t)、G_iso^′ (ω)、G_ani (t)、G_ani^′ (ω)及m_iso、m_ani，提高待反演参数稳定性，利用奇异值分解C_cov=UΣU^T得到U和特征值正定矩阵Σ，即

$ \begin{aligned} & \left(\begin{array}{ll} \widetilde{\boldsymbol{G}}_{\text {iso }}(t) & \widetilde{\boldsymbol{G}}_{\mathrm{ani}}(t) \\ \widetilde{\boldsymbol{G}}'_{\mathrm{iso}}(\omega) & \widetilde{\boldsymbol{G}'}_{\mathrm{ani}}(\omega) \end{array}\right)= \\ & \left(\begin{array}{cc} \boldsymbol{G}_{\text {iso }}(t) & \boldsymbol{G}_{\text {ani }}(t) \\ \boldsymbol{G}'_{\text {iso }}(\omega) & \boldsymbol{G}'_{\text {ani }}(\omega) \end{array}\right)\left(\operatorname{kron}\left(\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{-\frac{1}{2}} \boldsymbol{U}^{\mathrm{T}}, \boldsymbol{I}_L\right)\right)^{-1} \\ & \end{aligned} $

(14)

$ \left(\begin{array}{c} \widetilde{\boldsymbol{m}}_{\text {iso }} \\ \widetilde{\boldsymbol{m}}_{\text {ani }} \end{array}\right)=\left(\mathrm{kron}\left(\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{-\frac{1}{2}} \boldsymbol{U}^{\mathrm{T}}, \boldsymbol{I}_L\right)\right)\left(\begin{array}{c} \boldsymbol{m}_{\text {iso }} \\ \boldsymbol{m}_{\text {ani }} \end{array}\right) $

(15)

因此，地震记录可以写为

$ \left(\begin{array}{c}\boldsymbol{S}(t) \\ \boldsymbol{S}^{\prime}(\omega)\end{array}\right)=\left(\begin{array}{cc}\widetilde{\boldsymbol{G}}_{\text {iso }}(t) & \widetilde{\boldsymbol{G}}_{\text {ani }}(t) \\ \widetilde{\boldsymbol{G}}_{\text {iso }}^{\prime}(\omega) & \widetilde{\boldsymbol{G}}'_{\text {ani }}(\omega)\end{array}\right)\left(\begin{array}{l}\widetilde{\boldsymbol{m}}_{\text {iso }} \\ \widetilde{\boldsymbol{m}}_{\text {ani }}\end{array}\right) $

(16)

式中：kron(·)为Kronecker乘积算子；I_L为L阶单位矩阵，L为井曲线采样点数。构建协方差矩阵的方式有多种，本文由统计测井数据获得。在贝叶斯框架下，构建五维数据时频域分步反演目标泛函

$ \begin{aligned} & F(\tilde{\boldsymbol{m}})=(\boldsymbol{S}-\boldsymbol{G} \tilde{\boldsymbol{m}})^{\mathrm{T}}(\boldsymbol{S}-\boldsymbol{G} \tilde{\boldsymbol{m}})+(\boldsymbol{S}-\widetilde{\boldsymbol{G}} \tilde{\boldsymbol{m}})^{\mathrm{T}}(\boldsymbol{S}'-\boldsymbol{G} \tilde{\boldsymbol{m}})+ \\ & 2 \sigma^2 \sum\limits_{i=1}^6 \ln \left(1+\frac{\widetilde{\boldsymbol{m}}_i^2}{\sigma^2}\right)+\sum\limits_{i=1}^6 \lambda_i\left(\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_i-\widetilde{\boldsymbol{P}}_i \widetilde{\boldsymbol{m}}_i\right)^{\mathrm{T}}\left(\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_i-\widetilde{\boldsymbol{P}}_i \widetilde{\boldsymbol{m}}_i\right) \end{aligned} $

(17)

$\begin{aligned} & F(\tilde{\boldsymbol{m}})_{\mathrm{ani}}=\left(\boldsymbol{S}-\boldsymbol{S}_{\mathrm{iso}}-\boldsymbol{G} \widetilde{\boldsymbol{m}}_{\mathrm{ani}}\right)^{\mathrm{T}}\left(\boldsymbol{S}-\boldsymbol{S}_{\mathrm{iso}}-\boldsymbol{G} \tilde{\boldsymbol{m}}_{\mathrm{ani}}\right)+ \\ & \left(\boldsymbol{S}^{\prime}-\boldsymbol{S}_{\mathrm{iso}}^{\prime}-\widetilde{\boldsymbol{G}}' \widetilde{\boldsymbol{m}}_{\mathrm{ani}}\right)^{\mathrm{T}}\left(\boldsymbol{S}-\boldsymbol{S}_{\mathrm{iso}}^{\prime}-\boldsymbol{G}' \tilde{\boldsymbol{m}}_{\mathrm{ani}}\right)+ \\ & 2 \sigma^2 \sum\limits_{i=5}^6 \ln \left(1+\frac{\tilde{\boldsymbol{m}}_i^2}{\sigma^2}\right)+\sum\limits_{i=5}^6 \lambda_i\left(\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_i-\boldsymbol{P}_i \widetilde{\boldsymbol{m}}_i\right)^{\mathrm{T}}\left(\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_i-\boldsymbol{P}_i \widetilde{\boldsymbol{m}}_i\right) \end{aligned} $

(18)

式中：λ_i为模型参数的约束系数；P_i为积分算子；Λ_i=ln(R_i/R₀)/2，R₀为模型参数的初始值；σ²为背景噪声分布方差。反演过程首先将M_d、μ、ρ和f_ani四个参数作为第一步的重点优化参数，然后在此基础上，将δ_N和δ_T两个裂缝参数作为第二步的优化目标，从而最终求取全部待反演参数。

2 模型试算及实际资料处理 2.1 模型试算

选用A区的一口井的合成方位地震角度道集验证五维数据时频联合反演模型参数的可行性。首先，基于褶积模型及提取的地震子波，并添加高斯随机噪声，合成无噪及信噪比为5∶1、2∶1的不同方位角的地震角度道集(图 3、图 5、图 7)及其反演结果(图 4、图 6、图 8)。由图 4、图 6、图 8可见，模型参数反演结果与实际测井数据吻合较好(平均相对误差较小)，验证了五维数据时频联合反演的可行性。

2.2 实际资料处理

B区地层发育大量高角度裂缝，预测目标为裂缝型非常规页岩油气储层，可将其等效为HTI介质。选取45°、135°方位五维数据，每个方位的叠前炮检距数据均经过AVO保幅处理、从炮检距到角度道集的转换。图 9为45°、135°方位角地震数据小、中、大入射角部分角度叠加剖面，图 10与图 11分别为M_d、μ、ρ、f_ani、δ_N及δ_T初始模型剖面与反演剖面。由图 11可见：f_ani反演剖面的结构与图 9相似，并与测井解释结果一致，A井作为测试井并未参与初始模型的建立以及后期反演，f_ani反演剖面的低值异常对应测井解释的油气位置，f_ani解耦流体识别结果的油气界面清晰(红色区域)，且与钻井结果吻合，进一步验证了方法的有效性与可靠性。模型参数反演结果与原始测井数据的对比结果表明(图 12)，两者吻合较好，进一步验证了方法的有效性。

3 结论

本文方法旨在从五维数据中直接获取裂缝型储层参数。结合各向异性Gassmann孔隙弹性理论和线性滑动理论，以“固液缝”解耦流体因子为基础，推导了由“固液缝”解耦流体因子直接表征的新的纵波反射系数方程，消除了岩石骨架、裂缝、孔隙及油气水耦合对地震反射特征的影响；分析、讨论了新纵波反射系数方程的适用条件及敏感性，在此基础上，构建贝叶斯框架的地震反演目标函数，并基于时频资料研究五维数据流体识别方法，可直接获得裂缝型储层“固液缝”解耦流体因子，消除了累计误差的影响，提高了预测精度。模型试算及实际资料处理结果表明，所提方法可以合理、可靠地预测裂缝型储层油气分布。

附录A 式(10)的推导过程

将式(8)的扰度刚度矩阵代入式(9)，得

$ \begin{aligned} & R_{\mathrm{PP}}(\theta, \varphi)= \\ & \frac{\Delta M_{\mathrm{d}}+\Delta f_{\mathrm{ani}}-4 \Delta \delta_{\mathrm{T}} \mu \cos ^2 \varphi \sin ^2 \theta\left(\sin ^2 \theta \sin ^2 \varphi-\cos ^2 \theta\right)}{4 \rho V_{\mathrm{P}}^2 \cos ^2 \theta}- \\ & \frac{\Delta \delta_{\mathrm{N}} M_{\mathrm{d}}\left(1-\frac{2}{\gamma_{\mathrm{dry}}^2} \sin ^2 \theta \sin ^2 \varphi-\frac{2}{\gamma_{\mathrm{dry}}^2} \cos ^2 \theta\right)^2}{4 \rho V_{\mathrm{P}}^2 \cos ^2 \theta}+ \\ & \frac{\Delta \rho\left(2 \cos ^2 \theta-1\right)}{4 \rho \cos ^2 \theta}-\frac{2 \Delta \mu \sin ^2 \theta}{\rho V_{\mathrm{P}}^2} \end{aligned} $

(A-1)

整理可得

$ \begin{aligned} & R_{\mathrm{PP}}(\theta, \varphi)= \\ & \left(\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4}\right) \frac{\Delta M_{\mathrm{d}}}{M_{\mathrm{d}}}-\left(\frac{2}{\gamma_{\mathrm{sat}}^2} \sin ^2 \theta\right) \frac{\Delta \mu}{\mu}+ \\ & \left(\frac{1}{2}-\frac{\sec ^2 \theta}{4}\right) \frac{\Delta O}{\rho}+\frac{K_{\mathrm{f}} G(\phi)+2 F(\phi) K_{\mathrm{f}} \chi \delta_{\mathrm{N}}}{M_{\rm{s}}} \frac{\sec (\theta)^2}{4} \frac{\Delta f_{\mathrm{ani}}}{f_{\mathrm{ani}}}- \\ & \frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4}\left(1+\frac{2}{\gamma_{\mathrm{dry}}^2} \sin ^2 \theta \sin ^2 \varphi-\frac{2}{\gamma_{\mathrm{dry}}^2} \cos ^2 \theta\right)^2 \Delta \delta_{\mathrm{N}}+ \\ & \frac{1}{\gamma_{\mathrm{sat}}^2} \cos ^2 \varphi \sin ^2 \theta\left(1-\tan ^2 \theta \sin ^2 \varphi\right) \Delta \delta_{\mathrm{T}} \end{aligned} $

(A-2)

根据各向异性Gassmann流体项F(ϕ)K_f=M_s－M_d，而且$F(\phi) K_{\mathrm{f}} \Delta \delta_{\mathrm{N}} \ll M_{\mathrm{s}}$，因此

$ \begin{aligned} & R_{\mathrm{PP}}(\theta, \varphi)=a(\theta) \frac{\Delta M_{\mathrm{d}}}{M_{\mathrm{d}}}+b(\theta) \frac{\Delta \mu}{\mu}+c(\theta) \frac{\Delta \rho}{\rho}+ \\ & d(\theta) \frac{\Delta f_{\mathrm{ani}}}{f_{\mathrm{ani}}}+e(\theta, \varphi) \Delta \delta_{\mathrm{N}}+f(\theta, \varphi) \Delta \delta_{\mathrm{T}} \end{aligned} $

(A-3)

其中

$ \begin{aligned} & a(\theta)=\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4} ; b(\theta)=-\frac{2}{\gamma_{\mathrm{sat}}^2} \sin ^2 \theta \\ & c(\theta)=\frac{1}{2}-\frac{\sec ^2 \theta}{4} ; d(\theta)=\left(1-\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2}\right) \frac{\sec ^2 \theta}{4} \\ & e(\theta, \varphi)=-\frac{\gamma_{\mathrm{dry}}^2}{\gamma_{\mathrm{sat}}^2} \frac{\sec ^2 \theta}{4}\left[\frac{2}{\gamma_{\mathrm{sat}}^2}\left(1-\sin ^2 \theta \cos ^2 \varphi\right)-1\right]^2 \\ & f(\theta, \varphi)=\frac{1}{\gamma_{\mathrm{sat}}^2} \cos ^2 \varphi \sin ^2 \theta\left(1-\tan ^2 \theta \sin ^2 \varphi\right) \end{aligned} $