CN111948712A - 一种基于深度域地震记录的叠前线性反演方法 - Google Patents

Abstract

本发明涉及一种基于深度域地震记录的叠前线性反演方法，包括：步骤101：进行深度域井震标定；步骤102：对整个工区的实际深度域地震记录进行层位拾取并与深度构造解释结果结合，然后对深度域测井数据进行插值，并将沉积模式作为协同数据参与插值，最终获取整个工区的初始低频弹性参数数据体；步骤103：首先对整个工区的时间域地震记录进行时间域地震子波提取，然后基于点传播函数与时间域地震子波的对应关系进一步将时间转换为深度从而获取整个工区的点传播函数；步骤104：利用获取的点传播函数构建对应的非稳态褶积矩阵形式；步骤105：将非稳态褶积子波矩阵与反射系数线性近似公式结合，并利用贝叶斯推理线性反演方法进行整个工区的弹性参数反演。

Description

一种基于深度域地震记录的叠前线性反演方法

技术领域

本发明涉及一种叠前线性反演方法，具体是关于一种利用深度域地震记录直接进行叠前线性反演的方法。

背景技术

地震记录往往是时间域的，但仅从时间域地震记录解释地质问题在理论上存在不足且并未直接提供相关地质问题解释，主要是因为：1)从数学角度来看，时间域的地震记录与深度域的测井数据在坐标尺度上存在不统一，因此需要通过井震标定确定时间与深度关系。而在井震标定中，强烈的人机交互，如平移或拉伸等操作使得时深关系不确定性大大增强。2)由时间域地震记录获取的地震属性往往也是时间域的，并未与相关的地质问题直接联系，仍需要进行时深转换。随着计算机技术的发展，计算成本降低，有助于深度偏移算法高效运行，因此地质学家和地球物理学家将目光转向深度域偏移后的地震记录。

目前，深度域反演大致分为下列几类方法：1)通过数次时深转换将深度域地震记录转换到时间域再进行反演，但是几次时深转换使得有效信息大量丢失，甚至丢失反映薄层的地震信息，最终使得反演误差放大；2)忽略深度域子波的时变性，利用平均速度计算深度域子波，直接利用线性褶积模型直接进行地震记录反演，但该方法忽略深度域子波的时变性，在长时间窗内反演会造成估计的弹性参数准确度不高；3)利用滑动的窗口进行反演，在窗口内忽略子波的时变性，是方法2)的改良，使得弹性参数的准确度获得提升，但并未从根本上解决子波时变带来的反演误差大的问题；4)借助深度域偏移，将时间域合成的地震记录转换到深度域，再进行深度域反演，直接将偏移介入反演中，使得计算效率大大降低；5)利用一个等效速度将地震记录转换到满足褶积模型的“伪深度”域，进行线性褶积反演，再利用该等效速度将估计的弹性参数转换到真深度域，从而完成反演整个流程，但该方法仍然无法避免两次不同域之间的转换，使得反演误差放大；6)基于非确定性映射关系，利用统计原理和机器学习的模糊深度域反演方法，但如何提高标签数据的质量、如何扩充标签数据数量以及神经网络结构设计、效率等问题仍待解决，使得该类方法在工业界并未得到广泛应用。

综上所述，目前深度域反演方法研究主要存在以下问题：①理论研究不够充分，无法将深度域成像与反演有机的结合，未实现真正意义上的深度域反演；②忽略子波的时变性，采用长时窗或滑动时窗进行反演，均会使得反演结果准确度降低；③基于转换域的反演方法，其转换域物理意义不明确，且无法避免转换过程中的误差累计；④基于深度偏移进行时深转换的反演方法，计算效率较低，且转换过程有精度损失；⑤基于非确定性映射关系进行反演方法，由于标签数据的质量、数量以及神经网络结构设计、效率等问题，很难将其应用到实际生产应用中。因此亟需一种理论完善、计算效率高且可以直接对深度域地震记录进行反演的新方法。

发明内容

针对上述问题，本发明的目的是提供一种利用深度域地震记录直接进行叠前线性反演的方法，该方法利用深度域偏移成像理论建立点传播函数与深度域子波直接的非稳态褶积关系，并推导非稳态褶积的矩阵形式，最终与现有的线性反演理论结合实现深度域反演，有效地提高所反演的弹性参数精度，满足储层预测与流体识别要求。

