本發(fā)明涉及高速列車(chē)懸掛系統(tǒng)故障辨識(shí)領(lǐng)域。更具體地，涉及一種基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法。

背景技術(shù)：

鐵路是國(guó)民經(jīng)濟(jì)發(fā)展的大動(dòng)脈，是實(shí)現(xiàn)國(guó)家現(xiàn)代化建設(shè)的重要支撐。高速鐵路作為一種安全、高速、舒適的交通運(yùn)輸方式，對(duì)國(guó)家經(jīng)濟(jì)社會(huì)的發(fā)展起到巨大的推動(dòng)作用。2016年，調(diào)整后的《中長(zhǎng)期鐵路網(wǎng)規(guī)劃》明確提出建設(shè)以“八縱八橫”主通道為骨架、區(qū)域連接線(xiàn)銜接、城際鐵路補(bǔ)充的高速鐵路網(wǎng)。高速鐵路運(yùn)輸在我國(guó)交通運(yùn)輸尤其是鐵路運(yùn)輸中的地位變得越來(lái)越重要。安全性是高速鐵路的核心競(jìng)爭(zhēng)力，保障高速列車(chē)的安全運(yùn)行不僅關(guān)系到乘客的生命安全，更會(huì)對(duì)我國(guó)高速鐵路的發(fā)展和走出國(guó)門(mén)產(chǎn)生深遠(yuǎn)的影響。此外，隨著我國(guó)高速鐵路的快速發(fā)展，大量高速列車(chē)投入運(yùn)營(yíng)，由此產(chǎn)生的一系列維修保養(yǎng)問(wèn)題亟待解決。

作為保障高速列車(chē)安全運(yùn)行的關(guān)鍵設(shè)備，抗蛇行減震器能夠有效抑制車(chē)輛的蛇行運(yùn)動(dòng)，對(duì)車(chē)輛運(yùn)行的橫向穩(wěn)定性非常重要?，F(xiàn)有的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法僅通過(guò)間接分析車(chē)輛振動(dòng)數(shù)據(jù)判斷零部件的故障狀態(tài)，容易受到采集設(shè)備、數(shù)據(jù)處理方法和耦合故障的影響，辨識(shí)效率較低。

因此，需要提供一種直接對(duì)抗蛇行減震器的性能參數(shù)進(jìn)行辨識(shí)并基于性能參數(shù)辨識(shí)結(jié)果對(duì)抗蛇行減震器的故障進(jìn)行辨識(shí)的基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于提供一種基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法，解決現(xiàn)有的高速列車(chē)懸掛系統(tǒng)性能參數(shù)及故障辨識(shí)普遍采用半車(chē)模型進(jìn)行算法研究，模型簡(jiǎn)化程度較高，與實(shí)際情況差別較大的問(wèn)題，縮小模型與實(shí)際運(yùn)營(yíng)的高速列車(chē)間的差異，提高辨識(shí)結(jié)果的準(zhǔn)確性。

為達(dá)到上述目的，本發(fā)明采用下述技術(shù)方案：

一種基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法，包括如下步驟：

s1、建立車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型；

s2、根據(jù)車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型，建立懸掛系統(tǒng)的離散狀態(tài)方程和離散觀測(cè)方程；

s3、根據(jù)懸掛系統(tǒng)的離散狀態(tài)方程和離散觀測(cè)方程，基于邊緣化粒子濾波算法對(duì)懸掛系統(tǒng)的抗蛇行減震器進(jìn)行性能參數(shù)辨識(shí)和故障辨識(shí)。

優(yōu)選地，步驟s2中，

懸掛系統(tǒng)的離散狀態(tài)方程為：

xk+1＝axk+buk+ewk

其中，xk為第k個(gè)采樣點(diǎn)的狀態(tài)變量,uk為第k個(gè)采樣點(diǎn)的系統(tǒng)輸入，wk為第k個(gè)采樣點(diǎn)的過(guò)程噪聲，離散狀態(tài)方程中懸掛系統(tǒng)矩陣離散狀態(tài)方程中過(guò)程噪聲矩陣和為連續(xù)狀態(tài)方程中懸掛系統(tǒng)矩陣，為連續(xù)狀態(tài)方程中過(guò)程噪聲矩陣；

懸掛系統(tǒng)的離散觀測(cè)方程為：

yk＝cxk+duk+fvk

其中，離散狀態(tài)方程中懸掛系統(tǒng)矩陣和為連續(xù)狀態(tài)方程中懸掛系統(tǒng)矩陣，離散狀態(tài)方程中觀測(cè)噪聲矩陣為連續(xù)狀態(tài)方程中觀測(cè)噪聲矩陣，vk為第k個(gè)采樣點(diǎn)的觀測(cè)噪聲。

優(yōu)選地，步驟s3的具體過(guò)程為：

s3.1、建立懸掛系統(tǒng)的非線(xiàn)性系統(tǒng)模型：

θk～p(θk|θk-1)

xk＝f(xk-1,uk-1,wk-1)＝a(θk-1)xk-1+b(θk-1)uk-1+e(θk-1)wk-1

