CN107238862A - 基于贝叶斯反演框架的反射系数估计方法和装置 - Google Patents

基于贝叶斯反演框架的反射系数估计方法和装置 Download PDF

Info

Publication number
CN107238862A
CN107238862A CN201610190426.7A CN201610190426A CN107238862A CN 107238862 A CN107238862 A CN 107238862A CN 201610190426 A CN201610190426 A CN 201610190426A CN 107238862 A CN107238862 A CN 107238862A
Authority
CN
China
Prior art keywords
mrow
msup
reflection coefficient
coefficient sequence
sequence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201610190426.7A
Other languages
English (en)
Inventor
白英哲
徐春梅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201610190426.7A priority Critical patent/CN107238862A/zh
Publication of CN107238862A publication Critical patent/CN107238862A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6161—Seismic or acoustic, e.g. land or sea measurements

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了基于贝叶斯反演框架的反射系数估计方法和装置。反褶积是一种用来压缩地震子波、提高地震资料分辨率的手段，通常的反褶积方法往往要对地震子波和反射系数进行一系列假设从而获得反射系数的估计，而这些假设在很多情况下并不合理，不合理的假设阻碍了反褶积的应用效果。本发明中，基于贝叶斯反演框架进行稀疏盲反褶积，通过地质统计学，可以减少对子波和反射系数的不合理约束，从而获得更为精准的反射系数估计，提高地震记录的分辨率。

Description

