CN103018775B - 基于相位分解的混合相位子波反褶积方法 - Google Patents

Abstract

基于相位分解的混合相位子波反褶积方法，应用于地震勘探地震信号处理。特征：人工激发并接收地震记录；对地震记录进行噪声压制、反褶积、速度分析、偏移等常规地震资料处理；采用复赛谱分离方法由地震记录振幅谱估算地震子波振幅谱；按照不同的分解比例，得到一系列具有相同振幅谱，不同相位谱的混合相位子波集合；按照反褶积之后方差模的大小，确定实际混合相位地震子波；设定希望输出零相位地震子波；计算反褶积算子；反褶积算子与地震记录褶积；将地震记录绘制为能够反映地下结构特征的地震剖面图像。提高地震记录的分辨率。能够更加清晰地反映地下结构的内幕和细节。

Description

基于相位分解的混合相位子波反褶积方法

技术领域

本发明涉及勘探地震技术领域，特别涉及一种地震信号处理方法，是一种基于相位分解的混合相位子波反褶积方法。

背景技术

地震勘探是一种利用人工地震技术探测地下构造的勘探方法。它按照一定的方式在地表附近激发地震波，产生称之为地震子波的地震信号，地震子波由震源开始向地下传播，遇到地质界面之后，在界面处发生透射和反射，透射的地震子波继续向下传播，而反射子波在界面位置向上传播，其强度和极性依赖于界面反射系数的大小和符号。来自不同深度界面的地震子波以不同的时间到达地表，通过布设在地表的一种称之为检波器的接收装置，接收来自不同深度地质界面的反射波，接收到的数字信号称之为地震记录，因此，地震记录是不同深度地质界面反射子波的叠加。

地下的每一个地质界面都对应于地震记录中的一个反射子波，通过研究地震记录中不同时间的地震子波，可以研究地质界面的空间展布情况，进而预测和分析地下构造特征。但是，由于地震子波是具有一定延续时间的振动信号，当地层较薄时，来自不同界面的地震子波叠合在一起，使得从地震记录上无法分辨地质界面，这种利用地震信号分辨薄层的能力称之为地震信号的分辨率。

目前，基于地震记录褶积模型的反褶积技术是工业界提高地震勘探分辨率的主要方法，褶积模型表示为：

x(t)=w(t)*r(t)

其中，t为反射时间，单位是秒，x(t)是地震记录，若采用压电检波器接收地震信号，其单位是帕，若采用速度检波器接收地震信号，单位是米/秒，w(t)是地震子波，其单位与地震记录一致，r(t)是地层反射系数，无量纲，“*”表示褶积。

反褶积的目的就是从地震记录x(t)中恢复地层反射系数r(t)。由于地震子波也是未知的，褶积模型中同时存在两个未知量w(t)和r(t)，直接利用褶积模型无法从地震记录x(t)中恢复地层反射系数r(t)。

为了从地震记录x(t)中恢复地层反射系数r(t)，常规的反褶积方法对褶积模型进行了两个数学假设，首先假设地震子波w(t)是最小相位子波，其次假设反射系数r(t)具有高斯白噪的特点。在这两个假设条件下，常规反褶积方法可以得到反射系数r(t)在最下二乘意义下的解，并以次来近似实际地层反射系数。

勘探实践证明：实际地震勘探的子波绝大多数为混合相位子波，常规反褶积方法对地震子波所要求的最小相位假设与地震勘探实践相矛盾，这在很大程度上降低了常规反褶积方法提高地震勘探分辨率的能力。

如果能够从地震记录中估算出混合相位子波w(t)，那么褶积模型中只有反射系数r(t)是未知的，就能够利用维纳滤波方法由褶积模型计算反射系数r(t)，从而达到提高地震资料分辨率的目的。因此，从地震记录中估计混合相位子波成为提高地震资料分辨率的关键。

