0 引言

如今，随着世界各国对石油、天然气等能源的需求日益增长，以裂缝型储层为代表的非常规油气藏逐渐接替常规油气藏，成为油气勘探、开发的新目标。裂缝是裂缝型储层中流体的通道，连接了储层内的溶洞和溶孔^[1]，可为储层预测提供重要信息^[2]。因此裂缝对油气的分布具有重要影响。

目前，常规的裂缝预测方法主要包括地震属性分析方法和地震各向异性方法。基于叠后地震属性的分析技术发展相对成熟，如相干属性^[3]、曲率属性^[4]、蚂蚁追踪^[5]等，可表征长度大于1/4地震波长的大尺度裂缝。受地震数据纵向分辨率的影响，叠后属性技术还无法识别更小尺度的微裂缝。存在微裂缝的地下介质的各向异性特征显著^[6]，其中具有一组竖直裂缝且对称轴水平的横向各向同性（HTI）介质是一种典型的各向异性介质^[7]。在HTI介质中纵波振幅、纵波速度（走时）以及横波分裂特征随着入射角和裂缝方位角的改变而发生规律性变化。因此可以利用这三种信息反演裂缝方位角以及各向异性参数，进而估计裂缝的发育方向和相对发育密度^[8-9]。

然而，利用横波分裂预测裂缝需要采集高质量的横波数据，一般采用井下三分量数据^[10]或海底四分量数据^[11]。利用纵波速度（走时）特征只能预测大套储层中的裂缝，导致预测结果的纵向分辨率较低^[12]。利用纵波振幅各向异性特征预测裂缝的分辨率较高，稳定性强，是目前常用的方法^[13]，包括纵波振幅随方位角变化（AVAZ）^[14]、纵波振幅随炮检距变化（AVO）^[15]两类。大多数的各向异性反演方法需要建立初始模型，其中一种建模方法是将各向同性模型作为各向异性反演的初始模型^[16]，但建模结果与真实情况有较大偏差。此外，利用岩石物理信息^[17]或测井信息^[18]也可以建立各向异性参数模型，通过估算井旁道各向异性参数，并在层位约束下插值，从而得到各向异性参数初始模型。此种建模方法在井附近精度较高，但无法保证远井处的模型精度。由于三维初始模型的精度还取决于井的数量，受勘探过程中钻井数量的限制，特别是在勘探初始阶段，很难保证初始模型的准确性，从而影响最终各向异性参数反演结果^[19]。

虽然利用纵波走时信息预测裂缝的纵向分辨率低，但是在地震数据的每一个成像点位置（特别是在井间）均可提取各向异性信息，较插值方法可靠。针对地震各向异性反演的初始模型的建立问题，本文在HTI介质假设下，利用炮检距矢量片（OVT）数据的走时信息获取各向异性参数（$ \varepsilon $、$ \delta $），然后基于各向异性参数之间的经验关系估计另外一个各向异性参数$ \gamma $，从而建立初始各向异性参数模型。再利用各向同性反演获得的弹性参数建立初始弹性参数模型。最后，基于Rüger反射系数近似公式，利用高斯—牛顿法反演弹性参数和各向异性参数。合成地震数据和实际数据测试结果证明了方法的可行性和有效性。

1 方位走时各向异性建模 1.1 纵波走时拾取

利用叠前全方位角道集的纵波走时信息反演方位走时各向异性，首先要得到目的层纵波走时。在选取的时窗内，采用走时扫描函数获取纵波走时

$ Q\left(T\right)=W\left({\varphi }_{n}, T\right) \left|{}_{{T}_{0}}^{{T}_{1}}\right. $

(1)

式中：$ W\left({\varphi }_{n}, T\right)\mathrm{为} $时窗内$ T $时刻的振幅值，$ {\varphi }_{n}\left(n=\mathrm{1, 2}, \cdots , N\right) $为测线方位角，$ N $为方位测线数目；$ {T}_{0} $、$ {T}_{1} $分别为时窗的下限和上限。

在选定的时窗内，当振幅取最大$ \mathrm{值}W\left({\varphi }^{\mathrm{o}\mathrm{b}\mathrm{s}}, {T}^{\mathrm{o}\mathrm{b}\mathrm{s}}\right) $时，扫描得到的走时$ {T}^{\mathrm{o}\mathrm{b}\mathrm{s}} $即为测线方位角为$ {\varphi }^{\mathrm{o}\mathrm{b}\mathrm{s}} $的目的层走时。

