基于雙輸入單輸出ARX模型的塔機起重臂結構損傷診斷方法

王勝春，李文豪，安 宏，安增輝

（山東建筑大學機電工程學院，濟南 250101）

塔機作為建筑領域不可或缺的工程機械，其結構一旦出現問題會對人民的生命及財產安全產生極大的威脅。塔機金屬結構的主要部件為：起重臂、平衡臂、上下支承座和塔身等。塔機的金屬結構是其骨架，承受著塔機在運行過程中的自重和各種外荷載。組成塔機金屬結構的構件較多，其重量通常占整機重量的一半以上，而塔機起重臂長、自重大，因此，展開起重臂健康狀態的研究對塔機的安全運行尤為重要[1]。

對塔機起重臂健康狀態的實時診斷和定期維護，是保證塔機結構安全、正常運行的前提，而準確診斷起重臂結構健康狀態則是對塔機定期維護的基礎。對于結構損傷檢測的方法有很多種，其中基于時頻響應的識別方法需要在時域或頻域處理結構動力響應，而時間序列分析作為時域數據處理的重要工具，已被許多學者應用于結構損傷識別領域。Jayawardhana 等[2]提出了一種通過以Fisher 準則把用AR 模型的系數進行對比的損傷識別方法。Mosavi等[3]利用向量自回歸滑動平均模型求損傷指標，采用馬氏距離和Fisher準則為標準對連續鋼梁進行損傷識別，驗證了該方法的可能性。刁延松等[4]利用計量經濟學中的協整理論對時間序列的自回歸模型系數進行處理，提取損傷指數，識別結構損傷。Golnaz Shahidi等[5]用ARX模型對一個雙機架鋼管框架進行損傷檢測，結果表明該方法成功識別了損傷的發生，且對于損傷位置有精確的定位。Gul 等[6-7]提出基于自回歸滑動平均模型（Auto-regressive Moving Average Model，ARMA）的統計模式識別方法，提高了時間序列方法在損傷檢測過程的魯棒性。

上述方法在結構損傷識別中均取得一定效果，但存在一些不足。例如Gul等的研究僅對結構的特定損傷程度進行驗證，難以保證不同損傷程度的識別效果；基于時間序列模型結構損傷診斷方法的核心是時間序列建模，時間序列模型的精度對診斷效果至關重要，而模型參數求解的精確性是建立模型的關鍵一步，因此要盡可能減小模型參數的偏差。基于此，本文針對塔機起重臂結構的健康診斷，建立一種基于時序數據的雙輸入單輸出ARX 的起重臂結構損傷診斷模型，模型的輸入和輸出數據是臂架各節點的動態位移數據。接著求解模型參數的最小二乘估計，并且為進一步減小最小二乘法參數估計產生的偏差，提出利用花粉算法對其進行優化。然后定義模型的損傷特征參數和損傷程度參數，以完整狀態的動態響應作為狀態判斷的參考，最后分析損傷狀態響應信號與完整狀態響應信號的差異，從而診斷起重臂結構的健康狀態、損傷位置及損傷程度。

1 基于雙輸入單輸出ARX 模型的方法

1.1 雙輸入單輸出的ARX模型表達式

時間序列數據是通過數據的大小和相應的數據順序來包含信息的數據。時間序列動態數據建模可以很好地研究包含在動態數據中的特征信息，利用特征信息進行損傷診斷。根據動態時序數據建立一種雙輸入單輸出模型，其基本形式如式（1）所示。

式中：b—輸出函數的系數；a、c—輸入函數的系數；en—模型殘差。

BⅠC值的定義如式（2）所示。

式中：p—模型階次；N—數據段長度；σ2—殘差的方差。

定階圖如圖1 所示。由圖可知BⅠC 值從第5 階趨于穩定，由于當模型階次取得過大時會增加模型復雜性，且復雜模型帶來的巨大計算量會導致模型誤差的增大，因此根據系統的物理性質和BⅠC 準則定義模型的階次為5。

圖1 BⅠC定階圖

參數矩陣φ的最小二乘估計如式（3）所示。

