基于泛函序列時變自回歸滑動平均模型的彈箭時變模態(tài)參數(shù)遞推估計方法

余磊， 張永勵， 袁夢笛， 劉瑞卿

(西安現(xiàn)代控制技術(shù)研究所, 陜西 西安 710065)

0 引言

隨著現(xiàn)代戰(zhàn)爭任務(wù)需求的多樣化，火箭、導(dǎo)彈的設(shè)計也朝著智能化、高速化和大型化方向發(fā)展，具有大長細(xì)比的柔性制導(dǎo)火箭、導(dǎo)彈不斷出現(xiàn)。柔性彈體的固有頻率較低，彈體低階彈性振動容易受到外界干擾激勵，給剛性彈體姿態(tài)運動的穩(wěn)定性帶來不利影響[1]。因此辨識出彈體的固有頻率有助于提高控制的精度并保證任務(wù)順利完成。

一般在飛行試驗前，通常采用地面試驗方式來確定彈體的固有頻率。然而，制導(dǎo)火箭、導(dǎo)彈在飛行過程中，隨著構(gòu)型變化、模塊化結(jié)構(gòu)調(diào)整和載荷工況變化等過程，使得彈箭成為變結(jié)構(gòu)、變質(zhì)量、變參數(shù)的不確定性對象，傳統(tǒng)的時不變結(jié)構(gòu)模態(tài)參數(shù)辨識方法不再適用，需要開展時變結(jié)構(gòu)模態(tài)參數(shù)辨識方法研究。

對于彈箭時變結(jié)構(gòu)這類復(fù)雜細(xì)長柔性體的模態(tài)參數(shù)辨識方法研究，主要集中于時域法和頻域法。由于時域方法辨識的直接性、物理概念的清晰性，工程應(yīng)用中采用較多。時域法又可分為基于子空間類方法和基于時間序列類方法，其中子空間法即將動力學(xué)參數(shù)表征到航天器系統(tǒng)的狀態(tài)方程中，并結(jié)合高精度陀螺儀或加速度計等敏感器的測量信息，采用特征值或奇異值分解技術(shù)，得到系統(tǒng)矩陣的最小實現(xiàn)。基于子空間的辨識方法得到了廣泛應(yīng)用[2-4]，但是由于其不可避免地要使用特征值或奇異值分解技術(shù)，對大型工程結(jié)構(gòu)而言數(shù)值復(fù)雜性較高，該方法仍有較大的發(fā)展空間。

基于時間序列類的方法以描述時間序列的數(shù)學(xué)模型為研究對象，其中以20世紀(jì)80年代在信號處理和控制領(lǐng)域發(fā)展而來的時變自回歸滑動平均(TARMA)模型[5]最受關(guān)注，并成功應(yīng)用于生物醫(yī)學(xué)信號分析[6]、風(fēng)速數(shù)據(jù)建模[7]、橋梁健康監(jiān)測[8]、齒輪箱故障診斷[9]及風(fēng)車振動信號分析[10]等諸多領(lǐng)域。Poulimenos等[11]和Spiridonkos等[12]總結(jié)了一系列基于TARMA模型的方法，并將這些方法按照時變系數(shù)的演化形式分為3類，包括無結(jié)構(gòu)化、隨機(jī)結(jié)構(gòu)化和確定性結(jié)構(gòu)化。其中，確定性結(jié)構(gòu)化的時變系數(shù)演化形式即是將時變序列模型的時變系數(shù)用一組基函數(shù)展開，進(jìn)而將時變系數(shù)轉(zhuǎn)換為投影空間中的時不變投影系數(shù)，這類方法也被稱為泛函序列(FS)-TARMA模型。Rao等[13]在1970年基于時間多項式，對時變系數(shù)進(jìn)行展開，這也是FS思想的首次提出。1974年，Mendel等[14]對一般的泛函時間序列模型進(jìn)行了研究，討論了時變系數(shù)的快變與慢變，為泛函時間序列模型提供了基礎(chǔ)。之后，多種多樣的基函數(shù)被用于描述FS-TARMA模型的時變系數(shù)，如Chebyshev多項式[15]、Jacobi多項式[16]、Fourier級數(shù)[17]、B樣條、無網(wǎng)格形函數(shù)[18]等。FS-TARMA模型辨識精度較高，但是在估計過程中需要將總時間段內(nèi)的時變參數(shù)矩陣作為整體進(jìn)行估計，矩陣維度大，一般采用批量處理，計算效率較低，無法較好地適應(yīng)彈體飛行模態(tài)辨識的快速響應(yīng)要求。

