一种基于自适应复锐化扩散的沙漠地区地震勘探去噪方法与流程

本发明属于沙漠地区地震勘探随机噪声的压制方法，具体涉及的是沙漠地区地震勘探记录中随机噪声压制。本发明利用结构张量获取地震勘探数据纹理结构信息，进而构建复数域锐化扩散滤的自适应扩散系数，从而实现随地震勘探数据纹理结构特征自适应扩散滤波和信号保持。

背景技术：

随着易采易探的油气藏的减少和地震勘探技术的不断发展，地震勘探逐步向复杂恶劣地表地质地区推进，其勘探难度在不断加大，需要不断改进地震勘探数据处理方法来适应越来越复杂的数据和地质情况。沙漠地区蕴含大量的油气资源，近年来对沙漠地区的地震勘探活动日益增多。由于沙漠地区采集条件恶劣，沙层松散，地表受风速影响引起沙丘形状变化，噪声幅值大，因此，所收集到的地震勘探资料不同程度地被强随机噪声所污染，信噪比低。与其他地区地震勘探随机噪声相比，沙漠随机噪声具有非平稳，非线性，时空非均匀性等特征，且频率低，相关性强，和有效信号具有相似的时域波形和交叉的频带。沙漠随机噪声复杂的性质严重影响有效信号的辨识和拾取，不宜直接利用野外地震勘探资料做地质解释，需要对地震勘探数据进行噪声压制处理，提高地震勘探资料的信噪比，从中提取有效信息，为下一步的地质解释提供可靠的数据资料。

当前在地震勘探随机噪声压制处理方面有很多优秀的去噪技术，代表性的如小波变换，p-m扩散滤波，时频峰值滤波等。小波变换去噪方法将地震勘探资料分解到不同频率段的尺度上，并在各个层次设置相应的阈值，实现小波域地震勘探资料的随机噪声消减。但是因为沙漠随机噪声与有效信号频域范围有重叠部分，存在于相同尺度，且能量较强，信号和噪声很难分离。时频峰值滤波将地震勘探信号调制后进行时频域滤波，能够有效压制低信噪比地震勘探记录中的随机噪声。该方法是建立在高斯噪声基础上的，而沙漠低频随机噪声在时空滤波矢量空间呈非高斯分布，能量强，很难被有时频峰值滤波有效压制。p-m滤波作为经典非线性扩散滤波，该方法不受信号和噪声频带重叠的影响，同时可以在去除噪声的同时保留信号边缘等纹理结构部分。但是会出现阶跃效应，产生虚假边缘。复锐化扩散滤波是以偏微分方程为基础的滤波方法，其将复数域扩散滤波与shock滤波器结合，包含两部分，一部分为shock锐化，另一部分为复数域扩散平滑，两者同时作用。利用复锐化扩散滤波结果的虚部控制shock锐化器，解决了shock锐化器对噪声敏感的缺点，同时shock锐化器对信号的锐化增强作用改善了扩散部分易过分平滑信号的缺点。然而，当复锐化扩散滤波应用于沙漠地区地震勘探数据处理时，由于地震勘探数据具有丰富的纹理结构特征如同相轴结构，且沙漠随机噪声具有与同相轴结构相似的波形和频带，在消除沙漠随机噪声的同时保留这些纹理结构信息有助于地震资料的解释说明，而复数域扩散部分的扩散系数为定值，不能识别同相轴而做出相应改变，所以在扩散过程中会将这些有用的纹理结构信息同沙漠随机噪声一样平滑掉，无法有效恢复甚至增强同相轴信息。

通过分析上述方法表明，现有的地震勘探降噪方法在处理沙漠地区地震勘探随机噪声时面对的共同问题是沙漠随机噪声固有的复杂特性：低频率，强相关性，和有效信号具有相似的时域波形和交叉的频带。这严重影响地震勘探数据低频带信号的恢复，因此，需要探索有效压制沙漠地区地震勘探随机噪声，恢复弱有效信号，保持信号细节和结构地震勘探信号处理方法。

技术实现要素：

本发明提供一种基于自适应复锐化扩散的沙漠地区地震勘探去噪方法，以解决现有的地震勘探降噪方法在处理沙漠地区地震勘探随机噪声时存在的严重影响地震勘探数据低频带信号恢复的问题。

本发明采取的技术方案是，包括以下步骤：

