基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法.pdf

2023 年 8 月第 58 卷 第 4 期基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法曹静杰*1，2，许昌昊1，3，朱跃飞1，2，4（1.自然资源部京津冀城市群地下空间智能探测与装备重点实验室，河北石家庄 050031；2.河北省战略性关键矿产资源重点实验室，河北石家庄 050031；3.河北地质大学数理学院，河北石家庄 050031；4.中国矿业大学（北京）地球科学与测绘工程学院，北京 100083）摘要：高信噪比、高保真度和高分辨率地震数据是获得地下结构清晰成像的前提。在地震数据采集过程中，不利地形、坏道等因素造成地震数据不满足采样定理，需重建以获得完备的地震数据。多道奇异谱分析法是一种常用的地震数据重建与去噪方法，其关键是对划分的每个数据块确定有效奇异值个数。该参数需根据不同的数据特征判定，但目前的人工选择方法需花费大量人力和计算资源。针对每个数据块，基于频率域多道奇异谱分析框架，分析奇异值的离散点曲线和谱的规律，提出层次聚类法自动识别有效信号对应的奇异值个数，提高地震数据去噪效果和重建质量。对经过块 Hankel矩阵奇异值分解得到的奇异值序列，采用层次聚类法进行聚类，得到各频率的有效奇异值个数；再对有效频率内求取的奇异值个数求最大，获得数据块的有效奇异值个数；在阻尼多道奇异谱分析同时重建和去噪的框架下，进而提出一种改进的多道奇异谱分析（MSSA）方法以实现地震数据的同时重建和去噪。对模拟地震数据与实际地震资料进行测试，结果表明基于层次聚类的多道奇异谱分析法在同时重建和去噪方面具有优越性，能避免人工选择有效的奇异值个数，减小地震数据处理工作量，对大规模地震数据的重建与去噪具有现实意义。关键词：多道奇异谱分析，重建，去噪，层次聚类，奇异值分析中图分类号：P631 文献标志码：A doi：10.13810/ki.issn.1000-7210.2023.04.007Simultaneous reconstruction and denoising of seismic data using multichannel singular spectrum analysis based on hierarchical clusteringCAO Jingjie1，2，XU Changhao1，3，ZHU Yuefei1，2，4（1.Key Laboratory of Intelligent Detection and Equipment for Underground Space of BeijingTianjinHebei UrbanAgglomeration，Ministry of Natural Resources，Shijiazhuang，Hebei 050031，China；2.College of Mathematicsand Physics，Hebei GEO University，Shijiazhuang，Hebei 050031，China；3.Hebei Key Laboratoryof Strategic Critical Mineral Resources，Shijiazhuang，Hebei 050031，China；4.College of Geoscience and Surveying Engineering，China University of Mining and TechnologyBeijing，Beijing 100083，China）Abstract：High signaltonoise ratio，high fidelity，and highresolution seismic data are prerequisites for clear imaging of subsurface structures.During seismic data acquisition，factors such as unfavorable topography and bad channels lead to the phenomenon that seismic data cannot satisfy the sampling theorem，and reconstruction is required to obtain complete seismic data.The multichannel singular spectrum analysis method is a common seismic data reconstruction and denoising method，and its key is to determine the number of effective singular values for each divided 处理技术 文章编号：1000-7210（2023）04-0818-12*河北省石家庄市裕华区槐安东路 136号河北地质大学地球科学学院，050031。Email：本文于 2022年 8月 8日收到，最终修改稿于 2023年 5月 15日收到。本项研究受国家自然科学基金项目“面向城市地质的三维地震勘探压缩感知采集设计与数据重建研究”（41974166）、河北省自然科学基金项目“基于深度学习和模型驱动的地震数据重建方法研究”（D2021403010）、河北省在读研究生创新能力培养资助项目“基于矩阵完备理论的地震数据同时插值和去噪方法研究”（CXZZSS2022020）、河北地质大学科技创新团队项目（KJCXTD202106）、河北省自然资源厅项目“基于光纤传感的地下空间智能监测方法与应用”联合资助。第 58 卷 第 4 期曹静杰，等：基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法data block.This parameter needs to be determined based on different data characteristics，but current manual selection methods require a lot of labor and computational resources.For each data block，based on the framework of multichannel singular spectrum analysis in the frequency domain，the discrete point curves and spectral patterns of singular values are analyzed，and a hierarchical clustering method is proposed to automatically identify the number of singular values corresponding to the effective signals，which improves the denoising effect and reconstruction quality of seismic data.Then，the hierarchical clustering method is adopted to cluster the singular value sequence obtained by singular value decomposition of the block Hankel matrix and acquire the number of effective singular values of each frequency.Additionally，the number of singular values taken within the effective frequency is maximized to obtain the effective singular values for data blocks.In the framework of simultaneous reconstruction and denoising by damped multichannel singular spectrum analysis，an improved multichannel singular spectrum analysis（MSSA）method is put forward to realize simultaneous reconstruction and denoising of seismic data.Tests conducted on simulated and actual seismic data show that the MSSA method based on hierarchical clustering is superior in simultaneous reconstruction and denoising，and the accuracy of the obtained singular values is verified.The method can avoid the manual selection of the effective number of singular values and reduce the workload of seismic data processing，which is of practical significance for the reconstruction and denoising of largescale seismic data.Keywords：multichannel singular spectrum analysis，reconstruction，denoising，hierarchical clustering，singular value analysis曹静杰，许昌昊，朱跃飞.基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法 J.石油地球物理勘探，2023，58（4）：818829.CAO Jingjie，XU Changhao，ZHU Yuefei.Simultaneous reconstruction and denoising of seismic data using multichannel singular spectrum analysis based on hierarchical clusteringJ.Oil Geophysical Prospecting，2023，58（4）：818829.0引言地震勘探数据处理的最终目标是获得高信噪比、高保真度、高分辨率的反映地下结构的数据，而提高地震数据信噪比是获得高保真和高分辨率地震图像的基础。在地震数据采集过程中，由于环境及坏道等因素，经常会出现地震道缺失及随机噪声干扰等问题，降低了地震数据的信噪比。地震数据重建与去噪是获得高信噪比完备地震数据的主要手段，对于地震数据处理的后续环节，如AVO分析1、地震属性分析、波动方程偏移成像2 等，都有非常重要的影响。地震数据重建与去噪方法可分为四类。最早应用的方法是基于波动方程，通过速度模型正演地震道记录，代替缺失的地震数据34。这类方法的重建精度依赖于速度的准确性，且计算量大，难以在实际应用中普及。第二类是基于稀疏变换的方法。该类方法不依赖地下介质速度等信息，以地震数据在某些变换域的稀疏性为先验约束，构建反演模型并求解得到重建的数据。压缩感知理论出现后，该类方法迅速应用于地震勘探领域。常用的稀疏变换包括 Fourier变换5、Radon 变换67、Curvelet 变换810 及 Seislet 变换1112 等。Hennenfent等8 提出一种基于非规则采样Curvelet变换和阈值法对地震数据进行重建与去噪；仝中飞13 将 Curvelet 变换与阈值迭代法相结合，提出Curvelet变换阈值迭代法，把地震数据随机噪声衰减和地震缺失道重建问题描述为L1模最优化问题，利用阈值迭代法求解以达到去噪与重建的目的；Cao等14 将低冗余曲波变换引入地震数据重建，大幅度提高了曲波变换在数据重建中的计算速度；葛子建15 将Fourier变换法与凸集投影法相结合，建立了一套有效的针对不规则缺失、含假频地震数据和含噪数据的重建方法；曹静杰等 16 由迭代阈值法推导出加权凸集投影法，且证明它是解无约束优化问题的一种有效方法，加权因子可看作拟合误差项系数，进而提出一种改进的凸集投影法，与原始凸集投影法相比不需增加任何计算量，只需通过选取阈值进行重建与去噪。由于采用的变换不同，因此各种方法的特点也不同。基于Fourier变换的方法计算效率高，但基于线性同相轴假设，需将数据分块；基于Radon变换的方法需采用反演方法求解Radon变换系数，计算效率低；基于Curvelet变换的方法，冗余度高，计算效率低，但重建质量高。第三类是基于矩阵降秩的方法17。矩阵降秩方法是根据相邻地震道之间的相干性压制随机噪声，减少去噪时有效信号的损失，且通过迭代重建方法补全由于空间采样不均匀造成的地震道缺失1819。在缺失的地震数据矩阵中，缺失的列会导致解不稳定，一般819石 油 地 球 物 理 勘 探2023 年需将缺失地震数据排成 Hankel矩阵，再通过对 Hankel矩阵降秩实现数据重建与噪声消除。该类方法基于奇异谱分解（SVD），最早应用于核磁共振图像处理，称为 Cadzow 滤波20，后被 Trickett21 引入地震数据处理领域，并发展成多道奇异谱分析（MSSA）方法18，22。黄建平等23 提出一种基于SVD的联合去噪、规则化方法，迭代时在地震道缺失位置进行重建，在不缺失地震道位置压制噪声。Huang等24 将阻尼因子引入 MSSA 去噪方法，提出阻尼多道奇异谱分析（DMSSA）方法，以阻尼因子减弱随机噪声带来的奇异值的增量，将信号子空间与噪声子空间进行更彻底的分离，但并未讨论地震数据重建问题。马继涛等25基于频率域 SVD矩阵降秩运算，利用凸集投影迭代方法，实现了地震数据去噪和重建的同步处理。Chen等26将阻尼算子引入 5D 地震数据去噪和重建，使矩阵降秩方法在地震数据重建中更好地压制了随机噪声，但并未对奇异值数量的选取进行相关研究。矩阵降秩方法同样基于线性同相轴假设，一般需将数据分块，当重排的Hankel矩阵维数较大时，SVD需耗费较长时间，且须首先选取每个数据块的奇异值个数。朱跃飞等27 探讨了自动确定奇异值个数的方法，利用一种Akaike信息准则自动地确定地震信号的奇异值个数，然后基于 DMSSA 方法的框架去噪，避免了通过人工选取每个数据块的奇异值个数。第四类是基于人工智能的方法。由于近年来计算机硬件的不断更新，人工智能得到快速发展，特别是深度学习理论在各领域都产生了极大影响。字典学习方法2829 是其中的一个分支，该方法与基于稀疏变换法的区别在于，基函数不固定和对复杂构造有更好的表达能力，但其计算量大、泛化能力较差，难以实施大规模数据计算。深度学习方法，如卷积神经网络30、卷积自编码器31、残差学习32 等，在地震数据去噪、重建、属性分析、断层解释、初值拾取等领域，已得到一些应用。针对地震数据重建，江金生33 提出一种多层、多模块的卷积自编码器，改进后的卷积自编码器可更有效地提取地震数据特征，训练后的卷积自编码器深度神经网络可有效地完成地震数据中随机噪声、线性噪声、面波等的压制及地震数据重建。Liu 等34 提出基于部分卷积方法提高重建结果中的模糊问题，使同相轴更清晰。Siahkoohi等35 通过炮检互易原理，利用共检波点频率波数域随机缺失数据训练卷积网络，可实现很大比例缺失数据的重建。王峰36 在地震数据去噪中采用 DnCNN模型，让卷积神经网络学习数据中噪声的特征，建立统一用于二维和三维地震数据、同时去噪和重建的模型。除了卷积网络，Wang 等37将 ResNet 用于地震数据重建，生成对抗网络 GAN38和 Unet网络3940也被用于数据重建。MSSA法可通过矩阵降秩法有效压制噪声，并重建缺失地震道，其依据为随机噪声的存在会引起地震数据在频率域矩阵秩的增加，通过降秩法压制随机噪声。即首先将地震数据进行 Fourier变换，再将每一频率分量变为块 Hankel矩阵，然后对该 Hankel矩阵进行奇异值截断，通过多次迭代达到去噪及重建的效果。如何让计算机自动确定奇异值的个数，目前尚无很好的解决方案，但这对该类方法能否具有工业应用价值至关重要。针对上述问题，本文基于DMSSA的地震数据重建与去噪框架，提出一种层次聚类法自动确定奇异值的个数，即通过聚类法将奇异值分为两类，设定其中一类对应有效信号，另一类对应噪声，且在有效信号频率范围内取最大值以确定奇异值个数，从而实现地震数据同时重建和去噪。数值实验表明，该方法对模拟数据和实际数据均能取得较好效果。1方法原理1.1MSSA去噪原理设大小为M N O的三维含噪地震数据为D(x，y，t)，且x=x1，x2，xM；y=y1，y2，yN；t=t1，t2，tO。通过以下五步实现MSSA法去噪。第一步，对D(x，y，t)做离散 Fourier变换，使地震数据从时间域变换到频率域，即F(x，y，)，其中=1，2，J，J为信号的频率个数。第二步，对F(x，y，)每一频率分量做Hankel变换。显然F(x，y，i)=F()1，1F()1，2F()1，NF()2，1F()2，2F()2，NF()M，1F()M，2F()M，N（1）是F(x，y，)中频率为i的分量。将F(x，y，i)的820第 58 卷 第 4 期曹静杰，等：基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法每行构成Hankel矩阵Aq=F()q，1F()q，2F()q，KyF()q，2F()q，3F()q，Ky+1F()q，KxF()q，Kx+1F()q，N q=1，2，M（2）式中：Aq为F(x，y，i)第q行构成的 Hankel 矩阵；Kx=N2+1；Ky=N-Kx+1；符号 表示下取整运算。将所有Aq排列成一个块Hankel矩阵 Hi=A1A2ALyA2A3ALy+1ALxALx+1AM（3）式中：Lx=M2+1；Ly=M-Lx+1。第三步，对Hi进行SVD，得到Hi=Ui()000VHi（4）式中：Ui为Lx Kx阶酉矩阵；VHi为Ly Ky阶酉矩阵；=diag(1，2，s)为对角矩阵（s为非零奇异值个数），且1 2 s。第四步，对Hi做降秩处理，确定中合理的奇异值个数。假设地震信号中有 P个不同斜率的线性同相轴，对 中的奇异值做截断处理，保留前 P个奇异值，其余奇异值设为0，可达到降秩的效果。其过程可表示为=diag(1，2，P，0，0)（5）HMSSAi=(Ui1，Ui2)()000()Vi1H()Vi2H（6）第五步，将HMSSAi进行块Hankel矩阵反变换到频率域，再通过每个频率的 Fourier反变换完成地震数据的去噪。该过程如以下两式所示F(x，y，i)=-1(HMSSAi)（7）D(x，y，t)=-1FMSSA(x，y，)（8）式中：-1表示 Hankel矩阵变换的逆；-1为 Fourier反变换算子；FMSSA(x，y，)为F(x，y，i)的集合；D(x，y，t)为输出的去噪结果。1.2DMSSA去噪原理MSSA只能将数据分解为噪声子空间和信号加噪声子空间，并不能完全去除随机噪声。当地震数据中只有有效信号时，中包含与地震图像中线性同相轴数量相同个数的非零奇异值。当地震数据混入随机噪声，所有的奇异值大小都会发生改变，非零奇异值个数将会增加。传统的 MSSA 方法去噪保留前 P个奇异值，其余奇异值置0，但保留的前P个奇异值中依然残留随机噪声引起的变化量，其去噪结果中依然残存随机噪声。Huang等24 提出的DMSSA方法，将阻尼因子引入MSSA去噪方法，通过阻尼因子减弱随机噪声造成的奇异值增量，有效压制残存的随机噪声。式（5）和式（6）被替换为下列两式T=I-DP+1()-D（9）HDMSSAi=Ui1T(Vi1)H（10）式中：I是单位矩阵；D表示阻尼因子；T表示阻尼算子；P为预估的奇异值个数。通过式（9）和式（10）两式可利用阻尼因子D对第P+1个奇异值进行放大或缩小，用前P个奇异值与其相减，达到压制前P个奇异值中与随机噪声相关部分的目的。值得一提的是，D值越小，阻尼效果越强，其去噪效果越好，但会对有效信号造成轻微损害。当D 时，式（10）退化为式（6）。根据地震数据的复杂度及噪声的强弱，D的选取范围通常是15。本文研究重点在估计有效信号对应的奇异值，未讨论阻尼因子的选择，但经过实验选择阻尼因子为2，该阻尼因子能得821石 油 地 球 物 理 勘 探2023 年到可靠的结果。将上述过程用DMSSA表示。1.3DMSSA同时重建和去噪原理缺失的地震道通常在地震数据中以 0的形式存在，在对地震数据进行重建时，地震道缺失部分可看作有效信号与随机噪声相加为 0的情况，因此三维地震数据重建法的步骤与去噪原理的步骤相似，只是在重建空白地震道过程中采用了类似加权凸集投影的方式。具体重建流程如表1。表 1 DMSSA同时重建和去噪流程for =1：max for n=1：nmax Fn()=anF()+(1-an)SDMSSAFn-1()+(I-S)DMSSAFn-1()if Fn()-Fn-1()2F return Fn()endend表 1中：F0()=F()，an是一个随循环线性减少的标量，a1=1，anmax=0，S为采样算子。上述方法的含义是，对于频率1max范围内的每个频率，设定最大迭代次数为nmax，每次迭代都执行表 1 中Fn()的计算过程，如果迭代次数未达到nmax，但 是 满 足Fn()-Fn-1()2F，则 返 回Fn()为数据结果，将这些频率范围都执行上述运算，就得到频率域的同时重建和去噪结果，再通过Fourier反变换就能得到时间域的重建与去噪结果。该方法通过循环迭代的方式完成地震数据的重建与去噪，每次迭代中，在地震道缺失位置处由ADMSSA方法进行重建，在地震道没有缺失的位置由anF()+(1-an)SDMSSAFn-1()对数据进行去噪，F()为原始观测数据，Fn-1()为前一轮奇异谱分析后得到的数据，an控制两者所占的比率，由1线性递减至0。该方法可有效降低直接重建时由于原始观测数据的噪声对数据相干性造成的影响，实现地震数据的同时重建和去噪。1.4层次聚类法自动确定奇异值个数DMSSA类方法实现同时重建和去噪时，奇异值个数P的选取至关重要。若奇异值个数选取过大，则重建后残留的噪声过大。若奇异值个数选取过少，则会损害有效信号。在计算机运算过程中，由于实际地震资料的数据量过大，常需对地震数据进行分块处理；而每块的同相轴数量并不相同，因此需要人工确定 P，无法实现自动化重建与去噪。本节通过对奇异值的分布特点进行分析，采用层次聚类方法使计算机自动确定合适的奇异值个数 P，进而更好地完成地震数据的重建。图 1a是采用主频为 30 Hz的雷克子波正演得到一块大小为40 40300的模拟数据，时间采样率为2 ms，该数据有 3 个线性同相轴。图 1b 为随机抽取50%地震道且增加 10%随机噪声后的数据，信噪比为3.9 dB。图2为频率为40 Hz时块Hankel矩阵的奇异值分布，其中绿线为图 1a数据对应的奇异值，红线为待重建地震数据（图 1b）对应的奇异值。由图 2中看出，图1a地震数据第4项及之后的奇异值远小于前 3项奇异值，且绝大多数为 0。随机噪声的增加与地震道的缺失使奇异值序列中前 3项奇异值变小，非零项增多，但前 3项奇异值依然与 3项之后的奇异值有较大区别。上述含噪声地震数据的奇异值可分为两类：一类是与有效信号有关的奇异值，数值较大；另一类是与噪声和道缺失有关的奇异值，数值较小。本文引入一种层次聚类的方法，将多个元素按照相似程度自下而上的进行分类，从而确定奇异值个数。定义类与类之间的距离Dpq=min|dijzi Gp，zj Gq（11）式中：Gp、Gq为两个不同的类；zi和zj均为二维向量，分别属于Gp和Gq；dij为zi与zj之间的欧式距离。将式（4）中的奇异值排列成为二维数组(1，1)；(2，2)；(n，s)（12）使用(z1，z2，zs)进行表示。层次聚类的具体过程如下：（1）将(z1，z2，zs)中每个数组各自分为一类，计算各类之间的距离；（2）合并类间距最小的两个类，构成一个新类；（3）计算新类与当前各类的距离；（4）若 类 的 个 数 为 2，即终止计算；否 则，返 回822第 58 卷 第 4 期曹静杰，等：基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法步骤（2）。上述过程将奇异值分为两类，由于奇异值个数P远远小于奇异值的总数量s，故选取数量较少的一类，计算其包含奇异值的数量即可得到P。研究发现，使用DMSSA方法时只需对有效信号频率范围进行聚类分析，对有效频率范围以外可取值为零。在有效频率范围内，每个频率分量的聚类结果也有可能发生变化，为了使该方法结果稳定，选择有效信号频率范围内聚类，以结果的最大值作为最终的奇异值个数。为验证上述观点，图3给出了图1数据1150 Hz层次聚类法测得的奇异值数量。对图 1 数据而言，1090 Hz在信号的有效频率范围内，层次聚类法确定的奇异值个数稳定在3；超出该范围时测得的奇异值数量不断上下波动，其原因是超出该频率范围的为噪声。根据MSSA理论，可验证在1090 Hz内测得的奇异值数量是准确的。2模拟数据实验首先用模拟数据验证方法的有效性。图 4a为采用层次聚类法确定的奇异值个数P=3时图1b的重建结果，阻尼因子 D 取 2，重建后信噪比为 19.9 dB；图4b为重建误差。从图4可看出，重建后线性同相轴清晰可见，噪声残留较少，说明该方法预测的 P值准确，能重建出较高信噪比的地震图像。为了验证方法的稳定性，下面采用另一组模拟数据进行测试。图 5a是对主频为 30 Hz的雷克子波合成得到的一块维数为40 40 300的模拟地震数图 1模拟数据一（a）完整数据；（b）加噪且缺失的数据图 2图 1a和图 1b频率为 40 Hz时的奇异值分布图 4模拟数据一重建效果（a）模拟数据一层次聚类法重建与去噪结果；（b）图 1a与图 4a误差图 3模拟数据一层次聚类的结果823石 油 地 球 物 理 勘 探2023 年据，时间采样率为 2 ms，该数据有 5个线性同相轴。图 5b为随机抽取 50%地震道且增加 10%随机噪声后的数据，其信噪比为0.7 dB。图6为通过层次聚类测得的奇异值个数，可见在1090 Hz频率范围内，测得的P值为5。图7a为模拟数据二基于层次聚类法的重建结果，阻尼因子取2，重建后信噪比为18.6 dB，图7b为重建误差。图4与图7两组模拟地震数据的重建效果都相对较好，线性同相轴恢复清晰，噪声残留较少，证明了本文方法预测奇异值数量的可行性与准确性。3实际地震数据实验为了验证本文方法对实际数据的处理效果，选用了两块三维叠后地震数据。图8a是一块大小为40 40 276的三维实际地震数据，时间上有276个采样点，时间采样率为2 ms。主测线方向与联络测线方向上各有 40个采样点。对其加入 20%的随机噪声（图 8b），然后进行随机 30%地震道抽稀（图8c）。可见该数据噪声能量较强，信噪比低，地震道缺失严重。在计算大规模实际地震数据去噪和重建时需做分块处理。通过滑动分块将大块地震数据分成各个小块，对每个小块通过DMSSA方法进行重建。通过多次比较，发现空间分块采用10 10左右时能取得稳定的效果。图 9为对其中一个小块通过层次聚类法确定奇异值的结果。可见实际数据相比于模拟数据，有较大的差异性，奇异值波动更大；超出有效频率（1090 Hz）范围后，测得的奇异值个数出现严重波动；在有效频率范围内，奇异值测得的数量相对较为稳定，但与模拟地震数据不同的是，在有效频率范围内测得的奇异值数量并不完全一样。如何选择恰当的奇异值需进行细致探讨。若选取的奇异值过大，则会影响去噪效果；若选取的奇异值过小，则会损害有效信号。因此，本文对各频率测得的奇异值个数进行统计（表2），并按测得数量从大到小排列，取其前面 90%，选取其中最大值作为该小块最终确认的奇异值，这样选取奇异值可避免损伤大部分有效信号，且不会对去噪效果产生较大影响。对图 9在有效频率范围内（1090 Hz）的各频率所测奇异值的个数进行统计（表2）。可见测得的奇异值个数为 1、2、3的频率的累计总数为 76，占所有频率个数的 93.76%，因此取奇异值个数 1、2、3中的最大图 5模拟数据二（a）完整数据；（b）加噪且缺失的数据图 6模拟数据二层次聚类的结果图 7模拟数据二重建效果（a）模拟数据二基于层次聚类法的重建结果；（b）重建的误差824第 58 卷 第 4 期曹静杰，等：基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法值 3可避免对有效信号造成损伤，且会得到较好的去噪效果，显然此小块在有效频率范围内测得的奇异值的个数为3。对图8实际地震数据采用滑动分块重建方法，每小块大小为10 10，在主测线方向与联络测线方向均取步长为2，共分得256块。图10a为基于层次聚类法确定奇异值的重建结果，阻尼因子取2，其每块测得的奇异值数量如图 11所示。图 10b为去噪后残留的噪声。可见该实际数据的重建与去噪效果（图10a）良好，有效信号损失较少，剖面的同相轴轮廓恢复清晰，噪声清除干净，证明基于本文方法在实际地震数据中的有效性与可行性。对图8实际地震重建中，由图11可见大部分数据块中测得的奇异值在 13。为了验证层次聚类法所测奇异值的准确性，指定所有数据块的奇异值全部为图 10实际数据一重建效果（a）实际数据一同时重建和去噪结果；（b）图 10a与图 8b的差图9实际数据一某一分块基于层次聚类法确定的奇异值个数表 2图 9中测得奇异值个数统计奇异值12345610奇异值个数累计541752111比例/%66.6720.996.172.471.231.231.23图 8三维实际地震数据一及其加噪和缺失数据（a）三维原始实际数据；（b）增加 20%随机噪声后数据；（c）30%随机缺失且加噪数据825石 油 地 球 物 理 勘 探2023 年1、2、3，所得到的重建结果与本文每个数据块都自动确定奇异值个数的方法所得重建结果进行对比。从图12可见，当对每个数据块的奇异值直接取1时，其图像左下方的断层变得不再明显。对每个数据块的奇异值取 2和 3时有较好重建效果，但仍然残留少量噪声，且存在线性同相轴恢复不完整的情形。采用层次聚类确定奇异值的方法每个数据块自动确定适合的奇异值个数，重建图像的噪声去除干净，线性同相轴清晰可见，图像重建效果良好。图 13a 为三维实际地震数据实例二，大小为150 150 276，时间上有 276个采样点，时间采样率为2 ms，主测线方向与联络测线方向各有150个采样点。对其加入20%随机噪声（图13b），再对其随机图 11对所有分块测得的奇异值图 12指定奇异值与通过层次聚类法取确定奇异值重建结果对比（a）、（b）、（c）对应每个数据块的奇异值都取 1、2、3的效果图；（d）本文方法的重建效果；（e）、（f）、（g）、（h）为重建数据与原数据的残差。826第 58 卷 第 4 期曹静杰，等：基于层次聚类多道奇异谱分析的地震数据同时重建与去噪方法抽稀 30%地震道（图 13c）。该数据同样存在噪声能量强、信噪比低、地震道缺失严重等特点。图14为针对图13基于层次聚类法确定奇异值的重建结果，此时阻尼因子取值为 2，得到去噪后的残留噪声（图 14b）。可见该重建与去噪结果中，同相轴轮廓恢复清晰，噪声清除干净，有效信号损失较少。再次证明了本文方法对实际地震数据的有效性与可行性。4讨论本文方法需要对地震信号有效频率进行分析，选定有效信号的范围。对地震数据进行分块处理和拼接处理也是必须的，分块的大小对本文方法有一定影响。经验表明：当空间分块大约为10 10时，能取得稳定的结果；拼接方法的选择对于重建与去噪结果同样会产生影响；另外，阻尼因子也应根据噪声的强弱程度而选取。5结论矩阵降秩方法是地震数据去噪与重建的一种重要方法，秩的选择对于矩阵降秩效果的影响巨大。当秩选取过小时，将损害有效信息；但当秩选取过大时，会在去噪及重建后留下明显的残留噪声。面对海量地震数据，人工选取奇异值的方法具有一定的主观性，且不利于去噪和重建的工业化实现。本文通过对缺道且含有噪声的地震数据的奇异值进行分析，引入了层次聚类方法，将奇异值自动分为有效信号和噪声两类，选择在有效频带范围内实现聚类，使得多道奇异谱分析方法能自适应且稳定地确定奇异值的个数，并通过模拟数据与实际数据验证了方法的优越性。实验表明，利用层次聚类法确定的奇异值数量可靠稳定，重建后线性同相轴清晰可见，无太多残留噪声，地震数据重建效果良好。图 13三维实际地震数据二及其缺失的数据（a）三维实际地震数据二；（b）增加 20%随机噪声后的数据；（c）30%随机缺失且加噪的数据图 14三维实际地震数据二重建效果（a）实际数据二同时重建和去噪结果；（b）图 14a与图 13b的差827石 油 地 球 物 理 勘 探2023 年参 考 文 献1 付欣,张小龙,孙云强,等.基于双差反演策略的时移地震 AVO 反演J.石油地球物理勘探,2022,57(4):888896.FU Xin,ZHANG Xiaolong,SUN Yunqiang,et al.Timelapse seismic AVO inversion based on a doubledifference inversion strategyJ.Oil Geophysical Prospecting,2022,57(4):888896.2 陈生昌,刘亚楠,李代光.基于波动方程偏移的地震数据地层成像J.石油地球物理勘探,2022,57(5):11051113.CHEN Shengchang,LIU Yanan,LI Daiguang.Stratigraphic imaging with seismic data based on wave equation migrationJ.Oil Geophysical Prospecting,2022,57(5):11051113.3 RONEN J.Waveequation trace interpolationJ.Geophysics,1987,52(7):9739844 CANNING A,GARDNER G.Regularizing 3D data sets with DMOJ.Geophysics,1996,61(4):11031114.5 SACCHI M D,ULRYCH T J.Highresolution velo city gather and offset space reconstructionJ.Geophy sics,1995,60(4):11691177.6 薛亚茹,唐欢欢,陈小宏.高阶高分辨率 Radon变换地震数据重建方法J.石油地球物理勘探,2014,49(1):95100,131.XUE Yaru,TANG Huanhuan,CHEN Xiaohong.Seismic data reconstruction based on high order high resolution Radon transformJ.Oil Geophysical Prospec ting,2014,49(1):95100,131.7 唐欢欢,毛伟建,詹毅.3D 高阶抛物 Radon 变换在不规则地震数据保幅重建中的应用J.地球物理学报,2020,63(9):3