(1)摘要 基于桥墩均匀分布且高铁列车匀速行驶情况下, 每个桥墩源的震源时间函数仅存在时间延迟差异的假设, 分别通过构建均匀空间声波模型以及半无限空间弹性波模型, 计算格林函数, 给出反演所用的线性方程组。利用最小二乘法对高铁波场的等效震源时间函数进行反演, 测试在一定噪声水平下的反演效果, 探讨各种因素对反演结果的影响。
关键词 高铁; 震源时间函数; 反演; 格林函数
中国是世界上高铁线网规模最大的国家, 高铁列车的运行速度居世界前列。高速行驶的列车会造成车厢以及铁轨、桥墩的振动。高铁列车的振动可用于高铁故障诊断[1-4]。对高铁列车引起的桥体及其下路基振动的检测, 可用于桥体探伤或检测路基变化[5-7]。
除车身和桥梁路基的振动外, 该振动还会向更远处传播, 从而可以用作反演了解地下结构的震源。Chen等[8]对普通铁路列车激发的地震波进行研究，认为向下传播的地震波甚至有可能穿透地壳, 带来深层结构信息。徐善辉等[9]对京津城际铁路进行较长时间的观测, 发现高铁地震波场能够在几千米之外被检测到。在勘探地震学中, 利用单炮激发的多道地震记录进行偏移成像得到地下主要速度界面结构的理论和应用方法已经成熟。如果希望利用高铁激发的地震波进行地下结构成像, 则首先需要解决高铁震源并非点源的问题。可以将高铁地震震源视为若干有限个点源叠加(比如将桥墩视为震源), 或是无穷个点源叠加(比如将高铁列车或铁轨视为连续的源)。无论哪种震源模型, 都需要考虑高铁列车与铁轨、桥梁的耦合。若考虑铁轨上某一点受列车车轮的作用, 可以给出由若干个冲激函数组成的震源时间函数[10], 也有研究者通过正演给出列车经过时桥梁和路基的振动信息[11-13]。
本文探讨如何把天然地震震源的方法应用于高铁震源研究。天然地震震源研究主要通过波形反演震源机制和震源破裂过程。本文利用波形反演方法, 分别在声波介质以及半无限均匀弹性介质假设下, 给出反演高铁震源时间函数的方法，并进行理论测试。
在天然地震研究中, 震源时间函数的反演是反演地震破裂过程, 即震源破裂滑动的位移或速率随时间的变化。通常有两种描述震源的方法: 点源和有限震源。点源适于描述较小的地震, 研究大地震的破裂过程则通常需要将发震断层划分为多个子断层, 将子断层视为点源, 再反演每个点源的震源时间函数。在研究高铁地震问题中, 实际的震源是在地表或近地表一条很长的线上移动的, 涉及移动源叠加问题。大多数高铁线路都是架设在高架桥上的, 可以近似地将高铁桥梁的每一个桥墩视为点源。桥墩的宽度一般在 3m 左右, 在台站距离高铁线路超过 20m 的情况下, 满足空间点源近似条件。同时, 只要通过滤波保证地震波的波长远远大于桥墩的横向尺度, 也能满足时间点源近似条件。
反演的基础是正演, 可以将地震台站的记录视为震源时间函数与格林函数卷积的结果。因此, 格林函数在反演过程中尤其重要。对于近震的震源时间函数, 在反演中常用不同的方法[14-15]计算一维层状介质中的格林函数。考虑到桥墩源与观测台站距离较短, 浅层沉积层结构复杂且波速未知, 先采取较简单的介质情况, 分别是三维均匀空间声波方程的格林函数和半无限均匀弹性空间格林函数[16]。
对大地震的震源破裂过程进行有限震源模型反演时, 各子断层的震源时间函数甚至滑动角都可以不一致。对高铁震源做反演时, 给出如下假定。首先, 认为同一列高铁在通过反演涉及的所有桥墩的过程中匀速行驶, 这与大多数高铁线路的实际情况基本上一致。第二, 认为桥墩均匀分布, 也基本上符合实际情况。基于以上两个假定, 可以认为每个桥墩源的震源时间函数形式一致, 只是有一个由于高铁列车行驶到达时间差异引起的延时。反演得到的最终结果是高铁列车经过时的某一个桥墩处的震源时间函数, 以此作为高铁波场的等效震源时间函数。对于弹性波动方程, 根据相应格林函数的定义, 得到的结果是桥墩震源三分量力的震源时间函数。
求解高铁震源时间函数是反褶积问题:
(1)其中
是位于 m 处的观测点在 t 时刻 m方向的位移; K 为桥墩点源总数; 格林函数
为
处的桥墩点源在
时刻 n 方向的集中脉冲力作用下, 在观测点 x 处产生的位移; 
为桥墩点源 n方向力的震源时间函数。该式可矩阵化为
(2)
格林函数卷积矩阵 G 是每一个桥墩源的矩阵 Gi 的叠加:
。 (3)其中,
。(4)
该矩阵的行数 L 为台站记录地震信号的长度, 在实际计算中, 该矩阵的列数取震源时间函数的有效长度部分即可。
研究天然地震震源的过程中, 震源机制解反演通常需要利用多个台站的记录同时约束, 震源破裂过程反演则可以通过单台记录完成, 给出该台站的视震源时间函数。本文对高铁震源时间函数的反演也将从单台反演开始, 拓展到多台联合反演。
对式(2)给出的这一线性问题, 可以采用不同的反演方法。考虑到一般情况下震源时间函数的有效长度远小于记录的长度, 所以这是一个超定或混定问题, 可以采用最小二乘法直接求解。
计算中对时空坐标作如下定义。将 x 轴设置在高铁线路上, 并在原点处设置一个桥墩, 则各桥墩点源的位置得以确定。以高铁行驶方向为 x 轴正方向, 以垂直高铁线路的水平方向为 y 轴, 竖直向下为 z 轴, 且 x, y 和 z 满足右手系。在时间方面, 把高铁车头通过坐标原点的时刻记为t=0。
在声速为 c 的三维均匀空间中, 理想流体小振幅振动声波方程为
, (5)其中, p 为声压场, s 表示声源。对应的格林函数
满足
。(6)
由拉普拉斯变换可求出格林函数[17]为
。 (7)由于反演问题为线性, 且独立线性方程组方程数远大于未知数, 因此若利用理论计算得到的台站记录进行震源时间函数的反演, 则会得到与正演完全相同的震源时间函数。因此, 本文主要分析记录的信噪比以及格林函数的准确性对反演的影响。
给定作为参考模型的震源时间函数为正弦函数, 如图 1(a)所示。取高铁行驶速度 v=80m/s, 声速c=300 m/s, 计算 x1=(0, 100 m, 0)处台站的正演波形, 结果如图 1(b)所示。然后, 在得到的正演结果中加入 10%的随机噪声, 结果如图 1(c)所示。采用含噪声的记录进行最小二乘反演, 得到的震源时间函数如图 1(d)所示。可以看出, 添加噪声之后, 反演得到的震源时间函数与原震源时间函数较接近, 主要在振动幅度上有所差异, 同时产生部分较高频的振动成分。
下面讨论格林函数对反演结果的影响。在考虑声波方程的简单情况中, 格林函数主要受两个参数的影响: 波速和每个桥墩源的高铁列车到时。由于假设高铁匀速行驶, 因此高铁到时实际上由高铁列车行驶速度决定。
采用图 1(a)所示的震源时间函数, 在生成测试数据的过程中分别将车速取 0.5%, 1%和 2%的误差, 相应的最小二乘反演结果见图 2。可见, 速度偏差对震源时间函数的反演结果有较大影响, 主要体现在振幅差异上, 最大相对误差可超过 100%。在实际的高铁观测过程中, 通过两种方法计算车速: 1)利用高铁沿线上距离不太远的两个桥墩旁边的记录, 分析高铁到达时间, 得到平均车速; 2)通过分析近处台站记录的频谱, 得到高铁速度信息[18]。第一种方法计算的是平均速度, 只要高铁的到时测定精准, 得到的平均速度就比较准确, 但只能反映高铁在一定区间内行驶的平均情况。利用第二种方法计算车速时, 需要已知高铁列车车厢长度参数, 然后通过分析频谱得到非常精确的车速。
作为真实情况的高度简化, 声波方程能够帮助简化计算, 避免多分量的讨论, 便于直观地分析影响反演结果的因素。但是, 考虑到将来的实际应用, 则需要计算弹性波动方程。
波速对反演结果的影响如图 3 所示。可以看出,在异常程度相同(2%)的情况下, 波速异常对反演结果的影响远小于车速异常。波速的不准确会影响波传播的时间, 从而明显地影响震源时间函数的相位, 这在图 3(b)和(c)中均有体现。
(a)简谐震源时间函数, (b)声波方程正演结果, (c)为(b)中的记录添加 10%随机噪声后的结果, (d)震源时间函数反演结果
图1 利用简谐震源时间函数的正演波形以及添加噪声后反演得到的震源时间函数
Fig. 1 Waveform calculated from simple harmonic source time function and the inversion result obtained from data with noise
灰色虚线为原震源时间函数，黑色实线为反演结果，下同
图2 车速误差对震源时间函数反演结果的影响
Fig. 2 Influence of train speed error on the inversion results of source time function
图3 波速误差对震源时间函数反演结果的影响
Fig. 3 Effect of wave velocity error on the inversion result of source time function
半无限均匀空间格林函数计算使用 Johnson[16]给出的公式。该公式计算源在地表的格林函数时, 会产生振幅无穷大的瑞利波, 同时也会遇到数值问题。所以本文计算时认为桥墩源有一定的深度, 这与实际情况相符。
在反演过程中, 需要利用三分量的记录反演三分量的震源时间函数。记三分量记录为 U1, U2 和U3, 三分量震源时间函数为s1, s2和s3, 则反演公式[U]=[G][s]中 3 个矩阵的表达式变为
。 (8)由于矩阵比较复杂, 所以我们尝试通过单分量记录来反演三分量震源时间函数。合成地震记录所用的三分量震源时间函数为波形相同、振幅不同的正弦函数, 模拟高铁行驶过程中桥墩的真实受力, 即竖直的 z方向受力最大, 垂直于高铁的 y 方向受力最小。对合成记录添加 10%的随机噪声, 分别利用竖直分量、平行高铁分量和三分量记录同时反演震源时间函数, 结果如图 4 所示。结果表明, 由于水平分量的震源时间函数幅值远小于竖直分量, 若不添加新的约束, 在反演过程中将被竖直分量主导, 水平分量无法得到可靠的结果。由于实际计算波形的竖直分量能量大于水平分量, 因此利用竖直分量反演的结果更为理想。利用三分量同时反演的结果则非常接近原震源时间函数, 约束较好。
从地震记录保留的震源信息来看, 近场地震记录通常包含震源的高频信息, 而远场记录低频能量较高, 二者可分别反映震源时间函数的宏观特征和细节特征。因此, 天然地震震源研究中, 想要得到更精确的震源破裂过程, 通常采用近场强地面运动记录和远震记录联合反演[19]。
第1行: 利用平行高铁分量记录反演结果; 第2行: 利用竖直分量记录反演结果; 第3行: 利用三分量记录反演结果
图4 使用单分量或三分量记录反演三分量震源时间函数的结果
Fig. 4 Inversion of three-component source time function by using single component or three-component records
以左上标表示不同台站, 则N个台站联合反演的公式为
。 (9)利用带噪声的合成数据, 可以比较利用单个台站和多台联合反演得到的震源时间函数与正演所用震源时间函数的接近程度。为了计算方便, 采用声波的反演模型。单台站、两个台站和 5 个台站的反演结果如图 5 所示, 所有合成记录均添加了相同幅值的随机噪声。可以看出, 两个台站的约束效果明显好于单台。取 5 个台站时, 反演结果只比两个台站稍好, 但涉及的矩阵较大, 计算效率较低。
在大多数情况下, 可将高铁震源视为由若干个桥墩点源组成。在桥墩均匀分布且高铁匀速行驶的情况下, 近似地认为各桥墩的震源时间函数仅存在时间延时的差异。将各桥墩源高铁到时作为参数
, 代入格林函数中, 得到可直接用于反演震源时间函数的公式, 并用最小二乘法求解该线性问题, 获得高铁波场等效震源时间函数, 本文得到如下主要结论。
1)在基于声波模型的反演过程中, 由于列车行驶速度影响其到达每一个桥墩的时间, 因此对列车速度的微小估计偏差会造成反演结果的较大误差; 介质波速对反演结果的影响则比较小。
2)在考虑半无限空间的情况下, 可以使用单分量记录或多分量记录反演三分量的震源时间函数。若某一分量的震源时间函数振幅远大于其他分量, 则另外的分量无法得到正确的反演结果。采用单分量反演时, 该分量的信噪比越高, 反演效果越好, 三分量综合反演给出的结果最为理想。由于记录是三分量震源时间函数共同产生的, 无法单独对某一分量进行反演。
图5 单台反演以及多台联合反演结果
Fig. 5 Single inversion and multi-joint inversion results
3)多台站联合反演结果比单台站结果更可靠, 通常两个台站对震源时间函数就有较好的约束。若使用更多的台站数据, 则需要考虑计算效率。
参考文献
[1]赵晶晶, 杨燕, 李天瑞, 等. 基于近似熵及 EMD 的高铁故障诊断. 计算机科学, 2014, 41(1): 91-94
[2]李贵兵, 金炜东, 蒋鹏, 等. 面向大规模监测数据的高铁故障诊断技术研究. 系统仿真学报, 2014, 26(10): 2458-2464
[3]李智敏, 苟先太, 秦娜, 等. 高速列车振动监测信号的频率特征. 仪表技术与传感器, 2015(5): 99-103
[4]朱菲, 金炜东. 基本概率指派生成方法在高铁设备故障诊断中的应用研究. 中国铁路, 2018(4): 82-86
[5]姚京川, 杨宜谦, 王澜. 基于 Hilbert-Huang 变换的桥梁损伤预警. 中国铁道科学, 2010, 31(4): 46-52
[6]丁幼亮, 王超, 王景全, 等. 高速铁路钢桁拱桥吊杆振动长期监测与分析. 东南大学学报(自然科学版), 2016, 46(4): 848-852
[7]赵瀚玮, 丁幼亮, 李爱群, 等. 大跨多线高速铁路钢桁拱桥车——桥振动安全预警研究. 中国铁道科学, 2018, 39(2): 28-36
[8]Chen Q F, Li L, Li G, et al. Seismic features if vibration induces by train. Acta Seismologica Sinica, 2004, 17(6): 715-724
[9]徐善辉, 郭建, 李培培, 等. 京津高铁列车运行引起的地表振动观测与分析. 地球物理学进展, 2017, 32(1): 432-422
[10]Hung H H, Yang Y B. Elastic waves in visco-elastic half-space generated by various vehicle loads. Soil Dynamics and Earthquake Engineering, 2001, 21(1): 1-17
[11]和振兴, 翟婉明, 杨吉忠, 等. 铁路交通地面振动的列车-轨道-地基耦合数值方法研究. 振动工程学报, 2008, 21(5): 488-492
[12]Cheng Y S, Au F T K, Cheung Y K. Vibration of railway bridges under a moving train by using bridge-track-vehicle element. Engineering Structures, 2001, 23(12): 1597-1606
[13]Feng S J, Zhang X L , Zheng Q T , et al. Simulation and mitigation analysis of ground vibrations induced by high-speed train with three dimensional FEM. Soil Dynamics and Earthquake Engineering, 2017, 94: 204-214
[14]Bouchon M. A simple method to calculate Green’s function for elastic layered media. Bull Seism Soc Am, 1981, 71: 959- 971
[15]Kohketsu K. The extended reflectivity method for synthetic near-field seismograms. J Phys Earth, 1985, 33: 121-131
[16]Johnson L R. Green’s function for Lamb’s problem. Geophysical Journal International, 1974, 37(1): 99-131
[17]Aki K, Richards P G. Quantitative seismology. 2nd ed. Mill Valley, CA: University Science Books, 2002
[18]王晓凯, 陈文超, 温景充, 等. 高铁震源地震信号的挤压时频分析应用. 地球物理学报, 2019, 62(2): 2328-2335
[19]Yagi Y, Mikumo T, Pacheco J, et a1. Source rupture of the Tecoman, Colima, Mexico earthquake of 22 January 2003, determined by joint inversion of teleseismic body-wave and near-source data．Bull SeismSoc Am, 2004, 94(5): 1795-1807
Inversion of Effective Source Time Function of the High-Speed-Rail Wavefield
Abstract Based on the assumption that there is only time delay difference in the source time function of each pier when the piersare uniformly distributed and the high-speed train moves at a uniform speed, we calculate Green’s function by constructing uniform spatial acoustic wave model and semi-infinite spatial elastic wave model, respectively, and give the linear equations used for inversion. Thenthe least square method is used to invert the effective source time function of the high-speed rail. For checking availability of the inversion method, the effect of inversion at a certain noise level is tested, and the influence of various factors on the inversion results is discussed.
Key words high-speed train; source time function; inversion; Green’s function
doi: 10.13209/j.0479-8023.2019.075
收稿日期: 2019-07-29;
修回日期: 2019-08-14
国家自然科学基金(41874071)资助