基于最小二乘離散狀態空間的振級計算方法?

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

（1.南京航空航天大學機械結構力學及控制國家重點實驗室 南京，210016）

（2.江蘇環保產業技術研究院股份公司 南京，210036）

引 言

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

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

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

1 理論背景

1.1 計權網絡的狀態空間形式

文獻［7］給出了所關心頻帶內頻率計權傳遞函數的構造參數。對于Z 振級，采用Wk計權因子組成的計權網絡，其由4 種功能傳遞函數組成，分別為

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

將式（1）~（4）所代表的濾波器級聯，得到計權濾波器為

式（5）按系數展開為

聯立式（1）~（5），求得Wk計權傳遞函數的分子與分母系數如表1 所示。

表1 Wk計權傳遞函數的分子與分母系數Tab.1 The numerator and denominator coefficients of the transfer function of the weightings Wk

將這些系數轉化為連續SISO 系統的狀態空間，對應的系統矩陣分別為

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

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

表2 Wd計權傳遞函數的分子與分母系數Tab.2 The numerator and denominator coefficients of the transfer function of the weightings Wd

上述得到的連續系統不能直接將離散時域數據作為輸入，在實際應用中還需將其轉換為離散系統。

為實現該轉換，采用最小二乘向量擬合法，使連續系統和離散系統的頻響函數在關心頻帶內誤差最小。設連續系統頻響與離散系統頻響之間的誤差函數為

進一步得到所有譜線的誤差函數為

其中：Nf為譜線數；上標H 表示共軛轉置。

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

其中：J為Jacobian 矩陣。

P，L分別為

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

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

經過簡單化簡，可得到分子、分母多項式系數，進而得到離散系統的狀態空間矩陣Ad，Bd，Cd和Dd，其中Ad為分母多項式系數向量α的伴隨矩陣。

1.2 振級估計

得到離散狀態空間后，計權加速度可表示為SISO 系統的輸出，即

其中：y（k）為系統的狀態響應；x(k)為原始加速度；xw(k)為計權后的加速度。

基于計權后的加速度可采用線性平均或指數平均估計加速度的計權運行均方根值，即

其中：N為計算時間周期T內所包含的數據樣本數，通常T取1 s。

計算得到瞬時振級值為

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

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

2 算例分析

2.1 離散計權系統

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

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

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

2.2 仿真算例

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

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

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

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

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

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

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

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

2.3 實測算例

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

圖7 測試現場Fig.7 Test scene

設積分周期為1 s，步進時間間隔為0.1 s，基于圖8 所示的實測加速度信號，分別采用濾波法和本研究方法計權計算瞬時Z 振級，并取每趟列車經過時的最大值連成曲線，兩種方法估計出的最大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 結果顯示，本研究方法的計算速度約為濾波法的7 倍。實測算例進一步證明，所提出的計權計算Z 振級的方法，在兼顧了計算精度的同時，明顯提高了計算效率。

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

3 結 論

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

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

3） 所提出的計權方法由時域內離散序列的輸入-輸出模型實現，能夠實時進行加速度數據的計權計算。