高　偉， 于開平， 林　宏

( 1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱　150001；　2. 北京宇航系統(tǒng)工程研究所，北京　100076 )

?

基于商函數(shù)的動態(tài)載荷識別最優(yōu)正則化參數(shù)選取方法

高偉1， 于開平1， 林宏2

( 1. 哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱150001；2. 北京宇航系統(tǒng)工程研究所，北京100076 )

為解決載荷識別反問題，研究選取最優(yōu)正則化參數(shù)商函數(shù)方法。利用Tikhonov正則化方法的最優(yōu)化問題的最小二乘解，定義以正則化參數(shù)為自變量的商函數(shù)；根據(jù)不同的正則化參數(shù)，使用二次規(guī)劃理論，求解Tikhonov正則化方法的最優(yōu)化問題的最優(yōu)解；基于不同最優(yōu)解對應(yīng)商函數(shù)的不同特點，將最優(yōu)正則化參數(shù)的商函數(shù)方法，與廣義交叉檢驗(Generalized Cross-Validation, GCV)準(zhǔn)則所得載荷識別結(jié)果進行比較。數(shù)值仿真及試驗驗證結(jié)果表明，商函數(shù)方法對于共振區(qū)及非共振區(qū)下載荷識別問題具有較好的合理性和適用性。

載荷識別； 不適定問題； 正則化方法； 二次規(guī)劃； 形函數(shù)

0　引言

載荷識別和系統(tǒng)辨識是反問題理論中的兩類問題[1-3]。Tikhonov正則化方法[4]利用引入的正則化參數(shù)構(gòu)造正則化算子，能夠有效克服由系統(tǒng)矩陣的病態(tài)性引起的載荷識別不適定問題。最優(yōu)正則化參數(shù)的選取是載荷識別理論中主要的研究方向之一，主要的正則化參數(shù)選取方法有嶺跡法[5]、擬最優(yōu)準(zhǔn)則[6]、廣義交叉檢驗準(zhǔn)則(Generalized Cross-Validation，簡稱GCV)[7]、Morozov偏差原理[8]、啟發(fā)式規(guī)則(Heuristic rules)[9]和L曲線(L-curve)準(zhǔn)則[10]。正則化參數(shù)選取主要分為先驗策略和后驗策略兩類，后驗策略不需要事先知道最優(yōu)正則化參數(shù)的先驗信息。廣義交叉檢驗準(zhǔn)則和L曲線準(zhǔn)則作為后驗策略是工程上最常用的兩種最優(yōu)正則化參數(shù)選取方法，工程應(yīng)用結(jié)果[11]表明兩種方法具有更好的適應(yīng)性。L曲線準(zhǔn)則沒有收斂性證明，并且有人給出L曲線非收斂的反例[12]。L曲線方法需要對應(yīng)不同正則化參數(shù)的坐標(biāo)點畫出對數(shù)尺度下的L形曲線，坐標(biāo)點過少時易導(dǎo)致結(jié)果精度低，坐標(biāo)點過多時易導(dǎo)致計算量過大；在個別情況下，L曲線呈現(xiàn)凹形曲線[13]，給曲線的最大曲率位置確定帶來困難。雖然在理論上已經(jīng)證明GCV準(zhǔn)則的收斂性，但是有時用來確定最優(yōu)正則化參數(shù)的GCV函數(shù)曲線過于平坦，導(dǎo)致確定GCV函數(shù)的最小值出現(xiàn)困難[14]。

筆者研究一種基于二次規(guī)劃理論的確定最優(yōu)正則化參數(shù)的商函數(shù)方法。它屬于后驗策略方法，不需要事先知道結(jié)構(gòu)響應(yīng)中的噪聲水平。數(shù)值仿真及試驗驗證表明，商函數(shù)方法對于共振區(qū)及非共振區(qū)下載荷識別問題具有較好的適應(yīng)性。

1　正則化方法

離散系統(tǒng)方程的建立基于Duhamel積分：

(1)

