CN102096101B - 提高地震数据处理精度的方法及装置 - Google Patents

Abstract

本发明实施例提供一种提高地震数据处理精度的方法及装置，该方法包括：获取地震数据并计算地震数据的四阶累积量；通过三维滞后窗函数对地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；通过地震子波的三谱提取混合相位地震子波。通过本发明实施例，可以提高了高阶谱估计的精度和子波估计的精度。

Description

提高地震数据处理精度的方法及装置

技术领域

本发明涉及石油地球物理勘探领域，特别涉及一种利用地震记录的三谱实现提高地震数据处理精度的方法及装置。

背景技术

在地震正演模拟、常规反演以及地震反褶积方法中，都需要直接或者间接地利用地震子波。在地震正演问题中，需要利用波动方程或者褶积模型，结合地震子波正演形成地震数值模拟记录；在反演和反褶积问题中，也需要通过地震数据估计子波，地震子波提取的精度对最终反演结果有很大的影响。

地震子波的估计方法可以划分为两类，一类是确定性方法，另外一类是统计性方法。目前工业生产中广泛采用的确定性子波提取方法，就是利用测井资料首先计算出反射系数序列，然后结合井旁地震道由褶积模型理论求地震子波。确定性子波提取方法需要在勘探程度较高，测井资料丰富的勘探区域才可以得到良好的应用，在勘探程度低，井资料缺乏的勘探环境中，它的应用受到很大的限制。同时，井旁道提取的地震子波必定包含地层反射系数的响应信息，不同的井提取的子波往往具有较大的差异，具有很低的可信度。

在没有测井资料或者不能直接测量地震子波的情况下，一般使用统计性子波提取方法，这种方法一般需要对反射系数的分布做一些假设，然后利用地震数据的统计特征提取子波，其假设条件与实际情况的吻合程度对子波提取的精度有很大的影响，但它的优点是即使无测井资料，也可以得到比较准确的子波，可以有更广泛的应用，而实际地震子波往往是混合相位的特性也使得二阶统计量方法提取的子波是不准确的。从物理含义上看，二阶统计最大的缺点是只对振幅信息敏感，丢失了地震子波的相位信息。

为了解决混合相位子波提取问题，Lazear于1993年提出将高阶统计理论用于地震子波提取，Lazear提出了通过将初始地震子波的高阶累积量与地震记录的高阶累积量在最小平方条件下拟合，然后利用最速下降算法估计非最小相位地震子波的方法。但是最速下降法作为线性拟合方法有两个缺点：最终结果受给定的初值影响较大，并且有和初值相像的趋势；结果有可能陷入局部极值，如果初值给的不恰当，所得结果可能会有较大误差。

为了减小误差，Velis和Ulryc提出了一种非线性拟合方法，即模拟退火法，该方法能有效避免结果陷入局部极小的可能，但是模拟退火算法也有它自身的缺点，比如控制参数的选择比较困难，并且运算时间较长。

因此在实现本发明的过程中，发明人发现现有技术的缺陷在于：在地震子波非最小相位的情况下，二阶统计方法往往是不适用的；而高阶相关函数和高阶谱对地震子波的相位信息是敏感的，因此可以采用高阶谱方法实现对混合相位地震子波的估计，然而子波估计的结果往往受到高阶谱估计精度的影响，谱估计的精度不高。

发明内容

本发明实施例提供一种提高地震数据处理精度的方法及装置，目的在于提高高阶谱估计的精度和子波估计的精度。

为达到上述目的，本发明实施例提供一种提高地震数据处理精度的方法，所述方法包括：

获取地震数据并计算地震数据的四阶累积量；

通过三维滞后窗函数对所述地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；所述三维滞后窗函数通过标准的一维滞后窗函数构造；所述三维滞后窗函数为：

其中，所述一维滞后窗函数d(τ)为Parzen窗；

