基于贝叶斯反演框架的反射系数估计方法和装置与流程
本发明涉及地震资料处理领域，更具体地，涉及一种基于贝叶斯反演框架的反射系数估计方法和一种基于贝叶斯反演框架的反射系数估计装置。
背景技术：
众所周知，地震记录通常是带限的，由于接收和大地非弹性参数传播效应的影响，地震资料一般都缺少低频和高频分量。低频的缺失会导致反演的多解性；而高频的缺失使得资料的分辨率比较低。反褶积效果的好坏对整个处理的影响至关重要，然而不同地区不同条件下提高分辨率的难度是不同的。大地吸收、噪声干扰、几何扩散以及地震资料处理中的一系列问题都会给分辨率造成不同的影响。提高分辨率处理的依据主要有：不同频率成分有不同的信噪比；不同信噪比的频率成分对分辨率的贡献是不同的。对于地震资料而言，提高分辨率处理的重点在于高频段信号的可靠恢复，如反q滤波、反褶积、谱白化等方法。多道预测反褶积的主要目的是压制面波等干扰；地表一致性反褶积的主要目的是调整相位，使有效波一致化，频率适当提高。
传统的褶积模型把地震记录看作是地震子波与地下反射系数序列相褶积的结果；而反褶积则是通过压缩地震子波从地震记录中反演地下的反射系数序列。然而，地震反褶积实际上是一个“盲”过程，因为通常地表爆炸激发的地震子波是未知的，大地滤波效应(反射系数序列)也是未知的。而在地震子波和反射系数都未知时，常常要做统计性假设，如：地震子波是最小相位的或地层反射系数是高斯白噪声等，因此传统的反褶积也是一种统计性反褶积。
传统的反褶积方法通常根据这些假设条件，用地震记录的自相关代替子波的自相关，同时使用基于二阶统计量的线性滤波技术(如：预测反褶积，脉冲反褶积等)来实现子波估计与反褶积。这些假设和相应的方法在实际应用中取得了一定的效果，但是它们往往不符合实际情况，并且基于二阶统计量的反褶积方法不包含子波的相位信息。实际地下介质中传播的地震子波往往是混合相位的；地层反射系数也不是高斯白噪声。这些问题阻碍了反褶积估计地震反射系数的精度。
技术实现要素：
本发明提出了能够获得更为准确的反射系数估计的方法和装置。针对上述反褶积中沿用的一系列不确切的假设，从而获得更为准确的反射系数估计。本发明还提出了相应的装置。
根据本发明的一方面，提出了一种基于贝叶斯反演框架的反射系数估计方法，包括：
步骤101，基于下式求解子波序列w＝[w0,…,wp]t：
其中，
步骤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：
其中，
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服从均值为
其中
同样不妨假设先验概率依然服从高斯分布即：
其中
把式公式5、公式6代入式公式4，则r的后验概率可表示为：
使上式取最大值的解即为反射系数序列r的最优解，也就是所谓的最大后验概率估计(maximumaposterioriestimator，map)。
对公式7两边同时取对数：
并略去常数项-lnk1k2，定义新的目标函数为：
j＝(gr-d)t(gr-d)+μrtr公式9
其中
其中，
可以采用修正的柯西准则，减小强弱对比效应，使弱反射也能较好的体现出来，其表达式为：
其中：
在有井资料、或其他地质资料可利用的情况下，通常可运用一些简单的数学变换引入先验信息，来得到地球物理反演问题更适定、更有意义的解。由于波阻抗可看作是反射系数对时间的积分，因此我们把它引入到目标函数中作为反射系数的先验约束。
由反射系数的定义可知：
其中，i(ti)表示第i层波阻抗。当反射系数足够小的情况下(|r(i)|＜0.3)，有以下近似：
上式对时间的积分，可以得到相对波阻抗的表达式：
i(t0)表示初始波阻抗值，公式14的离散形式表示为：
把公式15写成矩阵形式：
ξ＝cr公式16
其中
则由最小二乘法可定义波阻抗约束项如下：
最终建立盲反褶积目标函数为：
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：
其中，
步骤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为预设参数。可先有共轭梯度迭代计算
其中，共轭梯度迭代的具体计算格式可以如下：①定义初始值：r←fx-d,β←0。②迭代：{δx←f′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：
其中，
反射系数求解单元，用于基于子波序列求解单元当前得到的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)，该电子电路可以执行计算机可读程序指令，从而实现本发明的各个方面。
这里参照根据本发明实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。
这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。
也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。
附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
- 还没有人留言评论。精彩留言会获得点赞！