基于商函數的動態載荷識別最優正則化參數選取方法

高　偉， 于開平， 林　宏

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

?

基于商函數的動態載荷識別最優正則化參數選取方法

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

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

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

載荷識別； 不適定問題； 正則化方法； 二次規劃； 形函數

0　引言

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

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

1　正則化方法

離散系統方程的建立基于Duhamel積分：

(1)

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

(2)

取Ω內中間時間單元近似載荷作為真實載荷的最優最小二乘近似，設S=2l，定義首個時間單元中形函數為

(3)

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

(4)

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

(5)

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

(6)

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

再令

(7)

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

(8)

則每個形函數對應的矩陣為

(9)

l-1 columns

(10)

(11)

(12)

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

y=Gf,

(13)

式(13)中系統矩陣往往是病態的。

正則化方法能夠有效克服系統矩陣的病態導致的不適定問題，主要思想為考察最優化問題：

(14)

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

(15)

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

(16)

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

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

2　最優正則化參數選取方法

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

(17)

商函數方法實現步驟：

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

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

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

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

3　數值仿真

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

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

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

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

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

(18)

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

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

(19)

3.1.1SISO系統非共振區載荷識別

圖2　SISO系統在2%噪聲水平下商函數Fig.2 QF under 2% noise level of SISO system

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

3.1.2SISO系統共振區載荷識別

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

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

圖4　SISO系統在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

噪聲水平/%商函數方法GCV方法最優參數RErr/%最優參數RErr/%20.01757.700.10208.7750.078610.150.102010.88

表2SISO系統共振區載荷識別結果

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

噪聲水平/%商函數方法GCV方法最優參數RErr/%最優參數RErr/%20.011111.930.374012.8050.027713.910.371016.25

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

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

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

3.2.1MIMO系統非共振區載荷識別

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

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

表3　MIMO系統非共振區載荷識別結果

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

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

3.2.2MIMO系統共振區載荷識別

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

表4　MIMO系統共振區載荷識別結果

4　試驗驗證

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

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

圖8　鋼制懸臂梁結構Fig.8 The cantilever beam structure

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

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

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

5　結論

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

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

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

[1]文祥榮,智浩,孫守光.結構動態載荷識別的精細逐步積分法[J].工程力學,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]王曉燕,黃維平,李華軍.地震動反演及結構參數識別的EKF算法[J].工程力學,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].工程力學,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；編輯：任志平

國家自然科學基金項目(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