Tikhonov方法在不適定模型修正中的應用

邱飛力，張立民，張衛華

（西南交通大學牽引動力國家重點實驗室，成都 610031）

Tikhonov方法在不適定模型修正中的應用

邱飛力，張立民，張衛華

（西南交通大學牽引動力國家重點實驗室，成都 610031）

數值建模和分析在結構動態設計中應用廣泛，為獲取準確的計算模型，基于參數靈敏度有限元修正技術得到迅速發展。然而，參數靈敏度矩陣病態性和修正目標函數方程組的不適定性，造成最小二乘直接法難以得到穩定的物理解。為此，對靈敏度矩陣的病態和模型優化方程組的不適定性進行研究，以鐵道車輛6自由度離散模型和有限元支架模型為載體，采用Tikhonov正則化方法分別完成了超定系統和欠定系統仿真模型的參數修正。從而，解決了病態不適定系統中最小二乘直接法無穩定物理解的缺陷。參數修正后的模型準確反應了實際結構的尺寸差別和質量，這表明該方法具有實際的應用價值。

模型修正；參數靈敏度；不適定系統；正則化方法

有限元結構修正方法隨著有限元理論及應用而飛速發展，其目的在于縮小有限元模型與實際結構的參數誤差，使得有限元模型能夠準確反應結構的特性［1－2］。目前，基于參數靈敏度分析的有限元修正技術得到了長足的發展和廣泛的應用［3－5］。模型修正通常轉換為線性系統的求解過程，該線性系統可能出現不適定特性，同時可能含有病態的靈敏度分析矩陣，利用最小二乘法難以獲取具有物理意義的參數值［6－7］。為避開靈敏度分析，對于理論模型可以采用Berman矩陣修正方法，僅涉及簡單矩陣運算無需迭代優化［8－9］。文獻［10］中對離散的理論模型進行了矩陣修正，但修正后的矩陣改變了原始質量和剛度矩陣的帶狀和稀疏性，物理意義不明顯。對于試驗模型，基于響應面方法概率統計理論的響應面方法，同樣可以避開求取靈敏度矩陣，然而多參數響應面方法的樣本空間數較大。文獻［11－12］中分別采用直接靈敏度方法和響應面法對欠定試驗模型進行了修正，結果表明修正后的響應面方法精度較靈敏度低且樣本空間數目大，但直接靈敏度方法獲取的修正參數改變過大，或超出材料屬性范圍。

針對待修正系統的解不唯一和解不連續依賴于系數矩陣等問題，引入正則化方法對待修正的線性系統進行約束（初始狀態、邊界條件等），平衡靈敏度矩陣的病態特性和模型修正系統不適定性，獲取模型修正的穩定解［13］。與Berman矩陣修正方法、響應面修正法、直接靈敏度修正方法比較，引入Tikhonov正則化方法進行模型修正，能夠保持矩陣帶狀和稀疏性，多參數參與修正也無需進行大樣本空間點計算，同時能夠確保修正后的參數具有明顯的物理意義。

1 正則化方法

1.1 系統不適定性

正則化方法出現后，不適定問題得到了廣大學者的研究［14］。不適定修正方程如（1）式所示。

式中：A∈Rm×n為參數靈敏度矩陣，b0∈Rm為修正目標殘差向量；ε為b0的擾動組成的誤差向量。

Hansen指出線性不適定系統特點為［15］：系數矩陣A秩虧損、奇異值逐步趨于零、條件數K（A）比較大。采用傳統的方法，如Gauss消元、LU分解、Cholesky分解、QR分解法等計算，A或b中存在的微小誤差會使得解x發生振蕩，與真實解相差極大，導致解喪失物理意義。

1.2 不適定問題可解性

對結構的靈敏度矩陣A作奇異值分解［16］：

式中：矩陣U＝（u1，u2，…，un）和Vnn＝（v1，v2，…，vn）的列向量正交。對角矩陣Λnn的對角元素滿足σ1＞σ2＞…＞σn＞0。靈敏度矩陣A的偽逆矩陣如式（3）。

Hansen給出了離散的Picard條件［17－18］：線性系統A x＝b的離散傅里葉系數趨于零的速度快于矩陣A的奇異值趨于零的速度，則表明該線性系統滿足離散Picard條件。離散Picard條件決定了系統是否可以通過正則化方法獲取具有物理意義的解。

離散傅里葉級數系數的計算式如式（5）所示。

1.3 正則化

20世紀60年代，Philips［19－20］提出了正則化方法。Tikhonov等［21－22］發表了正則化方法的專著，為求解不適定問題奠定了基礎。Tikhonov方法是目前應用最為廣泛，也是最著名的正則化方法，它將最小二乘法轉換為式子（6）的形式。

