0 引言

流体(含气性)预测对于裂缝性油气藏的勘探与开发具有重要意义。在成岩作用和后期构造应力作用下，储层常常存在一组或多组定向排列的裂缝^[1]。当地震波穿过裂缝时，一方面定向排列的裂缝会造成地震波速度的各向异性^[2]；另一方面，地震波对裂缝的挤压作用会引起流体流动，该过程伴随着地震波的频散和衰减^[3]。因此，裂缝流体的预测在考虑地震波形的振幅、各向异性等信息的同时，应考虑地震波的频散特征。

在裂缝型储层表征方面，振幅与角度和方位角(AVAZ)反演被广泛用于纵波地震数据的各向异性反演^[4-5]。具有垂直对称轴的横向各向同性(Transverse Isotropic with a Vertical Symmetry Axis, VTI)背景介质中嵌入一组平行排列的垂直裂缝可等效为正交各向异性介质^[6]。为提取正交各向异性介质中的各向异性特征，Bachrach根据Pšenčik等^[7]提出的公式和符号研究了基于正交各向异性OA(Orthorhombic Anisotropy, OA) 介质的线性AVAZ反演^[8]。潘新朋等^[9]研究了贝叶斯框架下基于OA介质的AVAZ反演方法。

裂缝流体预测通常依赖于与流体相关的异常特征的反演，这些异常特征被称为流体因子^[10-11]。流体识别的效率主要取决于两个方面：流体因子的敏感性和地震反演的可靠性^[12]。对于与频率无关的流体因子，应用最广的是以Zoeppritz方程为基础推导的弹性参数和AVO属性的组合^[13-16]。这种传统的流体识别因子在实际应用中取得了较好的效果，但未考虑到地震频散和衰减的影响，导致其分辨率和准确性较低。随着时频分析技术的进步，发展了与频率相关的流体因子^[17]。Taner等^[18]最先注意到含气层下方的“低频阴影”现象，并推测其由地震波衰减引起。但后续的研究表明，造成“低频阴影”的因素不一定与衰减有关^[19]。

在地震波频散和衰减的机制方面，Chapman等^[20]、张金伟等^[21]、Shuai等^[22]提出动态等效介质理论，从地震波诱导流体流动的角度解释了地震波的频散和衰减，即反射系数不仅与入射角和方位角有关，而且还随频率变化。该理论为FDAVO(Frequency-Dependent Amplitude Versus Offset)反演提供了理论基础。

在频散属性反演方面，对于各向同性介质，Wilson等^[23]基于Chapman等^[20]提出的岩石物理理论，结合时频分析方法和AVO技术初步开发了依赖频率的AVO反演框架，计算反映地震波频散程度的属性并将其作为流体因子。吴小羊^[24]将改进的Wigner-Ville分布(SPWVD)方法用于FDAVO反演。程冰洁等^[25]基于Smith-Gidlow近似式推导了截距、梯度和拟泊松比等频变AVO属性。罗鑫等^[26]根据f-μ-ρ反射系数公式^[16]进行了FDAVO反演。不过，实际储层中的孔隙空间具有多样性和复杂性，通过研究各向异性参数随频率的变化可进一步识别有效裂缝(流体饱和)。对于各向异性介质，Zhang等^[27]通过计算不同方位反射系数的差值以消除反射系数公式中的各向同性项，推导了表征各向异性参数频散程度的频散项。Ajaz等^[28]基于Rüger方程，推导出频变AVAZ反演的频域解析表达式，并在碳酸盐岩应用中进行了验证。上述研究表明，考虑地震波频散有利于裂缝性储层中流体信息的获取。

由于正交各向异性介质中地震响应的复杂性，准确刻画其中的流体分布特征仍然具有一定的挑战性。为此，本文提出了基于OA介质频率相关AVAZ(OA-FDAVAZ)反演的含气性预测方法。首先通过Chapman多组裂缝模型分析弹性参数的频变特性以及孔隙度、裂缝密度等参数对频散特征的影响；其次将OA介质纵波反射系数近似式拓展到频率域，构建OA-FDAVAZ反演的目标方程，给出频散组合因子的计算公式，并提出相应的反演流程；最后利用最小二乘法求解正交各向异性频散属性。模型数据和实际数据的计算结果均验证了方法的有效性和稳定性。

1 Chapman动态等效介质理论

