一种基于二维非均匀曲波变换的高精度地震数据重建方法与流程

本发明涉及的是空间非均匀采样下不规则缺失道的地震数据重建方法，具体是一种基于二维非均匀曲波变换的高精度地震数据重建方法。

技术背景

在野外数据采集过程中，由于采集设备、野外地形条件以及经济成本限制等原因，地震数据通常沿空间方向常进行不规则欠采样，从而导致采集到的地震数据不规则、不完整，出现空间假频，影响到后续其它处理方法的效果，降低了地震勘探的分辨率。为了克服这种情况，野外不规则缺失地震数据必须进行叠前数据重建，使得缺失道得到有效地恢复，然而目前大多数数据重建方法的前提条件是空间均匀采样下的不规则缺失地震道重建，而对于空间非均匀采样下的地震道缺失重建效果较差或无能为力。但由于野外复杂地形条件的限制或者海上电缆的水平偏移，很多情况下野外地震数据常进行空间非均匀采样，如不加处理则会引起覆盖次数的变化(地下不均匀照明)，在叠加成像时会形成扭曲的成像振幅(采集脚印)，加重了空间假频现象的出现，影响后续成像处理。

为了解决空间非均匀采样下地震道不规则缺失问题，地球物理领域常规处理方法为共面元叠加，从而将非均匀采样数据归位到均匀采样数据中来，满足后续其他处理方法的要求。然而共面元叠加处理方法忽略了每个面元内各道共中心点的真实位置，改变了部分地震道的振幅和相位，从而导致部分地震道位置出现严重偏差，降低了地震勘探资料的分辨率。而如果采用使用波动方程对非均匀采样数据进行重建时，则需要地下结构的先验信息，计算量非常巨大，对采样率要求也较高，从而也不能较好的解决该问题。因此，许多学者采用基于数学变换的重建方法对该问题进行处理，Duijndam(1999)提出基于傅立叶变换的二维非均匀采样重建技术，Hindriks(2000)利用傅立叶变换的方法实现三维非均匀采样重建技术。但是，Duijndam等人的傅立叶重建方法还是有其局限性，重建结果受最低速度和空间带宽的影响很大，随着采样间隔的逐渐增大，重建结果会逐渐变差。进而，Zwartjes(2007)等在此基础上进行了改进，达到压制假频和混淆能量的目的，Jin(2010)提出基于阻尼最小范数傅立叶反演下的五维地震数据重建，该方法能够重建空间非均匀采样下的不规则缺失地震数据，但不具有反假频能力。可见，以上方法都是采用傅里叶变换作为稀疏基，并没有采用曲波变换，而傅里叶变换作为全局变换，只适合同相轴近似线性或者平稳变化的地震信号，而曲波变换能够表征信号的局部细节特征，众多的研究结果也证明，基于曲波变换的数据重建方法效果显著。

尽管如此，以往基于二维曲波变换的重建方法前提条件是空间均匀采样，因为常规的曲波变换在计算过程中首先要应用到傅里叶变换，而傅里叶变换的前提条件是空间均匀采样，从而导致以往二维曲波变换只能处理空间均匀采样下的地震道缺失重建，而对于空间非均匀采样下的地震数据则不能直接重建，限制了该方法的进一步应用。

技术实现要素：

本发明的目的是为了能够高精度重建野外非均匀采样下不规则缺失道，并大幅度提高重建信号的保真度和信噪比，保护微弱的有效波信号，从而使反射波同相轴更加连续，而提出了一种基于二维非均匀曲波变换的高精度地震数据重建方法。

本发明提出了一种基于二维非均匀曲波变换的高精度地震数据重建方法，首先针对常规二维曲波变换方法难以对空间非均匀采样下地震缺失道进行重建问题，在多尺度多方向二维曲波正变换的基础上，采用非均匀傅里叶变换，建立曲波系数与空间非均匀采样下地震缺失数据之间的正反变换算子，在最小L1范数约束下，使用谱投影梯度法进行反演计算得到均匀曲波系数，并且通过对该曲波系数的标准二维曲波反变换处理，从而最终实现一种基于二维非均匀曲波变换的高精度地震数据重建方法。

