一种基于空间约束的地震反演方法与流程

本发明属于地球物理反演领域，具体涉及一种基于空间约束的地震反演方法。

背景技术：

随着勘探目标转向复杂型油气藏，非常规油气勘探开发已成为研究重点。通常通过反演算法处理后获得的地震数据预测油气储藏情况。目前大多数基于稀疏反演的算法在存在噪声的情况下不太稳定，不容易恢复原本构造，横向连续性较差的。对此特别需要一种提高横向连续性且容易恢复原本构造的反演方法。

技术实现要素：

本发明的目的是提出一种提高横向连续性且容易恢复原本构造的基于空间约束的地震反演方法。

为了实现上述目的，本发明提供一种基于空间约束的地震反演方法，包括：基于褶积模型和井旁地震道，提取地震子波；构建反演目标函数；基于所述反演目标函数，构建最终反演目标函数；对所述最终反演目标函数进行求解，获得反射系数。

优选的，所述基于褶积模型和井旁地震道，提取地震子波包括：基于所述褶积模型，利用测井资料中的速度信息和密度信息计算反射系数序列，结合所述反射系数序列和井旁地震道提取所述地震子波。

优选的，所述构建反演目标函数包括：基于贝叶斯定理，结合似然函数与先验分布，获得后验概率分布；使所述后验概率分布极大，获得所述反演目标函数。

优选的，所述后验概率分布为：

其中，p(m|d)是后验概率分布，d为观测数据，g为地震子波褶积矩阵，m为反射系数，cd为协方差，μ为正则化参数，f(m)是由先验分布得到的正则化项，其为柯西分布，p为待反演数据的长度，σ为尺度因子。

优选的，所述反演目标函数为：

其中，j(m)表示反演目标函数。

优选的，在所述反演目标函数中引入空间平滑约束，获得最终反演目标函数。

优选的，所述最终反演目标函数为：

其中，lfx为fx滤波算子。

优选的，所述对所述最终反演目标函数进行求解，获得反射系数包括：步骤401：设置对角加权矩阵；步骤402：基于所述对角加权矩阵求取反射系数；步骤403：对所述反射系数进行fx空间滤波，计算最终反演目标函数值；步骤404：更新所述对角加权矩阵，重复执行步骤402-步骤403，直到所述最终反演目标函数值满足收敛条件，所对应的反射系数即为最终反演的反射系数。

优选的，在所述步骤402中，通过以下公式求取反射系数：

m^j＝(g^tg+μq^j-1)^-1g^td^j，

其中，m为反射系数，d为观测数据，g为地震子波褶积矩阵，j为迭代次数，μ为正则化参数，q为对角加权矩阵。

优选的，通过以下公式更新对角加权矩阵：

qkk＝1/(1+(lfxm)k²/σ²)

其中，q为对角加权矩阵，k为数据列的第k个元素，lfx为fx滤波算子，m为反射系数，σ为尺度因子。

本发明的有益效果在于：本发明的基于空间约束的地震反演方法通过提取地震子波、建立最终反演目标函数，以及求取最终反演目标函数，与常规反演相比，在利用单道信息的基础上结合多道信息，使反演结果更真实可靠。

本发明在考虑实际资料存在噪音，常规基于稀疏反演的算法反演结果不太稳定，对此假设反射系数在空间上是连续的，在稀疏反演的基础上引入空间平滑约束，利用fx滤波来增强反演结果横向连续性，稳定反演结果，恢复弱小构造。

本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。

附图说明

通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。

图1示出了根据本发明的一个实施例的基于空间约束的地震反演方法的流程图。

图2示出了根据本发明的一个实施例的基于空间约束的地震反演方法的真实波阻抗楔形模型图。

图3示出了根据本发明的一个实施例的基于空间约束的地震反演方法的模型加入噪音后的地震数据。

图4示出了根据本发明的一个实施例的基于空间约束的地震反演方法的仅加入稀疏约束的反射系数。