步骤一：对沙漠地区含噪地震勘探数据i构建一致性判断算子d:

构建一致性判断算子d的具体步骤如下：

(1)对二维沙漠地区含噪地震勘探数据i求结构张量r，如下式：

其中为卷积算子，ix和iy表示分别对二维沙漠地区含噪地震勘探数据i的x方向和y方向求偏导数，之后采用标准差为ρ的高斯核gρ进行线性滤波，从而获得线性结构张量r，其中元素r12和r21相同，结构张量r是对称且正定的矩阵，可以获得该矩阵的两个标准正交的特征向量ω1和ω2,其中特征向量ω1平行于梯度方向，而特征向量ω2垂直于梯度方向，两个特征向量ω1和ω2分别对应特征值μ1和μ2，公式如下：

其中μ1代表地震勘探数据i在特征向量ω1上的变化强度；μ2代表地震勘探数据i在特征向量ω2上的变化强度；

(2)利用结构张量在梯度方向和一致性方向所对应的特征值μ1和μ2构建一致性判断算子d，用于分析地震勘探数据结构特征，一致性判断算子表达式如下：

d＝(μ1-μ2)^bμ1^q(5)

其中参数b和q随不同数据而做出调整，突出信号的一致性的同时将强噪声区域的一致性降到最低；

步骤二：利用一致性判断算子d来构建自适应阈值函数k，公式如下：

k(d)＝mexp(-d)(6)

其中函数exp(·)是以自然常数e为底的指数函数，参数m是调整阈值与复数域扩散结果虚部的尺度因子；

步骤三：对沙漠地区含噪地震勘探数据i进行自适应复数域锐化扩散滤波，其第n次迭代滤波的变化量如下：

其中表示第n次迭代过程中数据的变化量；a作为锐化系数；arctan(·)是反正切函数；im(·)表示虚部；θ为将实数域扩散引入到复数域的角度参数，一般设置为趋近于零的值；n为设定的迭代次数；表示对梯度求幅值大小；η和ξ分别表示数据梯度方向和一致性方向；为第n次迭代过程中的梯度方向扩散系数，为复数；为第n次迭代过程中的一致性方向扩散系数，为实数；和分别表示第n-1次迭代结果i^n-1在梯度方向和一致性方向上的二阶导；

第n次迭代滤波的结果：

其中δt表示每次迭代的时间步长；i^n-1表示第n-1次迭代结果；iⁿ表示第n次迭代结果；当n达到设定迭代次数n时，得到最终迭代滤波结果iⁿ为地震勘探数据的去噪结果。

本发明步骤三所述的借助自适应阈值函数k来进行构建自适应梯度方向扩散系数表示如下：

其中表示第n次迭代获得的自适应梯度方向扩散系数；w为自适应梯度方向扩散系数的权值。

本发明步骤三所述的借助自适应梯度方向扩散系数和结构张量的特征向量ω1构建自适应一致性方向扩散系数其中特征向量ω1＝(ω1x,ω1y)^t，并存在如下关系式：

自适应一致性方向扩散系数表示如下：

其中表示第n次迭代获得的自适应一致性方向扩散系数；λ1为自适应一致性方向扩散系数的权值；δ是用于判断界限的参数。

本发明从噪声传播特性出发，以复锐化扩散滤波为基础，提出自适应复锐化扩散作为压制沙漠地区随机噪声的方法，本发明借助结构张量构建复锐化扩散滤波在梯度方向上的扩散系数使其随地震勘探数据一致性程度非线性变化，从而控制对地震勘探数据信号区域和噪声区域的扩散程度。在此基础上，本发明进一步利用结构张量的特征向量调节一致方向上的扩散系数，使扩散强度随地震不同区域变化而变化，进而增强复锐化扩散对同相轴幅度和结构的保持能力。

本发明所解决的技术问题：复数域锐化扩散滤波因扩散部分在梯度方向和一致性方向上的扩散系数为定值，无法同时实现噪声压制和信号保护。针对这一问题,本发明提出梯度方向和一致性方向上的自适应扩散系数。基于构建的一致性算子，梯度方向上的扩散系数可以随着地震数据不同区域，即信号区域和噪声区域而变化：在噪声区域,梯度方向上的扩散系数设为最大值以进行强去噪；在信号区域，使梯度方向上的扩散系数趋近于零用于保护信号。针对大倾角同相轴衰减问题，本发明调节一致性方向扩散系数,使其在噪声区域处于最大值以增强去噪能力；在信号区域根据地震同相轴与竖直方向坐标轴的倾角大小自适应调整扩散系数，同相轴越陡峭其值越小，从而有效保护信号，改善去噪后的大倾角同相轴幅值衰减。本发明提出的自适应复锐化扩散能够根据地震数据结构自适应调整在梯度方向和一致性方向上的扩散系数，实现了对沙漠地区含噪地震勘探数据进行噪声压制处理的同时保护地震同相轴的目标。

