基于測點(diǎn)優(yōu)選和改進(jìn)L曲線法的動載荷識別

李曉旺， 趙海濤， 陳吉安

(上海交通大學(xué) 航空航天學(xué)院， 上海 200240)

動態(tài)載荷識別[1]是結(jié)構(gòu)動力學(xué)中的第2類反問題.目前，載荷識別方法主要有頻域法、時域法以及其他一些智能算法[2-4].對于每一種算法而言，其最終目的都是要建立起輸入載荷和輸出響應(yīng)以及系統(tǒng)特性之間的數(shù)學(xué)模型.由于反求輸入載荷的過程需要對數(shù)學(xué)模型中的傳遞函數(shù)矩陣求逆，并且被測響應(yīng)容易受到環(huán)境噪聲的干擾，往往會導(dǎo)致出現(xiàn)不適定問題.

為了盡量減小傳遞函數(shù)矩陣的病態(tài)性和響應(yīng)噪聲對載荷識別精度的負(fù)面影響，國內(nèi)外學(xué)者做了大量的研究工作.這些研究總體上可分為兩種思路：第1種思路為通過優(yōu)化響應(yīng)點(diǎn)的位置[5-6]獲得病態(tài)程度較低的傳遞函數(shù)矩陣，然后直接求廣義逆；第2種思路為避免對傳遞函數(shù)矩陣直接求逆，而是通過一定的技術(shù)手段對振動方程進(jìn)行特殊處理，間接反演動態(tài)激勵，比較典型的方法有奇異值分解(SVD)法[7]和 Tikhonov 正則化方法[8].對于正則化方法來說，如何選擇正則化參數(shù)是關(guān)鍵，而L曲線法[9-10]是最為常用的計(jì)算正則化參數(shù)的方法.

然而，目前的研究鮮見將兩種思路結(jié)合起來.對此，本文提出在時域中分2步解決不適定問題.第1步對不同測點(diǎn)對應(yīng)的傳遞函數(shù)矩陣進(jìn)行分析，找出條件數(shù)最小的傳遞函數(shù)矩陣，從而獲得最佳測點(diǎn).第2步采用Tikhonov正則化重構(gòu)載荷.同時，針對正則化方法精度不高的缺點(diǎn)，在確定正則化參數(shù)的過程中對L曲線法進(jìn)行改進(jìn)，即引入B樣條函數(shù)[11]對L曲線進(jìn)行插值，并通過計(jì)算B樣條曲線的各點(diǎn)曲率獲得更加準(zhǔn)確的正則化參數(shù).

1 基于Taylor多項(xiàng)式迭代的載荷識別方法

采用Taylor多項(xiàng)式迭代法[12]在時域內(nèi)識別載荷時間歷程.對于一個多輸入多輸出的振動系統(tǒng)，其振動微分方程可表示為

(1)

為了將振動方程解耦，總時間歷程被離散成n個時刻.對于任意相鄰的第i時刻t和第i+1時刻t+Δt，t+Δt時刻的振動響應(yīng)被近似地表示為t時刻響應(yīng)的Taylor多項(xiàng)式展開形式，即

(2)

(3)

(4)

(5)

式中：A1,A2,Ad,Av,Aa,B1,B2,Bd,Bv,Ba,D1,D2,Dd,Dv,Da均為m×m階定常矩陣；下標(biāo)a，d，v分別為加速度、位移和速度；F為從t1～tn+1時刻的全部激勵組合.

設(shè)結(jié)構(gòu)自由度數(shù)為m，則對式(5)作遞歸可得到總線性離散方程為

(6)

式中：Yl為從t1～tn時刻的全部響應(yīng)組合；Hl為總體傳遞函數(shù)矩陣.

2 載荷時間歷程重構(gòu)

2.1 響應(yīng)測點(diǎn)的優(yōu)選方法

(7)

通過遍歷全部組合情況可以找到條件數(shù)最小的傳遞函數(shù)矩陣H，此時的測點(diǎn)組合Y即為最佳選擇，從而實(shí)現(xiàn)了響應(yīng)測點(diǎn)的優(yōu)化配置.

2.2 改進(jìn)的L曲線法

在獲得最佳測點(diǎn)的響應(yīng)值并確定了條件數(shù)最小的傳遞函數(shù)矩陣之后，采用Tikhonov正則化技術(shù)重構(gòu)式(6)中的輸入載荷.設(shè)系統(tǒng)誤差e為

