混合概率模型驱动的叠前地震反演方法

李 坤 印兴耀

（①中国石油大学（华东）地球科学与技术学院，山东青岛266580；②海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛266071）

0 引言

叠前地震反问题作为储层地球物理学的核心内容，是解决复杂油气储层特征描述及油气判识的重要理论与方法，通常包括波动方程反演（层析成像、全波形反演等）、AVO（Amplitude Variation with Offsets，振幅随炮检距变化）、EI（Elastic Impedance，弹性阻抗）反演［1-4］。然而，由于带限地震数据缺失低频和高频信息，且受到噪声的干扰，叠前地震反问题表现出较强的不适定性，因此叠前地震反问题必须通过其他途径（如钻、测井数据、油藏建模数据等）补充缺失的模型信息。正则化和概率化反演（以贝叶斯理论为主）是缓解叠前地震反问题不适定性最重要的两种手段［5］。然而，随着油气藏类型日益复杂化，仅仅依靠单一的模型最优解，已经难以满足地球物理勘探人员对复杂油气储层岩性和含流体性质的精细描述的要求，故模型参数最优解、不确定度及置信区间等勘探风险评估逐渐引起储层地球物理学家的关注。

在地震反演的贝叶斯框架中，模型参数的先验信息可以通过随机变量的先验概率密度分布函数（Prior Probability Density Function，Prior PDF）定量表征。其认为地下介质模型参数是随机的，则地震反问题是指在多类型观测数据集的协同约束下匹配随机变量的随机实现过程［5］。目前，以贝叶斯推断为基础的反问题中，最大后验概率解（Maximum A Posteriori Estimation，MAP）［5-13］、后验PDF的显式解［14-20］及后验PDF的随机解（以马尔可夫链蒙特卡洛模型为主，MCMC）［21-27］是应用最广泛的三种概率解。贝叶斯MAP解是通过贝叶斯公式计算得到模型参数的后验概率密度分布，以最大化后验PDF获取模型参数的最优解。后验概率的显式解是指模型参数的后验统计特征（如均值和协方差等）可以通过显式表达式定量表征，以高斯型概率模型为应用典型。随机解是指通过随机采样的方式对后验PDF进行随机采样，利用大量随机样本点逼近理论后验PDF的求解方式。Hansen等［28］阐述了两类观测数据约束下高斯概率密度模型的地质统计模拟问题以及贝叶斯线性反问题的显式解。Grana等［14-15］在混合先验分布的基础上，采用序贯模拟方法对后验PDF的显式解进行单点模拟，实现了纵波阻抗、离散岩相的同步预测，该方法由于未引入模型参数的低频信息，需要在模拟相对纵波阻抗后补偿低频背景。Lang等［16］在高斯概率模型的基础上，通过后验概率模型的显式解实现了叠前地震纵、横波阻抗和密度的随机反演方法。Li等［27］联合马尔科夫链蒙特卡洛模型、混合高斯概率先验信息，提出了岩相驱动下的叠后地震随机反演方法，改善了叠后地震反演分辨率的同时，实现了“离散相态”与“弹性参数”的协同预测。

Hastie等［29］基于贝叶斯推断理论框架，系统阐述了混合概率模型在“离散变量”反演与分类领域的意义。混合概率模型在贝叶斯推断中是实现叠前地震“连续弹性参数”和“离散相态”协同预测的基础。基于此，本文依据岩性驱动的模型参数统计先验信息，把地层岩性（如砂、泥岩）与模型参数和观测地震数据联系在一起，提出了一种基于混合概率模型的时频联合域叠前地震概率化弹性阻抗反演方法；综合利用时域、频域地震、低频整合先验信息及已知模型点四类数据集实现对弹性阻抗、离散岩相的序贯协同预测；此外，将非线性边界修正算法引入叠前地震流体因子的直接提取过程，旨在提高叠前地震反演的稳定性，克服弹性参数提取过程中易出现“超界解”的问题。通过理论模型测试和实际勘探实例证实叠前地震概率化弹性阻抗反演方法的稳定性和可行性。

1 方法原理

1.1 时频联合域地震正演模型构建