的離散化。將載荷作用時間區(qū)域離散化為Q個長度為Δt的時間單元，將真實載荷f(t)用最小二乘近似函數(shù)表示[15]，設(shè)時間域內(nèi)局部支撐域Ω內(nèi)包含S-1時間單元(共S點)，則支撐域Ω內(nèi)形函數(shù)向量為

(2)

取Ω內(nèi)中間時間單元近似載荷作為真實載荷的最優(yōu)最小二乘近似，設(shè)S=2l，定義首個時間單元中形函數(shù)為

(3)

可以將所有時間單元上近似載荷表示為

(4)

則整個時間域上最小二乘近似載荷為

(5)

將式(5)代入式(1)并整理有

(6)

式中：fij為fi中第j分量。

再令

(7)

表示形函數(shù)Nj(t)的響應(yīng)，其中j=1,2,…S。對式(6)在時間域上離散化，由關(guān)于fij的向量長度為Q+2l-1的向量得

(8)

則每個形函數(shù)對應(yīng)的矩陣為

(9)

l-1 columns

(10)

(11)

(12)

由式(10-12)得離散系統(tǒng)方程模型：

y=Gf,

(13)

式(13)中系統(tǒng)矩陣往往是病態(tài)的。

正則化方法能夠有效克服系統(tǒng)矩陣的病態(tài)導(dǎo)致的不適定問題，主要思想為考察最優(yōu)化問題：

(14)

式中：α為常數(shù)，α>0；‖·‖為歐氏范數(shù)。式(14)的最小二乘解表達(dá)式為

(15)

式中：I為單位矩陣。將響應(yīng)y=yture+ynoise及系統(tǒng)矩陣G的SVD (Singular Value Decomposition)分解[12]并代入式(15)：

(16)

通過對正則化參數(shù)α進行合理取值，使得識別載荷fα兼顧數(shù)值反演的穩(wěn)定性和與真實載荷誤差的可控性。

在載荷識別問題的正則化求解過程中，采用不同的最優(yōu)正則化參數(shù)選取方法，將得到不同的最優(yōu)正則化參數(shù)，以及不同精度的識別結(jié)果。由式(16)可知，較小的正則化參數(shù)可以保留更多響應(yīng)ytrue的初始信息，并且更接近真實系統(tǒng)。在識別結(jié)果精度相同的情況下，能夠得到較小最優(yōu)正則化參數(shù)的最優(yōu)正則化參數(shù)選取方法在實際工程中具有更好的適用性。

2　最優(yōu)正則化參數(shù)選取方法

利用二次規(guī)劃算法[16-18]求解最優(yōu)化問題(14)，矩陣GTG+α2I在參數(shù)α較大時為非病態(tài)，求解得到的識別載荷具有非常高的精度且滿足式(15)，整理有：GTy-GTGfα=α2fα，等式兩邊取歐氏范數(shù)有‖GTy-GTGfα‖=α2‖fα‖。隨著常數(shù)α逐漸減小，矩陣GTG+α2I從非病態(tài)逐漸變?yōu)榫哂休^弱的病態(tài)性，識別載荷具有較好的精度，則等式退化為‖GTy-GTGfα‖≈α2‖fα‖；如果參數(shù)α繼續(xù)減小，則矩陣GTG+α2I的病態(tài)性較強，識別載荷精度非常差，則有‖GTy-GTGfα‖≠α2‖fα‖。設(shè)商函數(shù)H(α)為

(17)

商函數(shù)方法實現(xiàn)步驟：

(1)由SVD分解得到系統(tǒng)矩陣的最大奇異值σ1，給出初始參數(shù)向量α1=[σ1/50,2σ1/50,…σ1]；

(2)對于向量α1中每個分量，計算相應(yīng)H(α)函數(shù)值構(gòu)成的向量H=[H(σ1/50),H(2σ1/50),…H(σ1)]；

(3)給出閾值eop，如果H中所有分量與1的距離均不大于閾值eop，則更新參數(shù)向量α2=α1/2；

