一種運載火箭時變結構模態參數辨識的確定性演化方法

余 磊，劉 莉,2，崔 穎，岳振江，康 杰

(1. 北京理工大學宇航學院，北京 100081; 2. 北京理工大學飛行器動力學與控制教育部重點實驗室，北京 100081)

0 引 言

運載火箭的僅輸出時變結構模態參數辨識是航空航天工程領域的關鍵技術。運載火箭在工作狀態下由于構型變化、模塊化結構調整、載荷工況變化等過程，使得自身成為變結構、變質量、變參數的典型時變結構，獲取結構的模態參數等動力學特性對于實現高精度控制起著關鍵作用，并直接影響復雜空間任務的順利完成。

運載火箭的時變特性體現在微分方程系數中[1]，并進一步反映在振動響應信號中，即時變結構的振動響應信號具有非平穩特性[2]。針對非平穩信號的辨識方法主要分為時域和頻域，其中時域辨識方法直接以時域信號為基礎進行模態參數辨識，物理概念更為清晰。時域辨識方法又分為基于狀態空間的辨識方法[3-5]和基于時間序列的辨識方法，其中時變自回歸滑動平均模型(Time-varying autoregressive moving average, TARMA)是時域僅輸出辨識模型，作為一種參數化辨識模型，TARMA模型具有許多優點[6]：1)表達的簡約性；2)更高的精度；3)更好的時變動態特性的跟蹤能力；4)更高的分辨率；5)更好的靈活性。Spiridonakos等[7]根據TARMA模型的時變系數的演化形式不同，分為三類：無結構化、隨機結構化和確定性結構化。

確定性結構化的時變系數演化形式即是將TARMA模型的時變系數用一組基函數進行展開，進而將時變系數轉換為投影空間中的時不變投影系數，選擇合適的基函數可以提高辨識精度。這類方法統稱為泛函序列(Functional series，FS)TARMA模型。Rao[8]對TARMA模型的時變系數進行展開，基于時間多項式首次提出了FS-TAR模型。Mendel等[9]對一般的泛函時間序列模型進行了研究，討論了時變系數的快變與慢變，建立了相關的理論基礎。Liporace[10]基于任意階多項式對FS-TAR模型的時變系數進行展開，并將模型應用于語音分析。Grenier[11]將FS-TAR模型首次拓展到FS-TARMA模型，并從數學上證明了FS-TARMA模型的存在性和唯一性。Walter等[12]研究了模型參數估計的精度與基函數類型的關系，并指出足夠階數的任意函數簇能以任意精度逼近時變系數。Zou等[13]采用了Walsh函數，對分段時變系數進行了展開。目前，多種多樣的基函數被用于描述FS-TARMA模型的時變系數，如勒讓德多項式、切比雪夫多項式[14]、雅克比多項式[15]、傅里葉級數[16-17]、B樣條、無網格形函數[18]等。不同基函數所適應的問題各不相同，如勒讓德多項式、傅里葉級數適應于慢變/長周期信號，B樣條適應快變/慢變信號，無網格形函數適應快變信號。

小波方法因具有“數學放大鏡”的作用，近年來也被頻繁應用于頻域模態參數辨識方法中。文獻[19-20]將小波變換應用于脈沖響應函數提取和頻率響應的估計，給出了響應的模態參數辨識方法。Staszewski等[21]基于小波變換將傳統頻率響應函數拓展到時頻域，利用時間平均法提高了信噪比。Dziedziech等[22]進一步將該方法應用于頻率突變系統，結果表明，小波函數的多分辨率特性能夠在頻率突變處很好的擬合，得到了更為光滑的時頻分布圖。小波方法對快變信號的良好適應性，在時頻域/頻域模態參數辨識方法上得到了廣泛應用[23-24]，卻鮮有在時域僅輸出模態參數辨識方法上的應用。

本文針對運載火箭運行工況下的質量變化特性，在確定性演化方法上，利用小波函數能夠較好地反映信號局部細節的特點，選擇墨西哥帽小波函數作為TARMA模型時變系數泛函空間基底，構建了基于小波基的FS-TARMA模型。利用分布解耦思想以及TARMA模型與無限階截斷模型相等價的原理，發展了兩步最小二乘法對估計過程進行解耦，得到了時變系數的解耦估計。通過運載火箭時變有限元模型驗證了提出的墨西哥帽小波基FS-TARMA方法具有良好的辨識精度與時變跟蹤能力。

1 TARMA模型

TARMA模型是基于非平穩時間序列建立起來的離散模型，具有隨機差分方程的形式，不僅可揭示動態數據本身的結構和規律，還可定量地了解觀測數據之間的線性相關性，是解決動力學第一類反問題的有效工具。令na和nc分別表示AR階數和MA階數，TARMA(na,nc)模型可寫為如下形式。

