CN101354442B - 一种用于获取地层信息的混合相位反褶积方法及处理*** - Google Patents

Abstract

本发明提供一种用于获取地层信息的混合相位反褶积方法及处理***，为了适应复杂地区的地震资料中地震子波都是混合相位这一特点，使用本发明提供的混合相位反褶积的方法处理反射地震记录，首先求取与实际地震子波最匹配的混合相位子波，然后再求得混合相位子波的反子波滤波器，用这个反子波滤波器与实际地震道褶积，得到高分辨率的地震记录信息以提高地震勘探的精度。

Description

一种用于获取地层信息的混合相位反褶积方法及处理***

技术领域

本发明涉及石油地质勘探技术领域中的地震资料处理技术，尤其涉及一种用于获取地层信息的混合相位反褶积方法及处理***，以提高反射地震记录的分辨能力。

背景技术

在油气勘探过程中，了解地层信息非常重要，通过地层信息，可以得知一地区地下有没有油、多少米下地层有油、油有多厚等，以为是否开展地质勘探工作提供有效的指导帮助，因此，地层信息的获得是地质勘探的首要工作之一，而获得准确的地层信息一直是业界非常关注和亟待解决的。

目前，为了获得地层信息，常用的方法是采用反褶积，即反滤波方法对反射地震资料进行处理和分析，由于反射地震资料中不仅包含了对于地震勘探有用的地层信息，还包含了对于地质勘探无用的地震子波相关信息，因此，需要采用反褶积方法将地震子波相关信息去除掉，进而得到有用的地层信息。

然而，常规反褶积常要做一系列假设，其中最常用的就是假设地震子波是最小相位的，基于这种假设条件的方法我们称之为最小相位反褶积，也称为常规反褶积，而实际上这种假设条件既非严格理论推导出的结论，也非实测结果，只不过是在解决实际问题时使这一问题可以简单求解而做的一种近似假设，在地表与地下构造较简单的地区，如海洋、东部油田等地的地震勘探，这种假设尚可接受，然而在复杂区，由于地震勘探对地震子波相位影响很大，所以一般情况下实际地震子波是混合相位的，如此一来，上述假设地震子波是最小相位来进行反褶积就不再适用。

发明内容

为了弥补常规反褶积方法，也即最小相位反褶积方法对于复杂地区地震资料处理的不适用性，本发明从针对复杂区地震勘探的目标出发，对反褶积这一地震资料处理的基本方法进行了研究，期望用更符合复杂区实际地震资料的混合相位反褶积，来取代现在使用的基于地震子波最小相位假设的反褶积处理技术。

为了达成上述目的，本发明的技术方案是这样实现的：

一种用于获取地层信息的混合相位反褶积方法，该方法包括下列步骤：采集反射地震记录，选择该反射地震记录中的一条标准道的反射地震记录，该反射地震记录包括实际地震子波的频谱；应用谱模拟的方法，在空间窗口和时间窗口中求取与所述实际地震子波频谱相位对应的最小相位子波；采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波，该混合相位子波是时空变化的；采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波；求得所述选定的混合相位子波的反滤波器；采用多层可变时窗，将上述混合相位子波反滤波器与所述反射地震记录褶积，得到所述标准道的地层信息。

一种混合相位反褶积处理***，用于通过分析反射地震记录，获取反射地震记录中的地层信息，该***包括：采集单元，用于采集反射地震记录，选择该反射地震记录中的一条标准道的反射地震记录，该反射地震记录包括实际地震子波的频谱；最小相位子波计算单元，用于应用谱模拟的方法，在空间窗口和时间窗口中求取与所述实际地震子波频谱相位对应的最小相位子波；混合相位子波计算单元，用于采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波，该混合相位子波是时空变化的；混合相位子波选择单元，用于采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波；混合相位子波反滤波器计算单元，用于根据所述混合相位子波计算混合相位子波反滤波器；褶积单元，用于采用多层可变时窗，将上述混合相位子波反滤波器与所述反射地震记录褶积，得到所述标准道的地层信息。

本发明的有益效果包括：

本发明提供的一种用于获取地层信息的混合相位反褶积方法及处理***，应用混合相位反褶积方法处理反射地震记录，压缩了地震子波的长度，获得的地层信息比常规方法获得的地层信息更为精确，提高了反射地震记录的分辨能力。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：

图1为本发明的方法流程图；

图2A为本发明第一实施例的发射系数序列示意图；

图2B为本发明第一实施例的发射系数序列与原始混合相位子波的褶积地震道模型；

图3A为本发明第一实施例的地震子波振幅谱示意图；

图3B为本发明第一实施例的相位子波示意图；

图4A为本发明第一实施例的最小平方反褶积结果示意图；

图4B为本发明第一实施例的混合相位反褶积结果示意图；

图5为本发明第二实施例的模型反射系数序列示意图；

图6为本发明第二实施例的子波时变地震道模型示意图；

