脈沖風洞測力系統建模與載荷辨識方法研究

脈沖風洞測力系統建模與載荷辨識方法研究

王鋒1,2，賀偉1，毛鵬飛1，張小慶1

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

摘要：將載荷辨識技術應用于脈沖燃燒風洞模型測力。用子結構綜合法建立了測力試驗系統的動力學模型，在時域內將動力學方程進行離散，建立起天平測量信號與模型氣動載荷歷程之間的線性關系，作為載荷辨識的模型。采用Tikhonov正則化和子空間投影法相結合的混合正則化方法，將高維的、不適定的載荷辨識問題轉化為低維的適定問題，以利于快速求解。提出了一種新方法來確定合適的投影子空間維數，然后應用L曲線準則來尋找低維正則化問題的最優正則化參數。最后通過算例驗證了系統建模方法的精度和載荷辨識算法的有效性與穩定性。

關鍵詞：脈沖風洞;動力學建模;載荷辨識;正則化;反問題

中圖分類號:V231;O32文獻標志碼：A

Dynamicmodelingoftestingsysteminimpulsefacilitiesandloadidentificationmethod

WANG Feng1,2, HE Wei1,MAOPeng-fei1,ZHANGXiao-qing1(1.Air-breathingHypersonicTechnologyResearchCenter,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China;2.ScienceandTechnologyonScramjetLaboratory,ChinaAerodynamicsResearchandDevelopmentCenter,Mianyang621000,China)

Abstract:Load identification method was applied to force measurement in impulse combustion facilities. The dynamic model of the testing system was built using substructure synthesis method. Then, the differential equations were discretized to establish the linear relationship between the outputs of balance and the dynamic forces history acting on the aircraft model. A hybrid regularization method which combines the Tikhonov regularization and subspace projection method was proposed to solve the load identification problem. The hybrid method converts the ill-conditioned large scale problem into a well-posed small scale problem which can be solved easily. A new method was proposed to determine the size of the projection subspace. The L-curve criteria was employed to search the optimal regularization parameter of the dimension-reduced problem. The accuracy of the structural dynamic modeling method and the effective and stability of the load identification method were validated by a numerical example.

Keywords:impulsefacilities;structuraldynamicmodeling;loadidentification;regularization;inverseproblem

吸氣式高超聲速飛行器是當前國際航空航天技術的熱點研究領域，地面風洞試驗是最重要的研究手段，脈沖燃燒風洞是國內目前開展大尺度超燃發動機和機體-推進一體化高超聲速飛行器試驗研究的主力設備之一，具有口徑大、建造和運行成本較低的優點。由于其有效的試驗時間較短，只有約200ms～300ms，給模型測力試驗帶來一定的困難。目前采用常規的測力方法，為了在短時間內測得模型的氣動力，需盡量降低模型的質量，同時盡量提高天平的剛度，以便提高試驗系統的響應速度。吸氣式飛行器采用機身和推進系統一體化設計，內外流緊密耦合，為了反映真實的氣動-推進特性，希望采用全尺寸或接近全尺寸的試驗模型。所以，目前在風洞開展試驗的飛行器模型尺寸和質量越來越大，導致試驗系統的響應速度降低，彈性振動也難以衰減，基于穩態測量思想的常規測力方法面臨挑戰。因此，有必要研究新的測力方法，提高脈沖風洞測力的準確性，更好地發揮其優勢。

本文研究采用載荷辨識技術進行間接測力的方法。載荷辨識是根據測量得到的結構響應(應變、位移、速度或加速度等)，結合結構的動力學方程(或傳遞函數)反求結構所受的載荷歷程。作為結構動力學的第二類反問題，載荷辨識技術獲得了廣泛的研究和應用[1-5]，其難點在于問題本身的不適定性，因此在求解時普遍使用正則化方法，文獻[6]對此進行了較全面的總結。