其中：Ym、Yn如式（4）至式（5）所示。

式中：φ—模型系數矩陣；E—殘差矩陣。

為了進一步減小最小二乘法的參數估計產生的偏差，本文采用花朵授粉算法對模型參數進行優化[8]。所優化的目標函數如式（8）所示。

式中：N—數據段長度；p—模型階次；φ—模型系數。

花朵授粉算法模仿自然授粉的兩種機制。對于異花授粉由于其依靠傳粉者遠距離傳播花粉，將其對應為全局搜索過程；對于自花授粉，由于其授粉物理位置上距離較近，將其對應為局部搜索過程。引入切換概率p（p∈[0,1 ]），以此來權衡局部搜索過程與全局搜索過程之間的比重。

花朵授粉算法的實現步驟如下：

（1）參數初始化，包括花朵種群數n=25，轉換概率p=0.8，迭代次數Niter=4 000。

（2）根據目標函數計算各解的適應度值，求出當前最優解。

（3）如果轉換概率p＜rand 條件成立，執行異花授粉，對應全局搜索，按式（9）對當前解進行更新。

其中：Xt+1i，Xti分別是第t+1 代、第t代的解；g*是全局最優解；L是步長，L的計算公式參照文獻[9]。

如果轉換概率p＞rand 條件成立，執行自花授粉，對應局部搜索，按式（10）對當前解進行更新。

其中：?∈[0,1 ]上服從均勻分布的隨機數，Xjt、Xkt是從所求得解的內部隨機選取的異于Xit的解。

（4）計算式（3）得到的新解對應的適應度值，若新解的適應度值優，則用新解和新解對應的適應度值分別替換當前解和當前適應度值，否則保留當前解和當前適應度值。

（5）如果新解的適應值優于全局最優值，則更新全局最優解和全局最優值。

（6）判斷結束條件，滿足條件后退出程序，并輸出最優值和最優解。

1.2 損傷特征提取

（1）首先提取結構完好狀態相關節點的動態位移數據，將{ya1,ya2,…,yan}、{yc1,yc2,…,ycn}作為輸入數據，將{yb1,yb2,…,ybn}作為輸出數據代入式（1）中建立基于時域數據的雙輸入單輸出的ARX模型，此時e（n）為白噪聲序列。

（2）利用最小二乘法求出模型參數a、b和c，為減小最小二乘法對模型參數求解產生的偏差，進一步采用花朵授粉算法對最小二乘法求得的參數進行優化，進而建立基于時域數據的雙輸入單輸出的損傷診斷模型。

（3）將待檢測數據代入所建立的損傷診斷模型中，計算模型殘差的方差，如式（11）所示。

式中：σ2—殘差的方差，N—數據段長度，p—輸入階次。

（4）將計算所得的模型殘差的方差作為損傷特征參數來診斷結構是否出現損傷，當塔機起重臂結構出現損傷時表現為結構完好狀態的損傷特征參數與結構損傷狀態的損傷特征參數出現偏差[10]。

（5）定義損傷程度參數Dc來評價損傷工況的損傷特征參數與完好工況的損傷特征參數之間的偏差程度，其表達形式如式（12）所示。

式中：σ2a—完好工況殘差的方差；σ2b—損傷工況殘差的方差。

1.3 損傷診斷流程

（1）獲取塔機起重臂相關節點完好工況的動態位移數據和待檢測工況的動態位移數據；

（2）利用完好工況的動態位移數據建立基于雙輸入單輸出的ARX模型；

（3）定義損傷特征參數和損傷程度參數；

（4）將完好工況的動態位移數據和待檢測工況的動態位移數據分別輸入模型，求得對應的損傷特征參數和損傷程度參數；

（5）通過損傷特征參數和損傷程度參數對塔機起重臂健康狀態進行診斷。

2 有限元建模及動態信號提取

為驗證基于時域數據的雙輸入單輸出ARX 起重臂結構診斷模型的有效性，本文針對QTZ80塔機進行有限元建模及分析。

2.1 有限元模型的建立

