周昊天，趙 婕，曹為政，史曉鳴，于開平

(1.哈爾濱工業(yè)大學(xué) 航天學(xué)院，黑龍江 哈爾濱 150001； 2.黑龍江科技大學(xué) 計(jì)算機(jī)與信息工程學(xué)院，黑龍江 哈爾濱 150022； 3.上海機(jī)電工程研究所，上海 201109)

0 引言

近年來，高超聲速飛行器因其在速度上的巨大優(yōu)勢(shì)而受到研究和應(yīng)用領(lǐng)域的廣泛關(guān)注[1]，以高超聲速飛行器為背景的研究成為各領(lǐng)域的研究熱點(diǎn)。高超聲速飛行器的高速特性也給各領(lǐng)域的研究帶來了諸多挑戰(zhàn)[2]。如何辨識(shí)在氣動(dòng)加熱環(huán)境下結(jié)構(gòu)的模態(tài)參數(shù)是結(jié)構(gòu)動(dòng)力學(xué)領(lǐng)域遇到的挑戰(zhàn)之一。處于飛行狀態(tài)的飛行器模態(tài)參數(shù)的辨識(shí)問題一直備受關(guān)注[3]，因無法獲知其所受外載荷的詳細(xì)信息，故一般采用運(yùn)行模態(tài)分析(OMA)方法。此外，由于飛行器的熱環(huán)境隨飛行狀態(tài)的改變而改變[4]，結(jié)構(gòu)的動(dòng)力學(xué)特性也隨時(shí)間的變化而變化，因而需要用時(shí)變模態(tài)參數(shù)辨識(shí)方法解決這一問題。

子空間跟蹤方法是眾多時(shí)變模態(tài)參數(shù)辨識(shí)方法中常用的一類方法。文獻(xiàn)[5]提出了基于快速近似子空間跟蹤(FAST)的時(shí)變模態(tài)參數(shù)辨識(shí)方法，對(duì)一三連桿系統(tǒng)進(jìn)行了辨識(shí)；文獻(xiàn)[6-7]引入了投影估計(jì)子空間跟蹤(PAST)算法，對(duì)二連桿和移動(dòng)質(zhì)量-簡(jiǎn)支梁進(jìn)行了辨識(shí)；文獻(xiàn)[8]引入了自然冪迭代(NPI)和估計(jì)冪迭代(API)方法，對(duì)彈簧-質(zhì)量?jī)勺杂啥认到y(tǒng)和移動(dòng)質(zhì)量-懸臂梁系統(tǒng)進(jìn)行了辨識(shí)；文獻(xiàn)[9-10]利用改進(jìn)的API方法對(duì)一二連桿機(jī)構(gòu)和一個(gè)衛(wèi)星模型進(jìn)行了辨識(shí)；此外，還有基于新信息準(zhǔn)則(NIC)的子空間跟蹤算法[11]。近年來，也有一些學(xué)者進(jìn)行了熱環(huán)境下飛行器動(dòng)力學(xué)特性的試驗(yàn)研究。余艷輝等[12]研究了金屬熱防護(hù)結(jié)構(gòu)在熱載荷作用下的動(dòng)力學(xué)特性；蘇華昌等[13]提出了改進(jìn)的熱模態(tài)試驗(yàn)方法，利用振動(dòng)臺(tái)對(duì)試件施加猝發(fā)快掃形式的基礎(chǔ)激勵(lì)，并結(jié)合激光測(cè)振儀進(jìn)行非接觸式測(cè)量，最終獲得舵面一階頻率隨時(shí)間和溫度的變化規(guī)律；楊凱等[14]設(shè)計(jì)了一套溫度可控的實(shí)驗(yàn)設(shè)備，采用時(shí)頻分析方法對(duì)時(shí)變結(jié)構(gòu)進(jìn)行了辨識(shí)；李宇峰等[15]研究了氣動(dòng)加熱的模擬方法，獲得了不同溫度梯度條件下熱防護(hù)壁板結(jié)構(gòu)的振型、固有頻率等模態(tài)信息。

傳統(tǒng)的熱模態(tài)分析采用頻域方法，在某一短時(shí)間段內(nèi)將結(jié)構(gòu)視為時(shí)不變系統(tǒng)，并在此時(shí)間段內(nèi)進(jìn)行辨識(shí)，無法形成完整連續(xù)的辨識(shí)結(jié)果。針對(duì)此問題，本文總結(jié)了子空間跟蹤辨識(shí)模態(tài)參數(shù)的一般步驟，利用一個(gè)時(shí)變兩自由度系統(tǒng)模型的響應(yīng)數(shù)據(jù)，通過添加不同量級(jí)噪聲來對(duì)系統(tǒng)進(jìn)行再辨識(shí)，考察了3種應(yīng)用較為廣泛的子空間跟蹤算法和1種較新的基于新信息準(zhǔn)則的子空間跟蹤算法，并從中選出了抗噪能力最強(qiáng)的算法。最后，設(shè)計(jì)了飛行器舵面結(jié)構(gòu)的加熱試驗(yàn)，利用優(yōu)選出的子空間跟蹤算法對(duì)結(jié)構(gòu)升溫階段的模態(tài)參數(shù)進(jìn)行了辨識(shí)，得出了結(jié)構(gòu)完整的連續(xù)時(shí)間固有頻率變化規(guī)律。

