0 引言

地震勘探数据处理的最终目标是获得高信噪比、高保真度、高分辨率的反映地下结构的数据，而提高地震数据信噪比是获得高保真和高分辨率地震图像的基础。在地震数据采集过程中，由于环境及坏道等因素，经常会出现地震道缺失及随机噪声干扰等问题，降低了地震数据的信噪比。地震数据重建与去噪是获得高信噪比完备地震数据的主要手段，对于地震数据处理的后续环节，如AVO分析^[1]、地震属性分析、波动方程偏移成像^[2]等，都有非常重要的影响。

地震数据重建与去噪方法可分为四类。最早应用的方法是基于波动方程，通过速度模型正演地震道记录，代替缺失的地震数据^[3-4]。这类方法的重建精度依赖于速度的准确性，且计算量大，难以在实际应用中普及。

第二类是基于稀疏变换的方法。该类方法不依赖地下介质速度等信息，以地震数据在某些变换域的稀疏性为先验约束，构建反演模型并求解得到重建的数据。压缩感知理论出现后，该类方法迅速应用于地震勘探领域。常用的稀疏变换包括Fourier变换^[5]、Radon变换^[6-7]、Curvelet变换^[8-10]及Seislet变换^[11-12]等。Hennenfent等^[8]提出一种基于非规则采样Curvelet变换和阈值法对地震数据进行重建与去噪；仝中飞^[13]将Curvelet变换与阈值迭代法相结合，提出Curvelet变换阈值迭代法，把地震数据随机噪声衰减和地震缺失道重建问题描述为L₁模最优化问题，利用阈值迭代法求解以达到去噪与重建的目的；Cao等^[14]将低冗余曲波变换引入地震数据重建，大幅度提高了曲波变换在数据重建中的计算速度；葛子建^[15]将Fourier变换法与凸集投影法相结合，建立了一套有效的针对不规则缺失、含假频地震数据和含噪数据的重建方法；曹静杰等^[16]由迭代阈值法推导出加权凸集投影法，且证明它是解无约束优化问题的一种有效方法，加权因子可看作拟合误差项系数，进而提出一种改进的凸集投影法，与原始凸集投影法相比不需增加任何计算量，只需通过选取阈值进行重建与去噪。由于采用的变换不同，因此各种方法的特点也不同。基于Fourier变换的方法计算效率高，但基于线性同相轴假设，需将数据分块；基于Radon变换的方法需采用反演方法求解Radon变换系数，计算效率低；基于Curvelet变换的方法，冗余度高，计算效率低，但重建质量高。

第三类是基于矩阵降秩的方法^[17]。矩阵降秩方法是根据相邻地震道之间的相干性压制随机噪声，减少去噪时有效信号的损失，且通过迭代重建方法补全由于空间采样不均匀造成的地震道缺失^[18-19]。在缺失的地震数据矩阵中，缺失的列会导致解不稳定，一般需将缺失地震数据排成Hankel矩阵，再通过对Hankel矩阵降秩实现数据重建与噪声消除。该类方法基于奇异谱分解（SVD），最早应用于核磁共振图像处理，称为Cadzow滤波^[20]，后被Trickett^[21]引入地震数据处理领域，并发展成多道奇异谱分析（MSSA）方法^{[18, 22]}。黄建平等^[23]提出一种基于SVD的联合去噪、规则化方法，迭代时在地震道缺失位置进行重建，在不缺失地震道位置压制噪声。Huang等^[24]将阻尼因子引入MSSA去噪方法，提出阻尼多道奇异谱分析（DMSSA）方法，以阻尼因子减弱随机噪声带来的奇异值的增量，将信号子空间与噪声子空间进行更彻底的分离，但并未讨论地震数据重建问题。马继涛等^[25]基于频率域SVD矩阵降秩运算，利用凸集投影迭代方法，实现了地震数据去噪和重建的同步处理。Chen等^[26]将阻尼算子引入5D地震数据去噪和重建，使矩阵降秩方法在地震数据重建中更好地压制了随机噪声，但并未对奇异值数量的选取进行相关研究。