进一步，基于常规的二维曲波正变换的四个步骤：(1)对地震数据应用二维傅立叶变换，得到频率波数域系数；(2)在频率波数域形成角度楔形；(3)将每一个楔形围绕到原点进行重新装配；(4)对每一个楔形应用二维傅立叶反变换，得到曲波系数；定义反变换算子A：

$A\stackrel{def}{=}FT$

这里F代表二维傅立叶反变换，将频率波数域转换到时间空间域中，T表示曲波平铺算子，即将曲波系数变换到频率波数域的过程，定义正算子为A^H是：

A^H＝T^HF^H

该式F^H实现了常规的二维曲波正变换第(1)步，T^H则实现了第(2)～第(4)步；

由于二维傅立叶变换参与了曲波变换之中，导致常规的二维曲波变换不能处理非均匀采样数据，然而，用F_x来代表沿着空间轴的一维傅立叶反变换算子，F_t表示沿着时间轴的一维傅立叶反变换算子，因此，

$F=F_x\⊗F_t$

时间方向是均匀理想采样，不需要重建；因此，用非均匀傅立叶变换N_x代替一维傅立叶反变换算子F_x，因此新的非均匀曲波反变换算子可以定义为：

$B\stackrel{def}{=}(N_x\⊗F_t)T$

该算子将曲波系数与非均匀采样下不规则地震道建立相应的联系；

当野外地震道缺失严重，输入的不规则地震道一般少于傅里叶系数，属于欠定的情况，为此定义非均匀曲波正变换算子为：

从该式得知，正变换算子非线性映射非均匀采样地震数据y到曲波系数向量x是属于BP问题，由于该方程是欠定的，有无限多向量x满足y＝Bx，x是满足最小L1范数中的一个；方程两边相等约束保证了正反变换算子对是能量无损伤，满足称为二维非均匀曲波正反变换，而对于上述方程的求解，通过谱投影梯度法来求解此L1范数最小化问题；

在求解上式方程得到曲波系数后，重建后的地震波场f可以通过下式得到：

其中，A是标准的二维曲波反变换算子，是二维非均匀曲波变换正算子，通过该式重建出空间均匀采样下规则而完整的地震数据，实现一种基于二维非均匀曲波变换的高精度地震数据重建方法。

进一步，所述二维曲波正变换的定义为：

$c(j,l,k)=\<d,{\mathrm{\ψ}}_{j,l,k}\>={\∫}_{R^2}d\left(x\right)\stackrel{\‾}{{\mathrm{\ψ}}_{j,l,k}\left(x\right)}dx$

式中：ψ_j,l,k表示曲波函数，c(j,l,k)为曲波系数，j,l,k分别表示尺度，方向和位置参数，d(x)表示地震数据，其反变换为：

$d\left(x\right)=\<d,{\mathrm{\ψ}}_{j,l,k}\>{\mathrm{\φ}}_{j,l,k}={\∫}_{R^2}c(j,l,k){\mathrm{\ψ}}_{j,l,k}\left(x\right)dx.$

进一步，所述的非均匀傅里叶变换的定义为：

长度为N信号序列x[n]的一维非均匀傅里叶变换可表示如下：

$X\left(z_k\right)=X\left(z\right){|}_{z=z_k}=\underset{n=0}{\overset{N-1}{\Σ}}x\[n\]z_k^{-n},k=0,\mathrm{1...},N-1$

在这里X(z)是x[n]的Z变换，并且{z_i}_{i＝0,1…N-1}是z平面中的任意不同的点；上式如用矩阵表达上式，可得：

X＝Dx

其中，并且

如果Z变换是非奇异的，得到唯一的非均匀傅里叶反变换为：

x＝D^-1X

因此，给定了X和D，采用高斯消去法求解出x。

进一步，所述的谱投影梯度法算法如下：

