0 引言

我国自从1998年起开始即通过研究摆式列车在中国铁路既有线路实现提速至高速铁路的可行性，随后投入大量科研力量及资金到高铁的研究中，并于2008年8月1日开通运营我国首条商业运营高铁线路——京津城际铁路.截至2017年底，我国高铁营业里程达2.5万公里，占世界高铁总量的66%.每天有数千趟高铁列车高速运行在分布范围很广的高铁线路上，例如：仅广深港高速铁路在运营初期每日即开行127对列车.

如此庞大数量的高铁列车高速运行在高铁线路上，势必会引起高铁列车的振动.通过车载传感器大规模收集高铁列车振动数据，并通过数据分析可用于高铁故障诊断等(赵晶晶等，2014；朱菲和金炜东，2018).例如可通过在高铁上的传感器采集到列车运行过程中的振动信号，然后利用信息熵和经验模态分解(EMD)分析采集到的振动信号并提取振动数据特征，进而采用各种分类方法可实现对列车运行过程中出现的机械故障进行诊断(李贵兵等，2014).方松和曾京(2013)以郑州—济南线路为测试线路，截取车辆从启动、加速到200 km·h^-1并持续运行60 s、减速这个过程，选取动车组的轴箱和车体垂向振动数据作为研究对象，然后分别利用快速傅里叶变换、短时傅里叶变换、连续小波变换、希尔伯特黄变换进行分析并比较.

高铁列车高速运行在高铁线路上，除了自身振动以外，还会引起高铁路基的震动.一些学者进行过数值模拟、地质力学方面的分析，进行了轨道及高铁几十米范围内的小规模震动数据采集，并对相应数据进行了时域及频域分析(Forrest，1999；Degrande and Schillemans, 2001; Forrest and Hunt, 2006；Kouroussis et al., 2014；郑亚玮和陈俊岭，2015).例如：Degrande和Schillemans(2001)采集了布鲁塞尔—巴黎段高铁轨道及附近70 m范围内的震动信号，郑亚玮和陈俊岭(2015)采集了京津高铁轨道35 m范围内的震动信号，均发现了高铁轨道附近的震动呈现出时域周期性及频域窄带分立谱等特征，印证了高铁震源的周期性及窄带分立谱等特征.

高铁列车高速运行在高铁线路上，不但会引起车厢及路基的震动，而且会将震动以各种类型的地震波传播出去. Chen等(2004)利用9个检波器(包含有长周期检波器和短周期检波器)首次对大秦铁路上运行的重载低速列车进行观测，认为高铁激发的震动有潜力成为检测地壳结构的手段，但受限于分析手段，仅从时间域和频率域出发做了有益尝试.徐善辉等人于2013年在京津城际高铁廊坊段使用1000多个地震检波器在距离高铁1 km外长时间采集地震数据，发现高铁列车激发的振动信号可传播几公里远(徐善辉等，2017).

高铁列车高速运行引起的振动可能会威胁建筑物、工程结构的安全，工程界视之为一大危害.因此目前大部分研究聚焦于如何减少车厢及路基的振动以提高高铁列车及路基的安全系数.然而，高速运行的高铁列车也产生了一种全新的震源类型——高铁震源.高铁列车确定的长度和荷载以及均匀的运动速度使高铁震源具有地震学家所期望的可重复性；高铁震源分立谱特征为高铁附近大范围结构和物性的高精度探测提供了特有条件；大范围纵横交错的高速铁路分布及高密度运行的高铁列车为高质量地震成像提供了一直以来所需要的均布震源.因此对这种全新的可重复震源——高铁震源所激发出的地震波进行采集并分析，不但能分析高铁列车的运行状态，而且有望对高铁线路路基及近地表结构进行成像.然而，为充分发掘这种全新的高铁震源地震数据所蕴含的信息，采用何种信号分析处理手段进行预处理是极为关键的.直接在时间域观察波形无法获得信号中的频率成分及其频率位置，而常规的傅里叶变换将信号从时间域变到频率域也仅能获得能量沿频率的分布，因此单从时间域或者频率域出发无法刻画频率成分随时间的变化规律；受限于海森伯格不确定性原理的约束，短时傅里叶变换(Qian, 2001; Mallat, 2009)、连续小波变换(Goupillaud, 1984; Mallat，2009)、S变换(Stockwell et al., 1996；高静怀等，2003)等变换不能精确地刻画信号的时频特征；而常用的Cohen类时频分布(Cohen，1989；Baraniuk and Jones, 1993；Choi and Williams, 1989；Jones and Baraniuk, 1995；Zhao et al., 1990)只能对信号进行时频分析，但由于所对应的时频分布缺乏相位信息，无法重构信号.经验模式分解(Huang et al., 1998)可将信号分解为一系列本征模态函数，在勘探地震信号处理中取得了明显的应用效果(Zhou et al., 2012；Li et al., 2014).为了完善EMD的理论体系，Daubechies等(2011)提出了同步挤压小波变换并给出了一套完整的理论，随后在地震勘探等领域取得了诸多应用(Wang et al., 2014; 黄忠来等，2017).若信号是窄带的且振幅的变化率远低于频率的变化率，挤压类时频变换首先计算窄带信号的线性时频变换，然后沿频率方向将变换系数依据变换结果的相位偏导数据进行挤压，可显著地提高时频分布的频率分辨率.由于这种挤压类时频变换的理论基础就假设信号是窄带的，所以它极其适合分析窄带信号.前期的研究成果表明，高铁列车造成的轨道及列车振动呈现出窄带分立谱特性，所产生出的地震波也呈现出窄带特性，因此采用挤压时频分析极其适合分析高铁震源地震信号.因此，本文首次将挤压时频分析引入高铁震源地震信号处理，对北京大学在中国南方某地区采集到的高铁震源地震信号进行应用分析初探，利用挤压时频分析可发现高铁列车的运行状态(匀速、加速等)，并获得了关于高铁震源频谱特征的一些认识.本文主要安排如下：首先简单介绍挤压短时傅里叶变换的理论，然后着重从时间域、频率域分析实测高铁震源地震信号的时间域及频率域特征，最后以挤压短时傅里叶变换为工具对实测高铁震源地震信号进行分析并将其用于判断列车运行状态及重构部分频段信号.

