0 引言

我国于2008年8月1日开通首条商业运营高铁线路——京津城际铁路.根据2016年7月13日联合发布的《中长期铁路网规划》，到2025年我国的高速铁路将达到3.8×10⁴ km左右.截止2017年底，我国高铁营业里程达2.5×10⁴ km，占世界高铁总量的66%.每天有数千趟高铁列车高速运行在分布范围很广的高铁线路上.高速运行庞大数量的高速运行列车，不但会引起高铁列车的振动，还会引起铁路路基的振动.很多学者对高铁路基的振动进行过数值模拟及地质力学方面的分析，并对轨道及高铁线路几十米范围内的小规模做振动数据采集(Forrest，1999；Degrande and Schillemand, 2001; Forrest and Hunt, 2006；Kouroussis et al., 2014；郑亚玮和陈俊岭，2015)，并对相应数据进行了时域及频域分析.例如：Degrande和Schillemand(2001)采集了布鲁塞尔—巴黎间某段高铁轨道及附近70m范围内的振动信号，郑亚玮和陈俊岭(2015)采集了京津高铁某段轨道35m范围内的振动信号，均发现高铁轨道附近的振动呈现出时域周期性及频域窄带分立谱等特征，印证了高铁震源的周期性、窄带分立谱等特征.高铁列车高速运行在高铁线路上，除了会引起车厢及路基的振动，而且会将振动以各种类型的地震波传播出去.徐善辉等(2017)于2013年在京津城际高铁廊坊段使用1000多个地震检波器在距离高铁1 km外长时间采集地震数据，发现高铁列车激发的振动信号可传播几公里远.

高铁列车高速运行引起的振动可能会威胁建筑物、工程结构的安全，工程界视之为一大危害.然而，高铁列车的高速运行也会产生了一种全新的可重复震源类型——高铁震源.因此对这种全新的可重复震源-高铁震源所激的地震波进行采集并分析，不但可能分析高铁列车的运行状态，而且有望对高铁线路附近的地下结构进行成像.然而，实际采集到的数据不但包含有特殊的高铁震源激发出的地震波，还有各种各样的背景噪声.如何从检波器接收到的数据中分离出具有窄带分立谱特征的高铁震源地震信号和宽频带背景信号，对充分利用高铁震源信号是极为关键的.工业界常用陷波器来压制单频干扰信号，但考虑到高铁震源地震信号的能量集中在很多个窄频带上，采用陷波器不但需要知道频率成分的位置并以此设计很多个陷波器进行窄带信号滤除，而且在频带较多时会损伤宽频带信号(Xu et al., 2013)，因此需要发展新的方法以实现高铁震源地震信号和宽频带背景信号的分离.

随着信号稀疏表示理论发展(Chen et al., 2001; Mallat and Zhang, 1993)，Starck等(Starck et al., 2004; Starck et al., 2005; Bobin et al., 2017)提出形态成分分析(Morphological Component Analysis, MCA)理论，将具有不同形态特征的字典构成一组超完备冗余字典，可获得对复杂信号的更稀疏表示.根据MCA理论，实现复杂信号中各成分的稀疏表示和分离的两个关键为：(1)能否为复杂信号中不同的成分选择不同的稀疏表示字典；(2)每个稀疏表示字典仅能稀疏表示对应的成分，无法对其他成分进行稀疏表示.前期很多关于高铁振动的频谱分析已表明高铁震源地震信号的能量多集中在多个窄带上(Degrande and Schillemans, 2001；郑亚玮和陈俊岭，2015；王晓凯等，2019)，同时很多前期关于普通列车和高铁振动的研究也表明震源函数在频域具有明显的窄带分立谱特征(和振兴和翟婉明，2007；Sheng et al., 2004；Chen et al., 2011；Liu and Luo, 2019¹⁾).但是检波器接收到的信号不只包含窄带分立谱特征的高铁震源地震信号，还包含各种各样的宽频带背景信号.如何从复杂信号中分离出高铁震源地震信号是准确利用高铁震源地震信号的关键.考虑到两者在频域呈现出完全不同的形态(一个为窄带分立谱特征，而另一个为宽频带特征)，可使用局部余弦变换来稀疏表示高铁震源地震信号，而使用连续小波变换来稀疏表示宽频带背景信号.因此，根据MCA理论并结合分块坐标松弛法可实现两者的分离，为后续的分析和处理提供各自所需的数据.

1) Liu Y J, Luo Y. 2019. Explore High-Speed-Train-Induced vibration as a new seismic source. Geophysical Research Letters (Under review).

