吳佳靈 劉 溢 任大呈 李文峰

(北京航天計量測試技術(shù)研究所，北京 100076)

1 引 言

隨著現(xiàn)代工業(yè)的迅速發(fā)展、城市規(guī)模的日益擴大，振動對大都市生活環(huán)境和工作環(huán)境的影響引起了人們的普遍關(guān)注，國際上已把振動列為七大環(huán)境公害之一，并已開始著手研究振動污染規(guī)律、振動產(chǎn)生的原因、傳播路徑與控制方法等問題。引起環(huán)境振動有兩類：一類是人為的機械運動所引起的振動，一類是自然現(xiàn)象引起的振動，他們的危害程度隨著振動特性、振源布局和環(huán)境條件不同而異。

目前研究者們基于不同的用途提出了一些環(huán)境振動系統(tǒng)識別的方法，如：基于功率譜密度的峰值法[1]，基于離散時間數(shù)據(jù)的ARMA模型[2]，自然激勵技術(shù)(NExT)[3]，隨機子空間法[4]等。這些方法對數(shù)據(jù)進行減縮、方程進行求解、并進行矩陣運算，往往依賴于輸入模型或需要進行復雜的運算。本文首先構(gòu)建巴特沃斯(butterworth)帶通濾波器對原始數(shù)據(jù)進行濾波，提出了一種鄰近節(jié)點判定法來檢測環(huán)境振動，該方法使用鄰近節(jié)點構(gòu)建網(wǎng)絡進行分組計算，來減小計算量，再設(shè)定判據(jù)使用互相關(guān)方法判定環(huán)境振動。

2 數(shù)字濾波器設(shè)計

IIR 數(shù)字濾波器具有無限持續(xù)時間沖激響應，是遞歸型的線性時不變因果系統(tǒng)，其差分方程為

式中：y(n)——第n個時刻的輸出量；x(n)——第n個時刻的輸入量；ai——第i時刻輸入量系數(shù)；bi——第i時刻輸出量系數(shù)；M——輸入量階數(shù)；N——濾波器階數(shù)，進行Z變換后得到系統(tǒng)傳遞函數(shù)H(z)為

帶通濾波器可實現(xiàn)某一指定范圍的頻率成分可以順利通過，而不在該頻率范圍的頻率成分被抑制[5]。使用雙線性變換法設(shè)計出IIR的butterworth帶通濾波器，其中H(s)為模擬域傳遞函數(shù)的拉普拉斯變換，H(z)為數(shù)字域傳遞函數(shù)的Z變換，步驟如下：

(1)根據(jù)任務需求，明確性能指標，確定butterworth濾波器的輸入?yún)?shù)，包括邊界頻率和通帶最大衰減等；

(2)通過查表方法得到模擬butterworth濾波器的傳輸函數(shù)H(s)；

(3)通過如下公式模擬平面與數(shù)字平面建立對應關(guān)系

式中：Ω——模擬角頻率；T——取樣周期；ω——數(shù)字角頻率。

(4)利用雙線性變換法關(guān)系，將模擬butterworth濾波器H(s)轉(zhuǎn)換成數(shù)字butterworth濾波器H(z)；

3 鄰近節(jié)點判定

檢波器布置在不同的位置，振源信號以機械波的形式傳遞到檢波器的過程中，由于傳播介質(zhì)為各向異性，接收到的信號具有差異，如圖1所示，距離相近的節(jié)點接收到的信號由于傳播路徑相似，波形趨于一致。

圖1 鄰近節(jié)點接收信號示意圖Fig.1 The sketch how close nodes receive signal

本文所提出的鄰近節(jié)點判定法利用鄰近節(jié)點接收到的振動信號特征一致，而接收到的噪聲則表現(xiàn)出隨機性的特點，將檢波器陣列分組成為多個節(jié)點網(wǎng)絡，并使用滑動互相關(guān)方法，計算出每個網(wǎng)絡的平均互相關(guān)系數(shù)，通過設(shè)定判定閾值，最終判定環(huán)境振動，步驟如下：

(1)構(gòu)建鄰近節(jié)點網(wǎng)絡

節(jié)點總數(shù)為n，分別以每個節(jié)點為中心，建立節(jié)點網(wǎng)絡，第i個中心節(jié)點坐標為(xi,yi)，獲取節(jié)點i與其他節(jié)點間距離di。

獲取di中最小的k個節(jié)點min(di,k)，構(gòu)建成為節(jié)點i的鄰近節(jié)點波形數(shù)據(jù)集合Ci，即為節(jié)點i的鄰近節(jié)點網(wǎng)絡。該步驟使得每次只計算k個節(jié)點的數(shù)據(jù)，有效的減小了算法的空間復雜度。

(2)計算鄰近節(jié)點網(wǎng)絡的平均互相關(guān)系數(shù)Si