通过所述地震子波的三谱提取所述混合相位地震子波，以进行地震正演模拟、或常规反演、或地震反褶积；其中所述通过所述地震子波的三谱提取所述混合相位地震子波，具体包括：通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱；通过二阶统计方法获取地震子波的振幅谱；根据所述地震子波的相位谱和振幅谱生成所述混合相位地震子波；所述通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱包括：确定所述Pan-Nikias公式：其中，是所述地震子波的相位谱，n＝2,3,...,N，s(n)是所述地震子波的三谱的相位谱；构造线性方程组： $\begin{array}{c}s\left(2\right)=2\mathrm{\φ}\left(1\right)-\mathrm{\φ}\left(2\right)\\ s\left(3\right)=3\mathrm{\φ}\left(1\right)+2\mathrm{\φ}\left(2\right)-\frac{14}{6}\mathrm{\φ}\left(3\right)\\ .\\ .\\ .\\ s\left(N\right)=\mathrm{N\φ}\left(1\right)+(N-1)\mathrm{\φ}\left(2\right)+...+2\mathrm{\φ}(N-1)+\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)\end{array};$ 结合所述Pan-Nikias公式和所述线性方程组得到矩阵方程：Bφ＝s，其中，B是(N-1)×(N-1)维方阵， $s=[s\left(1\right),...,s(N-1),s\left(N\right)-\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)]$ 是所述地震子波的三谱的相位谱的矩阵，φ＝[φ(1),φ(2),...,φ(N-1)]^T是所述地震子波的相位谱的矩阵；计算地震子波的相位谱：φ＝B^-1s。

本发明实施例还提供一种提高地震数据处理精度的装置，其特征在于，所述装置包括：

计算单元，获取地震数据并计算地震数据的四阶累积量；

处理单元，通过三维滞后窗函数对所述地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；所述三维滞后窗函数通过标准的一维滞后窗函数构造；所述三维滞后窗函数为：

其中，所述一维滞后窗函数d(τ)为Parzen窗；

提取单元，通过所述地震子波的三谱提取所述混合相位地震子波，以进行地震正演模拟、或常规反演、或地震反褶积；其中，所述提取单元具体包括：递推单元，通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱；获取单元，通过二阶统计方法获取地震子波的振幅谱；生成单元，根据所述地震子波的相位谱和振幅谱生成所述混合相位地震子波；所述递推单元具体用于，确定所述Pan-Nikias公式：其中，是所述地震子波的相位谱，n＝2,3,...,N，s(n)是所述地震子波的三谱的相位谱；构造线性方程组： $\begin{array}{c}s\left(2\right)=2\mathrm{\φ}\left(1\right)-\mathrm{\φ}\left(2\right)\\ s\left(3\right)=3\mathrm{\φ}\left(1\right)+2\mathrm{\φ}\left(2\right)-\frac{14}{6}\mathrm{\φ}\left(3\right)\\ .\\ .\\ .\\ s\left(N\right)=\mathrm{N\φ}\left(1\right)+(N-1)\mathrm{\φ}\left(2\right)+...+2\mathrm{\φ}(N-1)+\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)\end{array};$ 结合所述Pan-Nikias公式和所述线性方程组得到矩阵方程：Bφ＝s，其中，B是(N-1)×(N-1)维方阵， $s=[s\left(1\right),...,s(N-1),s\left(N\right)-\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)]$ 是所述地震子波的三谱的相位谱的矩阵，φ＝[φ(1),φ(2),...,φ(N-1)]^T是所述地震子波的相位谱的矩阵；计算地震子波的相位谱：φ＝B^-1s。

本发明实施例的有益效果在于，通过三维滞后窗函数对地震记录的三谱进行平滑，从而提高了高阶谱估计的精度和子波估计的精度。

附图说明

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

图1是本发明实施例的混合相位地震子波的提取方法的流程图；

图2是本发明实施例的褶积模型中混合相位特征明显的子波、反射系数序列、合成地震记录的示意图；

图3是本发明实施例的反射系数序列的四阶累积量切片、反射系数三谱振幅谱切片、反射系数三谱相位谱切片的示意图；

图4是本发明实施例的混合相位地震子波的四阶累积量切片、三谱振幅谱切片、三谱相位谱切片的示意图；

图5是本发明实施例的合成地震记录的四阶累积量切片、三谱振幅谱切片、三谱相位谱切片的示意图；

图6是本发明实施例的三维滞后窗在τ₁＝5、τ₁＝40、τ₁＝60、τ₁＝80时的切片结果的示意图；

图7是本发明实施例的未加窗函数优化的地震记录的三谱振幅谱切片、三谱相位谱切片的示意图；

图8是本发明实施例的加窗函数优化的地震记录的三谱振幅谱切片、三谱相位谱切片的示意图；

图9是本发明实施例的混合相位子波提取结果的时移前、时移后的示意图；

图10是本发明实施例的未对地震记录的三谱进行优化直接进行子波提取、对地震记录的三谱进行优化的结果示意图；

