基于互相关算法的环境振动检测方法研究
2018-05-05吴佳灵任大呈李文峰
吴佳灵 刘 溢 任大呈 李文峰
(北京航天计量测试技术研究所，北京 100076)
1 引 言
随着现代工业的迅速发展、城市规模的日益扩大，振动对大都市生活环境和工作环境的影响引起了人们的普遍关注，国际上已把振动列为七大环境公害之一，并已开始着手研究振动污染规律、振动产生的原因、传播路径与控制方法等问题。引起环境振动有两类：一类是人为的机械运动所引起的振动，一类是自然现象引起的振动，他们的危害程度随着振动特性、振源布局和环境条件不同而异。
目前研究者们基于不同的用途提出了一些环境振动系统识别的方法，如：基于功率谱密度的峰值法[1]，基于离散时间数据的ARMA模型[2]，自然激励技术(NExT)[3]，随机子空间法[4]等。这些方法对数据进行减缩、方程进行求解、并进行矩阵运算，往往依赖于输入模型或需要进行复杂的运算。本文首先构建巴特沃斯(butterworth)带通滤波器对原始数据进行滤波，提出了一种邻近节点判定法来检测环境振动，该方法使用邻近节点构建网络进行分组计算，来减小计算量，再设定判据使用互相关方法判定环境振动。
2 数字滤波器设计
IIR 数字滤波器具有无限持续时间冲激响应，是递归型的线性时不变因果系统，其差分方程为
式中：y(n)——第n个时刻的输出量；x(n)——第n个时刻的输入量；ai——第i时刻输入量系数；bi——第i时刻输出量系数；M——输入量阶数；N——滤波器阶数，进行Z变换后得到系统传递函数H(z)为
带通滤波器可实现某一指定范围的频率成分可以顺利通过，而不在该频率范围的频率成分被抑制[5]。使用双线性变换法设计出IIR的butterworth带通滤波器，其中H(s)为模拟域传递函数的拉普拉斯变换，H(z)为数字域传递函数的Z变换，步骤如下：
(1)根据任务需求，明确性能指标，确定butterworth滤波器的输入参数，包括边界频率和通带最大衰减等；
(2)通过查表方法得到模拟butterworth滤波器的传输函数H(s)；
(3)通过如下公式模拟平面与数字平面建立对应关系
式中：Ω——模拟角频率；T——取样周期；ω——数字角频率。
(4)利用双线性变换法关系，将模拟butterworth滤波器H(s)转换成数字butterworth滤波器H(z)；
3 邻近节点判定
检波器布置在不同的位置，振源信号以机械波的形式传递到检波器的过程中，由于传播介质为各向异性，接收到的信号具有差异，如图1所示，距离相近的节点接收到的信号由于传播路径相似，波形趋于一致。
图1 邻近节点接收信号示意图Fig.1 The sketch how close nodes receive signal
本文所提出的邻近节点判定法利用邻近节点接收到的振动信号特征一致，而接收到的噪声则表现出随机性的特点，将检波器阵列分组成为多个节点网络，并使用滑动互相关方法，计算出每个网络的平均互相关系数，通过设定判定阈值，最终判定环境振动，步骤如下：
(1)构建邻近节点网络
节点总数为n，分别以每个节点为中心，建立节点网络，第i个中心节点坐标为(xi,yi)，获取节点i与其他节点间距离di。
获取di中最小的k个节点min(di,k)，构建成为节点i的邻近节点波形数据集合Ci，即为节点i的邻近节点网络。该步骤使得每次只计算k个节点的数据，有效的减小了算法的空间复杂度。
(2)计算邻近节点网络的平均互相关系数Si
该步骤计算每个邻近节点网络的平均互相关系数，以中心节点i的波形数据作为基准，其邻近节点j的波形数据对节点i的可移动量为[-Lδ，Lδ]，其中δ为采样间隔，±L确定可以移动的范围，以δ为时间间隔对节点j的波形进行移动，计算每个间隔中节点i与节点j的互相关系数，取得最大值作为两个节点的互相关系数sij，最终得到每个网络的平均互相关系数Si，如公式(6)、(7)。
式中：i=1,2,…,n；j=1,2,…,k；n——节点总数；k——临近节点的总数；δ——采样间隔；ui——第i个节点的数据；uj——第i个节点的数据，且uj∈Ci；L——平移的最大范围；l——平移量；M——数据的窗口长度；sij(t)时刻t，节点i与第j个邻近节点互相关系数最大值；Si——节点i与所有邻近节点的互相关系数的均值。
(3)获取环境振动数量
设定互相关阈值ρ∈(0，1)与数量阈值nρ
(4)遍历所有数据
使用(2M+1)δ为窗口长度，以(2M+1)δ为步长，从t=0时刻对全部数据进行扫描，为避免环境振动持续时间长引起重复计数，在p个步长时间(即p(2M+1)δ)内多次出现的环境振动计为一个振动，记录所有检测到的数据。
4 实验结果及分析
环境振动是一种宽频带的随机振动[6]，由于高频振动在地层中衰减很快，因此地表所感受到的主要是低频振动信号。环境振动的振源一般分为自然振源和人工振源。包括变压器、风机、压缩机等动力设备产生的稳态波振动；人员走动及城市轨道交通造成的非稳态随机振动；压缩空气气源、油泵等辅助设备激励引起的稳态波或冲击波振动等人工振源；以及大地脉动、浪涌和风力引起的非稳态随机振动等自然界振源。地震的振动频率在0.1Hz～30Hz；风激振频率范围大多在0.1Hz～2Hz。人工振源(城市轨道交通和动力设备等)振幅变化较大，振动频率在1Hz～150Hz；建筑物基础自然频率一般大于10Hz (与土壤特性相关)；地面垂直向的自然频率通常为6Hz～30Hz，水平向则更高些;实验室人员走动产生的振动频率一般在1Hz～3Hz范围内。环境振动的特点是不仅依赖于激励的大小，还依赖于土壤、基础、地板和建筑物其他结构部件所组成的动力系统对振动的滤波效果。
使用公开发布的San Jacinto密集阵列的检波器数据[7]来检验邻近节点判定方法。美国南加州大学的Yehuda Ben-Zion课题组在美国南加州地区的San Jacinto断层带密集地布置了1108个检波器，这些检波器被布置在600m×600m的区域内，共20行，每行不少于50个检波器，间隔约10m，每列间隔约30m，采样频率为500个采样/s，从2014年5月7日到2014年6月13日，记录环境中的振动信号。
图2 美国南加州San Jacinto密集阵列示意图Fig.2 San Jacinto dense-array in southern California
本实验下载了1108个检波器从2014年5月11日0点开始10个小时的数据(共77GB)，使用butterworth带通滤波器对原始数据进行滤波，由于环境中的振动信号一般为低频信号，所以重构采样率为δ=50Hz，带通滤波器的频率范围设为1Hz～10Hz，某一检波器的数据滤波前后对比图如图3所示，可看出通过构建带通滤波器有效的控制了噪声。
图3 滤波前后波形对比图Fig.3 Waveforms comparison between filtered and not filtered
设置邻近节点波形相对于中心节点波形可移动的范围为[-Lδ，Lδ]=[-0.5,0.5]s，步长与窗长(2M+1)δ=3s，10小时数据，p=5，即在15s内产生的环境振动计为一次，总时长为ltotal=10h，则最多可检测到环境振动个数为Nmax=ltotal/(p(2M+1)δ) =2 400个。为分析互相关阈值ρ与数量阈值nρ对监测环境振动数量的影响，实施以下实验：
(1)互相关阈值ρ对监测环境振动数量的影响
固定数量阈值nρ=n/2，互相关阈值ρ∈(0，1)与检测到环境振动数量的关系如图4所示。在互相关系数在(0，0.35]范围内，所有的信号都将被检测为环境振动；在(0.35，0.85)范围内，检测到环境振动与互相关阈值呈现反相关性；在[0.85，1)范围内，没有信号被检测为环境振动。互相关阈值范围在(0.4，0.6)较为合理。
图4 互相关阈值与检测环境振动数量关系Fig.4 The relationship between cross-correlation threshold and the number of ambient vibration
(2)数量阈值nρ对监测环境振动数量的影响
如图5所示，固定互相关阈值ρ=0.5，数量阈值nρ与检测到环境振动的个数呈现反相关。并且小于(0，150)范围内，曲线斜率最陡，(150，600)范围内逐渐平缓，在(600，1000)内检测到的环境振动较少，故数量阈值占节点总数的20%～60%较为合理。
图5 台站数阈值与检测环境振动数量关系Fig.5 The relationship between station number threshold and the number of ambient vibration
环境振动信号以机械波的形式进行传播，其频率取决于振源的频率，而波速仅与传播介质的性质相关，在固体中同一振源产生的机械波以横波与纵波的形式进行传播，横波波速按公式(8)计算，纵波波速按公式(9)计算。
式中：u——波速；G——介质的剪切模量。
式中：Y——介质的杨氏模量。
在液体和气体中，振源产生的机械波只能以纵波的形式传播，其波速按公式(10)计算。
式中：K——介质的容变模量。
由于检波器数量较多，对信号曲线进行抽样显示，每隔15个节点进行显示，挑取检测到的三个典型的环境振动，如图6所示，其中横坐标为采集数据的相对时间，纵坐标为检波器的台站序号。在同种固体介质中，剪切模量总小于其杨氏模量，所以对于检波器接收到固体中传播的机械波信号，先接收到纵波，后收到横波，在图6的(a)、(b)中可以清晰的看到两列以上的密集波形，其中最先到达的为速度最快的纵波，随后横波到达，不同台站间纵波的到时时间差小于0.1s，在600m×600m的空间中，检波器的最大距离为850m，故波速约为8km/s，符合纵波在地表的传播速度。图(a)纵波在9s左右到达，横波在18s左右到达，纵波与横波到达时间差为9s左右，图(b)纵波在10s左右到达，横波在13s左右到达，纵波与横波到达时间差为3s左右，纵波与横波时间差与振源到检波器的距离相关，距离越大时间差越大，图6(a)、(b)中，(a)振源的位置距离台站更远。在图6(c)中，检波器接收波形到达时间差为2.5s左右，检波器最大距离为850m，故波速为340m/s，符合纵波在空气中的传播速度。
(a)
(b)
(c)图6 检测到的典型环境振动Fig.6 Typical ambient vibration detected
5 结束语
本文通过构建butterworth带通滤波器对检波器采集的环境中的信号进行处理，建立了用于检测环境振动的邻近节点判定法，通过理论推导与实验得出以下结论：
(a)邻近节点判定法可以有效的检测出在地面或空气中传播的机械波，并可判断传播介质。
(b)邻近节点判定法可以用于大型节点网络，可降低算法的空间复杂度和时间复杂度。
(c)互相关阈值和台站数阈值均与检测环境振动的数量呈反比例关系，互相关阈值合理范围为0.4～0.6，台站数阈值的合理范围为台站总数的40%～60%。
[1] Bendat J S ,Piersol A G. Engineering applications of correlation and spectral analysis[M] . 2nd edition. New York : John Wiley &Sons , 1993.
[2] Andersen P , Brincker R , Kirkegaard P H. Theory of covariance equivalent ARMAV models of civil engineering structures[A] . Proceedings of IMAC14 , the 14th international modal analysis conference[C] . 1996. 518～524.
[3] James III G H , Carne T G, Lauffer J P. The natural excitation technique (NExT) for modal parameter extraction fromoperatingstructures[J] . International Journal of Analytical and Experimental Modal Analysis , 1995 , 10 (4) : 260～277.
[4] Van OverscheeP , De Moor B. Subspace identification for linear systems : theory , implementation and applications[M] . Dordrecht :Kluwer Academic Publishers , 1996.
[5] 华师韩，王青. 数字滤波在动态测量数据处理的应用[J]. 宇航计测技术, 2005 , 25 (2) :50～54.
[6] 朱石坚,楼京俊,何其伟,等. 振动理论与隔振技术[M ]. 北京:国防工业出版社, 2006.
[7] Ben-Zion, Yehuda, et al. "Basic data features and results from a spatially dense seismic array onthe San Jacinto fault zone." Geophysical Journal International 202.1 (2015): 370～380.