1 方法過程簡(jiǎn)介

考慮離散時(shí)變狀態(tài)方程

(1)

(2)

式中：k為時(shí)刻；M為選定的Hankel矩陣塊行數(shù)。

由此得到狀態(tài)量z(k)的更新方式為

(3)

式中：上標(biāo)“?”表示求對(duì)應(yīng)矩陣的偽逆。且

(4)

(5)

式(3)、(4)中Y(k)U?(k)和[U(k)UT(k)]-1的更新方式為

[Y(k)U?(k)]=[Y(k-1)U?(k-1)]-

z(k)wT(k)

(6)

[U(k)UT(k)]-1=[U(k-1)UT(k-1)]-1-

(7)

式(3)～(7)組成了狀態(tài)量z(k)的更新規(guī)則。對(duì)任意k時(shí)刻的輸入矩陣U(k)組建其對(duì)應(yīng)的投影矩陣U⊥(k)，使得U(k)U⊥(k)=0，則有

U⊥(k)=I-U?(k)U(k)=

I-UT(k)[U(k)UT(k)]-1U(k)

(8)

且

Y(k)U⊥(k)=

Y(k)-Y(k)UT(k)(U(k)UT(k))-1U(k)=

Γ(k)X(k)U⊥(k)

(9)

對(duì)Y(k)U⊥(k)進(jìn)行奇異值分解(SVD)，可得

Y(k)U⊥(k)=R(k)Σ(k)V(k)T

(10)

式中：由R(k)的前n列組成的矩陣即為k時(shí)刻的廣義能觀矩陣Γ(k)。

tr(WT(k)W(k))

(11)

式中：Rzz(k)=E[z(k)zT(k)]。此處：“E”表示期望；“tr”表示矩陣的跡。子空間跟蹤方法優(yōu)化問題具體求解過程可參考文獻(xiàn)[5-10]。

得到任一k時(shí)刻的信號(hào)子空間后，可通過式(12)估計(jì)離散系統(tǒng)矩陣，具體為

(12)

式中：W1(k)和W2(k)分別由W(k)的前(M-1)×m和后(M-1)×m行組成。對(duì)估得的離散系統(tǒng)矩陣進(jìn)行特征值分解，可得

(13)

式中：Λ(k)=diag(λ1(k),λ2(k),…,λn(k))表示由共軛復(fù)特征值組成的對(duì)角陣；ψ(k)為特征向量。k時(shí)刻第i階系統(tǒng)偽固有頻率ωi(k)和模態(tài)阻尼比ξi(k)可通過式(14)、(15)求得，即

ωi(k)=ln(λi(k))/(2Δtπ)

(14)

(15)

式中：上標(biāo)“R”表示實(shí)部；Δt為數(shù)據(jù)采樣周期。

綜上所述，基于子空間跟蹤的時(shí)變模態(tài)參數(shù)辨識(shí)方法的一般步驟可總結(jié)如下：

1)對(duì)激勵(lì)和響應(yīng)數(shù)據(jù)進(jìn)行預(yù)處理，使數(shù)據(jù)滿足遞推步驟的需要，得到狀態(tài)向量z(k)；

3)利用復(fù)模態(tài)分析理論得出時(shí)變頻率ωi(k)、時(shí)變模態(tài)阻尼比ξi(k)等模態(tài)參數(shù)；

4)重復(fù)步驟1)～3)，直至所有響應(yīng)數(shù)據(jù)被處理。

2 仿真計(jì)算

本節(jié)對(duì)一個(gè)有理論參考解的兩自由度彈簧-質(zhì)量系統(tǒng)進(jìn)行辨識(shí)，系統(tǒng)簡(jiǎn)圖如圖1所示。

圖1 兩自由度彈簧-質(zhì)量模型示意圖Fig.1 Schematic diagram of two degree-of-freedom spring-mass model

系統(tǒng)的運(yùn)動(dòng)控制方程為

(16)

(17)

系統(tǒng)剛度矩陣K(t)取值為

(18)

系統(tǒng)阻尼矩陣C(t)為

0.5M+10-4K(t)