不同变换域地震反演是通过考虑不同域内地震响应与合成地震数据间的匹配程度以搜索最优模型参数的过程［5］。频率域反演相比时间域反演具有频率分量自动解耦、频率选择自由、多分量逐级寻优及分辨率高的特点。因此，将混合域卷积模型应用于地震数据的正演

式中：S（ω，θi）、W（ω，θi）分别为入射角为θi的地震频谱和子波频谱，ω为角频率；rPP（t，θi）为基于Biot-Gassmann孔隙弹性理论的Russell线性AVO近似方程［30-33］；j为虚数单位；t为时间；N（ω）为噪声频谱。式（1）可以重组为矩阵形式

假设η为利用实际测井数据构建的弹性阻抗低频先验数据，λL为低频背景的正则化因子，D为正则化对角矩阵。联合式（3）与η，可将正演模型改写为

式中：H为输入的三种类型已知数据集；P为正演核矩阵；Cωω、Ctt和Cηη分别是和η的协方差矩阵；Cωt、Cωη分别是、η之间的协方差矩阵；Ctη是与η之间的协方差矩阵；CH为三种数据集H的联合协方差矩阵。

1.2 混合概率模型驱动的弹性阻抗反演方法

混合概率模型是多个“单峰”PDF的线性叠加，旨在量化待反演模型参数的“多峰”分布特征，是“单峰”PDF的推广形式［5］。待反演模型参数的先验概率密度分布往往受到储层岩性的影响，储层岩性不同将导致模型参数的统计特征差异［5，15，19，27］。具有“多峰”分布特征的混合概率模型是实现“离散型”—“连续型”模型参数协同预测的关键。高斯型PDF是数学、物理及工程等领域中最普遍的统计分布，具有线性变换性质，易于求解地球物理反问题的显式解［16-18，28-29］。因此，本文将高斯混合PDF作为弹性阻抗反演模型的先验PDF。其中，弹性阻抗的先验均值、方差及高斯混合PDF中的峰态数量由实际测井数据的解释结果和统计特征获得。

针对贝叶斯线性地震反问题，若目标工区中储层岩性类别的数量为M，将式（6）代入式（A-8）～式（A-9）（附录A），则模型最优解可用后验概率分布的显式解进行量化，推导得到三类数据集（即和η）协同约束下的第k个高斯概率分量pk（m|H）的显式解，其后验均值为

式中：λk＝p（zk）为第k种离散岩性zk的先验概率；表示均值为、协方差为的高斯PDF，分别为已知数据集H的训练均值和训练协方差。针对实际工区可以通过测井解释数据统计获取，其中Fk是第k种岩性出现的频数。

通过混合后验均值和后验协方差的显式表达式（式（7）、式（8）），即可计算三类已知数据集H约束下的待反演模型参数最优解。为了更好地分析待反演模型的不确定性，采用基于单点实现的序贯模拟算法（Sequential Simulation）（图1）对后验PDF进行随机采样。此外，将先前的模拟点集ms作为贝叶斯地震反演中下一个模拟点mi的约束条件。因此，针对第k个高斯概率分量pk（m|H），在待模拟点位置的后验均值可改写为

图1 利用序贯模拟算法的时频联合域地震反演示意图［16-17］

在贝叶斯推断中，混合先验概率模型搭建了储层弹性阻抗与离散岩性之间的关系，在弹性阻抗、观测数据等特征属性的协同约束下，即可获得储层岩性的后验概率密度分布；然后，依据最大后验概率解获取储层岩性的最优解。基于此，从被优选的目标分量中随机采样得到岩相约束下的弹性阻抗数据。Hansen等［28］和Lang等［16］提出了2类条件数据约束下的地震反问题的线性解。基于贝叶斯分类，可以得到多类数据集协同约束的第k个高斯概率分量的后验权重

为序贯模拟过程中条件数据的先验均值。待模拟点mi的可以用来判断、识别离散的储层岩性。当离散岩性分类（优选混合概率模型中的高斯分量）之后，待反演模型参数（即弹性阻抗）将从被选择的高斯分量中随机采样获得。值得注意的是，当访问下一个模拟点时，先前的模拟点被视为已知数据，且已知数据的邻域（矩形或圆形邻域）选择对于提高地震序贯模拟的计算效率具有重要意义，邻域半径通常设置为一个波长至两个波长之间。