為盡可能地保證有限元建模的準確性及對建模過程的簡化，本文在ANSYS 中參照文獻[11]對塔機進行建模。定義彈性模量E=210 GPa，泊松比u=0.3，密度ρ=7 800 kg/m2，建立模型如圖2 所示。本文以塔機回轉中心為坐標原點，將沿起重臂方向定義為X方向，將沿塔身方向定義為Y方向，將垂直于塔身和起重臂所在平面的方向定義為Z方向。

圖2 塔機有限元模型

2.2 動態信號的提取

本文選用的QTZ80 塔機起重臂的最大工作幅度為60 m，起重臂最大工作幅度處額定起重量1.1 t，最大起重量8 t，起升高度45 m，起升速度20 m/min。該塔機在起吊重物時，3.7 s 末重物離地，6 s 末達到額定起升速度。參考文獻[11]對塔機起升激勵動載荷進行求解，將求得的起升激勵載荷施加在塔機有限元模型中起重臂的最大工作幅度60 m處，采用完全法對重物離地過程進行瞬態分析，提取塔機起重臂上相關節點的動態位移數據。

2.3 損傷工況

本文所選用塔機的起重臂共七節，從起重臂最遠端至臂根處依次為臂節一至臂節七。針對塔機起重臂的臂節連接處銷軸松動這一損傷特征進行研究，在模型分析中可認為結構損傷只引起損傷單元剛度的減小，用彈性模量E的降低來模擬單元損傷[12]。設定臂節四與臂節五連接段的542 單元損傷位置如圖2 所示，損傷單元長度占起重臂長度的0.1%。本次研究選擇的損傷工況如表1所示。

表1 損傷工況表

經研究發現，當損傷程度為20%時，由于損傷單元長度小且損傷程度輕，該損傷對塔機起重臂宏觀位移產生的影響極小（位移平均變化量為4.5×10-5m），診斷效果不理想，因此本文設定損傷程度從30%開始。表中數據含義例如：臂節四與臂節五連接段的542 單元損傷90%（即該單元彈性模量減少90%）記為工況一。

采集起重臂臂節上節點的動態位移數據，采樣的時間間隔Δ=0.025 s，采樣點如表2所示。

表2 數據采樣表

因篇幅有限，本文僅列出Z方向的動態位移數據，以節點150為例，完好工況和工況四的動態位移數據如圖3所示。

圖3 節點150的Z方向位移對比

由圖3可知，節點150在完好工況與工況四的Z方向的動態位移數據非常接近，僅對動態位移數據進行觀察無法診斷結構的健康狀態，為方便后續的研究，后文所采用的動態位移數據均為Z方向的動態位移數據。

3 仿真結果與分析

針對塔機起升工況的起重臂損傷診斷，從起重臂的臂根開始：將節點150和節點231的動態位移數據作為輸入，將節點213的動態位移數據作為輸出，經過損傷診斷模型計算得到損傷特征參數，記為D1，具體的輸入輸出方式如表3所示。

表3 輸入輸出方式

3.1 起重臂結構健康狀態診斷

應用本文所述方法對工況一至工況四進行診斷，因篇幅限制原因，本文僅列出了損傷程度較輕的工況四的診斷結果，其他工況較工況四相比診斷效果更好。以200 個動態位移數據為一組，按時間順序取12組完好工況位移數據，并以相同的時間起點按時間順序取12組損傷工況位移數據，將每組數據輸入損傷診斷模型，分別得到12組完好工況對應的12 個損傷特征參數和12 組損傷工況對應的12 個損傷特征參數，將所提取的各組的動態位移數據的均值作為橫坐標，損傷特征參數作為縱坐標，對動態位移數據的均值和損傷特征參數分別進行歸一化處理并作圖。以D1為例，模型參數經花朵授粉算法優化前后的起重臂結構健康狀態診斷效果如圖4 所示。由圖4（a）可以看出，在損傷程度較輕時，完好工況的損傷特征參數（圓形符號）和損傷工況的損傷特征參數（三角形符號）不容易區分，因此若僅用最小二乘求解的參數建立診斷模型，會導致診斷效果不理想，而將最小二乘法求解的參數經花朵授粉算法優化以后建立的診斷模型精度更高，如圖4（b）所示，完好工況的損傷特征參數和損傷工況的損傷特征參數區分明顯，對于起重臂結構健康狀態的診斷效果更好。