相對于確定性結(jié)構(gòu)化方法的“笨重”，無結(jié)構(gòu)化方法則利用“短時時不變假設(shè)”，即假設(shè)在一小段特征時間內(nèi)系統(tǒng)參數(shù)不隨時間變化，類似于時間凍結(jié)的思想[19]，通過選取合適的特征時間間隔，該假設(shè)能夠很好地應(yīng)用于一般的時變結(jié)構(gòu)系統(tǒng)。采用加入時間窗和遺忘因子等[20-21]手段，可以實現(xiàn)高效的模態(tài)參數(shù)更新，計算效率較高，但是計算精度一般低于批量算法。

本文針對彈箭飛行狀態(tài)下的時變模態(tài)參數(shù)辨識問題，基于確定性結(jié)構(gòu)化FS-TARMA模型，提出了一種遞推的小波基FS-TARMA時變結(jié)構(gòu)模態(tài)參數(shù)估計方法。該方法借鑒了無結(jié)構(gòu)化TARMA模型的遞推估計思想，有效地改善了FS-TARMA模型的計算效率，并保留了FS-TARMA模型計算精度高的優(yōu)點。通過一個大型運載火箭時變有限元模型驗證了本文提出的算法具有良好的時變模態(tài)辨識精度、時變跟蹤能力以及較高的計算效率。

1 小波基FS-TARMA模型

1.1 TARMA模型描述

通過振動信號辨識固有頻率屬于動力學(xué)第一類反問題，TARMA模型是基于非平穩(wěn)時間序列建立起來的離散模型，是解決該類反問題的有效工具。令na和nc分別表示自回歸(AR)與滑動平均(MA)的階數(shù)，則TARMA(na,nc)可表示為

(1)

式中：t為離散化時間；ti為第i階離散化時間；x[t]∈k×1為非平穩(wěn)振動響應(yīng)(多維)信號，k為信號維度；e[t]∈k×1為滿足零均值和協(xié)方差Σe[t]∈k×k的不相關(guān)殘差序列；ts為第s階離散化時間；Ati[t]和Cts[t]分別表示第i階和第s階TARMA模型的時變AR和MA系數(shù)矩陣。

TARMA模型“時間凍結(jié)”模態(tài)參數(shù)可通過在每一時刻求解如下特征值問題[22]獲得，

(pz[t]I-D[t])vz[t]=0,z=1,2,…,kna，

(2)

(3)

式中：I為單位矩陣。

更進(jìn)一步，在Δt的時間間隔下，可得TARMA模型“時間凍結(jié)”固有頻率和阻尼比分別為

(4)

1.2 基函數(shù)選擇

采用墨西哥帽小波函數(shù)作為基本小波，即母小波，對其進(jìn)行二進(jìn)尺度變換和時間軸上的平移。為方便運算，記小波基最大尺度因子等于時間基函數(shù)最大維數(shù)，可得其小波基為

(5)

式中：l為尺度因子；td為平移因子，td=0,1,…,2j-1；j為時間基函數(shù)的展開維數(shù)，j=0,1,…,m，m為最大維數(shù)。當(dāng)尺度因子l=0,1,…,m時，ρl,td(t)的伸縮即能擬合慢變和快變的時變系數(shù)，從而能夠覆蓋低頻到高頻的整個頻帶。則TARMA模型的時變AR系數(shù)矩陣和時變MA系數(shù)矩陣可以用墨西哥帽小波基函數(shù)展開表示：

(6)