載荷辨識應用于模型測力，不需要結構響應達到穩態，可以放松對天平剛度的要求和對模型質量的限制，同時，辨識得到的是載荷隨時間的變化過程，能反映更豐富的試驗動態信息，如發動機點火時的推力建立過程，這對于精細地研究高超聲速飛行器的氣動-推力特性具有重要價值。Simmons等[7]最早將載荷辨識技術用于激波風洞模型阻力測量，使用的應力波天平，后來又將該技術發展用于多分量同時測量[8-10]。考慮到脈沖燃燒風洞的試驗時間遠大于激波風洞(5ms左右)，本文的方法在試驗系統建模、天平標定方面均和基于應力波天平的測力技術有所不同，并且提出了新的載荷辨識算法。

1試驗系統的動力學建模

載荷辨識需已知系統的動態特性，所以首先建立試驗系統的結構動力學方程，并希望方程的自由度盡量少，以降低后續的計算量。

圖1　試驗系統示意圖 Fig.1 Schematic of the testing system

1.1試驗系統結構組成

風洞測力試驗系統主要由飛行器模型、天平、支架以及數據采集系統組成，對于結構動力學建模，僅考慮前三部分。由于飛行器模型一般比較長、比較重，同時，為了避開其腹部的推進流道，通常采用背部支撐的方式固定在風洞中，系統的剖面如圖1所示。天平的浮動框與試驗模型固連，天平的固定框與支架上端固連，支架下端固定于地面。在來流作用下，試驗模型帶動天平浮動框運動，使天平的彈性元件產生應變，利用應變片產生輸出信號。

1.2建模方法

可以用有限元方法直接對整個試驗系統進行建模，然后通過模態降階得到低階模型，但是，為了在動力學方程中保留模型的剛體自由度和直接使用天平的標定數據，在此采用子結構方法來建模。將系統分為試驗飛行器、天平和支架三個子結構。由于試驗時天平的浮動框與試驗飛行器固連在一起，固定框與支架固連在一起，因此，在建模時把浮動框歸入試驗飛行器，把固定框歸入支架，而天平子系統本身只考慮連接浮動框與固定框的彈性測量元件。

子結構采用Craig-Bampton方法[11](后面簡稱C-B方法)建模。C-B方法的子結構模態包括界面約束模態(后面簡稱約束模態)和固定界面正規模態(后面簡稱正規模態)，因此，首先要確定子結構間的連接界面，并根據需要作必要的簡化。

前面已將天平浮動框歸入試驗飛行器、將天平固定框歸入支架，所以試驗飛行器與天平的連接面就是天平測量元件與浮動框的對接面，天平與支架的連接面就是測量元件與固定框的連接面。實際中，天平的浮動框和固定框的剛度遠遠大于測量元件的剛度，相對于彈性測量元件，兩個框可近似看作剛體，于是它們與測量元件的對接界面均可簡化為一個6自由度的點。將試驗飛行器(含浮動框)的界面點記作點Ⅰ、天平(測量元件)在浮動框一側的界面點記作點Ⅱ、天平(測量元件)在固定框一側的界面點記作點Ⅲ、支架的界面點記作點Ⅳ。

1.2.1試驗飛行器的模型

試驗飛行器用有限元方法建模。為了將其與天平的對接面簡化為一個點，在適當的位置，例如質心處，設置一節點作為主節點(即點Ⅰ)，將浮動框上與天平彈性元件對接面上的節點作為從節點，用剛性單元相連。在此基礎上，用成熟的商業有限元軟件可以直接進行子結構分析，輸出子結構的質量和剛度矩陣，從而建立其動力學方程

(1)

這里，下標B(Boundary)表示界面約束模態，下標N(Normal)表示正規模態。ηⅠ,B是試驗飛行器約束模態的位移向量，為6個剛體位移

(2)

ηⅠ,N為正規模態位移向量，假設取前nm階正規模態，則

(3)