(19)

為模擬舵面受熱時(shí)的剛度變化情況，令式(18)中剛度K1(t)=1 000 exp(-0.4t)，K2=2 000，仿真計(jì)算的時(shí)間步長(zhǎng)Δt=1 ms，時(shí)長(zhǎng)總計(jì)為5 s。辨識(shí)方法采用文獻(xiàn)[2-7]中所述方法。為模擬實(shí)測(cè)結(jié)構(gòu)響應(yīng)，在仿真響應(yīng)數(shù)據(jù)中添加一定量的噪聲，使最終的響應(yīng)數(shù)據(jù)信噪比(SNR)為10。各方法采用相同的計(jì)算參數(shù)，最終辨識(shí)結(jié)果如圖2所示。

圖2 信噪比為10時(shí)各算法辨識(shí)結(jié)果Fig.2 Identification results with differentalgorithms when SNR equals to 10

圖中：低于5 Hz的頻率為第一階固有頻率；10 Hz左右的頻率為第二階固有頻率；實(shí)線部分為仿真模型的理論參考值；點(diǎn)線部分為各算法辨識(shí)值。從圖中可以看出，各算法都能很好地辨識(shí)出兩階固有頻率，且能對(duì)固有頻率進(jìn)行跟蹤。其中：NPI算法的第二階辨識(shí)結(jié)果在某些時(shí)刻偏離理論參考值較遠(yuǎn)；API算法和PAST算法的辨識(shí)結(jié)果非常接近；NIC算法的辨識(shí)結(jié)果更為接近理論參考值，總體上看精度更高。

為對(duì)不同算法在不同量級(jí)噪聲情況下的辨識(shí)結(jié)果進(jìn)行量化評(píng)估，引入平均絕對(duì)誤差率(MAPE)，其定義為

(20)

通過以上4種算法并利用信噪比分別為25，15，5，2時(shí)的響應(yīng)數(shù)據(jù)進(jìn)行辨識(shí)，最終得到的辨識(shí)結(jié)果MAPE值如表1所示。

表1 不同算法辨識(shí)結(jié)果MAPE值比較

從表1中可以看出：在響應(yīng)信號(hào)質(zhì)量較好時(shí)，各算法均能很好地辨識(shí)結(jié)果，各階次相對(duì)理論參考值的誤差較??；隨著信噪比的降低，各算法辨識(shí)結(jié)果的誤差隨之增大，本例中，NPI算法和API算法在信噪比為5時(shí)的第二階辨識(shí)結(jié)果MAPE值超過了50%，說明辨識(shí)結(jié)果已不可信。從辨識(shí)結(jié)果的精度來看，在諸多算法中，PAST算法和NIC算法表現(xiàn)較好，在信號(hào)質(zhì)量較差時(shí)也能對(duì)結(jié)構(gòu)的固有頻率做出估計(jì)。在信噪比為5時(shí)，各算法辨識(shí)結(jié)果如圖3所示。

由圖可知：NPI算法辨識(shí)結(jié)果的第二階誤差較大，在2～3 s時(shí)間段內(nèi)失去了頻率跟蹤能力；從總體上看，相較于其他幾種算法，NIC算法的辨識(shí)結(jié)果更接近理論參考值，可認(rèn)為NIC算法的抗噪聲能力更強(qiáng)。

3 飛行器典型結(jié)構(gòu)熱模態(tài)辨識(shí)

飛行器中存在大量的梁、板類結(jié)構(gòu)，如飛行器舵面。在飛行器飛行過程中，舵面結(jié)構(gòu)所在的熱力學(xué)環(huán)境也會(huì)隨之變化，由此導(dǎo)致的結(jié)構(gòu)動(dòng)力學(xué)特性的變化一直是工業(yè)領(lǐng)域的關(guān)注熱點(diǎn)。因此，本節(jié)對(duì)一個(gè)鈦合金舵面結(jié)構(gòu)進(jìn)行辨識(shí)，該舵面結(jié)構(gòu)的三維尺寸如圖4所示。

圖4 舵面結(jié)構(gòu)尺寸以及傳感器布置圖(mm)Fig.4 Profile of rudder structure and sensor distribution(mm)

舵面結(jié)構(gòu)為鈦合金TA15材質(zhì)，舵面根部突出部分用于夾持實(shí)驗(yàn)夾具。圖中：“*”號(hào)標(biāo)記處為高溫傳感器布置點(diǎn)，其中，舵面左上方第二個(gè)位置處為激勵(lì)點(diǎn)，其余為加速度響應(yīng)測(cè)量點(diǎn)。