图7A为本发明第二实施例的单时窗提取整道地震子波频谱(中间频谱曲线)与给定的两个子波频谱对比示意图；

图7B为图6中地震道的实际频谱曲线示意图；

图8为分时窗提取的子波频谱示意图；

图9A为单时窗混合相位反褶积剖面示意图；

图9B为时变混合相位反褶积剖面示意图；

图10为本发明第三实施例的振幅谱示意图；

图11A为本发明第三实施例的最小平方反褶积得到的子波示意图；

图11B为本发明第三实施例的混合相位反褶积法得到的子波示意图；

图12A为本发明第三实施例的原始实际地震炮集示意图；

图12B为本发明第三实施例的最小平方反褶积结果示意图；

图12C为本发明第三实施例的混合相位反褶积结果示意图；

图13为本发明第四实施例的振幅谱示意图；

图14A为本发明第四实施例的混合相位反褶积法得到的子波示意图；

图14B为本发明第四实施例的最小平方反褶积法得到的子波示意图；

图15A为本发明第四实施例的原始实际地震剖面示意图；

图15B为本发明第四实施例的最小平方反褶积结果示意图；

图15C为本发明第四实施例的混合相位反褶积结果示意图；

图16为本发明第五实施例的海洋某地区实际叠后地震剖面示意图；

图17为本发明第五实施例的单时窗反褶积示意图；

图18为本发明第五实施例的另一单时窗反褶积示意图；

图19为本发明第五实施例的分时窗反褶积示意图；

图20A-图20D为本发明第五实施例的地震剖面混合相位反褶积处理的局部放大图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合实施例和附图，对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。

本发明的核心内容是，为了适应复杂地区的地震资料中地震子波都是混合相位这一特点，使用本发明提供的混合相位反褶积的方法处理反射地震记录，首先求取与实际地震子波最匹配的混合相位子波，然后再求得混合相位子波的反子波滤波器，用这个反子波滤波器与实际地震道褶积，得到高分辨率的地震记录信息以提高地震勘探的精度。

为了使本发明更加浅显易懂，同时也为了明确下面描述中的术语，在此，先对本领域的几个概念以及本发明涉及的术语进行简单说明：

相位子波是地震子波用最大相位、最小相位、零相位、混合相位区分的，而相位子波反滤波器可以看做是用来改变地震子波相位的一种“算法”。

混合相位反褶积的原理：

地震道在许多情形下，地震子波是未知的，而反射地震记录的自相关函数与地震子波的自相关函数相等，因此假定反射系数序列是白化的，并且与噪声不相关，则可从地震道中获取地震子波脉冲的自相关函数的估计，接着，通过这个自相关函数，计算出一个最小相位的反滤波器，如果子波是最小相位的，则这个反滤波器是最佳的，通过计算反滤波器的自相关函数及求解相应的方程组，来获得一个估计的最小相位地震子波。

混合相位反滤波器等于一个最小相位反滤波器和一个全通滤波器的褶积的原理及推导步骤：

将一个混合相位子波分解为最小相位和最大相位两个分量：

假设混合相位子波脉冲为p_t＝{p₀，.....，p_N}，那么它的Z变换可以表示为：

$P\left(Z\right)=\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{p}_j{Z}^j---(1-1)$

假定P(Z)在单位圆上没有零点，在单位圆内有β个零点，而在单位圆外有N-β个零点，于是，P(Z)可被分解为一个最小相位分量A(Z)和一个最大相位分量Z^-βB(Z^-1)，如下所示：

$P\left(Z\right)={\hat{p}}_{0}A\left(Z\right)Z^{\mathrm{\β}}B\left(Z^{-1}\right)---(1-2)$

其中，

是由下式定义的最小相位子波延迟脉冲

的第一个样点。因此可以得出：

A(Z)＝1+a₁Z+......+a_N-βZ^N-β             (1-3)

B(Z^-1)＝1+b₁Z^-1+....+b_βZ^-β              (1-4)

假设与混合相位子波脉冲p_t具有相同的振幅谱和自相关函数的最小相位子波延迟脉冲是 ${\hat{p}}_t=\{{\hat{p}}_0,{\hat{p}}_1,\·\·\·\·\·\·{\hat{p}}_N\},$ 它由下式(1-5)给出：

${\hat{p}}_t={\hat{p}}_0{a}_t*b_t---(1-5)$

那么它的Z变换可以表示为：

$\hat{P}\left(Z\right)={\hat{p}}_{0}A\left(Z\right)B\left(Z\right)---(1-6)$

将式(1—6)带入式(1—2)可得：

$P\left(Z\right)=\hat{P}\left(Z\right)\frac{Z^{\mathrm{\β}}B\left(Z^{-1}\right)}{B\left(Z\right)}---(1-7)$

由式(1—7)可知，混合相位子波脉冲是最小相位子波延迟脉冲与一个全通滤波器的褶积。

又由于最小相位子波延迟脉冲

的反向滤波器

的Z变换为：

