基于高斯混合模型的列車隨機振動加速度響應最大值分布統計分析

蒲珍華，吳夢雪，唐德發，朱金，李永樂

(1.西南石油大學 土木工程與測繪學院，成都 610500；2.西南交通大學 土木工程學院，成都 610031)

1.1 列車多體動力學模型

表1 列車模型主要參數表

圖1 列車多體動力學模型

1.2 軌道不平順數值模擬

基于中國高速鐵路總體技術條件，建議對列車進行平穩性分析時使用德國高速線路軌道譜[14]，且高速鐵路試驗段軌道譜的高低不平順在30 m波長以上的平順性基本與德國高速低干擾譜接近[15]。因此，選用德國高速低干擾譜，采用三角級數法對軌道不平順序列進行模擬。在假設軌道不平順為平穩遍歷的高斯白噪聲的前提下，軌道不平順可看作是不同三角級數的疊加，可通過式(1)得到。

(1)

式中：f(x)為模擬得到的軌道不平順序列；S(ωk)為功率譜密度函數，垂向不平順單位為m2/rad/m，水平不平順單位為1/rad/m；ωk為軌道不平順的空間頻率，rad/m；φk為第k個頻率的相位，是獨立均布于0～2π的隨機數。

用三角級數法模擬得到的軌道不平順功率譜密度(Power spectral density，PSD)與目標譜的吻合情況如圖2所示，由圖2可知，模擬的功率譜與目標譜吻合較好。

圖2 軌道垂向不平順模擬譜與目標譜的對比

1.3 列車加速度響應樣本

以隨機軌道不平順作為輸入激勵，得到列車以200 km/h速度行駛時的加速度響應時程曲線，如圖3所示。列車加速度測點位置根據《鐵道車輛動力學性能評定和實驗鑒定規范》(GB 5599—1985)規定設定于轉向架中心上方橫向1 m的車體地板上。

圖3 軌道不平順激勵作用下的列車加速度響應

2 高斯混合模型

高斯混合模型廣泛應用于統計分析領域[16]，其作為一種統計模型，多用于構建概率密度函數，能較好地刻畫參數空間中數據的分布及其特征，既具有非參數化方法的靈活性，又保持了參數化方法的精確性。高斯混合模型采用有限個特定概率分布密度函數的加權組合來擬合復雜的概率分布模型，通過選擇混合分量的類型和個數，可逼近任何連續的概率分布密度函數。高斯混合模型(GMM)由高斯(正態)分布的加權組合得到，其概率密度函數為

(2)

高斯混合模型可以逼近任何連續的概率分布函數，選擇合適的權重系數是得出模型類型數量與模型參數的關鍵。高斯混合模型是一種“軟分類聚類”，是基于假設數據集是由一個潛在的混合概率分布產生的，其中每個高斯分量都表示一個不同的聚類。首先通過估計樣本數據集的混合概率密度，然后計算樣本源中單個樣本之于各個高斯分量的后驗概率，最后將單個樣本分配到后驗概率最大的高斯分量所在的聚類組，從而得到樣本數據集所服從的高斯混合分布[18]。

2.1 參數估計

高斯混合模型(GMM)的期望為

(3)

式中：πk為隨機變量x取至第k個高斯成分的權重系數；N(xi;μk;∑k)為第k個類別的條件概率密度；μk、∑k分別為該高斯分量的均值和協方差矩陣。

對于高斯混合模型的參數，可用期望最大化(EM)算法進行迭代估計[19]。估計步驟為

1)E步

(4)

2)M步

(5)

(6)

(7)

3)收斂條件

不斷迭代E步與M步，至似然函數的變化量小于誤差值esp或迭代次數k≤K時，迭代結束，否則返回E步。隨著迭代次數的增加，迭代誤差越來越小，似然函數不斷收斂。可接受的迭代誤差esp=2×10-16，最大迭代次數K=500。似然函數為

(8)

綜上所述，EM是一種迭代算法，也是一種聚類算法，它可以通過迭代求出高斯混合模型的參數，并將樣本源中的單個樣本通過迭代收斂性進行聚類。高斯混合模型聚類通常采用貝葉斯信息準則(BIC)選擇模型，模型的BIC值越大，該模型就越符合實際。

2.2 擬合度檢驗

假設列車振動加速度響應的最大值服從高斯混合分布，采用期望最大化算法對該概率模型參數進行最大似然估計，再進行擬合度檢驗。

為了對比得到的高斯混合模型概率密度函數與由列車加速度響應最大值的樣本序列得到的頻率直方圖的擬合效果，采用擬合優度(AdjustedR2)與均方根誤差(Root Mean Squared Error，RMSE)兩項指標來對概率密度分布曲線的擬合優劣程度進行評價。擬合優度用于評價概率密度分布曲線與直方圖之間的相似程度，該值越接近于1，表示擬合程度越好；均方根誤差用于評價概率密度分布曲線與直方圖之間的偏離程度，該值越接近于0，表示偏離程度越小，擬合程度越好。

