微震（microseismicity, MS）是指振动幅值相对较小的地震事件^[1]。工程中所涉及的微震主要指在外界扰动下（如隧道开挖），围岩应力重分布过程中裂纹的生成和扩展将能量以应力波的形式释放和传播的工程地质现象^[2]。微震监测技术通过分析传感器所接受的信号，能够实时地监测岩体发生宏观破坏之前的岩石微破裂信息，为工程安全建设提供指导，并已广泛应用于矿山^[3-4]、水电^[5]等领域。实际工程中，受限于复杂的工程现场情况，各个微震通道的信号都会受到噪音干扰，影响微震初值到时拾取、震源定位、震源类型分析等一系列微震反演工作的开展。因此，对工程现场所接受到的信号进行降噪处理，能够最大程度上获取岩体破裂过程中的应力、变形情况^[6]。

微震信号是一种典型的随机、非平稳信号，分析这类信号的重要之处在于获得信号的时域局部性质^[7]。傅里叶变换（Fourier transform）作为一种经典的降噪方法，虽然能够在一定程度上抑制噪声，但由于其为一种全局变换方法，主要针对周期平稳信号有效，无法准确分析微震这类随机信号的时域局部性质^[8]，且对于含有尖峰的信号降噪效果不佳^[9-10]。小波变换（wavelet transform）在时频域上具有良好的局部性质，能够有效处理随机、非平稳信号，但是降噪过程中需要人为选择合适的小波基、小波分解层数、降噪阈值，如果选择不同，可能会使得结果出现偏差^[11]。此外，小波变换仍然延续傅里叶变换的思想，不能弥补一些傅里叶变换的缺陷^[12-13]（如小波变换和傅里叶变换都受限于海森堡测不准原则）。经验模态分解（empirical mode decomposition，EMD）^[14]作为一种自适应的信号分析方法，在无需人为选择基函数的条件下，能够将信号自适应分解为一系列不同尺度的固有模态函数（intrinsic mode functions，IMF），并通过选择特定的固有模态函数，实现特定的低通、带通、高通滤波^[15]。但EMD 分解过程会产生模态混叠问题，即同一个IMF分量当中出现不同频率或尺度的信号，或同一尺度或频率的信号被自行分解到多个不同的IMF分量中^[16]。虽然集成经验模态分解（ensemble empirical mode decomposition, EEMD）^[17]通过对信号多次加入不同的高斯白噪声，再进行EMD分解，将多个结果求平均获得更加可靠的IMF，但添加高斯白噪声会增加计算量且添加噪声幅值和迭代次数不恰当的时候会产生更多的伪IMF分量，造成较大的结果误差。

基于此，提出一种LMD–SVD联合降噪方法，进行较为高效准确的微震信号降噪。首先，为减缓EMD分解中存在的过度分解、模态混叠的影响，选用局部均值分解法（local mean decomposition，LMD）将带噪信号分解为一组频率由高到低的乘积函数（product function，PF）。通过计算各个PF分量与原信号的相关系数确定有效信号与噪音信号的分界PF分量，将分界处PF分量前的PF分量去除，实现第1步降噪；对于分界PF分量，由于其包含了部分有效信号和噪音信号，采用奇异值分解（singular value decomposition，SVD）对有效信号进行提取，实现第2次滤波；将第2步与第1步结果相叠加，实现LMD–SVD微震信号的联合降噪。

1 基本理论 1.1 局部均值分解

局部均值分解法（LMD）由Smith^[18]提出，是一种通过迭代过程实现的1维信号自适应时频分解方法。该方法由于在迭代次数、端点效应等方面的优势^[19]，可用于分析随机、非平稳信号，如视觉感知的脑电信号^[18]、机械故障信号^[20-21]、地震信号^[22]、微震与爆破信号识别^[23]等。与EMD分解的方法类似，该方法通过迭代的方法从原始信号中提取多个单一频率或者窄频带分布的信号，但从信号分解的角度，二者存在本质的区别。对EMD而言，每次迭代都需要针对上一次迭代中获得的残余信号的极大值、极小值进行3次样条插值（插值方法有多种，例如B样条插值，而3次样条只是在EMD中较为可靠常用的一种，插值方法的讨论可参考文献[13]），进而形成上一步的残余信号的上下包络线，通过求平均值的方式获得均值函数，然后从残余信号中减去，获得该步迭代的残余信号，重复该过程，直至达到停止条件。LMD的方法为更好地刻画随机信号的细节特点，通过构建局部均值和局部包络函数，采用滑动平局的方式进行平滑处理，直至相邻两点值不同^[22]。最后，LMD方法将从原始信号中分解出纯调频函数与纯包络函数，将二者相乘便可获得相应的一个PF分量，然后不断迭代实现所有PF分量的分解。基于文献[18，20，23-25]，将LMD方法的具体流程归纳如下：