矩阵降秩方法同样基于线性同相轴假设，一般需将数据分块，当重排的Hankel矩阵维数较大时，SVD需耗费较长时间，且须首先选取每个数据块的奇异值个数。朱跃飞等^[27]探讨了自动确定奇异值个数的方法，利用一种Akaike信息准则自动地确定地震信号的奇异值个数，然后基于DMSSA方法的框架去噪，避免了通过人工选取每个数据块的奇异值个数。

第四类是基于人工智能的方法。由于近年来计算机硬件的不断更新，人工智能得到快速发展，特别是深度学习理论在各领域都产生了极大影响。字典学习方法^[28-29]是其中的一个分支，该方法与基于稀疏变换法的区别在于，基函数不固定和对复杂构造有更好的表达能力，但其计算量大、泛化能力较差，难以实施大规模数据计算。深度学习方法，如卷积神经网络^[30]、卷积自编码器^[31]、残差学习^[32]等，在地震数据去噪、重建、属性分析、断层解释、初值拾取等领域，已得到一些应用。针对地震数据重建，江金生^[33]提出一种多层、多模块的卷积自编码器，改进后的卷积自编码器可更有效地提取地震数据特征，训练后的卷积自编码器深度神经网络可有效地完成地震数据中随机噪声、线性噪声、面波等的压制及地震数据重建。Liu等^[34]提出基于部分卷积方法提高重建结果中的模糊问题，使同相轴更清晰。Siahkoohi等^[35]通过炮—检互易原理，利用共检波点频率—波数域随机缺失数据训练卷积网络，可实现很大比例缺失数据的重建。王峰^[36]在地震数据去噪中采用DnCNN模型，让卷积神经网络学习数据中噪声的特征，建立统一用于二维和三维地震数据、同时去噪和重建的模型。除了卷积网络，Wang等^[37]将ResNet用于地震数据重建，生成对抗网络GAN^[38]和U-net网络^[39-40]也被用于数据重建。

MSSA法可通过矩阵降秩法有效压制噪声，并重建缺失地震道，其依据为随机噪声的存在会引起地震数据在频率域矩阵秩的增加，通过降秩法压制随机噪声。即首先将地震数据进行Fourier变换，再将每一频率分量变为块Hankel矩阵，然后对该Hankel矩阵进行奇异值截断，通过多次迭代达到去噪及重建的效果。如何让计算机自动确定奇异值的个数，目前尚无很好的解决方案，但这对该类方法能否具有工业应用价值至关重要。

针对上述问题，本文基于DMSSA的地震数据重建与去噪框架，提出一种层次聚类法自动确定奇异值的个数，即通过聚类法将奇异值分为两类，设定其中一类对应有效信号，另一类对应噪声，且在有效信号频率范围内取最大值以确定奇异值个数，从而实现地震数据同时重建和去噪。数值实验表明，该方法对模拟数据和实际数据均能取得较好效果。

1 方法原理 1.1 MSSA去噪原理

设大小为$ M\times N\times O $的三维含噪地震数据为$ {\boldsymbol{D}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{t}}) $，且$ {\boldsymbol{x}}={x}_{1}, {x}_{2}, \cdots , {x}_{M} $；$ {\boldsymbol{y}}={y}_{1}, {y}_{2}, \cdots , {y}_{N} $；$ {\boldsymbol{t}}={t}_{1}, {t}_{2}, \cdots , {t}_{O} $。通过以下五步实现MSSA法去噪。

第一步，对$ {\boldsymbol{D}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{t}}) $做离散Fourier变换，使地震数据从时间域变换到频率域，即$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{\omega}}) $，其中$ {\boldsymbol{\omega}}={\omega }_{1}, {\omega }_{2}, \cdots , {\omega }_{J} $，J为信号的频率个数。

第二步，对$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{\omega}}) $每一频率分量做Hankel变换。显然

$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\omega }_{i})= \\ \;\;\;\;\left(\begin{array}{cccc}F\left(\mathrm{1, 1}\right)& F\left(\mathrm{1, 2}\right)& \cdots & F\left(1, N\right)\\ F\left(\mathrm{2, 1}\right)& F\left(\mathrm{2, 2}\right)& \cdots & F\left(2, N\right)\\ & & ⋮& \\ F\left(M, 1\right)& F\left(M, 2\right)& \cdots & F\left(M, N\right)\end{array}\right) $

(1)

是$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{\omega}}) $中频率为$ {\omega }_{i}^{} $的分量。将$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\omega }_{i}) $的每行构成Hankel矩阵

