通过模拟计算天然气藏地震波场，获取地震波波动方程，可以明确不同各向异性介质中的地震波传播过程与规律，为地震响应特征分析和裂缝参数反演提供理论支撑。

采集地震数据、测井数据和地质勘探数据，通过时深转换将各个数据源归一化到统一的时间域和深度域，采用加权平均融合方法对多源数据进行整合与预处理^[5]。公式如下：

式（1）中，S、L、G分别表示地震数据、测井数据和地质构造数据；w_s、w_l、w_g分别表示对应数据源的权重系数；D_f表示加权融合后的数据集合。

在横波和纵波的地震勘探中，在应力的影响下，介质中的地震波传播路径是不规则的，无法准确获取地震波的加速度^[6]。因此，为明确地震波的传播规律，采用固体的非线性弹性理论，构造地震波在天然气藏介质中的声弹性刚度。即：

式（2）中， b表示受应力地层的有效刚度张量； u_k表示地层k受到的地震波扰动大小；α表示地震波相对反射系数随入射角的变化梯度；X_i表示反射界面i的单元面积；β表示纵横波速比；X_j表示入射界面j的单元面积；c₀表示纵横波的波阻抗。

利用微分方程构造二维各向异性介质的二阶声波波动方程式^[7]。并将储层表面处的检波地震数据加载为第二轮模拟的震源信号，形成包含表面多次波的传播方程^[8]。即：

P ( z ₁, z ₀)=W⁺( a, v²) Q⁺(f), e + E⁺(t), q （3）

式（3）中，W⁺表示上行波地震信号；v表示各向异性弹性张量；Q⁺表示上行波的多次波阶数；f表示声压；e表示地震波的透射损失；t表示地震波传播时间；q表示二阶中心差分的空间网格步长；P (z₁, z₀)表示天然气储层(z₁, z₀)位置处的地震波信号。

利用公式（5）所计算的天然气储层表面多次波模拟方程，在每一轮迭代中不断向表层加入上一层计算得到的地震波模拟数据^[9]，即可得到地震波的波场快照，以此为依据，对天然气藏的地震响应各向异性特征进行分析。

地震波在天然气藏介质中传播时，沿不同方位角传播的速度会有所不同，同时，地震波在反射或透射过程中，其振幅会随着方位角的变化而变化^[10]。在裂缝发育的天然气藏中，由于裂缝的存在改变了介质的弹性性质，导致地震波在不同方位上的反射系数或透射系数不同，从而引起振幅的方位变化，因此，本研究主要针对天然气藏地震响应的速度和振幅变化进行分析^[11]。

在弱散射假设条件下，基于水平地应力作用下天然气藏储层的介质线性化刚度张量矩阵，推导弹性界面满足应力连续和位移连续边界条件时地震波的入射系数与透射系数^[12]。公式为：式（4）中，A₂表示应力上限；A₁表示位移上限；v_p1、v_p2分别表示P波下行和S波上行的传播速度；r_p表示入射系数；ε表示横波粒子运动协方差矩阵元素；K表示速度各向异性幅度；P(z₁, z₀)表示天然气储层(z₁, z₀)位置处的地震波信号；h表示裂缝引起的扰动项；T表示传播方向与对称轴夹角；r_s表示透射系数。

根据测井和岩心数据构建裂缝型储层介质模型，具体包括水平裂缝（VTI）、垂直裂缝（HTI）和垂直于水平裂缝组合而成的正交各向异性介质裂缝（OA）模型^[13]。如图1所示。

图1 天然气储层裂缝型介质模型Fig.1 Fractured medium model of natural gas reservoir

图1中，x₁ x₂ x₃为介质模型坐标系，F₁为沿模型对称轴施加的水平单轴应力，以模拟介质的各向异性大小，该作用力的大小通常与裂缝的强度有关。表达式为：

式（5）中，λ表示介质的密度扰动；ν表示裂缝强度；θ表示方位角。