1）找出输入信号 $x(t)$ 的所有极值点 ${n_i}$ ，按照式（1）计算所有相邻局部极值点的平均值，并将所有的相邻极值点 ${n_i}$ 和 ${n_{i + 1}}$ 用直线连接，经滑动平均法进行平滑处理以后，获得对应的局部均值函数 ${m_{11}}(t)$ 。

${m_i} = \frac{{{n_i} + {n_{i + 1}}}}{2}$

(1)

2）利用极值点 ${n_i}$ ，按照式（2）获得原始信号的包络估计值，使用直线连接所有相邻的 ${a_i}$ ，采用滑动平局的方式进行平滑处理，获得相应的包络估计函数 ${a_{11}}(t)$ 。

${a_i} = \frac{{\left| {{n_i} - {n_{i + 1}}} \right|}}{2}$

(2)

3）按照式（3），利用原始信号 $x(t)$ 减去局部均值函数 ${m_{11}}(t)$ ，获得 ${h_{11}}(t)$ ；利用式（4），用 ${h_{11}}(t)$ 除以包络估计函数 ${a_{11}}(t)$ ，得到调频信号 ${s_{11}}(t)$ 。

${h_{11}}(t) = x(t) - {m_{11}}(t)$

(3)

${s_{11}}(t) = \frac{{{h_{11}}(t)}}{{{a_{11}}(t)}}$

(4)

4）将 ${s_{11}}(t)$ 视为一个新的输入信号，重复步骤1）～3），获得 ${a_{12}}(t)$ 。在理想情况下， ${s_{11}}(t)$ 作为一个纯调频函数，应该满足 ${a_{12}}(t) = 1$ ；若 ${a_{12}}(t) \ne 1$ ，则多次重复步骤1）～3），获得 ${a_{13}}(t)$ 、 ${a_{14}}(t)$ 、 $\cdots$ 、 ${a_{1n}}(t)$ ，直至 ${a_{1n}}(t) = 1$ ，此时 ${s_{1n}}(t)$ 为一纯调频信号。计算过程如式（5）～（6）所示：

