计及非稳态效应的地脉动功率谱密度包络修正及应用

2019-12-31陈学前范宣华王玉军胡绍全

振动与冲击 2019年24期

郝雨，陈学前，范宣华，王玉军，杜强，胡绍全

(中国工程物理研究院总体工程研究所，四川绵阳 621999)

地脉动作为一种随机载荷，对很多大型精密科学装置的稳定性会产生显著影响[1-2]，为评估地脉动对科学装置建筑地基及结构设计的影响，需要获得相应的地脉动载荷谱曲线[3-4]，而地脉动的测试数据一般是时域信号，需要对信号进行甄别和相应处理，得到地脉动载荷的功率谱密度包络，为建筑基础及结构载荷输入提供参考。

传统的功率谱密度计算一般针对平稳信号进行[5]。但是，地脉动信号是一个典型的非平稳随机过程，其功率谱密度是一个随时间变化的函数，传统的功率谱密度计算方法只能够反映一段时间内信号的平均能量特性，如此会使结构偏于危险，而且，功率谱密度计算的时间长度选取会对计算结果产生较大影响。

学者们针对非平稳随机信号的功率谱密度也展开了一定的研究。例如，王良曦等[6]应用Wiener过程模拟了不平路面的功率谱密度；孔凡等[7]应用局部平稳小波过程，提出了一种非平稳随机过程的演变功率谱密度估计方法。这两种方法虽然具有一定的通用性，并且能够重复体现功率谱密度曲线的时间变化历程，但是对于大量、长时非平稳时间序列而言，所需的时间和计算空间过大，难以在实践中应用。史康等[8]通过将风速时间序列分解为时变确定性趋势成分和平稳零均值随机成分，分析了实测风场的非平稳特性。但是，这种分解方法是针对某一类问题得到的，难以推广应用于地脉动信号中。林建生等[9]分析了泉州地震台的地脉动数据，给出了若干时段的功率谱密度，但是其中没有考虑功率谱密度的非平稳性。

本文以绵阳地震台台站的地脉动监测数据为例，分析了地脉动信号的数值特征,在此基础上提出了一种考虑地脉动非平稳特性的功率谱密度修正算法，计算了绵阳地区地脉动信号的功率谱密度包络。该算法消除了PSD(Power Spectral Denaity)计算中人为参数选取对计算结果的影响，同时避免了平均效应导致的地脉动PSD幅值估计偏低的现象。

1 研究对象

选取2015年—2016年绵阳地区花荄台站的数字地震观测资料作为研究对象，数据包括垂直(Z向)、东西(E向)、南北(N向)3个方向的地脉动数据。时间间隔为0.01 s，数据单位为counts/s，原始测量数据与加速度信号之间通过如下关系式联系：

时域关系式为

(1)

根据式(1)，原始测量数据与加速度信号的功率谱之间的关系式为

(2)

式中：x(t)为原始速度测量信号；a(t)为转换得到的加速度信号，g；Gx(f)和Ga(f)分别为它们的功率谱密度；η为测量灵敏度。

分析时一般选取无地震或其他干扰的时段[9]，去除烈度大于等于Ⅱ级的有感地震及非地震干扰引起的异常信号。具体方法如下：对第一类异常信号，根据文献[10]，烈度为Ⅴ的地震峰值速度为0.02～0.04 m/s，且烈度每增加一级，地震峰值速度约增加1倍，选取峰值速度2.5 mm/s作为无感地震和有感地震的分界线。计算中,对数据按小时分段，如果某段数据的最大速度峰值大于2.5 mm/s，则去掉整段数据。

第二类异常信号的特点是振动幅值明显高于常态地脉动，同时相比于无感地震信号，具有持续时间长的特点。根据地震局的测量数据，无感地震的持续时间一般在几秒甚至更短的数量级[11]。以每10 s为一个窗，如果连续3个窗内有超过5%的数据点大于该小时或当天数据均方根值的5倍，则去掉整段数据。事实上，由于非地震干扰在量级上和持续时间上与其他信号有明显差异，判断阈值的选取具有较好的鲁棒性。

2 地脉动信号的功率谱密度包络

2.1 平稳信号的功率谱密度

平稳随机信号的功率谱密度与时间无关，对零均值离散平稳时间序列xn=x(nΔt)，其单边功率谱密度

(3)

式中：Δt为离散时间序列的时间步长；N为信号总长度。

为避免产生过多的毛刺，通常采用滑动平均的办法，求解一段相对较长时间T内的平均功率谱密度。由于平稳随机信号的功率谱密度与时间无关，时间长度T的选择不影响功率谱密度计算的结果。

2.2 地脉动信号的非平稳特性

