動態載荷辨識問題的正則化求解新方法*

王　鋒， 武　龍， 張小慶， 賀　偉

(1.中國空氣動力研究與發展中心超高速空氣動力研究所高超聲速沖壓發動機技術重點實驗室， 四川 綿陽 261000；2.中國空氣動力研究與發展中心吸氣式高超聲速技術研究中心， 四川 綿陽 261000)

?

動態載荷辨識問題的正則化求解新方法*

王鋒1， 武龍2， 張小慶2， 賀偉2

(1.中國空氣動力研究與發展中心超高速空氣動力研究所高超聲速沖壓發動機技術重點實驗室， 四川 綿陽 261000；2.中國空氣動力研究與發展中心吸氣式高超聲速技術研究中心， 四川 綿陽 261000)

Tikhonov正則化方法是求解載荷辨識問題的有效方法，正則化參數的確定是影響求解準確性的關鍵問題。提出了一種新方法——虛擬載荷法，以虛擬載荷向量的范數最小作為正則化參數確定準則，將載荷辨識問題轉化為單參數優化問題，用一維搜索法求得最優解。通過將法方程的系數矩陣預先進行三對角化，建立了高效的求解算法。通過兩個數值算例對新方法進行了檢驗，結果表明新的參數確定準則是有效的，基于該準則得到的載荷辨識結果具有滿意的精度。

載荷辨識； 反問題； Tikhonov正則化； 正則化參數

引　言

結構的動態載荷辨識是根據測量的結構響應(位移、速度、加速度或應變等)信息，結合已知的結構動力學特性，反求結構所受激勵的過程，屬于反問題的一種，有重要的科學研究價值，也有著廣泛的工程應用。載荷辨識方法大致可分為頻域法和時域法，前者研究較早，方法較成熟，主要通過對關心頻段上的系統頻響函數矩陣求逆來實現載荷識別[1-2]，比較適用于穩態或平穩隨機載荷的識別[3-4]。時域方法直接從結構動力學方程(微分或積分形式)出發，建立整個待辨識時段上載荷輸入與結構輸出的關系方程，即辨識方程，通過求逆運算，一次得到整個時段上的載荷。然而，辨識方程通常是病態的，解對輸出測量噪聲十分敏感，幾乎不可能用常規的矩陣求(廣義)逆的方法直接求解，必須進行適當的正則化處理，基于原來不適定的問題構建一個參數化的適定問題，通過適當選擇控制參數來使適定問題的解落在原問題解的一個可接受的鄰域內。針對科學和工程中遇到的各種反問題，人們已提出了多種正則化方法[5-7]，如Tikhonov正則化方法、奇異值截斷方法、迭代正則化方法等等，其中Tikhonov方法是比較常用的。對每種正則化方法，其最優正則化參數的選取準則都是影響辨識結果的關鍵因素，是非常重要的研究內容。

本文采用Tikhonov正則化方法求解載荷辨識問題，提出了一種新的準則來選取正則化參數，新準則基于力學概念提出，物理意義明確，而且計算簡單。

1　載荷辨識的數學模型

考慮含測量傳感器的線性結構動力學系統，其一般形式的動力學方程為

(1)

y=Rx

(2)

此即載荷辨識問題的數學模型。為了保證解的唯一性并有利于提高辨識精度，通常要求結構的測點數不少于待辨識的載荷數目。方程(2)通常是病態的，即使y中僅含有微小的誤差，直接求解y=Rx得到的解也會嚴重偏離真值，從而失去意義。而實際中，測量誤差是不可避免的。

2　Tikhonov正則化方法

Tikhonov正則化方法將問題(2)變為參數化的變分問題，用如下最小化問題的解作為其近似解[10-12]

(3)

式中α(>0)是正則化參數。等號右邊的第一項使近似解xα盡量滿足等式(2)，但又不能使其嚴格成立，因為y中噪聲會使解x嚴重偏離真解，呈現高頻振蕩的形態，第二項的作用就是為了抑制x中的振蕩，使解盡量光滑，參數α用來在兩者之間進行折中、平衡。對式(3)中的目標函數關于x求駐值，可得