$\hat{H}\left(Z\right)=\frac{1}{\hat{P}\left(Z\right)}=\frac{1}{{\hat{p}}_0}\frac{1}{A\left(Z\right)}\frac{1}{B\left(Z\right)}---(1-8)$

由此可以得出混合相位子波反滤波器为：

$H\left(Z\right)=\hat{H}\left(Z\right)\frac{B\left(Z\right)}{Z^{\mathrm{\β}}B\left(Z^{-1}\right)}---(1-9)$

即，这个混合相位子波反滤波器等于一个最小相位反滤波器和一个全通滤波器的褶积。

本发明的方法主要包括下列步骤，如图1所示：

S101：采集反射地震记录，选择该反射地震记录中的一条标准道的反射地震记录，该反射地震记录包括实际地震子波的频谱。

为了探明实际地质情况，需要人工激发地震源，而在不同的地层，反射回来的地震记录是不同的，反射地震记录采集装置采集到的反射地震记录包含了各个地震道，也即各个层的地层信息以及地震子波，地震子波在各地震道都是相同的，但却是未知的，但是地层信息在不同地震道却并不相同，为了获取地层信息，本发明首先要从反射地震记录中选取一条标准道的反射地震记录，以便对其进行分析，去掉地震子波，留下地层信息，该反射地震记录就包含了该标准道的地层信息参数以及实际地震子波参数。由于实际地震子波与地层信息是混在一起的，所以很难将地层信息从反射地震记录中分离出来，本发明的目的，就是通过混合反褶积方法将地层信息从反射地震记录中分离出来，而且力求分离出来的地层信息最接近真实的地层环境。

其中，该反射地震记录包含了实际地震子波的参数以及地层信息的参数，如实际地震子波的参数频谱(包括振幅谱和相位谱)、震源等；地层信息的参数深度、反射系数、吸收衰减能力等。

S102：应用谱模拟的方法，在空间窗口和时间窗口中求取与所述实际地震子波频谱相位对应的最小相位子波。

由于本发明的最终目的是通过将求取的地震子波的混合相位反滤波器与反射地震记录进行褶积，进而去掉地震子波，得到高分辨率的地震记录信息，因此，求取地震子波的混合相位反滤波器是关键，而根据上述推导可知，混合相位反滤波器是一个最小相位反滤波器与一个全通滤波器的褶积生成的，因此，为得到混合相位反滤波器，首先要求取地震子波的最小相位反滤波器以及地震子波的全通滤波器，本发明的该步骤就是求取地震子波的最小相位反滤波器的关键步骤，是为了求得与实际地震子波最匹配的混合相位子波，首先提取与实际地震子波频谱相位对应的最小相位子波，进而求得该最小相位子波反滤波器。

根据常用的脉冲反褶积方法可以确定实际地震子波的最小相位反滤波器，同时也能求得与实际地震子波对应的最小相位子波。

谱模拟的方法是求取地震子波频谱相位对应的最小相位子波的常用方法。首先采用谱模拟方法求取地震子波的振幅谱，进而求取地震子波的自相关。过程如下：

对于公式利用最小二乘法模拟出待求的多项式的系数a_n，及其对应的频率ω。是利用最小二乘法求解一个a_n使得

达到最小。由此得到了地震子波的振幅谱，又由地震子波振幅谱和地震子波自相关的关系得到地震子波的自相关，然后在根据复变函数的理论，求出地震子波的相位谱，就可以求得地震子波。

由于地震子波在传播过程中，会由于各种可能的原因，例如阻力的作用，而导致信号衰减，例如随着传播时间的延长，信号越来越弱，或者随着传播深度的增加，信号越来越弱等，因此，实际的地震子波是随时间和空间的变化而不同的，本发明在求取最小相位子波的过程中，应用空间窗口和时间窗口对所求取的最小相位子波进行显示和追踪，更加接近实际地震子波，且更为直观。

S103：采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波，该混合相位子波是时空变化的。

为了求得混合相位子波的反子波滤波器，本步骤将上述步骤中采用谱模拟的方法求取的最小相位子波的最小相位反滤波器分解为一个长度很短(延迟时间长度为β)的滤波器，与一个长度较长的滤波器的褶积，采用Z变换和单位圆变换的方法，将短滤波器变换为最大相位，而长滤波器保持最小相位，那么结果为一个混合相位的反滤波器。对于不同的长度的β(β＝1，2，…)，可得出不同的混合相位反滤波器，所有这些反滤波器都具有相同的振幅谱。

这里首先需要说明的是，本步骤求取混合相位子波的目的也是为了获得混合相位子波的反滤波器，所以，本步骤通过上述具体的方法和步骤直接获得了混合相位子波的多个反滤波器。

