CN108663711B - 一种基于τ分布的贝叶斯地震反演方法 - Google Patents

Abstract

本发明公开了一种基于τ分布的贝叶斯地震反演方法，涉及地球物理反演及油气储层预测领域；其包括步骤：1：输入数据得到待反演参数的初始模型后确定τ分布的参数；2：选择第t道待反演参数的初始模型作为迭代初始模型，通过参数化的τ分布确定先验概率密度分布函数，并输入数据得到后验概率密度函数后进行MCMC反演得到反演结果；3：判断t是否大于输入数据中的地震道数据，若是，则结束反演得到基于反演结果的待反演波阻抗剖面，若否t累加后跳至步骤2继续反演；本发明解决了现有地震反演采用高斯分布模拟地震反射系数分布导致反演物理参数精度和分辨率低的问题，达到了提高反演物理参数的精度以及分辨率的效果。

Description

一种基于τ分布的贝叶斯地震反演方法

技术领域

本发明涉及地球物理反演及油气储层预测领域，尤其是一种基于τ分布的贝叶斯地震反演方法。

背景技术

地震反演是预测油气储层的重要步骤，它根据探测器已知的地震记录数据，通过地震记录数据与待求物理量的数学模型来建立最优化问题，并通过反演方法进行求解最优化问题，从而得到待求物理量最优估计的过程。

贝叶斯地震反演是一种地震反演技术，它通过对已知的井数据进行分析，建立先验信息来对最优化问题进行约束，并通过地震记录数据与待求物理量的关系建立似然函数，并且根据贝叶斯理论采用后验概率密度函数的形式来建立最优化问题的过程，Tarantola根据贝叶斯理论发展了后验概率密度函数评估法标志着贝叶斯反演的诞生；Lavielle将蒙特卡罗方法用于基于高斯分布的贝叶斯地震反演中，张广智将MCMC算法运用于基于高斯分布的贝叶斯反演中，证明了贝叶斯反演的可行性；理论研究及实际应用都证明了贝叶斯反演在地震勘探资料处理和反演中的实用性，能够得到精度和分辨率较好的地球物理参数。但现有技术中采用高斯分布模拟地震反射系数分布，其中高斯分布缺乏地震实际数据分析，导致反演物理参数精度不高且分辨率不高。

发明内容

本发明的目的在于：本发明提供了一种基于τ分布的贝叶斯地震反演方法，解决了现有地震反演采用高斯分布模拟地震反射系数分布导致反演物理参数精度和分辨率低的问题。

本发明采用的技术方案如下：

一种基于τ分布的贝叶斯地震反演方法，包括如下步骤：

步骤1：输入数据得到待反演参数的初始模型后确定τ分布的参数；

步骤2：选择第t道待反演参数的初始模型作为迭代初始模型，通过参数化的τ分布确定先验概率密度分布函数，并结合输入数据得到后验概率密度函数后进行MCMC反演得到反演结果；

步骤3：判断t是否大于输入数据中的地震道数据，若是，则结束反演得到基于反演结果的待反演波阻抗剖面，若否t累加后跳至步骤2继续反演。

优选地，所述步骤1包括如下步骤：

步骤1.1：输入井旁道地震数据、原始地震记录数据、子波数据和井的反射系数，对通过井旁道地震数据得到的井旁道待反演参数进行插值滤波得到待反演参数的初始模型

步骤1.2：基于井旁道待反演参数通过井旁道待反演参数和波阻抗的递推关系得到井旁道波阻抗z₀；

步骤1.3：从待反演参数的初始模型中选择井旁道对应的待反演参数的初始模型，通过递推得到井旁道波阻抗初始模型z_s，计算井旁道波阻抗z₀与井旁道波阻抗初始模型z_s的差值后并计算其均值μ(z₀-z_s)，通过公式1确定τ分布的未知参数λ，公式1如下：

λ＝α/μ(z₀-z_s)

其中，α表示已知参数，μ表示均值函数；

步骤1.4：基于步骤1.2通过井旁道波阻抗z₀计算反射系数r，计算如公式2所示：

其中，r_i表示反射系数r的第i项，z_0i表示井旁道波阻抗z₀的第i项；

基于井的反射系数r、输入数据中的子波数据w和原始地震记录数据计算噪声n，计算如公式3所示：

步骤1.5：通过噪声n求噪声方差δ_n，噪声方差δ_n计算如公式4所示：

其中，n_i表示噪声的第i项，m表示噪声中向量元素的个数，表示噪声的均值。

优选地，所述步骤2包括如下步骤：

步骤2.1：选择第t道的待反演参数的初始模型作为迭代初始模型x₁，通过反演参数与波阻抗递推关系求得第t道波阻抗z₁即公式5，公式5如下：

z₁＝f(x₁)

通过参数化的τ分布确定先验概率密度分布函数P(z₁)即公式6，公式6如下：

其中，z₁_min表示z₁-z₀中元素最小值，z₁_d表示z₁-z₀的最大值与最小值的差，z_1i表示第t道波阻抗z₁第i项，z_0i表示井旁道波阻抗z₀第i项；