Lazear（1993）提出利用高阶累积量估计子波的方法，由于二阶累量实际上是自相关函数，不包含任何相位信息，而三阶累量对于地震反射系数这样的随机过程而言近似为零，因此只能用四阶累量进行相位估计，而四阶累量的计算不仅算法复杂，其对随机噪声的敏感性也很难满足实际地震资料处理的要求。

Porsani和Ursin(1998,2000)给出了一种混合相位子波估计方法，该方法的基本做法是：假设子波的Z变换在单位园上无零点，在单位园内有多个个零点，通过求解扩展的Yule-Walker方程首先获得最小延迟分量的反滤波器，再通过解两次YW方程来获得混合相位子波的最大延迟分量。由最大、最小分量确定出混合相位子波，然后确定混合相位子波的反滤波器，将反滤波器应用到地震记录上。由于该方法多次求解YW方程，而每次求解都要引入白噪算子以增强算子的稳定性，白噪算子的累计误差严重影响了混合相位子波的估算精度和对实际地震资料的适应性。

发明内容

本发明的目的是：现有的地震记录反褶积方法假设地震子波为最小相位子波，而实际勘探中的地震子波多为混合相位子波，两者之间的矛盾不但制约了反褶积方法提高地震记录分辨率的能力，也很容易引发“剩余相位”的解释陷阱。本发明的目的是提供一种基于相位分解的混合相位子波反褶积方法，实现混合相位子波反褶积，提高地震记录分辨率。

本发明采用的技术方案是：基于相位分解的混合相位子波反褶积方法，基本原理:

混合相位子波反褶积的关键是从地震记录中提取混合相位子波，有了混合相位子波，则混合相位反褶积能够通过常规的确定性反褶积来实现。混合相位子波由振幅谱和相位谱两部分构成，其中振幅谱的估算比较稳定，能够通过传统的复赛谱分离法得到，因此，混合相位子波估算的难点是如何确定相位谱。相对于振幅谱，相位谱对误差比较敏感，现有的相位谱估算方法很难适应实际地震记录的复杂性。

混合相位子波由最大相位分量和最小相位分量两部分构成，最大相位分量和最小相位分量的不同分解比例决定了地震子波的相位特征。在所有的分解比例中，存在这样一种分解比例，利用该比例得到的混合相位子波使得反褶积之后地震记录的方差模最大，则该混合相位子波就是实际地震记录的子波。

基于相位分解的混合相位子波反褶积方法，步骤：

步骤A、人工激发地震波，利用检波器接收来自地下的反射地震信号，形成地震记录并记录在磁带上；

步骤B、从磁带上读取地震记录，对地震记录进行噪声压制、反褶积、速度分析、偏移等常规地震资料处理，得到常规处理之后的地震记录x(t)，其中，t是时间，单位是毫秒，若采用压电检波器接收地震信号，地震记录x(t)的单位是帕，若采用速度检波器接收地震信号，地震记录x(t)的单位是米/秒。对地震记录x(t)进行傅立叶变换，得到地震记录的振幅谱x(f)，其中，f是频率，单位是赫兹，地震记录振幅谱x(f)无量纲。将地震记录的振幅谱x(f)输出到步骤C所在的处理单元；

步骤C、接收步骤B输出的地震记录振幅谱x(f)，采用复赛谱分离方法由地震记录振幅谱x(f)估算地震子波振幅谱w(f)，步骤如下：

步骤C1、对地震记录的振幅谱x(f)取对数之后进行傅立叶反变换，得到复赛谱序列

步骤C2、对复赛谱序列进行低通滤之后再进行傅立叶变换，并取幂指数，得到地震子波的振幅谱w(f)，将其输出到步骤D所在的处理单元。

步骤D、接收步骤C2输出的地震子波振幅谱w(f)，按照不同的分解比例，将其分解为一系列具有相同振幅谱，不同相位谱的混合相位子波集合，过程如下：