2.3 樣本量

樣本量是指從總體中抽取的樣本元素的總個數，樣本量的大小是選擇檢驗統計量的一個重要要素。由抽樣分布理論可知，在大樣本條件下，如果總體為正態分布，則樣本統計量服從正態分布；如果總體為非正態分布，則樣本統計量漸近地服從正態分布[20]。

樣本量的計算公式為

(9)

式中：n為樣本量；α為顯著水平；Zα/2為置信區間對應的標準分數；E為估計誤差，一般小于0.1；σ為標準差，一般為0.5。

在確定樣本量時，取α為0.05，則置信度為95%，經查表，Zα/2為1.96；假定的估計誤差為0.05，則最小樣本量為n=384。在滿足樣本最小容量的情況下，增加樣本量有助于增加檢驗統計的精度，提高可靠性。

3 列車加速度響應最大值的概率模型

基于Monte-Carlo方法，分別取400、500、600、700、800、900和1 000個隨機生成的軌道不平順樣本，計算得到列車的加速度響應時程樣本，并將列車豎向和橫向加速度響應的最大值作為隨機變量，采用EM算法進行聚類分析，比較聚類為1～4類的BIC值，并選擇BIC值最大的一組參數，得出某個確定樣本容量下列車加速度響應最大值所服從的高斯混合模型。接著，對比不同樣本容量下列車加速度響應最大值的高斯混合模型與其相應頻率分布直方圖的擬合優度及均方根誤差，從中選取合適的高斯混合模型概率密度函數(對應的樣本容量記為Nrep)作為代表該車速下列車加速度響應最大值的概率密度模型，此時的樣本容量Nrep作為代表該車速下列車加速度響應最大值概率統計特征的最優樣本數量。

3.1 列車豎向加速度

當列車以200 km/h車速行駛時，以1 000個列車豎向加速度響應樣本為例，表2為列車豎向加速度最大值不同聚類個數的BIC值。由表2可得，當聚類個數N=2時，BIC值最大。表3為N=2時模型的參數估計值。將表3中的參數帶入到式(2)中，即可得到列車豎向加速度最大值所服從的高斯混合模型概率密度函數，如式(10)所示。圖4為得到的高斯混合模型概率密度函數與相應樣本數量下列車豎向加速度響應最大值的頻率分布直方圖的對比，從圖4可知，其擬合效果較好。

圖4 列車豎向加速度最大值的概率密度分布

表2 列車豎向加速度最大值不同聚類數的BIC值

表3 N=2時模型參數估計值

(10)

分別取400、500、600、700、800、900和1 000個隨機生成的軌道不平順樣本，計算得到車速為200 km/h時不同樣本容量下的列車豎向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度及均方根誤差，如表4和圖5所示。由表4和圖5可知，當樣本數量為400～700時，列車豎向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度隨著樣本數量的增加逐漸變大，在樣本數量達到700后趨于穩定；而均方根誤差隨著樣本數量的增加逐漸減小，在樣本數量達到700后，波動較小，趨于穩定。因而，車速為200 km/h時，代表列車豎向加速度響應最大值概率統計特征的最優樣本數量Nrep=700。圖6為不同樣本數量下的列車豎向加速度最大值所服從的高斯混合模型的概率密度曲線。

表4 列車豎向加速度最大值概率分布的擬合指標

圖5 不同樣本數量下列車豎向加速度最大值概率密度曲線擬合效果評判指標對比

圖6 不同樣本數量的列車豎向加速度最大值概率密度曲線

當車速分別為100、150、200 km/h時，不同樣本數量下列車豎向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度及均方根誤差如圖7所示。從圖7中可以看出，在樣本數量達到700后，3種不同車速下的列車豎向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度和均方根誤差均波動較小，趨于穩定。因而當車速分別為100、150、200 km/h時，代表列車豎向加速度響應最大值概率統計特征的最優樣本數量Nrep均可取為700。

圖7 不同車速和樣本數量下列車豎向加速度最大值概率密度曲線擬合效果評判指標對比

3.2 列車橫向加速度

當列車以200 km/h車速行駛時，仍以1 000個列車橫向加速度響應樣本為例，表5為列車橫向加速度最大值不同聚類個數的BIC值，由表5可得，當聚類個數N=3時，BIC值最大。表6為N=3時模型的參數估計值。將表6中的參數帶入式(2)中即可得到列車橫向加速度最大值所服從的高斯混合模型的概率密度函數，如式(11)所示。圖8為得到的高斯混合模型的概率密度函數與相應樣本數量下列車橫向加速度響應最大值的頻率分布直方圖的對比，由圖8可知其擬合效果較好。