(4) 重復(fù)步驟(2)和(3)，直到向量H中出現(xiàn)某個分量(假設(shè)為第q個)與1的距離大于閾值eop，并且接下來有連續(xù)多個分量(可取5個)與1的距離不大于eop時，假設(shè)參數(shù)向量αi=α1/2i-1，確定的最優(yōu)正則化參數(shù)即為向量αi中的第q+1個分量。

3　數(shù)值仿真

圖1　懸臂梁有限元模型Fig.1 The finite element model of cantilever beam structure

采用懸臂梁有限元模型模擬鋼制懸臂梁結(jié)構(gòu)(見圖1)。懸臂梁劃分為10個單元，從左到右依次為節(jié)點1到10。

鋼制懸臂梁結(jié)構(gòu)的一階固有頻率為9.3 Hz，采用數(shù)值仿真分別研究懸臂梁結(jié)構(gòu)在非共振區(qū)及共振區(qū)的載荷識別問題。取基函數(shù)向量為q(t)=[1,t,t2]T，即R=3，采樣點個數(shù)S=4[15]。確定式(15)形式的形函數(shù)，構(gòu)造y=Gf形式的離散系統(tǒng)方程，并利用正則化方法求解。分別使用商函數(shù)方法及GCV方法確定最優(yōu)正則化參數(shù)，并比較商函數(shù)及GCV方法的結(jié)果。

3.1單輸入單輸出(SISO)系統(tǒng)載荷識別

在節(jié)點8處施加載荷，采樣頻率Fs=1 000 Hz，使用節(jié)點5處的加速度響應(yīng)識別載荷，取閾值eop=10-4，對響應(yīng)數(shù)據(jù)添加噪聲：

(18)

式中：y為用于載荷識別的加速度響應(yīng)；ycal為通過有限元模型計算得到的加速度響應(yīng)；lp為噪聲水平；noise為與ycal相同長度的均值為0、方差為1的白噪聲序列；std(·)為標(biāo)準(zhǔn)差。

識別載荷與真實載荷freal之間相對誤差為

(19)

3.1.1SISO系統(tǒng)非共振區(qū)載荷識別

圖2　SISO系統(tǒng)在2%噪聲水平下商函數(shù)Fig.2 QF under 2% noise level of SISO system

施加正弦載荷f1(t)=40sin(40πt)，作用時間為0.30 s，分別在2%及5%噪聲水平下利用正則化方法識別載荷，利用商函數(shù)方法確定最優(yōu)正則化參數(shù)。在2%噪聲水平下商函數(shù)曲線見圖2，其中星號處為最優(yōu)正則化參數(shù)對應(yīng)坐標(biāo)點。由圖2可知，當(dāng)α≥0.017 5時，幾乎所有函數(shù)值為1；當(dāng)α<0.017 5時，函數(shù)值與1相差非常大，因此最優(yōu)正則化參數(shù)α=0.017 5。在2%噪聲水平下商函數(shù)及GCV方法的識別載荷見圖3，在2%和5%噪聲水平下商函數(shù)及GCV方法的非共振區(qū)識別結(jié)果見表1。由圖3可知，兩種方法的識別載荷與真實載荷擬合結(jié)果較好。

3.1.2SISO系統(tǒng)共振區(qū)載荷識別

施加正弦載荷f2(t)=40sin(20πt)，作用時間為0.40 s。分別在2%及5%噪聲水平下利用正則化方法識別載荷，利用商函數(shù)方法確定正則化參數(shù)。在2%噪聲水平下商函數(shù)及GCV方法的識別載荷見圖4。在2%和5%噪聲水平下商函數(shù)及GCV方法的共振區(qū)識別結(jié)果見表2。

圖3　SISO系統(tǒng)在2%噪聲水平下識別載荷f1Fig.3 Identified load f1 under 2% noise level of SISO system

圖4　SISO系統(tǒng)在2%噪聲水平下識別載荷f2Fig.4 Identified load f2 under 2% noise level of SISO system

Table 1 Results of load identification in the non-resonant region of SISO system

