CN108983287A - 一种基于凸集投影算法的曲波变换反假频地震数据重建方法 - Google Patents

Abstract

本发明公开了一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，针对规则缺失的假频地震数据或采样间隔较大的假频地震数据，首先将曲波变换分解成f‑k算子和曲波拼接算子，从f‑k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图并拾取分布图中的边界峰值，该边界代表有效波能量边界，再选择合适的阈值，根据边界峰值内的有效波能量建立相应的mask函数，最后将mask函数引入到传统凸集投影算法中，得到新的凸集投影算法；在迭代过程中，采用指数阈值参数公式，利用硬阈值算子重建出无假频的地震数据。本发明提高了重建信号的信噪比和计算速度，保护了微弱的有效波信号，使反射波同相轴更加连续、清晰。

Description

一种基于凸集投影算法的曲波变换反假频地震数据重建方法

技术领域

本发明涉及的是规则缺失和缺失间隔较大的地震数据重建方法，具体是一种基于凸集投影算法的曲波变换反假频地震数据重建方法。

技术背景

在野外数据采集过程中，为了得到相对完整而规则的地震数据，在数据采集之前必须进行野外观测系统的设计，但由于采集设备、野外地形条件以及经济成本限制等原因，地震炮点和检波点通常会偏离原始设计位置，甚至有些炮点检波点无法采集到有效的地震数据，从而导致地震数据沿空间方向常进行不规则欠采样，出现空间假频(Trad,2009；Zhaoet al.,2013)。这种空间欠采样的假频地震数据严重影响到后续的处理方法，比如多次波压制、速度分析、波动方程偏移等。为了解决这种问题，最直接也是最有效的方法就是在野外重新进行数据采集，但从经济角度出发，显然不可能重新进行数据采集来完美的解决该问题。因此，必须在室内提出高效率的反假频地震数据重建方法，使得缺失道数据得到有效地恢复。其实，从另外一个方面来讲，如果有了高效率的反假频重建方法，则可以在数据采集阶段主动缺失部分地震道，节省野外勘探成本，而将这种成本转移到室内数据重建中来。

目前已经有了很多数据重建方法，这些数据重建方法可以分为波动方程和数据处理两大类，其中数据处理方法主要包括预测滤波方法(Spitz,1991；Naghizadeh andSacchi,2007)，降秩方法(Oropeza and Sacchi,2011；Ma,2013)，傅里叶变换，曲波变换(Hennenfent,2008)，机器学习(Yu et al.,2016)等。然而这些方法尽管得到了广泛应用，但是只适合随机欠采样的地震数据重建，而对于规则欠采样数据重建，由于缺少反假频机制，因此会产生一个相对较差的重建结果。为了解决这个问题，Gulunay(2003)通过利用低频信息在f-k来重建假频地震数据。Zwartjes and Sacchi (2007)结合Gulunay方法，采用傅里叶反演实现反假频数据重建。基于这个思想。Naghizadeh(2007)也提出了多步自动回归的方法进行反假频重建。Curry(2010)根据沿着频率轴中同相轴的连续性，实现了基于傅立叶变换的径向自适应阈值策略重建方法。Naghizadeh(2012)通过利用整个频带信息，在f-k域进行角度扫描，建立一个面具函数进行约束假频，从而实现了反假频数据重建。Chiu(2014)在最小加权范数方法的基础上提出利用先验模型作为约束来解决反假频问题。然而以上所有反假频重建方法都是基于傅里叶变换，本质上来讲，这些方法都需要假设地震数据为线性或近似线性同相轴组成，尽管可以通过开设窗口函数的方法来重建非线性同相轴地震数据，但操作较为复杂，在窗口重叠部分重建效果较差，特别是对于那些缺失间隔较大的假频地震数据。

