CN110208856B - 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 - Google Patents
一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 Download PDFInfo
- Publication number
- CN110208856B CN110208856B CN201910489001.XA CN201910489001A CN110208856B CN 110208856 B CN110208856 B CN 110208856B CN 201910489001 A CN201910489001 A CN 201910489001A CN 110208856 B CN110208856 B CN 110208856B
- Authority
- CN
- China
- Prior art keywords
- data
- signal
- noise
- dimensional
- frequency
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 67
- 230000001629 suppression Effects 0.000 title claims abstract description 29
- 238000005192 partition Methods 0.000 title claims abstract description 10
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 27
- 230000008569 process Effects 0.000 claims abstract description 10
- 238000001228 spectrum Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 43
- 238000009826 distribution Methods 0.000 claims description 18
- 239000012634 fragment Substances 0.000 claims description 14
- 230000009467 reduction Effects 0.000 claims description 9
- 238000005457 optimization Methods 0.000 claims description 6
- 238000011084 recovery Methods 0.000 claims description 6
- 230000011218 segmentation Effects 0.000 claims description 6
- 230000002776 aggregation Effects 0.000 claims description 3
- 238000004220 aggregation Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 230000000877 morphologic effect Effects 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 claims description 3
- 238000003379 elimination reaction Methods 0.000 abstract description 20
- 230000008030 elimination Effects 0.000 abstract description 19
- 230000000694 effects Effects 0.000 abstract description 8
- 238000001914 filtration Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 230000006870 function Effects 0.000 description 6
- 238000005070 sampling Methods 0.000 description 5
- 230000002238 attenuated effect Effects 0.000 description 4
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 4
- 230000009466 transformation Effects 0.000 description 4
- 239000007789 gas Substances 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 239000003345 natural gas Substances 0.000 description 2
- 239000003208 petroleum Substances 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000004660 morphological change Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种基于流形分区2D‑VMD的沙漠复杂噪声压制方法,属于沙漠地区复杂噪声压制方法。根据沙漠噪声特点利用ISOMAP流形学习算法将含噪地震记录划分为含信号区域和不含信号区域,并将不含信号区域看成纯噪声清零;再对含信号区域内数据利用2D‑VMD算法进行多模态分解,实现不同形态、不同频率数据划分,最后将几个有效子模态信号求和恢复地震信号,实现噪声压制。有益效果是从二维角度对地震数据进行分析处理,考虑了相邻地震道间同相轴的相关性及同相轴形态对消噪效果的影响,同时在频域内可实现混叠较小的频谱分解过程,更适用于低频噪声压制,本发明能够有效实现复杂低频噪声的有效压制及有效地震信号的清晰连贯恢复。
Description
技术领域
本发明属于一种沙漠地区复杂噪声压制方法。
背景技术
地震勘探是地球物理勘探和油气勘探的重要途径之一。油气勘探是石油工业的基础,在国民经济中发挥着重要作用。近几十年来,大部分油田已进入高产阶段,主要老油田的产量正在减少,新油田的开发是一项紧迫而重要的任务。位于我国西北地区的塔里木油田是中国最大天然气产区和重要的石油生产基地。其所在的塔里木盆地是中国最大的含油气盆地,总面积56万平方公里。是中国最大天然气产区和重要的石油生产基地。
在地震勘探数据处理中,通常选用雷克子波来模拟地震子波,其有效信号主频通常设在30Hz到100Hz之间。一道地震数据中包含多个不同频率的子波,多道地震数据排列在一起就构成了多道地震记录,各个子波在多道地震记录中以同相轴形式呈现,反映了地下地层的分布情况。由于在地震数据采集过程中周围环境及地质条件不同,采集到的地震数据中会混入大量噪声,复杂的勘探环境和地质地表条件会造成含噪数据中噪声性质的复杂化。
与山地和平原地区相比,沙漠地区的地质地表条件相对复杂。在沙漠地区,蜂窝状沙丘波动很大,地震波的吸收和衰减更为严重,对地震资料的信噪比和分辨率有很大影响。造成沙漠地区噪声具有明显的非高斯、非线性和非平稳特性,且其频率集中在15Hz以下。噪声与信号在时域、频域及空域都重叠在一起,这些复杂特性给常规的噪声压制方法带来了很大的困难。
目前有多种噪声抑制算法主要包括Wiener滤波,卡尔曼滤波,自适应滤波,中值滤波,多项式拟合,小波变换,Curvelet变换,Shearlet变换和Contourlet变换,时频峰值滤波,经验模态分解(EMD)。虽然这些方法在多数条件下都能实现噪声压制,但都存在一定的应用条件或参数限制,如Weiner滤波在非平稳噪声条件下性能明显下降,小波变换等几种变换域方法中阈值选取方式的确定直接影响消噪效果,时频峰值滤波方法对于较低信噪比情况下及频率较低噪声压制效果较差,经验模态分解方法存在的模态混叠问题导致其对较低频噪声的压制效果也不够理想。综上所述,现有大多数消噪方法对复杂低频沙漠噪声都无法获得有效的消噪效果
发明内容
本发明提供一种基于流形分区2D-VMD的沙漠复杂噪声压制方法,针对沙漠地区噪声具有复杂且低频的特点,以解决现有消噪技术无法实现低频噪声的有效压制和有效信号恢复问题。
本发明采取的技术方案是,包括下列步骤:
步骤一:沙漠地区勘探数据中噪声具有非平稳、非线性、非高斯及低频的复杂特性,根据该地区噪声特点,通过对含噪地震数据记录YN×M进行多数据片段划分,将其与各个片段的位置信息合并构成高维空间的多维数据X;
步骤二:利用流形降维中的等距特征映射(ISOMAP)算法对高维空间的多维数据X进行非线性降维,获取含噪地震数据记录YN×M的低维空间表示,得到含噪地震数据记录YN×M的最优低维空间矩阵Topt,使含噪地震数据记录YN×M中的含信号数据片段与不含信号数据片段特征差异更突出更明显,在矩阵Topt中,含信号数据片段与不含信号的数据片段对应的特征点呈现不同的空间分布特点,利用密度空间聚类(DBSCAN)算法将最优低维空间矩阵Topt中的特征点聚为两类,得到含信号特征点类和不含信号特征点类,不含信号特征点直接清零,将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复,得到分区后的含信号数据矩阵S,该矩阵中纯噪声区域数据为零;将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复,得到分区后的含信号数据矩阵S;
步骤三:由于含信号数据矩阵S中噪声与信号在时域、频域和空域均存在重叠,噪声很难消减;利用二维变分模态分解(2D-VMD)算法对其进行模态分解,实现不同形态、不同频率数据划分,将有效模态子信号进行重构,实现噪声压制。
本发明步骤一具体包括:选取合适的参数n和m,将一张含噪地震数据记录YN×M均匀地分割成大小为n×m的子数据集合;
YN×M={Y1,Y2…YQ}, (1)
其中Yb为n×m的子数据记录,b=1,…,Q,Q=(N×M)/(n×m),再将每个子数据Yi转换成q×1的一维数据片段yi,q=n×m,于是含噪地震数据记录YN×M就被划分为许多数据片段的集合,设各个一维数据片段yi在含噪地震数据记录YN×M中对应的位置信息为L={l1,l2…lQ},将其与众多数据片段合并成高维空间的多维数据X={xi}。
本发明所述步骤二具体包括:
首先,计算高维数据中两个样本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。
本发明所述步骤三具体包括:
含信号数据矩阵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)
其中Yb(b=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={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范数达到最小值,从而构建优化函数如下:
其中为梯度,||·||2为L2范数,通过优化式(8)便可求出最优各模态子信号uk,不同的uk包含不同的频率成分及形态特征;由于沙漠地震记录中噪声频率较低,因此低频噪声被分解到低频模态信号中,并且地震记录中同相轴呈现不同的分布形态,因此不同的模态信号中包含的不同分布规律的地震数据,于是将包含有效频率成分并且有效分布形态的几个子信号进行相加便可实现有效地震信号的恢复,剩余几个子模态信号相加即为被消除的噪声。
选取模态阶数为,对图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方法不仅能够实现对信号频率上的划分,同时能够考虑相邻道间的相关性,在消噪过程中为同相轴连贯性的保持及恢复起到了重要作用。
Claims (4)
1.一种基于流形分区2D-VMD的沙漠复杂噪声压制方法,其特征在于:包括下列步骤:
步骤一:沙漠地区勘探数据中噪声具有非平稳、非线性、非高斯及低频的复杂特性,根据该地区噪声特点,通过对含噪地震数据记录YN×M进行多数据片段划分,将其与各个片段的位置信息合并构成高维空间的多维数据X;
步骤二:利用流形降维中的等距特征映射(ISOMAP)算法对高维空间的多维数据X进行非线性降维,获取含噪地震数据记录YN×M的低维空间表示,得到含噪地震数据记录YN×M的最优低维空间矩阵Topt,使含噪地震数据记录YN×M中的含信号数据片段与不含信号数据片段特征差异更突出更明显,在矩阵Topt中,含信号数据片段与不含信号的数据片段对应的特征点呈现不同的空间分布特点,利用密度空间聚类(DBSCAN)算法将最优低维空间矩阵Topt中的特征点聚为两类,得到含信号特征点类和不含信号特征点类,不含信号特征点直接清零,将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复,得到分区后的含信号数据矩阵S,该矩阵中纯噪声区域数据为零;将含信号特征点类中的数据点通过矩阵分割的逆过程进行恢复,得到分区后的含信号数据矩阵S;
步骤三:由于含信号数据矩阵S中噪声与信号在时域、频域和空域均存在重叠,噪声很难消减;利用二维变分模态分解(2D-VMD)算法对其进行模态分解,实现不同形态、不同频率数据划分,将有效模态子信号进行重构,实现噪声压制。
2.根据权利要求1所述的一种基于流形分区2D-VMD的沙漠复杂噪声压制方法,其特征在于:步骤一具体包括:选取合适的参数n和m,将一张含噪地震数据记录YN×M均匀地分割成大小为n×m的子数据集合;
YN×M={Y1,Y2…YQ}, (1)
其中Yb为n×m的子数据记录,b=1,…,Q,Q=(N×M)/(n×m),再将每个子数据Yb转换成q×1的一维数据片段yb,q=n×m,于是含噪地震数据记录YN×M就被划分为许多数据片段的集合,设各个一维数据片段yb在含噪地震数据记录YN×M中对应的位置信息为L={l1,l2…lQ},将其与众多数据片段合并成高维空间的多维数据X。
3.根据权利要求1所述的一种基于流形分区2D-VMD的沙漠复杂噪声压制方法,其特征在于,所述步骤二具体包括:
首先,计算高维数据中两个样本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。
4.根据权利要求1所述的一种基于流形分区2D-VMD的沙漠复杂噪声压制方法,其特征在于,所述步骤三具体包括:
含信号数据矩阵S经过2D-VMD分解后可得到多个谱带宽度有限的二维模态子信号,各模态子信号与其之间关系可表示为:
其中K为分解最大阶数,k为分解阶数,uk为二维模态子信号,各模态子信号均为频带有限信号,且具有不同的中心频率,频带之间虽有少部分重叠,但各频带内成分具有很好的聚集性,为了将信号频率调至基带,需将各模态子信号uk按下式变换为解析信号:
其中*表示卷积,<·>表示内积,δ(·)为狄拉克函数,⊥表示正交方向;接着,将解析信号乘以e指数函数进行调制,使其频谱调至基频,并令其沿频率方向梯度的L2范数达到最小值,从而构建优化函数如下:
由于沙漠地震记录中噪声频率较低,因此低频噪声被分解到低频模态信号中,并且地震记录中同相轴呈现不同的分布形态,因此不同的模态信号中包含的不同分布规律的地震数据,于是将包含有效频率成分并且有效分布形态的几个子信号进行相加便可实现有效地震信号的恢复,剩余几个子模态信号相加即为被消除的噪声。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910489001.XA CN110208856B (zh) | 2019-06-05 | 2019-06-05 | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910489001.XA CN110208856B (zh) | 2019-06-05 | 2019-06-05 | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110208856A CN110208856A (zh) | 2019-09-06 |
CN110208856B true CN110208856B (zh) | 2020-12-25 |
Family
ID=67791201
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910489001.XA Expired - Fee Related CN110208856B (zh) | 2019-06-05 | 2019-06-05 | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110208856B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110764147B (zh) * | 2019-11-06 | 2020-11-03 | 吉林大学 | 一种基于vmd局部f-x谱分解的沙漠勘探弱信号恢复方法 |
CN111723701B (zh) * | 2020-06-08 | 2022-05-20 | 西安交通大学 | 一种水中目标识别方法 |
CN112766044B (zh) * | 2020-12-28 | 2024-03-22 | 中海石油(中国)有限公司 | 疏松样品纵横波速度分析方法、装置及计算机存储介质 |
CN113156513B (zh) * | 2021-04-14 | 2024-01-30 | 吉林大学 | 一种基于注意力引导的卷积神经网络地震信号去噪方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013106720A1 (en) * | 2012-01-12 | 2013-07-18 | Schlumberger Canada Limited | Method for constrained history matching coupled with optimization |
CN108760316A (zh) * | 2018-08-16 | 2018-11-06 | 苏州大学 | 变分模态分解的变参信息融合方法 |
CN108845352A (zh) * | 2018-06-27 | 2018-11-20 | 吉林大学 | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 |
CN108919347A (zh) * | 2018-07-02 | 2018-11-30 | 东华理工大学 | 基于vmd的地震信号随机噪声压制方法 |
CN109657578A (zh) * | 2018-12-10 | 2019-04-19 | 北京航空航天大学 | 基于变分模态分解与拉普拉斯特征映射的特征提取方法 |
-
2019
- 2019-06-05 CN CN201910489001.XA patent/CN110208856B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013106720A1 (en) * | 2012-01-12 | 2013-07-18 | Schlumberger Canada Limited | Method for constrained history matching coupled with optimization |
CN108845352A (zh) * | 2018-06-27 | 2018-11-20 | 吉林大学 | 基于vmd近似熵与多层感知机的沙漠地震信号去噪方法 |
CN108919347A (zh) * | 2018-07-02 | 2018-11-30 | 东华理工大学 | 基于vmd的地震信号随机噪声压制方法 |
CN108760316A (zh) * | 2018-08-16 | 2018-11-06 | 苏州大学 | 变分模态分解的变参信息融合方法 |
CN109657578A (zh) * | 2018-12-10 | 2019-04-19 | 北京航空航天大学 | 基于变分模态分解与拉普拉斯特征映射的特征提取方法 |
Non-Patent Citations (8)
Title |
---|
Application of Variational Mode Decomposition in Random Noise Attenuation and Time-frequency Analysis of Seismic Data;W. Liu et al.;《78th EAGE Conference & Exhibition 2016》;20160602;全文 * |
Depositional sequence characterization based on seismic variational mode decomposition;Fangyu Li et al.;《Interpretation》;20170531;第SE97-SE106页 * |
Trace-transform-based time-frequency filtering for seismic signal enhancement in Northeast China;Ning Wu et al.;《Comptes Rendus Geoscience》;20160404;第360-367页 * |
变分模态分解及其在地震去噪中的应用;何元等;《中国地球科学联合学术年会 2014》;20141231;第1002页 * |
基于二维变分模态分解和自适应中值滤波的图像去噪方法;刘嘉敏等;《计算机应用研究》;20171031;第34卷(第10期);第3149-3152页 * |
基于低维流形的地震强噪声衰减;于四伟等;《中国地球科学联合学术年会 2017》;20171231;第728-731页 * |
基于变分模态分解与流形学习的滚动轴承故障特征提取方法;戚晓利等;《振动与冲击》;20181231;第37卷(第23期);第133-140页 * |
沙漠地区地震勘探随机噪声建模及其在噪声压制中的应用;李光辉等;《地球物理学报》;20160228;第59卷(第2期);第682-692页 * |
Also Published As
Publication number | Publication date |
---|---|
CN110208856A (zh) | 2019-09-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110208856B (zh) | 一种基于流形分区2d-vmd的沙漠复杂噪声压制方法 | |
CN107817527B (zh) | 基于块稀疏压缩感知的沙漠地震勘探随机噪声压制方法 | |
Tian et al. | Variable-eccentricity hyperbolic-trace TFPF for seismic random noise attenuation | |
US20080270033A1 (en) | Methods of hydrocarbon detection using spectral energy analysis | |
CN104849757A (zh) | 消除地震信号中随机噪声系统及方法 | |
CN103543469A (zh) | 一种基于小波变换的小尺度阈值去噪方法 | |
CN107144879A (zh) | 一种基于自适应滤波与小波变换结合的地震波降噪方法 | |
CN103364832A (zh) | 一种基于自适应最优核时频分布的地震衰减定性估计方法 | |
CN107179550B (zh) | 一种数据驱动的地震信号零相位反褶积方法 | |
CN110261910A (zh) | 基于自适应稀疏s变换的地震数据面波去除方法 | |
CN107783191A (zh) | 多维空间时空时频峰值滤波消减地震勘探随机噪声的方法 | |
Pecholcs et al. | A broadband full azimuth land seismic case study from Saudi Arabia using a 100,000 channel recording system at 6 terabytes per day: acquisition and processing lessons learned | |
CN105652322A (zh) | 多分量地震数据的t-f-k域极化滤波方法 | |
CN104570116A (zh) | 基于地质标志层的时差分析校正方法 | |
CN112327362A (zh) | 速度域的海底多次波预测与追踪衰减方法 | |
CN105319593A (zh) | 基于曲波变换和奇异值分解的联合去噪方法 | |
CN111708087A (zh) | 一种基于DnCNN神经网络对地震数据噪声压制的方法 | |
CN102338884B (zh) | 物探中的椭圆窗方向带通保幅滤波数据处理方法 | |
Zhang et al. | Seismic random noise attenuation and signal-preserving by multiple directional time-frequency peak filtering | |
CN104914471B (zh) | 适于黄土塬非纵测线的地滚波压制方法 | |
Wang et al. | Application of a new wavelet threshold method in unconventional oil and gas reservoir seismic data denoising | |
CN111257938A (zh) | 基于小波互相关时移地震虚拟震源波场重构方法和系统 | |
CN115840249A (zh) | 基于随机道矢量中值滤波及高分辨率Radon变换的黑三角噪音压制方法 | |
CN110764147B (zh) | 一种基于vmd局部f-x谱分解的沙漠勘探弱信号恢复方法 | |
CN112213785B (zh) | 一种基于特征增强去噪网络的地震资料沙漠噪声抑制方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201225 Termination date: 20210605 |
|
CF01 | Termination of patent right due to non-payment of annual fee |