基于短时互相关的地表横波波速研究

靳伯骜，高艳

（1.中国科学院噪声与振动重点实验室，北京 100190；2.中国科学院大学，北京 100049）

1 引言

随着城市的快速发展，城市地下空间的开发与利用受到了越来越多的关注。城市管线定位以及地铁轨道噪声相关研究中使用了弹性波理论，通过实验或仿真手段，对声源或反射体进行定位[1～3]。其中，地表波速是一个关键参数，其精度直接影响定位的精度[4]。同时，工程中也可利用纵波与横波波速确定土壤的等效体积模量以及剪切模量等相关参数。但由于城市环境和地下土况十分复杂，测试结果受周围环境和多模式声波影响很大。实际的波速计算中需要鲁棒性较强的算法。现有的计算方法有传递函数相位法和互相关方法，通过对两组一维时间信号进行处理，得到两个信号之间的时间延迟。

传递函数相位法和传统互相关方法适用于平稳信号的处理，对于时变或突发信号而言，需将一维时域信号映射到时频域进行处理。关于时频域的算法要追溯到20世纪60年代，语音编解码技术中已开始使用短时傅里叶变换[5、6]。随着计算技术的发展，时频域分析技术也逐渐演进[7、8]，先进的小波技术已应用于提取和计算湍流相干结构的研究[9]。对于互相关算法，在快速傅里叶变换之前已被大量使用，是一种十分常见的计算操作，目前已被广泛应用于计算时间延迟相关的研究中[10]。

本文提出利用短时加窗的互相关方法求取浅表横波波速。使用两个传感器进行水平横波的实验，分别通过传递函数梯度法、传统互相关法以及短时互相关方法计算地震波波速。实验发现地震波信号能量集中在55～250Hz内；传递函数梯度法求取的波速随频带选取变化；短时互相关方法比传统互相关方法的结果更接近传递函数梯度法中线性区间段的结果。

2 理论

现有的频率、时频域或尺度变换都是将时域信号通过一系列基函数变换到所对应的域中。例如，对于短时傅里叶变换有：

式①中：i为虚数单位，ω为角频率，t和τ分别代表时间和时间延迟，是时域信号，是时间窗函数，是短时傅里叶变换结果。由式①可发现，时频域变换通过两个算子和，可将一维的时域数据转到二维的时频域中。类似地，在小波变换中，可通过变化两个参数来计算时间尺度域的结果。

短时互相关算法与短时傅里叶变换和小波变换，有很多相似的地方。假设测量的两路信号为和，仿照短时傅里叶变换以及小波变换，短时互相关的计算表达式为：

式③中：δ代表时间延迟，τ代表时间窗起点。短时互相关通过移动时间窗函数截取信号片段与进行互相关，旨在对空间中多波速、多波数信号进行匹配，求解信号的时间延迟。基于以上算法进行数值仿真。仿真中假设信号传递路径是线性的，两路信号之间相差一个时间延迟，同时信噪比较高，忽略了噪声影响。仿真中的采样频率设为4096Hz，与分别使用了中心频率为200Hz的雷克子波（又称草帽函数），使用了32点的Hann窗作为短时窗函数，两个信号之间的延迟为38个数据点，对应的时间延迟为9.033ms。短时互相关的结果如图1所示。

图1 数值仿真中两个雷克子波的短时互相关结果

图1的灰度图中显示的是短时互相关结果，横轴是时间窗函数的起点时刻，纵轴对应的是互相关的时间延迟，灰度表示的是相应坐标下的互相关值；图1也绘制了时域的波形图，经过放缩与平移后绘在灰度图之上，其时间轴与灰度图的时间窗起点轴一致。由图1可知，两个雷克子波的短时互相关的峰值在实线峰值之前，表示该点的时间窗正好包含了实线信号的大部分能量，从而导致互相关得到最大峰值。时间延迟的坐标点与生成信号时输入的时间延迟一致。仿真结果表明，理想情况下短时互相关算法能够准确确定两个雷克子波之间的时间延迟。

