一种能够同时进行五维地震数据重建和噪声压制的方法

摘要

本发明提出一种能够同时进行五维地震数据重建和噪声压制的方法，首先针对低维数据重建精度的不足，采用傅里叶变换作为稀疏基，通过利用五维地震数据信息进行数据重建，在此过程中引入凸集投影算法，提出指数平方根衰减规律的阈值参数，采用硬阈值算子对每个时间切片单独进行重建，然后在传统凸集投影算法中引入加权因子，使得在的重建过程中减少噪声对重建结果的影响，最终实现了一种能够同时进行五维地震数据重建和噪声压制的方法。本发明提出的同时进行五维地震数据重建和噪声压制的方法，克服了传统二维或者三维地震数据重建精度不足，以及传统方法不能同时进行数据重建和噪声压制的缺点，大幅度提高了重建精度，降低了对计算机内存的要求。

说明书

技术领域

本发明涉及的是地震数据重建与噪声压制的方法，具体是一种基于傅里叶变换下能够同时进行五维地震数据重建和噪声压制的方法。

技术背景

在野外数据采集过程中，由于采集设备和经济成本的限制，地震数据沿空间方向上通常进行不规则采样，从而导致采集到的地震数据不规则、不完整，出现空间假频，这种稀疏分布的地震数据难以满足后续处理的要求，例如速度分析、自由表面多次波消除、偏移归位等。在这种情况下，需要发展较好的叠前数据重建方法。

目前，广泛使用的数据重建方法是基于某种数学变换，该方法不需要地下结构的先验信息，能够处理规则缺失和不规则缺失地震道的重建。同时，由于野外采集到的地震数据常常受到高频随机噪声干扰的影响，降低了地震记录的信噪比，影响到数据重建方法的效果，造成缺失道不能充分重建或者重建后精度较低，严重影响到信号的准确成像和波场归位，这就需要在数据重建过程中对随机干扰噪声进行适当处理。然而目前大部分有效的噪声压制方法与数据重建方法都是单独分开进行处理，缺少一种能够同时进行地震数据重建和噪声压制的方法。

对于数据重建算法，目前已经发展了各种有效的数据重建算法，其中凸集投影算法由于参数设置简单且效果显著而被广泛应用，其主要思想就是在每次迭代过程中，采用傅立叶变换将时间-空间域数据变换到频率－波数域，然后设置一个阈值保留较大振幅的有用信号，再对保留后的有用信号采用傅立叶反变换，最后将原始不需要重建的地震道置换到傅立叶反变换后的数据中去，通过多次迭代，从而将缺失地震道重建出来。该方法最早由Bregman(1965)提出，而后广泛应用于图像处理。Abma(2006)首次将POCS算法应用于不规则地震数据的重建，取得了较好的应用效果，但采用的是线性阈值参数，收敛速度慢，运算时间长。Gao(2010，2013)提出一种指数阈值参数和数据驱动阈值参数，实现了三维地震数据重建去和噪声压制，提高了该方法的收敛速度，而后Liu(2011)，Yang(2012)，Zhang(2013)分别采用曲波变换方法，在曲波域实现了基于POCS算法的二维和三维地震数据重建。但是地震数据采集已经发展到五维(x_s,y_s,x_r,y_r,t)，其中(x_s,y_s)和(x_r,y_r)分别代表偏移距和共中心点坐标，时间t为双程旅行时，因此地震数据重建和噪声压制方法也应该发展到五维，在有效的去除噪声同时高精度地重建出缺失道地震数据。由于时间方向不需要重建，本发明拟通过利用五维地震数据中的四维空间数据信息进行高精度地震数据重建，同时引入加权因子策略，最终实现一种能够同时进行五维地震数据重建和噪声压制的方法。

发明内容

本发明的目的是为了能够对野外缺失道地震数据重建，并且同时能够有效地压制噪声干扰，提高信噪比，而提出了一种能够同时进行五维地震数据重建和噪声压制的方法。

本发明提出了一种能够同时进行五维地震数据重建和噪声压制的方法，首先针对低维数据重建精度的不足，采用傅里叶变换作为稀疏基，通过利用五维地震数据信息进行数据重建，在此过程中引入凸集投影算法，提出指数平方根衰减规律的阈值参数，采用硬阈值算子对每个时间切片单独进行重建，然后在传统凸集投影算法中引入加权因子，使得在的重建过程中减少噪声对重建结果的影响，最终实现了一种能够同时进行五维地震数据重建和噪声压制的方法。