图5示出了根据本发明的一个实施例的基于空间约束的地震反演方法的反射系数。

图6示出了根据本发明的一个实施例的基于空间约束的地震反演方法的地震数据反演结果。

具体实施方式

下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。

根据本发明的一种基于空间约束的地震反演方法，包括：基于褶积模型和井旁地震道，提取地震子波；构建反演目标函数；基于反演目标函数，构建最终反演目标函数；对最终反演目标函数进行求解，获得反射系数。

具体地，基于褶积模型和井旁地震道提取地震子波；构建反演目标函数后，在反演目标函数的基础上再构建最终反演目标函数，对最终反演目标函数进行求解，获得反射系数。

根据示例性的实施方式，基于空间约束的地震反演方法基于空间约束的地震反演方法通过提取地震子波、建立最终反演目标函数，以及求取最终反演目标函数，与常规反演相比，在利用单道信息的基础上结合多道信息，使反演结果更真实可靠。

作为优选方案，基于褶积模型和井旁地震道，提取地震子波包括：基于褶积模型，利用测井资料中的速度信息和密度信息计算反射系数序列，结合反射系数序列和井旁地震道提取地震子波。

具体的，将提取的地震子波记为地震子波w。

作为优选方案，构建反演目标函数包括：基于贝叶斯定理，结合似然函数与先验分布，获得后验概率分布；使后验概率分布极大，获得反演目标函数。

作为优选方案，后验概率分布为：

其中，p(m|d)是后验概率分布，d为观测数据，g是子波褶积矩阵，m为反射系数，cd为协方差，μ为正则化参数，f(m)是由先验分布得到的正则化项，其为柯西分布，p为待反演数据的长度，σ为尺度因子。

作为优选方案，反演目标函数为：

其中，j(m)表示反演目标函数。

具体的，已知贝叶斯定理可表示为p(m|d)∝p(d|m)p(m)

p(m|d)是后验概率分布，p(d|m)是观测数据的似然函数，p(m)是模型参数的先验概率分布。

另外地震记录褶积模型的向量形式有

d＝gm+n

其中，d是观测数据，g为地震子波褶积矩阵，表示为：

w为地震子波，k为维数，实际地震数据反演过程中，子波矩阵g的维数k，需要根据实际地震数据采样点的个数进行确定。n表示数据噪音和一些理论误差，m为反射系数序列。假设噪音服从协方差为cd的高斯分布，则似然函数可表示为：

其中，n是待反演数据的长度。

假设先验分布给定，则后验分布有如下形式：

其中，u为正则化参数，f(m)是由先验分布得到的正则化项。

为使后验概率分布函数极大，则目标函数应有最小值：

其中，f(m)是关于所选择的先验分布的正则化项的形式。本发明采用的是柯西分布，特征在于有较窄的峰值宽度，向两侧的过渡也较平缓，具有“长尾巴”分布特征，能够平衡分辨率与弱反射信息，一定程度上压制噪音，其形式为σ为尺度因子。

作为优选方案，在反演目标函数中引入空间平滑约束，获得最终反演目标函数。

作为优选方案，最终反演目标函数为：

其中，lfx为fx滤波算子。

具体的，柯西约束虽然能保证反射系数在时间域的稀疏性，除此之外，本发明还引入空间平滑约束，利用地震道之间的信息来增加地震数据空间域的连续性，具体目标函数形式为：

其中，lfx为fx滤波算子。

本发明在考虑实际资料存在噪音，常规基于稀疏反演的算法反演结果不太稳定，对此假设反射系数在空间上是连续的，在稀疏反演的基础上引入空间平滑约束，利用fx滤波来增强反演结果横向连续性，稳定反演结果，恢复弱小构造。

作为优选方案，对最终反演目标函数进行求解，获得反射系数包括：步骤401：设置对角加权矩阵；步骤402：基于对角加权矩阵求取反射系数；步骤403：对反射系数进行fx空间滤波，计算最终反演目标函数值；步骤404：更新对角加权矩阵，重复执行步骤402-步骤403，直到最终反演目标函数值满足收敛条件，所对应的反射系数即为最终反演的反射系数。