因此，飛行器子系統共有6+nm個自由度。因為約束模態均為剛體模態，所以子陣KⅠ,BB為零矩陣，而子陣MⅠ,BB為飛行器的剛體質量矩陣(參考點為點Ⅰ)，而fⅠ,B正是飛行器所受氣動力和力矩向量(關于點Ⅰ的)。fⅠ,N是正規模態對應的廣義力，其第Ⅰ個分量的計算式如下

(4)

1.2.2天平子系統的模型

天平的浮動框和固定框分別歸屬飛行器和支架后，天平子系統僅剩下彈性測量元件，由于測量元件尺寸相對較小，質量也較輕，為簡化模型，忽略其慣性，從而將天平子系統建模為無質量的廣義彈簧，其與試驗飛行器和支架分別有一個對接點(點Ⅱ和點Ⅲ)，所以天平子系統共12個自由度，其模型為12×12的剛度矩陣。簡化模型示意圖見圖2。

雖然可以用有限元計算天平的剛度矩陣，但天平是試驗系統的關鍵環節，實際中要經過校準，顯然利用校準數據來得到其剛度矩陣將更為準確。

圖2　天平子系統的力學模型 Fig.2 Mechanical model of the balance

點Ⅱ、點Ⅲ的位移向量分別記作

天平子系統的力學模型為

(5)

式中，fⅡ、fⅢ分別是作用于點Ⅱ、點Ⅲ的載荷向量，每個載荷向量包括3個力和3個力矩。在校準時，是將天平固定框完全約束，通過剛性足夠的夾具與浮動框連接，然后在夾具上施加載荷，載荷施加點稱為校心，一般位于天平的幾何中心線上。在此，假設點Ⅱ即為校心。點Ⅲ選在同一中心線上位于固定框的一側，與固定框剛性連接，與點Ⅱ之間的距離為h。那么，校準時的邊界條件即為：點Ⅲ完全固定，在點Ⅱ施加校準載荷。于是在方程(5)中，ηⅢ=0，從而有

KⅡ,ⅡηⅡ=fⅡ

(6)

KⅢ,ⅡηⅡ=fⅢ

(7)

代入式(6)有

KⅡ,ⅡRⅡ=PⅡ

(8)

于是得到

(10)

利用式(7)和式(8)可得

KⅢ,ⅡRⅡ=-TB1PⅡ=-TB1KⅡ,ⅡRⅡ

(11)

所以有

KⅢ,Ⅱ=-TB1KⅡ,Ⅱ

(12)

因為剛度矩陣是對稱的，所以有

KⅡ,Ⅲ=KTⅢ,Ⅱ=-KTⅡ,ⅡTTB1=-KⅡ,ⅡTTB1

(13)

現在尚有KⅢ,Ⅲ待求。假設解除點Ⅲ的所有約束，允許天平產生剛體運動，在點Ⅱ沿6個自由度依次施加微位移Δx、Δy、Δz、Δθx、Δθy、Δθz，寫成位移向量矩陣

由于現在點Ⅲ自由，它要產生相協調的一組位移向量，記作ΔRⅢ，忽略二階小量，易得

(14)

因為點Ⅲ不受力，將ΔRⅡ、ΔRⅢ代入式(5)的第二組方程，有

KⅢ,ⅡΔRⅡ+KⅢ,ⅢΔRⅢ=0

(15)

再利用式(14)可得

(16)

KⅢ,Ⅲ=TB1KⅡ,ⅡTTB1

(17)

至此，利用天平校準時的加載及位移測量數據可以計算出其剛度矩陣。天平校準時還要獲取天平方程，即天平傳感器(一般是應變計)輸出與校心載荷的關系，這是常規天平校準的首要目的，具體過程在此不詳述，假設已獲得如下天平方程

v=NB,ⅡfⅡ

(18)

式中：v是天平輸出(一般是電壓)，其與載荷向量是等長的，系數矩陣NB,Ⅱ是方陣，且非奇異。現在需要推導天平輸出與點Ⅱ、Ⅲ位移的關系。因為校準時固定框(點Ⅲ)被完全約束，所以ηⅢ=0，式(6)成立，將其代入式(18)得