为实现上述目的，本发明采取以下技术方案：一种基于深度域地震记录的叠前线性反演方法，包括：

步骤101：利用基于点传播非稳态褶积运算，对深度域测井数据进行深度域正演模拟得到深度域地震记录，将模拟的深度域地震记录与实际深度域地震记录进行平移匹配，从而获得深度域测井数据在地震尺度上的准确位置，即得到准确匹配的深度域测井数据；

步骤102：首先对整个工区的实际深度域地震记录进行层位拾取，然后将拾取的层位信息与深度构造解释结果结合，对步骤101中已经准确匹配的深度域测井数据进行插值，并将沉积模式作为协同数据参与插值，最终获取整个工区的初始低频弹性参数数据体；

步骤103：首先对整个工区的时间域地震记录进行时间域地震子波提取，然后基于点传播函数与时间域地震子波的对应关系进一步将时间转换为深度从而获取整个工区的点传播函数；

步骤104：利用获取的点传播函数构建对应的非稳态褶积矩阵形式；

步骤105：将非稳态褶积子波矩阵与反射系数线性近似公式结合，并利用贝叶斯推理线性反演方法进行整个工区的弹性参数反演。

所述的叠前线性反演方法，优选地，在所述步骤101中，实际深度域地震记录即成像结果：

s＝M*d (1)

式中，M表示深度偏移算子；d表示偏移前时间域地震记录；

而偏移前时间域地震记录满足常规稳态褶积运算，其矩阵形式表达为：

d＝W*R (2)

式中，W表示常规稳态褶积运算构建的子波矩阵；R为反射系数；

由式(1)和式(2)得：

s＝(M*W)*R＝H*R (3)

式中，H表示正演运算与偏移运算结合的矩阵算子。

所述的叠前线性反演方法，优选地，在所述步骤103中，一维速度模型下，点传播函数与时间域地震子波满足以下对应关系：

式中，T表示时间域地震子波的周期；

表示点传播函数在一个波长λ内的慢度；h表示点传播函数的深度；v表示速度；

在均匀速度介质中，考虑双程旅行时，式(4)改写为：

式中，f＝1/T表示频率；k＝1/λ为波速。

所述的叠前线性反演方法，优选地，在所述步骤104具体为：将步骤103获取的整个工区的点传播函数利用对应的非稳态褶积矩阵形式进行排列，并在矩阵边界进行对应的点传播函数截取，以保持子波矩阵大小与常规稳态褶积一致。

所述的叠前线性反演方法，优选地，非稳态褶积表示为：

式中，s(t)表示地震记录；w(t,τ)表示变化的子波；r(t)为反射系数的向量形式；将式(6)转换为如下离散形式：

式中，s(i)、w(i,j)和r(i)分别表示地震记录、子波和反射系数的离散形式；i和j表示离散形式下的深度索引；

其中，将反射系数r(i)利用单位脉冲函数δ(i)分解为：

将式(8)代入式(7)得：

与式(2)中常规稳态褶积运算相同，此时非稳态褶积具备同样的矩阵形式：

式中，n表示子波序列的长度；m表示反射系数序列的长度；

在常规稳态褶积中：

其中，W表示稳态褶积子波矩阵，

表示非稳态褶积子波矩阵，其子波元素的上标表示不同深度下的点传播函数；此时，先通过时间域地震子波求取点传播函数，再按照稳态褶积子波矩阵对应位置截取子波元素以构建非稳态褶积子波矩阵。

所述的叠前线性反演方法，优选地，在所述步骤105具体包括：

①将非稳态褶积子波矩阵与反射系数线性近似公式结合，得到深度域正演矩阵为：

式中，s表示观测数据向量，即实际深度域地震记录；m为模型参数向量，即待反演的弹性参数；F为线性反射系数矩阵；

式(12)简记为：

s＝Gm (13)

式中，G表示深度域正演矩阵；

②基于贝叶斯推理框架构建反演目标函数为：

式中，σ_d表示数据的方差；R(m)表示模型的先验约束；