基于不同裂缝类型的介质模型，利用中心点道集提取地震波响应的速度谱，通过拟合各向异性速度函数确定地震响应的速度和振幅。公式为：

式（6）中，B表示速度函数；v'表示地震响应速度特征；δ₁、ε₁、γ₁均表示Thomsen参数，为各向异性系数；Γ表示各向异性参数差；y表示速度谱；G表示地震波振幅。

通过式（6）计算地震波作用下天然气储层的声波反射振幅，若某区域振幅出现极大值或极小值，则表明出现了波峰响应或波谷响应^[14]。基于以上计算过程，即可确定天然气藏地震响应在速度和振幅方面的各向异性特征，进而通过特征反演构建裂缝预测模型。

由于地震反射速度和振幅的变化特征中包含裂缝信息，所以利用天然气藏地震响应各向异性特征对裂缝参数进行反演以构建预测模型。在实际应用中，地震波的密度可以由叠前地震数据弹性反推得到，但叠前道集可能会存在随机噪声干扰、拉伸畸变和偏离等问题，因此，采用“去噪—拉平—去除”的方法对叠前道集进行优化处理。该过程可由下式表示：

式（7）中， n表示给定的时窗长度； δ_τ表示时窗τ对应的偏移距；-δ表示偏移距平均值；κ_τ表示时窗τ对应的噪声分量；-κ表示噪声分量平均值；z_ψ表示道ψ与初始道的相关系数。将与初始道相关系数最高的一道作为模型道，替代噪声分量，得到拉平道集。

将分析的地震响应各向异性特征作为反演参数，结合优化后的叠前道集对裂缝参数进行反演^[15]。公式为：

式（8）中， v'、 G分别表示地震波速度和振幅；ϕ(ϖ) 表示具有ϖ给入射角的地震叠前道集；Δρ表示应力增量；s₀表示地震子波向量；V表示反演得到的裂缝参数。

利用下式将反演得到的裂缝参数表示为裂缝密度，由此建立裂缝预测模型。即：

式（9）中，m₁、m₂分别表示背景介质横波和纵波的速度；Ε表示裂隙的长半轴；ω₀表示切向裂隙弱度；Ω表示裂缝密度预测值。

利用式（9）所示模型即可实现对天然气藏储层裂缝的预测。

基于上述研究数据，通过构建地震波的二维波动方程并利用天然气储层表面多次波模拟方法，模拟地震合成数据，从而获取地震波波场快照。如图4所示。

图4 三轮迭代轮次表面多次地震波波场快照Fig.4 Seismic wavefield snapshots of surface multiples after three iteration

图4中，（a）、（b）为第一轮到第二轮的上行波波场，（c）、（d）为第一轮到第三轮的下行波波场，每增加一轮正演模拟轮次，多次波阶数就增加一阶，形成多次波场。基于地震波波场与区域地层结构特征结合已有的测井资料与地震剖面解释结果，对研究区域天然气藏有利层系的两口测井进行地震响应各向异性特征，结果如图5和图6所示。

图5 A井地震响应各向异性特征Fig.5 Anisotropic characteristics of seismic response of well A

图6 B井地震响应各向异性特征Fig.6 Anisotropic characteristics of seismic response of well B

如图5（a）所示，目的层的地震波反射速度为4200 m/s，说明天然气藏储层呈现出近似垂直的条带状分布特征，使得天然气藏被相对致密的岩石部分封隔，从而导致天然气流体的弹性模量增加，进而提升了地震波的传播速度。如图5（b）所示，在深度-1000~-1200 m范围内，存在明显的弱波谷响应，波谷反射不连续，弱波谷的边界比较清晰，形成“波谷阴影”效应，主要是由于天然气具有较低的密度和较高的压缩性，会显著改变地层的弹性性质，使得目的层的声阻抗明显低于周围围岩，这种差异会导致地震波在含气层界面处发生反射和透射时，反射系数减小，从而使得记录到的振幅值降低，形成弱波谷响应。