作为优选方案，在步骤402中，通过以下公式求取反射系数：

m^j＝(g^tg+μq^j-1)^-1g^td^j，

其中，m为反射系数，d为观测数据，g为地震子波褶积矩阵，j为迭代次数，μ为正则化参数，q为对角加权矩阵。

作为优选方案，通过以下公式更新对角加权矩阵：

qkk＝1/(1+(lfxm)k²/σ²)

其中，q为对角加权矩阵，k为数据列的第k个元素，lfx为fx滤波算子，m为反射系数，σ为尺度因子。

具体的，(1)初始化对角加权矩阵q，设置对角加权矩阵q为单位矩阵，另外设置μ、σ和滤波器参数，例如，μ＝0.5，σ＝0.01，滤波器范围为10hz到100hz，权重为0.01，利用公式m^j＝(g^tg+μq^j-1)^-1g^td^j求取第一次迭代的初始反射系数m，j为迭代次数；

(2)对反射系数进行fx空间滤波，将m代入公式μf(lfxm)实现对反射系数进行fx空间滤波，进而计算获得目标函数值j(m)；

(3)利用公式qkk＝1/(1+(lfxm)k²/σ²)更新对角加权矩阵，k表示数据列的第k个元素，根据公式m^j＝(g^tg+μq^j-1)^-1g^td^j重新迭代获得反射系数，进而计算目标函数值；直至迭代满足收敛条件，此时的最终值即为反射系数反演结果。一般重复迭代4-5次即可得到相对满意的结果。

实施例

图1示出了根据本发明的一个实施例的基于空间约束的地震反演方法的流程图。图2示出了根据本发明的一个实施例的基于空间约束的地震反演方法的真实波阻抗楔形模型图。图3示出了根据本发明的一个实施例的基于空间约束的地震反演方法的模型加入噪音后的地震数据。图4示出了根据本发明的一个实施例的基于空间约束的地震反演方法的仅加入稀疏约束的反射系数。图5示出了根据本发明的一个实施例的基于空间约束的地震反演方法的反射系数。图6示出了根据本发明的一个实施例的基于空间约束的地震反演方法的地震数据反演结果。

如图1所示，一种基于空间约束的地震反演方法，包括：

s102：基于褶积模型和井旁地震道，提取地震子波；

其中，基于褶积模型和井旁地震道，提取地震子波包括：基于褶积模型，利用测井资料中的速度信息和密度信息计算反射系数序列，结合反射系数序列和井旁地震道提取地震子波；

具体的，将提取的地震子波记为地震子波w。

例如，褶积模型采用楔形模型，楔形模型模拟三层地层结构，从上至下地震速度依次为2000m/s、2500m/s、3000m/s,模型采用100个采样点，60道地震数据，地层密度默认都为2.1kg/m³。例如，直接采用主频为20hz的雷克子波。如图2所示为利用楔形模型速度和密度计算得到的波阻抗数据，图3为在图2的真实波阻抗楔形模型中加入信噪比为4的噪音后得到的地震数据。

s104：构建反演目标函数；

其中，构建反演目标函数包括：基于贝叶斯定理，结合似然函数与先验分布，获得后验概率分布；使后验概率分布极大，获得反演目标函数；

其中，后验概率分布为：

其中，p(m|d)是后验概率分布，d为观测数据，g为地震子波褶积矩阵，

m为反射系数，cd为协方差，μ为正则化参数，f(m)是由先验分布得到的正则化项，其为柯西分布，p为待反演数据的长度，σ为尺度因子；

其中，反演目标函数为：

其中，j(m)表示反演目标函数；

具体实现方式为：已知贝叶斯定理可表示为p(m|d)∝p(d|m)p(m)

p(m|d)是后验概率分布，p(d|m)是观测数据的似然函数，p(m)是模型参数的先验概率分布。

另外地震记录褶积模型的向量形式有