基于贝叶斯反演框架的反射系数估计方法和装置
技术领域
本发明涉及地震资料处理领域，更具体地，涉及一种基于贝叶斯反演框架的反射系数估计方法和一种基于贝叶斯反演框架的反射系数估计装置。
背景技术
众所周知，地震记录通常是带限的，由于接收和大地非弹性参数传播效应的影响，地震资料一般都缺少低频和高频分量。低频的缺失会导致反演的多解性；而高频的缺失使得资料的分辨率比较低。反褶积效果的好坏对整个处理的影响至关重要，然而不同地区不同条件下提高分辨率的难度是不同的。大地吸收、噪声干扰、几何扩散以及地震资料处理中的一系列问题都会给分辨率造成不同的影响。提高分辨率处理的依据主要有：不同频率成分有不同的信噪比；不同信噪比的频率成分对分辨率的贡献是不同的。对于地震资料而言，提高分辨率处理的重点在于高频段信号的可靠恢复，如反Q滤波、反褶积、谱白化等方法。多道预测反褶积的主要目的是压制面波等干扰；地表一致性反褶积的主要目的是调整相位，使有效波一致化，频率适当提高。
传统的褶积模型把地震记录看作是地震子波与地下反射系数序列相褶积的结果；而反褶积则是通过压缩地震子波从地震记录中反演地下的反射系数序列。然而，地震反褶积实际上是一个“盲”过程，因为通常地表爆炸激发的地震子波是未知的，大地滤波效应(反射系数序列)也是未知的。而在地震子波和反射系数都未知时，常常要做统计性假设，如：地震子波是最小相位的或地层反射系数是高斯白噪声等，因此传统的反褶积也是一种统计性反褶积。
传统的反褶积方法通常根据这些假设条件，用地震记录的自相关代替子波的自相关，同时使用基于二阶统计量的线性滤波技术(如：预测反褶积，脉冲反褶积等)来实现子波估计与反褶积。这些假设和相应的方法在实际应用中取得了一定的效果，但是它们往往不符合实际情况，并且基于二阶统计量的反褶积方法不包含子波的相位信息。实际地下介质中传播的地震子波往往是混合相位的；地层反射系数也不是高斯白噪声。这些问题阻碍了反褶积估计地震反射系数的精度。
发明内容
本发明提出了能够获得更为准确的反射系数估计的方法和装置。针对上述反褶积中沿用的一系列不确切的假设，从而获得更为准确的反射系数估计。本发明还提出了相应的装置。
根据本发明的一方面，提出了一种基于贝叶斯反演框架的反射系数估计方法，包括：
步骤101，基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,进入步骤102；
步骤102，基于当前得到的w重新计算反射系数序列r，进入步骤103；
步骤103，基于当前得到的w和r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前得到的r即为估计的反射系数序列；如果不收敛，则将r0更新为当前得到的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
根据本发明的另一方面，提出了一种基于贝叶斯反演框架的反射系数估计装置，包括：子波序列求解单元，用于基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,反射系数求解单元，用于基于子波序列求解单元当前得到的w重新计算反射系数序列r；控制单元，用于基于子波序列求解单元当前得到的w和反射系数求解单元当前得到r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前得到的r即为估计的反射系数序列；如果不收敛，则将r0更新为当前得到的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
本发明的各方面，基于贝叶斯反演框架的盲反褶积，通过统计学原理，减少反褶积算法中的各项假设，相比于现有技术中反褶积沿用的一系列不确切的假设，能显著提高得到的反射系数的准确度。
附图说明
通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的一个实施例的基于贝叶斯反演框架的反射系数估计方法的流程图。
图2示出了中国西部某工区的地震成像剖面。
图3示出了采用传统反褶积对图2所示地震成像剖面进行处理得到的结果。
图4示出了采用本发明对图2所示地震成像剖面进行处理得到的结果。
图5示出了图2、图3、图4的频谱分析对比示意图。
具体实施方式
下面将参照附图更详细地描述本发明的优选实施方式。虽然附图中显示了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
此处，先对本发明中涉及的原理进行简单说明。
在贝叶斯理论中，后验概率分布函数PPDF以p(m|d)表示，在给定数据向量d下，模型向量m的随机分布为：
其中p(d|m)为似然函数，p(m)为先验随机分布函数。如果仅关心PPDF的形状，分母p(d)是一个可被忽略的常数。则有：
p(m|d)∝p(m)p(d|m)， 公式2
考虑地震褶积模型：
d＝Gr+n， 公式3
其中d＝[d1,d2,...,dN]T是观测到的一道地震数据，N表示地震数据的样点个数，r＝[r1,r2,...,rM]T是反射系数序列，M是反射系数序列r中的元素个数,G是N×M维子波褶积矩阵，n＝[n1,n2,...,nN]T表示观测噪声。设子波序列为w＝[w0,…,wp]T，则有N＝M+P-1。由贝叶斯公式2可得以下近似：
p(r|d)∝p(r)p(d|r)， 公式4
式中：p(r|d)表示反射系数序列r的后验概率，p(r)表示反射系数序列r的先验概率，p(d|r)表示似然函数。
假定噪声向量n服从均值为(通常设)，方差为σn 2的高斯分布即：且独立，则似然函数可等价于噪声的概率密度分布，表示为：
其中
同样不妨假设先验概率依然服从高斯分布即：
其中
把式公式5、公式6代入式公式4，则r的后验概率可表示为：
使上式取最大值的解即为反射系数序列r的最优解，也就是所谓的最大后验概率估计(Maximum a Posteriori Estimator，MAP)。
对公式7两边同时取对数：
并略去常数项-lnK1K2，定义新的目标函数为：
J＝(Gr-d)T(Gr-d)+μrTr 公式9
其中显然，对公式9的最小化即等价于对公式10的最大化。因此可定义目标函数的一般形式为：
其中，ρr(ri)是表示反射系数序列r中元素ri的先验概率。公式10即为针对非最小相位子波和非白噪反射系列情况下的盲反褶积数学框架，其中J0(w,r)表示数据误差最小二乘拟合，其目的是为了寻求最佳子波和反射系数使噪声能量最小；Jr(r)表示反射系数的先验约束项，它体现了反射系数的随机性和统计学特征。
可以采用修正的柯西准则，减小强弱对比效应，使弱反射也能较好的体现出来，其表达式为：
其中：
在有井资料、或其他地质资料可利用的情况下，通常可运用一些简单的数学变换引入先验信息，来得到地球物理反演问题更适定、更有意义的解。由于波阻抗可看作是反射系数对时间的积分，因此我们把它引入到目标函数中作为反射系数的先验约束。
由反射系数的定义可知：
其中，I(ti)表示第i层波阻抗。当反射系数足够小的情况下(|r(i)|＜0.3)，有以下近似：
上式对时间的积分，可以得到相对波阻抗的表达式：
I(t0)表示初始波阻抗值，公式14的离散形式表示为：
把公式15写成矩阵形式：
ξ＝Cr 公式16
其中C为积分算子矩阵，其形式可表示为：
则由最小二乘法可定义波阻抗约束项如下：
最终建立盲反褶积目标函数为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r) 公式19
其中第一项、第二项如公式10中定义；第三项波阻抗约束，用来控制反演结果的准确性和稳定性。μ、α分别为稀疏约束和波阻抗约束因子，μ越大，反射系数越稀疏，α越大，反演结果越准确。
根据上述约束方程，可采用松弛交替迭代法同时求取反射系数和地震子波。具体步骤如下：
1)对任意反射系数r，目标函数可写为：
将上式写成矩阵形式：
其中R表示反射系数r的褶积矩阵，Qw是由[μJr(r)+αJI(r)]T构成的对角阵。为了进一步加快收敛速度，可假设一个预条件方程：
则原矩阵方程式21可改为：
设则有：
2)、把第1)步中得到的地震子波序列w代入公式19，并求r的偏导：
当时，
即：
把上式写成矩阵形式：
其中矩阵Q的对角元素为：
最终可以得到反射系数序列r的表达式为：
r＝(GTG+μQ+αCTC)-1(GTd+αCTξ) 公式31
实施例1
图1示出了根据本发明的一个实施例的基于贝叶斯反演框架的反射系数估计方法的流程图。在本实施例中，该方法包括：
步骤101，基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,进入步骤102；
步骤102，基于求解得到的w重新计算反射系数序列r，进入步骤103；
步骤103，基于当前的w和r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前的r即为估计得到的反射系数序列；如果不收敛，则将r0更新为当前的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
本实施例基于贝叶斯反演框架的盲反褶积，通过统计学原理，减少反褶积算法中的各项假设，相比于现有技术中反褶积沿用的一系列不确切的假设，能显著提高得到的反射系数的准确度。
在一个示例中，在上述步骤101中，可以采用共轭梯度迭代求解公式24以得到，迭代终止条件为：
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε1为预设参数。可先有共轭梯度迭代计算再运用得到地震子波序列w＝[w0,…,wp]T。
其中，共轭梯度迭代的具体计算格式可以如下：①定义初始值：r←Fx-d,β←0。②迭代：{Δx←F′r，s←Δx+βs，Δr←Fs，α←-γ/(Δr·Δr)，x←x+αs，r←r+αΔr}。式中F是过渡矩阵，x表示模型空间，d表示数据空间，r表示误差，β是调节系数。
在一个示例中，在步骤102中，可以通过直接求逆的方法重新计算反射系数序列r。具体地，可以基于下式重新计算反射系数序列r：
r＝(GTG+μQ+αCTC)-1(GTd+αCTξ)， 公式31
其中，Q为对角矩阵，其第i行第i列的元素q(m_cauchy)ii为：
σr为给定值，表示在反射系数序列r的先验概率分布符合高斯分布的情况下反射系数序列r的方差。
在另一个示例中，在步骤102中，参考公式25，可以基于下式重新计算反射系数序列r：
在这种情况下，可以采用共轭梯度迭代求解公式32以得到r，迭代终止条件为：
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε2为预设参数。
实施例2
在本发明的另一实施例中，公开了一种基于贝叶斯反演框架的反射系数估计装置。在本实施例中，该装置包括子波序列求解单元、反射系数求解单元和控制单元。
子波序列求解单元，用于基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,
反射系数求解单元，用于基于子波序列求解单元当前得到的w重新计算反射系数序列r。
控制单元，用于基于子波序列求解单元当前得到的w和反射系数求解单元当前得到r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前得到的r即为估计的反射系数序列；如果不收敛，则将r0更新为当前得到的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
在一个示例中，在子波序列求解单元中，可以采用共轭梯度迭代求解公式24以得到，迭代终止条件为：
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε1为预设参数。
在一个示例中，在反射系数求解单元中，可以基于下式重新计算反射系数序列r：
r＝(GTG+μQ+αCTC)-1(GTd+αCTξ)， 公式31其中，Q为对角矩阵，其第i行第i列的元素q(m_cauchy)ii为：
σr为给定值，表示在反射系数序列r的先验概率分布符合高斯分布的情况下反射系数序列r的方差。
在另一示例中，在反射系数求解单元中，可以基于下式重新计算反射系数序列r：
在这种情况下，可以采用共轭梯度迭代求解公式32以得到r，迭代终止条件为：
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε2为预设参数。
应用示例
为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
图2示出了中国西部某工区的地震成像剖面。
图3示出了采用传统预测反褶积对图2所示地震成像剖面进行处理得到的结果。这种反褶积方法中，假设地震子波为最小相位，地震反射系数为白噪音。对比图3和图2可以看出，采用传统预测反褶积可以使同相轴变细，使得整体剖面的分辨率得以提升，但其结果依然不能令人满意。
图4示出了采用本发明对图2所示地震成像剖面进行处理得到的结果。对比图4和图1，采用本发明可以使同相轴变细，提高剖面的分辨率。在对比图4和图3，可以看出，图4中示出的同相轴更细，意味着分辨率更高，同时整体剖面更为干净，意味着噪音更小，相比于图3，图4得到的反射系数估计更为准确。
图5示出了图2、图3、图4的频谱分析对比示意图。通过对比频谱可以看出，通过这两种反褶积，均可展宽频谱，因此可提升分辨率。相比而言，传统预测反褶积的展宽效果较差，而采用本发明的展宽效果更为显著，特别在低频，端，可有效展宽了低频部分，这对于当前的地震反演成像等应用都非常重要。
本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。
计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是――但不限于――电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。
这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。
用于执行本发明操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本发明的各个方面。
这里参照根据本发明实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。
这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。
也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。
附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。