該步驟計算每個鄰近節(jié)點網(wǎng)絡的平均互相關(guān)系數(shù)，以中心節(jié)點i的波形數(shù)據(jù)作為基準，其鄰近節(jié)點j的波形數(shù)據(jù)對節(jié)點i的可移動量為[-Lδ，Lδ]，其中δ為采樣間隔，±L確定可以移動的范圍，以δ為時間間隔對節(jié)點j的波形進行移動，計算每個間隔中節(jié)點i與節(jié)點j的互相關(guān)系數(shù)，取得最大值作為兩個節(jié)點的互相關(guān)系數(shù)sij，最終得到每個網(wǎng)絡的平均互相關(guān)系數(shù)Si，如公式(6)、(7)。

式中：i=1,2,…,n；j=1,2,…,k；n——節(jié)點總數(shù)；k——臨近節(jié)點的總數(shù)；δ——采樣間隔；ui——第i個節(jié)點的數(shù)據(jù)；uj——第i個節(jié)點的數(shù)據(jù)，且uj∈Ci；L——平移的最大范圍；l——平移量；M——數(shù)據(jù)的窗口長度；sij(t)時刻t，節(jié)點i與第j個鄰近節(jié)點互相關(guān)系數(shù)最大值；Si——節(jié)點i與所有鄰近節(jié)點的互相關(guān)系數(shù)的均值。

(3)獲取環(huán)境振動數(shù)量

設(shè)定互相關(guān)閾值ρ∈(0，1)與數(shù)量閾值nρρ的個數(shù)Nρ，當Nρ>nρ時判定為環(huán)境振動。

(4)遍歷所有數(shù)據(jù)

使用(2M+1)δ為窗口長度，以(2M+1)δ為步長，從t=0時刻對全部數(shù)據(jù)進行掃描，為避免環(huán)境振動持續(xù)時間長引起重復計數(shù)，在p個步長時間(即p(2M+1)δ)內(nèi)多次出現(xiàn)的環(huán)境振動計為一個振動，記錄所有檢測到的數(shù)據(jù)。

4 實驗結(jié)果及分析

環(huán)境振動是一種寬頻帶的隨機振動[6]，由于高頻振動在地層中衰減很快，因此地表所感受到的主要是低頻振動信號。環(huán)境振動的振源一般分為自然振源和人工振源。包括變壓器、風機、壓縮機等動力設(shè)備產(chǎn)生的穩(wěn)態(tài)波振動；人員走動及城市軌道交通造成的非穩(wěn)態(tài)隨機振動；壓縮空氣氣源、油泵等輔助設(shè)備激勵引起的穩(wěn)態(tài)波或沖擊波振動等人工振源；以及大地脈動、浪涌和風力引起的非穩(wěn)態(tài)隨機振動等自然界振源。地震的振動頻率在0.1Hz～30Hz；風激振頻率范圍大多在0.1Hz～2Hz。人工振源(城市軌道交通和動力設(shè)備等)振幅變化較大，振動頻率在1Hz～150Hz；建筑物基礎(chǔ)自然頻率一般大于10Hz (與土壤特性相關(guān))；地面垂直向的自然頻率通常為6Hz～30Hz，水平向則更高些;實驗室人員走動產(chǎn)生的振動頻率一般在1Hz～3Hz范圍內(nèi)。環(huán)境振動的特點是不僅依賴于激勵的大小，還依賴于土壤、基礎(chǔ)、地板和建筑物其他結(jié)構(gòu)部件所組成的動力系統(tǒng)對振動的濾波效果。

使用公開發(fā)布的San Jacinto密集陣列的檢波器數(shù)據(jù)[7]來檢驗鄰近節(jié)點判定方法。美國南加州大學的Yehuda Ben-Zion課題組在美國南加州地區(qū)的San Jacinto斷層帶密集地布置了1108個檢波器，這些檢波器被布置在600m×600m的區(qū)域內(nèi)，共20行，每行不少于50個檢波器，間隔約10m，每列間隔約30m，采樣頻率為500個采樣/s，從2014年5月7日到2014年6月13日，記錄環(huán)境中的振動信號。

圖2 美國南加州San Jacinto密集陣列示意圖Fig.2 San Jacinto dense-array in southern California

本實驗下載了1108個檢波器從2014年5月11日0點開始10個小時的數(shù)據(jù)(共77GB)，使用butterworth帶通濾波器對原始數(shù)據(jù)進行濾波，由于環(huán)境中的振動信號一般為低頻信號，所以重構(gòu)采樣率為δ=50Hz，帶通濾波器的頻率范圍設(shè)為1Hz～10Hz，某一檢波器的數(shù)據(jù)濾波前后對比圖如圖3所示，可看出通過構(gòu)建帶通濾波器有效的控制了噪聲。

圖3 濾波前后波形對比圖Fig.3 Waveforms comparison between filtered and not filtered

設(shè)置鄰近節(jié)點波形相對于中心節(jié)點波形可移動的范圍為[-Lδ，Lδ]=[-0.5,0.5]s，步長與窗長(2M+1)δ=3s，10小時數(shù)據(jù)，p=5，即在15s內(nèi)產(chǎn)生的環(huán)境振動計為一次，總時長為ltotal=10h，則最多可檢測到環(huán)境振動個數(shù)為Nmax=ltotal/(p(2M+1)δ) =2 400個。為分析互相關(guān)閾值ρ與數(shù)量閾值nρ對監(jiān)測環(huán)境振動數(shù)量的影響，實施以下實驗：