(4)

式中I為單位矩陣。由于正則化參數α的取值決定著解的好壞，因此如何選取α是個重要的研究課題。很多研究者從不同的觀點出發，已經建立了一些選取準則和算法[7,13]，常用的有L曲線方法、交叉檢驗方法、偏差準則等，各種方法都有一定的優缺點和適用條件，Bauer和Lukas對此作了比較全面的對比分析[14]。各種方法都是要盡量挖掘原問題可供利用的先驗知識或設立某些假設，從而導出確定正則化參數的某種準則。數學上，人們希望正則化解在真解的某個鄰域內，但在缺少真解信息的情況下，其鄰域也就無從談起。現有準則通常主要是發掘利用輸出的有關信息。針對載荷辨識問題，下面給出一種新的準則，該準則利用解的信息，具有明確的物理意義，而且計算簡單。

(5)

該準則還可以按下式的意義來理解

(6)

即將虛擬載荷x2與真實載荷x1的能量之比和x1本身同時穩定在盡可能小的水平。

不過‖x‖2隨α是單調減小的[13]，當α→∞時，x(α)→0，從而‖x2‖2→0。于是，可以猜想，若α從0開始增加，‖x2‖2隨α的變化趨勢應該是先減小(從大的振蕩趨向于真值0向量)，后增大(背離真值0向量)，再減小(因α的過正則化)，所以最佳正則化參數應該位于使‖x2‖2從遞減轉為遞增的那個點。

3　求解算法

在求解問題(5)的過程中，需要反復解方程(4)，有必要對其進行預處理，以降低計算量。RTR為實對稱矩陣，存在正交矩陣H將其變為對稱的三對角矩陣，具體算法可采用Householder變換[15]

HT(RTR)H=Q

(7)

令

(8)

帶入式(4)，并用HT左乘式(4)，再利用式(7)并注意H為正交矩陣，可得

(9)

上式中系數矩陣Q+αI為對稱三對角矩陣，可以用追趕法快速求解，比直接求解式(4)的計算量大幅下降。

問題(5)是個單參數優化問題，采用一維搜索算法來求解[16]。首先尋找目標函數的單谷區間，然后在此區間上用黃金分割法搜索最優的α。根據前面的討論，因為當α從0開始增大時，‖x2‖2先遞減再遞增，所以，在搜索單谷區間時應從一個充分小的α開始，以避免直接進入后面的單調遞減區間，導致搜索失敗。此外，因為在最優的正則化參數附近，正則解對參數的變化十分敏感，也就是說‖x2‖2隨α的變化曲線在最優解附近可能十分陡峭且狹窄，為了易于捕捉到包含最優解的單谷區間，應該在對數尺度下進行搜索，即令α=10β，將β作為搜索變量，這樣，‖x2‖2隨β的變化曲線在最優解附近會變得“平坦”和“寬闊”。

因為R=[R1R2]，將方程(4)分解

(10)

而原問題對應的方程是

(11)

當x2不為零時，由式(10)解出的x1中包含了x2產生的附加干擾，為了得到更好的解，在確定了最優的α之后，由式(11)來求出x1作為最終解。

綜上所述，問題的求解過程如下：

(1)將RTR變換為三對角矩陣，得到H和Q，計算式(9)的右端項；

(2)令α=10β，以β為自變量，搜索‖x2‖2隨β變化曲線的單谷區間，β的初始值為一個適當的負數，在搜索過程中需針對不同的β求解方程(9)，再根據式(8)得到目標函數值‖x2‖2；

(3)確定單谷區間[β1，β2]后，用黃金分割法在該區間上搜索最優的β，得到對應的正則化參數α；

在今年年初的時候，大部分人都無法準確地預料得到今年的行情會如此之“慘”。如果是從春節作為行情劃分的調整時點算起，那么截止到2018年末，市場整整調整了11個月，這在A股歷史上都是較為“罕見”的。所以，2018年，帶給投資者的幾乎都是損失，只是每個人虧多虧少而已。所有的公募股票型主動管理基金，都是以虧損告終。那只在今年6月30日才由原來保本型貨幣基金轉型而來的股票型基金，因為一直空倉到年末，反而拿了公募股票型基金的冠軍。這足以說明A股市場2018年的慘淡。