v=NB,ⅡKⅡ,ⅡηⅡ?SⅡηⅡ

(19)

這里，將SⅡ=NB,ⅡKⅡ,Ⅱ稱為點Ⅱ對應的傳感矩陣。尚需求點Ⅲ對應的傳感矩陣，將其記作SⅢ，則兩點都產生位移時，天平的輸出應該是

v=SⅡηⅡ+SⅢηⅢ

(20)

如果給定點Ⅱ的位移ηⅡ，而讓點Ⅲ自由運動，則根據式(14)有ηⅢ=TB2ηⅡ，此時天平僅有剛體位移，輸出應該為零，根據式(20)有

SⅡηⅡ+SⅢTB2ηⅡ=0

(21)

所以有

(22)

代入式(20)，有

此即天平子系統的傳感方程。

1.2.3支架子系統模型

支架子系統亦用有限元方法建模，其過程與試驗飛行器子系統一樣，需設置一節點作為主節點(即點Ⅳ)，選取固定框上與彈性元件對接處的所有節點為從節點，用剛性單元相連。在此基礎上，用商業有限元軟件執行子結構分析可以獲取子結構的質量和剛度矩

陣，得到其動力學方程為

(24)

支架約束模態的位移向量為

假設正規模態取前ns階，記為

則方程(24)共6+ns個自由度。

1.2.4系統模型綜合

圖3　點Ⅰ和點Ⅱ的相對位置 Fig.3 Position relationship of point Ⅰ and point Ⅱ

盡管各子結構之間的界面已簡化為一個點，但是各子結構的界面點是為建模方便而選取的，相鄰子結構的“界面”不一定重合，需要用剛性桿將其連接，以構成完整系統。

各子結構的界面點盡管不一致，但應該都是選在其對稱面(x-y平面，見圖1)內。假設點Ⅰ和點Ⅱ的相對位置如圖3所示，則用剛性桿連接后，其位移(微量)向量之間存在如下關系

ηⅡ=

(25)

類似地，可建立點Ⅲ和點Ⅳ位移向量的關系

(26)

利用式(25)、(26)，可將三個子結構的方程(1)、(5)和(24)綜合在一起，形成整個系統的動力學方程，其中獨立的自由度向量只有ηⅠ,B、ηⅠ,N、ηⅣ,B、ηⅣ,N，列在一起記作

系統的動力學方程為

(27)

其中

K=

實際中，天平置于飛行器內部，是不受外載荷作用的，因此fⅡ和fⅢ均為零。支架外部可以加整流外罩進行遮擋，也不受氣動載荷，因此fⅣ,B和fⅣ,N為零。非零載荷只有試驗飛行器的氣動載荷，其被分解為約束模態載荷和正規模態載荷，根據它們的定義，實際中前者是主要的，也是人們通常所指的氣動載荷，后者則相對較小，而且不直接引起天平輸出，因此可將fⅠ,N忽略，僅保留fⅠ,B。于是，試驗系統的力載荷向量簡化為

實際結構總是有阻尼的，考慮阻尼后系統動力學方程為

(28)

傳感方程(23)也需要用η來表示，易得

(29)

2載荷辨識模型

式(28)、(29)構成完整的試驗系統動力學輸入輸出方程。通過擴維降階將其變為狀態空間形式

(30)

y=Rx

(31)

式中：輸出y已知，輸入x待求。

3載荷辨識方法

3.1子空間投影法

載荷辨識是典型的反問題，方程(31)通常是病態的，輸出向量y中的噪聲會使得解x嚴重偏離真解，失去利用價值，而實際中測量噪聲(誤差)是不可避免的。為了得到有意義的解，必須對其進行正則化，Tikhonov正則化方法[15-16]是最常用的，它用如下最小化問題的解作為問題(31)的近似解[17]