步骤2.2：基于步骤1.5所得的噪声方差δ_n和输入数据中的原始地震记录和子波数据w进行贝叶斯反演得到似然函数计算如公式7所示：

其中，表示原始地震记录第i项，s_1i表示合成地震记录s₁第i项，r_1i表示地震反射系数r₁第i项，合成地震记录s₁和地震反射系数r₁计算如下：

将似然函数结合先验概率密度分布函数P(z₁)得到后验概率密度函数，计算如公式8所示：

其中，由可得；

步骤2.3：基于后验概率密度函数进行MCMC反演得到反演结果X_r。

优选地，所述步骤3包括如下步骤：

步骤3.1：判断t是否大于输入数据中的地震道数据，其中地震道数据为地震道数trace加1，若是，则结束反演得到根据公式9计算得到的待反演波阻抗剖面X_s，公式9如下：X_S(t,1)＝X_r，若否，t累加即t＝t+1后跳至步骤2继续反演。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1.本发明将传统的MCMC波阻抗反演中采用高斯分布模拟地震反射系数分布更改为采用τ分布模拟波阻抗分布，τ分布模拟波阻抗分布由实际数据得到，符合实际数据的分布，解决了现有地震反演采用高斯分布模拟地震反射系数分布导致反演物理参数精度和分辨率低的问题，达到了提高反演物理参数的精度以及分辨率的效果。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为方法流程图；

图2为含噪声的地震剖面；

图3为地震子波数据图；

图4为多道波阻抗初始模型示意图；

图5为单道井旁道波阻抗的直方图分布与对应的τ分布示意图；

图6为单道井旁道波阻抗反演结果及对比示意图；

图7为多道波阻抗反演结果示意图。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

下面结合图1-7对本发明作详细说明。

实施例1

反演参数为波阻抗时反演如下：

步骤1.1：输入井旁道地震数据、原始地震记录数据s～和子波数据w，对通过井旁道地震数据得到的井旁道待反演参数进行插值滤波得到待反演参数的初始模型

步骤1.2：基于井旁道待反演参数通过井旁道待反演参数和波阻抗的递推关系得到井旁道波阻抗z₀；

步骤1.3：从待反演参数的初始模型中选择井旁道对应的待反演参数的初始模型，通过递推得到井旁道波阻抗初始模型z_s，计算井旁道波阻抗z₀与井旁道波阻抗初始模型z_s的差值后并计算其均值μ(z₀-z_s)，α为已知参数，本实施例取4，通过公式1确定τ分布的未知参数λ，公式1如下：

λ＝α/μ(z₀-z_s)

其中，μ表示均值函数；

步骤1.4：基于步骤1.2通过井旁道波阻抗z₀计算反射系数r，计算如公式2所示：

其中，r_i表示反射系数r的第i项，z_0i表示井旁道波阻抗z₀的第i项；

基于井的反射系数r、输入数据中的子波数据w和原始地震记录数据计算噪声n，计算如公式3所示：

步骤1.5：通过噪声n求噪声方差δ_n，噪声方差δ_n计算如公式4所示：

其中，n_i表示噪声的第i项，m表示噪声中向量元素的个数，表示噪声的均值。

步骤2包括如下步骤：

步骤2.1：选择第t道的待反演参数的初始模型作为迭代初始模型通过反演参数与波阻抗递推关系如公式10求得第t道波阻抗z₁，公式10如下：

通过参数化的τ分布如公式6确定先验概率密度分布函数P(z₁)，公式6如下：

其中，z₁_min表示z₁-z₀中元素最小值，z₁_d表示z₁-z₀的最大值与最小值的差，z_1i表示第t道波阻抗z₁第i项，z_0i表示井旁道波阻抗z₀第i项；

步骤2.2：基于步骤1.5所得的噪声方差δ_n和输入数据中的原始地震记录和子波数据w进行贝叶斯反演得到似然函数计算如公式7所示：

其中，表示原始地震记录第i项，s_1i表示合成地震记录s₁第i项，r_1i表示地震反射系数r₁第i项，合成地震记录s₁和地震反射系数r₁计算如下：

将似然函数结合先验概率密度分布函数P(z₁)得到后验概率密度函数，计算如公式8所示：

步骤2.3：基于后验概率密度函数进行MCMC反演得到反演结果X_r。

步骤3包括如下步骤：

步骤3.1：判断t是否大于输入数据中的地震道数据，其中地震道数据为地震道数trace加1，若是，则结束反演得到根据公式9计算得到的待反演波阻抗剖面X_s，公式9如下：X_S(t,1)＝X_r，若否，t累加即t＝t+1后跳至步骤2继续反演。

迭代步骤如下：

S1：将原始的波阻抗模型z₁通过转移概率函数即采用均匀分布，得到新的波阻抗模型z₂，所用公式11如下：

z₂＝z₁+delta·(2·random(m,1)-1)

其中，delta表示步长选择；

S2：计算新的波阻抗的后验概率密度函数其具体步骤包括如下步骤：