(1)

式中:x[t]∈Rk×1表示離散的非平穩振動信號，e[t]∈Rk×1滿足均值為零且協方差為Σ[t]∈Rk×k的正態分布，Ai[t]和Ci[t]分別表示TARMA模型的時變AR和MA系數矩陣。

TARMA模型總是將輸入信號假設為白噪聲，因此保證了模型所表示的等價系統的固有特性包含了真實物理系統的固有特性，而沒有丟失信息。在實際應用中當輸入信號不可知時，白噪聲激勵假設是一種簡單易行的方式，解決了系統辨識問題中輸入信號未知時的困難，因此TARMA模型被廣泛地應用于僅輸出模態參數辨識領域。

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

(2)

(3)

TARMA模型“時間凍結”固有頻率和阻尼比為

(4)

2 TARMA模型的時變系數確定性演化方法

2.1 基于墨西哥帽小波基的TARMA模型

對于小波展開，可構造一個兩參數系統，使得

(5)

其中,j和d為整數指標，ρj,d(t)是小波基函數，通常形成一組正交基。展開的系數集aj,d稱為f(t)的離散小波變換。采用墨西哥帽小波函數作為基本小波，即母小波，對其進行二進尺度變換和時間軸上的平移，為方便運算，記小波基最大尺度因子等于時間基函數最大維數，可得其小波基為

(6)

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

(7)

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

e[t]?x[t]=

e[t]?x[t]=θTφ[t]+e[t]

(8)

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

(9)

(10)

其中，“?”為Kronecker積算子，基函數向量ρ[t]具有如下形式

(11)

由式(9)可知，參數矩陣θT為常數矩陣。通過將TARMA模型參數矩陣表示為一組以墨西哥帽小波基函數的線性組合，可將時變的TARMA模型轉化為時不變的FS-TARMA模型，從而把時變模型參數估計問題轉化為時不變模型參數估計問題。

2.2 時變系數的兩步最小二乘估計

根據參數估計原理，對式(8)中TARMA模型的殘差序列均方誤差(MSE)求和，得到估計的費用函數為

(12)

利用最小化MSE準則(MMSE)，獲得待估參數估計

(13)

式中:arg min(·)為最小化算子，即當費用函數取最小值時變量的取值。

從式(12)不難看出，由于殘差序列中存在MA部分，即回歸向量φ[t]中包含了殘差的過去時刻值，所以估計參數θT的過程是一個非線性優化問題。為解決該問題，常常需要將非線性優化過程轉化為線性求解過程，鑒于TARMA模型具有截斷形式的特點，本節采用兩步最小二乘法對該非線性優化問題進行線性化求解。

TARMA模型的一個基本性質是：任意階的TARMA模型都可以用一個無限階的TAR模型或TMA模型(也稱為逆函數表達)來等價表示[7]。兩步最小二乘法的基本思想是利用TARMA模型的等價性質，將TARMA模型轉化為兩個可以線性求解的TAR模型：第一步將TARMA模型相等價的無限階TAR模型進行截斷，求出該模型的殘差序列值；第二步利用該殘差序列值當做先驗已知的MA部分，代入到原TARMA模型中。TARMA模型的逆函數可表達為

C-1[B,t]°A[B,t]·x[t]=

C-1[B,t]°C[B,t]·e[t]?

H-1[B,t]x[t]=e[t]

(14)

兩步最小二乘的步驟可分為如下兩個部分：

1) 逆函數估計

對式(14)中H-1[t]取ng階截斷，截斷后的逆函數模型如式(15)所示。

x[t]=hT·φH[t]+e[t,h]?

(15)

(16)

(17)

2)AR和MA系數估計

(18)

同理利用最小化MSE準則，獲得待估參數估計

(19)

2.3 模型結構問題的討論

第2.2節討論了模型結構已知情況下的墨西哥帽小波基FS-TARMA時變模態參數估計方法，在實際情況下，模型結構往往是未知的，需要首先解決模型結構的選擇問題。對于確定性結構化TARMA模型而言，模型結構的選擇問題可以分為兩個子問題：一是模型的階數選擇，二是模型的結構參數選擇。

對于墨西哥帽小波基FS-TARMA時變模態參數估計方法的模型結構選擇問題，主要包括對逆函數截斷階數ng、AR階數na、MA階數nc、時變AR系數擬合階數pa和時變MA擬合系數pc的選取。

為了控制模型的過擬合程度，需要采取一些準則對擬合程度進行定量描述，常用的準則包括赤池信息準則(Akaike Information Criterion，AIC)、貝葉斯信息準則(Bayesian Information Criterion)和殘差平方和準則(Residual Sum of Squares，RSS)[27-28]，本文采用RSS準則對TARMA模型的模型階數進行選取。定義TARMA模型的殘差平方和