yk＝h(xk,uk,vk)＝c(θk)xk+d(θk)uk+f(θk)vk

其中，θk為第k個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)，p(·|θk-1)表示已知θk-1時(shí)的概率密度函數(shù)，xk為第k個(gè)采樣點(diǎn)的狀態(tài)變量,yk為第k個(gè)采樣點(diǎn)的觀測(cè)變量，uk為第k個(gè)采樣點(diǎn)的系統(tǒng)輸入，wk和vk分別為第k個(gè)采樣點(diǎn)的過(guò)程噪聲和觀測(cè)噪聲，f(·)為非線(xiàn)性狀態(tài)轉(zhuǎn)移方程，h(·)為非線(xiàn)性觀測(cè)方程；

系統(tǒng)待辨識(shí)變量ξk為：

ξk＝[xk^tθk^t]^t

由貝葉斯定理可得：

p(ξk|zk)＝p(xk,θk|zk)＝p(xk|θk,zk)p(θk|zk)

其中，zk為第k個(gè)采樣點(diǎn)的系統(tǒng)觀測(cè)值；

s3.2、初始化第1個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)θ1|0(i)：

如果已知初始概率密度函數(shù)為p(θ0|z0)，對(duì)其進(jìn)行采樣可以得到作為參數(shù)粒子初始值的第1個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)θ1|0(i)；如果沒(méi)有先驗(yàn)知識(shí)作為依據(jù)獲得p(θ0|z0)，則通過(guò)在參數(shù)取值范圍[θmin，θmax]內(nèi)均勻采樣獲得θ1|0(i)；

設(shè)系統(tǒng)初始狀態(tài)值為對(duì)應(yīng)協(xié)方差矩陣為p0，對(duì)進(jìn)行采樣，獲得系統(tǒng)狀態(tài)的預(yù)測(cè)初值x1|0(i)，并根據(jù)下式計(jì)算對(duì)應(yīng)的協(xié)方差矩陣p1|0(i)：

s3.3、依次對(duì)k＝2，k＝3,…,k＝t時(shí)，第k個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)進(jìn)行辨識(shí)，具體包括：

s3.3.1、權(quán)值更新及歸一化：

根據(jù)第k個(gè)采樣點(diǎn)的狀態(tài)預(yù)測(cè)粒子xk|k-1(i)計(jì)算對(duì)應(yīng)的觀測(cè)值yk|k-1(i)，根據(jù)yk|k-1(i)與實(shí)際觀測(cè)值yk的偏差計(jì)算粒子權(quán)重并進(jìn)行歸一化，得到αk(i)

yk|k-1(i)＝c(θk|k-1(i))xk|k-1(i)+d(θk|k-1(i))uk

rk(i)＝c(θk|k-1(i))pk|k-1(i)c^t(θk|k-1(i))+qv

其中，qv為系統(tǒng)觀測(cè)噪聲的協(xié)方差矩陣；

s3.3.2、參數(shù)辨識(shí)：

計(jì)算得到第k個(gè)采樣點(diǎn)的參數(shù)辨識(shí)值為

s3.3.3、重采樣：

對(duì)粒子{θk|k-1(i),xk|k-1(i),pk|k-1(i):i＝1,2,…,n}進(jìn)行重采樣，得到對(duì)應(yīng)的采樣結(jié)果為{θk(j),xk|k-1(j),pk|k-1(j):j＝1,2,…,n}，滿(mǎn)足pr{θk(j)＝θk|k-1(i)}＝αk(i)；

s3.3.4、卡爾曼濾波測(cè)量更新

計(jì)算第k個(gè)采樣點(diǎn)的卡爾曼濾波增益矩陣kk(i)、狀態(tài)變量的辨識(shí)結(jié)果xk(i)和均方誤差矩陣的辨識(shí)結(jié)果pk(i)：

rk(i)＝c(θk(i))pk|k-1(i)c^t(θk(i))+qv

xk(i)＝xk|k-1(i)+kk(i)(yk-c(θk(i))xk|k-1(i))

pk(i)＝pk|k-1(i)-kk(i)c(θk(i))pk|k-1(i)；

s3.3.5、粒子濾波預(yù)測(cè)更新：

粒子濾波預(yù)測(cè)更新過(guò)程應(yīng)用參數(shù)演化中心平滑模型實(shí)現(xiàn)，具體為

其中，a＝(3δ-1)/2δ；δ為折扣因子，取值范圍為(0,1]，一般情況下取值為0.95～0.99；為第k個(gè)采樣點(diǎn)的各參數(shù)粒子的蒙特卡洛均值，為噪聲，且有vk為第k個(gè)采樣點(diǎn)的參數(shù)粒子的方差矩陣，

s3.3.6、卡爾曼濾波預(yù)測(cè)更新：

根據(jù)卡爾曼濾波預(yù)測(cè)更新的結(jié)果以及粒子濾波預(yù)測(cè)更新的結(jié)果計(jì)算下一時(shí)刻的狀態(tài)預(yù)測(cè)值xk+1|k(i)以及均方誤差矩陣的預(yù)測(cè)值pk+1|k(i)：