(1)互相關(guān)閾值ρ對監(jiān)測環(huán)境振動數(shù)量的影響

固定數(shù)量閾值nρ=n/2，互相關(guān)閾值ρ∈(0，1)與檢測到環(huán)境振動數(shù)量的關(guān)系如圖4所示。在互相關(guān)系數(shù)在(0，0.35]范圍內(nèi)，所有的信號都將被檢測為環(huán)境振動；在(0.35，0.85)范圍內(nèi)，檢測到環(huán)境振動與互相關(guān)閾值呈現(xiàn)反相關(guān)性；在[0.85，1)范圍內(nèi)，沒有信號被檢測為環(huán)境振動?；ハ嚓P(guān)閾值范圍在(0.4，0.6)較為合理。

圖4 互相關(guān)閾值與檢測環(huán)境振動數(shù)量關(guān)系Fig.4 The relationship between cross-correlation threshold and the number of ambient vibration

(2)數(shù)量閾值nρ對監(jiān)測環(huán)境振動數(shù)量的影響

如圖5所示，固定互相關(guān)閾值ρ=0.5，數(shù)量閾值nρ與檢測到環(huán)境振動的個數(shù)呈現(xiàn)反相關(guān)。并且小于(0，150)范圍內(nèi)，曲線斜率最陡，(150，600)范圍內(nèi)逐漸平緩，在(600，1000)內(nèi)檢測到的環(huán)境振動較少，故數(shù)量閾值占節(jié)點總數(shù)的20%～60%較為合理。

圖5 臺站數(shù)閾值與檢測環(huán)境振動數(shù)量關(guān)系Fig.5 The relationship between station number threshold and the number of ambient vibration

環(huán)境振動信號以機械波的形式進行傳播，其頻率取決于振源的頻率，而波速僅與傳播介質(zhì)的性質(zhì)相關(guān)，在固體中同一振源產(chǎn)生的機械波以橫波與縱波的形式進行傳播，橫波波速按公式(8)計算，縱波波速按公式(9)計算。

式中：u——波速；G——介質(zhì)的剪切模量。

式中：Y——介質(zhì)的楊氏模量。

在液體和氣體中，振源產(chǎn)生的機械波只能以縱波的形式傳播，其波速按公式(10)計算。

式中：K——介質(zhì)的容變模量。

由于檢波器數(shù)量較多，對信號曲線進行抽樣顯示，每隔15個節(jié)點進行顯示，挑取檢測到的三個典型的環(huán)境振動，如圖6所示，其中橫坐標為采集數(shù)據(jù)的相對時間，縱坐標為檢波器的臺站序號。在同種固體介質(zhì)中，剪切模量總小于其楊氏模量，所以對于檢波器接收到固體中傳播的機械波信號，先接收到縱波，后收到橫波，在圖6的(a)、(b)中可以清晰的看到兩列以上的密集波形，其中最先到達的為速度最快的縱波，隨后橫波到達，不同臺站間縱波的到時時間差小于0.1s，在600m×600m的空間中，檢波器的最大距離為850m，故波速約為8km/s，符合縱波在地表的傳播速度。圖(a)縱波在9s左右到達，橫波在18s左右到達，縱波與橫波到達時間差為9s左右，圖(b)縱波在10s左右到達，橫波在13s左右到達，縱波與橫波到達時間差為3s左右，縱波與橫波時間差與振源到檢波器的距離相關(guān)，距離越大時間差越大，圖6(a)、(b)中，(a)振源的位置距離臺站更遠。在圖6(c)中，檢波器接收波形到達時間差為2.5s左右，檢波器最大距離為850m，故波速為340m/s，符合縱波在空氣中的傳播速度。

(a)

(b)

(c)圖6 檢測到的典型環(huán)境振動Fig.6 Typical ambient vibration detected

5 結(jié)束語

本文通過構(gòu)建butterworth帶通濾波器對檢波器采集的環(huán)境中的信號進行處理，建立了用于檢測環(huán)境振動的鄰近節(jié)點判定法，通過理論推導與實驗得出以下結(jié)論：

(a)鄰近節(jié)點判定法可以有效的檢測出在地面或空氣中傳播的機械波，并可判斷傳播介質(zhì)。

(b)鄰近節(jié)點判定法可以用于大型節(jié)點網(wǎng)絡，可降低算法的空間復雜度和時間復雜度。

(c)互相關(guān)閾值和臺站數(shù)閾值均與檢測環(huán)境振動的數(shù)量呈反比例關(guān)系，互相關(guān)閾值合理范圍為0.4～0.6，臺站數(shù)閾值的合理范圍為臺站總數(shù)的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] 華師韓，王青. 數(shù)字濾波在動態(tài)測量數(shù)據(jù)處理的應用[J]. 宇航計測技術(shù), 2005 , 25 (2) :50～54.

[6] 朱石堅,樓京俊,何其偉,等. 振動理論與隔振技術(shù)[M ]. 北京:國防工業(yè)出版社, 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.