基于模糊聚類的高速列車車速識別

翁 順 吳俐瀅 朱宏平 張 敏 曾永平 謝 毅 蔡劍橋 顏永逸

(1.華中科技大學，武漢 430074；2.中鐵二院工程集團有限責任公司，成都 610031)

隨著我國高速鐵路的快速發展，列車高速運行對橋梁動力性能提出了更高的要求。高速列車運行會對橋梁產生低周疲勞荷載作用，列車荷載是影響橋梁使用壽命的主要因素之一。列車的車速和荷載識別可為橋梁的安全評估與檢定提供重要依據。

目前，橋梁移動荷載識別研究的主要方法可以分為3類：

(1)在橋面上布設動態稱重系統(WIM)來實現對車輛荷載的在線監測。此種方法優點是適用性廣，可用于不同形式的橋梁結構，且可長時間采集數據[1]；缺點是費用昂貴且所得的測量結果與實際值偏差較大。

(2)利用車-橋耦合系統振動方程求解任意時刻車輛與橋梁接觸處的相互作用力，并將其用以識別移動荷載。如1988年O’Connor和Chan[2]將簡支梁簡化為由無質量彈性梁連接的質量集中在節點的模型，將橋上的移動荷載、慣性力、阻尼力等視為集中力，推導出的識別車橋相互作用力的計算方法，即解析法Ⅰ(IM Ⅰ)。Law和Chan基于系統識別理論提出了兩種荷載識別方法：時域法(TDM)[3]和頻時域法(FTDM)[4]；隨后與解析法Ⅰ相類似的解析法Ⅱ[5]也被提出。此類方法需要的系統參數多，識別結果(即動態接觸力)受路面不平順、行駛速度等因素影響大[6]。

(3)根據動力響應如動應變等直接識別橋上車輛的軸數、軸距、軸重和車速等。20世紀70年代Moses[7]等利用單一移動荷載的動響應和測點彎矩影響線成比例的特性，由跨中應變響應提出了最初的BWIM移動荷載識別方法；Peters[8]在已有BWIM的基礎上通過動應變曲線的面積得到總重量，并結合最小二乘原理迭代計算分配各軸重；王寧波等通過尋找應變曲線峰值點識別車輛行駛速度、軸數、軸距，并采用影響線擬合動應變響應的思路識別車輛各軸軸重。此類荷載識別方法在公路橋梁中廣泛應用，待識別信息量大，實時性有待提高。

與公路橋的車速識別不同，高速鐵路列車的行車速度有嚴格的要求，常固定不同的車速區間，因此判斷和識別高速車速時需要更高的精度。本文依據高速鐵路橋梁健康監測系統采集的橋梁動力響應監測數據，提出基于模糊聚類分析的高速鐵路列車車速識別方法，建立列車經過橋梁時的橋梁加速度響應自回歸模型(AR模型)，通過對自回歸系數的模糊聚類結果來判斷車速的浮動區間。

1 高速列車車速識別基本理論

1.1 加速度響應時間序列模型

基于時間序列中各時間點的觀測值間的相關關系，可依據已知的歷史數據對未來的發展進行預測，時間序列分析是研究此種相關關系的方法[9]。高速列車的加速度響應數據就是一種典型的時間序列。當高速列車通過鐵路橋梁時，列車車速和車輛荷載的激勵將會導致不同車速列車對橋梁產生不同的加速度響應，依據這種差異性便可由模糊聚類分析識別車速。

1.1.1 加速度時間序列自回歸模型

自回歸模型(Auto regressive，簡稱AR)是一種常用的處理時間序列的方法[10]。假設有一零均值的加速度時間序列為{at}，t=1，2，3，…，n，對其進行自回歸擬合，結果如下：

at=φ1at-1+φ2at-2+φ3at-3+…+φpat-p+et

(1)

式中：at——加速度時間序列在t時刻的觀測數據值；

φi(i=1,2,3,…,p)——自回歸系數；

p——AR模型的擬合階數；

式(1)建立的模型可稱為AR(p)模型。

1.2.2 AR(p)模型的階數p

若要建立AR(p)模型，首先需要確定AR(p)模型的階數p，合適的階數p將極大地提高模型的精確度。目前，確定AR(p)模型階數p的方法包括AIC準則、BIC準則和FPE準則等[11]。本文將綜合采用AIC準則和FPE準則確定階數p值。其中，AIC準則的表達式如下：

(2)

式中:N——用于建立AR(p)模型的加速度時間序列的長度；

