0 前言
以Zoeppritz方程[1]为基础，振幅随偏移距变化（AVO）技术可以研究岩石的岩性、物性和含油气性质[2]。为了便于AVO技术的研究和应用，学者们对Zoeppritz方程开展了众多线性与非线性近似的研究。在众多线性近似中，Aki和Richards线性近似方程[3]最为著名，也是AVO技术蓬勃发展的前身与最广泛应用的线性近似之一。在Aki-Richards近似的基础上，通过研究岩石物性参数之间的关系，学者们分别构建了波阻抗近似[4]、泊松比近似[5]和含杨氏模量的线性近似[6]等反射系数近似方程。在油气探测中，流体因子通常被用作储层含油气性的评价指标。此概念最早由Smith和Gidlow[7]引入。随后，关于流体因子的研究便得到了快速的发展。基于Biot-Gassmann方程，Russell等[8]提出了应用最广泛的反映孔隙流体类型的Russell流体因子，并将其应用于叠前反演识别含气储层[9]。然而，线性反演方法受制于固有假设条件，难以适应当今复杂地质条件下勘探精度的需要。随着优化算法的不断更新以及计算机计算效率的提升，高阶非线性反射系数近似方程被逐渐应用于叠前地震反演中[10-12]，并在储层含气性预测中具有较高的精度，且更适合于反射界面变化剧烈、地震波大角度入射的情况。Shi等[13]提出了基于贝叶斯理论的非线性反演方法，用于直接估算流体体积模量和孔隙度。此外，基于人工神经网络逆算子估计算法（ANNI）的叠前非线性反演算法也被众多学者引入到叠前地震反演方法中，取得了较好的储层弹性及物性参数预测结果，并进一步证实了非线性反演方法的稳定性和有效性[14-15]。因此，发展一种精度更高的非线性反演方法开展储层流体因子的预测，能够为碎屑岩含气储层预测精度的提升给予技术支持。
当地下多孔岩石中存在流体时，岩石将表现出黏弹性质，这将导致地震波在传播过程中产生衰减作用[16-18]。因此，传统的完全弹性介质理论难以准确描述储层的弹性性质。为表征黏弹性介质地震波传播机理，最常用的方法是通过引入近似常Q模型将弹性介质转变为黏弹性介质[19-20]。岩石孔隙中存在油气等黏性流体，不仅会引起地震波的衰减，同时利用地震波在黏弹性介质中的频散效应，基于线性近似推导出的黏弹性介质反射系数的近似表达式，利用振幅随频率变化（AVF）的反演方法可以同时兼顾振幅与频率信息，提高反演结果的分辨率，并得到了广泛发展[21-26]。然而，上述频率相关（频变）反演方法大多数基于线性反射系数近似方程，针对非线性反射系数方程的频率相关反演开展较少，使得复杂碎屑岩储层的流体检测存在一定的误差。为此，Chen等[27]基于近似常Q模型理论，推导了黏弹介质下的精确反射系数方程，发展了一种非线性反演方法对地层品质因子进行估计，将频率相关反演拓展到了基于精确方程的非线性领域，并对地下含气储层进行了有效预测。然而，在非线性精确方程反演中，考虑频率信息所构建的角度-频率联合非线性反演方法，虽然能够进一步提高储层参数的反演精度，但是精确非线性反射系数方程的固有复杂特征使得反演优化的计算效率有所下降，制约了其在大面积工区中的应用能力。为此，如何基于简化后的精确方程非线性近似式，发展一种频率相关的非线性反演算法，是目前碎屑岩储层含气性预测面临的重要难题。
本文针对碎屑岩储层流体识别精度的需求，基于Bortfeld非线性近似[28]，推导了黏弹性介质条件下含纵横波品质因子的频率相关非线性反射系数方程。该方程在线性近似的基础上提高了反射系数的精度，且相比于精确方程，降低了方程的复杂程度，为反演的大面积实施奠定了基础。之后，设计了AVO反演与AVF反演相结合的两步非线性频变反演方法，在提高地层品质因子预测精度的同时，保证反演结果的鲁棒性。合成记录测试验证了所提出方法的合理性与稳定性。实际工区应用进一步验证了该方法用于碎屑岩储层流体识别的可靠性与适用性。
1 理论与方法
1.1 黏弹性介质非线性反射系数方程推导
式中： 、 、 与 、 、 分别为反射界面上层与下层的纵波速度、横波速度、密度； 与 分别为纵波入射角和纵波透射角并满足：
由于复杂碎屑岩储层在孔隙流体赋存时具有的一定的非均质性与黏弹性，使得实际地震波在传播过程中其速度具有随频率的变化特征[29]，因此式(1)中的反射系数将不再随单一的入射角变化而变化，并可拓展写成随频率和入射角共同变化的形式如下：
式中： 为地震频率，在实际应用的取值将参考地震主频来设定。作为描述速度随频率变化的重要关系，Aki 与 Richards[3]给出了近似常Q模型用于解释地震波速度在黏弹介质中的频散关系：
式中： 为虚数单位， 为参考频率， 与 分别为参考频率处的纵波速度与横波速度， 与 分别为纵波品质因子与横波品质因子。将式(4)与式(5)中的品质因子相关项进行泰勒展开并保留一阶项，可以得到：
将式(6)与式(7)代入式(3)中，可得式（8）。
式(8)作为本次研究的核心反射系数方程，定义了纵波反射系数与速度、密度、品质因子以及振幅和频率的之间的非线性关系。该方程能够同时兼顾地震振幅中的角度信息与频率信息，为后续开展叠前频变反演，实现碎屑岩储层流体识别提供了理论基础。与线性方程[20]相比，其直接给出了品质因子而非逆品质因子与反射系数之间的非线性关系，同时，方程中的独立参数为绝对值而非相对值，在实际获取地层纵横波速度、密度及地层品质因子过程中可避免参数转换带来的累积误差。
表1 黏弹性双层介质模型参数Table 1 Model parameters of viscoelastic double-layer media |
| 模型2 | /(m·s-1) | /(m·s-1) | /(kg·m-3) | QP | QS |
|---|---|---|---|---|---|
| 上层介质 | 3 140 | 1 560 | 2 285 | 97 | 28 |
| 下层介质 | 3 035 | 1 480 | 2 260 | 74 | 18 |
可以看出，无论是非线性方程还是线性方程，其反射系数都随入射角的改变而产生较为明显变化，这为AVO反演纵横波速度与密度奠定了有力的基础。相比随入射角度变化幅度，非线性方程和线性方程反射系数随频率的变化幅度有所减小，但仍然具有一定的振幅随频率变化的关系，该关系为反演地层品质因子提供了理论保障。对比非线性方程与线性方程反射系数与精确方程之间的误差绝对值能够发现，非线性方程与精确方程之间的误差绝对值更小，精度更高，这为准确预测地层纵横波速度、密度以及品质因子提供了理论支撑。
1.2 非线性频变地震反演
碎屑岩地层黏弹性的存在，使得反射系数方程的独立参数增加到了5个，为此，提高多参数反演的鲁棒性则是本次研究需要解决的关键问题。品质因子的非线性频变反演分为两个步骤：第一步，基于完全弹性假设，式(7)退化为式(1)，利用振幅随入射角变化关系、采用非线性反演算法开展叠前AVO反演获取纵波速度、横波速度以及密度；第二步，以第一步所得到的反演结果作为已知，利用反射系数方程实部的振幅随频率变化关系开展频变非线性反演，对地层的纵波品质因子和横波品质因子进行预测。反演的详细实施方案描述如下：
第一步AVO非线性反演纵横波速度与密度。基于贝叶斯框架，假设似然函数服从高斯分布，先验分布服从具有长尾特征的修正柯西分布，可以得到后验概率分布 的表达式如下：
式中： 为地震数据； 为待反演参数； 与 分别为服从高斯分布下的似然函数和服从修正柯西分布下的先验分布，详细表达式展开如下：
式中： 为地震数据的方差， 与 分别为待反演参数的协方差矩阵以及其拓展矩阵， 为待反演参数的初始模型， 为入射角 下的子波矩阵， 为该入射角下的反射系数。将式(10)—(11)带入式(9)，考虑 (入射角的个数)个角度入射的情况，进一步添加L1范数的正则化项,提高反演结果对含油气碎屑岩储层的边界刻画能力，化简可得第一步非线性AVO反演纵横波速度和密度的目标函数如下：
式中： 为待反演参数序列的长度， 与 分别为第一步反演的待反演参数向量（纵波速度、横波速度以及密度）以及其初始模型。
第二步非线性AVF反演纵横波品质因子。与第一步反演的目标函数建立过程相同，可以建立含L1范数正则化算子的第二步非线性频变反演的目标函数如下：
式中： 为反演所用到的不同主频地震记录的总数； 为第二步反演的待反演的纵波品质因子和横波品质因子的参数向量， 则表示第二步反演的待反演的纵波品质因子和横波品质因子的参数向量的初始模型。
观察公式(12)与(13)不难发现，其似然函数所建立的最小二乘项耦合了不同角度与不同频率分量下的地震记录与待反演参数之间的非线性关系，与此同时，L1范数正则化项的加入增加了目标函数的复杂程度，因此，常规牛顿类算法难以满足目标函数稳定求解的需要。为此，本次研究采用马尔科夫链蒙特卡洛（Markov chain Monte Carlo, McMC）非线性优化算法[30-32]，分别对两步反演的目标函数进行求解，获取目标函数取得极小值时所对应的待反演参数的数值解。McMC作为一种随机优化算法，能够对L1范数约束的目标函数进行有效求解，同时在理论上具有避免参数求解陷入局部极小的优点，保障了反演结果的鲁棒性与可靠性。
本次研究所述的两步方法虽然在第一步骤的反演中采用了一定的假设，忽略了品质因子对反射系数的影响，但是，考虑到多参数的求解会在较大程度上影响反演结果的稳定性，为此，两步反演方法通过分步反演的策略降低每一步骤的独立待反演参数个数，在保证反演结果准确性的同时，提高反演结果的稳定性。
2 合成地震记录测试
合成地震记录可由褶积模型得到。考虑到褶积模型的褶积滤波效应，在模型测试中，不同主频的合成地震记录由不同主频的地震子波与反射系数褶积得到。采用两步反演策略，可以分别输出得到如图3a与图3b中蓝色曲线所示的新两步反演方法下的无噪声以及信噪比=10的非线性反演结果。经过AVO一步直接反演方法的无噪声与信噪比=10的反演结果分别由图3a与图3b中的青色曲线所示。可以看出，在合成地震记录无噪声情况下，两种方法都能够合理、准确地对地层纵横波速度、密度以及品质因子进行预测。当合成地震记录信噪比=10时，两种方法的预测结果的准确性均有所下降，但预测结果仍然能够反映地层的弹性参数的变化情况，具有较好的准确性，这为所推导的反射系数方程的准确性与合理性提供了理论支撑。
对比不同方法的反演结果可以看出，两步反演方法由于每一步骤反演中的独立参数更少，且能够充分利用丰富的地震频率信息，因此其反演结果的精度更高，特别是合成地震记录含信噪比=10的随机噪声时，其品质因子的反演结果的抗噪性在较大程度上优于一步直接反演的结果，两步反演方法能够为多参数反演提供更强的稳定性。为量化不同反演方法下的反演结果与模型之间的误差，图4a与图4b分别计算了图3a与图3b中无噪声与信噪比=10的反演结果与模型之间的误差绝对值。观察发现，新两步反演方法误差绝对值小于常规一步反演方法，从量化的角度验证了本次研究所提出方法的合理性与稳定性。以无噪声合成地震记录作为观测数据输入，图5与图6分别展示了选用不同入射角度的频变反演结果。可以看出，由小角度入射时方程误差更小，且小角度合成地震数据具有更丰富的频带范围，因此，小角度入射的反演结果相比于大角度入射具有更高的精度，特别是120 ms和500 ms位置左右的含气砂岩达到低品质因子异常区域，小角度反演结果的误差小于大角度反演。该结果为新两步反演方法在实际应用中地震入射角度的选取提供了参考依据。
图4 不同反演策略下模型测试反演结果与岩石弹性参数模型误差绝对值Fig. 4 The absolute errors between the inversion results and the models of rock elastic parameters under different inversion strategies |
图5 不同入射角度下模型测试反演结果Fig. 5 Inversion results from model tests under different incident angles |
合成地震记录的测试结果充分验证了本文所提出的反演策略对提高反演结果鲁棒性方面的优势，为碎屑岩储层的高精度流体识别提供了理论依据。
3 实际工区应用
在合成地震记录测试验证的理论可行性的基础上，选取了中国东部某碎屑岩工区实际地震资料开展叠前地震反演，以测试方法的实际应用能力。研究区发育河流相碎屑岩储层，物源位于西北方向，砂体自西北向东南逐渐尖灭；岩性以砂岩与泥岩互层为主，自上而下发育多套含气砂岩储层，储层非均质性强，地震波在传播过程中产生较强的衰减作用，为反演地层品质因子并用于流体识别提供了数据基础。图7展示了所选碎屑岩工区主力含气储层位置的3个角度（8°、22°和31°）部分叠加地震剖面及W-2井的井旁地震道集。可以看出，在2.6 s左右的多套含气砂体发育处，地震振幅随入射角的增大呈现增大趋势，该现象符合含气砂岩的AVO类型，进一步验证了地震数据的可靠性。
以图7所示的3个角度部分叠加地震剖面作为观测数据输入，采用第一步AVO非线性反演，得到如图8所示的纵横波速度和密度反演剖面。可以看出，本文所提出的方法在第一步的反演结果能够较好地描述地层的横向快速变化与纵向弹性性质变化状态，具有较高的分辨率。对比探井W-2井与盲井W-1井旁道反演结果与测井解释结果，两者吻合度较高，反演结果具有较好的准确性与合理性，为第二步AVF反演地层品质因子给予了数据保障。第二步为AVF反演，采用固定入射角度、依靠地震振幅随频率变化关系对地层品质因子进行反演。考虑到本文所推导的反射系数方程在小角度入射条件下精度更高，同时，小角度部分叠加地震数据具有较宽的频带范围，因此第二阶段的AVF反演将采用小角度部分叠加地震数据。图9a、图9b与图9c分别展示了图7a中的小角度部分叠加地震剖面的10 Hz、25 Hz及40 Hz主频处的分频后地震剖面。地震剖面中，随着地震主频的增加地震同相轴逐渐变细，符合物理学的认知。以图9所示的不同主频的地震剖面作为观测数据输入，图8所示的纵横波速度和密度反演结果作为先验已知，采用第二步AVF非线性反演，输出如图10a与图10b所示的纵波品质因子与横波品质因子反演结果。观察品质因子反演结果可以看出，纵横波品质因子均能较好地描述不同衰减特征地层的横向分布与纵向变化，具有较高的分辨率。对比探井W-2井与盲井W-1井旁道反演结果与实际测井品质因子解释结果可以看出，依赖于较高精度的非线性反射系数方程以及两步反演策略，井旁道反演结果与测井解释结果具有较高的吻合度，反演结果具有较高的准确性与合理性。考虑到含气碎屑岩储层具有较强的衰减特征，为此，利用纵横波品质因子的低值异常能够有效刻画含气储层的纵向分布与横向发育状态。观察图10中的品质因子反演结果，其低值异常能够对自上而下发育的多套含气砂岩进行有效刻画，刻画结果分辨率较高、横向变化较为剧烈，符合碎屑岩工区含气砂岩储层分布的地质认识。
图8 实际地震资料新方法第一步AVO反演结果Fig. 8 The first step's AVO inversion results for actual seismic data using the new method |
图9 小角度叠加地震剖面不同频带地震分量Fig. 9 The seismic components of different frequency bands in seismic profiles at small angles |
为进一步对比两步反演方法相较于一步直接AVO反演方法的优势，图11展示了基于一步AVO直接反演方法的5个待反演参数（纵波速度、横波速度、密度、纵波品质因子以及横波品质因子）的实际碎屑岩工区反演结果。可以看出，一步AVO直接反演结果能够较为合理、准确地对地层弹性参数和品质因子的纵向变化与横向分布进行了表征，其井旁道反演结果与实际测井解释结果具有一定的吻合程度，验证了本次研究所推导的反射系数方程的精度。但是，一步AVO直接反演受制于多参数直接反演稳定性较低以及未能利用地震数据丰富的频率信息，其反演结果的准确性低于本文所提出的两步反演方法，且对于含气砂岩的刻画能力较差，特别是图10与图11中反演结果的黑色椭圆所圈出的区域流体识别效果较差。
实际工区的应用结果验证了本文提出的方法的适用性和可靠性，为碎屑岩工区储层流体识别提供了新的方案。值得注意的是，该方法旨在定性层面阐述了含气储层与品质因子之间的关系，未能有效明确品质因子与孔隙流体性质（气/水）的定量关系，该不足也是未来进一步研究的重点方向。
4 结论
为提高碎屑岩储层流体识别精度，本次研究基于Zoeppritz方程的Bortfeld非线性近似，推导了黏弹性介质下的含纵横波品质因子的非线性反射系数方程，方程中的反射系数能够同时随入射角和地震频率变化，为开展频变反演获取碎屑岩地层品质因子提供了方程基础。
本文提出了一种两步频变非线性反演策略，相比于一步AVO直接反演，该反演策略提高了多参数优化的稳定性，为提高品质因子预测精度及碎屑岩储层流体预测可靠性提供了技术保障。
本次研究所提出的非线性频变反演方法能够对中国东部某实际工区的多套含气砂岩储层进行较为准确的预测，结果符合地质认识，具有较高的分辨率，为碎屑岩含气储层的横向展布刻画提供了新的方案。