曲波变换能够表征信号的局部细节特征，可以有效地重建非线性同相轴或者非平稳变化的地震数据，众多的研究结果也证明，基于曲波变换的数据重建方法效果显著(Naghizadeh and Sacchi,2010；张华等，2015)。然而在他们的文章中，都是采用随机欠采样策略来进行地震数据重建，而对于规则欠采样的地震数据则重建效果较差，即：不具有反假频能力。因此，需要发展有一种有效的反假频数据重建方法。而在数据重建过程中，凸集投影算法被广泛应用，但根据凸集投影算法的原理得知，该算法也不具有反假频功能，不能重建规则缺失的假频数据，为此，本发明采用曲波变换方法，利用全频段数据建立mask函数，提出一种基于凸集投影算法的曲波变换反假频地震数据重建方法，在保证重建精度和效率的同时，使得该方法能够有效地重建出规则缺失的严重假频地震数据。

发明内容

本发明的目的是为了能够高精度重建野外规则缺失或采样间隔较大所产生的假频地震数据，大幅度提高重建信号的信噪比和计算速度，保护微弱的有效波信号，从而使反射波同相轴更加连续、清晰，而提出了一种基于凸集投影算法的曲波变换反假频地震数据重建方法。

本发明提出了一种基于凸集投影算法的曲波变换反假频地震数据重建方法，针对规则缺失的假频地震数据或采样间隔较大的假频地震数据，首先将曲波变换分解成f-k算子和曲波拼接算子，从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取能量分布图中的边界峰值，该边界代表有效波能量边界，然后选择合适的阈值，根据边界峰值内的有效波能量建立相应的mask函数，最后将mask函数引入到传统凸集投影算法中，得到新的凸集投影算法。在迭代过程中，采用指数阈值参数公式，利用硬阈值算子重建出无假频的地震数据，从而形成了一种基于凸集投影算法的曲波变换反假频地震数据重建方法。

进一步，曲波变换C分解为：

C＝QF

其中，F表示f-k算子，Q表示曲波平铺算子。

进一步，从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取能量分布图中的边界峰值；

假设d(t,x)表示原始地震数据，D(f,k)表示其频率波数(f-k)域，利用整个频带信息，在f-k域上进行角度扫描。角度射线的起始位置位于f-k域原点，同时根据频率轴的对称特性，采用正频率及归一化频率和波数解释。沿角度射线对f-k系数求和得到能量函数表达式：

其中，p是f-k域路径求和的扫描斜率，参数n代表归一化频率的索引，N表示最大频率索引。算子代表最小方向的整数值。在扫描过程中，在相邻p值之间需要选择足够细的步长来获取能量的变化范围。

函数E(p)中的峰值是较强有效波能量的表示，需要寻找出最左边和最右边两个峰值，并将这两个峰值之内的能量都考虑为有效波能量，假设已经识别出了2个边界峰值所对应的斜率，如p_L1，p_L2，下标L₁和L₁表示p值的索引。

进一步，在确定了函数E(p)的有效波能量范围后，在有效波能量范围外，mask函数值为零；在有效波能量范围内，mask函数表达式定义为：

其中，τ为用户定义的阈值，该值可以保留大于τ值的有效波能量系数。通过这个mask函数可以消除假频干扰，保留原始信号。

进一步，传统凸集投影算法表达式为：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQFd^i-1(t,x),i＝1,2,···,M

其中M表示最大迭代次数，dⁱ(t,x)表示第i次迭代得到的重建数据，y^obs表示原始采集到的地震数据，满足y^obs＝d⁰(t,x)，I表示单位矩阵，C^T表示曲波反变换算子。S(t,x)表示采样矩阵，Tⁱ表示硬阈值算子，表达式为：

cⁱ表示第i次迭代得到的曲波系数，满足cⁱ＝Cdⁱ(t,x)，λⁱ表示第i次迭代阈值，其公式为：

式中Max为曲波变换系数绝对值的最大值，ε为接近零的小值。