Claims (10)

1.一种基于贝叶斯反演框架的反射系数估计方法，包括：
步骤101，基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,进入步骤102；
步骤102，基于当前得到的w重新计算反射系数序列r，进入步骤103；
步骤103，基于当前得到的w和r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前得到的r即为估计的反射系数序列；如果不收敛，则将r0更新为当前得到的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
2.根据权利要求1所述的方法，其中，在步骤101中，采用共轭梯度迭代求解公式24以得到，迭代终止条件为：
<mrow> <mfrac> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>+</mo> <mo>|</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> </mfrac> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mn>1</mn> </msub> <mo>,</mo> </mrow>
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε1为预设参数。
3.根据权利要求1所述的方法，其中，在步骤102中，基于下式重新计算反射系数序列r：
r＝(GTG+μQ+αCTC)-1(GTd+αCTξ)，公式31其中，Q为对角矩阵，其第i行第i列的元素q(m_cauchy)ii为：
<mrow> <msub> <mi>q</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>_</mo> <mi>c</mi> <mi>a</mi> <mi>u</mi> <mi>c</mi> <mi>h</mi> <mi>y</mi> <mo>)</mo> <mi>i</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <msubsup> <mi>r</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mo>/</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>r</mi> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow>
σr为给定值，表示在反射系数序列r的先验概率分布符合高斯分布的情况下反射系数序列r的方差。
4.根据权利要求1所述的方法，其中，在步骤102中，基于下式重新计算反射系数序列r：
5.根据权利要求4所述的方法，其中，在步骤102中，采用共轭梯度迭代求解公式32以得到r，迭代终止条件为：
<mrow> <mfrac> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>+</mo> <mo>|</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> </mfrac> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mn>2</mn> </msub> <mo>,</mo> </mrow>
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε2为预设参数。
6.一种基于贝叶斯反演框架的反射系数估计装置，包括：
子波序列求解单元，用于基于下式求解子波序列w＝[w0,…,wp]T：
其中，R表示反射系数序列r的褶积矩阵，r＝[r1,r2,...,rM]T,r0为初始反射系数序列，M是反射系数序列r中的元素个数，Qw是由[μJr(r)+αJI(r)]T构成的对角阵，ρr(ri)表示反射系数序列r中元素ri的先验概率，C为积分算子矩阵，表示为ξ＝Cr，μ和α均为给定的约束因子，d＝[d1,d2,...,dN]T是观测到的一道地震数据,d＝Gr+n,G是N×M维子波褶积矩阵,n＝[n1,n2,...,nN]T表示观测噪声,反射系数求解单元，用于基于子波序列求解单元当前得到的w重新计算反射系数序列r；
控制单元，用于基于子波序列求解单元当前得到的w和反射系数求解单元当前得到r判断盲反褶积目标函数J(w,r)是否收敛，如果收敛则确定当前得到的r即为估计的反射系数序列；如果不收敛，则将r0更新为当前得到的r，返回步骤101，进行下一轮迭代，直至盲反褶积目标函数J(w,r)收敛，其中，盲反褶积目标函数J(w,r)为：
J(w,r)＝J0(w,r)+μJr(r)+αJI(r)，
其中，
7.根据权利要求6所述的装置，其中，在子波序列求解单元中，采用共轭梯度迭代求解公式24以得到，迭代终止条件为：
<mrow> <mfrac> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>+</mo> <mo>|</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> </mfrac> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mn>1</mn> </msub> <mo>,</mo> </mrow>
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε1为预设参数。
8.根据权利要求6所述的装置，其中，在反射系数求解单元中，基于下式重新计算反射系数序列r：
r＝(GTG+μQ+αCTC)-1(GTd+αCTξ)， 公式31
其中，Q为对角矩阵，其第i行第i列的元素q(m_cauchy)ii为：
<mrow> <msub> <mi>q</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>_</mo> <mi>c</mi> <mi>a</mi> <mi>u</mi> <mi>c</mi> <mi>h</mi> <mi>y</mi> <mo>)</mo> <mi>i</mi> <mi>i</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <msubsup> <mi>r</mi> <mi>i</mi> <mn>2</mn> </msubsup> <mo>/</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>r</mi> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> <mo>,</mo> </mrow> 2
σr为给定值，表示在反射系数序列r的先验概率分布符合高斯分布的情况下反射系数序列r的方差。
9.根据权利要求6所述的装置，其中，在反射系数求解单元中，基于下式重新计算反射系数序列r：
10.根据权利要求9所述的装置，其中，在反射系数求解单元中，采用共轭梯度迭代求解公式32以得到r，迭代终止条件为：
<mrow> <mfrac> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <msup> <mi>J</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>+</mo> <mo>|</mo> <msup> <mi>J</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>w</mi> <mo>,</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> </mfrac> <mo>&amp;le;</mo> <msub> <mi>&amp;epsiv;</mi> <mn>2</mn> </msub> <mo>,</mo> </mrow>
其中，Jk(w,r)表示该共轭梯度迭代中第k次迭代后得到的J(w,r)，ε2为预设参数。
CN201610190426.7A 2016-03-29 2016-03-29 基于贝叶斯反演框架的反射系数估计方法和装置 Pending CN107238862A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610190426.7A CN107238862A (zh) 2016-03-29 2016-03-29 基于贝叶斯反演框架的反射系数估计方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610190426.7A CN107238862A (zh) 2016-03-29 2016-03-29 基于贝叶斯反演框架的反射系数估计方法和装置