(4)由方程(11)解出原問題的解x1。

4　算例分析

4.1平面桁架載荷辨識

考慮一平面桁架結構，如圖1所示，由79根長0.5 m、截面直徑15 mm的鋼質桿構成，在兩端簡支。假設在第31節點沿y方向作用一激勵載荷，如下

f(t)=100(1-cos2πω0t)sin6πω0t

(12)

其中，ω0=10 Hz。先積分求解系統的位移響應，然后對測點的位移計算結果疊加一定強度的高斯白噪聲,模擬實際測量結果，噪聲的均方差按下式確定

(13)

式中yj表示第j個測量輸出向量，l是其長度，ny是選取的測點個數，Δe是比例系數。

圖1　平面桁架結構Fig.1　Planar truss structure

在此，取第5和16節點的y向位移輸出來對f(t)進行辨識，而虛擬載荷假設沿y方向作用于第6個節點上，噪聲比例系數Δe取0.01。輸出采樣步長取0.001 s，時長0.5 s。載荷辨識結果與真實值的對比如圖2所示，二者吻合良好。圖3是計算得到的虛擬載荷歷程，其幅值相對于真實載荷較小，而且看上去類似零均值的噪聲，這是所期望的，說明真實載荷的信息沒有“泄露”到虛擬載荷上而影響到其辨識精度。

圖2　載荷辨識結果與真實值的對比Fig.2　Comparison between true load and identified load

圖3　虛擬載荷辨識結果Fig.3　Identified fictitious load

作為對比，不加虛擬載荷，而基于L曲線方法確定正則化參數，用Hanson開發的工具箱[17]來求解，得到最佳參數為6.8924×10-16，圖4中標出了其位置，比本文方法確定的參數值約小兩個量級，在該參數下得到的結果與真實載荷的對比如圖5所示，對噪聲的濾除效果稍差，辨識誤差比本文方法要大。

圖4　目標函數隨正則化參數的變化Fig.4　Curve of objective function vs. regularization parameter

圖5　基于L曲線方法的載荷辨識結果Fig.5　Load identification results based on L-Curve method

可以想象，取不同的虛擬載荷作用點，得到的辨識結果應該會有所不同。極端的情況是，虛擬載荷施加在真實載荷的作用點上并且方向相同時，兩者的作用效果相同，應該無法辨識出真實載荷。因此，虛擬載荷的施加點應避免與真實載荷作用點距離太近。為考察虛擬載荷作用點的影響，定義辨識誤差如下

(14)

圖6　虛擬載荷作用于不同節點時的辨識誤差Fig.6　Relative error history vs. action point of fictitious load

可見，當選擇距離實際載荷作用點較近的節點時，辨識誤差較大，如第30，32節點處約6%，第10，11節點處約9%，而選擇實際載荷作用點第31節點時，誤差為100%，辨識失敗。而其余節點對應的辨識誤差都差別不大，均在3.5%左右，這說明方法對虛擬載荷作用點的選擇并不敏感。即使將作用點選在較差的第10節點，辨識結果也是可以接受的，其與真實載荷的對比如圖7所示，甚至優于基于L曲線方法的結果。

圖7　虛擬載荷作用點選在節點10時的辨識結果Fig.7　Identification results when the 10th node is selected as action point of the fictitious load

4.2四邊固支板上的載荷辨識

考慮一矩形鋼板，長0.8 m，寬0.6 m，厚度2.5 mm，四邊固支。用有限元方法建模，單元劃分如圖8所示。假設板上同時作用2個載荷，選取4個節點的撓度測量值來對其進行辨識，使用1個虛擬載荷。圖8中，●表示真實載荷作用點，○表示虛擬載荷作用點，?為位移輸出測點。

圖8　矩形板網格及載荷作用點和測點分布Fig.8　Meshes of the plate, action points and measuring points