首先定义f(x)＝y-Bx，P(c)是向量c在一个闭凸集Ω上的正交投影，并且可以定义得到：

P_τ＝arg min||c-x||₂||x||≤τ

其中τ代表的是L1范数球的半径，此时有Ω＝{x|||x||₁≤τ}，将f(x)对x求导，得到其梯度方向：

g＝▽f(x)＝-B^T(y-Bx)，

由于残差则在第n次迭代中有gⁿ＝-Ω^Tr^n-1，给定P(c)和gⁿ，则：

$d^n=P_{\mathrm{\τ}}(\stackrel{\‾}{x}-a^{n-1}{g}^{n-1})-{\stackrel{\‾}{x}}^{n-1}$

同时其中d为每次更新方向，aⁿ是每次迭代时的步长，最终可以计算出曲波系数。

上述L1范数，采用通用表达式，表达式如下：

$\left|\right|x|{|}_1\stackrel{def}{=}{\Σ}_{i=1}^N\left|x\[i\]\right|$

其中x[i]表示向量x的第i个元素。

本发明在二维曲波变换过程中引入非均匀傅里叶变换，建立均匀曲波系数与空间非均匀采样下地震缺失道数据之间的正反变换算子，然后使用谱投影梯度算法进行反演计算得到均匀曲波系数，最后再进行二维均匀曲波反变换，最终实现一种基于二维非均匀曲波变换的高精度地震数据重建方法。

本发明创造性主要体现在：

1、相对于共面元叠加和波动方程方法，本发明技术可以在不改变振幅和相位的前提下直接重建出缺失道地震数据，具有较好的保真度，并且不需要地下结构先验信息，计算工作量较少。

2、相对于傅里叶变换，本发明技术更能够表征具有曲线状特征的地震信号局部细节特征，重建精度高，而且具有反假频能力。

3、相对常规二维曲波变换方法，本发明技术不仅能够处理空间均匀采样下的地震道缺失重建，而且能够处理空间非均匀采样下的地震道缺失重建。大幅度提高了重建信号的信噪比，保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰，

4、本发明技术对于其它非均匀信号处理领域也具有重要的借鉴意义。

综述，本发明克服了传统二维曲波变换方法不能重建非均匀采样下地震缺失道的缺点，并且该方法在重建过程能量中无损失，具有反假频功能，大幅度提高了重建信号的保真度和信噪比，保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰，同时对于其它非均匀信号处理领域也具有重要的借鉴意义。

附图说明

图1是本发明实施例中同时数据重建与噪声压制流程图。

图2是原始理论地震数据图。

图3是非均匀地震数据图。

图4是非均匀曲波变换规则化重建结果图。

图5是50％非均匀采样下不规则缺失地震数据图。

图6是非均匀采样下不规则缺失重建结果图。

图7是50％非均匀采样下规则缺失地震数据图。

图8是非均匀采样下重建前数据频谱分析图。

图9是非均匀采样下规则缺失重建结果图。

图10是非均匀采样下重建后数据频谱分析图。

具体实施方式

以下实施案例用于说明本发明，但不用来限制本发明的范围。

实施例1

实现该方法的步骤主要包括，常规二维曲波正反变换，非均匀傅里叶变换，建立非均匀曲波正反变换算子，谱梯度投影法求解等。具体步骤如下：

步骤1：二维曲波正反变换。为了得到高精度地震数据重建结果，需要在二维曲波变换的基础上进行数据重建，二维曲波正变换的定义为：

$c(j,l,k)=\<d,{\mathrm{\ψ}}_{j,l,k}\>={\∫}_{R^2}d\left(x\right)\stackrel{\‾}{{\mathrm{\ψ}}_{j,l,k}\left(x\right)}dx$

式中：ψ_j,l,k表示曲波函数，c(j,l,k)为曲波系数，j,l,k分别表示尺度，方向和位置参数，d(x)表示地震数据，其反变换为：