图11是本发明实施例的雷克子波、改进的Morlet子波和混合相位特征明显的子波的波形示意图；

图12是本发明实施例的数值模拟利用的反射系数序列的示意图；

图13是本发明实施例的子波估计长度等于真实子波长度时子波提取结果的示意图；

图14是本发明实施例的子波估计长度小于真实子波长度时子波提取结果的示意图；

图15是本发明实施例的子波估计长度大于真实子波长度时子波提取结果的示意图；

图16是本发明实施例的合成地震记录的信噪比为30dB的子波估计结果的示意图；

图17是本发明实施例的合成地震记录的信噪比为20dB的子波估计结果的示意图；

图18是本发明实施例的合成地震记录的信噪比为10dB的子波估计结果的示意图；

图19是本发明实施例的混合相位地震子波估计的实际地震数据的示意图；

图20是根据图19中5道地震数据利用三谱法提取的子波估计结果的示意图；

图21是本发明实施例的平均子波的频谱及其Z变换后根的分布示意图；

图22是本发明实施例的混合相位地震子波的提取装置的构成图；

图23是本发明实施例的提取单元的构成图。

具体实施方式

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

本发明实施例提供一种混合相位地震子波的提取方法，如图1所示，该方法包括：

步骤101，计算地震数据的四阶累积量；

步骤102，通过三维滞后窗函数对地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；

步骤103，通过地震子波的三谱提取混合相位地震子波。

在本实施例中，根据褶积模型，地震记录x(t)可以看作是地震子波w(t)与反射系数序列r(t)的褶积，即：

$x\left(t\right)=\underset{\mathrm{\τ}}{\mathrm{\Σ}}w\left(\mathrm{\τ}\right)r(t-\mathrm{\τ})=w\left(t\right)*r\left(t\right)---\left(1\right)$

其中，w(t)是地震子波，r(t)是反射系数序列。为表述方便，可以将式(1)改写成离散形式：

$x\left(n\right)=\underset{i=0}{\overset{q}{\mathrm{\Σ}}}w\left(i\right)r(n-i)---\left(2\right)$

其中，地震记录x(n)是在加性高斯噪声v(n)中观测的，因此：

y(n)＝x(n)+v(n)                 (3)

y(n)是实际记录的地震数据，假设满足以下条件：(1)反射系数序列r(n)是平稳的、零均值的独立同分布的非高斯过程；(2)地震子波允许是非最小相位的；(3)噪声v(n)与反射系数序列r(n)相互独立。

在步骤101实施时，根据高阶谱的定义，对于长度为N的不含噪声的零均值的地震信号x(n)，其三谱可以定义为：

T_x(ω₁,ω₂,ω₃)＝X(ω₁)X(ω₂)X(ω₃)X^*(ω₁+ω₂+ω₃)        (4)

其中，X(ω)是地震信号x(n)的傅里叶变换，上式又可以进一步表示为：

T_x(ω₁,ω₂,ω₃)＝|T_x(ω₁,ω₂,ω₃)|exp[jψ(ω₁,ω₂,ω₃)]        (5)

根据累积量的定义及性质，含有加性噪声v(n)的地震数据y(n)，其四阶累积量可由下式表示：

C_4y(τ₁,τ₂,τ₃)＝M_4w(τ₁,τ₂,τ₃)C_4r(τ₁,τ₂,τ₃)+C_4v(τ₁,τ₂,τ₃)   (6)

其中，M_4w是地震子波的四阶矩，其表示形式为：

$M_{4w}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2,{\mathrm{\τ}}_3)=\underset{n}{\mathrm{\Σ}}w\left(n\right)w(n+{\mathrm{\τ}}_1)w(n+{\mathrm{\τ}}_2)w(n+{\mathrm{\τ}}_3)---\left(7\right)$

由于噪声的四阶累积量C_4v趋于零，因此观测数据y(n)的四阶累积量与无噪声地震数据x(n)的四阶累积量相同。当反射系数为非高斯、独立同分布且数据趋于无穷长度时，反射系数的四阶累积量为常数，反射系数的峰度用δ_r表示，因此利用累积量的BBR公式，可以得到：

C_4y(τ₁,τ₂,τ₃)＝δ_rM_4w(τ₁,τ₂,τ₃)             (8)

其中，δ_r＝C_4r(0,0,0)。

根据(8)式，可得到地震记录的三谱与地震子波的三谱满足如下关系：

T_y(ω₁,ω₂,ω₃)＝δ_rW(ω₁)W(ω₂)W(ω₃)W(-ω₁-ω₂-ω₃)         (9)