圖8 列車橫向加速度最大值的概率密度分布

表5 列車橫向加速度最大值不同聚類數的BIC值

表6 N=3時模型參數估計值

(11)

分別取400、500、600、700、800、900和1 000個隨機生成的軌道不平順樣本，計算得到車速為200 km/h時不同樣本容量下列車橫向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度及均方根誤差，如表7和圖9所示。由表7和圖9可知，列車橫向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度和均方根誤差在樣本數量達到800后波動較小，趨于穩定。因而，車速為200 km/h時，代表列車橫向加速度響應最大值概率統計特征的最優樣本數量Nrep=800。圖10為不同樣本數量下列車橫向加速度最大值所服從的高斯混合模型的概率密度曲線。

圖9 不同樣本數量下列車橫向加速度最大值概率密度曲線擬合效果評判指標對比

圖10 不同樣本數量的列車橫向加速度最大值概率密度曲線

表7 列車橫向加速度最大值概率分布的擬合指標

當車速分別為100、150、200 km/h時，不同樣本數量下列車橫向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度及均方根誤差如圖11所示。從圖11中可以看出，在樣本數量達到800后，3種不同車速下列車橫向加速度響應最大值的高斯混合模型概率密度分布與其頻率分布直方圖的擬合優度和均方根誤差均波動較小，趨于穩定。因而當車速分別為100、150、200 km/h時，代表列車橫向加速度響應最大值概率統計特征的最優樣本數量Nrep均可取為800。

圖11 不同車速和樣本數量下列車橫向加速度最大值概率密度曲線擬合效果評判指標對比

4 不同車速下列車加速度響應最大值的統計規律

列車分別以車速100、150、200 km/h行駛1 000 m，基于Monte-Carlo方法，取1 000個隨機生成的軌道不平順樣本計算得到列車的加速度響應時程樣本。通過統計列車的豎向加速度最大值及橫向加速度最大值，得到列車不同車速下的加速度最大值箱型圖，如圖12所示。由圖12可知，列車的豎向加速度最大值和橫向加速度最大值均隨著車速的增加逐漸變得離散；列車的豎向加速度最大值和橫向加速度最大值的均值隨著車速的增加而增加。此外，列車豎向加速度最大值的均值都大于中位數，而列車橫向加速度最大值的均值在車速為100、150 km/h時大于中位數，在車速為200 km/h時的均值小于中位數。

圖12 不同車速下的列車加速度最大值箱型圖

通過對列車的豎向加速度最大值樣本及橫向加速度最大值樣本進行分析，進一步得到不同車速下列車加速度最大值的概率密度分布圖，如圖13所示。從圖13中可知，列車豎向加速度最大值和橫向加速度最大值的概率密度曲線均沿著橫坐標向右移動。另外，與車速為100、150 km/h相比，車速為200 km/h時列車豎向加速度與橫向加速度最大值的概率密度曲線的分布范圍均更大；而車速為100 km/h時列車豎向加速度與橫向加速度最大值的概率密度曲線最為高聳。這表明隨著車速的增加，列車加速度響應最大值分布的離散性增強。

圖13 不同運行車速下列車加速度最大值的概率密度曲線

由于列車振動會使乘車人員產生不適感或疲勞，因而引入平穩性指標來度量乘客乘車時的舒適程度。參考高速鐵路客車動力學性能評定資料，中國車體振動加速度的平穩性標準界限值可取為：豎向振動加速度av≤1.3 m/s2，橫向振動加速度ah≤1.0 m/s2。

在運行距離為1 000 m的情況下，選擇列車運行車速為100、150、200 km/h的工況，得到列車加速度最大值隨機變量的累積分布函數(Cumulative distribution function, CDF)曲線，如圖14～圖16所示。從圖中可以發現，列車在運行車速為100、150、200 km/h時的豎向加速度與橫向加速度均滿足平穩性要求。

圖14 100 km/h的列車加速度最大值的累計分布函數曲線

圖15 150 km/h的列車加速度最大值的累計分布函數曲線

圖16 200 km/h的列車加速度最大值的累計分布函數曲線

5 結論

1)采用高斯混合模型能夠有效地擬合列車加速度響應最大值的分布規律。

2)當車速分別為100、150、200 km/h時，列車豎向加速度響應的樣本數量達到700時可較好地代表列車豎向加速度響應最大值的概率統計特征；而列車橫向加速度響應的樣本數量達到800時能較好地代表列車橫向加速度響應最大值的概率統計特征。

3)列車的豎向加速度最大值和橫向加速度最大值均在車速為100 km/h時分布最為集中；整體來講，隨著車速的增加，列車加速度響應最大值分布的離散性增強。