1.2 HTI介质的走时

HTI介质（图 1）的纵波走时为^[20]

$ T=2\sum\limits_{k=1}^{K}\frac{{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(n\right)}+{a}_{1}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2} $

(2)

其中

$ \left\{\begin{array}{l}{a}_{0}^{\left(n\right)}={v}_{{}_{\mathrm{{\rm P}}0k}}^{-2}(1-2{\varepsilon }_{k})\\ {a}_{1}^{\left(n\right)}=2{v}_{{}_{\mathrm{{\rm P}}0k}}^{-2}(2{\varepsilon }_{k}-{\delta }_{k})\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}({\varphi }_{n}-{\phi }_{k})\\ {a}_{2}^{\left(n\right)}=2{v}_{{}_{\mathrm{{\rm P}}0k}}^{-2}({\delta }_{k}-{\varepsilon }_{k})\mathrm{c}\mathrm{o}{\mathrm{s}}^{4}(\varphi n-{\phi }_{k})\end{array}\right. $

(3)

式中：$ K $为地层层数；$ {v}_{\mathrm{{\rm P}}0k} $为垂直入射时第$ k $层的纵波速度；$ {\varepsilon }_{k} $、$ {\delta }_{k} $为第$ k $层的各向异性参数；$ {\theta }_{k} $为第$ k $层的纵波入射角；$ {d}_{k} $为第$ k $层的地层厚度；$ {\phi }_{k} $为第$ k $层的裂缝对称轴方位角。

式（2）是经过多次近似得到的，因此当各向异性较弱或者$ {\theta }_{k} $较小时，该式才有较好的适用性，因此该式适用于起伏较小的地层^[21]。

根据时距关系利用Dix公式^[22]计算$ {v}_{\mathrm{{\rm P}}0k} $，即

$ {v}_{\mathrm{{\rm P}}{0}_{k}}^{2}=\frac{{T}_{0, k}{v}_{\mathrm{R}, k}^{2}-{T}_{0, k-1}{v}_{\mathrm{R}, k-1}^{2}}{{T}_{0, k}-{T}_{0, k-1}} $

(4)

式中：$ {T}_{0, k} $为第1层到第$ k $层的自激自收时间；$ {v}_{\mathrm{R}, k} $为第1层到第$ k $层的均方根速度。

$ \mathrm{由}\mathrm{于}{\theta }_{k} $决定了地震波传播路径，并且满足费马原理，因此约束条件为

$ X=2\sum\limits_{j=1}^{k}\left({d}_{j}\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{j}\right) $

(5)

$ \frac{\partial }{\partial \left({d}_{j}\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{j}\right)}\left(T-\lambda X\right)=0 $

(6)

式中：$ X $为炮检距；$ \lambda $为拉格朗日算子。

式（2）表明目的层走时受上覆地层各向异性以及裂缝方位角的影响。因此为了计算目的层的各向异性参数，首先应该得到上覆地层的各向异性参数以及裂缝方位角。本文采用逐层反演方法，利用目的层走时和上覆地层的各向异性参数计算目的层的各向异性参数。首先，利用单层走时（式（2）中K=1）获取第一层的各向异性参数。然后再逐层获取下伏地层的各向异性参数。

根据式（2），得到第$ k $层的层间走时为

$ t=\frac{2{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(n\right)}+{a}_{1}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2} $

(7)

设第$ k $层上覆地层的各向异性参数以及裂缝方位角已知，根据式（7），不同测线方位的第$ k-1 $层的地震波走时为

$ {\mathit{\boldsymbol{{T}}}}_{N, k-1}^{\mathrm{c}\mathrm{a}\mathrm{l}}=\left(\begin{array}{c}{t}_{\mathrm{1, 1}}+{t}_{\mathrm{1, 2}}+\cdots +{t}_{1, i}+\dots +{t}_{1, k-1}\\ {t}_{\mathrm{2, 1}}+{t}_{\mathrm{2, 2}}+\cdots +{t}_{2, i}+\dots +{t}_{2, k-1}\\ \vdots\\ {t}_{n, 1}+{t}_{n, 2}+\cdots +{t}_{n, i}+\dots +{t}_{n, k-1}\\ \vdots\\ {t}_{N, 1}+{t}_{N, 2}+\cdots +{t}_{N, i}+\dots +{t}_{N, k-1}\end{array}\right) $