式中：pa為時變AR系數(shù)擬合階數(shù)；pc為時變MA系數(shù)擬合階數(shù)。將(6)式代入TARMA模型，可得FS-TARMA模型如下：

(7)

式中：參數(shù)矩陣θT和回歸向量φ[t]具有如下形式：

(8)

(9)

式中：?為Kronecker積算子；基函數(shù)向量ρpa[t]、ρpc[t]為

ρpa[t]=[ρ0,0[t],ρ1,0[t],…,ρpa,2j-1[t]]，
ρpc[t]=[ρ0,0[t],ρ1,0[t],…,ρpc,2j-1[t]].

(10)

由(8)式可知，參數(shù)矩陣θT為常數(shù)矩陣。通過將TARMA模型參數(shù)矩陣表示為一組以墨西哥帽小波基函數(shù)的線性組合，可將時變的TARMA模型轉(zhuǎn)化為時不變的FS-TARMA模型，從而把時變模型參數(shù)估計問題轉(zhuǎn)化為時不變模型參數(shù)估計問題。估計得到的θT代入(6)式中，即可得到系統(tǒng)的模態(tài)頻率。

2 時變系數(shù)的遞推估計方法

根據(jù)參數(shù)估計原理，對(7)式中的參數(shù)矩陣θT采用最小二乘準(zhǔn)則進(jìn)行估計，得到估計的費用函數(shù)為

(11)

參數(shù)矩陣θT的估計為

(12)

式中：arg min (·)為最小化算子，即當(dāng)費用函數(shù)取最小值時，變量的取值；N為辨識數(shù)據(jù)長度。

不難發(fā)現(xiàn)，由于殘差序列中存在MA部分，即回歸向量φ[t]中包含了殘差的過去時刻值，估計參數(shù)θT的過程是一個非線性優(yōu)化問題。針對該問題，Poulimenos等[11]和Spiridonakos等[12]就曾指出，可以采用兩步最小二乘方法，將非線性優(yōu)化過程轉(zhuǎn)化為線性過程進(jìn)行求解。但是，該估計方法需要將待估參數(shù)θT向量化處理，使得計算過程中矩陣維度大，計算耗時長。

為了改善這一問題，可以借鑒無結(jié)構(gòu)化TARMA模型的遞推辨識思想。注意到不同時間長度內(nèi)的θT是不同的，因此θT是關(guān)于N的變量，即θT可視為變量θ[t]，利用遞推方式，從初始時刻一步步更新θ[t]，最終得到θT的估計值。定義最小二乘費用函數(shù)[24]如下：

(13)

式中：τ為離散時間變量。

通過對費用函數(shù)進(jìn)行求導(dǎo)，可得

(14)

式中：

(15)

基于(14)式，可以定義協(xié)方差矩陣P[t]的形式為

(16)

它的逆矩陣可以表示為

P-1[t]P-1[t-1]+φ[t]φT[t].

(17)

采用矩陣逆定理[25]，(17)式可進(jìn)一步表示為

(18)

(18)式代入(14)式中，可得

(19)

定義增益矩陣K[t]，信號的一步超前預(yù)測[t|t-1]以及一步前向預(yù)測誤差T[t,t-1]具有如下形式：

(20)

[t|t-1]=T[t-1]φ[t]，

(21)

[t,t-1]=x[t]-[t|t-1] .

(22)

則(19)式可以表達(dá)為

[t]T=T[t-1]+[t,t-1]KT[t].

(23)

綜上，給出小波基FS-TARMA模型的遞推估計格式如下：

(24)

遞推估計初始化條件為θ[0]T=0，P[0]=αI，α為模型結(jié)構(gòu)參數(shù)。若振動信號的數(shù)據(jù)長度為N，則參數(shù)矩陣θT的估計值為T[N].

3 模型結(jié)構(gòu)選擇

第2節(jié)給出了模型結(jié)構(gòu)已知情況下的時變模態(tài)參數(shù)估計方法，在實際情況下，模型結(jié)構(gòu)往往是未知的，需要首先解決模型結(jié)構(gòu)的選擇問題。對于本文所提的方法，模型結(jié)構(gòu)可以分為兩個子問題：一是模型階數(shù)選擇；二是模型的結(jié)構(gòu)參數(shù)選擇。

