基于最小二乘離散狀態(tài)空間的振級(jí)計(jì)算方法?

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

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

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

引 言

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

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

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

1 理論背景

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

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

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

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

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

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

表1 Wk計(jì)權(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)空間，對(duì)應(yīng)的系統(tǒng)矩陣分別為

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

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

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

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

為實(shí)現(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）的等式兩端同時(shí)左乘Jacobian 矩陣的共軛轉(zhuǎn)置，得到正則方程為

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

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

1.2 振級(jí)估計(jì)

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

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

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

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

計(jì)算得到瞬時(shí)振級(jí)值為

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

取瞬時(shí)振級(jí)中的最大值即為最大振級(jí)。

2 算例分析

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

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

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

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

2.2 仿真算例

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

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

圖4 仿真振動(dòng)加速度信號(hào)Fig.4 The simulated acceleration signal

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

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

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

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

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

2.3 實(shí)測(cè)算例

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

圖7 測(cè)試現(xiàn)場(chǎng)Fig.7 Test scene

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

圖8 實(shí)測(cè)加速度信號(hào)Fig.8 Measured acceleration signal

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

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

圖10 兩種方法的計(jì)算時(shí)間對(duì)比Fig.10 Calculation time of two methods

3 結(jié) 論

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

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

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