3 地表横波波速实验

实验选址在中科院声学研究所院内的绿化带中，测试现场与周围建筑物、柏油路面以及高大乔木的距离大于3m，地表附有稀疏草本植被。在场地中心安装耙状声源底座，底座激励的方向平行于南北经线方向。在底座西侧沿东西纬线方向布放两个三轴加速度计（LC0161-2），加速度计、底座三点共线，间距1m。加速度计通过螺栓固定于特制金属底座上，底座有三个长为5cm的齿状物埋于地下，加速度计与地表紧密连接。加速度计的x轴、y轴、z轴分别平行于经线、纬线和铅垂线。测试时，为避免地表负载的影响，负责锤击操作的实验员站在底座上，数据采集设备通过信号线在5m外的柏油路面上进行，采样频率设为32 768Hz。由于手动锤击操作容易发生双击，并且受到现场多因素影响，一些文献中提出多次激励，对统计平均后的数据进行处理。鉴此，先后进行了10次锤击。

对采集到的时域信号进行初步处理，设置触发阈值为3m/s2，对1m处的加速度计的x轴数据进行逐点扫描。当满足触发条件时，判断此时刻为锤击起始点。截取每一个锤击起始点前3000点至起始点后12 000点的数据。通过数据的预处理剔除了存在明显差异的2组数据。

4 数据处理与分析

大多数的地表波速测量方法都是使用已知间距的传感器阵列，通过测量传感器之间的地震波到达时间差进行波速估计。通过对实验数据的预处理得到了8组锤击实验的数据，每组中包含两个三轴加速度计，共计六个通道的数据，每通道都有15 000个点，持续时间接近0.5s。由离散傅里叶变换特性可知，这些数据的频域分辨率接近2Hz。截取锤击脉冲片段的同时，也截取了每段脉冲之后相同长度的背景噪声片段。

首先，对每一段数据进行傅里叶变换，截取1kHz以内的频谱，幅值取对数后绘制于图2。图2中绘制了六个通道脉冲以及噪声的频域数据。其中（a）、（b）、（c）三个子图分对应近端传感器x轴、y轴、z轴三个通道，（d）、（e）、（f）对应远端传感器的三个通道；每个子图中都绘制了8次锤击实验的结果，其中实线是脉冲实验的数据，点图是脉冲后静息时的背景噪声。对比六个通道的数据发现，土壤中的地震波信号在55Hz至250Hz的频段内信噪较高。本文使用半透明灰色矩形窗对该频段进行了标记。

图2 两个三轴加速度计的脉冲与背景噪声片段的频谱

4.1 传递函数相位法的结果

由频域分析发现，8组锤击实验的一致性较好。由于锤击脉冲峰值幅度接近，且使用同一个判断阈值，故假设每组锤击实验起始点时刻的差异可忽略。对8组锤击实验结果求取平均，得到平均的锤击实验的时域波形。计算平均波形中x轴信号的传递函数相位，将其绘制于图3。

图3 平均时域信号的x轴互相关相位图以及三个频率区间段内计算得到的波速

从图3可发现，传递函数相位在55Hz至170Hz频段内呈线性变化。对照图2发现，该频段与地震波峰值频段对应。此外，在55Hz至250Hz内，传递函数相位梯度存在变化。实验中两个三轴加速度计的间距为1m，依此计算了三个频率区间段的等效波速，分别约为116.9m/s、82.2m/s、46.0m/s。初步分析猜测，波速变化的一种原因可能是由于地表土壤分层情况复杂，在200Hz附近存在水平横波（SH波）模式之外的波，造成波速结果偏移；另一种可能是埋地情况复杂，传递路径上存在中心频率接近200Hz的共振结构，导致两个传感器之间的相位发生偏移。