式中：a是正則化參數，C是正則化矩陣，d是對解x的限制向量（初始狀態、邊界條件等）。

式（6）可以改寫式（7）所示。

殘差向量J1（a），J2（a）的二范數，決定了正則化參數的取值和應用。

L－曲線法是利用對數尺度來描述殘差范數和解的限制范數的曲線對比，該方法的特征是對數尺度圖形中出現明顯的L形狀曲線，曲線拐點所對應的正則化參數作為優化參數［23］。

Hansen對L－曲線在正則化方法的應用作了深入研究，將L－曲線在對數尺度下最大曲率點作為其拐點［24］。

式中：k（a）取最大值時，對應的a為最優的正則化參數。

2 修正中的不適定問題

根據修正的目標模態階數、修正參數不同，靈敏度矩陣的維數不同。線性求解系統A x＝b，（A∈Rm×n）可以分為兩類：超定系統m＞n；欠定系統m＜n。

2.1 超定理論模型

以鐵道車輛的6自由度模型為對象，進行超定系統的模型修正分析，其組成見圖1［25］。

圖1 鐵道車輛6自由度模型Fig．1 The six-dofmodel of the vehicle system

鐵道車輛系統參數見表1和表2。

表1 車輛系統質量及轉動慣量參數表Tab.1 Themass andinertia moment of vehicle system

表2 車輛系統尺寸及懸掛剛度參數表Tab.2 The vehicle length and suspension stiffness

一系懸掛剛度k4和k6存在誤差，初始值為k4＝214 N／m，k6＝214 N／m，其理論值為k4＝2 140 040 N／m，k6＝2 140 040 N／m。

建立車輛系統垂向振動和垂平面轉動的六自由度微分方程，如式（15）所示。

矩陣［M］為對角矩陣，其主對角線元素分別為：m1，m2，m3，I1，I2，I3；剛度矩陣［K］如式（16）所示。

初始模型和理論模型的模態分析結果（見表3）。

表3 模態參數對比表（單位：Hz）Tab.3 Thecom parison of themodal parameters

從表3可知，初始參數車輛模型各階頻率均存在誤差，其最大誤差為98．0%。

選取1階、2階、4階作為修正目標，選取靈敏度方法進行模型修正，建立目標函數如式（17）。

式中：f，f0分別為計算模態頻率和理論模態頻率向量；b為計算模型與理論模型頻率殘差向量，ε為極小正數。

首先求解待修正模態頻率關于k1、k4的靈敏度，振動系統特征值對剛度參數ki的一階靈敏度公式為［26］：

對式（22）直接求解，修正后模態誤差見表4。

表4 優化后模態參數表Tab.4 The updated m odal parameters

優化后的剛度參數值為：k4＝214 N／m，k6＝245 634 714 N／m。由此可知，采用最小二乘直接法無法獲取穩定的物理解。

根據式（2）計算可知，［Aji］秩虧矩陣，奇異值σi逐步趨于零，其條件數為1．73×1014，式（22）具有不適定性。結合式（5）計算離散傅里葉系數 ri，繪制離散picard圖（見圖2）。

圖2 靈敏度矩陣離散Picard圖Fig．2 The picardcurve of the sensitivitymatrix

從圖2可知，靈敏度矩陣的傅里葉系數衰減速度快于奇異值，表明式（22）滿足離散picard條件。

按照式（6）構造目標優化式子，分析可知，系統剛度參數具有對稱性，存在k1＝k4關系，故取［C］＝｛1，－1｝，｛d｝＝｛0｝。

利用Tikhonov方法進行正則化參數的選取，根據式（12）繪制L－曲線見圖3。

圖3 L－曲線圖Fig．3 The L-curve

利用式（14）可以求出L－曲線的最大拐點位置，見圖3中小圓圈。根據式（10）獲取的最優正則化參數為a＝7．5×10－7。

將a代入式（6），對正則化約束的目標函數進行最小二乘求解，修正后各階模態對比見表5。

表5 修正后模態參數表Tab.5 The updated modal parameters

由表5知，修正后各階模態誤差較小，最大誤差僅為0．15%。

剛度參數隨迭代變化見圖4，修正后的一系懸掛剛度為k4＝k6＝2 150 214 N／m。與理論值比較，修正后剛度參數的相對誤差僅為0．48%。

圖4 剛度參數隨迭代次數變化Fig．4 The stiffness changing with iteration