(8)

式中$ {t}_{n, i} $（$ i= $1，2，$ \cdots $，$ k-1 $）为不同测线方位的某一层的地震波层间走时，当$ k=1 $时，$ {t}_{n, i}=0 $。由于第$ k $层的各向异性参数和裂缝方位角未知，则第$ k $层在不同测线方位的地震波层间走时为

$ {\mathit{\boldsymbol{{t}}}}_{N, k}^{\mathrm{c}\mathrm{a}\mathrm{l}}=\left(\begin{array}{c}\frac{2{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(1\right)}+{a}_{1}^{\left(1\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(1\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2}\\ \frac{2{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(2\right)}+{a}_{1}^{\left(2\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(2\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2}\\ \vdots\\ \frac{2{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(n\right)}+{a}_{1}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(n\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2}\\ \vdots\\ \frac{2{d}_{k}}{\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{k}}{\left[{a}_{0}^{\left(N\right)}+{a}_{1}^{\left(N\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}{\theta }_{k}+{a}_{2}^{\left(N\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}{\theta }_{k}\right]}^{1/2}\end{array}\right) $

(9)

在式（5）和式（6）的约束下，地震波在不同测线方位的炮点到检波点的总走时为

$ {\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{c}\mathrm{a}\mathrm{l}}=\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{ }({\mathit{\boldsymbol{{T}}}}_{N, k-1}^{\mathrm{c}\mathrm{a}\mathrm{l}}+{\mathit{\boldsymbol{{t}}}}_{N, k}^{\mathrm{c}\mathrm{a}\mathrm{l}}) $

(10)

根据式（1），从叠前全方位角道集中拾取第$ k $层在不同测线方位的实际走时信息为

$ {\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}={\left(\begin{array}{ccc}{T}_{1, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}& {T}_{2, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}& \begin{array}{cc}\cdots & {T}_{N, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}\end{array}\end{array}\right)}^{\mathrm{T}} $

(11)

根据式（10）和式（11），利用最小二乘法可反演$ {\varepsilon }_{k}\mathrm{、}{\delta }_{k} $以及$ {\phi }_{k}, $目标函数为

$ J\left(\mathit{\boldsymbol{{m}}}\right)=\left|\right|{\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{c}\mathrm{a}\mathrm{l}}|{|}_{2}^{2} $

(12)

式中$ \mathit{\boldsymbol{{m}}} $为反演参数矩阵。

通过纵波方位走时各向异性反演方法无法直接获取各向异性参数$ \gamma $。但是在一般情况下$ \varepsilon $和$ \gamma $的变化趋势一致，因此假设二者之间具有线性关系^[23]

$ \gamma =A\times \varepsilon +B $

(13)

式中$ A $、$ B $为拟合系数。

通过最小二乘方法迭代计算式（12），当合成走时$ {\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{c}\mathrm{a}\mathrm{l}} $与实际走时$ {\mathit{\boldsymbol{{T}}}}_{N, k}^{\mathrm{o}\mathrm{b}\mathrm{s}}\mathrm{的} $差最小时，输出$ \varepsilon \mathrm{、}\delta $以及φ。最后根据式（13）计算$ \gamma $。

在每一个成像点均可独立利用走时信息反演各向异性参数，不受层位和测井信息约束。然而，由于走时反演只能识别大套储层，并且由时距关系计算的地层速度为近似解，因此得到的参数精度以及分辨率不满足储层预测的要求^[24]。为此，本文将走时反演获取的各向异性参数结果作为AVAZ反演的初始模型。

2 各向异性反演

根据弱各向异性理论假设，Rüger等^[25-26]推导了HTI介质纵波反射系数$ R $随纵波入射角$ \theta $和裂缝对称轴方位角$ \phi $变化的公式

$ \begin{array}{c} R\left(\theta , \varphi \right)=\frac{1}{2}\frac{\Delta Z}{\stackrel{-}{Z}} +\frac{1}{2}\left\{\frac{\Delta {v}_{\mathrm{P}0}}{{\stackrel{-}{v}}_{\mathrm{P}0}}-{\left(\frac{2{\stackrel{-}{v}}_{\mathrm{S}0}}{{\stackrel{-}{v}}_{\mathrm{P}0}}\right)}^{2}\frac{\Delta G}{\stackrel{-}{G}}+\left[\Delta \delta +2{\left(\frac{2{\stackrel{-}{v}}_{\mathrm{S}0}}{{\stackrel{-}{v}}_{\mathrm{P}0}}\right)}^{2}\Delta \gamma \right]\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}(\varphi -\varphi )\right\}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta + \\ \frac{1}{2}\left[\frac{\Delta {v}_{\mathrm{P}0}}{{\stackrel{-}{v}}_{\mathrm{P}0}}+\Delta \varepsilon \mathrm{c}\mathrm{o}{\mathrm{s}}^{4}(\varphi -\varphi )+\Delta \delta \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}(\varphi -\varphi )\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}(\varphi -\varphi )\right]\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \end{array} $

(14)

式中：$ {v}_{\mathrm{P}0} $为纵波速度；$ {v}_{\mathrm{S}0} $为横波速度；$ \rho $为介质密度；$ Z=\rho {v}_{\mathrm{P}0} $为垂向纵波阻抗；$ G=\rho {v}_{\mathrm{S}0}^{2} $为垂向横波剪切模量；上划线“$ - $”代表反射界面上、下介质的物理量平均值；“$ \Delta $”代表反射界面上、下介质的物理量之差；$ \varphi $为测线方位角。

在各向异性反演中，采用高斯—牛顿法反演各向异性参数。为了提高反演精度以及解决高斯—牛顿法中的局部最小解问题，将各向异性反演分为两个阶段。同时考虑到各向异性随入射角的增加而越明显这一特性，在第一阶段使用小角度（θ < 20°）道集进行各向同性AVO反演，获得弹性参数（$ {v}_{\mathrm{P}0}\mathrm{、}{v}_{\mathrm{S}0}\mathrm{、}\rho $）初始模型。在第二阶段，在弹性参数初始模型和各向异性参数初始模型约束下，使用中、大角度（θ > 20°）道集进行各向异性AVO反演。上述过程往往是繁琐且耗时的。因此，为了优化迭代效率，Lu等^[27]根据阻尼最小二乘法对高斯—牛顿反演算法正则化，并利用欧氏距离相似性建立迭代反演停止条件，极大地提高了计算效率，且得到了更精确的反演结果。该优化方法的目标函数为

$ \Delta \mathit{\boldsymbol{{M}}}\cong {\left(\mathit{\boldsymbol{{H}}}+\mu \mathit{\boldsymbol{{I}}}\right)}^{-1}{\mathit{\boldsymbol{{J}}}}_{\mathrm{H}\mathrm{T}\mathrm{I}}^{\mathrm{T}}\left({\mathit{\boldsymbol{{R}}}}_{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathit{\boldsymbol{{R}}}}_{0}\right) $

(15)

其中

$ \mathit{\boldsymbol{J}}_{{\rm{HTI}}}^{\rm{T}} = \left( {\begin{array}{*{20}{c}} {\frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial {v_{{\rm{P}}0}}}}\quad \frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial {v_{{\rm{S}}0}}}}\quad \frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial \rho }}\quad \frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial \varepsilon }}\quad \frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial \delta }}\quad \frac{{\partial R\left( {{\theta _1}, \phi } \right)}}{{\partial \gamma }}}\\ {\frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial {v_{{\rm{P}}0}}}}\quad \frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial {v_{{\rm{S}}0}}}}\quad \frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial \rho }}\quad \frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial \varepsilon }}\quad \frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial \delta }}\quad \frac{{\partial R\left( {{\theta _2}, \phi } \right)}}{{\partial \gamma }}}\\ \vdots \\ {\frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial {v_{{\rm{P}}0}}}}\quad \frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial {v_{{\rm{S}}0}}}}\quad \frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial \rho }}\quad \frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial \varepsilon }}\quad \frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial \delta }}\quad \frac{{\partial R\left( {{\theta _m}, \phi } \right)}}{{\partial \gamma }}} \end{array}} \right) $