其中的Z变换和单位圆变换是为了根据最小相位子波求取全通滤波器而使用，因为根据上述推导可知，混合相位反滤波器等于一个最小相位反滤波器和一个全通滤波器的褶积，因此，得到实际地震子波的最小相位反滤波器和全通滤波器即可得到混合相位反滤波器。由于实际地震子波的最小相位反滤波器通过步骤S102中的脉冲反褶积方法和谱模拟方法已经求得，同时该方法也能求得与实际地震子波对应的最小相位子波。因此，根据(1-9)式，求得全通滤波器就可以求得实际地震子波的混合相位反滤波器。

根据上述推导可知，全通滤波器是与最大相位分量有关的，因此，只需将最小相位子波的单位圆外的零点依次变换到单位圆内，当这个变换关于单位圆对称时，它们的振幅谱保持不变，其变化的只是相位，每一个单位圆对称的零点变换都对应一个不同相位的混合相位子波。因为实际地震子波是有限长度的，所以零点变换有限次可以穷尽，通过本步骤的实施，可以得到多个混合相位子波的反滤波器，但只有一个是与实际地震子波匹配最好的。

其中，Z变换方法的原理如下：

Z变换是对离散序列进行的一种数学变换，将离散序列用级数的形式表示出来。此处，是将地震子波抽样成n个分量的离散序列，再将此离散序列用变量为Z，最高阶为n-1次的多项式表示出来。

其中，单位圆变换方法的原理如下：

根的单位圆变换是引入其他参数，使的在单位圆外(内)的根变换到单位圆内(外)。将根与单位圆的圆心连成一条射线，在射线上找一点，使得此点的值与原根的模是1，那么此点值就是所求根的单位圆对称变换。

例如：x-2＝0，此时根x＝2>1，为根在圆外，如果引入变换x＝1/y，则有1/y-2＝0，此时根y＝1/2<1，根在圆内。其中根x>1为单位圆外，x<1为单位圆内。

由于步骤S102采用时间窗口和空间窗口显示最小相位子波，因此，本步骤求得的混合相位子波也是时空变化的。

S104：采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波。

本步骤是把上述步骤中求得的那些混合相位反滤波器应用于地震数据，并根据P-范数和最小熵准则联合判定的方法，判断所得的不同长度β的混合相位子波反滤波器，以选出最佳的与实际地震子波最匹配的混合相位反滤波器。

选择P-范数和最小熵准则进行联合判定以选取与实际地震子波最匹配的混合相位子波的原因请见如下说明。

首先有两件事必须要考虑到，第一：实际地震子波的零点恰好落在单位圆上，它不在单位圆内又不在单位圆外。第二：与实际地震子波匹配最好的准则怎么定义，显然常用的最小方差是不能用的，因为所有的混合相位子波其振幅谱相同。对于第一件事解决的办法是给子波的自相关函数加一点扰动，使零点偏离单位圆上。第二件事的解决办法是用高阶范数定义匹配准则，porsani和ursin用的是5阶范数，为了让反褶积的结果能量向强反射系数集中，我们选择了最小熵与高阶范数交汇极值点作为判别准则。

根据反射系数序列是白化的，并且与噪声不相关的假设，地震数据(根据地震记录获得)的振幅谱是实际地震子波的振幅谱倍数，地震子波的自相关和地震子波的振幅谱可相互变换，所以求出地震子波的振幅谱，就可计算出地震子波的自相关，为进一步恢复实际地震子波打下基础。

根据本发明的较佳实施例，求取混合相位子波及混合相位子波反滤波器可以通过如下方式实现，即：首先以离散序列表示反射地震记录的振幅谱；然而，使用一个指数函数的多项式在最小平方意义下对所述以离散序列表示的反射地震记录的振幅谱进行拟合，得到的拟合曲线即为地震子波的振幅谱曲线，据此得到地震子波的振幅谱；而后，根据该地震子波的振幅谱求得地震子波的自相关；最后，再根据该地震子波的自相关，求得混合相位子波以及混合相位滤波器。

下面举例说明。需要说明的是，下述实现方式只是举例说明，而非以此作为限制。

用离散序列

来表示地震记录的振幅谱，假设地震子波的振幅谱是光滑的，且消除了仪器、检波器、震源和检波器组合以及震源和虚反射的影响，则可以构造一个指数函数的多项式，使它在最小平方意义下对

进行拟合，并认为这条拟合曲线就相当于地震子波振幅谱的曲线。该多项式的数学表达式为：

式中，L是多项式的阶数，

是频率轴上的坐标值，a_n为待求的多项式的系数。则令：

解这个最小二乘问题，可以求解出多项式的表达式，即地震子波的振幅谱。

地震子波的自相关和地震子波的振幅谱之间的关系为：

为地震子波自相关的频谱，经过傅立叶变换到时间域：

这样就得到了地震子波的自相关。

由于地震子波是混合相位子波，所以它可以写成一个最小相位子波和最大相位子波的褶积，这样上述地震子波的自相关R_p(t)可以写成如下的形式

R_p(t)＝(a_t*b_t)*(a_t*b_t)         (2-1)

b_t＝b(-t)

a_t＝a(-t)