Publications (1)

Publication Number Publication Date
CN107238862A true CN107238862A (zh) 2017-10-10

Family

ID=59983946

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610190426.7A Pending CN107238862A (zh) 2016-03-29 2016-03-29 基于贝叶斯反演框架的反射系数估计方法和装置

Country Status (1)

Country Link
CN (1) CN107238862A (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110308483A (zh) * 2019-05-23 2019-10-08 中国石油天然气股份有限公司 基于多任务贝叶斯压缩感知的反射系数求取方法及装置
CN110333534A (zh) * 2019-06-19 2019-10-15 中国地质大学（武汉） 一种基于Biot理论的贝叶斯时移AVO反演方法及系统
CN111505713A (zh) * 2020-01-21 2020-08-07 长江大学 基于多点地质统计的叠前地震反演方法
CN111694052A (zh) * 2019-03-14 2020-09-22 中国石油天然气股份有限公司 盲反演方法及装置
CN112578439A (zh) * 2019-09-29 2021-03-30 中国石油化工股份有限公司 一种基于空间约束的地震反演方法
CN113341463A (zh) * 2021-06-10 2021-09-03 中国石油大学(北京) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN113534239A (zh) * 2020-04-18 2021-10-22 中国石油化工股份有限公司 一种基于信号分离的波形压缩方法
CN113552631A (zh) * 2021-08-16 2021-10-26 中煤科工集团西安研究院有限公司 一种用于窄带信号的时频双域正则化稀疏反褶积方法及装置

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853930A (zh) * 2014-03-19 2014-06-11 中国科学院地质与地球物理研究所 地震矢量波场数值模拟方法和装置
US9075159B2 (en) * 2011-06-08 2015-07-07 Chevron U.S.A., Inc. System and method for seismic data inversion

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9075159B2 (en) * 2011-06-08 2015-07-07 Chevron U.S.A., Inc. System and method for seismic data inversion
CN103853930A (zh) * 2014-03-19 2014-06-11 中国科学院地质与地球物理研究所 地震矢量波场数值模拟方法和装置

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SAM ZANDONG SUN ET AL.: "PSO non-linear pre-stack inversion method and the application in reservoir prediction", 《SEG LAS VEGAS 2012 ANNUAL MEETING》 *
刘杰: "基于独立分量分析的地震盲反褶积方法及应用研究", 《中国优秀硕士学位论文全文数据库 基础科学辑》 *
张繁昌等: "修正柯西约束地震盲反褶积方法", 《石油地球物理勘探》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111694052A (zh) * 2019-03-14 2020-09-22 中国石油天然气股份有限公司 盲反演方法及装置
CN110308483A (zh) * 2019-05-23 2019-10-08 中国石油天然气股份有限公司 基于多任务贝叶斯压缩感知的反射系数求取方法及装置
CN110333534A (zh) * 2019-06-19 2019-10-15 中国地质大学（武汉） 一种基于Biot理论的贝叶斯时移AVO反演方法及系统
CN110333534B (zh) * 2019-06-19 2020-10-30 中国地质大学（武汉） 一种基于Biot理论的贝叶斯时移AVO反演方法及系统
CN112578439A (zh) * 2019-09-29 2021-03-30 中国石油化工股份有限公司 一种基于空间约束的地震反演方法
CN111505713A (zh) * 2020-01-21 2020-08-07 长江大学 基于多点地质统计的叠前地震反演方法
CN111505713B (zh) * 2020-01-21 2021-05-07 长江大学 基于多点地质统计的叠前地震反演方法
CN113534239A (zh) * 2020-04-18 2021-10-22 中国石油化工股份有限公司 一种基于信号分离的波形压缩方法
CN113341463A (zh) * 2021-06-10 2021-09-03 中国石油大学(北京) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN113341463B (zh) * 2021-06-10 2023-05-26 中国石油大学(北京) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN113552631A (zh) * 2021-08-16 2021-10-26 中煤科工集团西安研究院有限公司 一种用于窄带信号的时频双域正则化稀疏反褶积方法及装置
CN113552631B (zh) * 2021-08-16 2023-11-03 中煤科工集团西安研究院有限公司 一种用于窄带信号的时频双域正则化稀疏反褶积方法及装置

Similar Documents

Publication Publication Date Title
CN107238862A (zh) 基于贝叶斯反演框架的反射系数估计方法和装置
CN113127931B (zh) 基于瑞丽散度进行噪声添加的联邦学习差分隐私保护方法
Yue et al. Bayesian inference for additive mixed quantile regression models
CN110869814B (zh) 用于地震数据的全波形反演的系统和方法
Smyl et al. Data preprocessing and augmentation for multiple short time series forecasting with recurrent neural networks
CN109859120B (zh) 基于多尺度残差网络的图像去雾方法
US20090238406A1 (en) Dynamic state estimation
US11978146B2 (en) Apparatus and method for reconstructing three-dimensional image
CN109188410B (zh) 一种非视距场景下的距离校准方法、装置及设备
CN105427264A (zh) 一种基于群稀疏系数估计的图像重构方法
CN116451553A (zh) 一种改进的变分模态分解与BiGRU融合的水质预测方法
Dwivedi Quantifying predictability of Indian summer monsoon intraseasonal oscillations using nonlinear time series analysis
Sun et al. Full-waveform inversion using a learned regularization
Misra et al. Effect of using genetic algorithm to denoise MRI images corrupted with Rician Noise
CN107622476A (zh) 基于概率生成模型的图像超分辨处理方法
Wolter Frequency domain methods in recurrent neural networks for sequential data processing
Kåresen Evaluation and applications of the iterated window maximization method for sparse deconvolution
CN113537354B (zh) 一种基于深度学习的含水层结构阶段式随机反演识别方法
DE112015007217T5 (de) HMM-basiertes adaptives Spektrogrammnachverfolgungsverfahren
Drissi et al. Second-order optimization strategies for neural network quantum states
CN108805286A (zh) 基于高阶残留量化的高性能网络加速方法
Hanebeck et al. Progressive Bayesian estimation with deterministic particles
Yu et al. Modeling spatio-temporal extreme events using graphical models
Lee et al. Total variation-based image noise reduction with generalized fidelity function
Lyu Calibrating an earth system model using the adjoint method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20171010