(16)

$ \mathit{\boldsymbol{{H}}}=\left(\begin{array}{c}\frac{\partial R(\theta , \phi )}{\partial \mathit{\boldsymbol{{E}}}}\\ \frac{\partial R(\theta , \phi )}{\partial \mathit{\boldsymbol{{A}}}}\end{array}\right)\cdot \left(\begin{array}{cc}\frac{\partial R(\theta , \phi )}{\partial \mathit{\boldsymbol{{E}}}}& \frac{\partial R(\theta , \phi )}{\partial \mathit{\boldsymbol{{A}}}}\end{array}\right) $

(17)

式中：$ \mathit{\boldsymbol{{M}}} $为参数模型矩阵；μ为一个标量；$ \mathit{\boldsymbol{{I}}} $为单位矩阵；$ \mathit{\boldsymbol{{H}}} $为Hessian矩阵；$ \mathit{\boldsymbol{{J}}} $为Jacobin矩阵；$ {\mathit{\boldsymbol{{R}}}}_{\mathrm{o}\mathrm{b}\mathrm{s}} $和$ {\mathit{\boldsymbol{{R}}}}_{0} $分别为观测系数矩阵和初始反射系数矩阵；$ \mathit{\boldsymbol{{E}}} $为弹性参数矩阵；$ \mathit{\boldsymbol{{A}}} $为各向异性参数矩阵；m为入射角道集的个数。

