一种基于流形分区2D-VMD的沙漠复杂噪声压制方法与流程
本发明属于一种沙漠地区复杂噪声压制方法。
背景技术：
地震勘探是地球物理勘探和油气勘探的重要途径之一。油气勘探是石油工业的基础，在国民经济中发挥着重要作用。近几十年来，大部分油田已进入高产阶段，主要老油田的产量正在减少，新油田的开发是一项紧迫而重要的任务。位于我国西北地区的塔里木油田是中国最大天然气产区和重要的石油生产基地。其所在的塔里木盆地是中国最大的含油气盆地，总面积56万平方公里。是中国最大天然气产区和重要的石油生产基地。
在地震勘探数据处理中，通常选用雷克子波来模拟地震子波，其有效信号主频通常设在30hz到100hz之间。一道地震数据中包含多个不同频率的子波，多道地震数据排列在一起就构成了多道地震记录，各个子波在多道地震记录中以同相轴形式呈现，反映了地下地层的分布情况。由于在地震数据采集过程中周围环境及地质条件不同，采集到的地震数据中会混入大量噪声，复杂的勘探环境和地质地表条件会造成含噪数据中噪声性质的复杂化。
与山地和平原地区相比，沙漠地区的地质地表条件相对复杂。在沙漠地区，蜂窝状沙丘波动很大，地震波的吸收和衰减更为严重，对地震资料的信噪比和分辨率有很大影响。造成沙漠地区噪声具有明显的非高斯、非线性和非平稳特性，且其频率集中在15hz以下。噪声与信号在时域、频域及空域都重叠在一起，这些复杂特性给常规的噪声压制方法带来了很大的困难。
目前有多种噪声抑制算法主要包括wiener滤波，卡尔曼滤波，自适应滤波，中值滤波，多项式拟合，小波变换，curvelet变换，shearlet变换和contourlet变换，时频峰值滤波，经验模态分解(emd)。虽然这些方法在多数条件下都能实现噪声压制，但都存在一定的应用条件或参数限制，如weiner滤波在非平稳噪声条件下性能明显下降，小波变换等几种变换域方法中阈值选取方式的确定直接影响消噪效果，时频峰值滤波方法对于较低信噪比情况下及频率较低噪声压制效果较差，经验模态分解方法存在的模态混叠问题导致其对较低频噪声的压制效果也不够理想。综上所述，现有大多数消噪方法对复杂低频沙漠噪声都无法获得有效的消噪效果
技术实现要素：
本发明提供一种基于流形分区2d-vmd的沙漠复杂噪声压制方法，针对沙漠地区噪声具有复杂且低频的特点，以解决现有消噪技术无法实现低频噪声的有效压制和有效信号恢复问题。
本发明采取的技术方案是，包括下列步骤：
步骤一：沙漠地区勘探数据中噪声具有非平稳、非线性、非高斯及低频的复杂特性，根据该地区噪声特点，通过对含噪地震数据进行多数据片段划分，将其与各个片段的位置信息合并构成高维空间的多维数据x；
选取合适的参数n和m，将含噪地震数据记录yn×m均匀地分割成大小为n×m的子数据集合；
yn×m＝{y1,y2…yq},(1)
其中yi(i＝1,…,q)为n×m的子数据记录，q＝(n×m)/(n×m)，再将每个子数据yi转换成q×1(q＝n×m)的一维数据片段yi，于是含噪地震数据记录yn×m就被划分为许多数据片段的集合，设各个一维数据片段yi在含噪地震数据记录yn×m中对应的位置信息为l＝{l1,l2…lq}，将其与众多数据片段合并看成高维空间的多维数据x＝[yl]＝{xi}；
步骤二：利用流形降维中的等距特征映射(isomap)算法对高维空间的多维数据x进行非线性降维，获取含噪地震数据记录yn×m的低维空间表示，得到含噪地震数据记录yn×m的最优低维空间矩阵topt，使含噪地震数据记录yn×m中的含信号数据片段与不含信号数据片段特征差异更突出更明显，将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复，得到分区后的含信号数据矩阵s；
首先，计算高维数据中两个样本xi和xj的欧式距离dx(xi,xj)，设定半径ε，如果xj在xi的半径ε之内或者是xi的c个最近邻点，则连接xi和xj；否则将该点与其非近邻点之间的距离设为无穷大，从而构造高维空间的多维数据x的c点邻域图g；
接着，计算c点近邻图g中任意两个样本xi和xj之间的测地距离为，两个样本点xi和xj之间的最短路径dg(xi,xj)为：
dg(xi,xj)＝min{dg(xi,xj),dg(xi,xa)+dg(xa,xj)},a＝1,2,…,n×m，(2)
由最短路径得到测地距离矩阵为：
对测地距离矩阵进行变换可得：
h＝-(i-ppt)dg(i-ppt)/2,(4)
其中p＝{1,2,…,n×m}，i为单位阵；
选取低维维数为d，计算h的最大d个特征值λ1,…,λd和特征向量μ1,…,μd，并构成特征向量矩阵u＝[μ1,…,μd]，于是含噪地震数据记录yn×m的低维流形为：
计算不同维数对应的低维流形与含噪地震数据记录yn×m之间的残差值，将最小残差对应的维数dmin作为最佳低维维数，得到含噪地震数据记录yn×m的最优低维空间矩阵topt；
在矩阵topt中，含信号数据片段与不含信号的数据片段对应的特征点呈现不同的空间分布特点，利用密度空间聚类(dbscan)算法将最优低维空间矩阵topt中的特征点聚为两类，得到含信号特征点类和不含信号特征点类，不含信号特征点直接清零，将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复，得到分区后的含信号数据矩阵s，该矩阵中纯噪声区域数据为零；
步骤三：由于含信号数据矩阵s中噪声与信号在时域、频域和空域均存在重叠，噪声很难消减；利用二维变分模态分解(2d-vmd)算法对其进行模态分解，实现不同形态、不同频率数据划分，将有效模态子信号重构，实现噪声压制；
含信号数据矩阵s经过2d-vmd分解后可得到多个谱带宽度有限的二维模态子信号，各模态子信号与其之间关系可表示为：
其中k为分解最大阶数，k为分解阶数，uk为二维模态子信号，各模态子信号均为频带有限信号，且具有不同的中心频率，频带之间虽有少部分重叠，但各频带内成分具有很好的聚集性，为了将信号频率调至基带，需将各模态子信号uk按下式变换为解析信号：
其中*表示卷积，<·>表示内积，δ(·)为狄拉克函数，⊥表示正交方向；接着，将解析信号乘以e指数函数进行调制，使其频谱调至基频，并令其沿频率方向梯度的l2范数达到最小值，从而构建优化函数如下：
其中
本发明首先根据沙漠噪声特点利用isomap流形学习算法将含噪地震记录划分为含信号区域和不含信号区域，并将不含信号区域看成纯噪声清零；再对含信号区域内数据利用2d-vmd算法进行多模态分解，实现不同形态、不同频率数据划分，最后将几个有效子模态信号求和恢复地震信号，实现噪声压制。
本发明的有益效果是本发明根据沙漠地区噪声自身的复杂特点，提出了一种基于流形分区2d-vmd的沙漠地震数据噪声压制方法，该方法从二维角度对地震数据进行分析处理，考虑了相邻地震道间同相轴的相关性及同相轴形态对消噪效果的影响，同时在频域内可实现混叠较小的频谱分解过程，更适用于低频噪声压制。与现有消噪方法相比，本发明能够有效实现复杂低频噪声的有效压制及有效地震信号的清晰连贯恢复。
附图说明
图1(a)是流形降维实现含噪地震数据分区；纯净合成地震记录，包含两个同相轴，层速度分别为1500m/s和2000m/s；
图1(b)是含噪合成地震记录，其中包含噪声为实际沙漠噪声，可以看到噪声振动幅度较大；
图1(c)是流形降维后含信号数据和不含信号数据的位置分布图，其中红色表示含信号片段，黑色表示纯噪声片段，两类数据在三维低维空间表示中具有明显的位置分布差异；
图1(d)是划分出的含信号数据区域，每个数据块大小为64×5个数据点，可以看出划分出的含信号区域中同相轴保持较为完整；
图2(a)是图1(d)中含信号数据区域的wiggle图；
图2(b)是含信号区域数据进行4阶2d-vmd分解后的最低频段的模态子信号，主要包含低频噪声；
图2(c)是含信号区域数据进行4阶2d-vmd分解后的较低频段的模态子信号，主要包含有效信号；
图2(d)是含信号区域数据进行4阶2d-vmd分解后的较高频段的模态子信号，主要包含部分有信号及噪声；
图2(e)是含信号区域数据进行4阶2d-vmd分解后的最高频段的模态子信号，主要包含高频噪声；
图2(f)是将图2(c)和图2(d)中的模态子信号求和恢复出的有效地震信号；
图3(a)是合成含噪地震记录，共有五个地震同相轴，主频均为20hz，包含实际的沙漠噪声；
图3(b)是emd消噪结果，其中模态数选为4，消噪结果为模态2和模态3构成；部分噪声得到了有效压制，但仍有许多噪声残留下来；
图3(c)是vmd消噪结果，其中模态数选为4，消噪结果为模态2构成，消噪效果优于emd方法，但有些地震道的噪声压制效果不是很好；
图3(d)是本发明方法消噪结果，其中2d-vmd中模态数选为4，无信号处噪声被清零消除，同相轴得到了完整保持；
图4是三种方法单道波形对比图；其中emd方法消噪后噪声得到了压制，但信号幅度也发生的严重衰减；vmd方法消噪后噪声的压制效果比emd方法要好些，但是信号幅度也发生的衰减；本发明方法消噪后信号幅度保持较完整，纯噪声数据段数据直接清零，实现噪声消除；
图5(a)是原始含噪地震数据；可以看出初至时刻之前的随机噪声振动幅度较大，有效信号淹没在噪声当中；
图5(b)是emd算法消噪结果，高频噪声得到了压制，但是同相轴仍然没有得到清晰恢复；
图5(c)是vmd算法消噪结果，噪声压制效果比emd消噪结果要好一些，同相轴变得更清晰，但是连贯性仍然不是很好；
图5(d)是本发明中方法消噪结果，低频噪声也得到了有效压制，隐藏的有效同相轴显现出来，连贯性明显增强。
具体实施方式
包括下列步骤：
步骤一：沙漠地区勘探数据中噪声具有非平稳、非线性、非高斯及低频的复杂特性，根据该地区噪声特点，通过对含噪地震数据进行多数据片段划分，将其与各个片段的位置信息合并构成高维空间的多维数据x；
选取合适的参数n和m，将含噪地震数据记录yn×m均匀地分割成大小为n×m的子数据集合；
yn×m＝{y1,y2…yq},(1)
其中yi(i＝1,…,q)为n×m的子数据记录，q＝(n×m)/(n×m)，再将每个子数据yi转换成q×1(q＝n×m)的一维数据片段yi，于是含噪地震数据记录yn×m就被划分为许多数据片段的集合，设各个一维数据片段yi在含噪地震数据记录yn×m中对应的位置信息为l＝{l1,l2…lq}，将其与众多数据片段合并看成高维空间的多维数据x＝[yl]＝{xi}；
仿真生成一张40道，每道1024个采样点的合成纯净地震记录，如图1(a)所示，采用雷克子波模拟地震子波，采样间隔为1ms，道间距为30米。记录中共包含两个同相轴，主频均为20hz，层速度分别为1500m/s和2000m/s，对其加入实际沙漠地区噪声得到含噪地震记录如图1(b)所示。选取n＝64，m＝5将其分成了多个子块数据记录，每个子记录包含320个样本点，有的子记录为纯噪声，有的子记录既包含噪声又包含信号；接着，将每个子记录转换成320×1的一维向量，并与其位置信息相结合构成高维数据，图1(b)中分块后的各个子记录数据片段集的位置信息为其在记录中位置中心点的横纵坐标；
步骤二：利用流形降维中的等距特征映射(isomap)算法对高维空间的多维数据x进行非线性降维，获取含噪地震数据记录yn×m的低维空间表示，即低维流形，使含噪地震数据记录yn×m中的含信号数据片段与不含信号数据片段特征差异更突出更明显；
首先，计算高维数据中两个样本xi和xj的欧式距离dx(xi,xj)，设定半径ε，如果xj在xi的半径ε之内或者是xi的c个最近邻点，则连接xi和xj；否则将该点与其非近邻点之间的距离设为无穷大，从而构造高维空间的多维数据x的c点邻域图g；
接着，计算c点近邻图g中任意两个样本xi和xj之间的测地距离为，两个样本点xi和xj之间的最短路径dg(xi,xj)为：
dg(xi,xj)＝min{dg(xi,xj),dg(xi,xa)+dg(xa,xj)},a＝1,2,…,n×m，(2)
由最短路径得到测地距离矩阵为：
对测地距离矩阵进行变换可得：
h＝-(i-ppt)dg(i-ppt)/2,(4)
其中p＝{1,2,…,n×m}，i为单位阵；选取低维维数为d，计算h的最大d个特征值λ1,…,λd和特征向量μ1,…,μd，并构成特征向量矩阵u＝[μ1,…,μd]，于是含噪地震数据记录yn×m的低维流形为：
计算不同维数对应的低维流形与含噪地震数据记录yn×m之间的残差值，将最小残差对应的维数dmin作为最佳低维维数，得到含噪地震数据记录yn×m的最优低维空间矩阵topt；
在矩阵topt中，含信号数据片段与不含信号的数据片段对应的特征点呈现不同的空间分布特点，利用密度空间聚类(dbscan)算法将最优低维空间矩阵topt中的特征点聚为两类，得到含信号特征点类和不含信号特征点类，不含信号特征点直接清零，将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复，得到分区后的含信号数据矩阵s，该矩阵中纯噪声区域数据为零；
对图1(a)中的含噪数据进行流形降维后得到的含信号数据片段与不含信号数据片段的最优低维空间矩阵topt的分布图如图1(c)所示，其中最低维数选取为3，因此绘制的是三维分布图；利用密度空间聚类(dbscan)算法对三维分布图中数据点进行聚类；图1(c)中的数据点中空心圆点表示含信号数据特征点，实心圆点表示纯噪声数据特征点，将含信号数据特征点按照原始含噪数据分割的逆过程进行重构，恢复出含信号区域，纯噪声数据特征点对应数据可直接清零，恢复中的区域结果如图1(d)所示，从划分结果可以该算法对信号与噪声成分的识别与划分具有很好的准确性；
步骤三：由于含信号数据矩阵s中噪声与信号在时域、频域和空域均存在重叠，噪声很难消减；利用二维变分模态分解(2d-vmd)算法对其进行模态分解，实现不同形态、不同频率数据划分，将有效模态子信号重构，实现噪声压制；
含信号数据矩阵s经过2d-vmd分解后可得到多个谱带宽度有限的二维模态子信号，各模态子信号与其之间关系可表示为：
其中k为分解最大阶数，k为分解阶数，uk为二维模态子信号，各模态子信号均为频带有限信号，且具有不同的中心频率，频带之间虽有少部分重叠，但各频带内成分具有很好的聚集性，为了将信号频率调至基带，需将各模态子信号uk按下式变换为解析信号：
其中*表示卷积，<·>表示内积，δ(·)为狄拉克函数，⊥表示正交方向；接着，将解析信号乘以e指数函数进行调制，使其频谱调至基频，并令其沿频率方向梯度的l2范数达到最小值，从而构建优化函数如下：
其中
选取模态阶数为，对图1(d)中划分信号区域后的地震数据进行2d-vmd分解后得到了各个模态子信号如图2所示。图2(a)为图1(d)中数据的wiggle多道波形图，图2(b)到图2(e)为利用2d-vmd对图2(a)中含噪声数据进行4阶分解后得到的4个模态子信号，可以看出各个模态子信号的频率由低到高变化，且所包含的同相轴形态不同，低频噪声主要集中在图2(a)中，高频噪声主要集中在图2(e)中。将图2(c)和图2(d)中的两个模态子信号求和便可得到图2(f)中的重构有效地震信号。对比含噪信号与重构结果可以看出，低频噪声得到了有效压制，同时有效信号得到了完整保持。
实验例1合成地震记录
仿真生成一张80道，每道2048个采样点的合成地震记录，采样频率为1000hz，对其加入一些实际的沙漠噪声，得到的含噪记录如图3(a)所示。在该含噪记录中共有五个地震同相轴，主频均为20hz。分别利用emd、一维vmd和本发明中方法对其进行消噪处理，去噪结果如图3(b)，3(c)和3(d)所示。从三种方法的多道结果可以看出，emd方法消噪后，部分噪声得到了有效压制，但仍有许多噪声残留下来。一维vmd方法处理后，其消噪效果优于emd方法，但有些地震道的噪声压制效果不是很好。相比之下，本发明中方法处理之后，无信号处噪声被清零消除，同相轴区域噪声也得到了有效抑制，并且得到了完整保持。
接着，从含噪记录，emd消噪结果，vmd消噪结果及本发明方法消噪结果中分别提取出第15道地震数据，并截取第641到第1740个数据点绘制其时域波形如图4所示。图中分别用不同线条表示不同方法的消噪结果，可以看出emd方法消噪后噪声得到了压制，但信号幅度也发生的严重衰减。vmd方法消噪后噪声的压制效果比emd方法要好些，但是信号幅度也发生的衰减。本发明方法消噪后信号幅度保持较完整，此外在纯噪声数据段数据直接清零，实现噪声消除。
实验例2实际地震记录
为了验证本发明的实际应用效果，将该方法应用于图5(a)中的部分实际野外地震资料，该资料来自中国塔里木沙漠地区。该地震剖面包含216道，每道2250个样本点，采样频率为1000hz,含噪记录中存在较大的波动，同相轴形态改变，都是由低频噪声造成的。分别利用emd方法、一维vmd方法和本发明中方法对该含噪记录进行消噪处理，结果如图5(b)，5(c)和5(d)所示。为了方便观察消噪效果，我们圈出了3个代表性区域，其中区域a为初至前噪声部分，区域b和区域c内低频噪声造成了同相轴具有较大波动。
从消噪结果图可以看出，图5(b)中emd方法消噪后，区域a中初至前噪声得到了有效压制，区域b和区域c中部分高频噪声得到了压制，但是同相轴仍然没有得到清晰恢复。图5(c)中一维vmd方法消噪后，区域a中的初至前噪声也得到了一定程度的压制，区域b和区域c中的噪声压制效果比emd消噪结果要好一些，同相轴变得更清晰，但是连贯性仍然不是很好。图5(d)中本发明方法消噪后，不仅区域a中的初至噪声得到了有效压制，区域b和区域c中的低频噪声也得到了有效压制，隐藏的有效同相轴显现出来，连贯性明显增强。这种满意的效果主要归功于2d-vmd方法不仅能够实现对信号频率上的划分，同时能够考虑相邻道间的相关性，在消噪过程中为同相轴连贯性的保持及恢复起到了重要作用。
- 还没有人留言评论。精彩留言会获得点赞！