本文首次将形态成分分析引入到高铁震源地震信号处理中，实现高铁震源地震信号与宽频带背景信号的分离：首先介绍信号的稀疏表示理论及形态成分分析；然后分析实际高铁震源地震数据及理论震源函数的振幅谱特征进而构建相应的稀疏表示字典，并进一步构建超完备字典；最后将MCA理论并结合分块坐标松弛法应用于中国南方某高铁线路附近实际采集数据的分离，并用振幅谱和时频谱说明了方法的有效性.

1 形态成分分析 1.1 信号的稀疏表示理论

在信号处理领域，通过各种数学变换将信号转换到变换域进行分析，在实际应用中具有简单性和有效性.传统的基于非冗余正交基函数的数学变换，如小波变换(Daubechies, 1990)、短时傅里叶变换(Griffin and Lim, 1984)等，通常不能有效地表示复杂信号.近年来许多新的数学变换，如Ridgelet变换(Donoho, 2001)、Curvelet变换(Candès et al., 2006)、Shearlet变换(Easley et al., 2008)等，区别于传统的正交变换方法，能够更丰富地表示复杂信号.Mallat和Zhang(1993)提出了信号稀疏表示理论中最重要的概念，即字典.字典中的元素被称为原子，信号通过原子的线性组合来表示，其中字典内原子的数目远远高于信号的维数，并且Mallat通过对比说明了稀疏表示的方法与传统的时频分析方法相比，具有更好的优越性.诸如小波变换、短时傅里叶变换和Ridgelet变换等数学变换，都可以作为稀疏表示某种信号的字典.信号能够稀疏表示，是指信号能够用尽量少的变换系数来表示，即可以用尽可能少的变换原子的线性组合来表示信号.而使用多种字典构成的超完备字典，可以实现对复杂信号更稀疏的表示，其中选择的字典可以分别稀疏表示复杂信号的各个成分.自从超完备冗余字典的理论被广泛应用以来，构建冗余基函数稀疏表示信号越来越受到学者们的关注.

图 1是信号稀疏表示的示意图.假设给定任意信号s∈R^N与对应的超完备字典Φ，X为信号s使用超完备字典Φ的稀疏的表示系数，我们要求该超完备字典的维数L≥N.图中字典Φ中参与信号重构的原子以黑色标记，并在系数矩阵X中的对应位置也用黑色标记，表示此处系数不为零，没有参与重构的原子以灰色表示，并在系数矩阵X中的对应位置用白色标记，表示此处系数为零.从图 1可以看出，系数矩阵X中非零项的个数(即黑色块)占系数矩阵的总数的比例非常少，即字典Φ获得了信号s的比较稀疏的表示形式.

使用超完备字典Φ对任意的信号s∈R^N求取最稀疏表示，需要求解如下所示的优化问题：

(1)

式中，‖X‖₀是系数矩阵X的稀疏性度量，即系数矩阵X中的非零元素个数.

1.2 形态成分分析

基于信号稀疏表示理论的形态成分分析是解决图像或地震信号的多成分分离问题的一种方法.假设某种待分析复杂信号s∈R^N，由两种不同形态特征的信号分量和随机噪声组合而成：

(2)

其中，s为待分析信号，s₁、s₂为信号中的两种成分，n为随机噪声，而在实际应用时我们经常会忽略随机噪声，式(2)可写成如下的形式：

(3)

形态成分分析方法的假设条件，是对于混合信号中的每种信号分量，存在一种只对该信号分量稀疏表示而对另一信号分量不稀疏的字典.在式(3)模型中，s₁和s₂可以分别由字典Φ₁和Φ₂有效地稀疏表示，但是用Φ₂稀疏表示s₁以及用Φ₁稀疏表示s₂时稀疏性差.在形态分量分析理论中，确定两个分别对两种信号分量进行稀疏表示的字典Φ₁和Φ₂十分重要.在信号稀疏表示理论中，字典的选取构造有两类：一类是根据经验选择常用的数学变换作为字典，例如离散Ridgelet变换、Curvelet变换和Shearlet变换等；另一类是根据信号自适应计算能够稀疏表示信号的字典，即自适应字典，主要的算法有MOD算法(Engan et al., 1999)和K-SVD算法(Aharon et al., 2006)等.

基于以上的假设条件，对于含有两种信号s₁、s₂的混合信号s，分别选取只对其稀疏表示的字典Φ₁和Φ₂，并联合Φ₁和Φ₂构造一对超完备字典来稀疏表示复杂信号s，可得到信号s的最稀疏表示形式，即求解优化公式(4)：

(4)

式中，X₁为重构系数中与Φ₁对应的部分；X₂为重构系数中与Φ₂对应的部分.

