一種基于最優輸出跟蹤的多源動態載荷識別方法

陳 睿, 劉 杰, 張 正, 張 偉

(1. 湖南大學汽車車身先進設計制造國家重點實驗室，湖南 長沙 410082；2. 吉首大學物理與機電工程學院， 湖南 吉首 416000 )

引 言

在機械、結構系統的力學計算、結構設計、故障診斷中，載荷的確定是一個非常重要的問題，它為結構的設計、計算以及分析提供可靠的載荷依據，為減小振動、提高結構的可靠性、安全性，提供確切的環境條件。然而在很多實際工程中，如導彈及航天器在飛行時的載荷、發電站工作時機組及脈動水壓力、橋梁所受運動態載荷等，由于技術條件或經濟條件，對結構所受外載荷難以直接測量甚至根本無法實測。因此，近年來基于結構測取的響應識別結構所受外載荷的技術取得了很大的進步。

載荷識別在結構動力學中屬于第二類反問題，它根據已知系統的動態特性和實測的動力響應反求出結構所受的動態激勵。動態載荷識別的研究最早可追溯到20世紀70年代末，從軍事用途發展起來的[1]。目前動態載荷識別的方法主要有頻域法和時域法兩大類方法[2,3]。頻域法發展較早，是比較成熟的識別方法，其應用較為廣泛。該方法動態標定簡單、識別精度較高，但要求測試數據的樣本具有一定的長度，一般只適用于穩態或平穩隨機載荷識別，對沖擊型瞬態載荷的識別則受到限制。時域法提出較晚，它從系統動力學方程出發，根據響應的時間歷程直接確定動態載荷的時間歷程。時域法直觀，廣泛應用于工程中，其優點在于時域法可對非平穩動態載荷識別，對沖擊型載荷識別更是具有實際工程價值。目前國內外已開展了廣泛研究，并取得了一定的成果[4～6]。

本文從系統動力學方程出發，基于最優輸出跟蹤的基本思想，在載荷位置確定的基礎上，通過實測的動力響應反求出結構所受的動態載荷。該方法針對集中力載荷，實現了多源動態載荷識別，并采用L曲線法確定了最優輸出跟蹤中性能指標的關鍵參數。數值算例結果表明了本文提出的基于最優輸出跟蹤的動態載荷識別方法的有效性和可行性。該方法在載荷作用的全部時間域內整體進行反求，應用方便，計算量較小且效率較高，并具有較強的抑制噪聲能力。

1 載荷識別的正問題

結構動力響應的有限元離散形式為[7]

(1)

可將式(1)改寫為

(2)

假設一空間向量

(3)

建立相應的系統方程

(4)

該系統結構響應的輸出為

y=Dw

(5)

若實際測點響應量為z，(z=[z1z2…zn]，n為測點個數)。該響應可以為位移響應或速度響應。如果實測響應量z已知，需求結構的外載荷，就轉化為最優跟蹤問題，即尋求最優控制u，使得結構響應輸出y(t)盡量跟蹤上已知實測的響應輸出z(t)。

記e為跟蹤誤差，即

e(t)=y(t)-z(t)

(6)

2 基于最優輸出跟蹤的載荷識別

最優輸出跟蹤問題是尋求一種最優控制，既能保證跟蹤誤差盡可能小，又要避免使用取值“過大”的控制，故選用下式的最優化性能指標[8]

(7)

式中t0為起始時間，tf為終止時間，f表示末值型性能指標參數，Q表示平均誤差度量性能指標參數，R表示最小能量控制參數，且?t∈[t0,tf],f≥0,Q≥0,R>0。本文假設這些性能指標參數都為常數。

由于載荷識別問題中結構的外載荷是肯定存在的，所以肯定存在最優控制。本文的載荷識別方法基于載荷位置已知，并通過實測獲得測點的響應來反求載荷。此外，該載荷識別問題注重時間歷程，對終端誤差不關注，進而可忽略終端響應量的誤差，即令f=0。因此，最優性能指標式(7)可改寫為