由式(2)可知，AIC的值主要由兩部分組成，前部分主要與模型的殘差相關，后部分與模型的階數相關。模型越精確，則前部分的值越小，而后部分的值越大。因此，AIC值取得最小值時對應的p值為模型的合適階數。

FPE準則的表達式為：

(3)

FPE值取得最小值時對應的p值為模型的合適階數。

1.2.3 AR(p)模型的參數估計

求得模型的階數p后，由實測的加速度時間序列便可求得對應的AR(p)模型。對AR(p)模型的參數估計的方法可分為兩類：一類是對時間序列{at}進行直接估計，如Yule-Walker法[12]、Ulrych-Clayton法[13]、最小二乘估計[14]等；另一類是對時間序列進行遞推估計，如Burg法[15]、LUD法、Levinson法和遞推最小二乘法等。本文采用Yule-Walker法，具體的求解過程如下：

對式(1)兩邊同乘at-k，k≥0，得：

at-kat=φ1tt-kat-1+φ2tt-kat-2+…+φptt-kat-p+at-ket

(4)

同時對式(4)的兩邊取期望，得：

E[atat-k]=φ1E[at-1at-k]+φ2E[at-2at-k]+

…+φpE[at-pat-k]+E[etat-k]

(5)

其中E[atat-k]為加速度時間序列{at}的自協方差函數：

E[atat-k]=Rk=R-k=E[atat+k]

(6)

Rk=φ1Rk-1+φ2Rk-2+…+φpRk-p,k>0

(7)

兩邊同除R0，得到：

ρk=φ1ρk-1+φ2ρk-2+…+φpρk-p,k>0

(8)

其中ρk=Rk/R0，式(8)可表達為：

(9)

式(9)稱為Yule-Walker方程?？捎洖椋?/p>

(10)

則根據式(10)，式(9)可表達為：

φ=Γ-1ρ

(11)

根據式(11)可得到參數φi(i=1,2,3,…,p)的Yule-Walker估計。

1.2 車速聚類分析理論

聚類分析包括系統聚類法、模糊聚類法、有序樣品聚類法等多種分析算法[16-17]，其中模糊聚類法為其中應用范圍最廣、使用頻率最高。模糊聚類法引入了隸屬度的概念，隸屬度取值為0～1，每個聚類樣本根據隸屬度的大小被分到多個類中。

眾多的模糊聚類算法中使用最為廣泛的是基于目標函數的算法。1974年，由Dunn提出并經Bezdek推廣了模糊C-均值算法(FCM算法)[18-19]，FCM算法建立每個點到聚類中心的距離與隸屬度乘積之和的目標函數，通過對目標函數的迭代優化達到對樣本空間劃分。FCM算法優化的目標函數為：

(12)

式中：ci——第i類的加速度時間序列AR模型系數聚類中心；

rj——第j個AR模型系數聚類樣本；

m——模糊聚類參數，一般取為2[20-21]。

(13)

(14)

通過不斷迭代式(13)和式(14)來可優化目標函數J，并使其達到最優。

FCM算法的迭代過程如下：

(1)設置聚類的個數c和初始隸屬度矩陣，初始隸屬度矩陣的設置只需滿足元素值為0～1且每個樣本對各個聚類中心的隸屬度之和為1即可。

(2)根據設置的初始隸屬度矩陣，由式(13)求各個聚類中心ci，i=1,2,…,c。

(3)根據求得的聚類中心ci，由式(14)求隸屬度矩陣uij。

(4)如果前一次迭代和后一次迭代求得的目標函數的差值滿足條件|J(q)-J(q+1)|<ε則停止迭代，其中q為迭代次數。否則，重復第(2)～第(4)步，直到滿足要求為止，此時得到的聚類中心和隸屬度為滿足要求的聚類中心和隸屬度。

2 高速列車車速識別分析方法

2.1 列車過橋判定

高速列車以不同車速通過橋梁時，橋梁結構的加速度響應和環境激勵所產生的加速度響應間存在明顯差異，其加速度響應時間序列如圖1所示。在該段加速度時間序列中，測點序列為0～900時為無列車通過時橋梁的加速度響應時間序列(區段1)；測點序列為900～1 200 時為列車通過時橋梁的加速度響應時間序列(區段2)；測點序列為 1 200～2 000 時為列車通過后橋梁的加速度響應時間序列(區段3)。

圖1 列車通過時橋梁加速度響應時間序列圖

