一种基于线性Bregman算法的非均匀曲波三维地震数据重建方法与流程

本发明涉及的是空间非均匀采样下不规则缺失道的地震数据重建方法，具体是一种基于线性bregman算法的非均匀曲波三维地震数据重建方法。

技术背景

在野外数据采集过程中，为了得到相对完整而规则的地震数据，在数据采集之前必须进行野外观测系统的设计，但由于采集设备、野外地形条件以及经济成本限制等原因，地震炮点和检波点通常会偏离原始设计位置，甚至有些炮点检波点无法采集到有效的地震数据，从而导致地震数据沿空间方向常进行不规则欠采样，出现空间假频，影响到后续其它处理方法的效果，降低了地震勘探的分辨率(trad,2009；张华等，2013)。为了解决这种问题，最直接也是最有效的方法就是在野外重新进行数据采集，但从经济角度出发，显然不可能重新进行数据采集来完美的解决该问题。因此，必须在室内采取相应的不规则地震数据重建方法，使得缺失道数据得到有效地恢复。然而在目前已有的重建方法中，包括预测滤波方法(spitz,1991；naghizadehandsacchi,2007)，降秩方法(oropezaandsacchi,2011；ma,2013)，数学变换方法(白兰淑等，2014；唐欢欢等,2014)等，很少能够重建空间非均匀采样下的地震缺失道。但由于野外复杂地形条件的限制或者海上电缆的水平偏移，很多情况下野外地震数据常进行空间非均匀采样，如不加处理则会引起覆盖次数的变化(地下不均匀照明)，在叠加成像时会形成扭曲的成像振幅(采集脚印)，加重了空间假频现象的出现，影响后续成像处理。

为了解决空间非均匀采样下地震道不规则重建问题，地球物理领域常规处理方法为共面元叠加，从而将非均匀采样数据归位到均匀采样数据中来，满足后续其他处理方法的要求。然而共面元叠加处理方法忽略了每个面元内各道共中心点的真实位置，改变了部分地震道的振幅和相位，从而导致部分地震道位置出现严重偏差，降低了地震勘探资料的分辨率。另外一种方法就是基于波动方程的重建方法，然而该方法需要地下结构的先验信息，计算量非常巨大，对采样率要求也较高，从而也不能较好的解决该问题(ronen,1987)。尽管如此，许多学者采用基于数学变换的重建方法对该问题进行处理，duijndam等(1999)提出基于傅立叶变换的二维非均匀数据重建方法，hindriks和duijndam(2000)将其扩展到三维，实现了三维非均匀采样重建技术。但是，duijndam等人的傅立叶重建方法有其局限性，重建结果受最低速度和空间带宽的影响很大，随着采样间隔的逐渐增大，重建结果会逐渐变差。xu等(2005)采用重新正交化的过程，提出了基于抗泄露傅里叶变换的重建方法，但该方法对具有严重假频的地震数据则重建效果不好。zwartjes等(2007)也提出了反假频的非均匀地震数据重建方法，该方法采用无假频的低频信息来重建有假频的高频信息，达到压制假频和恢复缺失道的目的，但当低频也具有假频时，则该方法则会失效。jin(2010)提出基于阻尼最小范数傅立叶反演下的五维地震数据重建，该方法引入非均匀傅里叶变换方法，能够重建空间非均匀采样下的不规则缺失地震数据，但该方法抗假频能力不强。况且以上方法都是采用傅里叶变换作为稀疏基，尽管计算速度较快，但只适合处理近似线性同相轴或者平稳变化的地震信号，不能解决非线性同相轴或者非平稳地震数据的重建问题。

曲波变换能够表征信号的局部细节特征，可以有效地重建非线性同相轴或者非平稳变化的地震数据，众多的研究结果也证明，基于曲波变换的数据重建方法效果显著(naghizadehandsacchi,2010；刘国昌等，2011；张华等,2015)。尽管如此，以往基于曲波变换的重建方法前提条件仍然是空间均匀采样，而对于空间非均匀采样信号，由于缺乏空间连续性，曲波变换则不能有效地探测出地震波前特征，从而导致以往曲波变换不能有效地重建空间非均匀采样下的地震缺失道，从而限制了该方法的进一步应用。