1 挤压短时傅里叶变换

一个时间信号x(t)的短时傅里叶变换可通过如下公式得到：

(1)

式中ω表示角频率，g(t)表示窗函数，[]^*表示对括号内的复数取共轭.根据信号相关的性质，上述过程也可在频率域实现：

(2)

式中X(ξ)为信号x(t)的傅里叶变换，G(ω)为窗函数g(t)的傅里叶变换.若信号x(t)为一复正弦信号e^jω₀t，则信号加窗傅里叶变换STFT(t, ω)与ω₀的关系如下：

(3)

由(3)式可知，若在加窗傅里叶变换STFT(t, ω)不为零时，对(3)式求t的导数，可得到信号加窗傅里叶变换STFT(t, ω)沿时间方向的导数与信号x(t)的真实频率ω₀关系如下(Wu and Zhou, 2018)：

(4)

则构造新的变量 (t, ω)可反映信号的真实频率位置：

(5)

考虑到部分加窗傅里叶变换值为零，采用(5)式可能导致结果不稳定，可在(5)式中采用阻尼因子来提高稳定性.若信号x(t)为e^jω₀(t)t，其中ω₀(t)随时间缓慢变化，(t, ω)亦可反映信号频率成分的变化，则可将STFT(t, ω)挤压到新的时频位置(t, (t, ω))上得到分辨率更高的时频分布.若按照频率间隔Δω对STFT(t, ω)进行采样，且把频率按照间隔Δ划分得到一系列离散化频率_l，则挤压短时傅里叶变换为：

(6)

ε为一控制噪声的阈值，若短时傅里叶变换某些系数小于阈值ε，则不参与挤压.虽然零系数的时间导数并不一定为零，导致利用(5)式计算出来的频率不为零，但考虑(6)式为求和式，零系数对求和式没有任何贡献.若信号为多个窄带且频率成分随时间缓变的分量叠加，只要各个分量在短时傅里叶变换的结果上可以分辨，也可通过上述方法得到高精度的挤压短时傅里叶变换结果.与其他时频分布不同的是，通过挤压短时傅里叶变换的结果可以重构出原始信号：

(7)

若要重构[ω_min, ω_max]频带范围内的信号，可通过限制求和范围来实现：

(8)

2 实际高铁震源地震数据分析

2018年1月，北京大学地球与空间科学学院在中国南方某段高铁线路附近布置了51个三分量低频加速度检波器，检波器频带宽度为0.2~100 Hz，采样间隔为5 ms.我们选取其中的17个检波器进行分析处理，17个检波器(编号为T1至T17)与高铁线路的位置示意图如图 1所示，检波器阵列与高铁线路基本垂直.最近的检波器距离高铁大约75 m，最远检波器距离高铁大约187 m，检波器空间间隔大约为7 m.

2.1 时间域特征