$\left\{ \begin{gathered} \!\!\!\!\!\!{h_{11}}(t) = x(t) - {m_{11}}(t){\text{，}} \\ \!{h_{12}}(t) = {s_{11}}(t) - {m_{12}}(t){\text{，}} \\ \vdots \\ {h_{1n}}(t) = {s_{1(n - 1)}}(t) - {m_{1n}}(t) \\ \end{gathered} \right.$

(5)

$\left\{ \begin{gathered} {s_{11}}(t) = \frac{{{h_{11}}(t)}}{{{a_{11}}(t)}}{\text{，}} \\ {s_{12}}(t) = \frac{{{h_{12}}(t)}}{{{a_{12}}(t)}}{\text{，}}\\ \vdots \\ {s_{1n}}(t) = \frac{{{h_{1n}}(t)}}{{{a_{1n}}(t)}} \\ \end{gathered} \right.$

(6)

需要说明的是，理论上应该满足：

$\mathop {\lim }\limits_{n \to \infty } {a_{1n}}(t) = 1$

(7)

但在实际数值计算过程，式（7）是不能完全满足。将式（8）作为LMD迭代的停止条件。

$\left| {{a_{1n}}(t) - 1} \right| < 0.05$

(8)

5）将步骤4）迭代过程中产生的所有包络估计函数按照式（9）相乘获得包络信号（部分文献也将其称为瞬时幅值函数，instantaneous amplitude，IA）。

${a_1}(t) = {a_{11}}(t){a_{12}}(t) \cdot \cdot \cdot {a_{1n}}(t) = \prod\limits_{q = 1}^n {{a_{1q}}(t)} $

(9)

6）将包络信号 ${a_1}(t)$ 与纯调频信号 ${s_{1n}}(t)$ 按照式（10）相乘得到原始信号的第1个PF分量。

$P{F_1}(t) = {a_1}(t){s_{1n}}(t)$

(10)

$P{F_1}$ 作为一个调频–调幅信号，包含了原始信号中最高的频率成分，其瞬时幅值即为包络信号 ${a_1}(t)$ ，其瞬时频率 ${f_1}(t)$ 可通过纯调频信号 ${s_{1n}}(t)$ 按照式（11）计算获得：

${f_1}(t) = \frac{1}{{2{\text{π}} }}\frac{{{\rm{d}}[\arccos({s_{1n}}(t))]}}{{{\rm{d}}t}}$

(11)

7）将 $P{F_1}$ 分量从原始信号 $x(t)$ 中分离，获得信号 ${u_1}(t)$ 。将 ${u_1}(t)$ 作为新的原始信号，重复步骤1）～6），经过 $k$ 次循环，直至 ${u_k}(t)$ 为单调函数停止，如式（12）所示：

$\left\{ \!\! \begin{gathered} {u_1}(t) = x(t) - P{F_1}(t) {\text{，}}\\ {u_2}(t) = {u_1}(t) - P{F_2}(t){\text{，}} \\ \vdots \\ {u_k}(t) = {u_{k - 1}}(t) - P{F_k}(t) \\ \end{gathered} \right.$

(12)

最终原始信号 $x(t)$ 被分解为所有PF分量与 ${u_k}(t)$ 之和（ ${u_k}(t)$ 为LMD分解的余项），即：

$x(t) = \sum\limits_{p = 1}^k {P{F_p}(t) + u{}_k(t)} $

(13)

1.2 奇异值分解降噪 1.2.1 SVD降噪理论

对一个含噪信号 $X(N) = \{ x{}_1,x{}_2, \cdots ,x{}_N\} $ ，按式（14）进行相空间重构，形成 $m \times n$ 阶的Hankel矩阵。

${{{H}}_{m \times n}} = \left( {\begin{array}{*{20}{c}} {{x_1}}&{{x_2}}& \cdots &{{x_n}} \\ {{x_2}}&{{x_3}}& \cdots &{{x_{n + 1}}} \\ \vdots & \vdots & \vdots & \vdots \\ {{x_m}}&{{x_{m + 1}}}& \cdots &{{x_N}} \end{array}} \right) = {{{D}}_{m \times n}} + {{{W}}_{m \times n}}$

(14)

式中， $N = m + n - 1$ ， ${{{D}}_{m \times n}}$ 为含噪信号中有效信号子空间， ${{{W}}_{m \times n}}$ 为噪声信号子空间。

基于式（14），按照式（15）进行奇异值分解：

${{{H}}_{m \times n}} = {{{U}}_{m \times m}}{\varSigma}_{m \times n}{{V}}_{n \times n}^{\rm{T}}$

(15)

式中： ${{{U}}_{m \times m}}$ 和 ${{{V}}_{n \times n}}$ 为正交矩阵； ${{\varSigma}_{m \times n}}$ 为一个非负对角型矩阵，其秩为 ${{r}}$ ， ${{\varSigma}_{m \times n}}$ 可表示为 ${\varSigma} = \left( {\begin{array}{*{20}{c}} {{S}}&0 \\ 0&0 \end{array}} \right)$ ， ${{S}} = diag$ $ ({\sigma _1},{\sigma _2}, \cdots ,{\sigma _r})$ ， ${\sigma _i}$ 为矩阵 ${{{H}}_{m \times n}}$ 的奇异值，表示为 ${\sigma _i} = $ ${\varSigma} (i,i)$ ，通常满足 ${\sigma _1} \ge {\sigma _2} \ge \cdots \ge {\sigma _r} \ge 0$ 。根据奇异值分解定理和Frobeious范数下的矩阵最佳逼近定理^[26-27]可知：信号中的有效信号主要由前K个奇异值体现，而高频的噪声信号由后面部分的奇异值体现，因此，通过保留前面K个奇异值且令其他奇异值为0，形成新的对角矩阵 ${{\varSigma}_{m \times n}}$ ，再带入式（15）进行奇异值分解的逆变换获得降噪后的 ${\hat {{H}}_{m \times n}}$ ，并获得降噪后的信号 $\hat X(N)$ 。

1.2.2 Hankel矩阵维数的确定

对于相同的一个输入信号，可重构多种结构的Hankel矩阵，且不同结构会使得信号的SVD分解结果产生一定误差，因而会直接影响到信号的降噪结果^[28]。文献[29]对多组不同长度和不同频率的信号分析，确定构造矩阵 ${{{H}}_{m \times n}}$ 的最佳维数在 $m = N/2$ 左右，获得满足工程应用的降噪结果。故分别按照式（16）～（17）确定Hankel矩阵行数 $m$ 和列数 $n$ 。

$m = \left\{ \begin{array}{l} \!\!\!\!N/2{\text{，}}N{\text{为偶数}}{\text{；}}\\ \!\!\!\!(N + 1)/2{\text{，}}N{\text{为奇数}} \end{array} \right.$

(16)

$n = N + 1 - m$

(17)

1.2.3 奇异值分解阶数的确定

如第1.2.1节中所述，SVD降噪需要确定所保留奇异值个数K的大小，即奇异值有效阶数的确定。当K较小时，部分有效信号会被视为噪声去除；当K较大时，会有部分噪声未被去除，使得信号上出现“毛刺”，有较明显的波动。当前，在奇异值降噪过程中，应用最为广泛的奇异值有效阶数的确定方法是奇异值差分谱的方法，但该方法的缺点在于当相邻的两个奇异值相对于其他奇异值都大得多而二者之间也有较大差值时，奇异值差分谱也会在这两个奇异值的部分产生峰值，故会导致有效阶数选择不准确的问题。基于对文献[30]的分析，按照式（18），采用加权能量贡献率（percent of contribution to total energy，PCTE）作为奇异值有效阶数K的确定方法。值得注意的是，文献[30]给出的方法需要人为观察判断PCTE值是否大致为0。为提升判断效率，在经过对多组信号分析的基础上，采用保留PCTE值大于0.1%的奇异值以确定奇异值有效阶数。

$PCT{E_{{\sigma _i}}}\!=\! \left( {1\! -\! \frac{{\sqrt {\sigma _1^2 \!+\! \sigma _2^2\! + \! \cdots \! +\! \sigma _{i - 1}^2 \!+\! \sigma _{i + 1}^2 \!+\! \cdots \!+ \!\sigma _n^2} }}{{{{\left\| {{H}} \right\|}_F}}}} \right) \times 100{\text{%}}$

(18)

式中， ${\left\| {{H}} \right\|_F}$ 为矩阵 ${{{H}}_{m \times n}}$ 的Frobenious范数。

2 LMD–SVD联合降噪

LMD分解同EMD分解类似，会将原始信号分解为一系列从高频到低频的信号。早期的降噪思路是结合现场情况，以及有效信号与噪声信号特点，人为去掉某些特定的分量，将剩下的部分进行组合，获得降噪后的信号，如文献[15]所示。然而，目前对于LMD、EMD及其改进的相关算法，无论采用何种改进方式，都无法根本解决信号分解过程中的模态混叠的影响。因此，简单地舍弃某些特定分量会使得部分有效信号损失、部分噪音信号保留。此外，如果每次分解过程中都需要人为判断合适的分量，不但会受不同个体评价标准的影响，且会影响计算分析效率。针对这一问题，引入相关算法实现自动确定高频噪声与低频有效信号分界分量的位置，通过直接去除分界PF分量之前的高频部分实现第1步降噪；之后，引入SVD法从受模态混叠影响的分界PF分量信号中提取有效微震信号，实现第2步降噪。

2.1 分界PF分量的确定

通常而言，多数的工程噪声相对于有效信号为高频信号。针对LMD与EMD分解结果表现为从高频到低频分布的特点，相关学者引入信息论中的相关系数（correlation coefficient, CC）^[6，31]和K–L散度法（Kullback–Leibler divergence，K–L divergence）^[32]确定噪声信号、有效信号和二者的分界分量。一般而言，K–L散度法优点在于区分度较大，差别效果较为明显，但缺点在于适应性差，较为敏感；相关系数法适应性较好，但有时会出现区别效果不明显。结合文献[6,31]在微震信号分解及分量选取中的应用，再考虑到适用性问题，本文选择相关系数法作为分界分量的确定方式，即首先采用式（19）计算出 $ P{F_i}(i =1, $ $2, \cdots ,n)$ 与原信号 $x(t)$ 的相关系数 ${\,\rho _{P{F_i},x(t)}}$ ，然后将各个相关系数按照计算先后顺序排列，将相关系数序列中出现的第1个极值点位置确定为分界分量位置。

${\rho _{P{F_i},x(t)}} = \frac{{Cov(P{F_i},x(t))}}{{\sqrt {D(P{F_i})} \cdot \sqrt {D(x(t))} }}$

(19)

式中， $Cov(P{F_i},x(t))$ 为 $P{F_i}$ 与 $x(t)$ 的协方差， $D(P{F_i})$ 为 $P{F_i}$ 的方差， $D(x(t))$ 为 $x(t)$ 的方差。

2.2 联合降噪处理过程

在获得分界PF分量后，针对第2节指出直接去掉LMD分解的高频分量所存在的问题，引入SVD法提取分界PF分量中的有效信号（具体过程参考第1.2节）。将SVD处理结果与LMD分解结果保留的低频分量相加，可得到LMD–SVD联合降噪的结果。详细流程如下：

1）输入含噪微震信号，进行LMD分解，获得各个PF分量。

2）利用式（19）计算各个PF分量与输入含噪信号的相关系数。

3）寻找相关系数序列的第1个极值点，并将其作为噪声信号与有效信号的分界分量。

4）剔除分界分量之前的所有高频信号。

5）对分界信号按照第2节中的方式进行SVD降噪。

6）将第5）步获得的信号与剩下的PF分量进行叠加，获得LMD–SVD降噪结果。

3 仿真分析

为验证本文方法的合理性，引入模拟地震的Ricker子波进行仿真实验。Ricker子波的原始表达式为：

$f(t) = \left(1 - 2{{\text{π}} ^2}f_{\rm{p}}^2{t^2}\right)\exp\left( - {{\text{π}} ^2}f_{\rm{p}}^2{t^2}\right)$

(20)

式中： $f(t)$ 为振幅； ${f_{\rm{p}}}$ 为谱峰频率，本文取 ${f_{\rm{p}}} = 35$ Hz，采样频率为1 kHz，采样点数为 1 000。为了展示方便，将式（20）向右沿 $t$ 轴平移0.3 s，如式（21）所示，其余参数不变，其波形图和频谱图分别见图1（a）、（b）。向Ricker子波中添加高斯白噪声，获得加噪后信号的波形图和频谱图如图1（c）、（d）所示。

$f(t) = \left[1 - 2{{\text{π}} ^2}f_{\rm{p}}^2{(t - 0.3)^2}\right]\exp\left[ - {{\text{π}} ^2}f_{\rm{p}}^2{(t - 0.3)^2}\right]$

(21)

在获得加噪的Ricker子波之后，对该信号分别进行EMD和LMD分解，对分解后产生的各个分量分别求其频谱，分量的波形图和频谱图分别如图2和3所示。

为说明降噪效果，一方面，结合波形图和频谱图进行定性分析降噪效果；另一方面，引入信号学中的信噪比公式（22）作为定量分析指标：

$SNR = 10{\ln }\left( {\frac{{\displaystyle\sum\limits_{i = 1}^N {x{{(t)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^N {{{(x(t) - x(t)')}^2}} }}} \right)$

(22)

式中， $x(t)$ 为未受噪声污染的信号， $x(t)'$ 为采用某种手段降噪后的信号， $N$ 为信号采样点数目。仿真实验中，在已知未受噪声污染信号的条件下，信噪比越大，表明噪声的含量越小。

对带噪的Ricker子波进行EMD和LMD分解，结果分别如图2和3所示。图2的EMD分解中出现了明显的模态混叠现象，图3中的LMD分解结果相对减缓了模态混叠。

此外，由图2和3可知，IMF1与PF1分别代表EMD和LMD两种方法分解结果中的高频噪音，将其剔除，然后将剩下的分量合并形成降噪后的信号，降噪后的波形图和频谱图见图4。

由图2～4可知，对EMD和LMD算法而言，去掉第1个分量很大程度上减少了信号中的高频噪音，但是由于模态混叠的影响，高频部分会混合在其他分量当中，简单去掉这些分量会造成有效信号的损失。

针对上述问题，对LMD分解结果进行本文LMD–SVD方法降噪，将计算所得的各个分量与带噪信号之间的相关系数结果记录于表1中。结合表1可知，显然降噪过程是去掉PF1分量，然后对PF2分量进行SVD降噪处理，降噪结果的波形图和频谱图如图5所示。

利用式（22）计算各个方法降噪前后的信噪比，将其记录于表2中。由表2可知：由于LMD对模态混叠现象的减缓，使得其降噪效果明显优于EMD降噪后的结果；LMD–SVD由于引入了SVD方法处理PF2分量，使得降噪后的波形相对更平滑，拥有略高的信噪比。

4 应用实例

为验证本文方法在实际工程中运用的可行性，对布置于白鹤滩水电站左岸地下厂房的ESG微震监测系统所采集的含噪微震信号进行降噪处理。

白鹤滩水电站位于四川省宁南县和云南巧家县交界位置，是金沙江下游干流河段梯级开发的第2个水电站。左岸地下厂房为超大型地下厂房，围岩主要由隐晶质玄武岩、斜斑玄武岩、杏仁状玄武岩、角砾熔岩等组成，发育有原生构造和断裂构造^[33]。为研究地下洞室开挖影响下的围岩稳定性，于2014年10月10日将ESG微震监测系统引入白鹤滩左岸地下厂房。该监测系统由Hyperion数字信号处理系统、Palading数字信号采集系统、6通道单轴加速度的传感器及3维可视化软件等组成（传感器空间布置见图6）。在监测系统运行过程中，传感器采样频率为10 000 Hz，每个信号采样点数为15 000（具体参数见表3）。使用Hyperion数字信号处理系统中的WaveVis波形分析软件将待分析信号以.txt格式导出，在MATLAB平台上导入，并编写相关程序，进行含噪微震信号的降噪处理。

为简洁展示降噪效果，从中选出4组含噪微震信号进行降噪处理演示，结果如图7所示。由图7可知，微震信号中夹杂有部分高频随机噪声，尤其是信号4，受到较大噪声干扰，这对后续微震信号的初值时刻拾取、震源机制解分析、中心频率和主频范围分析等反演工作造成了一定干扰。经过LMD–SVD降噪处理后的信号，由频谱图中可知，高频的随机噪被很好地消除，几乎没有残留，且在微震信号频带范围与噪声频带范围的分界位置，没有出现频域截断现象，这使得降噪后的结果在时域波形上更加接近微震信号。由此可见，本文方法对微震信号的降噪处理具有一定的可行性。

5 结　论

1)结合相关文献及本文仿真实验可知，EMD分解虽然能够将原始信号分解为一系列由高频到低频的信号，但无法避免受到模态混叠现象的影响。LMD分解不同于EEMD等EMD的改进算法，能够在不添加高斯噪声的条件下，对模态混叠起到一定的抑制作用。

2)充分发挥LMD分解的上述优点，结合SVD方法从有效信号与噪声信号分界分量中提取有效信号分量，提出一种基于LMD–SVD的微震信号降噪方法。通过仿真实验表明该方法相比于EMD分解与LMD分解简单地剔除高频分量，本文方法能够更好地去除含噪信号中的高频信号。

3)通过分析白鹤滩左岸地下厂房的含噪微震信号，论证了本文方法在处理含噪微震信号时既不会在降噪过程中造成频域截断而使得时域波形失真，又能很好地去除信号中高频分量，进而从含噪信号中获得有效的微震信号。