技术实现要素：

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

本发明提出了一种基于线性bregman算法的非均匀曲波三维地震数据重建方法，首先针对常规二维曲波变换方法难以对空间非均匀采样下地震缺失道进行重建以及基于二维非均匀曲波变换的重建方法精度不高，应用不广泛等问题。沿着两个空间方向依次抽取非均匀三维地震数据的时间切片，在多尺度多方向二维曲波正变换的基础上，引入二维空间非均匀快速傅里叶变换，建立均匀曲波系数与空间非均匀采样下地震缺失数据之间的非均匀曲波反变换算子，然后采用线性bregman算法进行求解，通过选择合适的阈值因子和动态步长，并采用软阈值算子，精确地反演计算得到非均匀采样下不规则地震数据的均匀曲波系数，最后再进行常规曲波反变换，从而形成了一种基于线性bregman算法的非均匀曲波三维地震数据重建方法。

进一步，地震数据的重建问题就是从不完整的数据中恢复出完整的地震数据，假设如下线性正演模型

y＝md

这里y∈rm代表采集的不完整地震数据；d∈rⁿ，且n＞＞m，表示待重建的完整数据；m∈r^n×m表示随机采样矩阵。假设数据x是d在曲波变换域c中的稀疏表示，则上述方程可以写成：

y＝ax且

这里上标h代表共轭转置矩阵。

从上述方程可知，以往重建方法采用的算子c^h为常规均匀曲波反变换算子。事实上，该算子也可以为非均匀曲波反变换算子。为此，本发明采用和来分别表示常规均匀曲波反变换算子和非均匀曲波反变换算子。

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

式中：ψj,l,k表示曲波函数，c(j,l,k)为曲波系数，j,l,k分别表示尺度，方

向和位置参数，d(x)表示地震数据，其反变换为：

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

在这个方程中，f表示二维傅里叶变换，它实现了离散曲波变换第1步。t表示曲波拼接算子，即将频率波数域变换到曲波系数的过程，它实现了离散曲波变换第2～4步。

由于常规曲波正反变换算子满足因此可以定义均匀曲波反变换算子：

其中f^h表示二维傅里叶反变换，t^h表示曲波平铺算子。

在上述方程中，由于二维傅立叶变换参与了快速离散曲波变换之中，因此常规二维曲波变换不能处理非均匀采样数据，然而可以沿着时间切片的两个空间方向采用二维非均匀快速傅里叶反变换代替常规二维快速傅立叶反变换算子f^h。因此新的非均匀曲波反变换算子可以定义为：

该算子可以将离散曲波系数与非均匀采样下不规则地震道建立相应的联系。此时由于重建方程的稀疏解可以求解以下l1范数最优化问题得到：

subjecttoy＝ax

从而可以得出均匀曲波系数。在这个表达式中，x代表估计值，l1范数定义为x[i]是向量x中第i个元素。而对于该方程的求解，本发明选用线性bregman算法来求解此最小化问题；

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

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

进一步，所述非均匀快速傅里叶变换主要实现过程如下：首先将非均匀地震数据与某高斯短滤波器进行褶积，并对其结果进行均匀网格下密集采样，然后对密集采样后的数据进行快速傅里叶变换到频谱域，最后在频谱域进行反褶积校正，得到非均匀地震数据的频谱。

进一步，所述的线性bregman算法如下：

首先可以将l1范数下的最优化问题转为求解下述bp规则化问题

其中，λ是一个阈值权衡因子，在平衡l1范数和l2范数起到重要作用，算法1给出了线性bregman算法求解非均匀曲波变换过程的伪代码。

算法1.求解非均匀曲波变换的线性bregman算法

该函数用来处理数据中的噪声，动态步长tk被定义为：

算法1中的阈值为

算法1中的软阈值函数为

sλ＝sign(x)·max(|x|-λ,0)