反射系数的峰度δ_r是标量，它的存在不会影响地震子波的形态，只会影响估计子波的振幅大小。因此通过计算地震记录的三谱可以得到地震子波三谱的估计，进而可以估计地震子波。

下面以模型为例说明褶积模型中地震子波、反射系数以及合成地震记录的高阶统计量特征。

图2中(a)是混合相位特征明显的子波示意图，包括40个采样点；图2中(b)是反射系数序列示意图，满足贝努力-高斯分布，包括500个采样点；图2中(c)是合成地震记录示意图。

在本实施例中，反射系数的四阶累积量切片和相应的三谱切片(τ₁＝0)见图3。图3中(a)为图2中(b)反射系数序列的四阶累积量切片的示意图，其中，在计算累积量之前已对输入序列去均值；图3中(b)(c)分别为反射系数三谱振幅谱切片和相位谱切片(τ₁＝0)的示意图。

由图3中(a)可见，反射系数的四阶累积量在中心处有正尖峰，即δ_r＝C_4r(0,0,0)>0，表明反射系数序列是超高斯分布的；由图3中(b)(c)可见，反射系数的三谱在全空间都有展布。

在本实施例中，混合相位地震子波的四阶累积量切片(τ₁＝0)见图4中(a)，合成地震记录的四阶累积量切片(τ₁＝0)见图5中(a)。由图5中(a)可见，合成记录的四阶累积量是子波与反射系数的累积量的综合结果。

图4中(b)和图5中(b)分别是地震子波和合成记录的振幅谱切片(τ₁＝0)，图4中(c)和图5中(c)分别是地震子波和合成记录的相位谱切片(τ₁＝0)。通过观察发现，合成记录的三谱仍然是反射系数的三谱与地震子波的三谱综合的结果，合成记录的三谱与混合相位子波的三谱的空间展布基本一致，尽管频带宽度发生了变化，但主频是完全对应的，这表明合成记录的三谱能够比较准确地反映子波的频谱。

实际上，地震记录中的噪声的四阶累积量并不能总满足等于零这个条件，同时，反射系数的四阶累积量在零延迟时也不是多维尖脉冲。因此在地震数据较少或者地震记录含有加性高斯噪声的情况下，要提高三谱估计的精度，必须对谱估计过程进行优化处理。

在步骤102实施时，在滞后域构造三维滞后窗函数对地震记录的四阶累积量进行平滑处理，进而提高三谱估计的精度。

三维滞后窗函数可用标准的一维滞后窗d(τ)来构造，其中d(τ)应该满足：