进一步，由于规则欠采样会在频率域产生与有效波振幅能量相同的假频干扰，难以通过阈值法提取有效波信号。因此，传统凸集投影算法不能够重建规则缺失的假频地震数据，需要将mask函数引入到传统凸集投影迭代算法中，可得新的凸集投影算法表达式为：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQHFd^i-1(t,x),i＝1,2,···,M

mask函数能够消除由规则欠采样产生的假频干扰，使得每一次迭代过程中有效波能量聚焦，因此，新的表达式能够成功地重建出规则缺失的假频地震数据。

针对规则缺失的假频地震数据或采样间隔较大的假频地震数据，本发明首先将曲波变换分解成f-k算子和曲波拼接算子，从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取能量分布图中的边界峰值，该边界代表有效波能量边界，然后选择合适的阈值，根据边界峰值内的有效波能量建立相应的mask函数，最后将mask函数引入到传统凸集投影算法中，得到新的凸集投影算法。在迭代过程中，采用指数阈值参数公式，利用硬阈值算子重建出无假频的地震数据，从而形成了一种基于凸集投影算法的曲波变换反假频地震数据重建方法。

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

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

2、相对于傅里叶变换重建方法，本发明技术能够高精度地处理非线性同相轴地震数据以及缺失间隔较大的地震数据，并且重建精度高。

3、相对于基于凸集投影算法的曲波变换重建方法，本发明技术不仅能够重建规则缺失的假频地震数据，而且大幅度提高了计算效率，节省了运算时间，提高了重建后的信噪比，保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰，

综述，本发明克服了常规重建方法不能重建规则缺失假频地震数据的缺点以及以及常规重建方法重建精度低，应用不广泛的问题。而且该发明方法不仅可以重建规则缺失的假频地震数据，而且也可以重建缺失道间隔较大的假频地震数据，大幅度提高了重建信号的计算速度和信噪比，节省了运算时间，保护了微弱的有效波信号，从而使反射波同相轴更加连续、清晰。这对于指导复杂地区的地震数据采集、缺失道重建等方面具有重要的实用价值。

附图说明

图1是本发明实施例中反假频数据重建流程图。

图2是原始理论地震数据及其频谱图。

图3是规则缺失50％数据重建结果图。

图4是重建结果与信噪比关系曲线图。

图5是f-k域角度扫描能量分布图。

图6是f-k域mask函数图。

图7是规则缺失75％数据重建结果图。

具体实施方式

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

实施例1

实现该方法的步骤主要包括，曲波变换分解，角度扫描，mask函数的建立，传统凸集投影算法、新的凸集投影算法等。具体步骤如下：

步骤1：曲波变换分解。曲波变换C分解为：

C＝QF

其中F表示f-k算子，Q表示曲波平铺算子。

步骤2：角度扫描。从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取边界峰值。

假设d(t,x)表示原始地震数据，D(f,k)表示其频率波数(f-k)域，该发明利用整个频带信息，在f-k域上进行角度扫描。角度射线的起始位置位于f-k域原点，同时根据频率轴的对称特性，本发明仅用正频率及归一化频率和波数解释该方法。沿角度射线对f-k系数求和得到能量函数表达式：

其中，p是f-k域路径求和的扫描斜率，参数n代表归一化频率的索引，N表示最大频率索引。算子代表最小方向的整数值。在扫描过程中，在相邻p值之间需要选择足够细的步长来获取能量的变化范围。

函数E(p)中的峰值是较强有效波能量的表示，需要寻找出最左边和最右边两个峰值，并将这两个峰值之内的能量都考虑为有效波能量，假设已经识别出了2个边界峰值所对应的斜率，如p_L1，p_L2，下标L₁和L₁表示p值的索引。

步骤3：mask函数的建立。在确定了函数E(p)的有效波能量范围后，在有效波能量范围外，mask函数值为零；在有效波能量范围内，mask函数表达式定义为：

其中，τ为用户定义的阈值，该值可以保留大于τ值的有效波能量系数。通过这个mask函数可以消除假频干扰，保留原始信号。