實際應用中，區段2數據對計算結果的準確性影響最大，因此AR建模的時間序列必須包含區段2，而減少區段1和區段3的數據點數。列車以高速通過橋梁時在橋上的滯留時間很短，受限于采樣頻率區段2的時間序列長度很短，不一定能滿足模糊聚類的長度要求，且可能會降低聚類計算結果的精度。因此，本文選取全部區段2、部分區段1和部分區段3的數據進行分析。

橋梁健康監測系統采集加速度響應數據是實時連續進行的，而列車經過橋梁的時間很短且不定，在列車經過橋梁的起止時間未知的狀況下，如何提取區段2的加速度時間序列成為亟待解決的問題。本文提出了兩種方法用以提取區段2的加速度時間序列：閾值判斷法和樣本區間統計值法。

(1)閾值判斷法需對加速度設定某一閾值，當加速度響應超過閾值時將該時間點前后某一區間視為高速列車通過橋梁的時間。假設無列車通過時的加速度響應統計極值為amax，則設namax為所需閾值，其中n為依據實際情況采用的放大系數。假設有某一加速度時間序列{ai},i=1,2,…,n，經判斷at≥namax，認為{ak}為高速列車通過時橋梁的加速度響應時間序列，其中k=t-x,…,t-1,t,t+1,…,t+y，{ak}?{ai},x和y為根據實際情況采取的時間段距離參數，可令x=y。

由于環境荷載的影響和傳感器的采集誤差，可能出現某一個加速度響應數據點超過閾值，但對應時間點并無列車通過的情況，該點被稱為“壞點”。為區分所需測點與“壞點”，本文采用的方法為：在目標測點前后某一區間內重復判斷其他實測數據是否超過閾值，若存在多個點超過閾值的情況，則可判斷為所需的測點值，否則為“壞點”。

(2)樣本區間統計值法將某一加速度響應時間序列{ai},i=1,2,…,n劃分為多個樣本區間{ak}m，k=1,2,…,l，m=1,2,…,t，{ak}m?{ai}，分別求各個樣本區間的方差，當其中某一個區間的樣本方差與其他區間樣本方差相差較大時(一般相差100倍以上)可認為該區間的數據為所需的區段2中的數據。

與閾值判斷法相比，樣本區間統計值計算量較大，在實時采集數據的監測系統中不一定能適用，因此本文采用閾值判斷法。

2.2 監測數據標準化處理

在對加速度時間序列計算前，需對其進行標準化處理，即提取趨勢項、零均值化和標準化。由于本文加速度時間序列長度并不大，在環境溫度或者其他因素的影響下，不能產生明顯的趨勢項，因此只對其進行零均值化和標準化處理。

(1)零均值化

傳感器可能存在系統誤差，導致監測數據均值并不為0，為了減小這種誤差，剔除監測數據的均值。假設加速度時間序列{at},t=1,2,…,n的均值為μ，則零均值化后的時間序列為：

a′t=at-μ

(15)

(2)標準化

對于監測得到的加速度時間序列{at},t=1,2,…,n，其值過大或者過小都會影響AR建模的效果，因此一般需對其進行標準化處理。假設加速度時間序列{at},t=1,2,…,n的均值為μ，方差為σ2，則標準化處理后的時間序列為：

(16)

經標準化處理后的時間序列{a′t}服從正態分布，即a′t～N(0,1)。

2.3 加速度時間序列AR模型

假設高速列車在已知車速行駛過程中的某一測點的加速度時間序列樣本空間S，S={{a}1,{a}2,…,{a}m-1,{a}m}，其中{a}i為在已知車速Vi下橋梁的加速度時間序列，i=1，2，…，m。同一測點在未知車速下的加速度時間序列樣本空間為D,D={{a}m+1,{a}m+2,…,{a}m+l-1,{a}m+l}，其中{a}j為在未知車速Vj下橋梁的加速度時間序列，其中j=m+1，m+2，…，m+l。此時所有數據已經過預處理。

(1)確定AR模型的階數

由式(2)和式(3)分別計算AIC指標和FPE指標，令模型階數p從1～100變化，觀察兩個指標的變化趨勢，確定p所處的大概范圍，再由MATLAB計算每個p值下AR(p)模型的擬合匹配率，選取匹配率變化拐點時的p值作為本文所建立的AR模型的階數。

(2)建立AR(p)模型

分別建立加速度時間序列{a}i和{a}j的AR(p)模型如下：

(17)

(18)

式中：ai,t——時間序列{a}i中t時刻的值;

aj,t——時間序列{a}j中t時刻的值;