$\left\{\begin{array}{c}d\left(\mathrm{\τ}\right)=d(-\mathrm{\τ})\\ d\left(\mathrm{\τ}\right)0,\left|\mathrm{\τ}\right|>L\\ d\left(0\right)=1\\ D\left(\mathrm{\ω}\right)=\underset{\mathrm{\τ}=-L}{\overset{L}{\mathrm{\Σ}}}d\left(\mathrm{\τ}\right)\exp -(-\mathrm{j\ω\τ})\≥0,\left|\mathrm{\ω}\right|\≤\mathrm{\π}\end{array}\right.---\left(10\right)$

许多一维滞后窗都可用于三维滞后窗的构造，例如，Daniel窗、Hamming窗、Parzen窗和Sasaki窗等。但不限于此，可根据实际情况确定具体的实施方式。

其中Sasaki窗是一种常用的滞后窗，一维Sasaki窗的数学表达式为：

$d_s\left(\mathrm{\τ}\right)=\left\{\begin{array}{c}\frac{1}{\mathrm{\π}}\left|\sin \frac{\mathrm{\π\τ}}{L}\right|+(1-\frac{\left|\mathrm{\τ}\right|}{L})\cos \frac{\mathrm{\π\τ}}{L},\left|\mathrm{\τ}\right|\≤L\\ 0,\left|\mathrm{\τ}\right|>L\end{array}\right.---\left(11\right)$

另外一种常用的一维滞后窗为Parzen窗，其表达式为：

$d_p\left(\mathrm{\τ}\right)=\left\{\begin{array}{cc}1-6{\left(\frac{\left|\mathrm{\τ}\right|}{L}\right)}^2+6{\left(\frac{\left|\mathrm{\τ}\right|}{L}\right)}^3,& \left|\mathrm{\τ}\right|\≤\frac{L}{2}\\ 2{(1-\frac{\left|\mathrm{\τ}\right|}{L})}^3& \frac{L}{2}\≤\left|\mathrm{\τ}\right|\≤L\\ 0,& \left|\mathrm{\τ}\right|>L\end{array}\right.---\left(12\right)$

研究人员发现在许多一维窗函数中，Parzen窗得到的效果是最好的，因此，在本实施例中采用Parzen窗构造三维平滑窗函数。三维滞后窗的构造可以根据通用多维滞后窗构造公式完成，通用多维窗函数构造公式如式(13)所示：

在本实例中，构造的三维滞后窗函数为：

其中d(τ)为Parzen窗。

在本实施例中，利用上述公式设计的三维滞后窗如图6所示，其中图6中(a)(b)(c)(d)分别是三维滞后窗在τ₁＝5，τ₁＝40，τ₁＝60，τ₁＝80时的切片结果的示意图。

在本实施例中，为了验证了三维Parzen窗在谱估计过程中的作用，对图2中(c)合成地震记录的三谱进行了加窗平滑处理，优化前后合成地震记录的三谱(τ₁＝0)分别如图7和图8所示。

图7中(a)为未加窗函数优化的地震记录的三谱振幅谱切片的示意图；图7中(b)为未加窗函数优化的地震记录的三谱相位谱切片的示意图；图8中(a)为加窗函数优化的地震记录的三谱振幅谱切片的示意图；图8中(b)为加窗函数优化的地震记录的三谱相位谱切片的示意图。

比较上述附图可以明显发现，由于地震数据是有限长度的，反射系数不严格满足非高斯分布等原因引起的误差导致图7A中三谱振幅谱振动剧烈，而经过三维滞后窗的平滑优化处理，振幅谱已经得到明显的平滑，从而提高了谱估计的精度。

同时，通过观察优化前后的三谱相位谱切片也可以发现，相位谱相对未优化前已经变得更加光滑。因为三谱相位谱是高阶统计量，再加上相位表达和计算本身所受的影响因素较大，所以从图中无法观察出它包含的内在信息，虽然它表现的凌乱无章，但是它却比二阶统计量(自相关，功率谱)包含了更多的信息。

在步骤103实施时，通过地震子波的三谱得到混合相位地震子波，具体包括：通过Pan-Nikias公式由地震子波的三谱的相位谱递推出地震子波的相位谱；通过二阶统计方法获取地震子波的振幅谱；根据地震子波的相位谱和振幅谱生成混合相位地震子波。

下面具体介绍三谱域相位重构算法，三谱相位重构就是在已经求得地震子波的三谱T_w(ω₁,ω₂,ω₃)情况下，重构地震子波的相位谱。

根据三谱的定义，地震子波的三谱可定义为：

T_w(ω₁,ω₂,ω₃)＝W(ω₁)W(ω₂)W(ω₃)W^*(ω₁+ω₂+ω₃)        (15)

上式可进一步表示为：

地震子波的频谱可表示为：

W(ω)＝|W(ω)|exp[jφ(ω)]                    (17)

根据(15)，(16)和(17)式，可得到地震子波的三谱相位与地震子波相位之间的关系式：

利用上述基本方程，可得到递推的连续形式：

和离散形式：

其中n＝2,3,...,N。

上式(20)构成了三谱域内求解地震子波相位谱的MU递推算法，该算法的初始值φ(0)＝0和

$\mathrm{\φ}\left(1\right)=\underset{n=2}{\overset{N}{\mathrm{\Σ}}}\frac{Q\left(n\right)-Q(n-1)}{n(n-1)}+\frac{\mathrm{\φ}\left(N\right)}{N}-\frac{N-1}{N}\mathrm{\φ}\left(0\right)---\left(21\right)$

其中

φ(N)可以根据的值取0或者π。

Pan与Nikias通过改进上述算法，提出了三谱域相位重构的最小二乘算法，即Pan-Nikias公式，这种方法首先通过定义：

从而构造线性方程组：

$\begin{array}{c}s\left(2\right)=2\mathrm{\φ}\left(1\right)-\mathrm{\φ}\left(2\right)\\ s\left(3\right)=3\mathrm{\φ}\left(1\right)+2\mathrm{\φ}\left(2\right)-\frac{14}{6}\mathrm{\φ}\left(3\right)\\ .\\ .\\ .\\ s\left(N\right)=\mathrm{N\φ}\left(1\right)+(N-1)\mathrm{\φ}\left(2\right)+...+2\mathrm{\φ}(N-1)+\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)\end{array}---\left(24\right)$

结合公式(23)和(24)，可得到以下矩阵方程：

Bφ＝s                      (25)

其中

φ＝[φ(1),φ(2),...,φ(N-1)]^T               (26)

$s=[s\left(1\right),...,s(N-1),s\left(N\right)-\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)]---\left(27\right)$

在(28)式中，B是(N-1)×(N-1)维方阵。由于B是非奇异的，故可以通过直接求矩阵逆的方式计算地震子波的相位谱：

φ＝B^-1s                    (29)

在本实施例中，可利用图2所示褶积模型数据，进行混合相位子波提取。图9中(a)是混合相位子波提取结果的示意图；因为相位谱计算存在线性相移，所以在时间域需要进行时移，时移后的结果如图9中(b)所示。从中发现，提取的子波与原始子波非常相似，从而验证了这种方法的有效性。

在本实施例中，通过加多维滞后窗对地震记录的三谱进行优化，对最终的地震子波提取结果有多大的影响，也要通过子波提取进行验证。图10中(a)是未对地震记录的三谱进行优化，直接进行子波提取的结果示意图，提取结果与原始子波的相似度为0.85645；而如图10中(b)所示，经过优化处理后，子波提取结果与原始子波的相似度提高到了0.93967，提取的子波在波形前端和尾段的质量都有了明显提高，与真实子波在波形上更加逼近。

由上述实施例可知，利用三谱法估计混合相位地震子波的方法是建立在信号高阶统计理论之上的，高阶谱对地震子波的振幅和相位信息都是敏感的，它抛开了传统的地震子波最小相位、反射系数随机白噪假设，而二阶统计只对振幅信息敏感，因此通过计算地震记录的高阶谱，可以将真实子波的相位和振幅提取出来，从而达到混合相位子波恢复的目的。

在本实施例中，为了测试利用地震记录的三谱进行混合相位子波估计方法的有效性，可以分别利用雷克子波、改进的Morlet子波和混合相位特征明显的子波三种子波进行数值模拟试验，并将估计的子波与原始子波进行对比分析。

首先，介绍数值模拟采用的三种子波及反射系数序列。

(1)雷克子波，其解析表达式为：

$w\left(t\right)=[1-2{\left(\mathrm{\π}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\π}{f}_{p}t\right)}^2}---\left(30\right)$