本发明首先沿着两个空间方向依次抽取非均匀三维地震数据的时间切片，在多尺度多方向二维曲波正变换的基础上，引入二维空间非均匀傅里叶变换，建立均匀曲波系数与空间非均匀采样下地震缺失数据之间的非均匀曲波反变换算子，然后采用线性bregman算法进行求解，通过选择合适的阈值因子和动态步长，并采用软阈值算子，精确地反演计算得到非均匀采样下不规则地震数据的均匀曲波系数，最后再进行常规曲波反变换，从而形成了一种基于线性bregman算法的非均匀曲波三维地震数据重建方法。

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

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

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

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

4、相比于常规二维非均匀重建方法，本发明技术充分利用三维地震数据信息，从不同空间方向对缺失道进行重建，进一步提高重建精度，使得该发明在实际应用范围更为广泛、适用。

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

综述，本发明克服了常规二维曲波变换方法不能重建非均匀采样下地震缺失道的缺点以及常规二维非均匀重建方法重建精度低，应用不广泛的问题。并且该发明方法不仅可以重建非均匀带假频的缺失数据，而且也可以将非均匀网格数据归位到任意指定的均匀采样网格，大幅度提高了重建信号的分辨率和信噪比，保护了微弱的有效波信号，从而使反射波同相轴更加连续。这对于指导复杂地区非均匀地震数据采集、缺失道重建等方面具有重要的实用价值，同时对于其它非均匀信号处理领域也具有重要的借鉴意义。

附图说明

图1是本发明实施例中非均匀数据重建流程图。

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

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

图4是非均匀曲波变换三维地震数据规则化重建结果图。

图5是非均匀三维地震数据重建结果误差图。

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

图7是50％非均匀采样下三维不规则缺失地震数据重建结果图。

图8是非均匀采样下规则缺失三维地震数据图。

图9是非均匀采样下规则缺失三维地震数据重建结果图。

具体实施方式

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

实施例1

实现该方法的步骤主要包括，地震数据重建模型建立，常规二维曲波正反变换，非均匀快速傅里叶变换，建立非均匀曲波正反变换算子，地震波场重建，线性bregman算法求解等。具体步骤如下：

步骤1：数据重建模型建立。地震数据的重建问题就是从不完整的数据中恢复出完整的地震数据，假设如下线性正演模型

y＝md

这里y∈r^m代表采集的不完整地震数据；d∈rⁿ，且n＞＞m，表示待重建的完整数据；m∈r^n×m表示随机采样矩阵。假设数据x是d在曲波变换域c中的稀疏表示，则上述方程可以写成：

y＝ax且

这里上标h代表共轭转置矩阵。

从上述方程可知，以往重建方法采用的算子c^h为常规均匀曲波反变换算子，事实上，该算子也可以为非均匀曲波反变换算子。为此，本发明采用和来分别表示常规均匀曲波反变换算子和非均匀曲波反变换算子。

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

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

步骤3：非均匀傅里叶变换。首先将非均匀地震数据与某高斯短滤波器进行褶积，并对其结果进行均匀网格下密集采样，然后对密集采样后的数据进行快速傅里叶变换到频谱域，最后在频谱域进行反褶积校正，得到非均匀地震数据的频谱。

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

在这个方程中，f表示二维傅里叶变换，它实现了离散曲波变换第1步。t表示曲波拼接算子，即将频率波数域变换到曲波系数的过程，它实现了离散曲波变换第2～4步。

由于常规曲波正反变换算子满足因此可以定义均匀曲波反变换算子：

其中f^h表示二维傅里叶反变换，t^h表示曲波平铺算子。

在上述方程中，由于二维傅立叶变换参与了快速离散曲波变换之中，因此常规二维曲波变换不能处理非均匀采样数据，然而可以沿着时间切片的两个空间方向采用二维非均匀快速傅里叶反变换代替常规二维快速傅立叶反变换算子f^h。因此新的非均匀曲波反变换算子可以定义为：

该算子可以将离散曲波系数与非均匀采样下不规则地震道建立相应的联系。此时由于重建方程的稀疏解可以求解以下l1范数最优化问题得到：

subjecttoy＝ax

从而可以得出均匀曲波系数。在这个表达式中，x代表估计值，l1范数定义为x[i]是向量x中第i个元素。而对于该方程的求解，本发明选用线性bregman算法来求解此最小化问题；

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

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