(32)

式中α(>0)是正則化參數。上式與如下問題等價

(33)

實際中，由于結構頻率較高，測量的采樣率也很高，若待辨識的時段稍長則矩陣R的規模會很大，而在優化選擇正則化參數α的過程中，要反復求解式(33)，因此很難直接基于式(33)來求解。

為了使實際問題的求解更加可行，這里采用子空間投影法，將問題(32)的解限制在某個維數較低的解域上，從而降低求解規模。假設選定了一組基向量w1、w2、…、wk，將近似解x限制在由其張成的k維空間Wk上，即

x∈Wk=span{w1,w2,…,wk}

令

(34)

(35)

其等價問題為

(WTkRTRWk+αWTkWk)cα=(RWk)Ty

(36)

方程(36)的系數矩陣大小為k×k，通常k遠遠小于R的行數和列數，求解規模大大降低。

顯然，投影子空間的選擇，也就是投影基向量w1、w2、…、wk的選擇是影響近似解好壞的關鍵，為了逼近原問題的解，該子空間需盡量包含系數矩陣R前k個右奇異向量的信息。Krylov子空間具備這樣的特性[18]。k維Krylov子空間的定義如下[16, 18]

Krylov子空間的基向量可以按定義直接計算，但采用Golub-Kahan雙對角迭代法[19]可以得到性質更好的基向量，并為問題的求解帶來方便。基本算法如下：

w0=0，β1z1=y

對i=1,2,3,…,k，執行以下迭代：

αiwi=RTzi-βiwi-1

βi+1zi+1=Rwi-αkzi

式中，實數αi、βi的作用是分別使向量wi和zi單位化。迭代的結果得到如下關系式

RWk=Zk+1Bk

(37)

式中，

并且，Wk、Zk+1矩陣各自的列向量是相互正交的單位向量。將式(37)代入式(36)并利用Wk、Zk+1的正交性質，式(36)可簡化為

(BTkBk+αI)cα=BTkZTk+1y

(38)

根據前面的迭代算法，y=β1z1，所以BTkZTk+1y=β1BTkZTk+1z1=β1BTke1，其中e1標準歐氏空間的第1個向量。于是，式(38)進一步簡化為

(BTkBk+αI)cα=β1BTke1

(39)

注意BTkBk+αI為三對角矩陣，因此上式可以用追趕法快速求解。得到cα后，最終解由下式給出

(40)

3.2確定正則化參數的L曲線準則

(41)

(42)

式中zα是如下問題的解

(43)

當采用子空間投影后,有

(44)

(45)

為減小計算量，zα也用子空間投影法來近似計算，即令

zα≈Wkζ

(46)

從而將問題(43)變為如下容易求解的問題

(47)

于是，曲率的κ(α) 計算量大大降低，使得L曲線準則容易實施。

3.3子空間維數的確定方法

現在討論如何確定適當的子空間維數k。顯然，對同一個正則化參數α，取不同的子空間維數，得到的正則化解是不同的，因此，不同子空間維數下作出的L曲線是不重合的。但是，由于k維Krylov子空間的基向量具有逼近矩陣R前k個右奇異向量的性質，隨著k的增加，后面的基向量對正則解的貢獻逐漸減弱，所以L曲線本身隨k的增加將趨于收斂。在此，取如下指標作為判斷L曲線收斂的依據

J=‖xα‖2·‖rα‖2

(48)

3.4求解算法

根據前面的討論，對問題(31)完整的求解算法包括兩步，首先用3.3節的方法確定合適的子空間維數，然后在該子空間下用3.2節的L曲線準則確定最優的正則化參數，并得到對應的正則化解。第一步的詳細算法流程為：

(1)給定初始子空間維數k0和維數的增量Δk；給定i0和nα，計算出αi序列；給定參數ε；

(2)執行k0步Golub-Kahan雙對角迭代，得到式(37)的各矩陣；

