基于子結構的軌道系統基礎結構參數識別方法

葉 玲，陳華鵬，劉昌雨

(華東交通大學 土木建筑學院， 江西 南昌 330013)

隨著高鐵運輸的迅猛發展，軌道運輸在整個交通運輸系統中占有越來越重要的位置。開展車輛荷載作用下超長線狀結構參數識別方法研究，識別結構的損傷狀況，實現車輛參數和結構參數的實時監控和智能化管理，對保障城市交通正常運作至關重要。近年來，越來越多的學者開始關注由于車輛荷載引起的地鐵軌道和橋梁這一類超長線狀結構的損傷識別。Zhou等[1]提出了一種隧道結構損傷識別方法，確定了隧道結構楊氏模量與耦合共振頻率的關系。Zhou等[2]提出了一種基于波傳播的結構健康評估方法，以確定隧道結構的整體剛度，進一步評價隧道結構的服役性能。Liu等[3]提出了一種基于移動荷載作用下大跨度軌道橋梁位移響應靈敏度分析的損傷識別方法。Chen等[4]提出了一種基于有限測點響應同時識別結構損傷和力的方法。Ling等[5]提出了一種未知移動力作用下結構損傷識別方法，該方法以結構單元剛度變化為損傷指標，以動力響應靈敏度為優化方向。文獻[6-8]通過對隨機荷載作用下系統動力響應靈敏度特點的研究，提出了一種用于識別大型結構的損傷識別方法，并將其應用到大跨度軌道橋梁結構中，識別了車輛荷載作用下的橋梁結構損傷。Majumder等[9]將車輛系統簡化為單自由度移動質量塊，研究了基于車輛動力響應和橋梁結構動力響應對橋梁結構進行損傷識別的方法。

上述方法在對地鐵軌道和橋梁這一類超長線狀結構進行有限元建模時，通常只能選取一定長度的結構作為研究對象，這種選取一定長度模擬無限長軌道的方法，缺少物理意義，且在識別結構損傷時，也容易因為邊界問題引起識別誤差，同時，由于這類結構本身結構規模大、待識別參數多，極大地限制了識別效率。子結構參數識別法將子結構法和參數識別方法相結合，將大型結構劃分為若干個子結構，在保證相同條件下，有效減少研究對象的單元數目，提高算法的精度和效率，為超長線狀結構的參數識別提供了新思路。Weng等[10]在力和位移兼容性的約束下，分解整體結構的模態數據，提出了一種新的逆向結構損傷識別方法，可以高精高效的識別結構損傷。隨后Koh等[11]采用“分而治之”的戰略思想解決了當子結構和整體結構差別較大時，很難在數值上獲得合理準確結果的難題，給出了漸進結構識別方法。Yun等[12]采用子結構識別法和子矩陣縮放因子，克服了未知參數相關性等問題。Koh等[13]提出了一種不需要界面測量信息的子結構損傷識別方法。Law等[14]提出了一種從支撐激勵下結構的加速度響應中識別結構間耦合力的方法，并基于動態響應靈敏度分析對所識別的耦合力進行了局部結構損傷檢測，該方法在狀態空間域中建立了耦合力的識別方法，并用阻尼最小二乘法求解。

為了高精高效的識別軌道系統基礎結構參數，本文提出了一種基于子結構的超長線狀結構參數識別方法。在車輛荷載作用下，將子結構界面力作為待識別子結構的未知外荷載，表達為切比雪夫多項式，推導系統動力響應關于結構參數和界面力參數的靈敏度矩陣，根據車軌系統參數同步識別方法[15]，同時識別軌道子結構參數和界面力參數。通過一個車輛-軌道交互系統數值算例驗證了基于子結構的軌道系統基礎結構參數識別方法的精度和效率。

1 獨立子結構的動力響應分析

1.1 車輛-軌道交互系統的運動方程

車輛系統離散為附有二系彈簧的多剛體整車模型。假設車輛的各個組成部件都關于各自質心前后、左右對稱，取沿縱向線路對稱的半車結構進行研究。定義位移和力的方向以豎直向下為正，轉角和力矩的方向以逆時針方向為正，將每個車輛看作一個單元，則半車模型共有10個自由度，見圖1，車輛系統位移響應矩陣Xa可記為

(1)