(20)

殘差平方和在一定程度上描述了TARMA模型對信號的擬合誤差，因此可以控制模型的過擬合程度。基于RSS準則選取TARMA模型的模型階數的步驟如下。

值得說明的是，較高的模型階數可以更好地擬合信號，從而RSS準則一般相對較小；較低的模型階數意味著較高的計算效率，更低的計算復雜度。從信號擬合的角度來看，模型階數偏向于取大，然而在TARMA模型的模態參數辨識過程中，需要考慮實際的物理系統。當模型階數大于實際物理階數的時候，系統的特征根數量增多，辨識結果中會出現大量的虛假模態，虛假模態的剔除是十分困難的；當模型階數過低，極端情況下往往會無法辨識出正確模態。因此，在模態參數辨識領域，模型階數的選取不僅是信號擬合與計算效率之間的折中，更是物理系統先驗知識和計算效率之間的折中。

3 算例驗證

針對運載火箭的僅輸出模態分析，本質上是對運載火箭動力學響應的分析，運載火箭的動力學響應數據一方面來源于實驗/工作條件下的實測數據；另一方面來源于運載火箭物理模型的數值計算，即動力學的正問題。運載火箭的動力學響應實測數據對于研究者而言往往代價高昂，因此研究運載火箭的正問題不僅能獲得系統的動力學響應，基于力學原理建立的反映系統物理特性的動力學模型對運載火箭結構設計與驗證也具有重要意義。

為了驗證本文所提出的墨西哥帽小波基FS-TARMA時變模態參數估計方法的有效性，本節將對運載火箭時變動力學正問題進行建模，為方法的驗證提供數據基礎，時變結構動力學建模的流程如圖1所示。

圖1 運載火箭時變結構動力學建模流程圖Fig.1 The work flow of time-varying dynamic modelling for launch vehicle

3.1 運載火箭時變有限元模型

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

(21)

目前，商業軟件目前并不具備分析時變結構動力學問題的功能。因此，針對運載火箭的結構動力學正問題，以阿里安V-versatile芯級為原型，經過調研與合理假設，將芯級的各部段進行分組編號，并建立幾何曲線。在得到的每段幾何曲線上進行站點的劃分，站點劃分的數目應當滿足正確描述模型的需求：對于變截面梁，至少需要三個站點以上進行描述；對于等截面梁，則根據長度進行分配，本模型中對等截面梁單元以0.5 m的尺度進行等比例劃分，共得到了96個節點，95個梁單元。梁截面定義為為薄壁圓環，芯級端頭和尾段為變截面梁，將結構質量和液體質量分別分配到對應的站點上，結構質量為集中質量，液體質量則在相應的站點上添加對應耦合質量。表1給出了阿里安V號芯級的建模信息。

表1 阿里安V號芯級建模關鍵參數Table 1 The model structure parameters of Ariane V

考慮一種質量局部顯著變化工況，如圖2所示。

圖2 芯一級節點質量特性變化規律Fig.2 The time-varying mass feature of the first stage nodes

采用有限帶寬白噪聲激勵作為激勵源，低通截止頻率為80 Hz，作用位置為頭錐頂點，即整箭幾何原點，作用方向沿y軸正向。

這里引入“時間凍結”概念，即在忽略結構參數的變化率后，慢時變結構連續時間變化的模態參數近似等于“特征凍結”的結構模態參數。通過對時變系統進行“時間凍結”分析，以1000 Hz的采樣率獲得整箭的凍結結構集合。在凍結結構集合中，給定底部固支條件，提取阿里安V號的整箭凍結結構無阻尼橫向振動頻率，如圖3所示。

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

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

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

圖5 節點1加速度響應自功率譜Fig.5 Auto-power spectrum of node 1

取節點1的加速度響應信號作時頻分析，以80 Hz 進行重采樣，采用短時傅里葉變換對重采樣后的加速度響應信號進行分析，得到其自功率譜如圖5所示。由圖5可知，一階模態的能量與其他階相比很小，在功率譜中只能看到微弱的峰值。較弱的模態難以用辨識方法辨識出，因此，本文后續只針對運載火箭的第二、三、四階模態的辨識效果進行討論。為了更加全面地比較墨西哥帽小波基FS-TARMA方法辨識結構局部特征的有效性，本節將引入傳統傅里葉基FS-TARMA方法作為對比。

3.2 辨識結果

1)墨西哥帽小波基FS-TARMA辨識結果