φk——時間序列{a}i的AR模型自回歸系數;

φk——時間序列{a}j的AR模型自回歸系數。

對于樣本空間中的某一個時間序列，稱[φ1φ2…φp]或者[φ1φ2…φp]為一個聚類樣本，分別對應于已知車速狀態和未知車速狀態?；跇颖究臻gS={{a}1,{a}2,…,{a}m-1,{a}m}和樣本空間D={{a}m+1,{a}m+2,…,{a}m+l-1,{a}m+l}，定義兩個不同狀態下的聚類樣本空間(或稱為AR模型系數矩陣)，分別為：

(19)

(20)

為通過聚類識別車速，需將未知車速下的Φd的聚類結果與已知車速下的Φs的聚類結果進行對比，因此在實際的計算過程中，將Φd和Φs合并為總體樣本進行聚類計算，該總體樣本可寫為Φ=[Φs;Φd]的形式，Φ中的每一行即為式(13)中的一個樣本。

2.4 FCM聚類識別車速

FCM聚類識別車速的過程主要包括確定聚類中心個數和未知車速識別。根據求得的AR(p)模型系數矩陣Φs，該矩陣中的AR系數為已知車速條件下求得，每個AR(p)模型分別對應列車車速Vi，i=1,2,…,m，可先由Φs的聚類結果確定合適的聚類中心個數c。此時車速Vi將被劃分為c個區間，在識別未知車速時將以這c個車速區間為參考對象，c越大，即聚類中心個數越多，車速區間越多，識別結果越精確。然而，隨著c的增大，FCM聚類分析的結果越難判斷，因此選用合適的c值尤為重要，本文經研究確定合適的聚類中心個數c為3。采用FCM算法對總體聚類樣本Φ=[Φs;Φd]進行聚類劃分，得到隸屬度矩陣：

(21)

式中，uij表示第j個樣本對第i個聚類中心的隸屬度。由隸屬度矩陣可判斷車速所在區間。

3 某鐵路橋高速列車車速識別

3.1 算例及工況介紹

某鐵路橋梁設計行車速度 250 km/h，預留行車速度350 km/h，正線線間距5.0 m，鋪設無砟軌道。主橋結構全梁長572.1 m(含梁縫)。邊跨及部分中跨主梁為預應力混凝土箱梁，中跨主梁為箱形鋼-混凝土結合梁。主梁采用混合梁結構，梁端至中跨155.75 m范圍內采用混凝土箱梁結構，其余采用箱形鋼-混凝土結合梁，中間采用鋼混結合段過渡。索塔采用人字形混凝土塔，兩座索塔全高分別為124.5 m、127.6 m，橋面以上塔高88 m，為主跨的 1/3.409。主塔采用分離式樁基礎，兩承臺間設置系梁。斜拉索采用采用抗拉標準強度 1 670 MPa 鍍鋅平行鋼絲拉索，空間雙索面體系，扇形布置，全橋共 48 對斜拉索。

實測加速度數據采集頻率為30 Hz，共選取4個加速度測點，監測截面皆為中跨，編號分別為10-ZD02、10-ZD03、11-ZD02和11-ZD03。已知和未知車速下橋梁動力響應數據信息如表1所示。

表1 不同車速下加速度動力響應時間序列信息表

以表1中的不同車速建立工況，其中時間序列窗口數是指監測時間段內列車經過橋梁的次數(即可用于車速聚類的不同加速度響應時間序列窗口個數)。每一個時間序列窗口對應于一種車速工況，則共劃分為25種工況，其中已知車速工況17種，未知車速工況8種。已知和未知車速工況如表2和表3所示，未知車速用Vi表示。

表2 已知車速工況信息表

表3 未知車速工況信息表

3.2 監測數據選擇和預處理

目前橋梁暫處于試運行階段，經過的高速列車16節車廂編組，由此計算得到列車在橋梁上的運行時間為10～30 s，又因加速度數據的采集頻率為30 Hz，為取完所需數據，將時間序列的窗口長度定為 2 000。以工況17為例，根據2.1節介紹的閾值判斷法，4個測點所需的時間序列如圖2所示。

圖2 工況17時間序列圖

截取完所需數據點后，根據2.2節中的方法對數據進行預處理，得到的時間序列如圖3所示(以工況17為例)。

圖3 工況17預處理后的時間序列圖

3.3 加速度時間序列AR(p)模型建立

(1)確定AR(p)模型的階數