步骤4：传统凸集投影算法。传统凸集投影算法表达式为：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQFd^i-1(t,x),i＝1,2,···,M

其中M表示最大迭代次数，dⁱ(t,x)表示第i次迭代得到的重建数据，y^obs表示原始采集到的地震数据，满足y^obs＝d⁰(t,x)，I表示单位矩阵，C^T表示曲波反变换算子。S(t,x)表示采样矩阵，Tⁱ表示硬阈值算子，表达式为：

cⁱ表示第i次迭代得到的曲波系数，满足cⁱ＝Cdⁱ(t,x)，λⁱ表示第i次迭代阈值，其公式为：

式中Max为曲波变换系数绝对值的最大值，ε为接近零的小值。

步骤5：新的凸集投影算法。由于规则欠采样会在频率域产生与有效波振幅能量相同的假频干扰，难以通过阈值法提取有效波信号。因此，传统凸集投影算法不能够重建规则缺失的假频地震数据，需要将mask函数引入到传统凸集投影迭代算法中，可得新的凸集投影算法表达式：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQHFd^i-1(t,x),i＝1,2,···,M

mask函数能够消除由规则欠采样产生的假频干扰，使得每一次迭代过程中有效波能量聚焦，因此，新的表达式能够成功地重建出规则缺失的假频地震数据。

实现该方法具体操作为：

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

图2a为256道合成地震记录，该记录总共有4层地震反射波，每一层反射波能量有所差异，采样间隔为1ms，道距为4m，每道1024个采样点。然后对其进行50％规则欠采样，如图2b所示。图2c为原始理论数据的f-k频谱，图2d为规则欠采样后的f-k频谱，从中可以看出规则欠采样所引起的假频与真实频谱较为接近，不具有反假频的重建方法则很难对其进行处理，恢复出缺失地震数据，提高地震剖面的信噪比。首先采用本发明方法提出的反假频数据重建方法进行重建，迭代次数选择15次，重建结果如图3a所示，图3b为其f-k频谱，重建后的信噪比为21.34dB，明显可以看出本发明方法能够有效地将缺失道信息恢复出来，并且由于采用了多尺度多方向的曲波变换方法，重建后的频谱能量与原始数据频谱非常接近。

为了对比本发明方法与传统凸集投影方法的区别，分别采用不同的迭代次数进行重建，得出迭代次数与重建后信噪比关系曲线图，如图4所示，从图中可以看出在相同的迭代次数下，本发明方法重建后信噪比始终高于传统凸集投影方法，并且只要迭代12次就能够得到理想的重建效果，再增加迭代次数已经没有太大的意义，而传统方法至少需要迭代40次才能够得到较高的信噪比，并且最高的信噪比也只有18.75dB，而这个信噪比的取得，本发明方法只要迭代7次，从而也说明本发明提出的方法聚焦收敛速度快，重建精度高。

为了详细展示本发明方法如何得到图3a的重建结果，图5为从图1b计算出来的E(p)能量函数，其中-5≤p≤5，采样间隔为0.08，在能量分布图中可以清楚看出规则欠采样后有效信号的左右边界情况，该边界范围内的能量值表示有效信号能量，根据图5中有效信号能量范围，选择合适的阈值，建立mask函数，如图6所示，大于该信号阈值区域的值为1，小于阈值区域的值为0。然后将该mask函数引入到传统凸集投影方法中，使得在每次凸集投影算法迭代过程中，该mask函数可以很精确地拾取出原始信号的f-k频谱，再参与下一次迭代计算中，使得原始信号的有效波频谱能够很快聚焦，减少迭代次数，提高重建精度，并在保证精度的同时，其计算效率远高于传统凸集投影算法。