噪聲水平/%商函數(shù)方法GCV方法最優(yōu)參數(shù)RErr/%最優(yōu)參數(shù)RErr/%20.01757.700.10208.7750.078610.150.102010.88

表2SISO系統(tǒng)共振區(qū)載荷識別結(jié)果

Table 2 Results of load identification in the resonant region of SISO system

噪聲水平/%商函數(shù)方法GCV方法最優(yōu)參數(shù)RErr/%最優(yōu)參數(shù)RErr/%20.011111.930.374012.8050.027713.910.371016.25

由圖4可知，兩種方法的識別載荷能較好地反映真實載荷的時間歷程。由表1和表2可知，商函數(shù)方法對SISO系統(tǒng)載荷識別結(jié)果的噪聲有較好的適應(yīng)性。在零初始條件下，結(jié)構(gòu)在周期激勵作用下的響應(yīng)中包含穩(wěn)態(tài)周期項及過渡衰減項。在共振區(qū)響應(yīng)中，過渡衰減項初始時刻的幅值較大，到達(dá)穩(wěn)態(tài)需要的時間較長；與非共振區(qū)載荷識別結(jié)果比較，共振區(qū)識別載荷精度相對較差，同時造成對測量響應(yīng)中測量噪聲有較大敏感性。

3.2多輸入多輸出(MIMO)系統(tǒng)載荷識別

使用與SISO系統(tǒng)載荷識別情形下相同的初始條件，分別在節(jié)點8和節(jié)點6處施加正弦載荷及沖擊載荷，利用半正弦波模擬沖擊載荷[19]，采樣頻率Fs=1 000 Hz。高水平的噪聲對MIMO系統(tǒng)載荷識別結(jié)果有較大影響[20]，其中沖擊載荷的主要特征為最大幅值[21]，對應(yīng)的RErr指的是識別載荷最大幅值與真實載荷最大幅值之間的RErr。

3.2.1MIMO系統(tǒng)非共振區(qū)載荷識別

施加正弦載荷f1(t)=40sin(40πt)，作用時間為0.25 s；沖擊載荷為f3(t)=45sin(100πt)，取0.16～0.17 s時間段上長度為0.01 s的半正弦波模擬沖擊載荷，沖擊載荷時間歷程為0.25 s。使用節(jié)點7和節(jié)點9處的加速度響應(yīng)，分別在1%和2%噪聲水平下識別載荷。在1%噪聲水平下利用GCV方法確定最優(yōu)正則化參數(shù)，對應(yīng)識別載荷見圖5。由圖5可知，沖擊載荷識別效果很好；正弦載荷識別精度非常差，最優(yōu)正則化參數(shù)過大導(dǎo)致丟失過多的系統(tǒng)信息。

與SISO系統(tǒng)載荷識別結(jié)果比較，MIMO系統(tǒng)載荷識別對應(yīng)式(13)中方程的系統(tǒng)矩陣具有更高的階數(shù)及更強的病態(tài)性。降低對商函數(shù)值近似為1的標(biāo)準(zhǔn)，在MIMO系統(tǒng)載荷識別過程中閾值為eop=10-2。在1%噪聲水平下商函數(shù)方法識別載荷見圖6。在1%和2%噪聲水平下商函數(shù)及GCV方法的非共振區(qū)識別結(jié)果見表3。由圖6可知，識別載荷較好地反映真實載荷的時間歷程。

表3　MIMO系統(tǒng)非共振區(qū)載荷識別結(jié)果

圖5　MIMO系統(tǒng)非共振區(qū)在1%噪聲水平下GCV方法識別載荷Fig.5 Identified load of the GCV method under 1% noise level in the non-resonant region of MIMO system

圖6　MIMO系統(tǒng)非共振區(qū)在1%噪聲水平下商函數(shù)方法識別載荷Fig.6 Identified load of QFM in the non-resonant region under 1% noise level of MIMO system

3.2.2MIMO系統(tǒng)共振區(qū)載荷識別