d＝gm+n

其中，d是观测数据，g是子波褶积矩阵，表示为：

w为地震子波，k为维数，实际地震数据反演过程中，子波矩阵g的维数k，需要根据实际地震数据采样点的个数进行确定。n表示数据噪音和一些理论误差，m为反射系数序列。假设噪音服从协方差为cd的高斯分布，则似然函数可表示为：

其中，n是待反演数据的长度。

假设先验分布给定，则后验分布有如下形式：

其中，u为正则化参数，f(m)是由先验分布得到的正则化项。

为使后验概率分布函数极大，则目标函数应有最小值：

其中，f(m)是关于所选择的先验分布的正则化项的形式。本发明采用的是柯西分布，特征在于有较窄的峰值宽度，向两侧的过渡也较平缓，具有“长尾巴”分布特征，能够平衡分辨率与弱反射信息，一定程度上压制噪音，其形式为σ为尺度因子。

如图4所示，为采用柯西分布的正则化项得到的反射系数反演结果，与原始模型对比，可以看到地层界面模糊，基本识别不清，除此之外还存在很多条带状现象，严重影响地下构造识别。

s106：基于反演目标函数，构建最终反演目标函数；

其中，在反演目标函数中引入空间平滑约束，获得最终反演目标函数；

其中，最终反演目标函数为：

其中，lfx为fx滤波算子；

柯西约束虽然能保证反射系数在时间域的稀疏性，除此之外，本发明还引入空间平滑约束，利用地震道之间的信息来增加地震数据空间域的连续性，具体目标函数形式为其中lfx为fx滤波算子。

s108：对最终反演目标函数进行求解，获得反射系数；

其中，对最终反演目标函数进行求解，获得反射系数包括：

步骤401：设置对角加权矩阵；

步骤402：基于对角加权矩阵求取反射系数；

其中，在步骤402中，通过以下公式求取反射系数：

m^j＝(g^tg+μq^j-1)^-1g^td^j，

其中，m为反射系数，d为观测数据，g为地震子波褶积矩阵，j为迭代次数，μ为正则化参数，q为对角加权矩阵；

步骤403：对反射系数进行fx空间滤波，计算最终反演目标函数值；

步骤404：更新对角加权矩阵，重复执行步骤402-步骤403，直到最终反演目标函数值满足收敛条件，所对应的反射系数即为最终反演的反射系数。

其中，通过以下公式更新对角加权矩阵：

qkk＝1/(1+(lfxm)k²/σ²)

其中，q为对角加权矩阵，k为数据列的第k个元素，lfx为fx滤波算子，m为反射系数，σ为尺度因子。

(1)初始化对角加权矩阵q，设置对角加权矩阵q为单位矩阵，另外设置μ、σ和滤波器参数，例如，μ＝0.5，σ＝0.01，滤波器范围为10hz到100hz，权重为0.01，利用公式m^j＝(g^tg+μq^j-1)^-1g^td^j求取第一次迭代的初始反射系数m，j为迭代次数；

(2)对反射系数进行fx空间滤波，将m代入公式μf(lfxm)实现对反射系数进行fx空间滤波，进而计算获得目标函数值j(m)；

(3)利用公式qkk＝1/(1+(lfxm)k²/σ²)更新对角加权矩阵，k表示数据列的第k个元素，根据公式m^j＝(g^tg+μq^j-1)^-1g^td^j重新迭代获得反射系数，进而计算目标函数值；直至迭代满足收敛条件，此时的最终值即为反射系数反演结果。

如图5所示为采用本发明得到的反射系数反演结果，与只利用柯西约束的反演结果相比，可以看到反演结果基本能够反映地下地层结构，反射界面清晰，虽然仍存在一些条带状，但基本不影响界面识别。

如图6所示为利用图5反演结果计算得到的地震数据，与加入噪音的模型地震数据相比，可以看到地震数据得到了很好恢复，噪音去除较明显，最终也证明了该方法的有效性。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。