本发明的优点：研究针对沙漠地区随机噪声与信号频域范围部分重叠而导致的压噪保幅效果差的问题，提出以沙漠地区地震勘探含噪数据纹理结构信息为变量的自适应复数域锐化扩散方法，该方法通过构造一致性判断算子来获取沙漠地区地震勘探含噪数据的纹理结构信息，来区分数据中的有效信号区域和噪声区域，并进一步将此纹理结构信息融入到自适应阈值函数的构建中，从而使扩散系数可以自适应地随地震勘探数据纹理结构变化，从而实现沙漠地震勘探噪声有效压制的同时保持有效信号同相轴更加清晰和准确的目的，具有实用性和有效性。

本发明可以去除常规地震勘探数据去噪方法难以去除的沙漠噪声，提高沙漠噪声地震勘探数据的信噪比。

附图说明

图1是本发明的流程图；

图2是无噪声的模拟地震勘探记录图；

图3是加沙漠噪声后的模拟地震勘探记录图；

图4是沙漠噪声频率谱图；

图5是小波变换滤波处理后的地震勘探记录图；

图6是p-m滤波方法处理后的地震勘探记录图；

图7是斜坡保持复扩散滤波处理后的地震勘探记录图；

图8是复数域锐化扩散滤波处理后的地震勘探记录图；

图9是自适应复锐化扩散滤波处理后的地震勘探记录图；

图10是第一个同相轴第9道幅值对比图；

图11是第一个同相轴第9道频谱对比图；

图12是第二个同相轴第38道幅值对比图；

图13是第二个同相轴第38道频谱对比图；

图14是沙漠地区实际地震勘探含噪记录图；

图15是小波变换滤波处理后的实际地震勘探记录图；

图16是p-m滤波方法处理后的实际地震勘探记录图；

图17是斜坡保持复扩散滤波处理后的实际地震勘探记录图；

图18是复数域锐化扩散滤波处理后的实际地震勘探记录图；

图19是自适应复锐化扩散滤波处理后的实际地震勘探记录图。

具体实施方式

包括以下步骤：

步骤一：对沙漠地区含噪地震勘探数据i构建一致性判断算子d:

一致性判断算子的构建是基于对结构张量的分析为基础的，结构张量作为一种有效的数据结构分析工具，矩阵值的结构张量可以准确估计数据中的边缘方向及其上的变化强度，可以区分数据中的边缘、各向同性区域，允许同时进行方向估计和结构分析。此处我们利用结构张量来获得沙漠地震勘探数据的局部结构信息，并利用所获得的信息来控制复数域扩散过程，以实现效果更好的滤波结果；

构建一致性判断算子d的具体步骤如下：

(1)对二维沙漠地区含噪地震勘探数据i求结构张量r，如下式：

其中为卷积算子，ix和iy表示分别对二维沙漠地区含噪地震勘探数据i的x方向和y方向求偏导数，之后采用标准差为ρ的高斯核gρ进行线性滤波，从而获得线性结构张量r，其中元素r12和r21相同。结构张量r是对称且正定的矩阵，可以获得该矩阵的两个标准正交的特征向量ω1和ω2,其中特征向量ω1平行于梯度方向，而特征向量ω2垂直于梯度方向。两个特征向量ω1和ω2分别对应特征值μ1和μ2，公式如下：

其中μ1为结构张量在平行于梯度方向上的特征值，代表数据i在特征向量ω1上的变化强度；μ2为结构张量在垂直于梯度方向上的特征值，代表数据i在特征向量ω2上的变化强度；

(2)利用结构张量在梯度方向和一致性方向所对应的特征值μ1和μ2构建用于分析地震勘探数据结构特征的一致性判断算子d；