步骤6：线性bregman算法。首先可以将l1范数下的最优化问题转为求解下述bp规则化问题

其中，λ是一个阈值权衡因子，在平衡l1范数和l2范数起到重要作用，算法1给出了线性bregman算法求解非均匀曲波变换过程的伪代码。

算法1.求解非均匀曲波变换的线性bregman算法

该函数用来处理数据中的噪声，动态步长tk被定义为：

算法1中的阈值为

算法1中的软阈值函数为

sλ＝sign(x)·max(|x|-λ,0)

实现该方法具体操作为：

为了衡量数据重建方法的效果，定义信噪比公式snr＝20log10||d0||2/||d-d0||2来进行对比，单位为db，其中d0表示原始模型数据，d表示重建结果，信噪比越高，代表重建结果与模型数据越接近，处理效果越理想。

首先建立一个二维不均匀速度模型，设置256炮，每炮256道接收，炮距和道距都为12米，采样率4毫秒，采用声波有限差分算法得到该模型所有的二维地震记录，然后对所得到的正演地震数据沿检波器，炮点以及时间坐标排列成三维数据体。图2为原始理论三维地震数据图，由于不能全部显示理论的三维数据模型，只能从三个不同的方向进行显示，其中时间切片为0.44s，炮点和检波点距离都为1524m。首先为了验证本发明所提方法的重建效果，对理论三维地震数据的时间切片进行二维均匀傅里叶变换，然后再进行空间二维非均匀快速傅里叶反变换，得到新的空间非均匀采样下的256×256道三维地震数据，如图3所示(图3表示非均匀三维地震数据图)，此时信噪比为6.47db，尽管名义上的采样网格(道距和炮距)还是12×12米，但每道地震数据道距不均匀，其道距范围为0m～24m，显然，如果将非均匀采样下地震数据在均匀采样网格上进行显示，连续的地震波场则会被破坏，而常规重建方法的前提条件是均匀道距下不规则缺失的地震记录。为此采用本发明方法进行规则化重建，规则化重建后的采样网格为12m×12m，结果如图4所示(图4表示非均匀曲波变换三维地震数据规则化重建结果图)，规则化重建后的信噪比为42.23db，地震波场连续性显著提高，重建后的地震记录与原始记录非常接近。图5表示非均匀三维地震数据重建结果误差图，可以看出误差几乎忽略不计，重建后信噪比非常高，与原始地震记录相比较几乎没有视觉上的差异，从而说明本发明方法能够反映出地震波场的局部细节特征，因此重建方法精度高，保真度较好。

为了检验本发明方法在非均匀采样下不规则缺失道重建效果，同样采用空间非均匀傅里叶变换，并且对其进行50％非均匀采样，得到50％非均匀采样下的三维不规则地震数据，如图6所示(图6表示50％非均匀采样下的三维不规则缺失地震数据)，此时三维地震记录炮距或道距范围为0m～84m，可以看出地震数据道距极不均匀，同相轴不连续，不能直接应用于后续资料的处理，必须采用非均匀三维地震数据重建方法进行处理。为此利用本发明方法进行重建，重建后的采样网格为12m×12m，重建结果如图7所示(图7表示50％非均匀采样下三维不规则缺失地震数据重建结果图),重建后的信噪比为21.04db，可以看出尽管不均匀采样下地震数据缺失50％地震道，但是重建效果精度仍然较高，重建后同相轴更连续，能量损失较少。

为了进一步检验本发明方法的反假频能力，对原始理论非均匀三维地震数据进行规则欠采样，如图7所示(图7表示非均匀采样下规则缺失三维地震数据图)，规则欠采样带来较为严重的假频成分，与信号的真实频谱存在部分重叠，然后采用本发明方法进行重建，重建后道距为12m×12m，重建结果如图9所示(图9表示非均匀采样下规则缺失三维地震数据重建结果图)，重建后信噪比为18.55db，重建效果较好，缺失的地震道得到了有效的恢复，表明本发明方法具有较强的反假频能力，能够进行复杂地区非均匀采样下的不规则和规则缺失地震数据重建。

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