遞推小波基FS-TARMA模型的階數(shù)選擇包括AR階數(shù)na、MA階數(shù)nc、時變AR擬合階數(shù)pa和時變MA擬合系數(shù)pc. 一般地，為了控制模型的過擬合程度，需要采用一些準(zhǔn)則對擬合程度進(jìn)行定量描述，本文采用殘差平方和(RSS)準(zhǔn)則對模型的階數(shù)進(jìn)行選取。定義FS-TARMA模型的RSS為

(25)

步驟1在搜索區(qū)間內(nèi)對na=nc進(jìn)行搜索，當(dāng)RSS曲線取極小值或趨于平緩時，即可確定AR階數(shù)na.

步驟2在確定AR階數(shù)na后，以[0,na]為MA階數(shù)搜索區(qū)間，依次對nc進(jìn)行賦值，當(dāng)RSS曲線取極小值或趨于平緩時，即可確定MA階數(shù)nc.

步驟3在確定AR階數(shù)na和MA階數(shù)nc后，可假定pa=pc，給定搜索區(qū)間，最佳的估計值出現(xiàn)在RSS曲線取極小值或趨于平緩處。

在遞推過程中涉及到模型結(jié)構(gòu)參數(shù)α的選擇，協(xié)方差初始矩陣P[0]=αI，其中α一般選取較大的正數(shù)(例如α=104)，旨在保證初始階段的協(xié)方差矩陣足夠大。

4 數(shù)值算例

針對彈箭的僅輸出模態(tài)分析，本質(zhì)上是對彈箭動力學(xué)響應(yīng)的分析，彈箭的動力學(xué)響應(yīng)數(shù)據(jù)一方面來源于實驗、工作條件下的實測數(shù)據(jù)；另一方面來源于彈箭物理模型的數(shù)值計算，即動力學(xué)的正問題。開展彈箭動力學(xué)正問題的研究，不僅可以提供彈箭仿真條件下的動力學(xué)響應(yīng)數(shù)據(jù)，而且基于力學(xué)原理所建立的反映系統(tǒng)物理特性的動力學(xué)模型，也對彈箭結(jié)構(gòu)設(shè)計與驗證具有重要意義。本節(jié)以阿里安Ⅴ號運載火箭(以下簡稱阿里安Ⅴ號)為例進(jìn)行時變結(jié)構(gòu)動力學(xué)建模，建模流程如圖1所示。

圖1 阿里安Ⅴ號時變結(jié)構(gòu)動力學(xué)建模流程圖Fig.1 Flow chart of time-varying dynamic modelling for Ariane Ⅴ

4.1 阿里安Ⅴ號芯級時變有限元模型

根據(jù)振動理論，一個具有nd個自由度的線性振動系統(tǒng)，可以用一個nd維2階微分方程組來描述如下：

(26)

(27)

式中：

(28)

式中：Me、Ee、Ke、fe分別為單元的質(zhì)量矩陣、阻尼矩陣、剛度矩陣和載荷向量；ρ為質(zhì)量密度；μ為阻尼系數(shù)；N為形函數(shù)；D為彈性矩陣；B為應(yīng)變矩陣；Ve為分析體內(nèi)所有單元；dV為微元體積。