列车1经过时，检波器T1接收到三分量波形如图 2所示.图 2a、图 2b及2c分别为平行高铁线路分量、垂直高铁线路分量以及垂直地表分量的波形.列车经过时波形的振幅强度较大，振幅比无列车通过时高了两个数量级.图 2c所示波形均呈现出周期现象，如图中箭头所示.经过分析，发现此周期波形的周期大约为0.3 s(60个采样点).考虑到我国目前高铁运营速度大约为300 km·h^-1，即83.33 m·s^-1，波形中的周期现象对应的距离为大约25 m，而目前我国大部分高铁车厢长度约为25 m.因此，波形中的周期现象与车厢长度完全吻合.此外，图 2c所示波形中的稳定周期波形大约为16个，而我国高铁编组一般为16节车厢和8节车厢，因此周期波形的个数与车厢个数也是基本一致的.

2.2 频率域特征

列车1经过时，对检波器T1接收到三分量波形进行频域分析，三个分量的振幅谱分别如图 3a、图 3b及图 3c所示.图 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¹⁾)：

1) Liu Y J, Luo Y. 2019. Explore high-speed-train-induced vibration as a new seismic source. Geophysical Research Letters, under review.

(9)

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

(10)

若假定一节车厢长度为25 m，一节车厢四对轮子距该节车厢头部的距离分别为4 m、6.5 m、18.5 m及21 m，则对应A点受力函数的振幅谱如图 4c所示.该振幅谱也呈现出分立谱特征，且谱峰的位置与T1观测数据的振幅谱(图 4b)基本匹配.

2.3 时频域特征

仅从时间域和频率域来分析所接收到的震动信号，不能反映频率成分随时间的变化，也就不能反映列车的运行动态.例如：高铁列车采取加速或者减速操作时，频率成分随着时间是变化的，若采用时频分析能够刻画频率成分随时间的变化，就有可能反映列车的运行状态.考虑到对应的频率成分在短时傅里叶时频谱上基本可分辨且各频率成分基本为窄带的，非常适合采用挤压类变换来精确地获得高铁震动信号频率成分随时间的变化.采用挤压短时傅里叶变换计算检波器T1所接收到的列车1激发的震动信号，得到的三分量的时频谱分别如图 5a、图 5b及图 5c所示.由图 5可看出，经过挤压操作以后窄带信号的特征更加明显，频率成分随时间变化的规律更加清楚.由于列车经过时各分量的频率成分基本保持不变，说明列车在经过T1检波器是匀速运行的.由于挤压时频变换可以很精确地得到各频率成分的位置，我们还可以重构特定频率的分量.我们将列车1经过时各检波器收到的信号选择30 Hz进行重构，原始的平行高铁分量和30 Hz重构的平行高铁分量分别如图 6a和图 6b所示.从原始平行高铁分量(图 6a)中并不容易观察出波形的空间特征，但从30 Hz重构的平行高铁分量(图 6b)中可非常容易看到各检波器接收到信号的时差，能够清晰地观察出波形的空间特征.

我们选取了另外一趟高铁列车(列车2)经过时所收到的信号进行挤压时频分析，检波器T1接收到的三个分量的挤压时频谱图分别如7a、图 7b和图 7c所示.与列车1经过时各分量的挤压时频谱图(图 5)相比，列车2经过时挤压时频谱图中各窄带成分的频率随着时间的增大也逐渐增大，说明此时列车的运行状态有所改变.由式(10)可知，若速度增大导致各脉冲函数之间的间隔缩小，会导致频率升高.因此，由图 7可知，列车2在经过检波器阵列时明显处于加速状态.

3 结论

经过对中国南方某地区实际采集到的高铁震源地震信号进行时域和频域分析，发现高铁震源引起的地震信号呈现出窄带分立谱特征.我们首次将挤压时频变换这种极其适合分析窄带信号的时频分析工具引入到高铁震源地震信号分析中，并对实际采集到的高铁震源地震信号进行初步分析，得到了以下认识：(1)距离高铁线路较近的检波器所收到的震动信号具有周期特征，周期与车厢长度及列车运行速度有关；(2)高铁列车所引起的震动信号呈现出窄带分立谱特性，分立谱的形状与车厢参数(车厢长度、轮子间距等)及列车运行速度密切相关；(3)基于高铁列车基本具有相同的结构及经过某一路段时具有基本相同的运行速度，高铁列车所引起的震动具有很强的可重复性；(4)挤压时频变换极其适合分析高铁震源地震信号，利用挤压时频分析不但可判断高铁列车的运行状态(匀速或加速)而且可重构出某些感兴趣频带的信号.关于挤压时频变换在高铁震源地震信号处理中进一步应用，例如监测高铁安全隐患及探测地下结构等，还有待深入研究.

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