但是，式(4)所示的优化问题的求解有两大困难，首先是优化问题为l₀范数下的非凸问题，不容易直接求解，其次是实际信号受到噪声的干扰，不容易找到完全匹配信号分量的字典.因此，对于式(4)所描述的问题，我们将系数矩阵稀疏性度量的l₀范数转化为l₁范数，并适当松弛等式约束条件，转化为无约束条件，使不可解问题变得如下可以求解的最优化问题：

(5)

针对式(5)所示的优化问题，Bruce等(1998)提出分块坐标松弛(Block Coordinate Relaxation, BCR)算法进行求解.BCR算法的核心思想是设定一个合理的阈值策略，每次迭代对稀疏系数X₁和X₂进行交替更新，直到达到迭代终止条件.BCR算法的具体内容如下：

初始化：初始迭代步数k=0，初始系数解：

(6)

迭代：每步迭代步数k增加1，并且如式(7)所示的方式来交替计算稀疏系数X₁和X₂：

(7)

其中，Φ₁^*、Φ₂^*分别是Φ₁、Φ₂的伪逆.需要指出的是，稀疏系数X₁和X₂除了上述这种顺序更新的方法，还有另一种实现方式，即每步迭代时并行更新，如式(8)所示：

(8)

其中，T_λ1^k、T_λ2^k分别为变换Φ₁^*、Φ₂^*第k步的阈值参数，可以为硬阈值或者软阈值.

终止条件：当‖X^k－X^k－1‖₂²小于预设值，即继续迭代对最终结果影响足够小时，迭代终止.这里所指的预设的值也称为BCR算法迭代的最小阈值参数.

输出：得到最优解X₁^opt, X₂^opt：

(9)

得到最优解X₁^opt, X₂^opt之后，用字典Φ₁和稀疏表示系数中与Φ₁对应的部分重构s₁，用字典Φ₂和稀疏表示系数中与Φ₂对应的部分重构s₂，即

(10)

从而实现从混合信号s中分离出两种信号分量s₁和s₂的目标.

2 高铁震源地震信号的稀疏化建模

2018年1月，北京大学地球与空间科学学院在中国南方某段高铁线路附近布置了51个三分量低频检波器，检波器频带宽度为0.2~100 Hz，采样间隔为5 ms.我们选取其中的17个检波器进行分析处理，17个检波器(编号为T1至T17)与高铁线路的位置示意图见图 2，检波器阵列与高铁线路基本垂直.最近的检波器距离高铁大约75 m，最远检波器距离高铁线路大约187 m，检波器空间间隔大约为7 m.列车1经过时，检波器T1接收到三分量波形如图 3所示，将三个分量的总体振幅谱示于图 4a.由于采集时间范围内的通过列车较多，将多趟高铁列车经过时检波器T1所接收到的三分量信号振幅谱进行累加，结果如图 4b所示，大部分能量集中在少数几个窄带分立谱上且分立谱的位置与单趟列车经过基本相同，这也说明高铁列车震源的可重复性极好且呈现出明显的窄带分立谱特征.

前期诸多学者在研究高铁铁轨及路基振动时，将移动的列车对钢轨某处的作用力建模为多个延迟点力的累加(和振兴和翟婉明，2007；Sheng et al., 2004；Chen et al., 2011；Liu and Luo, 2019¹⁾)可得到理论上的震源函数.若一列列车中有M对车轮，高铁列车匀速通过铁轨且速度为v，假设列车车头通过铁轨上A点的时刻为0时刻，第m对车轮到车头的距离为d(m)，则A点受力为(Liu and Luo, 2019¹⁾)：

(11)

其中Q为一常数.A点受力f(t)的傅里叶变换为F(ω)：

(12)

若假定一节车厢长度为25 m，一节车厢四对轮子距该节车厢头部的距离分别为4 m、6.5 m、18.5 m及21 m，则对应A点受力的振幅谱如图 4c所示.该振幅谱也呈现出明显分立谱特征.

无论是从单趟列车经过时检波器所接收到信号的平均振幅谱(图 4a)上观察，还是从多趟列车经过时检波器所接收到信号的平均振幅谱(图 4b)上观察，大部分能量都集中在少数几个很窄的频带上.但除此之外，振幅谱上还存在一些宽频带信号.因此将窄带高铁震源地震信号与宽频带信号进行分离有利于充分利用高铁震源地震信号.

基于窄带高铁震源地震信号与宽频带信号的不同频域特性，本文选取离散余弦变换来稀疏表示窄带高铁震源地震信号，选取连续小波变换来稀疏表示宽频带背景信号.其中离散余弦变换的定义为

(13)