e=Y-HF

(8)

引入罰函數(shù)J的概念

J=(eTe)+λ(FTF)

(9)

當(dāng)J對F的1階導(dǎo)數(shù)為0，誤差e達(dá)到最小值.此時激勵F可以表示為

F=(HTH+λI)-1HTY

(10)

式中：I為單位矩陣；λ為正則化參數(shù).

正則化技術(shù)的核心在于如何選擇合適的正則化參數(shù)λ，L曲線法是常用的一種計(jì)算λ的方法.該方法的思路是選取一系列不同的λ值，分別計(jì)算每一個λ所對應(yīng)的lg‖Y-HF‖和lg‖F(xiàn)‖值.然后以lg‖Y-HF‖為橫坐標(biāo)，lg‖F(xiàn)‖為縱坐標(biāo)作圖，其曲線形狀類似于字母“L”.在L曲線的拐點(diǎn)位置可找出曲率最大的一個點(diǎn)，該點(diǎn)對應(yīng)的λ值即為最合適的正則化參數(shù).

然而，傳統(tǒng)的L曲線法存在一定的缺陷.首先，該方法的第1步對λ值的選取具有一定的盲目性，很可能會遺漏最準(zhǔn)確的λ值.其次，由于受到矩陣H奇異性的影響，對lg‖Y-HF‖和lg‖F(xiàn)‖求1階和2階導(dǎo)數(shù)很容易偏離真實(shí)值，從而直接造成L曲線的曲率計(jì)算結(jié)果不可靠.

為彌補(bǔ)上述不足，本文建立了基于B樣條插值的改進(jìn)L曲線算法.首先依據(jù)常規(guī)方法獲得L曲線，然后截取L曲線的拐點(diǎn)部分(L-Co)進(jìn)行插值操作，最后計(jì)算B樣條各點(diǎn)曲率從而獲得最佳正則化參數(shù).

2.2.1L-Co的B樣條插值 設(shè)L-Co共有Nc個坐標(biāo)，第j個坐標(biāo)為(xj，yj)，其對應(yīng)的正則化參數(shù)為λj，則共可以組成Nc個B樣條插值函數(shù)的控制頂點(diǎn)，任意一個控制頂點(diǎn)Pi可表示為

(11)

(i=0,1,…,Nc-1)

對L-Co進(jìn)行k次B樣條插值后，曲線可寫為

(12)

0≤u≤1

式中：Ni,k(u)為定義在節(jié)點(diǎn)矢量U=[u0u1…uNc+k]上的k次B樣條基函數(shù)，Ni,k(u)的計(jì)算公式為

(13)

(14)

經(jīng)過B樣條插值操作后，L-Co涵蓋的λ值遠(yuǎn)遠(yuǎn)超過預(yù)先選取的控制點(diǎn)的數(shù)量，避免了遺漏最佳正則化參數(shù)的可能.與選取同樣數(shù)量λ值計(jì)算L曲線坐標(biāo)的傳統(tǒng)方法相比，B樣條插值操作的計(jì)算量要小的多.

2.2.2L-Co的各點(diǎn)曲率求解 經(jīng)過插值操作后，計(jì)算L-Co各點(diǎn)的曲率問題可以轉(zhuǎn)化為計(jì)算B樣條函數(shù)的曲率問題.對k次B樣條曲線求1階導(dǎo)數(shù)的公式為

(15)

令

則式(15)可以改寫為

(16)

將式(16)和(12)比較可知，k次B樣條曲線的1階導(dǎo)數(shù)恰好是一條新的k-1次B樣條曲線.因此，求S(u)的2階導(dǎo)數(shù)只需根據(jù)式(15)計(jì)算S′(u)的1階導(dǎo)數(shù)即可求出S″(u)的值.插值后的L-Co是一條二維平面曲線，所以S′(u)和S″(u)分別對應(yīng)坐標(biāo)(x′，y′)和(x″，y″)，則L-Co各點(diǎn)的曲率Cu的計(jì)算公式為

(17)

曲率最小值所對應(yīng)的λ值即為最優(yōu)正則化參數(shù)，代入式(10)可以重構(gòu)輸入載荷.

3 數(shù)值算例研究