施加正弦載荷f2(t)=40sin(20πt)，作用時間0.25 s；其余條件與3.2.1相同。使用節(jié)點7和節(jié)點9處的加速度響應(yīng)，分別在1%和2%噪聲水平下識別載荷。在1%噪聲水平下利用商函數(shù)方法確定最優(yōu)正則化參數(shù)，對應(yīng)識別載荷見圖7。在1%和2%噪聲水平下商函數(shù)及GCV方法的共振區(qū)識別結(jié)果見表4。由圖7可知，商函數(shù)方法的識別載荷與真實載荷的擬合結(jié)果較好。由表4可知，GCV方法確定的最優(yōu)正則化參數(shù)較大，丟失過多系統(tǒng)信息而導(dǎo)致識別載荷精度較差。由表3和表4可知，商函數(shù)方法可以通過調(diào)整閾值得到滿足不同系統(tǒng)狀態(tài)的最優(yōu)正則化參數(shù)，并得到精度較好的識別載荷。

表4　MIMO系統(tǒng)共振區(qū)載荷識別結(jié)果

4　試驗驗證

以鋼制懸臂梁結(jié)構(gòu)(見圖8)作為試驗?zāi)Ｐ?，幾何尺寸?.900 m×0.050 m×0.009 m，彈性模量為200 GPa，密度為7.8×103kg/m3。初始條件為零，將懸臂梁結(jié)構(gòu)劃分為10個單元，從左至右依次為節(jié)點1至10。在節(jié)點8處施加頻率為20 Hz的正弦形式載荷，并同時測量真實載荷數(shù)據(jù)，載荷作用時間為0.50 s。在節(jié)點3-6處放置加速度傳感器測量結(jié)構(gòu)響應(yīng)，采樣頻率為1024 Hz。

圖7　MIMO系統(tǒng)共振區(qū)在1%噪聲水平下商函數(shù)方法識別載荷Fig.7 Identified load of QFM in the resonant region under 1% noise level of MIMO system

圖8　鋼制懸臂梁結(jié)構(gòu)Fig.8 The cantilever beam structure

圖9　試驗驗證識別載荷Fig.9 The identified load in the experimental verification

試驗驗證使用與數(shù)值仿真離散系統(tǒng)方程形式相同的式(13)，使用節(jié)點6處加速度響應(yīng)，分別利用商函數(shù)及GCV方法進行載荷識別，并與真實載荷進行比較(見圖9)。

商函數(shù)方法確定的最優(yōu)正則化參數(shù)及識別載荷RErr分別為α=0.109 0及15.74%，GCV方法確定的最優(yōu)正則化參數(shù)及識別載荷RErr分別為α=0.202 0及19.98%。由于GCV方法確定的最優(yōu)正則化參數(shù)相對較大，在載荷識別過程中丟失較多系統(tǒng)信息；商函數(shù)方法確定的識別載荷與GCV方法確定的識別載荷比較，具有較高的精度，并且與真實載荷吻合較好。

5　結(jié)論

(1)商函數(shù)方法屬于后驗策略型方法，在共振區(qū)及非共振區(qū)下載荷識別問題中能有效確定最優(yōu)正則化參數(shù)，并得到高精度的識別載荷，說明商函數(shù)方法在工程上具有更好的適應(yīng)性。

(2)商函數(shù)方法的判定標(biāo)準(zhǔn)基于正則化方法的最小二乘解，不會出現(xiàn)陷入局部最優(yōu)解的問題。將商函數(shù)方法的函數(shù)值是否為1作為選取最優(yōu)正則化參數(shù)的判定標(biāo)準(zhǔn)，商函數(shù)值與閾值之間容易區(qū)分，能夠有效克服GCV方法最小值附近函數(shù)值過于接近而導(dǎo)致的確定最小值困難的問題。

(3)對于不同的實際工程需要，商函數(shù)方法可以通過調(diào)整閾值得到滿足需要的最優(yōu)正則化參數(shù)。

[1]文祥榮,智浩,孫守光.結(jié)構(gòu)動態(tài)載荷識別的精細(xì)逐步積分法[J].工程力學(xué),2001,18(4):117-122.