(3)針對每個αi求解方程(39)，由(44)、(45)式得到‖xα‖2、‖rα‖2，計算指標函數值序列J(αi,k)，找出最小值對應的序號ik和最小值Jk；

(4)將當前子空間下的結果ik、Jk與上一步的結果ik-Δk、Jk-Δk作比較，若滿足條件

終止，得到子空間維數k；否則往下執行；

(5)k←k+Δk，繼續執行Δk步Golub-Kahan，擴充子空間，轉到第3步。

確定了子空間維數后，對應低維正則化問題的L曲線就確定了，現在需要尋找最大曲率點的正則化參數。這時，L曲線上各點的曲率κ是參數α的函數，即式(41)，可以采用一維搜索算法[23]求函數κ(α)的近似極大值解，起點可以取前面確定子空間維數時得到的αik，具體算法在此不再詳細列出。得到問題(39)的正則解cα后代入式(40)，最終得到問題(31)的正則化解。

4數值算例

考慮一個虛擬的試驗飛行器及其測力系統，其模型、天平及支架如圖4所示。飛行器長2.1m，材料為超硬鋁，重255.24kg，天平和支架為鋼材。

圖4　試驗系統三維圖 Fig.4 View of the testing system

用有限元法直接對全系統建模，稱其為全階模型，作為簡化模型的對比基準。天平也用有限元建模，通過靜力計算，模擬實際校準過程，用得到的數據建立簡化模型的剛度矩陣。每個測量通道取四個點的應變組成惠氏電橋，根據經驗設置應變片的靈敏度，將應變轉化為電壓，據此計算天平傳感矩陣。

4.1模態計算結果

用有限元軟件分別對試驗飛行器和支架進行子結構分析并輸出質量、剛度矩陣，得到子結構動力學模型。其中試驗飛行器取7個正規模態，支架取6個正規模態，加上各自的6個約束模態，系統簡化模型共25個自由度，而全階模型的節點數超過30萬，自由度超過100萬。對全階模型和簡化模型進行模態分析，得到的前10階固有頻率列于表1，可見簡化模型縱向和橫側向的前兩階頻率精度都是比較高的。

表1　固有頻率計算結果對比

4.2瞬態響應結果

某脈沖燃燒風洞典型的來流動壓變化曲線如圖5所示(對原始數據作低通濾波的結果)。現在計算試驗系統在其作用下的天平輸出。對于全階模型，用工程方法[24]計算試驗飛行器表面氣流壓力分布，模擬真實的受力狀態。對于簡化模型，先將飛行器的壓力分布系數積分，得到合力和力矩(關于質心)系數，再作為集中力加載。計算時，兩種模型的阻尼均以模態阻尼方式施加，各階阻尼比均取0.005。

假設側滑角為零，則氣動載荷只有阻力、升力和俯仰力矩。在圖5所示來流作用下，全階模型和簡化模型得到的天平三個通道的輸出對比如圖6所示，低頻信號吻合良好。高頻響應的誤差來自兩個方面，一是簡化模型本身的高頻特性誤差較大，二是兩種模型的加載方式不同，全階模型加載的是分布力，而簡化模型加載的是等效集中力，忽略了式(1)中的fⅠ,N。

由模態分析和瞬態響應結果的對比，說明低階簡化模型反映了試驗系統主要的動力學特性。

圖5　典型的風洞來流動壓曲線 Fig.5 The history of dynamic pressure of the flow

圖6　瞬態響應對比 Fig.6 Comparison of transient responses

圖7　載荷辨識結果對比 Fig.7 Comparison of the identified loads

4.3載荷辨識結果