1.3 叠前地震Gassmann流体参数反演方法

复杂油气储层的饱和岩石为干岩石骨架与多种流体填充物构成的多孔双相介质，宏观岩石模量受干岩石骨架和孔隙流体的综合影响。Russell等［30-31］基于Biot-Gassmann孔隙弹性理论推导了岩石饱含流体情况下的线性地震AVO反射系数方程

式中：θ为入射角；f是Gassmann流体项，反映岩石孔隙流体的作用；μ为岩石剪切模量，反映岩石骨架的作用；ρ为饱和岩石密度；为干岩石骨架的纵横波速度平方比；是饱和岩石的纵横波速度平方比；Δlnf、Δlnμ和Δlnρ分别为Gassmann流体项、剪切模量和密度的地层反射率。

由式（15）可知，弹性阻抗反演比叠前地震AVO反演具有更高的计算效率和稳定性［1，4］，因此提出了基于模型正则化和边界约束的改进弹性阻抗反演方法，以预测表征油气储层含流体性质的Gassmann流体项等弹性参数。与式（15）对应的流体弹性阻抗方程可以改写为

式中：F、μ和ρ分别表示井旁道实测Gassmann流体项、剪切模量和密度构成的向量；n为重采样后的测井数据样点数。

采用阻尼最小二乘法（或共轭梯度法）求解式（18），获得a′（θi）、b′（θi）和c′（θi）。令弹性参数的相对值表示为Rf、Rμ、Rρ和REI，即

式中N为单道地震数据的样点数。为了实现非线性方程的线性化，对式（16）等号两边取自然对数，并将式（19）代入，得

式中：Aθi＝diag［a（θi）］N；Bθi＝diag［b（θi）］N；Cθi＝diag［c（θi）］N；diag［·］表示对角矩阵运算；Q为入射角的数量。

考虑到Gassmann流体项、剪切模量及密度的大尺度先验背景可以从实际测井数据中获得，故将弹性参数的已知先验信息用L2范数表示，并引入到Gassmann流体项等参数的反演中，提高弹性参数反演的可靠性，则弹性参数反演的目标泛函和更新量分别为

式中：Δξf、Δξμ和Δξρ均为先验数据的迭代残差；先验信息的正则化算子Df、Dμ、Dρ是单位矩阵；ξf、ξμ和ξρ是弹性参数的先验背景；λf、λμ和λρ分别是Gassmann流体项、剪切模量和密度的先验约束权重。式（21）～式（22）未考虑待反演模型的下限和上限信息，然而实际地震数据受到噪声的干扰，反演结果往往会超出模型边界的解（固有不稳定性引起）。为了改善该问题，引入非线性边界转化算法。假定当前迭代次数为y，则中间变量的更新量为

2 理论模型测试

为了验证混合概率模型驱动的叠前地震概率化反演可行性，利用重采样后的测井数据设计了一维理论模型（图2）。图2展示了Gassmann流体项、剪切模量、密度、小角度弹性阻抗的先验概率密度分布及理论离散岩相。采用30Hz零相位Ricker小波、2ms采样间隔和精确Zoeppritz方程，合成了信噪比S/N＝1情况下的AVO地震数据，AVO地震数据的入射角分别为10°、20°和30°（图3）。

第一步是基于贝叶斯混合概率模型的地震弹性阻抗反演处理。不同角度（10°、20°和30°）的弹性阻抗与离散岩相的估计结果如图4～图6所示。从图4～图6可以看出，50次模拟实现（灰线）的均值解（红线）与理论模型（绿线）的拟合度较高。不同角度情况下，50次离散岩相判识的最大后验概率解（图4c、图5c和图6c的左半部分）均与理论岩相（图4c、图5c和图6c的右半部分）具有高度一致性。与之对应的，50次模拟实现的岩相误差主要出现在地层界面上及先验均值之间的模型参数过渡区，该现象产生的根本原因是这些区域的模型参数在不同后验PDF中的概率是相近的。因此，基于贝叶斯混合概率模型的地震弹性阻抗反演方法可以有效地实现弹性阻抗和离散岩相的概率化预测，为后续的边界约束叠前地震概率化弹性阻抗反演奠定了数据基础。