考虑到实际储层通常发育多组裂缝，Chapman^[29]将单裂缝模型扩展到包含两组不同方向、大小和连通性的裂缝集。双裂缝模型证明了由两个特征频率(速度随频率变化最大时对应的频率)控制了频率相关各向异性。这两个特征频率可由流体流动性和裂缝的尺度定义。包含两组裂缝的频变刚度张量为

$ \mathit{\boldsymbol{C}}={\mathit{\boldsymbol{C}}}_{0}-{\phi }_{\mathrm{p}}{\mathit{\boldsymbol{C}}}_{{\mathrm{p}}}-{\varepsilon }_{1}{\mathit{\boldsymbol{C}}}_{1}-{\varepsilon }_{2}{\mathit{\boldsymbol{C}}}_{2} $

(1)

式中：$ {\phi }_{\mathrm{p}} $为孔隙度；$ {\varepsilon }_{1} $和$ {\varepsilon }_{2} $分别是两组裂缝的密度；$ {\mathit{\boldsymbol{C}}}_{0} $表示背景介质的各向同性弹性张量；$ {\mathit{\boldsymbol{C}}}_{{\mathrm{p}}} $是球形孔隙的校正张量；$ {\mathit{\boldsymbol{C}}}_{1} $和$ {\mathit{\boldsymbol{C}}}_{2} $分别为两组裂缝的校正张量。

在Chapman模型中，最重要的两个参数分别是两种尺度下的弛豫时间。微观尺度下的弛豫时间($ {\tau }_{\mathrm{m}}) $由孔隙与微裂隙之间、不同方向微裂隙之间的流体流动引起，对应着传统的喷射流频率；介观尺度下的弛豫时间($ {\tau }_{\mathrm{f}} $)由裂缝与孔隙、微裂隙之间的流体流动引起，对应着较低的特征频率。两种弛豫时间有如下关系

$ {\tau }_{\mathrm{f}}=\frac{r}{\varsigma }{\tau }_{\mathrm{m}} $

(2)

式中：$ r $为裂缝半径；$ \varsigma $为岩石颗粒尺度(微观尺度)。

在长波长假设条件下，周期性薄互层中发育一组垂直裂缝可等效为OA介质(图 1)。为了解OA介质中弹性参数的频变特性和流体敏感性，通过Chapman模型进行数值模拟。如图 2所示，设置各向同性背景介质中的拉梅系数$ \lambda =31\ \mathrm{G}\mathrm{P}\mathrm{a} $、$ \mu =13\ \mathrm{G}\mathrm{P}\mathrm{a} $，骨架密度为$ 2.5\ \mathrm{g}/\mathrm{c}{\mathrm{m}}^{3} $，孔隙度为0.1。为了模拟正交特性，将两组裂缝的极化角分别设置为$ 0 $°和$ 90 $°，方位角保持一致(均为0°)。此外，两组裂缝密度均为0.02，平均裂缝半径为1 m。为区分不同流体饱和下的频散程度，考虑孔隙空间中分别充填水和气两种情况。假设水饱和时的流体体积模量为3.2 GPa，流体密度为$ 1.1\ \mathrm{g}/\mathrm{c}{\mathrm{m}}^{3} $，$ {\tau }_{\mathrm{m}}=3\times {10}^{-5}\mathrm{s} $；气饱和时的流体体积模量为0.024 GPa，流体密度为$ 65\times {10}^{-3}\mathrm{g}/\mathrm{c}{\mathrm{m}}^{3} $，$ {\tau }_{\mathrm{m}}=9\times {10}^{-7}\mathrm{s} $。

通过计算式(1)得到矩阵张量$ \mathit{\boldsymbol{C}} $后，求解Christoffel方程^[30]可得到具有频率依赖性的复相速度$ V\left(f\right) $。根据以下公式可进一步得到实相速度$ {V}_{\mathrm{P}}\left(f\right) $和逆品质因子$ {Q}_{\mathrm{P}}^{-1}\left(f\right) $

$ {V}_{\mathrm{P}}\left(f\right)={\left[\mathrm{R}\mathrm{e}\left(\frac{1}{V\left(f\right)}\right)\right]}^{-1} $

(3)

$ {Q}_{\mathrm{P}}^{-1}\left(f\right)=\frac{\mathrm{I}\mathrm{m}\left[V{\left(f\right)}^{2}\right]}{\mathrm{R}\mathrm{e}\left[V{\left(f\right)}^{2}\right]} $