S2.1：计算先验概率密度函数，所用公式12如下：

其中，z₂_min表示中元素最小值，z₂_d表示的最大值与最小值的差；

S2.2：计算似然函数，所用公式13如下：

其中，

S2.3：计算后验概率密度函数，所用公式14如下：

S3：计算判别函数a(z₁,z₂)，所用公式15如下：

S4：随机选择u＝random(0，1)，比较u与a的大小，若u<a，则新的波阻抗替换z₂原始波阻抗z₁，否则，原始波阻抗z₁不变；将每一次结果储存在矩阵Z的每一列中，直到达到预先规定的迭代次数int；取Z中后面k*delta列，求均值作为每一道最后反演结果；如图6，7所示，设定为k＝0。

效果分析：通过图4的初始模型与图7的反演结果比较，图7反演结果中有更多的细节信息，有更多的细节信息，提高了反演物理参数的精度；相同times初始模型，相对于反演结果更平滑；如图6所示initial表示初始模型中井旁道对应波阻抗，inversion表示反演的结果中井旁道对应波阻抗，true表示对应的井旁道波阻抗，单道反演结果比初始模型更接近真实数据。本发明解决了现有地震反演采用高斯分布模拟地震反射系数分布导致反演物理参数精度和分辨率低的问题，达到了提高反演物理参数的精度以及分辨率的效果。

Claims (3)

1.一种基于τ分布的贝叶斯地震反演方法，其特征在于：包括如下步骤：

步骤1：输入数据得到待反演参数的初始模型后确定τ分布的参数；

步骤2：选择第t道待反演参数的初始模型作为迭代初始模型，通过参数化的τ分布确定先验概率密度分布函数，并结合输入数据得到后验概率密度函数后进行MCMC反演得到反演结果；

步骤3：判断t是否大于输入数据中的地震道数据，若是，则结束反演得到基于反演结果的待反演波阻抗剖面，若否t累加后跳至步骤2继续反演；

所述步骤1包括如下步骤：

步骤1.1：输入井旁道地震数据、原始地震记录数据、子波数据和井的反射系数，对通过井旁道地震数据得到的井旁道待反演参数进行插值滤波得到待反演参数的初始模型

步骤1.2：基于井旁道待反演参数通过井旁道待反演参数和波阻抗的递推关系得到井旁道波阻抗z₀；

步骤1.3：从待反演参数的初始模型中选择井旁道对应的待反演参数的初始模型，通过递推得到井旁道波阻抗初始模型z_s，计算井旁道波阻抗z₀与井旁道波阻抗初始模型z_s的差值后并计算其均值μ(z₀-z_s)，通过公式1确定τ分布的未知参数λ，公式1如下：

λ＝α/μ(z₀-z_s)

其中，α表示已知参数，μ表示均值函数；

步骤1.4：基于步骤1.2通过井旁道波阻抗z₀计算反射系数r，计算如公式2所示：

其中，r_i表示反射系数r的第i项，z_0i表示井旁道波阻抗z₀的第i项；

基于井的反射系数r、输入数据中的子波数据w和原始地震记录数据计算噪声n，计算如公式3所示：

步骤1.5：通过噪声n求噪声方差δ_n，噪声方差δ_n计算如公式4所示：

其中，n_i表示噪声的第i项，m表示噪声中向量元素的个数，表示噪声的均值。

2.根据权利要求1所述的一种基于τ分布的贝叶斯地震反演方法，其特征在于：所述步骤2包括如下步骤：

步骤2.1：选择第t道的待反演参数的初始模型作为迭代初始模型x₁，通过反演参数与波阻抗递推关系求得第t道波阻抗z₁即公式5，公式5如下：

z₁＝f(x₁)

通过参数化的τ分布确定先验概率密度分布函数P(z₁)即公式6，公式6如下：

其中，z₁_min表示z₁-z₀中元素最小值，z₁_d表示z₁-z₀的最大值与最小值的差，z_1i表示第t道波阻抗z₁第i项，z_0i表示井旁道波阻抗z₀第i项；

步骤2.2：基于步骤1.5所得的噪声方差δ_n和输入数据中的原始地震记录和子波数据w进行贝叶斯反演得到似然函数计算如公式7所示：

其中，表示原始地震记录第i项，s_1i表示合成地震记录s₁第i项，r_1i表示地震反射系数r₁第i项，合成地震记录s₁和地震反射系数r₁计算如下：

将似然函数结合先验概率密度分布函数P(z₁)得到后验概率密度函数，计算如公式8所示：

其中，由可得；

步骤2.3：基于后验概率密度函数进行MCMC反演得到反演结果X_r。

3.根据权利要求1所述的一种基于τ分布的贝叶斯地震反演方法，其特征在于：所述步骤3包括如下步骤：

步骤3.1：判断t是否大于输入数据中的地震道数据，其中地震道数据为地震道数trace加1，若是，则结束反演得到根据公式9计算得到的待反演波阻抗剖面X_s，公式9如下：X_S(t,1)＝X_r，若否，t累加即t＝t+1后跳至步骤2继续反演。