一种能够同时进行五维地震数据重建和噪声压制的方法，地震数据的重建问题描述为由一组不完整数据通过线性算子的作用恢复出完整的数据，假设如下线性正演模型：

y_obs＝Md

这里y_obs∈Rⁿ代表采集的地震数据；d∈R^N，且N≥n，表示待重建的无假频数据；M∈R^n×N表示对角矩阵，其元素1和0分别表示已知地震道和未知地震道。

假设系数x是d在傅里叶域F中的稀疏表示，则上述方程为：

y_obs＝Md＝MF^Hx＝Ax

其中A＝MF^H，常称为测量矩阵，上式同样可以写成

>

在这个表达式中，代表估计值，L₁范数定义为x[i]是向量x中第i个元素，通过求解上述方程，原始无假频的数据就重建出来。

在五维地震数据重建和噪声压制过程，采用凸集投影算法，具体重建和噪声压制步骤如下：

步骤1：在时间域中输入具有缺失道的五维含噪地震数据，然后采用傅里叶变换对时间域五维地震数据进行稀疏变换，得到傅里叶系数，根据傅里叶系数的大小选择合适的阈值参数公式；

步骤2：在傅里叶域中，采用硬阈值算子进行处理，也即大于阈值λ_i的曲波系数保留，而其它的傅里叶系数置零；

步骤3：将阈值化后的傅里叶系数做反变换得到时间域地震数据，选择合适的加权因子，减少迭代过程中的噪声影响，然后再将未缺失含噪地震道填充到反变换后的地震数据中去；

步骤4：最后将得到的地震数据代入步骤(1)，重新进行迭代，直到运行N次结束，即可得到最终的重建和噪声压制结果。

进一步，所述的凸集投影算法，其迭代表达式为：

d_k(t,x_s,y_s,x_r,y_r)＝y_obs(t,x_s,y_s,x_r,y_r)+[I-M(x_s,y_s,x_r,y_r)]F_t^-1T_k(t,x_s,y_s,x_r,y_r)

×F_td_k-1(t,x_s,y_s,x_r,y_r),k＝1,2,…,N.

其中，d_k表示第k次迭代得到的(t,x_s,y_s,x_r,y_r)域重建数据，d₀表示原始采集到数据y_obs(t,x_s,y_s,x_r,y_r)，满足d₀＝y_obs，F_t和F_t^-1表示关于时间变量t的四维正反傅立叶变换，T_k(t,x_s,y_s,x_r,y_r)表示硬阈值算子。

进一步，所述硬阈值算子，其表达式如下：