(4)

式中：Im表示取复数的虚部；Re表示取复数的实部；f为频率。

各向异性参数$ \delta $的分量表示为

$ {\delta }_{x}=\frac{{\left({c}_{13}+{c}_{55}\right)}^{2}-{\left({c}_{33}-{c}_{55}\right)}^{2}}{2{c}_{33}\left({c}_{33}-{c}_{55}\right)} $

(5)

$ {\delta }_{y}=\frac{{\left({c}_{23}+{c}_{44}\right)}^{2}-{\left({c}_{33}-{c}_{44}\right)}^{2}}{2{c}_{33}\left({c}_{33}-{c}_{44}\right)} $

(6)

式中：$ {\delta }_{x} $为$ \delta $在对称面$ \left[{x}_{1}, {x}_{3}\right] $内的分量；$ {\delta }_{y} $为$ \delta $在对称面$ \left[{x}_{2}, {x}_{3}\right] $内的分量；$ {c}_{ij} $(i、j=1, 2, …, 5)为对应的矩阵张量元素。

计算不同含气饱和度$ {S}_{\mathrm{g}} $下纵波速度$ {V}_{\mathrm{P}}\left(f\right) $ (图 3a)和逆品质因子$ {Q}_{\mathrm{P}}^{-1}\left(f\right) $ (图 3b)随频率的变化(地震波入射角为0°)。结果表明：纵波速度和逆品质因子在地震频带范围(1~100 Hz)内随频率变化明显，不同含气饱和度下频散和衰减的差异较大；随含气饱和度的增加，纵波速度降低，逆品质因子和特征频率均增大。计算图 3a中纵波速度关于频率的导数如图 3c所示，可见：随含气饱和度的增加，纵波速度的频散程度增大。

此外，两个各向异性参数分量$ {\delta }_{x} $和$ {\delta }_{y} $随频率的变化(图 4)和其关于频率的导数的绝对值随频率变化的规律(图 5)均与纵波速度类似。因此，可利用与含气饱和度有关的各向异性频散和衰减信息预测有效裂缝。

为进一步探究孔隙和裂缝发育程度对地震波频散的影响，在其他参数保持一致的情况下，选择频散现象较明显的饱和气模型，计算地震频带内纵波速度随频率、孔隙度$ {\phi }_{\mathrm{p}} $和裂缝密度$ \varepsilon $(两组裂缝密度相同)的变化规律，并定义地震频带范围内的速度变化量$ \mathrm{\Delta }M={V}_{\mathrm{P}}\left({f}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)-{V}_{\mathrm{P}}\left({f}_{\mathrm{m}\mathrm{i}\mathrm{n}}\right) $以衡量纵波速度的频散程度，其中$ {f}_{\mathrm{m}\mathrm{a}\mathrm{x}}=100\ \mathrm{H}\mathrm{z} $，$ {f}_{\mathrm{m}\mathrm{i}\mathrm{n}}=1\ \mathrm{H}\mathrm{z} $。

由图 6a可知，纵波速度随孔隙度和频率而变化，且频散程度随孔隙度的增加先急剧增大、随后增速放缓(图 6b)，这表明孔隙度对地震波频散有一定影响，但这种影响难以反映孔隙度的变化。相比之下，纵波速度在低频时随裂缝密度变化更明显(图 7a)，频散程度与裂缝密度近似呈线性正相关(图 7b)，因此可以利用频散程度表征裂缝发育情况。

反演各向异性参数时，需要用到OVT域叠前地震数据，其中包含地震波入射角($ \theta ) $和方位角($ \varphi ) $信息，因此有必要研究不同入射角和方位角对弹性参数频散程度的影响。在前述模型的基础上分别计算频散速度随入射角和方位角的变化关系(图 8~图 9)。结果表明，在正交各向异性介质中，纵波速度同时随频率和地震波入射角而变化；在低频范围内，纵波速度随入射角的变化更明显(图 8a)。纵波速度频散程度随入射角的增加而减小(图 8b)，因此各向异性频散属性反演时选择的入射角不宜过大。