③对反演目标函数求导可得：

式中，

μ表示模型均值可用初始模型代替；当选取高斯分布为先验分布时，

C_m表示弹性参数的协方差矩阵；

④令导数等于0，此时式(15)整理得到：

式中，μ_h表示先验信息权重系数，且

⑤对于式(16)采用迭代重加权最小二乘算法求解最终得到准确的弹性参数，该弹性参数即可用于储层预测与流体的准确识别。

本发明由于采取以上技术方案，其具有以下优点：1、本发明直接在深度域进行测井数据与地震记录的匹配，不需要获取时深转换关系，简化了井震标定这一步骤。2、本发明不需要进行时深转换，避免了转换误差累积，使得反演结果准确度提高。3、本发明基于完善点传播函数偏移成像理论，整个反演过程并未忽略相应的物理意义。4、本发明利用时间域子波与低频纵波速度获取随深度变化的点传播函数的，降低了点传播函数获取难度。5、本发明可以与现有的线性反演理论结合，计算效率高，且获取的弹性参数稳定、连续，满足储层预测与流体识别的要求。6、从程序化角度出发，本发明在现有成熟的时间域线性反演的基础上修改，替代相应的褶积矩阵即可进行相应计算，因此编码易于实现，满足工业需求。

附图说明

为了更清楚的说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍，下面描述中的附图仅仅是本发明的一些实施例。

图1是本发明叠前线性反演方法的流程图；

图2(a)是时间域稳态褶积矩阵的示意图；

图2(b)是深度域非稳态褶积点传播函数矩阵的示意图；

图3是本发明实施例合成的深度域地震记录的示意图；

图4是本发明实施例无噪声情况下采用深度域叠前线性反演得到的纵波速度、横波速度、密度的示意图；

图5是本发明实施例信噪比为3:1情况下采用深度域叠前线性反演得到的纵波速度、横波速度、密度的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明实施例提供的基于深度域地震记录的叠前线性反演方法，具体包括以下步骤：

步骤101：测井数据与地震记录的匹配即井震标定，是测井数据作为硬数据参与地震反演的关键。本发明利用基于点传播非稳态褶积运算，对深度域测井数据进行深度域正演模拟得到深度域地震记录，将模拟的深度域地震记录与实际深度域地震记录进行平移匹配，从而获得深度域测井数据在地震尺度上的准确位置，即得到准确匹配的深度域测井数据，以为后续的操作提供技术支持。由于该步骤均在深度域进行，仅需要稍微的平移操作即可，因此可以避免常规井震标定中拉伸操作。

其中，实际深度域地震记录即成像结果：

s＝M*d (1)

式中，M表示深度偏移算子；d表示偏移前时间域地震记录。

而偏移前时间域地震记录满足常规稳态褶积运算，其矩阵形式可以表达为：

d＝W*R (2)

式中，W表示常规稳态褶积运算构建的子波矩阵；R为反射系数。

由式(1)和式(2)可得：

s＝(M*W)*R＝H*R (3)

式中，H表示正演运算与偏移运算结合的矩阵算子，可反映速度变化(与速度有关)，该运算操作为任意散射点正演和偏移运算后的点传播函数脉冲响应。因此，矩阵H可以由离散后的点传播函数构建。可见，深度域地震记录是点传播函数与反射系数的非稳态褶积(即点传播函数随速度变化)。

步骤102：首先对整个工区的实际深度域地震记录进行层位拾取，由于层位拾取的人机交互问题，可利用机器学习进行训练，减少人机交互时间；然后将拾取的层位信息与深度构造解释结果结合，对步骤101中已经准确匹配的深度域测井数据进行插值(插值方法可利用距离关系的克里金插值)，并将沉积模式(即对一定环境中的现代沉积物的物理、化学、生物特征综合研究的基础上概括出的沉积环境及其沉积物的物化模型)作为协同数据参与插值，最终获取整个工区的初始低频弹性参数数据体。

步骤103：经测试表明低频纵波速度可以代替真实纵波速度计算点传播函数，因此可利用步骤102获取的初始低频弹性参数数据体获取点传播函数。但在此之前需先对整个工区的时间域地震记录进行时间域地震子波提取，然后基于点传播函数与时间域地震子波的对应关系进一步将时间转换为深度从而获取整个工区的点传播函数。