当t>β时：

g_t＝b_t*a_t*b_t＝0                 (2-2)

因为b_t和a_t是非因果的，而b_t是因果的，其长度为β，例如，当β＝1时：

a(m)＝{0 a(2) a(1) 1 0 0}

b(m)＝{0 0 b(1) 1 0 0}

b(m)＝{0 0 0 1 b(1) 0}

                 ↑

                t＝0

上面三项褶积的结果为：

a_t*b_t*b_t＝g_t＝{g_-3，g_-2，g_-1，g₀，g₁，0}                  (2-3)

在方程(2-1)两边都同时乘以a^-1 _t，则：

$R_p\left(t\right)*a_t^{-1}={\stackrel{\&OverBar;}{b}}_t*({\stackrel{\&OverBar;}{a}}_t*b_t)---(2-4)$

当t>β时：

$R_p\left(t\right)*a_t^{-1}={\stackrel{\&OverBar;}{b}}_t*({\stackrel{\&OverBar;}{a}}_t*b_t)=0---(2-5)$

上式表明，最小相位子波的逆和原始子波自相关的褶积只含有β个非零项(t>0)。

令： $c_t=a_t^{-1},$ 对于不同的β：

$\underset{j=0}{\overset{\mathrm{Nc}}{\mathrm{\&Sigma;}}}{R}_p(\mathrm{\&tau;}-j)c_j=0$ τ＝β+1，......β+Nc                  (2-6)

通过解这个矩阵方程(2-6)，就可以得到不同β时的最小相位子波的反子波的一个近似解，理由是，将(2-3)代入(2-6)式得：

$\underset{j=0}{\overset{\mathrm{Nc}}{\mathrm{\&Sigma;}}}g(\mathrm{\&tau;}-j)f\left(j\right)=0$ τ＝β+1，β+2....，β+Nc               (2-7)