纵波速度同时随频率和方位角而变化(图 9a)，纵波速度随方位角呈现周期性变化，且在低频时这种变化更剧烈。纵波频散程度分别在方位角为$ 90 $°、$ 270 $°和$ 360 $°时取得最大值，因为此时地震波传播方向与垂直裂缝的法向平行(图 9b)。这说明反演的方位角范围只需选择在一个周期内(如$ 0°~180 $°)的数据即可提高反演效率。

2 OA-FDAVAZ反演方法 2.1 原理

2001年，Pšenčik等^[7]推导的任意各向异性和对称轴背景下的正交各向异性介质纵波反射系数近似方程为

$ R\left(\theta , \varphi \right)={R}^{\mathrm{i}\mathrm{s}\mathrm{o}}\left(\theta \right)+{R}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta , \varphi -{\varphi }_{\mathrm{s}}\right) $

(7)

$ {R}^{\mathrm{i}\mathrm{s}\mathrm{o}}\left(\theta \right)=P+G\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta $

(8)

$ \begin{array}{l}{R}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta , \varphi \right)=\left[{{\mathit{\Gamma}} }_{x}\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)+\right.\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.{{\mathit{\Gamma}} }_{y}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\right]\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \end{array} $

(9)

上述式中：$ \theta $为纵波入射角；$ \varphi $为OVT道集在观测坐标系下的方位角；$ {\varphi }_{\mathrm{s}} $是裂缝对称轴在观测坐标系下的方位角；上标iso和ani分别表示各向同性和各向异性介质；截距$ P=0.5\mathrm{\Delta }{I}_{\mathrm{P}}/\overline{{I}_{\mathrm{P}}} $，$ \mathrm{\Delta }\left(\cdot \right) $表示上、下层介质参数的差值，一般为下层参数减去上层参数，$ \overline{\left(\cdot \right)} $表示上、下层介质参数的平均值，$ {I}_{\mathrm{P}}=\rho {V}_{\mathrm{P}} $，$ \rho $为密度，$ {V}_{\mathrm{P}} $为背景介质的纵波速度；梯度$ G=0.5\left[\mathrm{\Delta }{V}_{\mathrm{P}}/\overline{{V}_{\mathrm{P}}}-{\left(2\overline{{V}_{\mathrm{S}}}/\overline{{V}_{\mathrm{P}}}\right)}^{2}\mathrm{\Delta }\left(\rho {V}_{\mathrm{S}}^{2}\right)/\left(\overline{\rho }\overline{{V}_{\mathrm{S}}^{2}}\right)\right] $，$ {V}_{\mathrm{S}} $为背景介质的横波速度；$ {{\mathit{\Gamma}} }_{x} $为各向异性参数在对称面$ \left[{x}_{1}, {x}_{3}\right] $内的分量；$ {{\mathit{\Gamma}} }_{y} $为各向异性参数在对称面$ \left[{x}_{2}, {x}_{3}\right] $内的分量。各向异性部分的参数定义为

