夏遵平， 王 彤， 王 彧， 范新亮

（1.南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點實驗室 南京，210016）

（2.江蘇環(huán)保產(chǎn)業(yè)技術(shù)研究院股份公司 南京，210036）

引 言

近年來，由交通干擾、機(jī)械運(yùn)轉(zhuǎn)等人為因素產(chǎn)生的環(huán)境振動問題尤為突出。以交通干擾為例，軌道交通（高鐵、地鐵、輕軌等）為人們的出行提供了便捷，但列車運(yùn)行產(chǎn)生的環(huán)境振動將嚴(yán)重干擾軌道線路周圍居民的生活，甚至損毀周邊敏感建筑［1-2］。因此，環(huán)保部門將環(huán)境振動問題作為環(huán)境綜合治理的重要內(nèi)容，并為此新修訂了環(huán)境振動的測試與評價準(zhǔn)則［3-4］，規(guī)定采用人體全身振動評價指標(biāo)中的計權(quán)振動加速度級（簡稱振級）作為環(huán)境振動的評價指標(biāo)［5-6］。

我國新版人體全身振動評價標(biāo)準(zhǔn)［7］直接引用國際同名標(biāo)準(zhǔn)ISO 2631-1—1997［8］，該標(biāo)準(zhǔn)詳細(xì)規(guī)定了加速度計權(quán)的方式和振級的計算方法，要求加速度需按1/3 倍頻程進(jìn)行計權(quán)［9］。因此，對測試得到的離散加速度數(shù)據(jù)按1/3 倍頻程的方式進(jìn)行計權(quán)是振級計算過程中的重要步驟［10］。常見的加速度計權(quán)方法有頻譜方法［11］、百分比帶寬濾波法［12-13］等。頻譜方法首先估計加速度信號功率譜密度，然后根據(jù)Parseval 定理計算每個1/3 倍頻程帶內(nèi)的總功率并乘以對應(yīng)的計權(quán)因子。該方法的優(yōu)點是利用了快速Fourier 變換，計算效率很高，但對于要求以1 s 為積分周期的振級估計，其頻率分辨為1 Hz，故無法基于該Fourier 譜對頻程帶小于1 Hz 的1/3 倍頻程進(jìn)行劃分，只能采用近似估計，因而可能引入較大的誤差［14］。等百分比帶寬濾波法首先根據(jù)1/3 倍頻程帶進(jìn)行數(shù)字濾波，得到每個頻程帶內(nèi)的時域數(shù)據(jù)，然后求其均方根并乘以對應(yīng)的計權(quán)因子。該方法避免了低頻帶寬無法劃分的問題，但需要采用濾波器組對每一頻程帶進(jìn)行時域濾波，計算效率較低。此外，文獻(xiàn)［7］給出了一種用于擬合計權(quán)網(wǎng)絡(luò)曲線的模擬濾波器，采用該濾波器計權(quán)不再直接基于1/3 倍頻程的劃分，但需要采用合適的方法將其轉(zhuǎn)化為數(shù)字濾波器，如文獻(xiàn)［10］中采用的雙線性變換法。然而，欲使變換后的系統(tǒng)在所關(guān)心頻帶內(nèi)具有可靠的精度，該方法所需的采樣頻率要遠(yuǎn)大于采樣定理的要求。在采樣時間固定時，較高的采樣頻率意味著更大的數(shù)據(jù)樣本數(shù)，將耗費更多的數(shù)據(jù)存儲和計算開銷。

筆者提出一種基于最小二乘-離散狀態(tài)空間的加速度計權(quán)方法。首先，由文獻(xiàn)［7］給定的計權(quán)參數(shù)構(gòu)建出單輸入-單輸出（single input single output， 簡稱SISO）系統(tǒng)的傳遞函數(shù)；然后，采用最小二乘擬合法，將該傳遞函數(shù)轉(zhuǎn)化為離散的SISO 狀態(tài)空間模型；最后，以原始加速度數(shù)據(jù)作為系統(tǒng)輸入，得到的系統(tǒng)響應(yīng)即為計權(quán)加速度。本方法避免了1/3倍頻程的劃分和對應(yīng)頻帶數(shù)據(jù)的轉(zhuǎn)換，有效提高了計算效率。通過仿真和實測算例驗證了所提方法具有兼顧計算精度和計算效率的優(yōu)點。

1 理論背景

1.1 計權(quán)網(wǎng)絡(luò)的狀態(tài)空間形式

文獻(xiàn)［7］給出了所關(guān)心頻帶內(nèi)頻率計權(quán)傳遞函數(shù)的構(gòu)造參數(shù)。對于Z 振級，采用Wk計權(quán)因子組成的計權(quán)網(wǎng)絡(luò)，其由4 種功能傳遞函數(shù)組成，分別為