$ {{\boldsymbol{A}}}_{q}=\left(\begin{array}{cccc}F\left(q, 1\right)& F\left(q, 2\right)& \cdots & F\left(q, {K}_{y}\right)\\ F\left(q, 2\right)& F\left(q, 3\right)& \cdots & F\left(q, {K}_{y+1}\right)\\ & & ⋮& \\ F\left(q, {K}_{x}\right)& F\left(q, {K}_{x+1}\right)& \cdots & F\left(q, N\right)\end{array}\right)q=\mathrm{1, 2}, \cdots , M $

(2)

式中：$ {{\boldsymbol{A}}}_{q} $为$ {\boldsymbol{F}}({\boldsymbol{x}}, {\boldsymbol{y}}, {\omega }_{i}) $第$ q $行构成的Hankel矩阵；$ {K}_{x}=⌊\frac{N}{2}⌋+1 $；$ {K}_{y}=N-{K}_{x}+1 $；符号$ ⌊\cdot ⌋ $表示下取整运算。将所有$ {{\boldsymbol{A}}}_{q} $排列成一个块Hankel矩阵

$ {{\boldsymbol{H}}}_{{\omega }_{i}}=\left(\begin{array}{cccc}{{\boldsymbol{A}}}_{1}& {{\boldsymbol{A}}}_{2}& \cdots & {{\boldsymbol{A}}}_{{L}_{y}}\\ {{\boldsymbol{A}}}_{2}& {{\boldsymbol{A}}}_{3}& \cdots & {{\boldsymbol{A}}}_{{L}_{y}+1}\\ & & ⋮& \\ {{\boldsymbol{A}}}_{{L}_{x}}& {{\boldsymbol{A}}}_{{L}_{x}+1}& \cdots & {{\boldsymbol{A}}}_{M}\end{array}\right) $

(3)

式中：$ {L}_{x}=⌊\frac{M}{2}⌋+1 $；$ {L}_{y}=M-{L}_{x}+1 $。

第三步，对$ {{\boldsymbol{H}}}_{{\omega }_{i}} $进行SVD，得到

$ {{\boldsymbol{H}}}_{{\omega }_{i}}={{\boldsymbol{U}}}_{{\omega }_{i}}\left(\begin{array}{cc}{\boldsymbol{\varSigma}}& 0\\ 0& 0\end{array}\right){{\boldsymbol{V}}}_{{\omega }_{i}}^{\mathrm{H}} $

(4)

式中：$ {{\boldsymbol{U}}}_{{\omega }_{i}} $为$ {L}_{x}\times {K}_{x} $阶酉矩阵；$ {{\boldsymbol{V}}}_{{\omega }_{i}}^{\mathrm{H}} $为$ {L}_{y}\times {K}_{y} $阶酉矩阵；$ {\boldsymbol{\varSigma}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\sigma }_{1}, {\sigma }_{2}, \cdots , {\sigma }_{s}\right) $为对角矩阵（s为非零奇异值个数），且$ {\sigma }_{1} > {\sigma }_{2} > \cdots > {\sigma }_{s} $。

第四步，对$ {{\boldsymbol{H}}}_{{\omega }_{i}} $做降秩处理，确定Σ中合理的奇异值个数。假设地震信号中有P个不同斜率的线性同相轴，对Σ中的奇异值做截断处理，保留前P个奇异值，其余奇异值设为0，可达到降秩的效果。其过程可表示为

$ {{\boldsymbol{\varSigma}}}^{\text{'}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\sigma }_{1}, {\sigma }_{2}, \cdots , {\sigma }_{P}, 0, \cdots , 0\right) $

(5)

$ {{\boldsymbol{H}}}_{{\omega }_{i}}^{\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}}=\left({{\boldsymbol{U}}}_{1}^{{\omega }_{i}}, {{\boldsymbol{U}}}_{2}^{{\omega }_{i}}\right)\left(\begin{array}{cc}{{\boldsymbol{\varSigma}}}^{\text{'}}& 0\\ 0& 0\end{array}\right)\left[\begin{array}{c}{\left({{\boldsymbol{V}}}_{1}^{{\omega }_{i}}\right)}^{\mathrm{H}}\\ {\left({{\boldsymbol{V}}}_{2}^{{\omega }_{i}}\right)}^{\mathrm{H}}\end{array}\right] $

(6)

第五步，将$ {{\boldsymbol{H}}}_{{\omega }_{i}}^{\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}} $进行块Hankel矩阵反变换到频率域，再通过每个频率的Fourier反变换完成地震数据的去噪。该过程如以下两式所示

$ \widehat{{\boldsymbol{F}}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\omega }_{i}\right)={{\boldsymbol{\varGamma }}}^{-1}\left({{\boldsymbol{H}}}_{{\omega }_{i}}^{\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}}\right) $

(7)

$ \widehat{{\boldsymbol{D}}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{t}}\right)={{\boldsymbol{\varPsi }}}^{-1}\left[{{\boldsymbol{F}}}^{\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{\omega}}\right)\right] $