應用Tikhonov方法修正后的不適定車輛模型，各階模態值與理論模態高度吻合，且保證了剛度參數k4、k6的物理意義。

2.2 欠定試驗模型

某支架結構模型（見圖5）。由三種類型鋼材組成，槽鋼、方鋼和角鋼，分別對應圖中編號1、2、3。

利用沖擊激勵法［27］，對支架進行實驗模態分析，其頻響函數（見圖6）。

由圖6可知，支架結構的主要模態頻率為30．10 Hz、61．84 Hz，分別對應支架的1階、2階彈性模態。

圖5 支架幾何結構及型材編號圖Fig．5 The geometry and component numbers of frame

圖6 支架頻響函數圖Fig．6 The frequencyresponse function of the frame

支架結構有限元模型初始參數（見表6）。

表6 支架模型初始參數表Tab.6 Theinitial parameters of the framemodel

應用有限元軟件ANSYS進行模態分析，提取0～100 Hz內的彈性模態，將其與試驗模態對比見表7。

表7 有限元與實驗模態對表Tab.7 The comparison ofmodal parameters

由表7可知，有限元計算的1階、2階彈性模態誤差較大，需要進行修正。考慮到支架結構材料一致，型材厚度不同，故將彈性模態EX、密度ρ、槽鋼厚度R1、方鋼厚度R2、角鋼厚度R3共5個參數作為修正參數。

利用ANSYS軟件計算1階、2階彈性模態的差分靈敏度［28］矩陣［A］，結合式（21），采用最小二乘法直接求解，1階、2階彈性模態頻率隨迭代變化見圖7。

圖7 彈性模態隨迭代振蕩圖Fig．7 The elastic modal changing with iteration

由圖7可知，模態頻率隨迭代呈現振蕩現象。且優化中密度參數取值低至3 900 kg／m3，彈性模量取值低至103 GPa。優化參數值遠超出材料屬性誤差范圍，已喪失了物理意義。

根據式（2）求得靈敏度矩陣的奇異值σi，結合式（5）計算系統離散傅里葉系數ri，其離散Picard圖（見圖8），由圖8可知系統滿足離散Picard條件。

圖8 靈敏度度矩陣離散picard曲線Fig．8 The picardcurve of the sensitivitymatrix

為約束修正參數穩定性，C矩陣選取單元矩陣。結合表7有限元模態誤差，建立優化表達式（23）。min（‖［A］2×5｛x｝5×1－｛b｝2×1‖+a‖C x‖） （23）根據式（12）繪制L－曲線（見圖9）。

圖9 L－曲線圖Fig．9 The L-curve

圖9中小圓圈即為L－曲線最大拐點位置，正則化參數取值為a＝1．55×10－4。將a代入式（23）進行模型修正，1階、2階彈性模態隨迭代次數變化（見圖10和圖11）。

圖10 第1階彈性模態頻率迭代變化圖Fig．10 The first ordermodal frequency changing with iteration

表8 修正后模態與實驗模態對表Tab.8 Thecom parison of the updated and experim entalmodal

圖11 第2階彈性模態頻率迭代變化圖Fig．11 The second ordermodal frequency changing with iteration

修正后模型模態參數值（見表8）。

分析表8，修正后的1階和2階彈性模態頻率誤差極小，與實驗模態值保持高度一致。

修正各參數值見表9。

表9 修正后各參數值Tab.9 The parameters value of updated model

從表9可以看出，彈性模量EX、密度變化量比較小，參數值具有物理意義。方鋼厚度為0．003 8 m，較角鋼、槽鋼厚，這與支架實際結構十分吻合。

稱取支架實際結構質量，并計算修正前后支架有限元模型質量，對比見表10。

表10 支架質量對比表Tab.10 The com parison of the framemass

由表10可知，修正后支架質量更接近真實值（誤差4%），驗證了修正后模型的合理性和準確性。

3 結 論

應用正則化方法修正不適定性的計算模型，得出以下結論：

（1）直接采用最小二乘法進行不適定模型修正，修正后的鐵道離散車輛模型剛度參數k4、k6與理論值相差甚遠；支架有限元模型修正中出現振蕩不收斂，且彈性模量EX、密度ρ超出了材料屬性范圍，喪失物理意義。因此，在不適定模型修正中，不建議直接采用最小二乘法。

（2）鐵道車輛超定理論模型和支架欠定試驗模型修正中，引入了較小的正則化參數a和對應的正則化矩陣C對模態頻率目標函數進行約束，從而獲取了穩定且具有物理意義的剛度k、彈性模量EX、密度ρ值，解決了傳統的直接修正法在不適定模型修正中的振蕩和無物理解的缺陷。