取工況5下測點11-ZD03的加速度時間序列來確定模型階數，根據AIC準則和FPE準則計算模型階數，結果如圖4和圖5所示。

由圖4和圖5可知，隨著AR模型階數的增加，AIC和FPE的值都在不停的減小，但均存在明顯的拐點。當模型的階數從1～10和10～100變化時，AIC的值和FPE的值的遞減速率存在較大的差異，模型階數大于10階后兩個指標基本沒有變化。由圖6可知，大于10階后模型匹配率變化較小。綜合考慮，取AR(p)模型的合適階數為10階。

圖4 AIC曲線圖

圖5 FPE曲線圖

圖6 模型匹配率曲線圖

(2)建立AR(p)模型

為了驗證AR模型的準確性和可靠性，計算對應時間窗口長度的時間序列，將原時間序列與求得的時間序列進行比較，由此判斷AR模型是否符合要求。以工況5為例，AR模型的自回歸系數如表4所示。各測點的原始時間序列與模型求得的時間序列比較曲線如圖7所示。由圖7可知，AR建模預測得到的加速度時間序列與實測的加速度時間序列的變化趨勢基本一致，可用于聚類分析。

圖7 工況5各測點數據對比圖

表4 工況5下各測點AR模型系數表

同理，對表2和3中列出的所有工況下的各個測點時間序列進行AR建模并求解AR系數。

3.4 AR系數聚類識別車速

選取聚類中心個數為3進行算例分析。對已知車速的17種工況下的加速度時間序列的AR系數矩陣進行FCM聚類分析，得到已知車速工況下的隸屬度信息。以測點10-ZD03為對象，計算得到的聚類中心個數為3下的隸屬度信息，結果如表5所示，聚類效果如圖8所示。

圖8 聚類中心個數為3時測點10-ZD03聚類效果圖

表5 測點10-ZD02各待測工況隸屬度表

由表5可知，當聚類中心個數為3時，工況1～工況4對類別1的隸屬度均在0.7以上，對類別2和3的隸屬度均在0.18以下，隸屬度的差別較大，可將工況1～工況4劃分為類別1；工況5～工況8對類別2的隸屬度均在0.56以上，大部分在0.67以上，對類別1和3的隸屬度均在0.25以下，大部分在0.20以下，可將工況5～工況8劃分為類別2；工況11～工況17對類別3的隸屬度均在0.45以上，大部分在0.60以上，對類別1和3的隸屬度均在0.35以下，大部分在0.20以下，可將工況11～工況17劃分為類別3。此時類別1對應的車速區間為160～200 km/h，類別2對應的車速區間為318～340 km/h，類別3對應的車速區間為340～380 km/h。

將未知車速和已知車速的所有AR系數組成的AR系數矩陣再次進行FCM聚類分析，得到各未知車速工況下的隸屬度信息如表6～表9所示，聚類效果圖如圖9所示。

表6 測點10-ZD02各待測工況隸屬度表

表7 測點10-ZD03各待測工況隸屬度表

表8 測點11-ZD02各待測工況隸屬度表

圖9 聚類效果圖

表9 測點11-ZD03各待測工況隸屬度表

由表6可知得測點10-ZD02的聚類結果：待測工況19和待測工況24對3個聚類中心的隸屬度皆在0.3左右，因此，待測工況19、工況24可能為3個類別中的任意一類；待測工況21、工況23、工況25對類別2的隸屬度均在0.76以上，對類別1、類別3的隸屬度均在0.17以下，可視為類別2；待測工況18、工況20、工況22對類別3的隸屬度均在0.56以上，對類別1和類別2的隸屬度均在0.27以下，可視為類別3。

由聚類分析結果可知，除了少數幾個工況，其余工況類別劃分基本一致。待測工況19、工況24為類別1時對應的車速區間為160～200 km/h；待測工況21、工況23、工況25為類別2時對應的車速區間為318～340 km/h；待測工況18、工況20、工況22為類別3時對應的車速區間為340～380 km/h。

4 結束語

本文采用時間序列AR系數聚類的方法，實現了對高速列車通過橋梁時車速的識別。通過計算加速度時間序列AR系數并進行模糊聚類，分析了不同車速列車導致橋梁產生的加速度響應的差異性，進而識別了列車車速。

以列車過橋時的橋梁加速度響應時間序列為研究對象，對不同車速下的加速度時間序列建立AR模型，確定了AR模型的階數，對各個工況下的AR系數進行模糊聚類分析，將已知車速劃分為3個聚類區間，并以此為參考識別了未知車速所處的區間。