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

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

(西安現代控制技術研究所, 陜西 西安 710065)

0 引言

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

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

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

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

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

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

1 小波基FS-TARMA模型

1.1 TARMA模型描述

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

(1)

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

TARMA模型“時間凍結”模態參數可通過在每一時刻求解如下特征值問題[22]獲得，

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

(2)

(3)

式中：I為單位矩陣。

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

(4)

1.2 基函數選擇

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

(5)

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

(6)

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

(7)

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

(8)

(9)

式中：?為Kronecker積算子；基函數向量ρ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)式可知，參數矩陣θT為常數矩陣。通過將TARMA模型參數矩陣表示為一組以墨西哥帽小波基函數的線性組合，可將時變的TARMA模型轉化為時不變的FS-TARMA模型，從而把時變模型參數估計問題轉化為時不變模型參數估計問題。估計得到的θT代入(6)式中，即可得到系統的模態頻率。

2 時變系數的遞推估計方法

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

(11)

參數矩陣θT的估計為

(12)

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

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

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

(13)

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

通過對費用函數進行求導，可得

(14)

式中：

(15)

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

(16)

它的逆矩陣可以表示為

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

(17)

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

(18)

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

(19)

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

(20)

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

(21)

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

(22)

則(19)式可以表達為

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

(23)

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

(24)

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

3 模型結構選擇

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

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

(25)

步驟1在搜索區間內對na=nc進行搜索，當RSS曲線取極小值或趨于平緩時，即可確定AR階數na.

步驟2在確定AR階數na后，以[0,na]為MA階數搜索區間，依次對nc進行賦值，當RSS曲線取極小值或趨于平緩時，即可確定MA階數nc.

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

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

4 數值算例

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

圖1 阿里安Ⅴ號時變結構動力學建模流程圖Fig.1 Flow chart of time-varying dynamic modelling for Ariane Ⅴ

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

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

(26)

(27)

式中：

(28)

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

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

表1 阿里安Ⅴ號芯級建模關鍵參數Tab.1 Modal structure parameters of Ariane Ⅴ

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

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

圖3 芯級節點質量特性突變規律Fig.3 Time-varying mass feature of the first stage node

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

圖4 阿里安Ⅴ號“凍結結構”橫向模態頻率Fig.4 Frozen structure modal frequency of Ariane Ⅴ

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

圖5 加速度響應Fig.5 Acceleration responses

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

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

4.2 模型參數選擇與計算結果

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

圖8 遞推小波基FS-TARMA模型凍結結構在不同 階數模態頻率的辨識結果Fig.8 Identified results of “frozen structure” modal frequencies with different orders using recursive wavelet basis FS-TARMA moclel

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

圖9 小波基FS-TARMA模型凍結結構在不同階數 模態頻率的辨識結果Fig.9 Identified results of “frozen structure” modal frequencies with different orders using wavelet basis FS-TARMA model

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

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

(29)

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

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

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

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

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

5 結論

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