其中f_p为雷克子波的主频。

(2)改进的Morlet子波，其解析表达式为：

$w\left(t\right)=e^{\mathrm{imt}}{e}^{-\frac{1}{2}{\left(\mathrm{ct}\right)}^2}---\left(31\right)$

其中m＝2πf为子波的中心角频率，i为虚数单位，c为常数，用以控制高斯函数从而调制小波函数，控制子波波形，公式(31)给出的是复Morlet小波的解析表达式，本模拟中只利用了其实部。

(3)混合相位特征明显的子波，该子波无明确的解析表达式。

上述三种子波的波形如图11所示，其中横坐标为子波的延续时间，纵坐标为归一化振幅。

如图11所示，(a)给出的雷克子波主频为40Hz，采样间隔为2ms，38个采样点，延续时间为76ms；(b)给出的Morlet子波的中心角频率为4π，采样间隔为2ms，40个采样点，延续时间为80ms；(c)给出的混合相位子波主频为50Hz，采样间隔为2ms，40个采样点，延续时间为80ms。

数值模拟利用的反射系数序列如图12所示，反射系数序列服从[-0.4,0.4]之间的均匀分布，采样间隔为2ms，1000个采样点，延续时间为2s。

在本实施例中，还进行无噪声子波估计试验及分析。在实际地震资料处理中，由于地震子波的延续时间是未知的，因此往往需要进行大量的试验工作，以测试地震子波的长度。不同的子波估计长度，对提取的地震子波以及后期的地震反褶积处理具有很大的影响。

因此，利用无噪声条件下的合成地震数据，对上述介绍的三种子波分别进行子波提取试验，并测试估计子波长度与真实子波的长度不同时，最终子波估计结果的准确性。

首先，测试当子波估计长度等于真实子波长度时，子波提取的准确性。将雷克子波、改进的Morlet子波和混合相位特征明显的子波的估计长度设定为真实的子波长度，即将子波估计长度参数设为76ms、80ms和80ms。在此条件下子波提取结果如图13所示，其中点线为估计子波，实线为原始子波。从中可以发现，当子波估计长度与真实长度一致时，利用该方法得到的结果与原始子波在形态上有较好的一致性。