在迭代过程中，将每次迭代增量$ \Delta \mathit{\boldsymbol{{M}}} $添加到初始模型$ {\mathit{\boldsymbol{{M}}}}_{0} $中。在给定的迭代约束下，即当$ {\mathit{\boldsymbol{{R}}}}_{\mathrm{o}\mathrm{b}\mathrm{s}} $接近$ {\mathit{\boldsymbol{{R}}}}_{0} $时，算法迭代停止。最后便可得到精确的$ {v}_{\mathrm{P}0}\mathrm{、}{v}_{\mathrm{S}0}\mathrm{、}\rho \mathrm{、}\varepsilon \mathrm{、}\delta $和$ \gamma $结果。图 2为基于HTI介质的方位走时各向异性反演流程。

3 合成数据测试

建立一个水平层状的HTI模型（图 3），其弹性参数和各向异性参数如表 1所示。Tsuneyama等^[28]通过分析盐水饱和砂岩和页岩的测井数据得到式（13）中$ \mathrm{的}A=1.2006 $，$ B=-0.0282\mathrm{。} $本文据此设置$ \gamma $，并且在实际数据测试中也采用上述值。Thomsen^[29]对各向异性介质的研究表明，大多数沉积岩呈弱各向异性，其各向异性参数取值范围一般为0~0.2。

对于叠前全方位角道集，设置X=280 m以及φ间隔为15°；对于入射角道集，设置φ为0°、60°和120°以及θ为5°~40°。利用式（14）计算图 3在时域的PP波反射系数。然后，将反射系数序列与主频为30 Hz的雷克子波卷积，生成叠前全方位角道集（图 4）和入射角道集（图 5）。

对每一层地层底界面对应的同相轴提取走时（图 4中红线）。为了验证走时反演效果，对第1层~第4层的各向异性参数的所有可能值进行最小二乘扫描，得到反演结果（表 2~表 5）。反演结果表明，走时反演的最大相对误差在10%以内，该精度对于建立初始模型是可以接受的。方位走时各向异性反演的模型测试结果表明，在没有任何约束条件下，方位走时各向异性反演依然能得到较准确的结果，说明反演方法的稳定性较高。

利用式（13），根据反演得到的$ \varepsilon $计算$ \gamma $，并利用最终反演结果（$ \varepsilon \mathrm{、}\delta \mathrm{、}\gamma $）建立各向异性参数初始模型（图 6）。

各向异性反演分两个阶段。第一阶段为各向同性反演，第二阶段为各向异性反演。对小角度（5°~20°）的入射角道集（图 5）进行各向同性反演，共迭代20次。当迭代停止时，各向同性反演残差（$ {\mathit{\boldsymbol{{R}}}}_{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathit{\boldsymbol{{R}}}}_{0} $）很小，约为0.16。反演结果（图 7）表明，反演结果与真实模型匹配较好。因此，可利用该反演结果建立弹性参数初始模型。

在图 6和图 7的约束下，对中、大角度（20°~40°）的入射角道集（图 5）进行各向异性AVAZ反演，共迭代了28次。当迭代停止时，各向异性AVAZ反演残差很小，约为0.32。各向异性反演结果（图 8）表明，在初始模型值误差较小的情况下，反演值始终与理论值吻合较好。因此，若能提供较精确的各向异性参数初始模型，基于HTI介质的各向异性AVAZ反演可以得到更准确的Thomsen参数。