2個激勵載荷分別取間斷的三角形脈沖和連續的正弦信號，作用時長0.4 s。4個測點的位移響應如圖9所示，采樣步長取0.001 s。可見其主要由正弦激勵主導，兩次脈沖激勵僅引起相對較小的擾動。

圖9　4個測點的位移輸出Fig.9　History of outputs from the four measurement points

圖10　載荷辨識結果Fig.10　Load identification result

虛擬測量噪聲強度的定義仍為式(13)，當取Δe等于0.01時，2個載荷的辨識結果分別如圖10所示，與真實載荷吻合很好。此時，目標函數‖x2‖2隨α的變化曲線如圖11所示，變化規律依然符合預期。而采用L曲線準則確定的最佳正則化參數為8.06×10-15，標于圖11中，比本文方法確定的參數約小兩個數量級。基于該參數得到的辨識結果如圖12所示，與圖10相比，精度要差一些。

圖11　目標函數隨正則化參數的變化Fig.11　Curve of objective function vs. regularization parameter

圖12　基于L曲線方法的載荷辨識結果Fig.12　Load identification results based on L-Curve method

現考察噪聲強度對辨識結果的影響。取不同大小的噪聲，辨識結果的誤差和對應的最優正則化參數列于表1，隨著噪聲的增大，辨識誤差隨之增大，最優正則化參數也在增大，在對數標度下均與噪聲大小近似成線性關系，如圖13所示。

表1　噪聲水平對辨識結果的影響

圖13　辨識誤差和最優正則化參數隨噪聲水平的變化Fig.13　Variations of relative error and the optimal regularization parameter vs. level of noise

5　結　論

本文采用Tikhonov正則化方法求解線性結構的載荷辨識問題，通過增加一個虛擬載荷，使增廣問題的解部分已知，從而給出確定最優正則化參數的準則，該準則物理意義明確，計算簡單。在此基礎上給出了優化求解算法，算例表明，本文方法能夠穩健地辨識結構所受動態載荷，具有一定的應用價值。

不過，目前的工作仍屬于初步探索，尚需開展進一步的理論研究，為新方法的合理性提供支撐，至少包括兩方的深入工作：(1)從數學上證明虛擬載荷范數隨正則化參數變化曲線極小點的存在性和唯一性，或者給出極小點存在的條件，即什么情況下本文提出的準則是適用的；(2)建立虛擬載荷作用點優劣的判定函數，用于選擇最佳的虛擬載荷作用點。

本文的工作也為其他反問題求解提供了新的思路，即，設法將問題進行增廣，注入已知信息，然后在求解過程中充分利用已知信息來約束增廣問題的解，使其逼近真解。

[1]胡杰,張希農.一種頻域載荷識別的優化方法[J].噪聲與振動控制, 2009,(6):34—36.

HU Jie, ZHANG Xi-nong. An optimization method of load identification in frequency domain[J]. Noise and Vibration Control, 2009,(6):34—36.

[2]智浩,文祥榮,繆龍秀,等.動態載荷的頻域識別方法[J].北方交通大學學報,2000,24(4):5—10.

ZHI Hao, WEN Xiang-rong, MIAO Long-xiu, et al. Dynamic loading identification in frequency domain[J]. Jounal of North Jiaotong University, 2000,24(4):5—10.

[3]周林,鄭四發,王彬星,等.動態載荷識別位置優化的傳遞函數相干法[J].振動工程學報,2011,24(1):14—19.

ZHOU Lin, ZHENG Si-fa, WANG Bin-xing, et al. Coherence analysis method for dynamic force identif ication[J]. Journal of Vibration Engineering, 2011,24(1):14—19.

[4]朱濤,肖守訥,陽光武.載荷識別研究進展及其運用于鐵道輪-軌載荷研究概述[J].鐵道學報,2011,33(10):29—36.

ZHU Tao, XIAO Shou-ne, YANG Guang-wu. State-of-the-art development of load identification and it′s application in study on wheel-rail forces[J]. Journal of the China Railway Society, 2011,33(10):29—36.