对地震子波振幅谱w(f)取对数之后进行傅立叶反变化，得到零相位子波的复赛谱序列设定最小分解比例λ₀＝0.0，最大分解比例λ_m＝1.0，分解比例增量Δλ=0.01，分别计算不同分解比例λ_j=λ₀+(j-1)Δλ，j＝1,2…,101所对应的混合相位子波w_j(t)，具体步骤如下：

（a）按照分解比例λ_i，依据下式得到最大相位的复赛谱序列和最小相位的复赛谱序列

${\stackrel{\&OverBar;}{v}}_j\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}\left(t\right),& t\&GreaterEqual;0\\ 0,& t<0\end{array}\right.---\left(1\right)$

${\stackrel{\&OverBar;}{u}}_j\left(t\right)=\left\{\begin{array}{cc}0,& t\&GreaterEqual;0\\ \stackrel{\&OverBar;}{w}\left(t\right)-{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}(-t),& t<0\end{array}\right.---\left(2\right)$

式中，是最小相位的复赛谱序列，是最大相位的复赛谱序列，零相位子波的复赛谱序列，和的单位与地震记录x(t)的单位一致，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，λ_i是分解比例因子，无量纲，t是时间，单位是毫秒。

（b）分别对和进行傅立叶变换，取幂指数，再进行傅立叶反变换，得到分解比例λ_j所对应的最小相位子波分量v_j(t)和最大相位子波分量u_j(t)；

（c）最小相位子波分量v_j(t)和最大相位子波分量u_j(t)进行褶积，得到分解比例λ_j所对应的混合相位子波w_j(t)；

（d）重复步骤（a）到步骤（c），得到101个具有相同振幅谱，不同相位谱的混合相位子波集合w_j(t),j＝1,2,…,101，将所有子波输出到步骤E所在的处理单元。

步骤E、接收步骤D输出的所有混合相位子波w_j(t),j=1,2,…,101，按照反褶积之后方差模的大小，从101个混合相位子波集合中确定实际混合相位地震子波w(t)。具体步骤为：

步骤E1、计算每个混合相位子波w_j(t)的反滤波器a_j(t)；

步骤E2、反滤波器a_j(t)与地震记录x(t)褶积，得到反褶积之后的地震记录y_j(t)；

步骤E3、采用下式计算反褶积之后地震记录y_j(t)的方差模e_j，

$e_j=\frac{\underset{t}{\mathrm{\&Sigma;}}{y}_j^4\left(t\right)}{{\left(\underset{t}{\mathrm{\&Sigma;}}{y}_j^2\left(t\right)\right)}^2}---\left(3\right)$

式中，e_j为方差模，无量纲，y_j(t)是反褶积之后的地震记录，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，t为反射时间，单位是毫秒，“∑”表示求和运算。

步骤E4、比较所有方差模的大小，最大方差模所对应的混合相位子波即为实际混合相位地震子波w(t)，将该子波输出到步骤G所在的处理单元；

步骤F、接收步骤B输出的地震记录振幅谱x(f)，利用下式确定希望输出零相位子波b(t)，并将其输出到步骤G所在的处理单元。

$b\left(t\right)=\frac{\sin 2\mathrm{\&pi;}{f}_{c}t}{\mathrm{\&pi;t}}---\left(4\right)$

式中，b(t)是希望输出零相位子波，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，t为反射时间，单位是毫秒，f_c是由振幅谱x(f)所确定的地震记录最高有效频率，单位是赫兹。

步骤G、接收步骤E4所输出的实际混合相位地震子波w(t)和步骤F输出的希望输出子波b(t)，按照维纳滤波方法，由这两个子波计算反褶积算子a(t)，并输出到步骤H所在的处理单元；

步骤H、接收步骤G所输出的反褶积算子a(t)，与步骤B中常规处理之后的地震记录x(t)进行褶积运算，得到混合相位子波反褶积之后的地震记录y(t)；

步骤I、将上述流程处理之后的地震记录绘制为能够反映地下结构特征的地震剖面图像。

本发明的核心内容是：克服传统反褶积方法关于地震子波为最小相位的数学假设，直接从地震记录中提取混合相位子波，再利用维纳滤波方法实现混合相位子波反褶积，提高地震记录的分辨率。其中，直接从地震记录中估算混合相位子波是本发明的关键。本发明提取混合相位子波的做法是，首先利用复赛谱分离技术由地震记录的振幅谱估算地震子波的振幅谱，然后，在复赛谱域进行不同比例的相位分解，得到一系列相同振幅谱、不同相位谱的子波集合，在这些混合相位子波集合中，反褶积之后方差模最大的子波即为实际地震记录的混合相位子波。

本发明的有益效果：本发明基于相位分解的混合相位子波反褶积方法，与传统的反褶积方法相比，本方法克服了地震子波为最小相位的数学假设，改进了混合相位反褶积方法的有效性。利用本方法处理之后的地震记录，能够更加清晰地反映地下结构的内幕和细节。本方法已经在多个勘探区块取得了实际试用，提高了地震勘探的精度。

附图说明

图1是某区块常规处理之后的地震记录剖面图像。

图2是地震记录振幅谱。

图3是地震子波振幅谱。

图4是实际混合相位地震子波，在101个相同振幅谱、不同相位谱子波集合中，该子波反褶积之后的方差模最大，即为实际地震子波。

图5是最高有效频率为80Hz的希望输出零相位子波。

图6是利用维纳滤波方法估算的反褶积算子。

图7是经过本发明处理之后的地震记录剖面图像，与图1所示的常规处理地震剖面图像相比，分辨率得到明显改善，更加精细地刻画了地下结构的沉积特征和层序关系。

图8是塔里木油田X工区常规处理之后的地震记录剖面图像。

图9是塔里木油田X工区提取的实际混合相位地震子波。

图10是塔里木油田X工区计算的混合相位反褶积算子。

图11是塔里木油田X工区经过本发明处理之后的地震记录剖面图像,与图8所示的常规处理地震记录图像相比，分辨率得到明显提高，更好好地揭示了强反射界面之间的内幕结构和层序关系。

图12是塔里木油田Y工区常规处理之后的地震记录剖面图像。

图13是塔里木油田Y工区经过本发明处理之后的地震记录剖面图像,与图12所示的常规处理地震图像相比，地震记录的分辨能力明显提高。两图中虚线圆圈所标注的部位是一规模较大的缝洞群反射，在常规处理图像上，看似一个规模较大的缝洞体，在本发明处理之后的地震图像上，分辨为由诸多小洞构成的缝洞群。这一地质认识对于井位部署和油藏开发具有重要的指导意义。

具体实施方式

实施例1：为了使本发明的技术方案和技术优势更加清晰，下面结合某油田的实施例和附图。对本发明基于相位分解的混合相位子波反褶积方法做进一步的详细说明。

步骤A、人工激发地震波，利用检波器探测来自地下的反射地震信号，并记录在磁带上。本实施例具体的野外采集配置为：6条测线同时接收，每条测线布设180个检波器，线间距100米，道间距50米，激发点与检波点的最小距离为50米，记录长度为6000毫秒，采样间隔2毫秒，共激发并接收了667炮地震记录。

步骤B、从磁带上读取地震信号，对地震信号进行噪声压制、反褶积、速度分析、偏移等常规地震信号处理，得到常规处理之后的地震记录x(t)，对地震记录x(t)进行傅立叶变换，得到地震记录的振幅谱x(f)，将地震记录的振幅谱x(f)输出到后面的处理单元。图1是常规处理之后地震记录的剖面图像，图2是地震记录的振幅谱。

步骤C、接收步骤B输出的地震记录振幅谱x(f)，采用复赛谱分离方法由地震记录振幅谱x(f)估算地震子波振幅谱w(f)，步骤如下：

步骤C1、对地震记录的振幅谱x(f)取对数之后进行傅立叶反变换，得到复赛谱序列

步骤C2、对复赛谱序列进行低通滤之后再进行傅立叶变换，并取幂指数，得到地震子波的振幅谱w(f)，将其输出到步骤D所在的处理单元。图3是地震子波的振幅谱。

步骤D、接收步骤C2输出的地震子波振幅谱w(f)，按照不同的分解比例，将其分解为一系列具有相同振幅谱，不同相位谱的混合相位子波集合，过程如下：

对地震子波振幅谱w(f)取对数之后进行傅立叶反变化，得到零相位子波的复赛谱序列设定最小分解比例λ₀＝0.0，最大分解比例λ_m＝1.0，分解比例增量Δλ=0.01，分别计算不同分解比例λ_j＝λ₀+(j-1)Δλ，j＝1,2…,101所对应的混合相位子波w_j(t)，具体步骤如下：

（a）按照分解比例λ_i，依据下式得到最大相位的复赛谱序列和最小相位的复赛谱序列

${\stackrel{\&OverBar;}{v}}_j\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}\left(t\right),& t\&GreaterEqual;0\\ 0,& t<0\end{array}\right.$