對于墨西哥帽小波基FS-TARMA方法，根據第2.3節中提到的模型結構選擇方法，首先給定一組較為完備的泛函空間，取pa=pc=35，逆函數截斷階數ng=20，使得模型結構選擇問題簡化，只需根據RSS準則選取AR階數na與MA階數nc，由圖5可知，在感興趣的頻帶范圍內，具有4階橫向頻率，因此na至少需要為6，設置初始搜索區間[6,12]，以節點1的振動響應信號選取模型階數如圖6所示，na=10和nc=8時RSS曲線變得平緩，因此選擇其作為模型階數的推薦值。最終確定的模型結構為pa=pc=35、ng=20、na=10和nc=8，利用節點1～96共96組振動響應數據進行分析，基于墨西哥帽小波基FS-TARMA方法的辨識結果如圖7所示。

圖6 墨西哥帽小波基FS-TARMA模型階數選擇Fig.6 Model selection of Mexican hat wavelet basis FS-TARMA

圖7 墨西哥帽小波基FS-TARMA模態參數辨識結果Fig.7 Modal frequency estimation results of Mexican hat wavelet basis FS-TARMA

2)傅里葉基FS-TARMA辨識結果

對于傅里葉基FS-TARMA方法，其模型結構與墨西哥帽小波基FS-TARMA方法的模型結構類似，區別在于其擬合時變系數的泛函空間的基函數不同，本方法采用傅里葉基作為泛函空間的基底，其形式如下所示[7]:

(22)

傅里葉基FS-TARMA方法的模型結構包括截斷階數ng、AR階數na、MA階數nc、時變AR系數擬合階數pa和時變MA擬合系數pc。同理，根據第2.3節中提到的模型結構選擇方法，首先給定一組較為完備的泛函空間，取pa=pc=35和逆函數截斷階數ng=20，使得模型結構選擇問題進一步簡化，只需要根據RSS準則選取關鍵的AR階數na與MA階數nc，以節點1的振動響應信號選取模型階數如圖8所示。最終確定的模型結構為pa=pc=35、ng=20、na=10和nc=8，利用節點1～96共96組振動響應數據進行分析，基于傅里葉基FS-TARMA方法的辨識結果如圖9所示。

圖8 傅里葉基FS-TARMA模型階數選擇Fig.8 Model selection of Fourier basis FS-TARMA

圖9 傅里葉基FS-TARMA模態參數辨識結果Fig.9 Modal frequency estimation results of Fourier basis FS-TARMA

3)方法的對比分析

由圖7和圖9可知，兩種方法均能較好地辨識出時變結構的第二、三、四階橫向固有頻率。此外，與傅里葉基FS-TARMA的模態參數辨識結果相比，墨西哥帽小波基FS-TARMA的虛假極點較少，辨識結果較為“干凈”。

為了定量比較不同模型的辨識精度與時變跟蹤能力，定義均值絕對誤差(Mean absolute error，MAE)為[29]

(23)

表2 模態頻率估計誤差Table 2 MAE values for different methods

由表2可知，本文所提的小波基FS-TARMA方法，在辨識精度上要略好于基于傳統的傅里葉基FS-TARMA方法。為了能夠在局部模態頻率迅速變化處更好地比較兩種方法的優劣，這里對辨識結果進行局部放大，如圖10所示。

圖10 模態參數辨識結果局部放大圖Fig.10 The local detail of modal frequency estimation results for different methods

由圖10可知，對于模態頻率拐點處，墨西哥帽小波基FS-TARMA方法可以更好地擬合局部細節特征，而傳統的傅里葉基FS-TARMA方法則只能平滑過渡。因此，對于具有局部細節特征的非平穩振動信號的辨識問題，與傳統傅里葉基FS-TARMA相比，本文所提的小波基FS-TARMA時變結構模態參數估計方法虛假極點更少，估計精度更高，并且能夠準確反映模態局部特征。

4 結 論

本文針對運載火箭時變結構模態參數辨識問題，基于TARMA模型，在時變系數的確定性演化方法上展開了研究。利用小波函數能夠較好反映信號局部細節的特點，選擇墨西哥帽小波函數作為TARMA模型時變系數泛函空間的基底，構建了基于墨西哥帽小波基的FS-TARMA模型。利用分步解耦思想以及TARMA模型與無限階截斷模型相等價的原理，采用兩步最小二乘法對時變系數實現了解耦估計。通過建立的阿里安V號芯級運載火箭時變有限元模型，對所提辨識方法進行了驗證，結果表明：針對非平穩振動響應信號，本文所提墨西哥帽小波基FS-TARMA方法可以準確地對時變系統的模態頻率進行估計，與傳統傅里葉基FS-TARMA相比，具有更好的辨識精度，并能夠較好地反映出模態局部細節特征。