為驗(yàn)證所提方法的可行性和有效性，建立了3個不同的數(shù)值算例.3個算例所施加載荷的數(shù)量和復(fù)雜程度逐漸增加，擾動噪聲的水平也依次升高.采用本文算法分別對3個算例的輸入載荷進(jìn)行重構(gòu)，并與最差測點(diǎn)識別結(jié)果以及奇異值分解法識別結(jié)果進(jìn)行對比.

3.1 正弦載荷識別算例

算例(T)1建立了一個懸臂梁模型，懸臂梁各項(xiàng)參數(shù)分別為彈性模量E=70 GPa，泊松比μ=0.35，密度ρ=2 700×103kg/m3，長寬高尺寸為600 mm×60 mm×30 mm.受正弦激勵的懸臂梁模型如圖1所示，對懸臂梁模型施加2個正弦時變載荷，2個動態(tài)激勵可表示為

(18)

圖1 受正弦激勵的懸臂梁模型

圖2 T1中不同測點(diǎn)組合對應(yīng)的條件數(shù)

根據(jù)圖2可知，第209種組合情況可獲得最小條件數(shù)為97，對應(yīng)最佳測點(diǎn)為7、10、12；第45種組合情況可獲得最大條件數(shù)為441，對應(yīng)最差測點(diǎn)為1、7、12.

在分別獲得最佳測點(diǎn)和最差測點(diǎn)的響應(yīng)值前提下，對被測響應(yīng)混入5%的Gaussian白噪聲.采用改進(jìn)的L曲線法計(jì)算最佳測點(diǎn)下的最佳正則化參數(shù)λ.T1中的整體L曲線如圖3所示，對L-Co進(jìn)行B樣條插值的結(jié)果如圖4所示，所有λ對應(yīng)的L-Co曲率如圖5所示，其中Cu為曲率.

由圖5可知，當(dāng)λ=9.87×10-6時，L-Co的曲率達(dá)到最大值.將該λ值代入式(10)，即可反求輸入載荷的時間歷程.本文方法和兩種對照方法的識別結(jié)果如圖6所示，其識別誤差如表1所示.

由圖6可知，本文算法對正弦載荷的識別結(jié)果與真實(shí)值吻合得較好.由表1可知，本文算法的識別誤差低于另外兩種算法.同時，SVD法的識別誤差也保持在較低水平，低于選用最差測點(diǎn)的識別誤差.這說明對于載荷數(shù)量較少且形式簡單、噪聲影響水平較低的載荷識別問題來說，SVD法也可以適用，測點(diǎn)的選取對識別精度的影響要大于不同重構(gòu)算法的影響.

圖3 T1的L曲線圖

圖4 T1中L-Co的B樣條插值結(jié)果

圖5 T1中L-Co各點(diǎn)曲率計(jì)算結(jié)果

表1 T1載荷識別誤差

圖6 T1載荷時間歷程重構(gòu)結(jié)果

3.2 正弦疊加載荷識別算例

T2選用和T1相同的懸臂梁模型，該結(jié)構(gòu)受到3個正弦疊加載荷的激勵，受力情況如圖7所示.

3個動態(tài)激勵的表達(dá)式為

(19)

圖7 受正弦疊加激勵的懸臂梁模型

圖8 T2中不同測點(diǎn)組合對應(yīng)的條件數(shù)

依據(jù)圖8所示的計(jì)算結(jié)果，當(dāng)組合數(shù)為271時，得到最小條件數(shù)為112，對應(yīng)的最佳測點(diǎn)為2、7、9、11；當(dāng)組合數(shù)為1時，得到最大條件數(shù)為811，對應(yīng)的最差測點(diǎn)為1、2、3、4.

接下來通過仿真獲得最佳測點(diǎn)和最差測點(diǎn)的響應(yīng)值，對響應(yīng)混入10%的Gaussian白噪聲.采用改進(jìn)的L曲線法計(jì)算最佳正則化參數(shù)λ.圖9～11分別為T2的L曲線圖、B樣條插值后的L-Co圖，以及L-Co的曲率圖.

根據(jù)圖11中的曲率計(jì)算結(jié)果，當(dāng)λ=2.06×10-5時，L-Co的曲率獲得最大值.將該λ值代入式(10)重構(gòu)3個正弦疊加激勵的時間歷程.動載荷的識別結(jié)果如圖12所示，其識別誤差e如表2所示.

表2 T2載荷識別誤差

圖9 T2的L曲線圖

圖10 T2中L-Co的B樣條插值結(jié)果