数据如若沿梯度方向的变化强度远大于一致性方向的变化强度，则会表现出明显的边缘结构特征，而当梯度方向的变化强度与一致性方向的变化强度相差不大时,则表现出平坦区域，所以结构张量的特征值可以很好地将数据中结构特征突出的部分与结构特征弱的部分区分开，通过结构张量锁定边缘位置和边缘强度，然后根据所提取结构信息有针对性地配合复数域扩散进行去噪。

一致性判断算子d表达式如下：

d＝(μ1-μ2)^bμ1^q(5)

其中参数b和q随不同数据而做出调整，突出信号的一致性的同时将强噪声区域的一致性降到最低。随着信噪比降低，参数b和q设置相应升高；信噪比升高时，参数b和q设置相应降低；

步骤二：利用一致性判断算子d来构建自适应阈值函数k，以此实现在噪声区域阈值较大进行扩散去噪，而在信号区域阈值较小从而保护信号纹理结构不被破坏，公式如下：

k(d)＝mexp(-d)(6)

其中函数exp(·)是以自然常数e为底的指数函数，参数m是调整阈值与复数域扩散结果虚部的尺度因子，从而有效利用阈值函数中包含的纹理结构信息。

步骤三：对沙漠地区含噪地震勘探数据i进行自适应复数域锐化扩散滤波，其第n次迭代滤波的变化量如下：

其中扩散方程的第一项是shock锐化器，用于锐化增强；后两项是扩散部分，用于平滑去噪；表示第n次迭代过程中数据的变化量；a作为锐化系数；arctan(·)是反正切函数；im(·)表示虚部；θ为将实数域扩散引入到复数域的角度参数，一般设置为趋近于零的值；n为设定的迭代次数；表示对梯度求幅值大小；η和ξ分别表示数据梯度方向和一致性方向；为第n次迭代过程中的梯度方向扩散系数，为复数；为第n次迭代过程中的一致性方向扩散系数，为实数；和分别表示第n-1次迭代结果i^n-1在梯度方向和一致性方向上的二阶导；

在复数域锐化扩散滤波方法中，梯度方向和一致性方向上的两个扩散系数均为常数，在地震勘探数据的噪声区域和信号区域保持相同的扩散强度。一般为了避免信号被过分平滑而丢失细节，梯度方向上的扩散系数取值较小，而同时为了保有压噪的能力，在一致性方向上的扩散系数取值较大。但是这样为了保护信号，会大大衰减噪声压制效果。所以为了同时实现保护信号和压制噪声的目的，我们提出了可以随噪声区域和信号区域而自适应变化的梯度方向扩散系数，使其在噪声区域保持最大值1，以此实现强去噪，而在信号区域使其取值接近于零，以此实现保护信号的目的。在此基础上，我们又提出了一致性方向上的扩散系数，在噪声区域取最大值1以此实现强去噪，而在信号区域可以随着地震同相轴的倾斜角度大小而变化。因为地震同相轴越陡峭其幅值衰减越严重，所以为了保护陡峭同相轴的幅值，我们利用同相轴与竖直方向坐标轴的夹角来控制一致性方向上的扩散系数，使其在越陡峭的同相轴处扩散系数越小以此实现保幅的目的。

第n次迭代滤波的结果iⁿ为：

其中δt表示每次迭代的时间步长；i^n-1表示第n-1次迭代结果；iⁿ表示第n次迭代结果；当n达到设定迭代次数n时，得到最终迭代滤波结果iⁿ为地震勘探数据的去噪结果；

所述步骤三中借助自适应阈值函数k来进行构建自适应梯度方向扩散系数表示如下：

其中表示第n次迭代获得的自适应梯度方向扩散系数；w为自适应梯度方向扩散系数的权值。

所述步骤三中借助自适应梯度方向扩散系数和结构张量的特征向量ω1构建自适应一致性方向扩散系数其中特征向量ω1＝(ω1x,ω1y)^t，并存在如下关系式：

自适应一致性方向扩散系数表示如下：

其中表示第n次迭代获得的自适应一致性方向扩散系数；λ1为自适应一致性方向扩散系数的权值；δ是用于判断界限的参数。因为梯度方向扩散系数在噪声区域和信号区域有明显的分界，所以利用梯度方向扩散系数划分区域；在信号区域，利用上式中的反正切函数求解出地震同相轴与竖直坐标轴的夹角，以此实现在同相轴由平缓变为陡峭的过程中一致性方向扩散系数由大变小，从而解决同相轴在陡峭区域被严重衰减的问题，而在噪声区域，一致性方向扩散系数被赋予最大值1，有效压制噪声。