根据图6（a）可以看出，在-1400~-1200 m的深度带，出现了明显的速度异常带，速度为5600 m/s，说明B井天然气大部分以游离态聚集在裂缝中，降低了地层的有效应力，使岩石颗粒之间的结合力减弱，从而导致地震波传播速度升高。根据图6 （b）可知，B井天然气藏地震剖面的异常响应为强波峰反射，说明其振幅强度较大，但影响范围较为局限，未向更深或更浅层大规模延伸，分析原因可知，B井天然气层还具有较高的饱和度且赋存于不同岩性地层界面处，导致地震波会发生强烈的反射，增强了反射能量，形成了强波峰。

为验证本文对于地震响应各向异性特征分析的正确性，引入文献[1]SV/P波斜入射、文献[2]基于真实地形、文献[3]基于地震数值模拟和文献[4]地质统计学四种方法与本文方法进行对比，以全波形反演（FWI）得到的各向异性参数为理论值，将本文方法和文献方法的计算结果与理论结果相比较，评估不同方法的分析准确性。结果如表3所示。

表3 地震响应各向异性参数计算偏差Table 3 Calculation deviation of anisotropic parameters of seismic response

从表3中的数据可以明显看出，本文方法在计算各向异性参数时，偏差均控制在较低水平，三个参数的计算偏差值均低于2.0，表明本文方法通过多源数据整合和地震波场反演能够更准确地捕捉介质的各向异性特征，计算出与理论值更为接近的各向异性参数，从而提高地震响应各向异性特征分析的准确性。而文献[1]方法过于简化各向异性介质模型，忽略了介质的非均质性，导致计算结果偏差较大；文献[2]方法只使用了部分地震波的传播信息，而忽略了其他模式波和频段信息，导致计算结果不准确；文献[3]在处理地震数据时，对噪声的抑制能力较弱，降低了数据质量，会使数据质量下降，导致计算出的各向异性参数偏离理论值；文献[4]由于无法准确描述地震波与各向异性介质的相互作用，从而影响了参数的计算精度。综合对比分析可知，本文方法能够更真实地模拟地震波在各向异性介质中的传播行为，提高各向异性参数的计算精度，进而确保地震响应各向异性特征分析的准确性。

通过对地震响应各向异性特征进行反演，并结合地震叠前道集预测裂缝密度。首先应对地震叠前道集进行优化处理，以去除噪声与拉伸畸变，结果如图7所示。

图7 叠前道集优化前后对比Fig.7 Comparison of pre-stack gathers before and after optimization

根据图7可以看出，优化处理前，地震叠前道集存在随机噪声和畸变，影响数据质量，经过优化后，有效抑制了噪声分量，同时将目的层相轴拉平归位，提升了数据的保真度。

以图7叠前道集为基础，对研究区域天然气储层裂缝进行预测，同时将本文方法与文献[1]、文献[2]、文献[3]和文献[4]进行对比，得到的结果如表4所示。

表4 裂缝预测结果对比Table 4 Comparison of fracture prediction results

通过分析表4中的数据可知，本文方法针对VTI、HTI（水平横向各向同性）、OA三种裂缝类型，通过构建不同的裂缝介质模型，并结合地震响应各向异性特征反演，在500~2000 m的范围深度，本文方法预测出的裂缝密度与实测密度基本相等，通过计算得到平均预测误差仅为2.1%；而文献[1]仅使用纵波速度反演裂缝密度，忽略了横波分裂特征，预测误差较大，为17.3%；文献[2]采用统一的裂缝密度约束，未考虑深度相关变化，对不同类型的裂缝均未取得较优的预测结果，平均误差为14.6%；文献[3]通过局部迭代优化反演参数，易陷入次优解，深度地层预测偏差较大，平均误差为8.9%，文献[4]未对地震数据进行噪声抑制，导致预测结果不准确，平均误差为15.2%，进一步证明了本文方法在预测储层裂缝中的有效性，可为复杂裂缝型储层开发提供更可靠的裂缝参数。