其中，一维速度模型下，点传播函数与时间域地震子波满足以下对应关系：

式中，T表示时间域地震子波的周期；

表示点传播函数在一个波长λ内的慢度；h表示点传播函数的深度；v表示速度。

在均匀速度介质中，考虑双程旅行时，式(4)可以改写为：

式中，f＝1/T表示频率；k＝1/λ为波速。

步骤104：将步骤103获取的整个工区的点传播函数利用对应的非稳态褶积矩阵形式进行排列，并在矩阵边界进行对应的点传播函数截取，以保持子波矩阵大小与常规稳态褶积一致。

其中，非稳态褶积表示为：

式中，s(t)表示地震记录；w(t,τ)表示变化的子波；r(t)为反射系数的向量形式。

将式(6)转换为如下离散形式：

式中，s(i)、w(i,j)和r(i)分别表示地震记录、子波和反射系数的离散形式；i和j表示离散形式下的深度索引。

其中，将反射系数r(i)利用单位脉冲函数δ(i)分解为：

将式(8)代入式(7)可得：

与式(2)中常规稳态褶积运算相同，此时非稳态褶积具备同样的矩阵形式：

式中，n表示子波序列的长度；m表示反射系数序列的长度。

在常规稳态褶积中：

其中，W表示稳态褶积子波矩阵，

表示非稳态褶积子波矩阵，其子波元素的上标表示不同深度下的点传播函数，因此构建非稳态褶积子波矩阵时需要先通过时间域地震子波求取点传播函数，再按照稳态褶积子波矩阵对应位置截取子波元素以构建非稳态褶积子波矩阵。

步骤105：将非稳态褶积子波矩阵与反射系数线性近似公式结合，并利用贝叶斯推理线性反演方法进行整个工区的弹性参数反演，具体包括：

①将非稳态褶积子波矩阵与反射系数线性近似公式结合，得到深度域正演矩阵为：

式中，s表示观测数据向量，即实际深度域地震记录；m为模型参数向量，即待反演的弹性参数；F为线性反射系数矩阵。

式(12)可简记为：

s＝Gm (13)

式中，G表示深度域正演矩阵。

②基于贝叶斯推理框架构建反演目标函数为：

式中，σ_d表示数据的方差；R(m)表示模型的先验约束。

③对反演目标函数求导可得：

式中，

μ表示模型均值可用初始模型代替；当选取高斯分布为先验分布时，

C_m表示弹性参数的协方差矩阵。

④令导数等于0，此时式(15)可整理得到：

式中，μ_h表示先验信息权重系数，且

⑤对于式(16)可采用迭代重加权最小二乘算法求解最终得到准确的弹性参数，该弹性参数即可用于储层预测与流体的准确识别。

图2(a)和图2(b)分别展示了时间域稳态褶积矩阵和深度域非稳态褶积点传播函数矩阵，二者对比可以发现点传播函数随深度变化，与波的传播效应实际吻合。图3展示了本发明实施例合成的深度域地震记录，该图中子波不再对称，发生相位畸变，且随着深度增加，子波拉伸效应越明显，该现象与实际深度域地震记录相匹配。可见，图2与图3均验证了本发明该方法的理论合理性。图4展示了本发明实施例无噪声情况下采用深度域叠前线性反演得到的纵波速度、横波速度、密度，结果表明本发明方法反演的弹性参数具备较高的准确度。图5展示了本发明实施例信噪比为3:1情况下采用深度域叠前线性反演得到的纵波速度、横波速度、密度，结果表明本发明方法在噪音干扰下，仍可获取稳定、准确的结果。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (6)

1.一种基于深度域地震记录的叠前线性反演方法，其特征在于，包括：

步骤101：利用基于点传播非稳态褶积运算，对深度域测井数据进行深度域正演模拟得到深度域地震记录，将模拟的深度域地震记录与实际深度域地震记录进行平移匹配，从而获得深度域测井数据在地震尺度上的准确位置，即得到准确匹配的深度域测井数据；