目前，隨著計算機(jī)技術(shù)的發(fā)展，大型復(fù)雜有限元模型的動力學(xué)離散化一般在商業(yè)軟件內(nèi)完成。對于彈箭的時變過程，若每一個瞬時的有限元模型代表當(dāng)前時刻系統(tǒng)的物理特征，就需要根據(jù)時間步長建立大量有限元模型，而商業(yè)軟件目前并不具備分析時變結(jié)構(gòu)動力學(xué)問題的功能。因此，針對彈箭的結(jié)構(gòu)動力學(xué)正問題，以阿里安V-versatile芯級為原型，經(jīng)過調(diào)研與合理假設(shè)，采用集中質(zhì)量- 梁模型建立反映運載火箭橫向振動特性的有限元模型[26-27]，按照如下步驟進(jìn)行建模：1)將芯級各部段進(jìn)行分組編號，并建立幾何曲線。2)在得到的每段幾何曲線上進(jìn)行站點劃分，站點劃分?jǐn)?shù)目應(yīng)當(dāng)滿足正確描述模型的需求：對于變截面梁，至少需要3個站點以上進(jìn)行描述；對于等截面梁，則根據(jù)長度進(jìn)行分配，本文模型中對等截面梁單元以0.5 m的尺度進(jìn)行等比例劃分，共得到96個節(jié)點、95個梁單元。3)梁截面定義為薄壁圓環(huán)，芯級端頭和尾段為變截面梁，將結(jié)構(gòu)質(zhì)量和液體質(zhì)量分別分配到對應(yīng)的站點上，結(jié)構(gòu)質(zhì)量為集中質(zhì)量，液體質(zhì)量則在相應(yīng)的站點上添加對應(yīng)耦合質(zhì)量。表1給出了阿里安Ⅴ號芯級的建模信息。圖形界面的編寫不在本文研究范圍內(nèi)，并聚焦于結(jié)構(gòu)特性上的變化對模態(tài)參數(shù)的影響，最終建立的阿里安Ⅴ號芯級模型示意圖可參考其在有限元軟件中的顯示結(jié)果，如圖2所示。

表1 阿里安Ⅴ號芯級建模關(guān)鍵參數(shù)Tab.1 Modal structure parameters of Ariane Ⅴ

圖2 阿里安Ⅴ號芯級有限元模型Fig.2 Finite element model of Ariane Ⅴ

考慮一種質(zhì)量消耗工況，如圖3所示。

圖3 芯級節(jié)點質(zhì)量特性突變規(guī)律Fig.3 Time-varying mass feature of the first stage node

采用有限帶寬白噪聲激勵作為激勵源，低通截止頻率為80 Hz，作用位置為頭錐頂點，即整箭幾何原點，作用方向沿y軸正向。通過對時變系統(tǒng)進(jìn)行“時間凍結(jié)”分析，以1 000 Hz的采樣頻率獲得整箭的凍結(jié)結(jié)構(gòu)集合。在凍結(jié)結(jié)構(gòu)集合中，給定底部固支條件，提取阿里安Ⅴ號的整箭凍結(jié)結(jié)構(gòu)無阻尼橫向振動頻率，如圖4所示。

圖4 阿里安Ⅴ號“凍結(jié)結(jié)構(gòu)”橫向模態(tài)頻率Fig.4 Frozen structure modal frequency of Ariane Ⅴ

使用Newmark-beta數(shù)值積分法求解阿里安Ⅴ號時變結(jié)構(gòu)的動響應(yīng)，并設(shè)置比例阻尼E(t)=a0M(t)+a1K(t)，其中a0=0.1,a1=0.001. 提取1號、40號和80號節(jié)點的動響應(yīng)信號如圖5所示。

圖5 加速度響應(yīng)Fig.5 Acceleration responses

圖7 小波基FS-TARMA模型的階數(shù)選擇Fig.7 Order selection of wavelet basis FS-TARMA model

取節(jié)點1的加速度響應(yīng)信號作時頻分析，以80 Hz的頻率進(jìn)行重采樣，采用短時傅里葉變換對重采樣后的加速度響應(yīng)信號進(jìn)行分析，得到其自功率譜如圖6所示。由圖6可知，1階模態(tài)的能量與其他階相比很小，在功率譜中只能看到微弱的峰值。較弱的模態(tài)難以用辨識方法辨識出，因此，本文后續(xù)只針對運載火箭的第2、3、4階模態(tài)的辨識效果進(jìn)行討論。

4.2 模型參數(shù)選擇與計算結(jié)果

