基于互相關(guān)算法的環(huán)境振動(dòng)檢測(cè)方法研究

吳佳靈 劉 溢 任大呈 李文峰

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

3 鄰近節(jié)點(diǎn)判定

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實(shí)驗(yàn)結(jié)果及分析

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

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

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

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

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

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

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

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

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

(2)數(shù)量閾值nρ對(duì)監(jiān)測(cè)環(huán)境振動(dòng)數(shù)量的影響

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

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

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

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

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

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

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

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

(a)

(b)

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

5 結(jié)束語(yǔ)

本文通過(guò)構(gòu)建butterworth帶通濾波器對(duì)檢波器采集的環(huán)境中的信號(hào)進(jìn)行處理，建立了用于檢測(cè)環(huán)境振動(dòng)的鄰近節(jié)點(diǎn)判定法，通過(guò)理論推導(dǎo)與實(shí)驗(yàn)得出以下結(jié)論：

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

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

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

[6] 朱石堅(jiān),樓京俊,何其偉,等. 振動(dòng)理論與隔振技術(shù)[M ]. 北京:國(guó)防工業(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.