${\stackrel{\&OverBar;}{u}}_j\left(t\right)=\left\{\begin{array}{cc}0,& t\&GreaterEqual;0\\ \stackrel{\&OverBar;}{w}\left(t\right)-{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}(-t),& t<0\end{array}\right.$

（b）分别对和进行傅立叶变换，取幂指数，再进行傅立叶反变换，得到分解比例λ_j所对应的最小相位子波分量v_j(t)和最大相位子波分量u_j(t)；

（c）最小相位子波分量v_j(t)和最大相位子波分量u_j(t)进行褶积，得到分解比例λ_j所对应的混合相位子波w_j(t)；

（d）重复步骤（a）到步骤（c），得到101个具有相同振幅谱，不同相位谱的混合相位子波集合w_j(t),j=1,2,…,101，将所有子波输出到步骤E所在的处理单元。

步骤E、接收步骤D输出的所有101个混合相位子波w_j(t)，按照反褶积之后方差模的大小，从101个混合相位子波集合中确定实际混合相位地震子波w(t)。具体步骤为：

步骤E1、计算每个混合相位子波w_j(t)的反滤波器a_j(t)；

步骤E2、反滤波器a_j(t)与地震记录x(t)褶积，得到反褶积之后的地震记录y_j(t)；

步骤E3、采用下式计算反褶积之后地震记录y_j(t)的方差模e_j，

$e_j=\frac{\underset{t}{\mathrm{\&Sigma;}}{y}_j^4\left(t\right)}{{\left(\underset{t}{\mathrm{\&Sigma;}}{y}_j^2\left(t\right)\right)}^2}$

步骤E4、比较所有方差模的大小，最大方差模所对应的混合相位子波即为实际混合相位地震子波w(t)，将该子波输出到步骤G所在的处理单元；

步骤F、接收步骤B输出的地震记录振幅谱x(f)，利用下式确定希望输出零相位子波b(t)，并将其输出到步骤G所在的处理单元。

$b\left(t\right)=\frac{\sin 2\mathrm{\&pi;}{f}_{c}t}{\mathrm{\&pi;t}}$

式中，f_c是由振幅谱x(f)所确定的地震记录最高有效频率，单位是赫兹，本实施例中f_c＝80Hz，图5为本实施例的希望输出零相位子波。

步骤G、接收步骤E4所输出的实际混合相位地震子波w(t)和步骤F输出的希望输出零相位子波b(t)，按照维纳滤波方法，由这两个子波计算反褶积算子a(t)，并输出到步骤H所在的处理单元，图6是计算得到的反褶积算子。

步骤H、接收步骤G所输出的反褶积算子a(t)，与步骤B中常规处理之后的地震记录x(t)进行褶积运算，得到混合相位子波反褶积之后的地震记录y(t)；

步骤I、将上述流程处理之后的地震记录绘制为能够反映地下结构特征的地震剖面图像。图7是按照上述流程处理之后的地震记录剖面图像，与图1所示的常规处理地震剖面图像相比，分辨率得到明显改善，更加精细地刻画了地下结构的沉积特征和层序关系。

实施例2

本实施例为塔里木油田X工区的应用实例，该地区地表为沙丘覆盖，所采集的地震记录经历了地表沙丘的强烈吸收，高频损失严重，地震子波相位呈现典型的非最小相位特征。为此，利用本发明进行了混合相位子波反褶积处理，提高地震资料的分辨率，增强地震信号反映地下结构的能力。

图8是常规处理的地震记录剖面图像，分辨率较低，不能反映强反射界面之间的内幕结构细节。图9是提取的实际地震子波，地震子波呈现典型的混合相位特征。图10是混合相位子波反褶积算子，常规反褶积算子的能量集中在正半轴上，而该反褶积算子的能量在正轴和负轴方向上均有分布，具有更强的压缩子波提高分辨率能力。图11是本发明处理之后的地震记录剖面图像，与图8所示的常规处理地震记录剖面图像相比，分辨率得到明显提高，更好地揭示了强反射界面之间的内幕结构和层序关系。

实施例3

本实施例为塔里木油田Y工区的应用实例，该地区油气储层为碳酸盐岩缝洞型储层，埋深大，非均质性强，地层吸收严重，现有的地震分辨率不能满足碳酸盐岩储层精细解释的地质需求。

图12是常规处理之后地震记录剖面图像，图13是本发明处理之后的地震记录剖面图像，两者相比，本发明处理之后地震记录的分辨能力明显提高。两图中虚线圆圈所标注的部位是一规模较大的缝洞群反射，在常规处理图像上，看似一个规模较大的缝洞体，在本发明处理之后的地震图像上，分辨为由诸多小洞构成的缝洞群。这一地质认识对于井位部署和油藏开发具有重要的指导意义。

Claims (1)

1.一种基于相位分解的混合相位子波反褶积方法，其特征是：

步骤A、人工激发地震波，利用检波器接收来自地下的反射地震信号，形成地震记录并记录在磁带上；

步骤B、从磁带上读取地震记录，对地震记录进行噪声压制、反褶积、速度分析、偏移常规地震资料处理，得到常规处理之后的地震记录x(t)，其中，t是时间，单位是毫秒，若采用压电检波器接收地震信号，地震记录x(t)的单位是帕，若采用速度检波器接收地震信号，地震记录x(t)的单位是米/秒；对地震记录x(t)进行傅立叶变换，得到地震记录的振幅谱x(f)，其中，f是频率，单位是赫兹，地震记录振幅谱x(f)无量纲；将地震记录的振幅谱x(f)输出到步骤C所在的处理单元；

步骤C、接收步骤B输出的地震记录振幅谱x(f)，采用复赛谱分离方法由地震记录振幅谱x(f)估算地震子波振幅谱w(f)，步骤如下：

步骤C1、对地震记录的振幅谱x(f)取对数之后进行傅立叶反变换，得到复赛谱序列

步骤C2、对复赛谱序列进行低通滤之后再进行傅立叶变换，并取幂指数，得到地震子波的振幅谱w(f)，将其输出到步骤D所在的处理单元；

步骤D、接收步骤C2输出的地震子波振幅谱w(f)，按照不同的分解比例，将其分解为一系列具有相同振幅谱，不同相位谱的混合相位子波集合，过程如下：

对地震子波振幅谱w(f)取对数之后进行傅立叶反变化，得到零相位子波的复赛谱序列设定最小分解比例λ₀＝0.0，最大分解比例λ_m＝1.0，分解比例增量Δλ＝0.01，分别计算不同分解比例λ_j＝λ₀+(j-1)Δλ，j＝1,2…,101所对应的混合相位子波w_j(t)，具体步骤如下：

(a)按照分解比例λ_i，依据下式得到最大相位的复赛谱序列和最小相位的复赛谱序列

${\stackrel{\&OverBar;}{v}}_j\left(t\right)=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}\left(t\right),& t\&GreaterEqual;0\\ 0,& t<0\end{array}\right.---\left(1\right)$