[5]Sun X, Liu J, Han X, et al. A new improved regularization method for dynamic load identification[J]. Inverse Problems in Science and Engineering, 2014,22(7):1062—1076.

[6]Wang L, Han X, Liu J, et al. An improved iteration regularization method and application to reconstruction of dynamic loads on a plate[J]. Journal of Computational and Applied Mathematics, 2011,235:4083—4094.

[7]Vogel C R. Computational Methods for Inverse Problems[M]. Philadelphia: SIAM, 2002.

[8]郭杏林,毛玉明,趙巖,等.基于Markov參數精細積分法的載荷識別研究[J].振動與沖擊,2009,28(3):27—31.

GUO Xing-lin, MAO Yu-ming, ZHAO Yan, et al. Load identification based on precise time-step integration for Markov parameters[J]. Journal of Vibration and Shock, 2009,28(3):27—31.

[9]韓旭,劉杰,李偉杰,等.時域內多源動態載荷的一種計算反求技術[J].力學學報,2009,41(4):595—602.

HAN Xu, LIU Jie, LI Wei-jie, 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):595—602.

[10]Lampe J, Reichel L, Voss H. Large-scale Tikhonov regularization via reduction by orthogonal projection[J]. Linear Algebra and its Applications, 2012,436(8):2845—2865.

[11]Reichel L, Rodriguez G. Old and new parameter choice rules for discrete ill-posed problems[J]. Numerical Algorithms, 2012:1—23.

[12]梅立泉,崔維庚.面載荷識別的TSVD正則化方法[J].應用力學學報,2010,27(1):140—146.

MEI Li-quan, CUI Wei-geng. TSVD regularization method for area load reconstruction[J]. Chinese Journal of Applied Mechanics, 2010,27(1):140—146.

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

[14]Bauer F, Lukas M A. Comparing parameter choice methods for regularization of ill-posed problems[J]. Mathematics and Computers in Simulation, 2011,81(9):1795—1841.

[15]張賢達.矩陣分析與應用[M].北京:清華大學出版社,2004:248—250.

ZAHNG Xian-da. Matrix Analysis and Applications[M]. Beijing: Tsinghua University Press, 2004:248—250.

[16]袁亞湘,孫文瑜.最優化理論與方法[M].北京:科學出版社,1997:56—75.

YUAN Ya-xiang, SUN Wen-yu. Optimization Theory and Method[M]. Beijing: Science Press, 1997:56—75.

[17]Hansen P C. Regularization tools version 4.0 for Matlab 7.3[J]. Numerical Algorithms, 2007,46:189—194.

A technique for solving dynamical force identification problems by Tikhonov regularization method

WANGFeng1,WULong2,ZHANGXiao-qing2,HEWei2

(1.Science and Technology on Scramjet Laboratory, Hypervelocity Aerodynamics Institute of CARDC, Mianyang 621000, China；2.Air-breathing Hypersonic Technology Research Center of CARDC, Mianyang 621000, China)

Tikhonov regularization method is effective in solving load identification problems, and proper selection of the regularization parameter is a key point that determines the solution accuracy. In this work, a new method which is named fictitious load method is proposed. Taking the minimization of the norm of the fictitious load is as the criterion of regularization parameter selection, the load identification problem isthus converted into a one-variable optimization problem which can be solved by one dimension search technique. By transforming the coefficient matrix of the normal equation into tridiagonal matrix, an efficient algorithm is developed. Two numerical examples are presented to validate the new method, and the results show that the new regularization parameter selection criterion is effective to revover the loads on structures with good accuracy.

load identification; inverse problem; Tikhonov regularization; regularization parameter

2014-06-24；

2015-11-16

國家自然科學基金資助項目(11372339)；高超聲速沖壓發動機技術重點實驗室基金資助項目(STSKFKT2012001)

O327

A

1004-4523(2016)01-0031-07

10.16385/j.cnki.issn.1004-4523.2016.01.005

王鋒(1976—),男，博士，副研究員。電話： 15908208358； E-mail： fwang_cardc@163.com