(8)

式中：$ {{\boldsymbol{\varGamma }}}^{-1} $表示Hankel矩阵变换的逆；$ {{\boldsymbol{\varPsi }}}^{-1} $为Fourier反变换算子；$ {{\boldsymbol{F}}}^{\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{\omega}}\right) $为$ \widehat{{\boldsymbol{F}}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\omega }_{i}\right) $的集合；$ \widehat{{\boldsymbol{D}}}\left({\boldsymbol{x}}, {\boldsymbol{y}}, {\boldsymbol{t}}\right) $为输出的去噪结果。

1.2 DMSSA去噪原理

MSSA只能将数据分解为噪声子空间和信号加噪声子空间，并不能完全去除随机噪声。当地震数据中只有有效信号时，$ {{\boldsymbol{\varSigma}}}^{\text{'}} $中包含与地震图像中线性同相轴数量相同个数的非零奇异值。当地震数据混入随机噪声，所有的奇异值大小都会发生改变，非零奇异值个数将会增加。传统的MSSA方法去噪保留前P个奇异值，其余奇异值置0，但保留的前P个奇异值中依然残留随机噪声引起的变化量，其去噪结果中依然残存随机噪声。

Huang等^[24]提出的DMSSA方法，将阻尼因子引入MSSA去噪方法，通过阻尼因子减弱随机噪声造成的奇异值增量，有效压制残存的随机噪声。

式（5）和式（6）被替换为下列两式

$ {\boldsymbol{T}}={\boldsymbol{I}}-{\sigma }_{P+1}^{D}{\left({{\boldsymbol{\varSigma}}}^{\text{'}}\right)}^{-D} $

(9)

$ {{\boldsymbol{H}}}_{{\omega }_{i}}^{\mathrm{D}\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}}={{\boldsymbol{U}}}_{1}^{{\omega }_{i}}{{\boldsymbol{\varSigma}}}^{\text{'}}{\boldsymbol{T}}{\left({{\boldsymbol{V}}}_{1}^{{\omega }_{i}}\right)}^{\mathrm{H}} $

(10)

式中：I是单位矩阵；D表示阻尼因子；T表示阻尼算子；P为预估的奇异值个数。通过式（9）和式（10）两式可利用阻尼因子D对第P+1个奇异值进行放大或缩小，用前P个奇异值与其相减，达到压制前P个奇异值中与随机噪声相关部分的目的。

值得一提的是，D值越小，阻尼效果越强，其去噪效果越好，但会对有效信号造成轻微损害。当D$ \to \mathrm{\infty } $时，式（10）退化为式（6）。根据地震数据的复杂度及噪声的强弱，D的选取范围通常是1~5。本文研究重点在估计有效信号对应的奇异值，未讨论阻尼因子的选择，但经过实验选择阻尼因子为2，该阻尼因子能得到可靠的结果。将上述过程用DMSSA表示。

1.3 DMSSA同时重建和去噪原理

缺失的地震道通常在地震数据中以0的形式存在，在对地震数据进行重建时，地震道缺失部分可看作有效信号与随机噪声相加为0的情况，因此三维地震数据重建法的步骤与去噪原理的步骤相似，只是在重建空白地震道过程中采用了类似加权凸集投影的方式。具体重建流程如表 1。

表 1中：$ {{\boldsymbol{F}}}_{0}\left(\omega \right)={\boldsymbol{F}}\left(\omega \right) $，$ {a}_{n} $是一个随循环线性减少的标量，$ {a}_{1}=1 $，$ {a}_{{n}_{\mathrm{m}\mathrm{a}\mathrm{x}}}=0 $，$ \widehat{S} $为采样算子。

上述方法的含义是，对于频率$ {\omega }_{1} $~$ {\omega }_{\mathrm{m}\mathrm{a}\mathrm{x}} $范围内的每个频率$ \omega $，设定最大迭代次数为$ {n}_{\mathrm{m}\mathrm{a}\mathrm{x}} $，每次迭代都执行表 1中$ {{\boldsymbol{F}}}_{n}\left(\omega \right) $的计算过程，如果迭代次数未达到$ {n}_{\mathrm{m}\mathrm{a}\mathrm{x}} $，但是满足$ {‖{{\boldsymbol{F}}}_{n}^{}\left(\omega \right)-{{\boldsymbol{F}}}_{n-1}^{}\left(\omega \right)‖}_{F}^{2} < \varepsilon $，则返回$ {{\boldsymbol{F}}}_{n}\left(\omega \right) $为数据结果，将这些频率范围都执行上述运算，就得到频率域的同时重建和去噪结果，再通过Fourier反变换就能得到时间域的重建与去噪结果。

该方法通过循环迭代的方式完成地震数据的重建与去噪，每次迭代中，在地震道缺失位置处由ADMSSA方法进行重建，在地震道没有缺失的位置由$ {a}_{n}{\boldsymbol{F}}\left(\omega \right)+\left(1-{a}_{n}\right)\widehat{S}\cdot \mathrm{D}\mathrm{M}\mathrm{S}\mathrm{S}\mathrm{A}\left[{{\boldsymbol{F}}}_{n-1}\left(\omega \right)\right] $对数据进行去噪，$ {\boldsymbol{F}}\left(\omega \right) $为原始观测数据，$ {{\boldsymbol{F}}}_{n-1}\left(\omega \right) $为前一轮奇异谱分析后得到的数据，$ {a}_{n} $控制两者所占的比率，由1线性递减至0。该方法可有效降低直接重建时由于原始观测数据的噪声对数据相干性造成的影响，实现地震数据的同时重建和去噪。

1.4 层次聚类法自动确定奇异值个数

DMSSA类方法实现同时重建和去噪时，奇异值个数P的选取至关重要。若奇异值个数选取过大，则重建后残留的噪声过大。若奇异值个数选取过少，则会损害有效信号。在计算机运算过程中，由于实际地震资料的数据量过大，常需对地震数据进行分块处理；而每块的同相轴数量并不相同，因此需要人工确定P，无法实现自动化重建与去噪。本节通过对奇异值的分布特点进行分析，采用层次聚类方法使计算机自动确定合适的奇异值个数P，进而更好地完成地震数据的重建。

图 1a是采用主频为30 Hz的雷克子波正演得到一块大小为$ 40\times 40\times $$ 300 $的模拟数据，时间采样率为2 ms，该数据有3个线性同相轴。图 1b为随机抽取50%地震道且增加10%随机噪声后的数据，信噪比为-3.9 dB。图 2为频率为40 Hz时块Hankel矩阵的奇异值分布，其中绿线为图 1a数据对应的奇异值，红线为待重建地震数据（图 1b）对应的奇异值。由图 2中看出，图 1a地震数据第4项及之后的奇异值远小于前3项奇异值，且绝大多数为0。随机噪声的增加与地震道的缺失使奇异值序列中前3项奇异值变小，非零项增多，但前3项奇异值依然与3项之后的奇异值有较大区别。上述含噪声地震数据的奇异值可分为两类：一类是与有效信号有关的奇异值，数值较大；另一类是与噪声和道缺失有关的奇异值，数值较小。本文引入一种层次聚类的方法，将多个元素按照相似程度自下而上的进行分类，从而确定奇异值个数。

定义类与类之间的距离

$ {D}_{pq}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\left.{d}_{ij}\right|{\boldsymbol{z}}_{i}\in {{\boldsymbol{G}}}_{p}, {\boldsymbol{z}}_{j}\in {{\boldsymbol{G}}}_{q}\right\} $

(11)

式中：$ {{\boldsymbol{G}}}_{p} $、$ {{\boldsymbol{G}}}_{q} $为两个不同的类；$ {\boldsymbol{z}}_{i} $和$ {\boldsymbol{z}}_{j} $均为二维向量，分别属于$ {{\boldsymbol{G}}}_{p} $和$ {{\boldsymbol{G}}}_{q} $；$ {d}_{ij} $为$ {\boldsymbol{z}}_{i} $与$ {\boldsymbol{z}}_{j} $之间的欧式距离。

将式（4）中$ {\boldsymbol{\varSigma}} $的奇异值排列成为二维数组

$ \left[(1, {\sigma }_{1});(2, {\sigma }_{2});\cdots ;(n, {\sigma }_{s})\right] $

(12)

使用$ \left({\boldsymbol{z}}_{1}, {\boldsymbol{z}}_{2}, \cdots , {\boldsymbol{z}}_{s}\right) $进行表示。

层次聚类的具体过程如下：

（1）将$ \left({\boldsymbol{z}}_{1}, {\boldsymbol{z}}_{2}, \cdots , {\boldsymbol{z}}_{s}\right) $中每个数组各自分为一类，计算各类之间的距离；

（2）合并类间距最小的两个类，构成一个新类；

（3）计算新类与当前各类的距离；

（4）若类的个数为2，即终止计算；否则，返回步骤（2）。

上述过程将奇异值分为两类，由于奇异值个数P远远小于奇异值的总数量s，故选取数量较少的一类，计算其包含奇异值的数量即可得到P。

研究发现，使用DMSSA方法时只需对有效信号频率范围进行聚类分析，对有效频率范围以外可取值为零。在有效频率范围内，每个频率分量的聚类结果也有可能发生变化，为了使该方法结果稳定，选择有效信号频率范围内聚类，以结果的最大值作为最终的奇异值个数。

为验证上述观点，图 3给出了图 1数据1~150 Hz层次聚类法测得的奇异值数量。对图 1数据而言，10~90 Hz在信号的有效频率范围内，层次聚类法确定的奇异值个数稳定在3；超出该范围时测得的奇异值数量不断上下波动，其原因是超出该频率范围的为噪声。根据MSSA理论，可验证在10~90 Hz内测得的奇异值数量是准确的。

2 模拟数据实验

首先用模拟数据验证方法的有效性。图 4a为采用层次聚类法确定的奇异值个数P=3时图 1b的重建结果，阻尼因子D取2，重建后信噪比为19.9 dB；图 4b为重建误差。从图 4可看出，重建后线性同相轴清晰可见，噪声残留较少，说明该方法预测的P值准确，能重建出较高信噪比的地震图像。

为了验证方法的稳定性，下面采用另一组模拟数据进行测试。图 5a是对主频为30 Hz的雷克子波合成得到的一块维数为$ 40\times 40 $$ \times 300 $的模拟地震数据，时间采样率为2 ms，该数据有5个线性同相轴。图 5b为随机抽取50%地震道且增加10%随机噪声后的数据，其信噪比为-0.7 dB。

图 6为通过层次聚类测得的奇异值个数，可见在10~90 Hz频率范围内，测得的P值为5。

图 7a为模拟数据二基于层次聚类法的重建结果，阻尼因子取2，重建后信噪比为18.6 dB，图 7b为重建误差。图 4与图 7两组模拟地震数据的重建效果都相对较好，线性同相轴恢复清晰，噪声残留较少，证明了本文方法预测奇异值数量的可行性与准确性。

3 实际地震数据实验

为了验证本文方法对实际数据的处理效果，选用了两块三维叠后地震数据。

图 8a是一块大小为$ 40\times 40\times 276 $的三维实际地震数据，时间上有276个采样点，时间采样率为2 ms。主测线方向与联络测线方向上各有40个采样点。对其加入20%的随机噪声（图 8b），然后进行随机30%地震道抽稀（图 8c）。可见该数据噪声能量较强，信噪比低，地震道缺失严重。

在计算大规模实际地震数据去噪和重建时需做分块处理。通过滑动分块将大块地震数据分成各个小块，对每个小块通过DMSSA方法进行重建。通过多次比较，发现空间分块采用$ 10\times 10 $左右时能取得稳定的效果。

图 9为对其中一个小块通过层次聚类法确定奇异值的结果。可见实际数据相比于模拟数据，有较大的差异性，奇异值波动更大；超出有效频率（10~90 Hz）范围后，测得的奇异值个数出现严重波动；在有效频率范围内，奇异值测得的数量相对较为稳定，但与模拟地震数据不同的是，在有效频率范围内测得的奇异值数量并不完全一样。如何选择恰当的奇异值需进行细致探讨。若选取的奇异值过大，则会影响去噪效果；若选取的奇异值过小，则会损害有效信号。因此，本文对各频率测得的奇异值个数进行统计（表 2），并按测得数量从大到小排列，取其前面90%，选取其中最大值作为该小块最终确认的奇异值，这样选取奇异值可避免损伤大部分有效信号，且不会对去噪效果产生较大影响。

对图 9在有效频率范围内（10~90 Hz）的各频率所测奇异值的个数进行统计（表 2）。可见测得的奇异值个数为1、2、3的频率的累计总数为76，占所有频率个数的93.76%，因此取奇异值个数1、2、3中的最大值3可避免对有效信号造成损伤，且会得到较好的去噪效果，显然此小块在有效频率范围内测得的奇异值的个数为3。

对图 8实际地震数据采用滑动分块重建方法，每小块大小为$ 10\times 10 $，在主测线方向与联络测线方向均取步长为2，共分得256块。图 10a为基于层次聚类法确定奇异值的重建结果，阻尼因子取2，其每块测得的奇异值数量如图 11所示。图 10b为去噪后残留的噪声。可见该实际数据的重建与去噪效果（图 10a）良好，有效信号损失较少，剖面的同相轴轮廓恢复清晰，噪声清除干净，证明基于本文方法在实际地震数据中的有效性与可行性。

对图 8实际地震重建中，由图 11可见大部分数据块中测得的奇异值在1~3。为了验证层次聚类法所测奇异值的准确性，指定所有数据块的奇异值全部为1、2、3，所得到的重建结果与本文每个数据块都自动确定奇异值个数的方法所得重建结果进行对比。

从图 12可见，当对每个数据块的奇异值直接取1时，其图像左下方的断层变得不再明显。对每个数据块的奇异值取2和3时有较好重建效果，但仍然残留少量噪声，且存在线性同相轴恢复不完整的情形。采用层次聚类确定奇异值的方法每个数据块自动确定适合的奇异值个数，重建图像的噪声去除干净，线性同相轴清晰可见，图像重建效果良好。

图 13a为三维实际地震数据实例二，大小为$ 150\times 150\times $$ 276 $，时间上有276个采样点，时间采样率为2 ms，主测线方向与联络测线方向各有150个采样点。对其加入20%随机噪声（图 13b），再对其随机抽稀30%地震道（图 13c）。该数据同样存在噪声能量强、信噪比低、地震道缺失严重等特点。

图 14为针对图 13基于层次聚类法确定奇异值的重建结果，此时阻尼因子取值为2，得到去噪后的残留噪声（图 14b）。可见该重建与去噪结果中，同相轴轮廓恢复清晰，噪声清除干净，有效信号损失较少。再次证明了本文方法对实际地震数据的有效性与可行性。

4 讨论

本文方法需要对地震信号有效频率进行分析，选定有效信号的范围。对地震数据进行分块处理和拼接处理也是必须的，分块的大小对本文方法有一定影响。经验表明：当空间分块大约为$ 10\times 10 $时，能取得稳定的结果；拼接方法的选择对于重建与去噪结果同样会产生影响；另外，阻尼因子也应根据噪声的强弱程度而选取。

5 结论

矩阵降秩方法是地震数据去噪与重建的一种重要方法，秩的选择对于矩阵降秩效果的影响巨大。当秩选取过小时，将损害有效信息；但当秩选取过大时，会在去噪及重建后留下明显的残留噪声。面对海量地震数据，人工选取奇异值的方法具有一定的主观性，且不利于去噪和重建的工业化实现。本文通过对缺道且含有噪声的地震数据的奇异值进行分析，引入了层次聚类方法，将奇异值自动分为有效信号和噪声两类，选择在有效频带范围内实现聚类，使得多道奇异谱分析方法能自适应且稳定地确定奇异值的个数，并通过模拟数据与实际数据验证了方法的优越性。实验表明，利用层次聚类法确定的奇异值数量可靠稳定，重建后线性同相轴清晰可见，无太多残留噪声，地震数据重建效果良好。