(8)

式中r=R/Q。

從式(4)和(8)易知哈密爾頓函數為

λTAw+λTBu

(9)

u*(t)=R-1BTλ

(10)

由極大值原理知共軛方程和橫截條件分別為

λ(tf)=0

(11)

將式(10)代入式(4)并結合式(11)得正則方程和端點條件

(12)

根據式(12)，可以假設

λ(t)=-p(t)w(t)+g(t)

(13)

微分式(13)并結合式(12)得

(14)

將式(13)代入式(11)，并聯合式(14)得

p(t)BR-1BTp(t)]w(t)+[p(t)BR-1BT-

(15)

由于式(15)對任意w(t)均成立，利用w(t)的任意性可知，p(t)和g(t)分別滿足

(16)

由式(11)及w(tf)的任意性得

p(tf)=0，g(tf)=0

(17)

從而可以得到最優跟蹤器為

(18)

進而得到最優跟蹤器相應的響應

(19)

因此，要最終求得最優跟蹤器u*(t)，即結構所受的外載荷，就要求解式(16)和(17)，即帶終端條件的Riccati微分方程以及帶初始條件的微分方程(19)。

通過上述方法求解微分方程求得p,g,x，即可求出最優跟蹤器，即結構所受的外載荷u*(t)。

3 L曲線法確定關鍵性能指標參數

利用上述方法進行載荷識別過程中，需要確定最優化性能指標中的指標參數，即式(8)中的權重系數r。當測量響應中存在誤差時，跟蹤誤差和所獲得的取值都會不同程度的變化，因此該權重系數r=R/Q的選取至關重要，既要保證跟蹤誤差盡可能小，又要避免取值“過大”。為了能有效地確定r值，本文采用適應性強的L曲線法以獲取最佳的r值。

L曲線法是指用對數尺度來描述殘差范數和解的范數的曲線對比[10]，該方法的典型特征是尺度圖形中出現一條明顯的L曲線，曲線的拐點所對應的參數當作優化參數。式(8)中eT(t)e(t)和uT(t)u(t)都是權重系數r的函數，選擇不同的r值，以eT(t)e(t)的對數為橫坐標，uT(t)u(t)的對數為縱坐標，得到許多點(lg(eT(t)e(t))，lg(uT(t)u(t)))，經過曲線擬合得到一條曲線。再利用L曲線法選擇最佳的r值，其中本文采用最短距離來定位L曲線的最優點。

4 算 例

為了驗證上述基于最優輸出跟蹤的載荷識別方法的正確性和有效性，下面給出3個多源載荷識別的數值算例。在數值算例中，測點的位移響應通過有限元數值仿真計算得到。此外，測點的選擇至關重要，該測點的響應要對外載荷具有較強的敏感性，則一般通過敏感性分析來選取合適的測點，且測點的個數要大于等于外載荷的數目，以避免載荷識別時出現不適定性問題。載荷識別前，在仿真計算得到的位移響應中加入一定水平的隨機噪聲以模擬實驗測量的響應，此時帶噪聲的位移響應可用下式表示

Yerr=Ycal+lnoise·std(Ycal)·rand(-1,1)

(20)

式中Ycal為計算得到的位移響應；std(Ycal)為計算的位移響應的標準差；lnoise為百分數，代表噪聲水平；rand(-1,1)為[-1,1]之間的隨機數。

4.1 算例 1

如圖1所示一均勻的四層剪切樓板結構[11]，每一層樓板有m=175.126 kg，k=57.328 kN/m，并在第4層樓板和第2層樓板分別施加指數衰減載荷F1和三角形載荷F2的剪切力模擬瞬態沖擊載荷，以及在第1層和第3層分別布置測點。

圖1 結構簡圖

按結構動力學的方法可寫出結構對應的質量陣和剛度陣：

當位移響應中未加入隨機噪聲時，利用上述載荷識別的方法進行載荷反求重構，其結果如圖2所示。從圖中可以看出，在測量響應沒有噪聲污染的情況下，本方法可以很好地實現動態載荷時程的重構。同時，確定最佳r值的L曲線如圖3所示。

當在測點位移響應加入5%水平的隨機噪聲時，利用本文所述方法進行載荷反求重構的結果如圖4所示。從圖中可以看出，在響應數據受噪聲影響的情況下，本方法仍可較好地進行動態載荷識別。圖5給出了確定最佳r值的L曲線。

表1給出了當測點位移含有5%水平的隨機噪聲時，9個不同時間點下識別的載荷值及其相對誤差。從表中可以看出，本方法在測量響應受到較低水平噪聲影響下識別載荷結果的相對誤差均在10%以內。從上述結果可看出，本文所述的載荷識別方法能有效地在時域上重構出施加在結構上的多源載荷。

圖2 多源載荷識別的結果(零噪聲水平)

圖3 確定最佳r值的L曲線(0%噪聲水平)

圖4 多源載荷識別的結果(5%噪聲水平)

圖5 確定最佳r值的L曲線(5%噪聲水平)

表1 不同時刻識別的動態載荷及其相對誤差(5%噪聲水平)

4.2 算例2

圖6所示為一個高空索道塔架結構有限元模型。該結構立柱下端材料為等邊角鋼∠200×16，上端為等邊角鋼∠200×14，橫桿、斜桿、撐桿下端為不等邊角鋼∠160×100×12，上端為不等邊角鋼∠140×90×10，其他用不等邊角鋼∠100×80×8，材料特性參數為密度7 860 kg/m3，彈性模量210 GPa，泊松比0.33，阻尼假設為比例阻尼，與質量陣相關的系數為10，與剛度陣相關的系數為2×10-4。整個塔架結構采用梁單元建模，塔架的四個支腳固支。在圖6中所示的箭頭處分別施變頻正弦載荷F1和隨機動態載荷F2，這些載荷模擬了高空索道在工作過程中鋼索對塔架的垂向作用力。

圖6 高空索道塔架結構的有限元模型

圖6中圓點處標出了響應測量的位置，通過仿真計算，可以得到測點的位移響應，在仿真得到的位移響應中加入10%水平的隨機噪聲來模擬實驗測量的響應。利用本文提出的載荷識別方法進行載荷時間歷程重構，其結果如圖7所示。從圖中可以看出，在測量響應數據受到較高噪聲水平影響下，本文方法能夠正確有效地實現動態載荷的識別。圖8給出了確定最佳r值的L曲線。

圖7 多源載荷識別的結果(10%噪聲水平)

圖8 確定最佳r值的L曲線(10%噪聲水平)

4.3 算例3

圖9所示一個25桿桁架結構，其中桿(1)～(4)有相同的橫截面積A1=400 mm2，桿(16)～(25)、桿(11)～(15)和桿(5)～(10)的橫截面積分別為A2=500 mm2，A3=600 mm2和A4=800 mm2。橫向和縱向桿的長度L=15.24 m，桿的彈性模量為200 GPa，泊松比為0.33，密度為2 800 kg/m3。阻尼假設為比例阻尼，與質量陣相關的系數為0，與剛度陣相關的系數為0.003。連接點12為餃接支座，6，8和10為滾動支座。在3個不同的節點位置分別施加不同時間歷程的載荷，其中一種為兩個周期的正弦載荷，另外兩種分別為一個周期的三角載荷和半個周期的三角載荷，其作用位置分別如圖9中箭頭F1，F2和F3所示。

建立結構的有限元模型，采用桿單元劃分網格，共25個單元和12個節點。選擇測點響應為連接點2處的橫向位移、連接點3處的縱向位移和連接點7的橫向位移。

圖9 25桿桁架結構

為了測試該算法對噪聲的適應能力，在位移響應中加入15%水平的隨機噪聲作為測量響應，該動態載荷識別的結果如圖10所示。從圖中可以看出，在測量響應數據受到較高噪聲水平影響下，本文載荷識別的方法可以較準確地實現動態載荷時程的重構。圖11給出了確定最佳r值的L曲線。

圖10 多源載荷識別的結果(15%噪聲水平)

圖11 確定最佳r值的L曲線(15%噪聲水平)

表2給出了在15%噪聲水平下3個不同時刻識別載荷的相對誤差。識別載荷的相對誤差都在12%以內，這說明本文所述的多源載荷識別方法能有效地抑制噪聲對識別結果的影響，具有較好的穩健性。

表2 15%噪聲水平下不同時刻識別載荷的相對誤差

5 結 論

工程實際中動態載荷的識別是一個難度較高且較為復雜的反問題，有許多共性問題需要亟待解決。本文基于最優輸出跟蹤的思想，提出了一種多源動態載荷識別的方法，并采用L曲線法確定了性能指標中的關鍵參數。數值算例表明該方法能有效抑制噪聲對載荷識別結果的影響，能較為準確、穩健地實現載荷的反演識別。該方法計算簡單快捷、抗噪聲能力較強，具有一定的工程實際實用價值。

參考文獻：

[1] Bartlett F D, Flannelly W D. Modal verification of force determination for measuring vibration loads[J]. Journal of the American Helicopter Society， 1979: 10—18.

[2] Liu Y, Shepard WS. Dynamic force identification based on enhanced least squares and total-squares schemes in the frequency domain[J]. Journal of Sound and Vibration, 2005, 282: 37—60.

[3] 劉杰. 動態載荷識別的計算反求技術研究 [D]. 長沙：湖南大學, 2011.Liu Jie. Research on computational inverse techniques in dynamic load identification [D]. Changsha:Hunan University, 2011.

[4] 韓旭, 劉杰, 李偉杰, 等. 時域內多源動態載荷的一種計算反求技術 [J]. 力學學報, 2009, 41(4): 598—605.Han Xu, Liu Jie, Li Weijie, et al. A computational inverse technique for reconstruction of multisource loads in time domain [J]. Chinese Journal of Theoretical and Applied Mechanics, 2009, 41(4): 598—605.

[5] 伍乾坤, 韓旭, 劉杰, 等. 一種直接求解的動態載荷識別方法[J]. 應用力學學報, 2011, 28(2): 201—205.Wu Qiankun, Han Xu, Liu Jie, et al. A dynamic load identification method based on direct solution [J]. Chinese Journal of Applied Mechanics, 2011, 28(2): 201—205.

[6] 章繼峰, 張博明, 杜善義. 最優化方法在動態載荷時程重構中的應用研究[J]. 振動與沖擊, 2006, 25(4): 5—8.Zhang Jifeng, Zhang Boming, Du Shanyi. Application of optimization method to kinetic loading time history reconstruction [J]. Journal of Vibration and Shock, 2006, 25(4): 5—8.

[7] Paultre P. Dynamics of Structures [M]. Wiley-ISTE, 2010.

[8] Anderson B D O, Moore J B. Linear Optimal Control[M]. Prentice-Hall Inc, 1971.

[9] 李榮華, 劉播. 微分方程數值解法[M]. 北京: 高等教育出版社, 2009.Li Ronghua, Liu Bo. Numerical Solution of Differential Equations [M]. Beijing: Higher Education Press, 2009.

[10] Hansen P C, Leary D P. The use of the L-curve in the regularization of discrete ill-posed problems [J]. SIAM Journal on Scientific and Statistical Computing, 1993, 14:1 487—1 503.

[11] 帕茲 M. (美)著．結構動力學——理論與計算[M]. 李裕澈, 劉勇生，譯. 北京: 地震出版社, 1993.Paz M. Structural Dynamics-Theory and Computation [M]. Beijing: Publish House of Seismological Bureau, 1993.