式中：xc和φc分別為車體豎向位移和轉角；xeu和φeu(u=1,2)分別為車輛第u個轉向架豎向位移和轉角；xwi(i=1～4)為車輛第i個車輪豎向位移。

圖1 車輛模型

由Hamilton原理[16]，可得車輛系統運動方程

(2)

列車的質量、阻尼和剛度矩陣分別可表示為

(3)

(4)

(5)

式(3)～式(5)中：Mc和Me分別為列車車體質量和轉向架質量；Jc和Je分別為列車車體點頭慣量和轉向架點頭慣量；Cs1和Cs2分別為一系懸掛裝置和二系懸掛裝置的阻尼；Ks1和Ks2分別為一系懸掛裝置和二系懸掛裝置的剛度；L1為單個轉向架上兩車輪中心之間的半距；L2為車體上兩個轉向架中心之間的半距。

本論文僅考慮車輛荷載對超長線狀結構的豎向作用，將隧道-土體簡化為彈性塊式支承，用彈性支承塊式無砟軌道模型模擬軌道結構，彈性支承塊式無砟軌道模型由鋼軌、混凝土道床板、隔離層及混凝土底座等組成，見圖2。由于軌道結構沿縱向線路橫向對稱，故取對半軌道結構研究。鋼軌用Bernoulli-Euler梁結構模擬；軌枕簡化為僅在豎向發生振動的質量塊Ms，且沿軌道長度方向等間距布設；鋼軌與軌枕之間，軌枕與道床之間均通過彈簧-阻尼器連接，兩組彈簧-阻尼器剛度系數和阻尼系數分別記為kp,cp和kb,cb。

圖2 彈性支承塊式無砟軌道模型

因軌道模型具有周期性結構特性，故可根據軌枕所在位置劃分鋼軌結構，形成軌道單元，則軌道系統轉化為由多個軌道單元首尾相連而成的結構。其中軌道單元見圖3。彈性支承塊式無砟軌道的振動主要體現在鋼軌和混凝土支承塊上。鋼軌只考慮豎向位移和轉動，軌枕只考慮豎向位移，定義位移和力的方向以豎直向下為正，則一個軌道單元包含5個自由度。

圖3 軌道單元模型

軌道單元的質量矩陣可看作鋼軌質量矩陣與軌枕質量矩陣的疊加，具體表達式為

Mcell=Mr+Ms

(6)

式中：

(7)

(8)

其中，mr為單位長度(1 m)鋼軌質量；ms為軌枕的質量；l為一個軌道單元中鋼軌的長度。

同理，軌道單元的剛度矩陣表達式為

Kcell=Kr+Ks

(9)

式中：

(10)

(11)

其中，I為鋼軌轉動慣量；E為鋼軌彈性模量；kp和kb分別為鋼軌墊片和道床的剛度。

軌道單元的阻尼矩陣表達式為

Ccell=Cr+Cs

(12)

式中：

(13)

式中：cp和cb分別為鋼軌墊片和道床的阻尼。

鋼軌阻尼矩陣采用瑞利阻尼矩陣，可以表示為鋼軌質量矩陣和剛度矩陣的線性組合，即Cr=αMr+βKr，α和β分別為瑞利阻尼系數，根據結構動力學原理，α、β分別為

(14)

其中，η1和η2分別為鋼軌的一階和二階頻率；ξ1和ξ2分別為鋼軌相關的一階和二階阻尼比。

軌道系統總的質量矩陣、阻尼矩陣和剛度矩陣可以通過組集軌道單元的質量矩陣、阻尼矩陣和剛度矩陣得到，軌道系統的運動方程最終可表示為

(15)

Nbj,i=

(16)

式中：xh,i(t)為t時刻第i個車輪經過鋼軌的距離；j為t時刻第i個車輪已通過的鋼軌單元個數。

t時刻車輪通過的距離xh,i(t)可以表示為

(17)

式中：x0為起始時刻第4個車輪所在位置到鋼軌最左端的距離；V為車輛行駛速度，本文車速定義為恒速。

本文采用余弦函數表示鋼軌焊接接頭和鋼軌波磨的形狀，即軌道位移不平順表達為

(18)

式中：μ和分別為軌道位移不平順的統計波深和波長，根據軌道線路實地測量獲得；L為實際的軌道位移不平順長度。