应用实例：

将本发明应用于模拟沙漠地区含噪地震勘探数据i的噪声压制处理，图3是由无噪声地震勘探数据图2通过加入沙漠噪声获得，道数为40道，每道包含600采样点，采样频率为1000hz。无噪声地震勘探数据由三条同相轴组成，从上往下的频率分别为30hz,25hz,20h。对图3中所加沙漠噪声进行频率谱分析，给出图4沙漠噪声频谱分析图，可以发现同相轴的频谱范围与沙漠噪声的低频频谱范围有重叠。本发明对具有这一特征的模拟沙漠地区含噪地震勘探数据i进行噪声压制，同时与经典小波变换算法，p-m模型算法，斜坡保持复扩散滤波以及复数域锐化扩散滤波进行比较。

五种方法的去噪后的数据见图5-图9。从图5可以观察到，小波变换算法虽然信号保持效果较好，但是存在噪声压制问题，其不能有效去除沙漠噪声，甚至产生一些畸变，如箭头所指向的区域。在图6-图8中，从箭头指示的区域可以观察到，如第一条同相轴的1-14道左右和26-36道左右，以及第二条同相轴的34-40道左右，p-m模型，斜坡保持复扩散滤波以及复数域锐化扩散滤波都存在信号衰减的问题，在陡峭同相轴的下降区域信号有明显的衰减，不能得到有效的恢复。而从图9可以看出经过自适应复数域锐化扩散滤波处理后的数据，在第一条同相轴的1-14道左右和26-36道左右，以及第二条同相轴的34-40道左右区域取得了较好的信号恢复结果，相较于p-m模型，斜坡保持复扩散滤波以及复数域锐化扩散滤波，明显降低了信号衰减，于此同时也获得了令人较满意的噪声压制的结果。

在给出整体去噪结果图的基础上，我们给出了以上5种方法的单道幅值对比图和单道频谱对比图，我们取了第一条同相轴的第9道(图10和图11)和第二条同相轴的第37道(图12和图13)，分别是整体去噪结果图中的箭头所指区域。从单道幅值对比图(图10和图12)我们可以观察到，经过p-m模型，斜坡保持复扩散滤波以及复数域锐化扩散滤波处理后的第一个同相轴和第二个同相轴的地震波，其波峰和波谷的幅值都有明显衰减，而经过自适应复数域锐化扩散滤波处理后，第一个同相轴和第二个同相轴的地震波波峰和波谷的幅值得到有效恢复，衰减得到明显改善。而从单道频谱对比图(图11和图13)我们可以观察到，经过p-m模型，斜坡保持复扩散滤波，复数域锐化扩散滤波以及自适应复数域锐化扩散滤波后，大部分低频沙漠噪声被有效压制。相较于p-m模型，斜坡保持复扩散滤波以及复数域锐化扩散滤波，经过自适应复数域锐化扩散滤波处理后，有效信号区域的频谱也得到有效保护。而经过小波变换处理后的数据，低频沙漠噪声频谱仍有大量残余，噪声没有得到有效清除。

将本发明应用于沙漠地区实际地震勘探数据的噪声压制处理，图14是沙漠地区实际地震勘探数据，道数为100道，每条道的采样点为560，采样频率为1000hz。本发明对该沙漠地区实际地震勘探数据进行噪声压制，同时与小波变换算法，经典p-m模型算法，斜坡保持复扩散滤波以及复数域锐化扩散滤波进行比较。

五种方法的去噪后的数据见图15-图19。从箭头指示的区域可以观察到，经过小波变换处理后的数据，沙漠地震勘探噪声仍有大量残余，噪声没有得到有效清除。p-m模型，斜坡保持复扩散滤波以及复数域锐化扩散三者能够压制大部分的沙漠地区地震勘探噪声，后两者的信号恢复效果相差不大，且比p-m模型效果更好。而经过自适应复数域锐化扩散滤波处理后的数据，相较于其他四种方法，在箭头指示的区域信号有更明显的恢复甚至增强，且明显压制了沙漠地震勘探噪声。

合成记录和实际沙漠地区共炮点记录表明,本发明提出的去噪方法在含有低频沙漠噪声的地震数据的处理中，较传统的去噪模型具有优越性，在有效清除沙漠噪声的同时，恢复甚至增强被沙漠噪声污染的有用信息。