以全階模型計算的瞬態響應當作試驗測量結果，對飛行器所受載荷歷程進行辨識。輸出采樣率取5 000Hz，采樣時長0.8s，有3個輸出，3個待辨識載荷，原始問題式(31)中矩陣R的大小為12 000×12 003。根據對真實試驗測量信號的統計分析，測量噪聲的標準差約為0.002 3mV，在此按0.005mV給計算輸出加白噪聲。計算時，初始子空間維數取550，增量取10，子空間收斂判定參數ε取0.001，正則化參數序列取10-30～10-20之間15個對數均布的點。最終確定的子空間維數為920，遠低于原始問題的維數，最優正則化參數為α=1.453 8×10-24。3個載荷的辨識結果與真實值的對比見圖7，二者吻合良好。定義辨識誤差為

(49)

對于辨識結果中的噪聲，可以通過適當的平滑處理來有效地降低。圖8是對阻力曲線用20點平均法進行平滑處理后的結果，噪聲大大降低，但載荷本身的低頻信息仍得以保留。

圖8　平滑后的阻力曲線 Fig.8 Smoothed drag force

現在考察測量信號中不同噪聲強度對辨識結果的影響。分別取噪聲的標準差為0.002mV、0.004mV、0.008mV、0.016mV、0.032mV和0.064mV時得到的辨識結果如表2所示。其中，辨識誤差取的是3個載荷的平均誤差。可見，隨著噪聲強度的增加，子空間維數在減小、正則化參數在增大，即對問題的正則化作用不斷加強，相應地，辨識誤差也在逐漸變大。均符合預期的規律。

表2　噪聲對辨識結果的影響

表3　不同計算條件下的結果

4.4算法的穩定性分析

輸出測量信號中的噪聲是隨機的，因此算法每次執行得到的結果亦可能有所不同，包括子空間大小、正則化參數以及辨識結果的精度。但是噪聲的統計特性基本不變，所以，對于一個穩定的算法，計算結果不應有大幅波動。另外，子空間維數增量Δk和正則化參數序列αi的設置也會對計算過程和結果產生影響。

現分別取Δk等于5、10、15，每種情況執行5次辨識算法，每次的噪聲隨機生成，正則化參數序列也在10-30～10-20之間隨機生成15個點。計算結果如表3所示，可見盡管子空間維數有所波動，但變化不大，而辨識結果的精度也變化很小，表明算法是穩定、魯棒的。

5結論

本文應用子結構綜合法，并充分利用天平校準數據，建立了試驗系統的低維動力學模型，可以比較準確地反映系統的低頻動態特性。通過對動力學方程進行時域離散，建立起試驗飛行器氣動力(力矩)與天平響應之間的線性關系方程，用來進行氣動載荷的辨識。

將Tikhonov正則化方法與基于Krylov子空間的投影法結合，把高維的不適定問題變為低維的適定問題，并充分利用Golub-Kahan雙對角迭代算法的性質，大大提高了求解效率。

提出了確定投影子空間大小的方法，該方法僅需計算低維問題的正則解及剩余向量的范數，計算簡單。子空間維數確定后，問題的規模大大減小，再應用L曲線準則確定最優的正則化參數，計算量大為降低，用一維搜索算法可快速得到最優參數。數值試驗表明，算法具有良好的穩定性。

在對試驗系統進行建模時，結構的阻尼參數是假定已知的，而實際中通常難以準確知道結構的阻尼參數，因此，利用測量響應，對載荷和阻尼參數同時進行辨識是值得進一步研究的課題。

參考文獻

[1]HollandsworthPE,BusbyHR.Impactforceidentificationusingthegeneralinversetechnique[J].InternationalJournalofImpactEngineering, 1989, 8(4): 315-322.

[2]ChoiK,ChangFK.Identificationofimpactforceandlocationusingdistributedsensors[J].AIAAJournal, 1996, 34(1): 136-142.

[3]MaoYM,GuoXL,ZhaoY.Astatespaceforceidentificationmethodbasedonmarkovparametersprecisecomputationandregularizationtechnique[J].JournalofSoundandVibration, 2010, 329(15): 3008-3019.