式中f(j)＝a(j)*c(j)，f(0)＝1，f(j)是因果序列，g(m)由(2-3)式定义，当m>β时，g(m)＝0，所以当Nc很大时，f(j)逼近于 $\mathrm{\&sigma;}\left(j\right)=\left\{\begin{array}{cc}1& j=0\\ 0& j\&NotEqual;0\end{array}\right.,$ 这就证明了c(j)逼近a(j)的逆a^-1(j)。由c(j)可求得最小相位子波A(Z)，于是：

$\hat{P}\left(Z\right)={\hat{p}}_{0}A\left(Z\right)B\left(Z\right)---(2-8)$

为β＝0时的最小相位子波，所以B(Z)可以(2-8)求得，这样混合相位的子波为：

$P\left(Z\right)=\hat{P}\left(Z\right)\frac{Z^{\mathrm{\&beta;}}B\left(Z^{-1}\right)}{B\left(Z\right)}---(2-9)$

知道了混合相位的子波，可以得到混合相位子波的反滤波器为：

$H\left(Z\right)=\frac{1}{A\left(Z\right)}\frac{1}{Z^{\mathrm{\&beta;}}B\left(Z^{-1}\right)}---(2-10)$

对每一个不同的β，就会有一个不同的混合相位子波反滤波器，我们将这个混合相位子波反滤波器同所选取的标准道的反射地震记录进行褶积：

e_β＝x*h_β+n*h_β

  ＝w*r*h_β+n*h_β                         (2-11)

  ＝μ_β*r+n*h_β

实际上，地震子波为有限长度的，因而其零点个数也是有限的，在单位圆内的零点也是有限的，所以只有有限个β，对于所得的e_β是不同的，本发明采取最小熵判别准则进行选取，以从上述众多的混合相位滤波器中选取最佳的混合相位滤波器，下面为推导判别准则的过程，考虑包含相位信息的高阶累计量：

MC(ββ)＝E(e⁴(m)-3E²(e²(m))               (2-12)

对于一非高斯白噪序列e(m)，即k，l充分大时有 $E\left(e_k^2\left(m\right)\right)=E\left(e_l^2\left(m\right)\right),$ 在样本无穷多的条件下，当h_β为实际子波的逆，也就是当μ(m)＝δ(m)时，(2-12)式将取得最大值。在实际中不可能得到无穷多的样本，常用有限截断的估计值代替，这样(2-12)式就变为：

$\mathrm{MC}\left(\mathrm{\&beta;}\right)=\frac{1}{M}[\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^4\left(j\right)-3{\left(\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^2\left(j\right)\right)}^2]---(2-13)$

上式两边都同时除以 $\frac{1}{M}{\left(\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^2\left(j\right)\right)}^2,$ 得：

$\mathrm{MC}\&prime;\left(\mathrm{\&beta;}\right)=\frac{\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^4\left(j\right)}{{\left(\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^2\left(j\right)\right)}^2}-3---(2-14)$

最大值与(2-14)等价的判别准则为：

$\mathrm{ME}\left(\mathrm{\&beta;}\right)=\frac{\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^4\left(j\right)}{{\left(\underset{j=0}{\overset{M}{\mathrm{\&Sigma;}}}{e}_{\mathrm{\&beta;}}^2\left(j\right)\right)}^2}---(2-15)$

将每一个不同的β所得到的ME(β)和β＝0时相比较：

$\mathrm{\&Phi;}\left(\mathrm{\&beta;}\right)=\frac{\mathrm{ME}\left(\mathrm{\&beta;}\right)}{\mathrm{ME}\left(0\right)}---(2-16)$

得到相对最小熵判别准则(2-16)。对每一个不同的β，计算其混合相位地震子波和反子波滤波器，计算出它们的Φ(β)，选择一个使Φ(β)为最大的滤波器，作为最佳滤波器，其对应的子波作为原地震道的子波。

其中，P-范数判定方法是：

对于反褶积后的地震道e_β，其中Lp范数为 $\mathrm{\&phi;}\&prime;(\mathrm{\&beta;},p)={\left\{\mathrm{\&Sigma;}\right|e_{\mathrm{\&beta;}}{|}^p\}}^{\frac{1}{p}},$ 将此式标准化为： $\mathrm{\&phi;}(\mathrm{\&beta;},p)=\frac{\mathrm{\&phi;}\&prime;(\mathrm{\&beta;},p)}{\mathrm{\&phi;}\&prime;(0,p)},$ 对于给定的p值，我们求出φ(β，p)的值，最大值对应的β延迟滤波器即为最佳滤波器。

其中，最小熵判定方法是：

考虑公式(2-14)(2-15)(2-16)，对每一个不同的β，计算其混合相位地震子波和反子波滤波器，计算出它们的Φ(β)，选择一个使Φ(β)值为最大的滤波器，作为最佳滤波器，其对应的子波作为原地震道的子波。

采用P-范数和最小熵联合判断，可以将上述两种判别准则中不同β所对应的φ(β，p)和Φ(β)的值在同坐标中绘制成曲线图，两条曲线的交汇点的极值点所对应的β就是我们选择的长度，与β相对应的滤波器就是最佳滤波器。

实际上要穷尽所有的混合相位子波的做法要花费大量的计算时间，利用现有技术中普遍使用的Yule-Walker双循环解Zoeppritz矩阵方法可快速求解，在此不再赘述。

S105：采用多层可变时窗，将上述选定的混合相位子波反滤波器与所述反射地震记录褶积，得到所述标准道的地层信息。

本步骤中，将上述步骤中根据P-范数与最小熵判别准则选定的混合相位子波反滤波器与所述反射地震记录褶积，是常规的反褶积方法，用以得到所述标准道的地层信息。

而根据本发明，本步骤在进行常规反褶积的过程中，采用了多层可变时窗进行显示，适应了地震子波时空变化的特点，更符合实际情况。

本发明可以带来如下有益效果：

1、本发明首先在地震道中选择一个标准道，然后应用谱模拟的方法，计算出地震子波的自相关，在这一步中，参数的选择很重要，子波子相关拟合的好坏直接关系到以后褶积的效果。

2、在计算反滤波器的过程中，在默认值附近选取适当的反滤波器的算子长度，以期得到最好的反褶积效果。

3、在计算过程中，可能会遇上子波Z变换多项式的根恰好落在单位圆上的情况，这时双循环求解Zoeppritz矩阵方程组将发散，通常在矩阵的主对角线加上一点扰动，使多项式的根偏离单位圆，将立即会得到稳定的解。

总之，混合相位反褶积方法的研究表明，这种方法是提高地震资料分辨率有效和实用的处理技术，理论模型资料和叠前及叠后的实际资料的处理，都取得了很好的效果。

下面以混合相位反褶积的数值模型以及实际地震资料混合相位反褶积举例来说明应用本发明的有益效果。

混合相位反褶积的数值模型试例：

范例一：理论地震道模型

请参照图2A所示的反射系数序列图，本实施例共有六个反射系数，图2B所示即为图2A的反射系数序列与图3B所示的原始混合相位子波的褶积地震道模型，其纵轴为时间轴，单位ms；横轴为地震道，单位10m。

应用谱模拟的方法，根据地震道的振幅谱，模拟出地震子波的振幅谱，如图3A所示，其为模拟出来的地震子波和实际给出的混合相位子波的振幅谱。其中，A为实际子波的振幅谱，B为模拟子波振幅谱。

再请参照图3B，其中，A为原始混合相位子波，B为求取的混合相位子波，C为求取的最小相位子波，从图3B可以看出，本发明的混合相位反褶积法恢复出来的子波和原始子波基本上是同一个子波，而最小相位法恢复出来的子波则和原始子波误差较大。

请参照表1—1，其显示了对本实施例不同β所求得的Φ(β)，当β取6时Φ(β)得到最大值。

表1-1：不同β所求得的Φ(β)

 

β	0	1	2	3	4	5	6	7	8
										Φ(β)	1.000	0.997	0.682	0.414	0.985	0.300	1.66	0.945	1.650

请参照图4A、图4B，其为反褶积后的输出地震道示意图，其中图4A为最小平方反褶积结果，而图4B为混合相位反褶积结果，纵轴为时间轴，单位ms，横轴为地震道，单位10m。从图中可以看到，利用本发明的混合相位反褶积后的地震道的结果优于最小相位反褶积后的结果。

在理想情况下，本发明的混合相位反褶积将地震子波压缩成尖脉冲，尖脉冲时间位置与地层界面位置一致，没有任何子波旁瓣产生。

反之，由于地震子波是混合相位的，用最小相位子波反褶积产生的旁瓣形成新的规则噪音，并且反射系数越大，产生的旁瓣就越多。

范例二：时变地震道模型

请参照图5，其是给定的地下反射系数序列示意图，形成图6所示的模型地震道。浅部3层是用主频为50Hz的雷克子波与反射系数褶积成。深部3层是用主频为30Hz的雷克子波与反射系数合成的。地震道的子波是时变的。

再请参照图7A所示的单时窗提取子波的频谱图，以及图7B所示的时变地震子波地震道的实际频谱示意图。从图中可以看出，尽管模型地震剖面图6没有加随机噪音，但整个地震道的频谱曲线仍然变化复杂，很不光滑。拟合后的子波频谱见图7A中间的一条频谱曲线，虽然它含有高频成分，但其主要频率部分还是低频成分，并且其主频与相对较低主频(此处为30Hz)的雷克子波的主频非常接近。从此子波的频谱特征来看，用此单时窗提取的子波合成地震记录的分辨率必然接近于用相对较低主频(此处为30Hz)的雷克子波合成的地震记录的分辨率，也就是说分辨率相对较低。

而用分时窗从模型剖面中提取地震子波(见图8)与实际给定的子波更接近。此处分两个时窗，时窗1包含相对浅层的3个反射界面；时窗2包含相对深层的3个反射界面。图8是分时窗提取子波的振幅谱图。时窗1提取的子波频谱与主频为50Hz的雷克子波的频谱吻合较好；时窗2提取的子波频谱与主频为30Hz的雷克子波的频谱吻合较好。用这两个子波合成的地震记录与模型吻合较好。从图9A及图9B可见，单时窗相位反褶积的地震剖面分辨率明显比时变混合相位反褶积后的地震剖面低，因此，分时窗提取的子波更适用于时变子波的地震记录的反褶积。

实际地震资料混合相位反褶积试例：

范例一：实际叠前炮集地震资料的混合相位反褶积。

请参照图12A，其是某地区的一个单炮记录，该记录是2ms采样，40道接收，1501个采样点，根据本发明的方法流程，首先通过地震道的振幅谱拟合出子波的振幅谱，图10是拟合出来的子波的振幅谱和原始地震道的振幅谱。

利用最小相位反褶积法和混合相位反褶积法得到的反褶积算子，恢复出来的子波如图11A、11B所示，其中，图11A为用最小相位反褶积法得到最小相位子波，图11B为用最小相位反褶积法得到混合相位子波，用最小相位反褶积(Promax地震处理软件)和混合相位反褶积这两个不同的反褶积算子对整个地震道进行褶积，得到如图12B和图12C所示的效果，从图中可以看出，利用混合相位反褶积法得到的地震道高频成分明显增多，同相轴变细，连续性提高，能够将复合波分离，分辨率较原始地震道和利用最小相位反褶积后的地震道都有改进。其中，图12A—图12C的纵轴为时间轴，单位为ms，横轴为地震道，单位为25m。

范例二：实际叠后地震剖面的混合相位反褶积。

请参照图15A，其为某地区的叠后资料，该剖面共200道，道间距25米，时间长度2000毫秒。经过反褶积处理后，对褶积后地震道进行了分析，地震道的频谱有所拓宽，提取的实际地震子波见图14A以及图14B。图13为振幅谱示意图，其中A为地震道的振幅谱，B为模拟出来的地震子波的振幅谱。

从图15A-图15C比较实际资料和通过最小相位(Promax地震处理软件)和混合相位滤波器处理后的资料，可以看到，经过处理后，在1100-120ms处，存在一个透镜体状的构造，而在原始资料上，不能分析出来。另外在1300ms和1400ms附近，经过最小相位和混合相位处理过后的资料同相轴明显增多，但是利用混合相位反褶积处理出来的剖面更清晰，分辨率更高。其中，图15A—图15C的纵轴为时间轴，单位为ms，横轴为地震道，单位为25m。

范例三：海洋实际地震资料的混合相位反褶积。

由图16的海洋某地区实际叠后地震剖面图可以看到在3000ms的地方，有一反射界面。如果整道提取一个子波进行反褶积，就会出现以下的问题。在进行谱模拟提取地震子波的时候，如果给的频带范围很宽，浅层的分辨率提高了，但3000ms处的反射界面就很不清楚，见图17所示的单时窗反褶积图，在该图中，谱模拟频率范围为1-150，由图中可以看出，浅层的分辨率有提高，看到3000ms的反射界面。如果频带给的较窄，深层反射界面较清楚，但浅层的分辨率却提高的不够，见图18所示的单时窗反褶积图，在该图中，谱模拟频带范围为1-100，由图中可以看到3000ms的反射界面，但是浅层的分辨率没有提高。这说明原地震剖面子波时变现象较严重。我们用分时窗解决这个时变问题。在浅层，可以给较宽的频带范围，在深层，给相对窄一些的频带范围，这样，浅层分辨率可以提高很多，深层反射界面也可看的较清楚，处理效果见图19所示的分时窗反褶积图，在该图中，浅层谱模拟时频带范围为1-150Hz，深层为1-100Hz，从图中可以看出，与单时窗相比，不仅浅层分辨率提的较高，而且深层反射轴也较清楚了。具体请参照图20A-图20D所示的地震剖面混合相位反褶积处理的局部放大图，其中，图20A为原剖面；图20B的谱模拟频率范围为1—150，由图中可以看出，浅层好，深部效果不大；图20C的谱模拟频带范围为1—100，由图中可以看出，浅层的分辨率没有提高；图20D的频带范围，在浅层谱模拟时为1—150Hz，深层1—100Hz，从图中可以看出，浅、深层分辨率都提高了。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种用于获取地层信息的混合相位反褶积方法，其特征在于，该方法包括下列步骤：

采集反射地震记录，选择该反射地震记录中的一条标准道的反射地震记录，该反射地震记录包括实际地震子波的频谱；

应用谱模拟的方法，在空间窗口和时间窗口中求取与所述实际地震子波频谱相位对应的最小相位子波；

采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波，该混合相位子波是时空变化的；

采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波；

求得所述选定的混合相位子波的反滤波器；

采用多层可变时窗，将上述混合相位子波反滤波器与所述反射地震记录褶积，得到所述标准道的地层信息；

其中，采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波的步骤包括：

计算所述最小相位子波的反滤波器；

将所述最小相位子波反滤波器分解为一个短滤波器与一个长滤波器的褶积，其中，短滤波器的长度为β；

采用Z变换和Z多项式根的单位圆对称投影变换的方法，将短滤波器变换为最大相位，而长滤波器保持最小相位，得到一个混合相位的反滤波器；

选取β为不同的值，得到多个混合相位反滤波器，所有这些反滤波器都具有相同的频谱。

2.根据权利要求1所述的方法，其特征在于，所述采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波的步骤包括；

根据反射地震记录的振幅谱，通过谱模拟求得地震子波的振幅谱；

根据所述地震子波的振幅谱求得最小相位地震子波；

根据所述最小相位地震子波的相位的单位圆对称投影变换，求得混合相位子波以及混合相位滤波器。

3.根据权利要求2所述的方法，其特征在于，所述根据反射地震记录的振幅谱，求得地震子波的振幅谱的步骤包括：

将反射地震记录的振幅谱以离散序列表示；

使用一个指数函数的多项式在最小平方意义下对所述离散序列进行拟合，得到一条拟合指数多项式曲线；

根据所述拟合指数多项式曲线得到地震子波的振幅谱曲线。

4.一种混合相位反褶积处理***，用于通过分析反射地震记录，获取反射地震记录中的地层信息，其特征在于，该***包括：

采集单元，用于采集反射地震记录，选择该反射地震记录中的一条标准道的反射地震记录，该反射地震记录包括实际地震子波的频谱；

最小相位子波计算单元，用于应用谱模拟的方法，在空间窗口和时间窗口中求取与所述实际地震子波频谱相位对应的最小相位子波；

混合相位子波计算单元，用于采用Z变换和单位圆变换的方法，根据所述最小相位子波求得多个混合相位子波，该混合相位子波是时空变化的；

混合相位子波选择单元，用于采用P-范数和最小熵准则联合判定的方法，从所述多个混合相位子波中选择与所述实际地震子波最匹配的混合相位子波；

选定的混合相位子波反滤波器计算单元，用于根据所述选定的混合相位子波计算混合相位子波反滤波器；

褶积单元，用于采用多层可变时窗，将上述混合相位子波反滤波器与所述反射地震记录褶积，得到所述标准道的地层信息；

其中，所述混合相位子波计算单元包括：

计算模块，用于计算所述最小相位子波的反滤波器；

分解模块，用于将所述最小相位子波反滤波器分解为一个短滤波器与一个长滤波器的褶积，其中，短滤波器的长度为β；

变换模块，用于采用Z变换和Z多项式根的单位圆对称投影变换的方法，将短滤波器变换为最大相位，而长滤波器保持最小相位，得到一个混合相位的反滤波器；

选择模块，用于选取β为不同的值，得到多个混合相位反滤波器；

其中，所有这些反滤波器都具有相同的频谱。

Images