同时为了进一步检验本发明方法的反假频能力，将原始数据进行75％规则缺失，使缺失间隔进一步加大，其缺失数据如图7a所示，图7c为其f-k频谱，可以看出由于缺失道较多，使得在频率域中有效波能量与假频能量相互干扰，非常严重。然后采用本发明方法进行重建，重建结果如图7b所示，重建后的信噪比分别为17.58dB，重建效果较好，缺失道得到了有效的恢复，并且从频率域(图7d)也可以看出，重建后的有效波能量与原始信号能量几乎相同，表明重建误差较小，从而也进一步说明本发明方法重建效果较好，能够满足实际应用要求。

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

Claims (6)

1.一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，针对规则缺失的假频地震数据或采样间隔较大的假频地震数据，首先将曲波变换C分解成f-k算子和曲波拼接算子，从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取能量分布图中的边界峰值，该边界代表有效波能量边界，然后选择合适的阈值，根据边界峰值内的有效波能量建立相应的mask函数，最后将mask函数引入到传统凸集投影算法中，得到新的凸集投影算法；在迭代过程中，采用指数阈值参数公式，利用硬阈值算子重建出无假频的地震数据。

2.根据权利要求1所述的一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，曲波变换C分解为：

C＝QF

其中F表示f-k算子，Q表示曲波平铺算子。

3.根据权利要求1所述的一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，从f-k域原点出发，利用全频段信息，采用角度扫描策略得出能量分布图，并拾取能量分布图中的边界峰值；

假设d(t,x)表示原始地震数据，D(f,k)表示其频率波数(f-k)域，利用整个频带信息，在f-k域上进行角度扫描；角度射线的起始位置位于f-k域原点，同时根据频率轴的对称特性，采用正频率及归一化频率和波数解释；沿角度射线对f-k系数求和得到能量函数表达式：

其中，p是f-k域路径求和的扫描斜率，参数n代表归一化频率的索引，N表示最大频率索引；算子代表最小方向的整数值；在扫描过程中，在相邻p值之间需要选择足够细的步长来获取能量的变化范围；

函数E(p)中的峰值是较强有效波能量的表示，需要寻找出最左边和最右边两个峰值，并将这两个峰值之内的能量都考虑为有效波能量，假设已经识别出了2个边界峰值所对应的斜率，如下标L₁和L₁表示p值的索引。

4.根据权利要求3所述的一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，在确定了函数E(p)的有效波能量范围后，在有效波能量范围外，mask函数值为零；在有效波能量范围内，mask函数表达式定义为：

其中，τ为用户定义的阈值，该值可以保留大于τ值的有效波能量系数；通过这个mask函数可以消除假频干扰，保留原始信号。

5.根据权利要求1所述的一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，传统凸集投影算法表达式为：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQFd^i-1(t,x),i＝1,2,···,M

其中M表示最大迭代次数，dⁱ(t,x)表示第i次迭代得到的重建数据，y^obs表示原始采集到的地震数据，满足y^obs＝d⁰(t,x)，I表示单位矩阵，C^T表示曲波反变换算子；S(t,x)表示采样矩阵，Tⁱ表示硬阈值算子，表达式为：

cⁱ表示第i次迭代得到的曲波系数，满足cⁱ＝Cdⁱ(t,x)，λⁱ表示第i次迭代阈值，其公式为：

式中Max为曲波变换系数绝对值的最大值，ε为接近零的小值；

由于规则欠采样会在频率域产生与有效波振幅能量相同的假频干扰，难以通过阈值法提取有效波信号；因此，传统凸集投影算法不能够重建规则缺失的假频地震数据，需要将mask函数引入到传统凸集投影迭代算法中。

6.根据权利要求1所述的一种基于凸集投影算法的曲波变换反假频地震数据重建方法，其特征在于，新的凸集投影算法表达式为：

dⁱ(t,x)＝y^obs(t,x)+[I-S(t,x)]C^TTⁱQHFd^i-1(t,x),i＝1,2,···,M

mask函数能够消除由规则欠采样产生的假频干扰，使得每一次迭代过程中有效波能量聚焦，因此，新的表达式能够成功地重建出规则缺失的假频地震数据。