其次，当子波估计长度不等于真实子波长度时，模拟结果如图14和图15所示，其中图14中三种子波的预设长度均小于真实子波的长度，三种子波的估计长度参数分别设为60ms、76ms和72ms。图15中三种子波的估计长度均大于真实子波的长度，分别设置为80ms、84ms和88ms，点线为估计子波，实线为原始子波。

观察发现，在子波估计长度参数与子波真实长度不同时得到的地震子波与原始子波在波形上仍能保持总体一致，子波估计长度参数对子波提取结果的影响不大。

在本实施例中，利用地震记录的高阶谱进行混合相位地震子波估计，除了能解决地震子波非最小相位问题外，这种方法的优势还在于它具有较强的抗高斯噪声的能力，能够从高斯背景的噪声中分离出非高斯信号。

但是，由于实际地震数据是有限长度的，并且噪声也不完全满足高斯分布，因此噪声的四阶累积量在很多情况下并不等于零。为了验证这种方法的抗噪性能，通过在人工合成记录中加入不同信噪比的高斯噪声，分析了噪声对子波提取的影响。

在不同的噪声水平下得到的结果如图16、图17和图18所示，三幅图中合成地震记录的信噪比分别为30dB、20dB和10dB，噪声强度依次增加。其中，点线为估计子波，实线为原始子波。在子波提取过程中，子波估计参数的长度均等于真实子波的长度。

分析发现，在噪声强度较低时，子波估计结果比较准确，估计子波与原始子波在波形上保持了较好的一致性；但随着噪声水平的增加，估计结果变差。在噪声水平很高时，估计子波的波形与原始子波相比发生了很大的变化，估计结果将不可信。地震资料的信噪比越高，子波估计结果将越准确。

在本实施例中，图19是从某地区实际地震数据中截取的一部分的示意图，包含80道地震数据，时间采样间隔为1ms，记录长度为1s。图20是根据其中的5道地震数据利用三谱法提取的子波，其中图20中(a)(b)(c)(d)(e)的对应的道号分别为5，15，25，35，45，图20中(f)为提取的5个混合相位子波的平均子波。通过观察图20中根据各道提取的子波波形可以发现，各道提取的子波具有很好的一致性和稳定性。

图21是平均子波的频谱及其Z变换后根的分布示意图。如图21所示，子波的根在单位圆内外均有分布，混合相位特征非常明显，从而验证了高阶谱方法对混合相位地震子波提取的有效性。

由上述实施例可知，本发明实施例将数字信号处理领域中的高阶谱估计优化算法引入地震勘探领域中，提出了三谱法混合相位地震子波估计技术。利用地震记录的高阶谱进行混合相位地震子波估计，除了能解决地震子波非最小相位问题外，这种方法的优势还在于它具有较强的抗高斯噪声的能力，能够从高斯背景的噪声中分离出非高斯信号。

但是，由于实际地震数据是有限长度的，并且噪声也不完全满足高斯分布，因此噪声的四阶累积量在很多情况下并不等于零，但是通过加窗函数优化，能够降低噪声对子波估计结果的影响，提高子波估计的精度，模型地震数据和实际资料均得到了较好的效果，从而验证这种技术具有较好的应用前景。

本发明实施例还提供一种混合相位地震子波的提取装置，如图22所示，所述装置包括：计算单元2201、处理单元2202和提取单元2203；其中，

计算单元2201计算地震数据的四阶累积量；

处理单元2202通过三维滞后窗函数对地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；

提取单元2203通过地震子波的三谱提取混合相位地震子波。

进一步地，如图23所示，提取单元2203具体可包括：递推单元2301、获取单元2302和生成单元2303；其中，

递推单元2301通过Pan-Nikias公式由地震子波的三谱的相位谱递推出地震子波的相位谱；

获取单元2302通过二阶统计方法获取地震子波的振幅谱；

生成单元2303根据地震子波的相位谱和振幅谱生成混合相位地震子波。

优选地，三维滞后窗函数通过标准的一维滞后窗函数构造。所述一维滞后窗函数包括：Daniel窗、Hamming窗、Parzen窗或者Sasaki窗。

优选地，所述三维滞后窗函数为：

其中d(τ)为Parzen窗。

本实施例的装置的各组成部分分别用于实现前述实施例的方法的各步骤，由于在方法实施例中，已经对各步骤进行了详细说明，在此不再赘述。