图2 实际测井数据

图3 S/N＝1情况下入射角为10°（a）、20°（b）和30°（c）的合成地震数据

图4 S/N＝1情况下小角度（10°）弹性阻抗与离散岩相的估计结果

图5 S/N＝1情况下中角度（20°）弹性阻抗与离散岩相的估计结果

图6 S/N＝1情况下大角度（30°）弹性阻抗与离散岩相的估计结果

第二步是基于弹性阻抗的Gassmann流体项、剪切模量等弹性参数的同步反演。图7展示了S/N＝1情况下反演得到的f、μ和ρ。由图可见，本方法的50次模拟实现（灰线）的弹性参数预测结果被有效地限制在上、下限（棕色虚线）范围内，与理论模型高度一致；f预测结果的相对误差为12.8%，μ的相对误差为17.9%，证明了该边界约束弹性阻抗反演算法在参数提取方面的可行性；存在噪声干扰的情况下，f的估计结果受随机噪声的影响程度较小，μ的反演精度比f低，与理论分析相一致，主要是由于μ对地震AVO反射系数的贡献度低导致的。

为了阐述基于弹性阻抗的Gassmann流体项、剪切模量等弹性参数的同步反演方法与常规反演方法之间的差异，图8展示了不同反演方法得到的Gassmann流体项反射率的频谱计算结果。由图可见，最小二乘确定性反演方法仅能获取模型参数的“最优解”，反演结果的频带范围主要与输入地震数据的频宽有关，高频信息不能得到有效恢复（图8b）；本文方法反演结果的频谱（图8c、图8d）比常规最小二乘反演方法具有更宽的频带，尤其是在模型参数的高频信息恢复方面，该概率化反演方法可以更有效地恢复地下介质高频信息（图8d）。

图7 S/N＝1情况下50次随机模拟实现的弹性参数及其相对误差估计结果

图8 S/N＝1情况下不同反演方法得到的Gassmann流体项反射率的频谱分析对比

3 实际资料处理

以中国华南地区Well-1井含油储层的地震勘探为例，验证本文混合概率模型驱动的叠前地震弹性阻抗反演方法的实用性。

研究区砂岩储层主要为含油层与油水同层。由图9可以看出，小角度弹性阻抗可以较好地区分优质砂岩储层（含油砂岩、含油水砂岩，图9a椭圆框）与泥岩（干层）（图9a）；相比中角度、大角度弹性阻抗，小角度弹性阻抗可以更好地区分砂岩储层（图9b）。不同岩性情况下，小角度弹性阻抗的统计特征（先验均值、协方差）差异更大。然而，中角度与大角度弹性阻抗在砂岩储层和泥岩干层位置的先验均值却比较接近，故在概率化推断过程中将存在较强的不确定性。由于小、中、大角度弹性阻抗具有相似的统计特征“拐点”，因此本文将弹性阻抗岩相判识结果的“交集”作为目标工区岩相预测的依据。

图9 弹性阻抗的岩石物理参数交会与统计

图10 储层含流体性的岩石物理参数交会与统计

图11 小角度（5°～16°）叠加的地震剖面、弹性阻抗及离散岩相的反演结果

由图10可见，Gassmann流体项较好地区分了油层、油水同层及泥岩干层，同时Gassmann流体项、密度等参数服从混合高斯概率密度分布，表现为“多峰”的概率统计特征。基于此，可将Russell近似公式作为“流体相”驱动叠前地震反演的地球物理映射关系，即模型参数为表征孔隙流体的Gassmann流体项、干岩石骨架的剪切模量、岩石密度及孔隙含流体类型（流体相）。