圖11 T2中L-Co各點(diǎn)曲率的計(jì)算結(jié)果

圖12 T2載荷時間歷程重構(gòu)結(jié)果

由圖12可知，本文算法可以較為準(zhǔn)確地反演出正弦疊加激勵的時間歷程，其識別結(jié)果相比于另外兩種算法更靠近真實(shí)載荷.從表2可以看出，由于T2的載荷復(fù)雜度、載荷數(shù)量以及噪聲水平均高于T1，導(dǎo)致算法的識別誤差比T1有所提高.而從數(shù)據(jù)值來看，3個載荷的識別誤差均在6%以下，仍屬于較低水平且低于兩種對照算法的誤差，體現(xiàn)了所提算法的優(yōu)越性.

3.3 Gaussian白噪聲載荷識別算例

T3建立了一個桁架模型，其各項(xiàng)參數(shù)為E=200 GPa，μ=0.3，ρ=7 800×103kg/m3，截面積為20 cm2，所有水平桁架和豎直桁架的長度均為1 m.對桁架模型施加4個Gaussian白噪聲載荷，載荷幅值分別為80 N、70 N、60 N、50 N，如圖13所示.

由圖14可知，當(dāng)測點(diǎn)組合數(shù)為595時，條件數(shù)取到最小條件數(shù)為194，對應(yīng)的最佳測點(diǎn)為3、4、9、11、12；當(dāng)組合數(shù)為109時，條件數(shù)取到最大條件數(shù)為 1 215，其對應(yīng)的最差測點(diǎn)為1、2、7、10、12.

圖13 受Gaussian白噪聲激勵的桁架模型

圖14 T3中不同測點(diǎn)組合對應(yīng)的條件數(shù)

對被測響應(yīng)混入15%的Gaussian白噪聲，然后采用改進(jìn)的L曲線法求解正則化參數(shù)λ的最佳值.圖15～17分別為T3的L曲線圖、B樣條插值后的L-Co圖，以及L-Co的曲率圖.

由圖17可知，當(dāng)L-Co的曲率達(dá)到最大值時，對應(yīng)的最佳λ值為4.91×10-5.根據(jù)式(10)重構(gòu)輸入激勵，T3的載荷重構(gòu)結(jié)果如圖18所示，其識別誤差如表3所示.

與前兩個算例相比，T3的Gaussian白噪聲激勵是時變載荷中復(fù)雜程度最高的一種.從圖18中的被識別激勵時間歷程曲線可以看出， 本文算法重構(gòu)的載荷與真實(shí)值的吻合程度顯著高于其他兩種對照算法，說明該算法具有優(yōu)越的穩(wěn)健性和抗噪性.由表3可知，本文算法的識別誤差雖然比T1和T2稍高，但仍處于較低水平且明顯低于對照算法.由于載荷數(shù)量較多，導(dǎo)致傳遞函數(shù)矩陣維數(shù)較大，即使選取最佳測點(diǎn)也有較強(qiáng)的病態(tài)性，從而引起SVD法的識別誤差迅速增大，此時不同識別算法對精度的影響超過了測點(diǎn)的選取對精度的影響.

圖15 T3的L曲線圖

圖16 T3中L-Co的B樣條插值結(jié)果

圖17 T3中L-Co各點(diǎn)曲率計(jì)算結(jié)果

圖18 T3載荷時間歷程重構(gòu)結(jié)果

表3 T3載荷識別誤差

4 結(jié)語

針對載荷識別中的不適定問題，在時域內(nèi)建立了測點(diǎn)優(yōu)選和改進(jìn)的L曲線法相結(jié)合的算法.首先通過遍歷所有測點(diǎn)組合找出條件數(shù)最小的傳遞函數(shù)矩陣，進(jìn)而獲得對應(yīng)的最佳測點(diǎn).然后在使用正則化方法反求輸入載荷的過程中引入二次B樣條函數(shù)對L曲線進(jìn)行插值，從而獲得了更加精確的正則化參數(shù)，并將該算法的載荷重構(gòu)結(jié)果分別與最差測點(diǎn)識別結(jié)果以及SVD法識別結(jié)果作對比.數(shù)值算例結(jié)果表明，該方法可以有效減弱傳遞函數(shù)矩陣的病態(tài)性，并具有較強(qiáng)的穩(wěn)健性和抗噪性，從而提高了載荷識別的精度.