車輛與軌道之間采用非線性Hertz接觸彈簧模擬輪軌法向接觸力，則輪軌法向接觸力[17]可以表達為

(19)

式中：G為輪軌接觸常數,m/N2/3；δX為輪軌間的彈性壓縮量，其值根據車輪接觸面形狀有不同取值，m。

對于磨耗型踏面車輪計算式為

G=3.86R-0.115×10-8

(20)

對于錐形踏面車輪計算式為

G=4.57R-0.149×10-8

(21)

式中：R為車輪半徑,m。

輪軌間的相對位移由該時刻，車輪的豎向位移和車輪所在鋼軌接觸點處鋼軌的豎向位移確定，即

δX=xw,i-xr,ii=1～4

(22)

式中：xw,i為第i個車輪的豎向位移,m；xr,i為第i個車輪所在鋼軌處鋼軌的豎向位移,m。

當輪軌接觸界面存在軌道不平順Xirr時，輪軌法向接觸力可以表達為

(23)

式中：x=xw,i-xr,i-xirr,i；xirr,i為第i個車輪所在鋼軌處的軌道位移不平順，m；xirr,i=Xirr(xh,i(t))。

令Kw=G-3/2，則輪軌法向接觸力可以表達為

(24)

(25)

將式(2)、式(15)、式(25)組合到一起，即可得模擬車輛荷載作用下超長線狀結構的車輛-軌道交互系統運動方程，表示為

(26)

由式(3)～式(5)，可獲得車輛系統質量矩陣Ma、剛度矩陣Ka和阻尼矩陣Ca，通過式(6)、式(9)、式(12)，可組集得到軌道系統質量矩陣Mb、剛度矩陣Kb和阻尼矩陣Cb。在給定的時間序列中，可寫出所有指示矩陣D和形函數矩陣N。該運動方程可通過基于Newmark的時間步內交叉迭代算法[15]計算得到。

1.2 子結構方法的應用

應用子結構方法[18]劃分軌道結構，根據是否有車輛經過，可以將軌道結構劃分有車輛荷載作用的軌道子結構和沒有車輛荷載作用的軌道子結構。軌道整體結構劃分為軌道子結構的示意圖見圖4，根據是否有車輛通過，可化分為有車輛荷載作用的軌道子結構1和沒有車輛荷載作用的軌道子結構2。Fi1和Fi2分別為軌道子結構1和軌道子結構2的界面力，車輛以速度V在軌道子結構1上運動。

圖4 軌道整體結構劃分子結構示意

對于軌道子結構1，同時受到車輛荷載和子結構界面力作用，軌道子結構1的運動方程可以表達為

(27)

由式(27)第一個方程式可得

(28)

界面力Fi(1)用切比雪夫正交多項式表示，即

(29)

式中：m為正交多項式的項數；c為界面力正交系數；T為切比雪夫正交基，計算式為

(30)

式中：Dt為車輛荷載的持續時間。

軌道子結構1上同時受到車輛荷載和界面力作用Q(1)，計算式為

Q(1)=Ds(1)Ns(1)fwr(1)-Fi(1)

(31)

將式(31)代入式(28)中，軌道子結構1運動方程可表示為

(32)

結合式(26)、式(32)，即可得到軌道子結構1和車輛系統組成的車輛-軌道交互系統運動方程為

(33)

同理，對于軌道子結構2，其運動方程可以表示為

(34)

式中：上標(2)表示第二個軌道子結構。

由式(34)第一個方程式可得

(35)

式(35)右邊為相鄰軌道子結構對軌道子結構2的界面力，該界面力同樣可以用切比雪夫正交多項式表示，同式(29)，作用在軌道子結構1和軌道子結構2上的界面力為一對大小相等，方向相反的作用力和反作用力，軌道子結構2的運動方程可以表達為

(36)

同理，式(33)和式(36)均可以通過時間步內交叉迭代算法[15]計算。

2 基于子結構的軌道系統基礎結構參數識別方法

軌道結構損傷體現在結構剛度的弱化，故將結構參數定義為軌道系統中鋼軌的單元剛度的變化，對結構參數的識別即為對鋼軌單元剛度的識別。定義為

(37)

(38)