高溫實(shí)驗(yàn)裝置原理如圖5所示。圖中：(1)為與試驗(yàn)件相匹配的夾具；(2)為水冷組件；(3)為紅外石英燈陣；(4)為試驗(yàn)件；(5)為高溫加速度傳感器；(6)為測(cè)溫?zé)犭娕迹?7)為激振桿；(8)為力傳感器；(9)為激振器。數(shù)據(jù)采集系統(tǒng)集成了信號(hào)發(fā)生器的功能。激振器通過功率放大器與數(shù)據(jù)采集器相連，其余傳感器通過電荷放大器與數(shù)據(jù)采集器相連。溫控PLC內(nèi)置了溫控程序，結(jié)合測(cè)溫?zé)犭娕及凑赵O(shè)定好的溫度曲線對(duì)石英燈陣進(jìn)行功率控制，從而實(shí)現(xiàn)控制試驗(yàn)件溫度的功能。

圖5 實(shí)驗(yàn)裝置原理圖Fig.5 Schematic diagram of experiment setup

實(shí)驗(yàn)過程中對(duì)試件施加均方根值為1 N的白噪聲激勵(lì)。未進(jìn)行加熱前首先采集30 s的結(jié)構(gòu)振動(dòng)響應(yīng)數(shù)據(jù)；而后通過石英燈陣加熱使舵面溫度在24 s內(nèi)由常溫升高至500 ℃，此為快速升溫階段；當(dāng)舵面溫度達(dá)到500 ℃后保持恒溫狀態(tài)100 s；之后關(guān)閉加熱系統(tǒng)，此為自然降溫過程。試件溫度變化曲線如圖6所示。在以上加熱、恒溫和降溫過程中，數(shù)據(jù)采集系統(tǒng)同時(shí)對(duì)舵面結(jié)構(gòu)進(jìn)行激勵(lì)，并采集振動(dòng)響應(yīng)數(shù)據(jù)。

圖6 試件溫度變化曲線Fig.6 Temperature curve of test structure

從上一節(jié)對(duì)仿真模型的辨識(shí)結(jié)果來看，NIC算法的抗噪聲能力更強(qiáng)，因此采用該算法對(duì)舵面結(jié)構(gòu)進(jìn)行辨識(shí)。由于結(jié)構(gòu)表面溫度迅速升至500 ℃，而結(jié)構(gòu)自身參數(shù)的變化要滯后于這樣的快速變化，因而選取前50 s進(jìn)行辨識(shí)。利用NIC算法辨識(shí)出的前四階結(jié)果如圖7所示。從圖中可以看出，隨著溫度升高，各階固有頻率均出現(xiàn)了不同程度的下降。第一階固有頻率在初始時(shí)刻和跟蹤結(jié)束時(shí)擾動(dòng)較大，這主要與本試驗(yàn)中低頻加速度信號(hào)幅值較小有關(guān)；更高階頻率的辨識(shí)值相對(duì)穩(wěn)定，能反映結(jié)構(gòu)固有頻率隨溫度變化情況。實(shí)驗(yàn)結(jié)果表明：NIC算法可用于熱模態(tài)分析，對(duì)加速度響應(yīng)幅值較高的頻段辨識(shí)效果較好；基于NIC子空間跟蹤的模態(tài)參數(shù)辨識(shí)方法可用于有激勵(lì)信號(hào)的情況，當(dāng)激勵(lì)信號(hào)為白噪聲時(shí)，該方法可轉(zhuǎn)換為僅利用響應(yīng)數(shù)據(jù)進(jìn)行辨識(shí)的運(yùn)行模態(tài)分析方法。

圖7 前四階辨識(shí)結(jié)果Fig.7 Identification results of first four orders

4 結(jié)束語

本文總結(jié)了子空間算法的一般步驟，結(jié)合一個(gè)有理論參考解的時(shí)變兩自由度彈簧-質(zhì)量模型考察了4種子空間算法的辨識(shí)精度，結(jié)果發(fā)現(xiàn)NIC算法的抗噪能力較強(qiáng)。利用能夠?qū)︼w行器舵面結(jié)構(gòu)進(jìn)行受控加熱的實(shí)驗(yàn)裝置采集了舵面結(jié)構(gòu)在快速升溫過程中的加速度響應(yīng)數(shù)據(jù)，并僅利用響應(yīng)數(shù)據(jù)進(jìn)行了時(shí)變模態(tài)參數(shù)辨識(shí)。辨識(shí)結(jié)果表明：基于NIC子空間跟蹤的模態(tài)參數(shù)辨識(shí)方法適用于熱模態(tài)分析，對(duì)振動(dòng)響應(yīng)幅值較高的頻段辨識(shí)效果較好，具有一定的工程應(yīng)用價(jià)值。但受子空間類算法本身所限，模型定階、計(jì)算參數(shù)選擇等問題仍待解決。