由上述实施例可知，通过三维滞后窗函数对地震记录的三谱进行平滑，从而提高了高阶谱估计的精度和子波估计的精度。

本领域普通技术人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

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

Claims (2)

1.一种提高地震数据处理精度的方法，其特征在于，所述方法包括：

获取地震数据并计算地震数据的四阶累积量；

通过三维滞后窗函数对所述地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；所述三维滞后窗函数通过标准的一维滞后窗函数构造；所述三维滞后窗函数为：

其中，所述一维滞后窗函数d(τ)为Parzen窗；

通过所述地震子波的三谱提取混合相位地震子波，以进行地震正演模拟、或常规反演、或地震反褶积；其中所述通过所述地震子波的三谱提取所述混合相位地震子波，具体包括：通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱；通过二阶统计方法获取地震子波的振幅谱；根据所述地震子波的相位谱和振幅谱生成所述混合相位地震子波；所述通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱包括：确定所述Pan-Nikias公式：

其中，是所述地震子波的相位谱，n＝2,3,...,N，s(n)是所述地震子波的三谱的相位谱；

构造线性方程组：

$\begin{array}{c}s\left(2\right)=2\mathrm{\φ}\left(1\right)-\mathrm{\φ}\left(2\right)\\ s\left(3\right)=3\mathrm{\φ}\left(1\right)+2\mathrm{\φ}\left(2\right)-\frac{14}{6}\mathrm{\φ}\left(3\right)\\ .\\ .\\ .\\ s\left(N\right)=\mathrm{N\φ}\left(1\right)+(N-1)\mathrm{\φ}\left(2\right)+...+2\mathrm{\φ}(N-1)+\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)\end{array},;$

结合所述Pan-Nikias公式和所述线性方程组得到矩阵方程：Bφ＝s，

其中，B是(N-1)×(N-1)维方阵， $s=[s\left(1\right),...,s(N-1),s\left(N\right)-\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)]$ 是所述地震子波的三谱的相位谱的矩阵，φ＝[φ(1),φ(2),...,φ(N-1)]^T是所述地震子波的相位谱的矩阵；

计算地震子波的相位谱：φ＝B^-1s。

2.一种提高地震数据处理精度的装置，其特征在于，所述装置包括：

计算单元，获取地震数据并计算地震数据的四阶累积量；

处理单元，通过三维滞后窗函数对所述地震数据的四阶累积量进行平滑处理，并估计地震子波的三谱；所述三维滞后窗函数通过标准的一维滞后窗函数构造；所述三维滞后窗函数为：

其中，所述一维滞后窗函数d(τ)为Parzen窗；

提取单元，通过所述地震子波的三谱提取混合相位地震子波，以进行地震正演模拟、或常规反演、或地震反褶积；其中，所述提取单元具体包括：递推单元，通过Pan-Nikias公式由所述地震子波的三谱的相位谱递推出地震子波的相位谱；获取单元，通过二阶统计方法获取地震子波的振幅谱；生成单元，根据所述地震子波的相位谱和振幅谱生成所述混合相位地震子波；所述递推单元具体用于，确定所述Pan-Nikias公式：

其中，是所述地震子波的相位谱，n＝2,3,...,N，s(n)是所述地震子波的三谱的相位谱；

构造线性方程组：

$\begin{array}{c}s\left(2\right)=2\mathrm{\φ}\left(1\right)-\mathrm{\φ}\left(2\right)\\ s\left(3\right)=3\mathrm{\φ}\left(1\right)+2\mathrm{\φ}\left(2\right)-\frac{14}{6}\mathrm{\φ}\left(3\right)\\ .\\ .\\ .\\ s\left(N\right)=\mathrm{N\φ}\left(1\right)+(N-1)\mathrm{\φ}\left(2\right)+...+2\mathrm{\φ}(N-1)+\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)\end{array},;$

结合所述Pan-Nikias公式和所述线性方程组得到矩阵方程：Bφ＝s，

其中，B是(N-1)×(N-1)维方阵， $s=[s\left(1\right),...,s(N-1),s\left(N\right)-\frac{6-(N+1)(N+2)}{6}\mathrm{\φ}\left(N\right)]$ 是所述地震子波的三谱的相位谱的矩阵，φ＝[φ(1),φ(2),...,φ(N-1)]^T是所述地震子波的相位谱的矩阵；

计算地震子波的相位谱：φ＝B^-1s。