步骤102：首先对整个工区的实际深度域地震记录进行层位拾取，然后将拾取的层位信息与深度构造解释结果结合，对步骤101中已经准确匹配的深度域测井数据进行插值，并将沉积模式作为协同数据参与插值，最终获取整个工区的初始低频弹性参数数据体；

步骤103：首先对整个工区的时间域地震记录进行时间域地震子波提取，然后基于点传播函数与时间域地震子波的对应关系进一步将时间转换为深度从而获取整个工区的点传播函数；

步骤104：利用获取的点传播函数构建对应的非稳态褶积矩阵形式；

步骤105：将非稳态褶积子波矩阵与反射系数线性近似公式结合，并利用贝叶斯推理线性反演方法进行整个工区的弹性参数反演。

2.根据权利要求1所述的叠前线性反演方法，其特征在于，在所述步骤101中，实际深度域地震记录即成像结果：

s＝M*d (1)

式中，M表示深度偏移算子；d表示偏移前时间域地震记录；

而偏移前时间域地震记录满足常规稳态褶积运算，其矩阵形式表达为：

d＝W*R (2)

式中，W表示常规稳态褶积运算构建的子波矩阵；R为反射系数；

由式(1)和式(2)得：

s＝(M*W)*R＝H*R (3)

式中，H表示正演运算与偏移运算结合的矩阵算子。

3.根据权利要求2所述的叠前线性反演方法，其特征在于，在所述步骤103中，一维速度模型下，点传播函数与时间域地震子波满足以下对应关系：

式中，T表示时间域地震子波的周期；

表示点传播函数在一个波长λ内的慢度；h表示点传播函数的深度；v表示速度；

在均匀速度介质中，考虑双程旅行时，式(4)改写为：

式中，f＝1/T表示频率；k＝1/λ为波速。

4.根据权利要求3所述的叠前线性反演方法，其特征在于，在所述步骤104具体为：将步骤103获取的整个工区的点传播函数利用对应的非稳态褶积矩阵形式进行排列，并在矩阵边界进行对应的点传播函数截取，以保持子波矩阵大小与常规稳态褶积一致。

5.根据权利要求4所述的叠前线性反演方法，其特征在于，非稳态褶积表示为：

式中，s(t)表示地震记录；w(t,τ)表示变化的子波；r(t)为反射系数的向量形式；

将式(6)转换为如下离散形式：

式中，s(i)、w(i,j)和r(i)分别表示地震记录、子波和反射系数的离散形式；i和j表示离散形式下的深度索引；

其中，将反射系数r(i)利用单位脉冲函数δ(i)分解为：

将式(8)代入式(7)得：

与式(2)中常规稳态褶积运算相同，此时非稳态褶积具备同样的矩阵形式：

式中，n表示子波序列的长度；m表示反射系数序列的长度；

在常规稳态褶积中：

其中，W表示稳态褶积子波矩阵，

表示非稳态褶积子波矩阵，其子波元素的上标表示不同深度下的点传播函数；此时，先通过时间域地震子波求取点传播函数，再按照稳态褶积子波矩阵对应位置截取子波元素以构建非稳态褶积子波矩阵。

6.根据权利要求5所述的叠前线性反演方法，其特征在于，在所述步骤105具体包括：

①将非稳态褶积子波矩阵与反射系数线性近似公式结合，得到深度域正演矩阵为：

式中，s表示观测数据向量，即实际深度域地震记录；m为模型参数向量，即待反演的弹性参数；F为线性反射系数矩阵；

式(12)简记为：

s＝Gm (13)

式中，G表示深度域正演矩阵；

②基于贝叶斯推理框架构建反演目标函数为：

式中，σ_d表示数据的方差；R(m)表示模型的先验约束；

③对反演目标函数求导可得：

式中，

μ表示模型均值可用初始模型代替；当选取高斯分布为先验分布时，

C_m表示弹性参数的协方差矩阵；

④令导数等于0，此时式(15)整理得到：

式中，μ_h表示先验信息权重系数，且

⑤对于式(16)采用迭代重加权最小二乘算法求解最终得到准确的弹性参数，该弹性参数即可用于储层预测与流体的准确识别。

Images