根據(jù)第3節(jié)給出的模型結(jié)構(gòu)選擇方法，為簡化選擇流程，首先給定一組較為完備的泛函空間，取pa=pc=40. 根據(jù)RSS準(zhǔn)則選取模型的AR和MA階數(shù)如圖7所示，即na=10,nc=8. 在初始化參數(shù)α=104時，利用節(jié)點1到節(jié)點96共96組振動響應(yīng)數(shù)據(jù)進(jìn)行分析得到模態(tài)參數(shù)的辨識結(jié)果，如圖8所示。

圖8 遞推小波基FS-TARMA模型凍結(jié)結(jié)構(gòu)在不同 階數(shù)模態(tài)頻率的辨識結(jié)果Fig.8 Identified results of “frozen structure” modal frequencies with different orders using recursive wavelet basis FS-TARMA moclel

為了比較遞推算法與批量算法的辨識效果，給出在相同模型結(jié)構(gòu)下的小波基FS-TARMA辨識結(jié)果，如圖9所示。

圖9 小波基FS-TARMA模型凍結(jié)結(jié)構(gòu)在不同階數(shù) 模態(tài)頻率的辨識結(jié)果Fig.9 Identified results of “frozen structure” modal frequencies with different orders using wavelet basis FS-TARMA model

由圖8可知，本文所提的遞推小波基FS-TARMA方法能夠有效地識別時變系統(tǒng)的動態(tài)特征，在3階模態(tài)頻率上均能較好地跟蹤系統(tǒng)的時變特征。通過圖8和圖9的對比可知：在辨識效果上，本文所提遞推算法的虛假極點比非遞推算法少，并且跟蹤精度更高；在辨識精度上，二者在圖中的差距并不明顯，都體現(xiàn)出了良好的辨識精度。

為了定量比較不同方法的辨識精度與時變跟蹤能力，定義均值絕對誤差(MAE)[28]為

(29)

式中：f(t)為“凍結(jié)結(jié)構(gòu)”模態(tài)頻率值；l(t)為模態(tài)頻率辨識值；L為數(shù)據(jù)長度。引入MAE計算方法，則兩種方法的模態(tài)頻率辨識結(jié)果的定量比較如表2所示。

表2 固有頻率估計誤差Tab.2 MAE values for identification methods Hz

由表2可知，本文所提遞推小波基FS-TARMA方法的辨識精度略低于批量算法，但是辨識精度上相當(dāng)接近，3階模態(tài)辨識頻率最大相對誤差在5%以內(nèi)。為了進(jìn)一步比較不同計算方法之間的計算效率，將估計過程所花費的CPU計算時間用于表征不同方法的計算效率，使用的計算機(jī)CPU主頻為3.5 GHz，內(nèi)存16 GB. 兩種方法的計算時間如圖10所示。

圖10 CPU計算時間Fig.10 CPU computation time

由圖10可知，在相同的時變系統(tǒng)、模型結(jié)構(gòu)和計算條件下，本文所提遞推算法計算效率與原方法相比，得到了較大的提升(提升約為9.38倍)。可見，本文所提算法具有更好的在線辨識能力，通過對結(jié)構(gòu)化TARMA模型的遞推化，不僅保留了結(jié)構(gòu)化TARMA模型辨識精度高的優(yōu)點，而且大大提升了其計算效率，為彈箭飛行模態(tài)的辨識提供了新的思路。

5 結(jié)論

針對彈箭在飛行狀態(tài)下的時變模態(tài)參數(shù)辨識問題，本文基于FS-TARMA模型，利用小波基函數(shù)作為泛函基底構(gòu)建時變參數(shù)矩陣，并借鑒于無結(jié)構(gòu)化TARMA模型的遞推思想，將投影參數(shù)矩陣視為數(shù)據(jù)長度的變量，提出了一種新的遞推估計算法，實現(xiàn)了投影參數(shù)矩陣的遞推估計。通過建立的阿里安Ⅴ號芯級運載火箭時變有限元模型，對所提辨識方法進(jìn)行了驗證，辨識結(jié)果表明：本文所提算法保留了批量算法辨識精度高的優(yōu)點，同時大大提升了計算效率。