Wen X R, Zhi H, Sun S G. Highly precise direct integration scheme for structural dynamic load identification [J]. Engineering Mechanics, 2001,18(4):117-122.

[2]王曉燕,黃維平,李華軍.地震動反演及結(jié)構(gòu)參數(shù)識別的EKF算法[J].工程力學(xué),2005,22(4):20-23.

Wang X Y, Huang W P, Li H J. Inversion of ground motion and identification of structural parameters by EFK [J]. Engineering Mechanics, 2005,22(4):20-23.

[3]李東升,郭杏林.逆虛擬激勵法隨機載荷識別試驗研究[J].工程力學(xué),2004,21(2):134-139.

Li D S, Guo X L. Experimental random loading identification using inverse pseudo excitation method [J]. Engineering Mechanics, 2004,21(2):134-139.

[4]Tikhonov A N. Solution of incorrectly formulated problems and the regularization method [J]. Soviet Mathematics Doklady, 1963,4:1035-1038.

[5]Hoerl A E, Kennard R W. Ridge regression: Biased estimation for non-orthogonal problems [J]. Technometrics, Special 40th Anniversary Issue, 2000,42(1):80-86.

[6]Kindermann S, Neubauer A. On the convergence of the quasi optimality criterion for (iterated) Tikhonov regularization [J]. Inverse Problem, 2008,2(2):291-299.

[7]Golub G H, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter [J]. Technometrics, 1979,21:215-223.

[8]Bonesky T. Morozov's discrepancy principle and Tikhonov-type functionals [J]. Inverse Problems, 2009,25(1):015015.

[9]Jin B, Lorenz D A. Heuristic parameter-choice rules for convex variational regularization based on error estimates [J]. SIAM J. Numer. Anal, 2010,48(3):1208-1229.

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

[11]Hanke M, Hansen P C. Regularization method for large-scale problems [J]. Surv. Math. 1nd, 1993,3:253-315.

[12]Vogel C R. Non-convergence of the L-curve regularization parameter selection method [J]. Inverse Problems, 1996,12:535-547.

[13]Hansen P C. The L-curve and its use in the numerical treatment of inverse problems [M]. Computational Inverse Problems in Electrocardiology, WIT Press, Southampton, 2001:119-142.

[14]Varah J M. Pitfalls in the numerical solution of linear Ⅲ-posed problems [J]. SIAM Journal on Scientific and Statistical Computing, 1983,4:164-176.

[15]Liu J, Sun X S, Han X, et al. A novel computational inverse technique for load identification using the shape function method of moving least square fitting [J]. Computers and Structures, 2014,144:127-137.

[16]Coleman T F, Li Y. A reflective Newton method for minimizing a quadratic function subject to bounds on some of the variables [J]. Siam Journal on Optimization, 1993,6(4):1040-1058.

[17]Gill P E, Murray W, Wright M H. Practical optimization [M]. London: Academic Press Inc., 1981.

[18]Gould N, Toint P L. Preprocessing for quadratic programming [J]. Mathematical Programming, 2004,100(1):95-132.

[19]Choi K, Chang F K. Identification of impact force and location using distributed sensors [J]. AIAA Journal, 1996,34(1):136-142.

[20]Li Z, Feng Z P, Chu F L. A load identification method based on wavelet multi-resolution analysis [J]. Journal of Sound and Vibration, 2014,333:381-391.

[21]Mao B Y, Xie S L, Xu ML, et al. Simulated and experimental studies on identification of impact load with the transient statistical energy analysis method [J]. Mechanical Systems and Signal Processing, 2014,46:307-324.

2015-11-13；編輯：任志平

國家自然科學(xué)基金項目(11372084)

高偉(1981-)，男，博士研究生，副教授，主要從事載荷識別理論方面的研究。

于開平,E-mail: yukp@hit.edu.cn

10.3969/j.issn.2095-4107.2016.02.014

O327

A

2095-4107(2016)02-0105-07