4 实际数据测试

本文利用四川盆地某页岩油储层的OVT地震数据测试所提方法的效果。页岩储层在构造应力作用下，易产生裂缝。岩心以及地层微电阻率扫描成像（FMI）测井中均发现高角度裂缝，整体具有明显的HTI介质特征。

以某一个OVT为例，炮检距为1100~1150 m的方位覆盖次数较高，因此抽取炮检距为1127 m的方位角道集（图 9）进行方位走时各向异性反演。

与常规的叠前地震数据相比，OVT数据不仅包含炮检距信息，还包含方位角信息。因此，为了减小走时拾取误差，首先根据地质资料的岩性信息以及地震记录中同相轴的连续性，将目的层细分为6层（图 9a）。其次，拾取OVT数据中同一炮检距的多个方位的走时，并且剔除异常值。即使其中有几道数据的误差较大，但对于拾取的大量走时数据而言，其影响可以忽略不计。同时，为了增加走时反演的稳定性，对拾取的走时进行方位插值，从而得到任意方位走时（图 9b）。通过图 9a以及时距关系计算分层地层参数，得到地层参数分层结果（表 6）。

图 10为裂缝方位预测结果。可以看出，基于走时反演得到的裂缝方位（图 10b）与实际裂缝发育方位（图 10a）基本一致，为AVAZ反演提供了良好的模型基础。利用计算得到的地层参数以及走时拾取结果进行方位走时各向异性反演得到各向异性参数$ \varepsilon $、$ \delta $，并利用式（13）计算各向异性参数$ \gamma $，得到方位走时各向异性反演建模结果（图 11）。然后，抽取井点附近θ为8°~35°且φ为0°、45°、90°、135°的入射角道集（图 12）。对小角度（8°~20°）的地震数据进行各向同性反演，对中、大角度（20°~35°）的地震数据进行各向异性反演。各向同性反演结果（图 13）表明，弹性参数（$ {v}_{\mathrm{P}0} $、$ {v}_{\mathrm{S}0} $、$ \rho $）反演结果与实际测井曲线变化趋势基本吻合。因此可以将各向同性反演结果作为各向异性反演的弹性参数初始模型。

在图 11、图 13的约束下进行各向异性反演，得到各向异性反演结果（图 14）。由图可见：①弹性参数反演结果（$ {v}_{\mathrm{P}0} $、$ {v}_{\mathrm{S}0} $、$ \rho $）与井点处的实际结果较吻合。②1220~1420 ms井段的FMI数据指示目的层上部的裂缝发育程度大于下部；各向异性参数曲线指示裂缝发育位置的各向异性程度相对较大，且整体上目的层段上部（层位Hor_1与Hor_2之间的地层）的各向异性程度大于下部（Hor_2与Hor_3之间的地层），与FMI数据基本吻合。

根据图 2，抽取CDP为450~600的过井叠前地震数据（叠后剖面如图 15所示）进行反演。其中，抽取炮检距为1127 m的叠前方位道集进行方位走时各向异性反演，抽取φ为0°、45°、90°、135°的入射角道集进行各向异性反演。根据方位走时各向异性反演以及各向同性反演得到初始模型剖面（图 16），根据各向异性反演得到反演剖面（图 17）。由图 17可见，各向异性参数高值区（红色区域）裂缝较发育，整体上目的层上部裂缝较下部发育，与井上裂缝发育情况吻合，证实了所提方法的可行性。

5 结束语

传统反演方法利用测井资料和层位插值构建初始模型的做法导致初始模型精度较低。为此，本文充分考虑了OVT数据中的走时和振幅信息，提出了一种基于HTI介质的方位走时各向异性建模的反演方法。首先，获取各个地层的方位纵波走时。然后，利用得到的方位走时信息进行方位走时各向异性反演，并利用反演结果建立各向异性参数初始模型。最后，基于Rüger反射系数近似公式，采用高斯—牛顿法进行各向异性AVAZ反演获得各向异性参数。模型测试表明，方位走时各向异性反演可以在无测井约束的情况下获得精度较高的各向异性参数。在获取各层的各向异性参数的基础上，通过各向异性AVAZ反演得到了较为精确且分辨率更高的反演结果。四川盆地实际数据测试结果表明，利用所提方法成功地获得了储层的弹性和各向异性信息。