通过对锤击信号进行频谱分析，从图2可发现，x轴的信号的能量高于其余两轴。考虑到加速度计安装时存在倾斜误差，初步判定实验中SH波模式的能量占主导。同时，图3的线性频段与图2的峰值频段较为吻合，均在55Hz至170Hz之间。由上述分析初步判定，实验中的水平锤击主要激发出了SH波，能量集中在55Hz至170Hz之间，且波速约为116.9m/s。

4.2 传统互相关的结果

传统互相关方法通常用于描述两个时间序列之间的相关程度。在传递路径分析中，当激励信号为宽带脉冲或白噪声时，互相关峰值所对应的时间延迟即为信号到达两个传感器的时间差。

本文计算了三个轴向的8次锤击实验的互相关函数，同时计算了x轴时域平均波形的互相关，绘制于图4。由于每次锤击的力度不同，同一个轴向的结果不同锤击的幅度偏差较明显，不易于观察，因此图4对x轴、y轴、z轴数据分别进行了能量归一化处理。

图4 x轴、y轴、z三轴的传统互相关结果以及x轴时域平均波形的互相关结果

由图4可发现，x轴互相关峰值远高于z轴互相关峰值，且x轴与z轴互相关零点、峰值点位置十分一致。初步判定传感器安装时存在俯仰角方向的偏差，z轴加速度计采集到了水平方向的分量。对于压缩波、Rayleigh波以及竖直横波（SV波）模式而言，均存在竖直方向的运动分量，而实验结果中竖直振动十分微弱，因此可推测实际测试中主要激发出了SH模式的地震波，与频谱分析的结论一致。图4实线为x轴时域平均波形的互相关，通过时间延迟以及已知的传感器间距计算得到地震波波速约为110.0m/s，与传递函数相位法得到的波速相差接近6%。

4.3 短时互相关的结果

类似于短时傅里叶变换，短时互相关算法使用了时频域局域化的时间窗函数对数据进行截断。实际处理时，首先对x轴时域平均波形进行降采样至4096Hz。短时截断时，时间窗函数设为256点的哈恩窗（Hann Window）。时间窗起点从锤击信号数据片段起点开始，步长为6，计算每个时间窗起点所对应的互相关函数，绘制于图5。

图5 x轴时域平均波形的短时互相关结果

图5的灰度图中，每一条平行于y轴的灰度数据对应一个时刻的互相关，纵轴坐标即为互相关的时间延迟。得到短时互相关的灰度图后，又计算了每一个时间窗起点下的互相关峰值所对应的时间延迟，用黑色圆点绘制于图5中。可发现，在时间窗起点为0.04s至0.1s的区间段内，互相关峰值较为明显，峰值坐标较为一致。这是因为在这个区间段内，时间窗能将锤击信号的大部分能量包含进去，使得这个区间内的互相关峰值较高。求取所有区间段的峰值中的最大值，用黑色叉号标记于图5中。记录该最大值的时间延迟，得到短时互相关方法下的地震波波速约为113.8m/s。该结果介于传递函数相位法和传统互相关方法的结果。

对比三种方法求取互相关的过程发现，传递函数相位法受到地表复杂性以及所选取频率区间的影响，其结果存疑；传统互相关方法中，8次锤击实验的结果一致性较差，x轴时域平均波形的结果与时域数据的结果存在不一致，有待进一步研究；短时互相关方法物理图景较为清晰，能够分析多模式、时变信号，但对于其窗函数选取以及算法鲁棒性等问题仍有待进一步研究。

5 结论

通过水平锤击声源在地表激发水平横波，使用两个三轴加速度计对地表振动进行测量。通过数据分析初步判定，该方法所产生的地震波主要以SH波为主。通过传递函数相位方法、传统互相关法以及短时互相关方法计算地震波波速，分别约为116.9m/s、110.0m/s及113.8m/s。对结果进行分析发现，实验现场的地表水平横波波速在110m/s至120m/s附近，但更高精度的结果仍有待进一步研究。同时，短时互相关方法拥有分析多模式、时变信号的能力，物理含义清晰，存在进一步研究的价值。