由图11～图13可以看出：①由于概率化地震反演的不确定性，单次模拟结果中仍存在序贯随机采样的干扰（图11c～图11d、图12c～图12d、图13c～图13d）。因此，地震概率化反演往往需要多次随机实现的均解或最大条件概率解（图11d、图12d、图13d）。对比图11c和图11d，10次实现的均值解可以有效地抑制单次实现中的随机干扰。②井旁道的弹性阻抗估计结果（红线）和砂泥岩模型的岩性判识结果与实际测井弹性阻抗计算结果、砂岩储层解释结果高度一致，可以较准确地识别出砂岩油藏的发育位置（图11e～图11h、图12e～图12h、图13e～图13h）。

根据井旁道弹性阻抗的岩石物理参数交会与统计特征分析（图9），角度弹性阻抗可以有效地识别砂岩储层，且将小、中、大3个角度弹性阻抗的岩相判识结果的“交集”作为砂泥岩中储层的最终判识结果（图14a），可以降低岩相判识的不确定性，与实际砂岩油层解释结果吻合较好。图14b展示了边界约束叠前弹性阻抗反演预测的Gassmann流体项f，从图中可以看出，Gassmann流体项在油砂储层位置存在明显低值异常，且砂泥岩－岩相与流体识别结果在预测位置上保持了较高的一致性，证明了本文混合概率模型驱动的叠前地震概率化反演方法在岩性预测和储层流体识别中具较好的实用性，应用前景广泛。

图12 中角度（16°～26°）叠加的地震剖面、弹性阻抗及离散岩相的反演结果

图13 大角度（26°～30°）叠加的地震剖面、弹性阻抗及离散岩相的反演结果

图14 叠前地震概率化弹性阻抗反演的岩相和Gassmann流体项f的估算结果

4 结论

本文依据地下介质待反演模型参数的先验概率服从混合型概率密度模型，提出了基于时域、频域地震、低频整合先验信息及已知模型数据点四类条件数据集协同约束下的混合概率模型驱动的叠前地震时频联合域弹性阻抗反演方法。主要结论如下。

（1）混合概率模型是实现“离散型”—“连续型”模型参数协同预测的关键，只有当不同离散相的待反演参数具有不同的统计特征时，才能够实现地震反演的离散相态判识。高斯混合PDF与高斯PDF的乘积仍可以表征为标准的混合高斯PDF，此时，“连续弹性参数”、“离散相态”的后验均值与后验协方差具有显式表达式，易于采用序贯模拟随机采样算法求解。

（2）相比常规确定性反演方法，序贯模拟算法是通过逐点遍历的方式实现模型参数的随机采样，先前模拟点将作为下一个模拟点的约束条件，有助于实现“连续弹性参数”与“离散岩相”的概率化协同判识，便于实现和理解，使多数据协同约束的地球物理反问题变得清楚易懂。

（3）理论测试和实际应用验证了该方法的可行性及稳定性。应用实例表明，混合概率模型驱动的叠前地震弹性阻抗反演的岩相和概率多重实现，与实际测井数据和油藏解释结果具有较高的一致性，有助于复杂岩性油气藏的储层刻画和流体直接检测，应用前景广泛。

附录A

贝叶斯理论是在当前观测样本数据的情况下，结合对模型参数的先验认知，计算得到模型参数后验概率密度分布（PDF）的统计学习方法。

式中：p（m|H）表示后验概率密度分布；p（H|m）为观测数据的似然概率密度分布；m为表征该概率模型的未知参数；p（m）被称为未知参数的先验信息；p（H）为已知观测数据的PDF，与θ的取值无关［34］。若待反演模型参数服从高斯混合PDF

式中：M为高斯混合概率模型中所包含概率分量的数量；；n为单道地震数据的样点数；为模型均值；为模型参数协方差。在高斯似然函数的前提下

式中：m为被选择频率分量的个数；H为观测地震数据；P为正演核矩阵；CH为协方差矩阵。利用贝叶斯公式（式（A-1））推导后验PDF（贝叶斯公式的分母与待估计模型参数无关，可忽略分母项）

将式（A-6）改写为标准的高斯型PDF形式

则式（A-7）所示的第k个标准型后验PDF分量的后验均值的显式解为

式（A-8）、式（A-9）即描述了高斯混合概率模型驱动下的标准后验PDF的显式解。