其中：Hh(s)，Hl(s)，Ht(s)和Hs(s)分別為高通、低通、加速度-速度轉(zhuǎn)換和高階濾波器的傳遞函數(shù)；s為拉氏變量。

將式（1）~（4）所代表的濾波器級聯(lián)，得到計權(quán)濾波器為

式（5）按系數(shù)展開為

聯(lián)立式（1）~（5），求得Wk計權(quán)傳遞函數(shù)的分子與分母系數(shù)如表1 所示。

表1 Wk計權(quán)傳遞函數(shù)的分子與分母系數(shù)Tab.1 The numerator and denominator coefficients of the transfer function of the weightings Wk

將這些系數(shù)轉(zhuǎn)化為連續(xù)SISO 系統(tǒng)的狀態(tài)空間，對應(yīng)的系統(tǒng)矩陣分別為

其中：A∈R8×8為狀態(tài)矩陣；B∈R8×1為輸入矩陣；C∈R1×8為輸出矩陣；D∈R1×1為前饋矩陣。

同理，根據(jù)Wd計權(quán)值亦可求得用于橫向（y）和縱向（x）振級計算的狀態(tài)空間矩陣，其計權(quán)傳遞函數(shù)的分子與分母系數(shù)如表2 所示。

表2 Wd計權(quán)傳遞函數(shù)的分子與分母系數(shù)Tab.2 The numerator and denominator coefficients of the transfer function of the weightings Wd

上述得到的連續(xù)系統(tǒng)不能直接將離散時域數(shù)據(jù)作為輸入，在實際應(yīng)用中還需將其轉(zhuǎn)換為離散系統(tǒng)。

為實現(xiàn)該轉(zhuǎn)換，采用最小二乘向量擬合法，使連續(xù)系統(tǒng)和離散系統(tǒng)的頻響函數(shù)在關(guān)心頻帶內(nèi)誤差最小。設(shè)連續(xù)系統(tǒng)頻響與離散系統(tǒng)頻響之間的誤差函數(shù)為

進(jìn)一步得到所有譜線的誤差函數(shù)為

其中：Nf為譜線數(shù)；上標(biāo)H 表示共軛轉(zhuǎn)置。

要使式（11）取極小值，可由最小二乘得

其中：J為Jacobian 矩陣。

P，L分別為

將式（13）的等式兩端同時左乘Jacobian 矩陣的共軛轉(zhuǎn)置，得到正則方程為

其中：R=Re(PHP)；S=Re(PHL)；T=Re(LHL)。

經(jīng)過簡單化簡，可得到分子、分母多項式系數(shù)，進(jìn)而得到離散系統(tǒng)的狀態(tài)空間矩陣Ad，Bd，Cd和Dd，其中Ad為分母多項式系數(shù)向量α的伴隨矩陣。

1.2 振級估計

得到離散狀態(tài)空間后，計權(quán)加速度可表示為SISO 系統(tǒng)的輸出，即

其中：y（k）為系統(tǒng)的狀態(tài)響應(yīng)；x(k)為原始加速度；xw(k)為計權(quán)后的加速度。

基于計權(quán)后的加速度可采用線性平均或指數(shù)平均估計加速度的計權(quán)運(yùn)行均方根值，即

其中：N為計算時間周期T內(nèi)所包含的數(shù)據(jù)樣本數(shù)，通常T取1 s。

計算得到瞬時振級值為

其中：a0=10-6m/s2，為基準(zhǔn)加速度有效值。

取瞬時振級中的最大值即為最大振級。

2 算例分析

2.1 離散計權(quán)系統(tǒng)

當(dāng)采樣頻率為400 Hz 時，分別采用雙線性變換法［6］和本研究的最小二乘狀態(tài)空間法計算離散系統(tǒng)，得到1~100 Hz 內(nèi)Wk和Wd計權(quán)因子組成的頻響函數(shù)，不同方法得到的Wk和Wd計權(quán)因子曲線分別如圖1，2 所示。圖中理論值取自于文獻(xiàn)［7］中的1/3 倍頻程頻率對應(yīng)的計權(quán)因子。結(jié)果顯示，在400 Hz（等于4fmax）的采樣率下，通過雙線性變換方法獲得的離散系統(tǒng)的頻響曲線在頻率大于31.5 Hz后呈現(xiàn)明顯的偏差，而采用本研究方法獲得的頻響曲線在整個頻帶內(nèi)都與理論值匹配較好。

圖1 不同方法得到的Wk計權(quán)因子曲線Fig.1 Curves of Wk obtained by different methods