$>T_k(t,x_s,y_s,x_r,y_r)=\left(\begin{array}{cc}1,& \left|S_{k-1}(t,x_s,y_s,x_r,y_r)\right|\ge {\lambda }_k\\ 0,& \left|S_{k-1}\left(\right(t,x_s,y_s,x_r,y_r)\right|<{\lambda }_k\end{array}\right),{\lambda }_k\in \lambda$>

S_k-1表示第k-1次迭代得到的傅里叶系数，满足S_k-1＝F_td_k-1(t,x_s,y_s,x_r,y_r)，λ表示N维阈值集合，λ＝{λ₁,λ₂,…,λ_N}，且满足|Cd|_max＝λ₁>λ₂>…>λ_N≥ε，N表示最大迭代次数，其中式中ε为接近零的小值，与五维地震数据中噪声的能量有关，不同数据ε值有所差别。

进一步，所述在传统凸集投影算法中引入加权因子，其表达式如下：

d_k(t,x_s,y_s,x_r,y_r)＝a×y_obs(t,x_s,y_s,x_r,y_r)+[I-a×M(x_s,y_s,x_r,y_r)]F_t^-1T_k(t,x_s,y_s,x_r,y_r)×F_td_k-1(t,x_s,y_s,x_r,y_r),k＝1,2,…,N.

其中a为加权因子，其范围在0<a≤1，如果a＝1，则方程与常规的凸集投影算法一样，此时原始噪声数据被带入到重建后的地震数据中，只能进行数据重建，不能进行噪声压制。不同的数据加权因子的选择不一样，与噪声能量的强度有关。

进一步，所述阈值参数公式，其表达式如下：

$>{\lambda }_k=Max·e^{\left(\ln \right(ϵ)-ln(Max\left)\right]\sqrt{\frac{k-1}{N-1}}}$>

其中Max为傅里叶变换系数绝对值的最大值；N为总的迭代次数，ε为接近零的小值，与五维地震数据中噪声的能量有关。

本发明的优点：本发明采用傅里叶变换作为稀疏基，通过利用五维地震数据信息进行数据重建，在此过程中引入凸集投影算法，提出指数平方根衰减规律的阈值参数，采用硬阈值算子对每个时间切片单独进行重建，然后在传统凸集投影算法中引入加权因子，实现了一种能够同时进行五维地震数据重建和噪声压制的方法，从而克服了传统二维或者三维地震数据重建精度不足，以及传统方法不能同时进行数据重建和噪声压制的缺点，降低了对计算机内存的要求，大幅度提高了重建精度，保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰。

附图说明

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

图2是原始模拟地震数据图。

图3是原始加噪数据图。

图4是20％随机欠采样图。

图5是五维地震数据重建结果图。

图6是阈值法五维地震数据噪声压制结果图。

图7是同时五维地震数据重建与噪声压制结果图。

图8是三维地震数据同时重建与噪声压制结果图。

具体实施方式

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

实施例1

实现该方法的步骤主要包括，重建方程的构建，凸集投影重建算法，采用加权因子进行噪声压制，阈值参数公式等。具体步骤如下：

步骤1：重建方程的构建。地震数据的重建问题描述为由一组不完整数据通过线性算子的作用恢复出完整的数据，假设如下线性正演模型：

y_obs＝Md

这里y_obs∈Rⁿ代表采集的地震数据；d∈R^N，且N≥n，表示待重建的无假频数据，M∈R^n×N表示对角矩阵，其元素1和0分别表示已知地震道和未知地震道。

假设系数x是d在傅里叶域F中的稀疏表示，则上述方程为：

y_obs＝Md＝MF^Hx＝Ax

其中A＝MF^H，常称为测量矩阵，上式同样可以写成

>

在这个表达式中，代表估计值，L₁范数定义为x[i]是向量x中

第i个元素，通过求解上述方程，原始无假频的数据就重建出来。

步骤2：凸集投影重建算法。针对传统重建算法只能利用二维或者三维地震信息的缺点，提出五维地震数据重建方法，并采用硬阈值算子和新的按指数规律衰减的阈值参数公式，其算法迭代表达式为：

d_k(t,x_s,y_s,x_r,y_r)＝y_obs(t,x_s,y_s,x_r,y_r)+[I-M(x_s,y_s,x_r,y_r)]F_t^-1T_k(t,x_s,y_s,x_r,y_r)

×F_td_k-1(t,x_s,y_s,x_r,y_r),k＝1,2,…,N.

其中，d_k表示第k次迭代得到的(t,x_s,y_s,x_r,y_r)域重建数据，d₀表示原始采集到数据y_obs(t,x_s,y_s,x_r,y_r)，满足d₀＝y_obs，F_t和F_t^-1表示关于时间变量t的四维正反傅立叶变换，T_k(t,x_s,y_s,x_r,y_r)表示硬阈值算子，其元素满足：

$>T_k(t,x_s,y_s,x_r,y_r)=\left(\begin{array}{cc}1,& \left|S_{k-1}(t,x_s,y_s,x_r,y_r)\right|\ge {\lambda }_k\\ 0,& \left|S_{k-1}\left(\right(t,x_s,y_s,x_r,y_r)\right|<{\lambda }_k\end{array}\right),{\lambda }_k\in \lambda$>

S_k-1表示第k-1次迭代得到的傅里叶系数，满足S_k-1＝F_td_k-1(t,x_s,y_s,x_r,y_r)，λ表示N维阈值集合，λ＝{λ₁,λ₂,…,λ_N}，且满足|Cd|_max＝λ₁>λ₂>…>λ_N≥ε，N表示最大迭代次数，其中式中ε为接近零的小值，与五维地震数据中噪声的能量有关，不同数据ε值有所差别。

步骤3：采用加权因子进行噪声压制。针对传统重建算法不能同时进行数据重建和噪声压制的缺点，对被重新置入的原始数据加入一个加权因子，其算法迭代表达式为：

d_k(t,x_s,y_s,x_r,y_r)＝a×y_obs(t,x_s,y_s,x_r,y_r)+[I-a×M(x_s,y_s,x_r,y_r)]F_t^-1T_k(t,x_s,y_s,x_r,y_r)×F_td_k-1(t,x_s,y_s,x_r,y_r),k＝1,2,…,N.

其中a为加权因子，其范围在0<a≤1，如果a＝1，则方程与常规的POCS算法一样，此时原始噪声数据被带入到重建后的地震数据中，只能进行数据重建，不能进行噪声压制。不同的数据加权因子的选择不一样，与噪声能量的强度有关。

步骤4：阈值参数公式。不同的阈值参数会获得不同的重建效果，而合适的阈值参数在满足精度要求下，可以减少迭代次数并节省计算成本，因此，阈值参数的选取工作尤为重要。本发明在前人基础上，结合傅里叶变换的特点，选用按规律衰减的阈值参数，也即按照指数平方根规律衰减的阈值参数，在保证重建精度的前提下，可以更快提高收敛速度，节省部分计算时间，该阈值参数公式为：

$>{\lambda }_k=Max·e^{\left(\ln \right(ϵ)-ln(Max\left)\right]\sqrt{\frac{k-1}{N-1}}}$>

式中Max为|Cd|的最大值，即傅里叶变换系数绝对值的最大值。

实现该方法具体操作为：

为了更详细地阐述五维地震数据重建和噪声压制的效果，本发明定义信噪比SNR＝20log₁₀||f₀||₂/||f-f₀||₂来对比处理数据后的质量，单位为dB，其中f₀表示原始模型数据，f表示重建结果，信噪比越高，代表重建结果与模型数据越接近，效果越理想。

首先将本发明算法应用于五维数据理论模型，假设五维地震数据(x_s,y_s,x_r,y_r,t)包括三个双曲线同相轴，其中(x_s,y_s)和(x_r,y_r)代表偏移距(offset_x,offset_y)和共中心(CMP_x,CMP_y)点坐标，时间t为双程旅行时。为了提高运算速度，减少空间内存，该理论模型的空间大小设置为21道×21道×21道×21道，采样点为256个，采样率为1ms，空间采样间隔为10m，采用30Hz雷克子波进行模拟，由于很难显示五维数据图像，为此从中选择了多个共中点道集进行显示，图2表示CMP_y＝10，offset_y＝20时，三个共中心点道集(图2表示原始模拟地震数据图)，可以看出三条曲线同相轴较为连续，数据质量较好。图3表示为加了一定能量的高斯随机噪声(图3表示原始加噪数据)，模拟野外噪声干扰数据，然后对其进行20％随机欠采样，欠采样结果如图4所示(图4表示20％随机欠采样图)，此时信噪比为0.892dB，然后采用本发明方法进行数据重建与噪声压制，在重建和噪声压制的迭代过程中，仍然只对时间切片进行处理，减少一维正反傅立叶变换，使得在运算过程中降低数据重建的维数，节省内存空间。首先进行单独的五维地震数据重建，结果如图5所示(图5表示五维地震数据重建结果)，可以看出重建效果较好，同相轴较为连续，重建后信噪比6.82dB，其信噪比较低的原因是没有对噪声干扰进行压制，此时相当于加权因子a＝1。为了对比同时数据重建与噪声压制的效果，对重建后的含噪数据进行五维阈值迭代法单独噪声压制，噪声压制结果如图6所示(图6表示阈值法五维地震数据噪声压制结果)，可以看出噪声压制效果较好，随机噪声几乎可以去除干净，同相轴非常连续，信噪比为15.31dB，但是以上处理都是数据重建与噪声压制都是分开进行，没有统一。最后利用本发明方法进行同时数据重建与噪声压制，此时的加权因子a＝0.4，处理结果如图7所示(图7表示同时五维地震数据重建与噪声压制结果)，可以看出尽管缺失80％地震道，在某一维方向地震道缺失严重，但是由于五维地震数据重建方法可以充分利用四维空间信息对某一空间方向进行重建，缺失道地震数据恢复较好，并且噪声干扰也得到了压制，信噪比为16.07dB，处理后的效果优于先重建后噪声压制的效果，从而也说明五维地震数据同时重建与噪声压制的优越性，为了进一步体现其效果，也采用本发明方法对三维地震数据进行同时重建与噪声压制，其结果如图8所示(图8表示三维地震数据同时重建与噪声压制结果)，信噪比为8.84dB，可以看出三维重建算法在采样率非常低的情况下同时重建和噪声压制效果较差，而五维数据重建由于多利用了二维地震数据的空间信息，因此重建和噪声压制效果更好，这也说明通过对时间切片进行处理，采用POCS算法同时进行五维地震数据重建和噪声压制的方法效果显著，能够应用于实际资料处理。

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