$d\left(x\right)=\<d,{\mathrm{\ψ}}_{j,l,k}\>{\mathrm{\φ}}_{j,l,k}={\∫}_{R^2}c(j,l,k){\mathrm{\ψ}}_{j,l,k}\left(x\right)dx$

步骤2：非均匀傅里叶变换。针对传统基于二维曲波变换的重建方法只能重建均匀采样下地震道缺失的现象，引入非均匀傅里叶变换，从而可以解决非均匀采样下地震道缺失重建问题，其定义表达式为：

长度为N信号序列x[n]的一维非均匀傅里叶变换可表示如下

$X\left(z_k\right)=X\left(z\right){|}_{z=z_k}=\underset{n=0}{\overset{N-1}{\Σ}}x\[n\]z_k^{-n},k=0,\mathrm{1...},N-1$

在这里X(z)是x[n]的Z变换，并且{z_i}_{i＝0,1…N-1}是z平面中的任意不同的点。上式也可以用矩阵表达上式，可得

X＝Dx

其中，并且如果Z变换是非奇异的，可以得到唯一的非均匀傅里叶反变换为

x＝D^-1X

因此给定了X和D，可以采用高斯消去法求解出x。

步骤3：建立非均匀曲波正反变换算子。常规的二维曲波变换主要进行了四个步骤，(1)对地震数据应用二维傅立叶变换，得到频率波数域系数；(2)在频率波数域形成角度楔形；(3)将每一个楔形围绕到原点进行重新装配；(4)对每一个楔形应用二维傅立叶反变换，得到曲波系数。为此，定义反变换算子A：

$A\stackrel{def}{=}FT$

这里F代表二维傅立叶反变换，将频率波数域转换到时间空间域中，T表示曲波平铺算子，即将曲波系数变换到频率波数域的过程，定义正算子为A^H是：

A^H＝T^HF^H

该式F^H实现了曲波变换第(1)步，T^H则实现了第(2)～第(4)步。

由于二维傅立叶变换参与了曲波变换之中，导致常规的二维曲波变换不能处理非均匀采样数据，然而，可以用F_x来代表沿着空间轴的一维傅立叶反变换算子，F_t表示沿着时间轴的一维傅立叶反变换算子，因此，

$F=F_x\⊗F_t$

时间方向是均匀理想采样，不需要重建。因此，可以用非均匀傅立叶变换N_x代替一维傅立叶反变换算子F_x，因此新的非均匀曲波反变换算子可以定义为：

$B\stackrel{def}{=}(N_x\⊗F_t)T$

该算子可以将曲波系数与非均匀采样下不规则地震道建立相应的联系。

当野外地震道缺失严重，输入的不规则地震道一般少于傅里叶系数，属于欠定的情况，为此定义非均匀曲波正变换算子为：

从该式得知，正变换算子非线性映射非均匀采样地震数据y到曲波系数向量x是属于BP问题，由于该方程是欠定的，有无限多向量x满足y＝Bx，x是满足最小L1范数中的一个。方程两边相等约束保证了正反变换算子对是能量无损伤，满足称为二维非均匀曲波正反变换算子，而对于上述方程的求解，本发明选用谱投影梯度法来求解此L1范数最小化问题。

在求解上式方程得到曲波系数后，重建后的地震波场f可以通过下式得到

其中，A是标准的二维曲波反变换算子，是二维非均匀曲波变换正算子，通过该式可以重建出空间均匀采样下规则而完整的地震数据，实现一种基于二维非均匀曲波变换的高精度地震数据重建方法。

步骤4：谱梯度投影算法。首先定义f(x)＝y-Bx，P(c)是向量c在一个闭凸集Ω上的正交投影，并且可以定义得到

P_τ＝arg min||c-x||₂||x||≤τ

其中τ代表的是L1范数球的半径，此时有Ω＝{x|||x||₁≤τ}，将f(x)对x求导，得到其梯度方向：

g＝▽f(x)＝-B^T(y-Bx)，