由于軌道剛度矩陣Kb是由軌道結構中鋼軌的單元剛度EI組成的，將式(38)帶入軌道剛度矩陣表達式中，即可得到軌道剛度矩陣關于鋼軌單元剛度相對變化量的一階導?Kb/?λp。

子結構界面力用切比雪夫正交多項式表示，則對界面力的識別即為對正交多項式系數cm的識別，需要識別的正交多項式系數的數量為m。計算式(33)關于結構參數的一階導數，可獲得軌道子結構1下系統動力響應關于結構參數λp的靈敏度矩陣

計算式(33)關于正交系數cm的一階導數，可獲得軌道子結構1下系統動力響應關于界面力正交系數cm的靈敏度矩陣

(40)

計算式(36)關于結構參數的一階導數，可獲得軌道子結構2下系統動力響應關于結構參數λp的靈敏度矩陣

(41)

計算式(36)關于正交系數cm的一階導數，可獲得軌道子結構2下系統動力響應關于界面力正交系數cm的靈敏度矩陣

(42)

式(39)～式(42)可以通過交叉迭代算法計算。

本章采用基于靈敏度分析的有限元模型修正方法，將通過子結構模型計算得到的系統加速度響應與整體結構同條件下測量得到的系統加速度響應的殘差作為目標函數J，可以表示為

(43)

參數識別過程中需要修正的參數γ(1)由結構參數λ以及界面力正交系數c組成，其中，λ=[λ1λ2…λp]T，c=[c1c2…cm]T，γ(1)=[λ1λ2…λpcp+1cp+2…cp+m]T。基于子結構法的參數識別方程可以表示為

(44)

基于子結構的軌道系統基礎結構參數識別方法的具體流程見圖5。

圖5 基于子結構的軌道系統基礎結構參數識別流程圖

基于子結構的軌道系統基礎結構參數識別方法具體步驟如下：

Step2當迭代次數k=1時，利用式(39)、式(41)計算系統加速度響應關于待識別結構參數λ的靈敏度矩陣Sλ；利用式(40)、式(42)計算系統加速度響應關于界面力參數c的靈敏度矩陣Sc。

Step3根據參數同步識別方法，利用式(44)計算出當前迭代步的待識別結構參數λk和ck。

Step4判斷待識別參數λk和ck是否滿足優化過程的目標函數收斂準則，即|J|≤Tolerance，滿足條件，則迭代停止；若不滿足，令k=k+1，重復Step1—Step4，直到目標函數滿足收斂條件后，迭代停止，得到參數最終的識別結果λk和ck。

3 車輛-軌道交互系統參數識別數值算例

本文用于模擬軌道交通系統的車輛模型和軌道模型分別見圖1和圖2。車輛在軌道上運行距離為98.1 m，軌道單元長度l=0.545 m，經過180個軌道單元，車輛從軌道左側勻速通過軌道，車速為恒速V=200 km/h，考慮波長和波深分別為3.3 cm、25 μm的軌道不平順，時間步長取為1×10-4s，時程曲線見圖6。軌道各參數見表1，車輛各參數見表2。以車輛荷載下的軌道結構為研究對象，將軌道結構在節點91處劃分成兩個軌道子結構，軌道子結構有限元模型見圖4。選取時間歷程前0.5 s的測量動力響應作為子結構結構參數識別的測量響應，前0.5 s內車輛僅在軌道子結構1上運行。采用基于子結構的參數識別方法計算結構參數和界面力參數。

圖6 軌道不平順時程曲線

表1 軌道參數

表2 車輛參數

車輛荷載作用下軌道子結構單元剛度發生變化的各種工況見表3，工況1為軌道結構單元81發生20%的單元剛度折減，此時發生單元剛度變化的單元處于軌道子結構1上，故可僅以軌道子結構1為研究對象，同步識別軌道子結構1未知界面力作用下的車軌交互系統參數和界面力參數；工況2為在工況1的基礎上，對測量響應加入5%噪聲的情況；工況3為軌道單元81和單元122都出現單元剛度折減，變化程度分別為20%和30%的情況，單元81位于軌道子結構1，單元122位于軌道子結構2中，同時選取兩個軌道子結構作為研究對象，同步識別結構參數和界面力參數；工況4為在工況3的基礎上，測量響應加入5%噪聲的情況。

表3 車輛荷載下子結構參數識別的各種工況

本文測量噪聲用一組正態隨機分布數據進行模擬，加了噪聲之后的測量加速度響應可表示為