${\stackrel{\&OverBar;}{u}}_j\left(t\right)=\left\{\begin{array}{cc}0,& t\&GreaterEqual;0\\ \stackrel{\&OverBar;}{w}\left(t\right)-{\mathrm{\&lambda;}}_j\stackrel{\&OverBar;}{w}(-t),& t<0\end{array}\right.---\left(2\right)$

式中，是最小相位的复赛谱序列，是最大相位的复赛谱序列，零相位子波的复赛谱序列，和的单位与地震记录x(t)的单位一致，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，λ_i是分解比例因子，无量纲，t是时间，单位是毫秒；

(b)分别对和进行傅立叶变换，取幂指数，再进行傅立叶反变换，得到分解比例λ_j所对应的最小相位子波分量v_j(t)和最大相位子波分量u_j(t)；

(c)最小相位子波分量v_j(t)和最大相位子波分量u_j(t)进行褶积，得到分解比例λ_j所对应的混合相位子波w_j(t)；

(d)重复步骤(a)到步骤(c)，得到101个具有相同振幅谱，不同相位谱的混合相位子波集合w_j(t),j＝1,2,…,101，将所有子波输出到步骤E所在的处理单元；

步骤E、接收步骤D输出的所有混合相位子波w_j(t),j＝1,2,…,101，按照反褶积之后方差模的大小，从101个混合相位子波集合中确定实际混合相位地震子波w(t)；具体步骤为：

步骤E1、计算每个混合相位子波w_j(t)的反滤波器a_j(t)；

步骤E2、反滤波器a_j(t)与地震记录x(t)褶积，得到反褶积之后的地震记录y_j(t)；

步骤E3、采用下式计算反褶积之后地震记录y_j(t)的方差模e_j，

$e_j=\frac{\underset{t}{\mathrm{\&Sigma;}}{y}_j^4\left(t\right)}{{\left(\underset{t}{\mathrm{\&Sigma;}}{y}_j^2\right(t\left)\right)}^2}---\left(3\right)$

式中，e_j为方差模，无量纲，y_j(t)是反褶积之后的地震记录，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，t为反射时间，单位是毫秒，“∑”表示求和运算；

步骤E4、比较所有方差模的大小，最大方差模所对应的混合相位子波即为实际混合相位地震子波w(t)，将该子波输出到步骤G所在的处理单元；

步骤F、接收步骤B输出的地震记录振幅谱x(f)，利用下式确定希望输出零相位子波b(t)，并将其输出到步骤G所在的处理单元；

$b\left(t\right)=\frac{\sin 2{\mathrm{\&pi;f}}_{c}t}{\mathrm{\&pi;}t}---\left(4\right)$

式中，b(t)是希望输出零相位子波，若采用压电检波器接收地震信号，单位是帕，若采用速度检波器接收地震信号，单位是米/秒，t为反射时间，单位是毫秒，f_c是由振幅谱x(f)所确定的地震记录最高有效频率，单位是赫兹；

步骤G、接收步骤E4所输出的实际混合相位地震子波w(t)和步骤F输出的希望输出零相位子波b(t)，按照维纳滤波方法，由这两个子波计算反褶积算子a(t)，并输出到步骤H所在的处理单元；

步骤H、接收步骤G所输出的反褶积算子a(t)，与步骤B中常规处理之后的地震记录x(t)进行褶积运算，得到混合相位子波反褶积之后的地震记录y(t)；

步骤I、将上述流程处理之后的地震记录绘制为能够反映地下结构特征的地震剖面图像。