圖2 不同方法得到的Wd計權(quán)因子曲線Fig.2 Curves of Wd obtained by different methods

2.2 仿真算例

在30~80 Hz 頻帶內(nèi)按2 Hz 等間隔取頻率，生成幅值為m/s2、初始相位隨機(jī)的26 組單頻加速度正弦信號。各頻率成分的均方根為1 m/s2，計權(quán)均方根即為計權(quán)因子數(shù)值，如圖3 所示。將上述26 組信號加法混合，并按400 Hz 采樣率生成仿真算例，以驗證本方法的計算精度和計算效率，仿真振動加速度信號如圖4 所示。

圖3 各頻率成分的計權(quán)均方根Fig.3 The weighted-rms of each frequency component

圖4 仿真振動加速度信號Fig.4 The simulated acceleration signal

以Z 振級（Wk計權(quán)）的計算為例，采用不同方法估計出的瞬時Z 振級曲線如圖5 所示。其中，理論值由26 組頻率分布在1/3 倍頻程帶內(nèi)的數(shù)據(jù)乘以對應(yīng)的計權(quán)因子求得，濾波法計算值采用等百分比帶寬濾波器分別濾出各1/3 倍頻程帶內(nèi)的時域數(shù)據(jù)并乘以對應(yīng)的計權(quán)因子求得。

圖5 不同方法估計出的瞬時Z 振級曲線Fig.5 Instantaneous VLz values computed by different methods

由圖5 可知，本研究方法的計算精度與濾波法相近，均明顯高于雙線性變換法，這是由于采用雙線性變換法進(jìn)行連續(xù)-離散轉(zhuǎn)化時產(chǎn)生誤差造成的。為檢驗各方法的計算效率，對5 組分別為200，400，800，1 600 和3 200 s 采樣長度的數(shù)據(jù)進(jìn)行振級計算，不同數(shù)據(jù)長度下各方法的計算效率對比如圖6 所示。

圖6 不同數(shù)據(jù)長度下各方法的計算效率對比Fig.6 Calculation time of the methods under different data length

由圖6 可知，本研究方法具有與雙線性變換方法相近的計算效率，即處理長樣本數(shù)據(jù)的速度明顯高于傳統(tǒng)的濾波法。

2.3 實測算例

根據(jù)《環(huán)境影響評價技術(shù)導(dǎo)則-城市軌道交通》［3］中的規(guī)定，實測南京某地鐵隧洞壁的振動源強(qiáng)，測試現(xiàn)場如圖7 所示。采用無人值守采集設(shè)備采集數(shù)據(jù)，連續(xù)測試17 趟列車經(jīng)過時的垂向加速度信號。

圖7 測試現(xiàn)場Fig.7 Test scene

設(shè)積分周期為1 s，步進(jìn)時間間隔為0.1 s，基于圖8 所示的實測加速度信號，分別采用濾波法和本研究方法計權(quán)計算瞬時Z 振級，并取每趟列車經(jīng)過時的最大值連成曲線，兩種方法估計出的最大Z 振級和計算時間對比分別如圖9，10 所示。

圖8 實測加速度信號Fig.8 Measured acceleration signal

圖9 兩種方法估計出的最大Z 振級Fig.9 Maximum Z-Vibration level computed by two methods

圖9 中2 條最大Z 振級曲線近乎重合，表明兩種方法計算的精度相近。圖10 結(jié)果顯示，本研究方法的計算速度約為濾波法的7 倍。實測算例進(jìn)一步證明，所提出的計權(quán)計算Z 振級的方法，在兼顧了計算精度的同時，明顯提高了計算效率。

圖10 兩種方法的計算時間對比Fig.10 Calculation time of two methods

3 結(jié) 論

1） 通過單輸入-單輸出離散狀態(tài)空間模型完成加速度數(shù)據(jù)的計權(quán)，過程中避免了1/3 倍頻程的轉(zhuǎn)換，有效提高了計算效率，并為實時加速度數(shù)據(jù)計權(quán)提供了算法依據(jù)。

2） 本研究方法使用最小二乘復(fù)頻域擬合的方法獲得離散狀態(tài)空間系統(tǒng)的參數(shù)，在不增加采樣頻率的前提下，使離散系統(tǒng)和連續(xù)系統(tǒng)的頻響幅值在關(guān)心頻帶內(nèi)的誤差最小，提高了計算精度。

3） 所提出的計權(quán)方法由時域內(nèi)離散序列的輸入-輸出模型實現(xiàn)，能夠?qū)崟r進(jìn)行加速度數(shù)據(jù)的計權(quán)計算。