（3）超定理論模型和欠定的試驗模型，分別經過12次和30次迭代修正，修正后的模態頻率誤差極小（＜1%），說明Tikhonov方法精度高且收斂快。

（4）對比修正前后的欠定試驗模型，修正后的模型準確地反應了型材的結構差別（方鋼較厚）；同時，修正后模型質量更接近于真實結構質量。這表明修正后參數真實可靠，正則化方法在實測結構的模型修正中具有可行性和適用性。

［1］吳曉菊．有限元結構模型修正綜述［J］．特種結構，2009，26（1）：39－45．

WU Xiao-ju．The summary of the finite element structure updating［J］．Special Structure，2009，26（1）：39－45．

［2］芮強，王紅巖，歐陽華江．機械結構動力學模型修正技術的現狀和發展［J］．裝甲兵工程學報，2012，26（2）：1－8．

RUIQiang，WANG Hong-yan，OUYANG Hua-jiang．Status and development of model updating in mechanical structural dynamics［J］．Journal of Academy of Armored Force Engineering，2012，26（2）：1－8．

［3］姚春柱，王紅巖，滕永超，等．基于靈敏度分析的車輛點焊結構有限元模型修正研究［J］．中國工程機械學報，2012，10（3）：325－328．

YAO Chun-zhu，WANG Hong-yan，TENG Yong-chao，et al．Sensitivity analysis-based FEMrevision for vehicle spotwelding structures［J］．Chinese Journal of Construction Machinery，2012，10（3）：325－328．

［4］劉小川，張凌霞，牟讓科．基于靈敏度分析和響應面方法的有限元模型修正［J］．振動工程學報，2008，21（S）：102－105．

LIU Xiao-chuan，ZHANG Ling-xia，MOU Rang-ke．Finite element model updating based on sensitivity analysis and response surface method［J］．Journal of Vibration Engineering，2008，21（S）：102－105．

［5］林恰．基于敏感性分析的懸索橋有限元模型修正［D］．成都：西南交通大學，2007．

［6］李紹新，許曉安，徐軍發，等．基于L曲線法正則化參數優化反演光子相關光譜［J］．華中科技大學學報：自然科學版，2010，38（6）：102－106．

LI Shao-xin，XU Xiao-an，XU Jun-fa，et al．Inversion of photon correlation spectroscopy by means of L-curve optimizing regular parameters［J］．J．Huazhong Univ．of Sci．＆T ech：Natural Science Edition，2010，38（6）：102－106．

［7］Miller G F．Fredholm equation of equation of the first kind in numerical solution of integral equations［M］．Oxford Garden Press，1974．

［8］祁泉泉，辛克貴．改進的Berman-Baruch法在非比例阻尼有限元模型修正中的應用［J］．工程力學，2012，7（29）：1－5．

QI Quan-quan，XIN Ke-gui．Application of the improved Berman-Bruch method for model updating in the nonproportional damping structure［J］．Engineering Mechanics，2012，7（29）：1－5．

［9］李龍梅，賈昌樂．有限元動力模型的修正——改進Berman修正法［J］．煙臺大學學報：自然科學，1994，4：35－41．

LILong-mei，JIA Chang-le．Finite-elementmodel updating—an improve Berman method［J］．Journal of Yantai University：Natural Science and Engineering，1994，4：35－41．

［10］邱飛力．基于優化算法的車輛懸掛系統物理參數識別研究［D］．成都：西南交通大學．2009．

［11］邱飛力，張立民，張衛華，等．支架結構建模中設計參數的修正與優化［J］．噪聲與振動控制，2014，34（1）：36－40．

QIU Fei-li，ZHANG Li-min，ZHANG Wei-hua．Updating and optimization of design parameters in frame structuremodeling［J］．Noise and Vibration Control，2014，34（1）：36－40．

［12］邱飛力，張立民，張衛華，等．基于響應面方法的支架結構模型修正研究［J］．噪聲與振動控制，2014，34（3）：139－143．

QIU Fei-li，ZHANG Li-min，ZHANGWei-hua，et al．Frame model updating study based on response surface method［J］．Noise and Vibration Control，2014，34（3）：139－143．

［13］Ahmadian H，Mottershead J E，Fris-well MI．Regulation methods for finite element model updating［J］．Mechanical Systems and Signal Processing，1998，12（1）：47－64．

［14］Hansen P C．Analysis of discrete ill-posed problems bymean of the L-curves［J］．Society for Industrial and Applied Mathematics，1992，34（4）：561－580．

［15］Hansen P C．The discrete picard condition for ill-posed problems［J］BIT，1990，30：658－672．