式中，x[n]表示待分析信号的时域采样点，CT_x(k)表示离散余弦变换系数，采样长度为N，ω(k)如下

(14)

而离散余弦变换的反变换为

(15)

此外，离散余弦变换是针对实信号所定义的一种变换，实信号正变换后得到的还是一个实信号，而且一维地震数据做离散余弦变换得到的系数也是一维的，因此运算量很小，在数据处理速度上有很大的优势.

连续小波变换定义为

(16)

式中，x(t)为待分析信号，WT_x(a, τ)为变换系数，a为尺度因子，ψ(t)为Morlet母小波.而连续小波变换反变换为

(17)

式中，常数为容许条件.

3 实际数据处理

将上节选择的局部余弦变换和连续小波变换联合构成超完备字典，然后利用BCR算法对中国南方某地的高铁震源数据进行分离以提取窄带高铁震源地震信号和宽频带信号.图 5a为列车1经过时T1检波器所接收的垂直地表分量，图 5b和图 5c分别为分离出的窄带高铁震源地震信号和剩余宽频带信号.从列车1经过时T1检波器所接收的垂直地表分量(图 5a)上可以观察到，从大约4.5 s开始有16个较为稳定周期波形，但比较模糊.经过估计，发现此波形的周期大约为0.3 s.考虑到我国目前高铁运营速度大约为300 km·h^-1(83.33 m·s^-1)，波形中的周期现象对应的距离为大约25 m，而目前我国高铁车厢长度大约为25 m，编组为16节车厢.因此，波形中的周期现象及个数与车厢长度及个数吻合，应为车厢引起的周期现象.但图 5a中从4.5 s开始的周期现象比较模糊.将列车1经过时T1检波器所接收的垂直地表分量进行分离，得到的窄带高铁震源地震信号如图 5b所示，可以明显地观察到从4.5 s以后的周期现象变得更加清晰，而且16个周期波形的一致性也变好，说明提取出的窄带信号能够更加清晰地反映列车的运行状态等.

进一步地，在频域分析高铁震源地震信号的分离结果.图 6a为列车1经过T1检波器时所接收的垂直地表分量振幅谱，可以看出大部分能量集中在少数几个窄带内，但有明显的宽频带背景信号.图 6b为提取出的列车1经过T1检波器时垂直地表分量的窄带高铁震源地震信号的平均的振幅谱，窄带特征极为明显且宽频带背景信号的能量大幅削弱.图 6c为列车1经过T1检波器时垂直地表分量的剩余宽频带信号的平均振幅谱，呈现出很宽的频带.频域的分析结果表明，采用本文方法可以较好地分离出窄带高铁震源地震信号和宽频带背景信号.

图 7为高铁震源地震信号分离结果的时频域分析.图 7a为列车1经过T1检波器时所接收的垂直高铁分量的时频分布(采用加窗傅里叶变换)，频率成分基本不随时间变化，但也能观察到列车引起的宽频带信号.图 7b为提取出的列车1经过时T1检波器时垂直高铁分量的窄带高铁震源地震信号时频分布，频率成分随时间基本不变化这个特征十分清晰，且不受宽频带背景信号的影响.图 7c为列车1经过T1检波器时垂直高铁分量的宽频带背景信号时频分布，呈现出很宽的频带.因此，时频域的分析结果也表明，采用本文方法可以较好地分离出具有窄带分立谱特征的高铁震源地震信号和宽频带背景信号.

4 结论

前期的实际数据分析表明高铁旁的检波器不但能够接收到高铁震源地震信号，同时也接收到宽频带信号，而频谱分析结果及理论震源函数振幅谱均表明高铁震源地震信号的能量在频率域集中在多个窄带内，具有窄带分立谱特征.本文首次将形态成分分析引入到高铁震源地震信号处理中，并根据高铁震源地震信号与宽频带背景信号在频率域的形态成分差异分别选择局部余弦变换和连续小波变换作为两个稀疏表示字典，并将两个变换相结合形成超完备字典，然后利用分块坐标松弛法实现高铁震源地震信号与宽频带背景信号的分离.利用本文方法对北大在中国南方某高铁线路附近实际采集到的高铁震源地震数据进行处理，结果表明：无论是在时间域、频率域还是时频域，利用本文方法可实现带高铁震源地震信号与宽频带背景信号的分离.然而，本文只是利用形态成分分析理论及分块坐标松弛法实现了高铁震源地震信号和宽频带背景信号的分离，但对如何利用高铁震源地震信号和宽频带背景信号并未做探讨，关于后续数据的应用还有待进一步研究.

致谢  感谢中国科学院地质与地球物理研究所李幼铭研究员的有益讨论及他的大力支持和鼓舞.