圖4 花粉優化前后診斷效果對比圖

工況四的完整診斷結果如圖5所示。

圖5的損傷狀態診斷結果表明：

結構完好工況下的損傷特征參數和損傷工況下的損傷特征參數明顯劃分為兩類：結構完好工況下的損傷特征參數均分布在0 附近，損傷工況下的損傷特征參數均分布在1 附近，兩者之間有明顯的差異。因此在各損傷工況下，用本文所述方法均可較為準確地診斷出此時塔機起重臂結構出現損傷。

由圖5 的（a）至圖5（i）可知，利用起重臂上的任意一組相鄰三個節點的動態位移數據經過本文所述的結構損傷診斷模型的處理均可診斷出起重臂結構出現損傷，因此在實際應用中，可用較少的傳感器，例如三個傳感器即可對塔機起重臂的健康狀態進行診斷。

圖5 工況四損傷狀態診斷結果

圖5 工況四損傷狀態診斷結果

3.2 起重臂結構損傷位置及損傷程度診斷

應用上述方法對工況一至工況四的起重臂結構損傷位置及損傷程度進行診斷。應用式（12）求得D1至D9對應的損傷程度參數Dc，將損傷程度參數Dc進行歸一化處理，作損傷位置及損傷程度診斷圖，如圖6所示。

圖6 損傷位置及損傷程度診斷圖

圖6損傷位置及損傷程度診斷結果表明：由圖6中對工況一至工況四的Dc最大值的標記可以看出，工況一至工況四的D4處的Dc均最大，表明D4受結構損傷的影響最大。由表3可知，D4為由節點332和節點455的動態位移數據作輸入，由431節點的動態位移數據作輸出所得的損傷特征參數，預設的損傷部位542 單元由圖2 可知位于節點332 和節點455 之間。因此，可通過本文所述方法將起重臂損傷位置定位在一個確定的范圍，若想獲得損傷位置更精確的定位，可考慮通過在該范圍內布置更密的傳感器來實現。

不同工況下的損傷特征參數分布基本相同，即使損傷程度較小，也能準確地診斷出損傷位置。隨著損傷程度的增加，損傷程度參數Dc不斷增大。因此，可利用該規律在對起重臂結構損傷程度進行診斷時，將某一特定的損傷程度設為基準損傷工況，通過有限元分析計算該損傷程度下對應的各完好工況損傷特征參數與損傷工況損傷特征參數的Dc，并通過對比待測工況與基準損傷工況的Dc的分布情況。對起重臂結構損傷位置及損傷程度進行較為準確的診斷。

4 結語

本文根據時間序列模型建立了基于雙輸入單輸出ARX的起重臂結構損傷診斷模型，并且為減小最小二乘法對模型參數估計產生的偏差，進一步提出用花朵授粉算法對所求參數進行了優化。介紹了損傷特征參數的提取方法，并且定義了損傷程度參數Dc，以此對塔機起重臂結構健康狀態、損傷位置及損傷程度進行診斷。對塔機起升狀態的動態激勵載荷進行求解，并建立塔機有限元模型，通過瞬態動力學求解以獲取塔機起重臂上節點的完好工況及損傷工況的動態位移數據。仿真結果分析表明：

（1）應用本文所述的基于雙輸入單輸出ARX起重臂結構損傷診斷模型可以利用少量傳感器（例如三個傳感器），即可對塔機起重臂結構的健康狀態進行有效的診斷。

（2）在塔機起重臂結構出現損傷時，可用本文所述方法通過布置多個傳感器對起重臂結構損傷的位置及損傷的程度進行較為準確的診斷。

（3）本文所述的基于雙輸入單輸出ARX起重臂結構損傷診斷模型對塔機起重臂健康狀態的實時診斷及損傷預警具有一定的參考價值。