由于残差则在第n次迭代中有gⁿ＝-Ω^Tr^n-1，给定P(c)和gⁿ，则

$d^n=P_{\mathrm{\τ}}(\stackrel{\‾}{x}-a^{n-1}{g}^{n-1})-{\stackrel{\‾}{x}}^{n-1}$

同时其中d为每次更新方向，aⁿ是每次迭代时的步长，最终可以计算出曲波系数。

实现该方法具体操作为：

为了详细阐述本发明所提的一种基于二维非均匀曲波变换的高精度地震数据重建方法的处理效果，定义信噪比SNR＝20log₁₀||f₀||₂/||f-f₀||₂来对比处理数据后的质量，单位为dB，其中f₀表示原始模型数据，f表示重建结果，信噪比越高，代表重建结果与模型数据越接近，效果越理想。

图2为采用40Hz雷克子波合成的256道二维理论地震数据图(图2表示原始理论地震数据图)，该记录总共有4层地震反射波，每一层反射波能量有所差异，采样间隔为1ms，道距为4m，每道1024个采样点，横坐标表示距离，纵坐标表示时间。首先为了验证本发明所提方法的重建效果，对理论地震数据空间域进行均匀傅里叶变换，然后再进行空间非均匀傅里叶反变换，得到新的空间非均匀采样下的256道地震数据，如图3所示(图3表示非均匀地震数据图)，此时信噪比为14.75dB，尽管名义上的道距还是4米，但每道地震数据道距不均匀，其道距范围为1.03m～7.16m，显然，如果将非均匀采样下地震数据在均匀采样网格上进行显示，连续的地震波场则会被破坏，而常规重建方法的前提条件是均匀道距下不规则缺失的地震记录。为此采用本发明方法进行规则化重建，规则化重建后的道距为4m，结果如图4所示(图4表示非均匀曲波变换规则化重建结果图)，规则化重建后的信噪比为46.17dB，地震波场连续性显著提高，重建后的地震记录与原始记录非常接近，误差几乎忽略不计，重建后信噪比非常高，与原始地震记录相比较几乎没有视觉上的差异，从而说明本发明方法能够反映出地震波场的局部细节特征，因此重建方法精度高，保真度较好。

为了检验本发明方法在非均匀采样下不规则缺失道重建效果，同样采用空间非均匀傅里叶变换，得到128道新的地震数据，如图5所示(图5表示50％非均匀采样下不规则缺失地震数据图)，相当于平均道距为8m左右，但是由于各检波器道距不均匀，此时128道地震记录道距范围为1.25m～18.26m，可以看出地震数据道距极不均匀，同相轴不连续，不能直接应用于后续资料的处理，必须采用非均匀地震数据重建方法进行处理。为此利用本发明方法进行重建，重建后的道距为4米，重建结果如图6所示(图6表示非均匀采样下不规则缺失重建结果图),重建后的信噪比为20.01dB，可以看出尽管不均匀采样下地震数据缺失50％地震道，但是重建效果精度仍然较高，重建后同相轴更连续，能量损失较少。

为了进一步检验本发明方法的反假频能力，对原始非均匀采样地震记录进行50％规则欠采样，如图7所示(图7表示50％非均匀采样下规则缺失地震数据图)，图8为其二维频谱图(图8表示非均匀采样下重建前数据频谱分析图)，从中可以看出规则欠采样带来较为严重的假频成分，与信号的真实频谱存在部分重叠，然后采用本发明方法进行重建，重建后道距为4m，重建结果如图9所示(图9表示非均匀采样下规则缺失重建结果图)，重建后信噪比为16.02dB，重建效果较好，缺失的地震道得到了有效的恢复，图10为其二维频谱分析图(图10表示非均匀采样下重建后数据频谱分析图)，可以看出假频几乎消失，其二维频谱与地震记录的真实频谱更加接近，表明本发明方法具有较强的反假频能力，能够进行复杂地区非均匀采样下的不规则和规则缺失地震数据重建。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式来实现本发明的实用功能。