(45)

選取工況1和2為研究對象，即以軌道子結構1為目標子結構，界面力Fi1和Fi2作為軌道子結構1的等效輸入力，識別軌道子結構1的結構參數和界面力參數，其中，界面力表達為切比雪夫多項式，階數取為30。子結構參數識別方法的目標函數為通過子結構模型計算得到的系統加速度響應與整體結構同條件下測量得到的系統加速度響應的殘差，故整體結構條件下測量得到的系統加速度響應為已知條件，理論界面力可通過將整體軌道動力響應帶入式(32)計算得到。無測量噪聲(工況1)和5%測量噪聲(工況2)下軌道子結構1的界面力時程曲線分別見圖7、圖8。工況1和工況2中所識別界面力與理論界面力的相對誤差e計算式為

(46)

工況1和工況2中所識別界面力的相對誤差見表4。表4可知，在無測量噪聲條件下，界面力Fi1和Fi2識別結果與相應理論界面力的相對誤差值為3.36%和3.26%，兩者的時程曲線非常吻合；在5%測量噪聲條件下，界面力Fi1和Fi2識別結果與相應理論界面力的相對誤差值為10.68%和10.77%，相對誤差值變大，雖然出現了微小的偏差，但是兩則的時程曲線基本一致。

圖7 工況1理論界面力與識別界面力時程曲線

圖8 工況2理論界面力與識別界面力時程曲線

表4 工況1和工況2中所識別界面力的相對誤差值 %

圖9 結構參數識別結果

不考慮噪聲和考慮5%噪聲條件下，以軌道子結構1為研究對象時，識別的結構參數識別見圖9。由圖9(a)、圖9(b)可知，不考慮噪聲條件時，軌道結構單元剛度發生變化的位置和程度可以準確識別出來；考慮5%測量噪聲時，出現了識別誤差，但是誤差較小，即仍然可以準確的識別出軌道結構單元剛度發生變化的位置和程度。以工況3和工況4為研究對象，即軌道子結構1和軌道子結構2上均存在單元剛度發生變化時，需要將軌道子結構1和軌道子結構2同時為識別對象，將軌道子結構1識別出的界面力作為軌道子結構2的輸入荷載，分析工況3和工況4，識別在車輛荷載作用下的系統參數和界面力參數。無測量噪聲(工況3)和5%測量噪聲(工況4)下軌道結構單元剛度相對變化量的識別結果見圖9(c)～9(f)，軌道子結構1和軌道子結構2中發生單元剛度變化的軌道單元位置和程度均能準確識別出來，考慮測量噪聲條件時，軌道子結構的邊界單元以及發生單元剛度變化的單元附近的單元出現了識別誤差。

為了驗證基于子結構的軌道系統基礎結構參數識別方法的高效性，與通過整體參數識別法計算同一車軌交互系統模型參數的結果對比。整體參數識別法和子結構參數識別法計算車軌交互系統模型參數所用的時間見表5。由表5可知，子結構參數識別法迭代一步僅需44.59 s，需迭代8次，所用總時間為356.7 s；整體參數識別法迭代一步需129.78 s，需迭代13次，所用總時間為1.687×103s，約為子結構參數識別法計算所用時間的4.73倍，這是由于軌道整體結構的自由度個數為905個，軌道子結構1的自由度個數只有455個，在識別過程中每一次迭代的動力響應和動力響應靈敏度矩陣的維數由整體參數識別法的905×(5×103)減少至子結構參數識別法455×(5×103)，同時所需識別的參數數量也由整體參數識別法的180個減少至子結構參數識別法的120個。與整體參數識別法對比，子結構參數識別法在計算過程中包含更少的自由度個數和待識別參數個數，大幅減少了計算時間，有效提高了計算效率。

表5 整體參數識別法和子結構參數識別法計算效率比較

4 結論

針對軌道系統基礎結構參數識別效率極低的難題，本文提出了一種基于子結構的軌道系統基礎結構參數識別方法。從對車輛-軌道交互系統的數值分析結果，可以得到以下結論：

(1)基于子結構的參數識別方法可以同步識別子結構參數和荷載參數。

(2)由于子結構方法減少了計算過程中各矩陣的維度和待識別參數個數，極大地提高了參數識別的效率。