地脉动信号是一个典型的非平稳随机过程。应用自相关函数(Augmented Dickey-Fuller,ADF)法[12]和游程法[13]对一段地脉动信号进行了平稳性检验，结果表明：地脉动信号在从几十秒～1天的任意一段时间长度内都是非平稳的。经典的做法是将地脉动信号视作局部平稳的，在每时间长度T内视作平稳随机过程，求出一条功率谱密度曲线，再将所有的曲线合并取包络。这样，时间长度T的选择会对结果产生较大影响，时间长度T越大，求出的功率谱密度包络幅值越低，原因在于振动量级较高的时间段会被振动量级较低的时间段平均。但是，时间长度T选择过小时，又会影响功率谱密度的求解精度和频率分辨率。

以花荄监测点E方向2016-03-01的地脉动信号为例，当时间长度取不同值时，得到的PSD包络曲线如图1所示。以不同时间长度计算出的功率谱密度谱形基本一致，但是幅值不同，T=30 s时的功率谱密度曲线比T=10 min时高约10 dB，并且没有发现明显的收敛趋势。这提示我们必须考虑地脉动信号的非平稳特性，对经典的功率谱密度包络求解方法进行修正。

图1 时间长度T对功率谱密度包络的影响Fig.1 Effect of time length T to the PSD envelope

2.3 地脉动信号的功率谱密度修正

在地震分析中，通常可以把非平稳地震信号近似为一个平稳随机过程和一个反映振幅强度随时间变化的函数的乘积[14]，即

x(t)=A(t)z(t)

(4)

式中：x(t)为非平稳地震信号；z(t)为平稳随机过程；A(t)为反映振幅强度的量(见图2)。

图2 非平稳地脉动信号的分解Fig.2 Decomposition of non-stationary ground microtremor signal

由于A(t)相对于z(t)是缓变的，式(4)可以写成频域的形式

Gx(t,f)=A(t)·Gz(f)

(5)

式中：Gx和Gz分别为x(t)和z(t)的功率谱密度，由于z(t)是平稳随机过程，Gz与时间无关。图1中不同曲线的形状相似性验证了这一近似的合理性。

将一段时间T分成若干子区间，根据以上假设，在每一子区间内，其功率谱密度是成比例的，有

Gi(f)=AiGz(f)

(6)

式中：Gi为第i个子区间内的地脉动信号功率谱密度；Ai在每个子区间内是一个常数。

考虑到功率谱密度的积分与时域信号均方根值的关系，式(6)可以写成

(7)

式中：RMS(xi)和RMS(x)分别为第i个子区间内和整段地脉动信号的均方根值。功率谱密度的包络即可写为

(8)

式中：Env为该段信号的功率谱密度包络。

为防止高振幅信号被低振幅信号平均，导致遗漏危险载荷，子区间时间长度应当尽可能小。考虑地脉动信号相邻两次穿越0值之间的部分，称之为一个游程，以每一个游程作为一个子区间，这可以作为子区间的最细划分。在一个游程内，地脉动时域信号为一个半波，如图3所示。假设半波的形状为半正弦波，则该游程内信号的均方根值与最大值之间的关系为

(9)

为验证这一假设，随机选取若干段游程，统计其均方根值和峰值的比例，如图4所示。结果表明，绝大多数游程的均方根值都位于0.5～0.7，式(9)偏于保守，且误差小于3 dB。

综合式(8)、式(9)，修正后的功率谱密度包络可以写成

(10)

图3 地脉动信号游程示意图Fig.3 Sketch of run-length of the ground microtremor signal

这一公式建立起功率谱密度包络与平均功率谱密度的关系，公式中只需要计算整段信号的均方根值及最大峰值，不需要实际进行游程的划分和稳态量/非稳态量的分解，对于大批量数据的处理显著节约了计算量，同时避免了参数和游程选取的任意性。

图5显示了时间长度T不同时根据式(10)计算得到的典型地脉动修正功率谱密度包络曲线，同时与未经修正的包络曲线进行了对比。结果显示，经过式(10)修正，功率谱密度包络的计算与时间长度T的选取基本无关。包络曲线与未修正的曲线形状一致，但幅值更高，有效避免了局部高量级振动被平均化的风险。

3 分析结果

应用上述方法，计算了绵阳花荄监测站2015年—2016年地脉动数据的功率谱密度，分析频率范围为0～50 Hz，分析频率间隔为0.048 8 Hz，对应的滑动平均窗长度为2 048个数据点，时间长度T选取为10 min。得到的功率谱密度包络如图6所示。在Intel Core i7-4770 CPU@3.4 GHz处理器上，使用MATLAB进行4线程并行计算，处理全部2年共计189亿个数据点所需的时间约3.5 h。

从图6中可见，地脉动信号的能量大部分集中在25～40 Hz的范围内，最大峰值约为3.1×10-12g2/Hz，其中横向(E向和N向)的振动幅度略大于垂直方向(Z向)。