xk+1|k(i)＝a(θk+1|k(i))xk(i)+b(θk+1|k(i))uk

pk+1|k(i)＝a(θk+1|k(i))pk^t(i)a^t(θk+1|k(i))+e(θk+1|k(i))qwe^t(θk+1|k(i))

其中，qw為系統(tǒng)過(guò)程噪聲的協(xié)方差矩陣；

s3.4、在各采樣點(diǎn)，根據(jù)參數(shù)辨識(shí)結(jié)果跟蹤抗蛇行減震器的性能參數(shù)衰減，進(jìn)而實(shí)現(xiàn)故障辨識(shí)。

優(yōu)選地，步驟s3.4的具體方法為：設(shè)置故障閾值，將各采樣點(diǎn)對(duì)應(yīng)的性能參數(shù)值與故障閾值比較，若采樣點(diǎn)對(duì)應(yīng)的性能參數(shù)值小于等于故障閾值則判斷抗蛇行減震器在該采樣點(diǎn)的采樣時(shí)刻出現(xiàn)故障，其中，設(shè)性能參數(shù)的正常值為θnormal，則故障閾值θf(wàn)ault設(shè)置為θf(wàn)ault＝0.5θnormal。

優(yōu)選地，步驟s3.3還包括：將系統(tǒng)觀測(cè)輸出等間隔分成若干段，在完成基于第k`-1段觀測(cè)數(shù)據(jù)的參數(shù)辨識(shí)后，以第k`-1段觀測(cè)數(shù)據(jù)的參數(shù)辨識(shí)中最后一次參數(shù)辨識(shí)的結(jié)果作為參考值，在范圍內(nèi)均勻采樣重新生成參數(shù)粒子，作為第k`段參數(shù)辨識(shí)過(guò)程的參數(shù)粒子初始值。

本發(fā)明的有益效果如下：

本發(fā)明所述技術(shù)方案將基于非線(xiàn)性濾波的參數(shù)辨識(shí)算法應(yīng)用于高速列車(chē)抗蛇行減震器性能參數(shù)及故障的辨識(shí)中。通過(guò)建立高速列車(chē)懸掛系統(tǒng)的整車(chē)模型，克服了模型簡(jiǎn)化程度過(guò)高的問(wèn)題，使得參數(shù)辨識(shí)過(guò)程更加貼近高速列車(chē)運(yùn)行的實(shí)際過(guò)程。在此基礎(chǔ)上，本發(fā)明將再次均勻采樣策略引入邊緣化粒子濾波算法，使參數(shù)粒子重新獲得多樣性，有效克服了算法收斂后的粒子貧化問(wèn)題。通過(guò)對(duì)突發(fā)故障情況下的抗蛇行減震器性能參數(shù)進(jìn)行辨識(shí)，可以實(shí)現(xiàn)性能參數(shù)的跟蹤監(jiān)測(cè)，從而為車(chē)輛維修部門(mén)提供運(yùn)維保障的決策支持。

附圖說(shuō)明

下面結(jié)合附圖對(duì)本發(fā)明的具體實(shí)施方式作進(jìn)一步詳細(xì)的說(shuō)明；

圖1示出基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法的流程圖。

圖2示出高速列車(chē)懸掛系統(tǒng)結(jié)構(gòu)示意圖，其中，①為軸箱彈簧、②為一系垂向減震器、③為空氣彈簧、④為抗蛇行減震器、⑤為二系橫向減震器。

圖3示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型俯視圖。

圖4示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型俯視圖。

圖5示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型各自由度運(yùn)動(dòng)的正方向示意圖。

圖6示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型的仿真模型示意圖。

圖7示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型零極點(diǎn)分布示意圖。

圖8示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型車(chē)體子系統(tǒng)的零極點(diǎn)分布示意圖。

圖9示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型轉(zhuǎn)向架子系統(tǒng)的零極點(diǎn)分布示意圖。

圖10示出高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型轉(zhuǎn)向架子系統(tǒng)的零極點(diǎn)分布示意圖。

圖11示出本發(fā)明中的邊緣化粒子濾波算法流程流程圖。

圖12示出再次均勻采樣策略原理示意圖。

圖13示出抗蛇行減震器故障設(shè)置圖。

圖14示出抗蛇行減震器阻尼突發(fā)故障情況下的性能參數(shù)辨識(shí)結(jié)果示意圖。

具體實(shí)施方式

為了更清楚地說(shuō)明本發(fā)明，下面結(jié)合優(yōu)選實(shí)施例和附圖對(duì)本發(fā)明做進(jìn)一步的說(shuō)明。附圖中相似的部件以相同的附圖標(biāo)記進(jìn)行表示。本領(lǐng)域技術(shù)人員應(yīng)當(dāng)理解，下面所具體描述的內(nèi)容是說(shuō)明性的而非限制性的，不應(yīng)以此限制本發(fā)明的保護(hù)范圍。

本發(fā)明公開(kāi)的基于非線(xiàn)性濾波的抗蛇行減震器性能參數(shù)及故障辨識(shí)方法(以下將抗蛇行減震器性能參數(shù)簡(jiǎn)述為參數(shù))采用基于非線(xiàn)性濾波的參數(shù)辨識(shí)算法和再次均勻采樣策略，實(shí)現(xiàn)對(duì)抗蛇行減震器進(jìn)行性能參數(shù)辨識(shí)，進(jìn)而實(shí)現(xiàn)對(duì)抗蛇行減震器進(jìn)行故障辨識(shí)。如圖1所示，該方法包括如下步驟：

s1、建立車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型；

s2、根據(jù)車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型，建立懸掛系統(tǒng)的離散狀態(tài)方程和離散觀測(cè)方程；

s3、根據(jù)懸掛系統(tǒng)的離散狀態(tài)方程和離散觀測(cè)方程，基于邊緣化粒子濾波算法對(duì)懸掛系統(tǒng)的抗蛇行減震器進(jìn)行性能參數(shù)辨識(shí)和故障辨識(shí)。

其中，步驟s1的具體過(guò)程為：

如圖2所示，車(chē)輛懸掛系統(tǒng)是指位于車(chē)體與轉(zhuǎn)向架及轉(zhuǎn)向架與輪對(duì)之間，起到支撐、緩沖和引導(dǎo)作用的一系列彈簧和阻尼元件，一般包括兩系懸掛系統(tǒng)。一系懸掛系統(tǒng)位于轉(zhuǎn)向架構(gòu)架與輪對(duì)之間，起到緩沖輪軌沖擊、支撐轉(zhuǎn)向架和引導(dǎo)車(chē)輛運(yùn)行的作用，主要包括軸箱、軸箱彈簧和垂向減震器等；二系懸掛系統(tǒng)位于車(chē)體與轉(zhuǎn)向架構(gòu)架之間，起到緩沖高頻振動(dòng)、支撐車(chē)體和引導(dǎo)車(chē)輛運(yùn)行的作用，主要包括空氣彈簧、橫向減震器和抗蛇行減震器等。

本發(fā)明對(duì)二系懸掛系統(tǒng)抗蛇行減震器的故障辨識(shí)問(wèn)題開(kāi)展研究，首先要建立懸掛系統(tǒng)的整體模型。車(chē)輛系統(tǒng)是復(fù)雜的非線(xiàn)性多體系統(tǒng)，車(chē)體、轉(zhuǎn)向架構(gòu)架和輪對(duì)的剛度要比懸掛系統(tǒng)大很多，因此可以不考慮其彈性，并將車(chē)輛系統(tǒng)簡(jiǎn)化為多剛體系統(tǒng)。分別考慮車(chē)體和轉(zhuǎn)向架構(gòu)架的橫移、搖頭和側(cè)滾運(yùn)動(dòng)以及輪對(duì)的橫移和搖頭運(yùn)動(dòng)，建立車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型：

首先，建立車(chē)輛懸掛系統(tǒng)狀態(tài)空間模型：

y＝cx+du

其中，

狀態(tài)變量x為

觀測(cè)變量y為

其中，yc、φc和ρc分別表示車(chē)體的橫向運(yùn)動(dòng)位移、搖頭角和側(cè)滾角，yb1、φb1和ρb1分別表示前轉(zhuǎn)向架的橫向運(yùn)動(dòng)位移、搖頭角和側(cè)滾角，yb2、φb2和ρb2分別表示后轉(zhuǎn)向架的橫向運(yùn)動(dòng)位移、搖頭角和側(cè)滾角，yw1和φw1分別表示前轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)位移和搖頭角，yw2和φw2分別表示前轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)位移和搖頭角，yw3和φw3分別表示后轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)位移和搖頭角，yw4和φw4分別表示后轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)位移和搖頭角，和分別表示車(chē)體的橫向運(yùn)動(dòng)速度、搖頭運(yùn)動(dòng)角速度和側(cè)滾運(yùn)動(dòng)角速度，和分別表示前轉(zhuǎn)向架的橫向運(yùn)動(dòng)速度、搖頭運(yùn)動(dòng)角速度和側(cè)滾運(yùn)動(dòng)角速度，和分別表示后轉(zhuǎn)向架的橫向運(yùn)動(dòng)速度、搖頭運(yùn)動(dòng)角速度和側(cè)滾運(yùn)動(dòng)角速度，和分別表示前轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)速度和搖頭運(yùn)動(dòng)角速度，和分別表示前轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)速度和搖頭運(yùn)動(dòng)角速度，和分別表示后轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)速度和搖頭運(yùn)動(dòng)角速度，和分別表示后轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)速度和搖頭運(yùn)動(dòng)角速度，和分別表示車(chē)體的橫向運(yùn)動(dòng)加速度、搖頭運(yùn)動(dòng)角加速度和側(cè)滾運(yùn)動(dòng)角加速度，和分別表示前轉(zhuǎn)向架的橫向運(yùn)動(dòng)加速度、搖頭運(yùn)動(dòng)角加速度和側(cè)滾運(yùn)動(dòng)角加速度，和分別表示后轉(zhuǎn)向架的橫向運(yùn)動(dòng)加速度、搖頭運(yùn)動(dòng)角加速度和側(cè)滾運(yùn)動(dòng)角加速度，和分別表示前轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)加速度和搖頭運(yùn)動(dòng)角加速度，和分別表示前轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)加速度和搖頭運(yùn)動(dòng)角加速度，和分別表示后轉(zhuǎn)向架前側(cè)輪對(duì)的橫向運(yùn)動(dòng)加速度和搖頭運(yùn)動(dòng)角加速度，和分別表示后轉(zhuǎn)向架后側(cè)輪對(duì)的橫向運(yùn)動(dòng)加速度和搖頭運(yùn)動(dòng)角加速度；

u為軌道橫向不平順激勵(lì)，作用于車(chē)輛的四組輪對(duì)，u＝[ua1ua2ua3ua4]^t，ua1表示前轉(zhuǎn)向架前側(cè)輪對(duì)的軌道橫向不平順激勵(lì)，ua2表示前轉(zhuǎn)向架后側(cè)輪對(duì)的軌道橫向不平順激勵(lì)，ua3表示后轉(zhuǎn)向架前側(cè)輪對(duì)的軌道橫向不平順激勵(lì)，ua4表示后轉(zhuǎn)向架后側(cè)輪對(duì)的軌道橫向不平順激勵(lì)；

a、b、c和d分別為懸掛系統(tǒng)矩陣，可由運(yùn)動(dòng)體各自由度運(yùn)動(dòng)的微分方程得到；

引入系統(tǒng)過(guò)程噪聲w和觀測(cè)噪聲v，對(duì)車(chē)輛懸掛系統(tǒng)狀態(tài)空間模型進(jìn)行變換，得到車(chē)輛懸掛系統(tǒng)橫向動(dòng)力學(xué)模型：

y＝cx+du+fv＝cx+d(u+b^-1ew)+fv-db^-1ew

其中，e和f分別為過(guò)程噪聲矩陣和觀測(cè)噪聲矩陣，u+b^-1ew作為新的系統(tǒng)輸入，實(shí)現(xiàn)過(guò)程噪聲和觀測(cè)噪聲的引入。

圖3和圖4分別示出了懸掛系統(tǒng)橫向動(dòng)力學(xué)模型的俯視圖和后視圖。為方便建立模型，指定運(yùn)動(dòng)體各自由度運(yùn)動(dòng)的正方向，如圖5所示。

如圖6所示，可在simulink中搭建高速列車(chē)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型的仿真模型，并對(duì)整車(chē)系統(tǒng)、車(chē)體子系統(tǒng)、前/后轉(zhuǎn)向架構(gòu)架子系統(tǒng)和各位輪對(duì)子系統(tǒng)的零極點(diǎn)分布進(jìn)行分析。

步驟s2的具體過(guò)程為：

基于非線(xiàn)性濾波的參數(shù)辨識(shí)是在離散域進(jìn)行的，因此在進(jìn)行參數(shù)辨識(shí)前需要先對(duì)連續(xù)系統(tǒng)進(jìn)行離散化。

懸掛系統(tǒng)的連續(xù)狀態(tài)方程為

其中，x(t)為狀態(tài)變量，y(t)為輸出變量，u(t)為系統(tǒng)輸入，w(t)和v(t)分別為過(guò)程噪聲和觀測(cè)噪聲，和分別為連續(xù)狀態(tài)方程中懸掛系統(tǒng)矩陣，和分別為連續(xù)狀態(tài)方程中過(guò)程噪聲矩陣和觀測(cè)噪聲矩陣。

連續(xù)系統(tǒng)離散化主要是對(duì)描述系統(tǒng)動(dòng)態(tài)特性的狀態(tài)方程而言，輸出方程為靜態(tài)的代數(shù)方程，離散化后應(yīng)保持不變，即

基于連續(xù)系統(tǒng)狀態(tài)方程的解，可得到懸掛系統(tǒng)的離散狀態(tài)方程為

xk+1＝axk+buk+ewk

其中，xk為第k個(gè)采樣點(diǎn)的狀態(tài)變量,uk為第k個(gè)采樣點(diǎn)的系統(tǒng)輸入，wk為第k個(gè)采樣點(diǎn)的過(guò)程噪聲，離散狀態(tài)方程中懸掛系統(tǒng)矩陣離散狀態(tài)方程中過(guò)程噪聲矩陣

懸掛系統(tǒng)的離散觀測(cè)方程為：

yk＝cxk+duk+fvk

其中，離散狀態(tài)方程中懸掛系統(tǒng)矩陣和為連續(xù)狀態(tài)方程中懸掛系統(tǒng)矩陣，離散狀態(tài)方程中觀測(cè)噪聲矩陣為連續(xù)狀態(tài)方程中觀測(cè)噪聲矩陣，vk為第k個(gè)采樣點(diǎn)的觀測(cè)噪聲。

步驟s3的具體過(guò)程為：

s3.1、建立懸掛系統(tǒng)的非線(xiàn)性系統(tǒng)模型：

θk～p(θk|θk-1)

xk＝f(xk-1,uk-1,wk-1)＝a(θk-1)xk-1+b(θk-1)uk-1+e(θk-1)wk-1

yk＝h(xk,uk,vk)＝c(θk)xk+d(θk)uk+f(θk)vk

其中，θk為第k個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)，p(·|θk-1)表示已知θk-1時(shí)的概率密度函數(shù)，xk為第k個(gè)采樣點(diǎn)的狀態(tài)變量,yk為第k個(gè)采樣點(diǎn)的觀測(cè)變量，uk為第k個(gè)采樣點(diǎn)的系統(tǒng)輸入，wk和vk分別為第k個(gè)采樣點(diǎn)的過(guò)程噪聲和觀測(cè)噪聲，f(·)為非線(xiàn)性狀態(tài)轉(zhuǎn)移方程，h(·)為非線(xiàn)性觀測(cè)方程。

將系統(tǒng)待辨識(shí)變量ξk分為兩部分：

ξk＝[xk^tθk^t]^t

由貝葉斯定理可得

p(ξk|zk)＝p(xk,θk|zk)＝p(xk|θk,zk)p(θk|zk)

其中，zk為第k個(gè)采樣點(diǎn)的系統(tǒng)觀測(cè)值，p(xk|θk,zk)易于求出解析解，可通過(guò)卡爾曼濾波器實(shí)現(xiàn)。而p(θk|zk)無(wú)法直接求得解析解，本發(fā)明中通過(guò)粒子濾波器實(shí)現(xiàn)。則系統(tǒng)待辨識(shí)變量ξk的最小均方誤差估計(jì)為

s3.2、初始化第1個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)θ1|0(i)：

如果已知初始概率密度函數(shù)為p(θ0|z0)，對(duì)其進(jìn)行采樣可以得到作為參數(shù)粒子初始值的第1個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)θ1|0(i)；如果沒(méi)有先驗(yàn)知識(shí)作為依據(jù)獲得p(θ0|z0)，則通過(guò)在參數(shù)取值范圍[θmin，θmax]內(nèi)均勻采樣獲得θ1|0(i)；

設(shè)系統(tǒng)初始狀態(tài)值為對(duì)應(yīng)協(xié)方差矩陣為p0，對(duì)進(jìn)行采樣，獲得系統(tǒng)狀態(tài)的預(yù)測(cè)初值x1|0(i)，并根據(jù)下式計(jì)算對(duì)應(yīng)的協(xié)方差矩陣p1|0(i)：

s3.3、依次對(duì)k＝2，k＝3,…,k＝t時(shí)，第k個(gè)采樣點(diǎn)的懸掛系統(tǒng)的抗蛇行減震器性能參數(shù)進(jìn)行辨識(shí)，具體包括：

s3.3.1、權(quán)值更新及歸一化：

根據(jù)第k個(gè)采樣點(diǎn)的狀態(tài)預(yù)測(cè)粒子xk|k-1(i)計(jì)算對(duì)應(yīng)的觀測(cè)值yk|k-1(i)，根據(jù)yk|k-1(i)與實(shí)際觀測(cè)值yk的偏差計(jì)算粒子權(quán)重并進(jìn)行歸一化，得到αk(i)

yk|k-1(i)＝c(θk|k-1(i))xk|k-1(i)+d(θk|k-1(i))uk

rk(i)＝c(θk|k-1(i))pk|k-1(i)c^t(θk|k-1(i))+qv

其中，qv為系統(tǒng)觀測(cè)噪聲的協(xié)方差矩陣；

s3.3.2、參數(shù)辨識(shí)：

計(jì)算得到第k個(gè)采樣點(diǎn)的參數(shù)辨識(shí)值為

s3.3.3、重采樣(粒子濾波測(cè)量更新)：

對(duì)粒子{θk|k-1(i),xk|k-1(i),pk|k-1(i):i＝1,2,…,n}進(jìn)行重采樣，得到對(duì)應(yīng)的采樣結(jié)果為{θk(j),xk|k-1(j),pk|k-1(j):j＝1,2,…,n}，滿(mǎn)足pr{θk(j)＝θk|k-1(i)}＝αk(i)；

s3.3.4、卡爾曼濾波測(cè)量更新

計(jì)算第k個(gè)采樣點(diǎn)的卡爾曼濾波增益矩陣kk(i)、狀態(tài)變量的辨識(shí)結(jié)果xk(i)和均方誤差矩陣的辨識(shí)結(jié)果pk(i)：

rk(i)＝c(θk(i))pk|k-1(i)c^t(θk(i))+qv

xk(i)＝xk|k-1(i)+kk(i)(yk-c(θk(i))xk|k-1(i))

pk(i)＝pk|k-1(i)-kk(i)c(θk(i))pk|k-1(i)；

s3.3.5、粒子濾波預(yù)測(cè)更新：

粒子濾波預(yù)測(cè)更新過(guò)程應(yīng)用參數(shù)演化中心平滑模型實(shí)現(xiàn)，具體為

其中，a＝(3δ-1)/2δ；δ為折扣因子，取值范圍為(0,1]，一般情況下取值為0.95～0.99；為第k個(gè)采樣點(diǎn)的各參數(shù)粒子的蒙特卡洛均值，為噪聲，且有vk為第k個(gè)采樣點(diǎn)的參數(shù)粒子的方差矩陣，

s3.3.6、卡爾曼濾波預(yù)測(cè)更新：

根據(jù)卡爾曼濾波預(yù)測(cè)更新的結(jié)果以及粒子濾波預(yù)測(cè)更新的結(jié)果計(jì)算下一時(shí)刻的狀態(tài)預(yù)測(cè)值xk+1|k(i)以及均方誤差矩陣的預(yù)測(cè)值pk+1|k(i)：

xk+1|k(i)＝a(θk+1|k(i))xk(i)+b(θk+1|k(i))uk

其中，qw為系統(tǒng)過(guò)程噪聲的協(xié)方差矩陣。

s3.4、在各采樣點(diǎn)，根據(jù)參數(shù)辨識(shí)結(jié)果跟蹤抗蛇行減震器的性能參數(shù)衰減，進(jìn)而實(shí)現(xiàn)故障辨識(shí)，具體方法為：設(shè)置故障閾值，將各采樣點(diǎn)對(duì)應(yīng)的性能參數(shù)值與故障閾值比較，若采樣點(diǎn)對(duì)應(yīng)的性能參數(shù)值小于等于故障閾值則判斷抗蛇行減震器在該采樣點(diǎn)的采樣時(shí)刻出現(xiàn)故障，其中，設(shè)性能參數(shù)的正常值為θnormal，則故障閾值θf(wàn)ault可設(shè)置為θf(wàn)ault＝0.5θnormal。

本發(fā)明采用的非線(xiàn)性濾波算法為邊緣化粒子濾波算法，將狀態(tài)辨識(shí)過(guò)程和參數(shù)辨識(shí)過(guò)程進(jìn)行有效分離，分別采用卡爾曼濾波和粒子濾波進(jìn)行辨識(shí)。當(dāng)邊緣化粒子濾波算法的參數(shù)辨識(shí)結(jié)果收斂后，所有參數(shù)粒子的權(quán)重基本相同，并且粒子值接近參數(shù)的正常值。此時(shí)即使抗蛇行減震器發(fā)生突發(fā)故障，由于參數(shù)粒子已經(jīng)貧化，無(wú)法繼續(xù)逼近參數(shù)的故障值，導(dǎo)致算法不能及時(shí)對(duì)突發(fā)故障做出反應(yīng)。因此，實(shí)現(xiàn)故障辨識(shí)的關(guān)鍵是使參數(shù)粒子重新獲得多樣性，并在系統(tǒng)觀測(cè)輸出的修正下持續(xù)逼近故障值。因此本發(fā)明引入再次均勻采樣策略，將參數(shù)辨識(shí)過(guò)程分割為若干段，每一段參數(shù)辨識(shí)過(guò)程結(jié)束后，對(duì)參數(shù)粒子進(jìn)行均勻采樣，使參數(shù)粒子重新獲得多樣性，以解決粒子貧化問(wèn)題，步驟s3.3還包括：將系統(tǒng)觀測(cè)輸出等間隔分成若干段(即，將總采樣點(diǎn)分為多組)，在完成基于第k`-1段觀測(cè)數(shù)據(jù)的參數(shù)辨識(shí)后(即，在完成第k`-1組的參數(shù)辨識(shí)后)，以第k`-1段觀測(cè)數(shù)據(jù)的參數(shù)辨識(shí)中最后一次參數(shù)辨識(shí)的結(jié)果作為參考值(即，以第k`-1組的參數(shù)辨識(shí)中最后一次參數(shù)辨識(shí)的結(jié)果作為參考值)，在范圍內(nèi)均勻采樣重新生成參數(shù)粒子，作為第k`段參數(shù)辨識(shí)過(guò)程的參數(shù)粒子初始值。

下面代入具體設(shè)置作進(jìn)一步說(shuō)明：

設(shè)置系統(tǒng)觀測(cè)輸出分為4段，分別為階段1、階段2、階段3和階段4，系統(tǒng)觀測(cè)輸出的分段長(zhǎng)度t＝250，粒子數(shù)ns＝1000，采樣間隔ts＝0.5ms，被辨識(shí)參數(shù)的正常值為θnormal，故障值為θf(wàn)ault＝0.5θnormal，均勻采樣的參考值為通過(guò)在范圍內(nèi)均勻采樣獲得初始參數(shù)粒子。假設(shè)車(chē)體前側(cè)抗蛇行減震器發(fā)生突發(fā)故障，后側(cè)抗蛇行減震器正常。同時(shí)對(duì)故障和工況正常抗蛇行減震器進(jìn)行參數(shù)辨識(shí)：

在simulink中搭建懸掛系統(tǒng)橫向動(dòng)力學(xué)模型的仿真模型后，設(shè)置車(chē)輛運(yùn)行速度為80m/s，系統(tǒng)采樣頻率為2000hz，系統(tǒng)采樣間隔為0.5ms。對(duì)懸掛系統(tǒng)橫向動(dòng)力學(xué)模型的零極點(diǎn)分布進(jìn)行分析，整車(chē)系統(tǒng)、車(chē)體子系統(tǒng)、前/后轉(zhuǎn)向架構(gòu)架子系統(tǒng)和各位輪對(duì)子系統(tǒng)的零極點(diǎn)分布分別如圖7至圖10所示。如果各級(jí)系統(tǒng)的零極點(diǎn)均位于虛軸左側(cè)，說(shuō)明系統(tǒng)是穩(wěn)定的。

本發(fā)明采用的非線(xiàn)性濾波算法是邊緣化粒子濾波算法，該算法的基本流程如圖11所示。在此基礎(chǔ)上，引入在此均勻采樣策略，基本原理如圖12所示。

再次均勻采樣策略的基本過(guò)程是：將系統(tǒng)觀測(cè)輸出等間隔分成若干段，在完成基于第k`-1段觀測(cè)數(shù)據(jù)的參數(shù)辨識(shí)后，以參數(shù)辨識(shí)結(jié)果作為參考值，并在范圍內(nèi)均勻采樣重新生成參數(shù)粒子，作為第k`段參數(shù)辨識(shí)過(guò)程的初始參數(shù)粒子。由圖12可知，算法對(duì)突發(fā)故障的反應(yīng)時(shí)間由兩方面因素決定：一是算法本身的辨識(shí)速度，二是系統(tǒng)觀測(cè)輸出的分段長(zhǎng)度。當(dāng)算法辨識(shí)速度較快時(shí)，需要的迭代次數(shù)較少，因此可以將系統(tǒng)觀測(cè)輸出的分段長(zhǎng)度設(shè)置得盡量短，使算法可以很快辨識(shí)出突變后的參數(shù)。

各階段的具體情況如表二所示：

表二故障辨識(shí)各階段的設(shè)置情況

具體說(shuō)明如下：

(1)階段2在第126個(gè)采樣點(diǎn)處發(fā)生突發(fā)故障，參數(shù)由正常值θnormal跳變?yōu)楣收祥撝郸萬(wàn)ault；

(2)階段2開(kāi)始于正常工況，結(jié)束于故障工況，如果辨識(shí)結(jié)果在參數(shù)突變前已經(jīng)收斂，則粒子貧化現(xiàn)象將導(dǎo)致算法不能及時(shí)對(duì)突發(fā)故障做出反應(yīng)，因此階段2的參數(shù)辨識(shí)結(jié)果理論上應(yīng)收斂于正常值θnormal；

(3)階段3均勻采樣的參考值理論上應(yīng)接近參數(shù)的正常值θnormal，并且如果算法有效，參數(shù)辨識(shí)結(jié)果理論上應(yīng)該收斂于參數(shù)的故障值θf(wàn)ault。

如圖13所示，對(duì)抗蛇行減震器進(jìn)行故障仿真，設(shè)置車(chē)體前側(cè)的兩個(gè)抗蛇行減震器發(fā)生突發(fā)故障，車(chē)體后側(cè)的兩個(gè)抗蛇行減震器處于正常狀態(tài)。同時(shí)對(duì)故障和正常的抗蛇行減震器進(jìn)行參數(shù)辨識(shí)，結(jié)果如圖14所示。由圖可知，對(duì)正常的抗蛇行減震器而言，參數(shù)辨識(shí)結(jié)果能夠收斂在正常值附近；對(duì)故障的抗蛇行減震器而言，參數(shù)辨識(shí)結(jié)果能夠跟蹤參數(shù)跳變，并且辨識(shí)值能夠收斂在真實(shí)值附近。綜上所述，基于再次均勻采樣策略的邊緣化粒子濾波算法能夠?qū)崿F(xiàn)突發(fā)故障情況下性能參數(shù)的跟蹤監(jiān)測(cè)，從而實(shí)現(xiàn)故障辨識(shí)。

顯然，本發(fā)明的上述實(shí)施例僅僅是為清楚地說(shuō)明本發(fā)明所作的舉例，而并非是對(duì)本發(fā)明的實(shí)施方式的限定，對(duì)于所屬領(lǐng)域的普通技術(shù)人員來(lái)說(shuō)，在上述說(shuō)明的基礎(chǔ)上還可以做出其它不同形式的變化或變動(dòng)，這里無(wú)法對(duì)所有的實(shí)施方式予以窮舉，凡是屬于本發(fā)明的技術(shù)方案所引伸出的顯而易見(jiàn)的變化或變動(dòng)仍處于本發(fā)明的保護(hù)范圍之列。