$ \left\{\begin{array}{l}{{\mathit{\Gamma}} }_{x}=0.5\left(\mathrm{\Delta }{\delta }_{x}-8{\overline{{V}_{\mathrm{S}}}}^{2}\mathrm{\Delta }{\gamma }_{x}/{\overline{{V}_{\mathrm{P}}}}^{2}\right)\\ {{\mathit{\Gamma}} }_{y}=0.5\left(\mathrm{\Delta }{\delta }_{y}-8{\overline{{V}_{\mathrm{S}}}}^{2}\mathrm{\Delta }{\gamma }_{y}/{\overline{{V}_{\mathrm{P}}}}^{2}\right)\\ {\delta }_{x}=\left({A}_{13}+2{A}_{55}-{V}_{\mathrm{P}}^{2}\right)/{V}_{\mathrm{P}}^{2}\\ {\delta }_{y}=\left({A}_{23}+2{A}_{44}-{V}_{\mathrm{P}}^{2}\right)/{V}_{\mathrm{P}}^{2}\\ {\gamma }_{x}=\left({A}_{55}-{V}_{\mathrm{S}}^{2}\right)/\left(2{V}_{\mathrm{S}}^{2}\right.)\\ {\gamma }_{y}=\left({A}_{44}-{V}_{\mathrm{S}}^{2}\right)/\left(2{V}_{\mathrm{S}}^{2}\right.)\end{array}\right. $

(10)

式中$ {A}_{ij}={C}_{ij}/\rho $(i, j=1，2，…，5)是密度归一化的正交弹性刚度张量的元素。

Chapman^[31]认为纵波速度、横波速度及各向异性参数均与频率有关，当界面两侧频散性质存在差异时，反射系数会随频率而变化，即反射系数可视为入射角、方位角和频率的函数。因此，假设式(7)~式(9)中的$ R\left(\theta , \varphi \right) $、$ P $、$ G $、$ {{\mathit{\Gamma}} }_{x} $和$ {{\mathit{\Gamma}} }_{y} $是频率的函数，则可以将方程(7)拓展到频率域，即

$ \begin{array}{l}R\left(\theta , \varphi , f\right)=P\left(f\right)+G\left(f\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\mathit{\Gamma}} }_{x}\left(f\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\mathit{\Gamma}} }_{y}\left(f\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \end{array} $

(11)

将式(11)关于某一参考频率$ {f}_{0} $对$ P $、$ G $、$ {{\mathit{\Gamma}} }_{x} $和$ {{\mathit{\Gamma}} }_{y} $进行泰勒级数展开，并省略高阶项，整理得到

$ \begin{array}{l}R\left(\theta , \varphi , f\right)\approx R\left(\theta , \varphi , {f}_{0}\right)+\left(f-{f}_{0}\right)\left[{I}_{\mathrm{a}}+\right.\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {I}_{\mathrm{b}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +{I}_{\mathrm{c}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.{I}_{\mathrm{d}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right]\end{array} $

(12)

式中：$ {I}_{\mathrm{a}} $为纵波截距频散属性；$ {I}_{\mathrm{b}} $为各向同性梯度频散属性；$ {I}_{\mathrm{c}} $和$ {I}_{\mathrm{d}} $分别是不同对称面内的各向异性梯度项$ {{\mathit{\Gamma}} }_{x} $和$ {{\mathit{\Gamma}} }_{y} $关于频率的导数，代表了各向异性参数的频散程度，统称为各向异性梯度频散项。则有

$ \begin{array}{l}R\left(\theta , \varphi , {f}_{0}\right)\approx P\left({f}_{0}\right)+G\left({f}_{0}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\mathit{\Gamma}} }_{x}\left({f}_{0}\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{\mathit{\Gamma}} }_{y}\left({f}_{0}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \end{array} $

(13)

$ \left\{\begin{array}{l}{I}_{\mathrm{a}}=\frac{\mathrm{d}}{\mathrm{d}f}\left[P\left(f\right)\right]\\ {I}_{\mathrm{b}}=\frac{\mathrm{d}}{\mathrm{d}f}\left[G\left(f\right)\right]\\ {I}_{\mathrm{c}}=\frac{\mathrm{d}}{\mathrm{d}f}\left[{{\mathit{\Gamma}} }_{x}\left(f\right)\right]\\ {I}_{\mathrm{d}}=\frac{\mathrm{d}}{\mathrm{d}f}\left[{{\mathit{\Gamma}} }_{y}\left(f\right)\right]\end{array}\right. $

(14)

式(13)可进一步简化为

$ \begin{array}{l}\mathrm{\Delta }R\left(\theta , \varphi , f\right)={I}_{\mathrm{a}}\left(f-{f}_{0}\right)+{I}_{\mathrm{b}}\left(f-{f}_{0}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {I}_{\mathrm{c}}\left(f-{f}_{0}\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta +\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {I}_{\mathrm{d}}\left(f-{f}_{0}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left(\varphi -{\varphi }_{\mathrm{s}}\right)\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \end{array} $

(15)

式中$ \mathrm{\Delta }R\left(\theta , \varphi , f\right)=R\left(\theta , \varphi , f\right)-R\left(\theta , \varphi , {f}_{0}\right) $。式(15)是反演的目标方程，包含了参考频率$ {f}_{0} $和各频率分量$ {f}_{i}\left(i=\mathrm{1, 2}, \cdots , K\right) $等，K为频率分量的总个数。其中，方位角和入射角为已知量，而$ {I}_{\mathrm{a}} $、$ {I}_{\mathrm{b}} $、$ {I}_{\mathrm{c}} $、$ {I}_{\mathrm{d}} $是待反演参数。$ {f}_{0} $和$ {f}_{i} $可以通过地震数据频谱分解得到。

2.2 技术流程

OA-FDAVAZ反演技术流程(图 10)如下。

(1) 获取OVT域宽方位道集，将OVT道集划分为N个方位角扇区的角度道集，并抽取各方位的部分叠加角道集$ \mathit{\boldsymbol{S}}\left(t, \theta , \varphi \right) $，其中t为时间。

(2) 通过常规AVAZ反演得到裂缝方位角$ {\varphi }_{\mathrm{s}} $。

(3) 对叠前地震记录进行频谱分析，获取地震记录的主频作为参考频率$ {f}_{0} $。

(4) 采用平滑伪Wigner-Ville分布(SPWVD)将$ \mathit{\boldsymbol{S}}\left(t, \theta , \varphi \right) $转换为一系列不同频率的振幅谱$ \mathit{\boldsymbol{S}}\left(t, \theta , \varphi , f\right) $，其中包括参考频率$ {f}_{0} $，即

$ \mathit{\boldsymbol{S}}\left(t, \theta , \varphi \right)\to \mathit{\boldsymbol{S}}\left(t, \theta , \varphi , f\right) $

(16)

(5) 由于地震子波的影响，不同频率之间的能量分布并不均衡，此现象称为“子波叠印”。为了消除“子波叠印”，设计一种权重函数消除子波的叠加，以平衡不同频率分量的振幅谱，也就是谱均衡，即

$ \mathit{\boldsymbol{r}}\left(t, \theta , \varphi , f\right)=\mathit{\boldsymbol{S}}\left(t, \theta , \varphi , f\right)\mathit{\boldsymbol{W}}\left(\theta , \varphi , f\right) $

(17)

式中：$ \mathit{\boldsymbol{r}}\left(t, \theta , \varphi , f\right) $是反射系数的振幅谱；$ \mathit{\boldsymbol{W}}\left(\theta , \varphi , f\right) $是谱均衡的权重函数，其计算公式为

$ \mathit{\boldsymbol{W}}\left(\theta , \varphi , f\right)=\frac{\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathit{\boldsymbol{S}}\left(\theta , \varphi , {f}_{0}\right)\right]}{\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathit{\boldsymbol{S}}\left(\theta , \varphi , f\right)\right]} $

(18)

式中：$ \mathrm{m}\mathrm{a}\mathrm{x}\left[\mathit{\boldsymbol{S}}\left(\theta , \varphi , {f}_{0}\right)\right] $是参考频率$ {f}_{0} $下振幅谱的最大值；$ \mathrm{m}\mathrm{a}\mathrm{x}\left[\mathit{\boldsymbol{S}}\left(\theta , \varphi , f\right)\right] $为分解频率$ f $下振幅谱的最大值。

(6) 计算不同频率反射系数振幅谱的差值$ \mathrm{\Delta }R\left(\theta , \varphi , f\right) $，然后将$ {\varphi }_{\mathrm{s}} $、$ {f}_{0} $和$ \mathrm{\Delta }R\left(\theta , \varphi , f\right) $代入式(15)，利用最小二乘法求解得到各向异性梯度频散属性$ {I}_{\mathrm{c}} $和$ {I}_{\mathrm{d}} $。

(7) 计算频散组合因子。为进一步降低反演的多解性，降低背景干扰，借鉴Russell流体识别因子的思路^[16]，频散属性的二次及二次以上的幂量纲组合形式有助于突出由孔隙流体导致的储层频散异常。因此，本文定义$ {I}_{\mathrm{c}} $和$ {I}_{\mathrm{d}} $的组合因子$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $用于裂缝气藏含气性预测，即

$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}}^{}={I}_{\mathrm{c}}^{2}+{I}_{\mathrm{d}}^{2} $

(19)

3 模型测试

设计三层理论模型(图 11a)测试本文方法的可靠性。第一层为各向同性介质。第二层分为三部分，由左往右依次为各向同性层、饱和水裂缝层和饱和气裂缝层，其中裂缝层的参数设置见图 2模型。第三层的性质与第一层相同。具体参数如表 1所示。

基于Chapman模型和反射率法计算一系列不同频率、方位角和入射角的反射系数后，与Ricker子波褶积合成叠前方位地震数据(方位角为$ 15 $°$ \sim 165 $°，间隔为$ 30 $°；入射角为$ 5^\circ \sim 40^\circ $，间隔为$ 5^\circ $)。图 11b为方位角$ 15^\circ $、入射角$ 30^\circ $时模型的叠前方位地震数据，其中裂缝层部分出现“低频阴影”现象，且在饱和气的部分更为明显(图 11b中红框部分)。图 12展示了含气模型的合成方位角道集，地震道在不同方位、不同角度的差异明显。

对上述数据进行OA-FDAVAZ反演，其中分解频率分别为10、20、30、40、50 Hz。由于所选Ricker子波的主频为35 Hz，故选35 Hz作为反演的参考频率，得到模型的正交各向异性频散属性$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $(图 13)。由图可见，各向同性地层不包含裂缝，各向异性频散程度最弱。中间地层饱和水的裂缝层发育两组正交裂缝，但其中流体流动较慢，弛豫时间较长，因此在地震频带的频散程度较低；而饱和气的裂缝层弛豫时间较短，在地震频带的频散程度较高。由此可以验证OA-FDAVAZ反演方法预测裂缝含气性具有可靠性。

4 实际资料应用 4.1 研究区概况

选择四川盆地东部五百梯气田实际资料测试本文方法效果。研究区目标储层为志留系龙马溪组($ {\mathrm{S}}_{1}\mathrm{l} $)超深层页岩。$ {\mathrm{S}}_{1}\mathrm{l} $发育陆棚相沉积，顶部由泥岩和粉砂质泥岩组成，伽马(GR)值较低，纵波速度相对较高；底部由富含有机质、硅质的黑色页岩组成，伽马值较高，速度和密度较低，为优质烃源岩。在上覆层的压实作用和成岩胶结作用下，$ {\mathrm{S}}_{1}\mathrm{l} $页岩基质的孔隙度较低(5%左右)。裂缝为页岩气提供了聚集空间，因此可采用OA-FDAVAZ反演方法预测裂缝储层中的含气性。

4.2 单道数据测试

采用研究区A井旁的OVT域叠前地震数据进行测试。通过分析各方位覆盖次数和炮检距范围，将处理好的OVT道集划分为6个方位角扇形区($ 0^\circ \sim 30 ^\circ$、$ 30^\circ \sim 60^\circ $、$ 60^\circ \sim 90^\circ $、$ 90^\circ \sim 120^\circ $、$ 120^\circ \sim 150^\circ $和$ 150^\circ \sim 180^\circ $)，然后在每个方位角扇区内对三个入射角范围($ 6^\circ \sim 14^\circ $、$ 14^\circ \sim 22^\circ $和$ 22^\circ \sim 30^\circ $)内的数据进行部分叠加。图 14a为经过方位时差校正的各方位地震叠加道集。当入射角相同时，振幅随方位角变化显著，表明该区域方位各向异性较强。图 14b是方位角固定为$ 15^\circ $，入射角依次为$ 10^\circ $、$ 18^\circ $和$ 26^\circ $的地震记录, 通过SPWVD方法得到的时频分析结果。对比不同入射角的地震道能量，可以看出最大值均位于26 Hz附近。据此确定井旁OVT道集的主频为26 Hz，选择该频率作为OA-FDAVAZ反演的参考频率。

应用本文方法得到各向异性频散属性，并与对应的测井曲线进行对比(图 15)。由图可见，储层(红色阴影部分)自然伽马为高值异常，而密度和纵波速度均为低值异常，与页岩储层特征一致；同时，正交各向异性频散属性$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $曲线在对应位置也出现峰值，与测井曲线储层特征具有较好的相关性。上述说明了本文方法提出的$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $属性对裂缝含气性的敏感性较强。OA-FDAVAZ反演方法结合了地震数据中的方位与频率信息，考虑了流体在频域引起的异常和裂缝引起的各向异性特征，可以对高含气饱和度的裂缝性储层进行有效识别。

4.3 页岩储层含气性识别

研究区实际地震数据信噪比较高，横向连续性强，垂向分辨率较高且具有足够的方位覆盖率。图 16为不同方位的大、中、小入射角部分叠加剖面，剖面反演的道集范围、反演参数与单道测试一致。

研究区内A井是一口高产井，日产气量为$ 11.03\times {10}^{4}{\mathrm{m}}^{3} $；D井的产气量较低，日产气量仅为$ 2.45\times {10}^{4}{\mathrm{m}}^{3} $。过井叠后剖面如图 17a所示。

由图 17可见，A井处目标储层的正交各向异性频散幅值$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $较高(图 17b，红色椭圆内)，表明其为高含气区域；D井处正交各向异性频散属性响应较弱(图 17b，白色椭圆内)。本文方法反演结果与钻探结果相吻合。由此进一步说明了正交各向异性频散属性$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $对裂缝储层的含气饱和度较敏感，强各向异性频散的储层中含气饱和度较高，可以用于识别储层含气性。

各向异性梯度$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $表征裂缝密度，从常规AVAZ反演结果(图 17c)可见，A井目标层位附近的各向异性梯度$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $大于D井，表明A井目标层的各向异性更强，裂缝密度更大。但相比于图 17b，各向异性梯度$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $在非储层位置也有较明显的异常。这是由于$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $与方位各向异性有关，而方位各向异性不仅包含了有效裂缝的影响，还可能与地应力等因素引起各向异性有关。因此，$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $在识别有效裂缝过程中可能存在较多背景干扰，具有一定的局限性。

FDAVO反演是一种广泛用于流体识别的方法，可根据不同的Zoeppritz方程近似式估算多种频散属性。本文采用基于Shuey近似式^[32]推导的纵波速度频散属性$ {I}_{\mathrm{a}} $，它反映储层内与空腔流体饱和度有关的速度频散。对比图 17b与图 17d可知，A井处储层纵波速度频散属性$ {I}_{\mathrm{a}} $幅度明显弱于正交各向异性梯度频散属性$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $，这是由于纵波速度频散属性$ {I}_{\mathrm{a}} $是根据FDAVO反演方法计算所得，该方法所依据的近似式方程是假设地下储层为各向同性介质而推导所得。它没有考虑地震数据中的方位信息，忽略了裂缝引起的各向异性的影响。与此同时，该区域的各向异性较强，与FDAVO的假设条件不匹配，导致$ {I}_{\mathrm{a}} $的反演结果难以反映目标层的实际情况。另外，高含气区和低含气区在$ {I}_{\mathrm{a}} $响应方面的差异较小，不利于区分不同含气饱和度的储层。相比之下，OA-FDAVAZ方法的反演结果与气井的生产数据一致。

从反演结果平面图(图 18)看，$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $在A井处较强，在D井处较弱，这与实际生产结果一致。另外，背景岩层中的能量被压制，裂缝储层与背景岩层分离，这说明$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $准确反映了高含气饱和度裂缝储层的分布范围。A井和D井处$ {B}_{\mathrm{a}\mathrm{n}\mathrm{i}} $的幅值均较强(图 18b)，两井差异较小，干扰较多。A井和D井处$ {I}_{\mathrm{a}} $差别不大(图 18c)，区分高产井和低产井的效果较差。

综上所述，$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $主要由裂缝含气饱和度决定，与实际生产数据更吻合，且准确区分了高产井与低产井。同时，表明了基于正交各向异性介质的$ {I}_{\mathrm{a}\mathrm{n}\mathrm{i}} $可以有效地识别储层含气性，并且OA-FDAVAZ反演方法对于裂缝性储层中的含气性预测是可靠的。

5 结论

本文基于Chapman双裂缝模型，分析了不同饱和度情况下流体的频率、入射角和方位角等参数与弹性参数的关系，提出了正交各向异性介质频变AVAZ反演方法，并将其应用于四川盆地五百梯页岩储层的裂缝含气性预测，得出以下结论。

(1) 当其他参数保持一致时，含水储层弹性参数的频散程度远小于含气储层，说明利用各向异性频散属性可识别流体类型。此外，弹性参数的频散程度随地震波入射角的增加而减小，随方位角的增加呈现周期性变化，当地震波传播方向与垂直裂缝的法向平行时，频散程度出现最大值。

(2) 本文OA-FDAVAZ反演方法考虑了地震数据的方位和频率信息，通过高分辨率谱分解技术和最小二乘法，得到了正交各向异性频散属性，并提出频散组合因子的计算方法。

(3) 理论模型测试表明含水和含气裂缝层均表现出高值异常，但含气层的频散值要明显高于含水层，可见正交各向异性频散属性不仅对流体敏感，而且易于区分不同类型的流体。

(4) 相比于常规AVAZ反演和FDAVO反演，本文方法在实际资料应用中更能准确地刻画含气储层的分布。