［16］鄒紅星，王殿軍，戴瓊海，等．延拓矩陣的奇異值分解［J］．電子學報，2001，23（3）：289－292．

ZOU Hong-xing，WANG Dian-jun，DAIQiong-hai．Singular value decomposition for extendedmatrix［J］．Acta Electronica Sinica，2001，23（3）：289－292．

［17］Hansen P C．Rank-defcient and discrete Ill-posed problems［M］．SIAM，Philadelphia，1998．

［18］CalvettiD，Reichel L，Sgallari F，etal．A regularizing lanczos iteration method for underdetermined linear systems，Journal of Computational and Applied Mathematics［M］．Namerical aspects of linear inversion，2000，115：101－120．

［19］Phillips D L．A technique for the numerical solution of certain integral equations of the firstkind［J］．J．ACM，1962，9：84～97．

［20］Tikhonov A N．Solution of incorrectly formulated problems and the regularization method［J］．Soviet ath，Dokl，1963，4：1035－1038．

［21］吉洪諾夫．不適定問題的解法［M］．北京：地質出版社，1979．

［22］Tikhonov A N，Arsenin V Y，Solution of ill-posed problems［M］．Wiley，New York，1977．

［23］霍智勇，朱秀昌．用于圖象復原的離散L－曲線自適應算法［J］．南京郵電大學學報：自然科學版，2008，28（4）：27－33．

HUO Zhi-yong，ZHU Xiu-chang．An adaptive discrete L-curve algorithm for image restoration［J］．Journal of Nanjing University of Posts and Telecommunications：Natural Science，2008，28（4）：27－33．［24］Hansen P C，Leary D P．The use of L-curves in the regularizations of discrete ill-posed problems［J］．Society For Industry And Applied Mathmatics，1993，14（6）：1487－1503．

［25］張立民，董鐵軍．鐵道車輛懸掛系統振動特征頻率靈敏度分析［J］．振動與沖擊，2009，128（1）：28－31．

ZHANG Li-min，DONG Tie-Jun．Analysis of eigenvalue sensitivity for a railway vehicle suspension system［J］．Journal of vibration and shock，2009，128（1）：28－31．

［26］王北京，王紅巖．車輛典型部件結構的有限元模型修正方法［J］．裝甲兵工程學報，2011，25（6）：40－44．

WANG Bei-jing，WANG Hong-yan．Updatingmethod of finite elementmodel for structures of typical components of vehicle［J］．Journal of Academy of Armored Force Engineering， 2011，25（6）：40－44．

［27］于金朋，邱飛力，馬紀軍，等．基于沖擊激勵的支架結構參數研究［J］．北京交通大學學報，2011，35（6）：98－101．

YU Jing-peng，QIU Fei-li，MA Ji-jun，etal．Research on stents structure parameters by impactive incentive［J］．Juornal of Beijing Jiaotong University，2011，35（6）：98－101．

［28］韓建平，董小軍，錢炯．基于環境振動測試的鋼筋混凝土渡槽空腹桁架拱有限元模型修正［J］．蘭州理工大學學報，2009，35（5）：122－124．

HAN Jiang-ping，DONG Xiao-jun，QIAN Jiong．Finite element model updating of arched Vierendeel truss of RC aqueduct based on ambient vibration test［J］．Journal of Lanzhou University of Technology，2009，35（5）：122－124．

Parameters updating of simulation models with ill-posed characteristics

QIU Fei-li，ZHANG Li-min，ZHANGWei-hua
（The National State Key laboratory，Southwest Jiaotong University，Chengdu 610031，China）

Along with the wide application of numerical analysis and modeling，getting a correct simulation model becomes an urgent requirement and consequently the parameter-sensitivity updatingmethod has been developed rapidly．The direct least square method can't always get the steady physical solution in the cases of ill-posed target function equations and ill-conditioned sensitivity matrixes．The ill characteristics of the sensitivity matrixes and target function equationswere investigated．A six-DOF discrete vehicle and a frame finite elementmodelwere updated with the Tikhonov regulation method by using the over determined and under determined simulation model respectively．The defect of the direct least squaremethod was solved．The updated models reflect exactly the realmass and the size difference．of the structure．It's proved that themethod is applicable in engineering practice．

model updating；parameter-sensitivity；ill-posed system；regulation method

TH113．1；O327

A

10．13465／j．cnki．jvs．2015．12．021

2014－04－15 修改稿收到日期：2014－05－26

邱飛力男，博士生，1987年生

張立民男，研究員，1960年生

郵箱：zhang＿lm01＠163．com