[4]HwangJS,KareemA,KimH.Windloadidentificationusingwindtunneltestdatabyinverseanalysis[J].JournalofWindEngineeringandIndustrialAero-dynamics, 2011, 99(1): 18-26.

[5]RonasiH,JohanssonH,LarssonF.Anumericalframeworkforloadidentificationandregularizationwithapplicationtorollingdiscproblem[J].Compu-ters&Structures, 2011, 89: 38-47.

[6]王林軍. 正則化方法及其在動態載荷辨識中的應用[D]. 長沙：湖南大學, 2011.

[7]SandersonSR,SimmonsJM.Dragbalanceforhypervelocityimpulsefacilities[J].AIAAJournal, 1991, 29(12): 2185-2191.

[8]SmithAL,MeeDJ,DanielWJT,etal.Design,modellingandanalysisofasixcomponentforcebalanceforhypervelocitywindtunneltesting[J].ComputersandStructures, 2001, 79: 1077-1088.

[9]RobinsonMJ,MeeDJ,TsaiCY,etal.Three-componentforcemeasurementsonalargescramjetinashocktunnel[J].JournalofSpacecraftandRockets, 2004, 41(3): 416-425.

[10]RobinsonMJ,HannemannK,SchrammJM.Designandimplementationofaninternalstresswaveforcebalanceinashocktunnel[J].CEASSpaceJournal, 2011, (1): 45-57.

[11]CraigRR,Jr,BamptonMCC.Couplingofsubstructuresfordynamicanalysis[J].AIAAJournal, 1968, 6: 1313-1319.

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

GUOXin-lin,MAOYu-ming,ZHAOYan,etal.Loadidentificationbasedonprecisetime-stepintegrationforMarkovparameters[J].JournalofVibrationandShock, 2009, 28(3): 27-31.

[13]ZhongWX,WilliamsFW.Aprecisetimeintegrationmethod[J].JournalofMechanicalEngineeringScience, 1994, 208: 427-430.

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

HANXu,LIUJie,LIWei-jie,etal.Acomputationalinversetechniqueforreconstructionofmultisourceloadsintimedomain[J].ChineseJournalofTheoreticalandAppliedMechanics, 2009, 41(4): 595-602.

[15]CalvettiD,MorigiS,ReichelL,etal.Tikhonovregulariza-tionandtheL-curveforlargediscreteill-posedproblems[J].JournalofComputationalandAppliedMathematics, 2000, 123(1-2): 423-446.

[16]LampeJ,ReichelL,VossH.Large-scaletikhonovregulari-zationviareductionbyorthogonalprojection[J].LinearAlgebraanditsApplications, 2012, 436(8): 2845-2865.

[17]BauerF,LukasMA.Comparingparameterchoicemethodsforregularizationofill-posedproblems[J].MathematicsandComputersinSimulation, 2011, 81(9): 1795-1841.

[18]HansenPC.Discreteinverseproblems:Insightandalgorithms[M].Philadelphia:SocietyforIndustrialandAppliedMathematics, 2010.

[19]GolubG,KahanW.Calculatingthesingularvaluesandpseudo-inverseofamatrix[J].JournaloftheSocietyforIndustrial&AppliedMathematics,SeriesB:NumericalAnalysis, 1965, 2(2): 205-224.

[20]HansenPC,O’learyD.Theuseofl-curveintheregulariza-tionofdiscreteill-posedproblems[J].SIAMJournalonScientificComputing, 1993, 14: 1487-1503.

[21]VogelCR.Computationalmethodsforinverseproblems[M].Philadelphia:SIAM, 2002.

[22]RegińskaT.Aregularizationparameterindiscreteill-posedproblems[J].SIAMJournalonScientificComputing, 1996, 17(3): 740-749.

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

[24]BonelliF,CutroneL,VottaR,etal.Preliminarydesignofahypersonicair-breathingvehicle[C]. 17thAIAAInternationalSpacePlanesandHypersonicSystemsandTechnologiesConference,SanFrancisco,California, 11-14April, 2011.