谷神星一號(hào)運(yùn)載火箭的時(shí)變模態(tài)參數(shù)辨識(shí)

馬志賽 丁千 劉百奇 劉建設(shè) 張軍鋒 董威利

摘要： 針對(duì)運(yùn)載火箭這一大型復(fù)雜系統(tǒng)，單純依靠機(jī)理分析難以建立能夠精確描述其時(shí)變特征的動(dòng)力學(xué)模型，借助模態(tài)辨識(shí)技術(shù)獲取其在飛行狀態(tài)下的時(shí)變動(dòng)力學(xué)特性尤為必要。面向運(yùn)載火箭的飛行模態(tài)辨識(shí)需求，發(fā)展基于時(shí)變自回歸滑動(dòng)平均模型的僅輸出遞推辨識(shí)方法，引入指數(shù)加權(quán)遺忘機(jī)制進(jìn)行時(shí)變特性跟蹤，能夠在激勵(lì)未知的情況下僅基于響應(yīng)數(shù)據(jù)準(zhǔn)確快速辨識(shí)運(yùn)載火箭的時(shí)變模態(tài)。以谷神星一號(hào)運(yùn)載火箭為例，開展飛行遙測(cè)數(shù)據(jù)的快速處理，準(zhǔn)確辨識(shí)發(fā)射準(zhǔn)備階段和飛行階段的關(guān)鍵低階時(shí)變模態(tài)，辨識(shí)結(jié)果與有限元分析結(jié)果的變化規(guī)律相互吻合，驗(yàn)證了所提時(shí)變模態(tài)僅輸出遞推辨識(shí)方法是有效可行的，服務(wù)了運(yùn)載火箭有限元模型修正和姿控系統(tǒng)設(shè)計(jì)的實(shí)際工程需求。

關(guān)鍵詞： 運(yùn)載火箭； 時(shí)變模態(tài)參數(shù)； 谷神星一號(hào)； 遞推辨識(shí)

中圖分類號(hào)： V475.1； O327??? 文獻(xiàn)標(biāo)志碼： A??? 文章編號(hào)： 1004-4523（2024）05-0729-08

DOI： 10.16385/j.cnki.issn.1004-4523.2024.05.001

引? 言

太空探索是世界強(qiáng)國(guó)為提升綜合國(guó)力、搶占科技制高點(diǎn)、保持核心競(jìng)爭(zhēng)力而大力發(fā)展的戰(zhàn)略性新興產(chǎn)業(yè)。航天器進(jìn)入太空的過程需要借助運(yùn)載器來克服地球引力，新一代運(yùn)載火箭追求更高的運(yùn)載效率，其尺寸與推力不斷增大，面臨的發(fā)射與運(yùn)行環(huán)境愈加惡劣，結(jié)構(gòu)系統(tǒng)設(shè)計(jì)難度日益增加。例如，針對(duì)運(yùn)載火箭這一大型復(fù)雜系統(tǒng)，單純依靠機(jī)理分析難以建立能夠精確描述其時(shí)變特征的動(dòng)力學(xué)模型，仍需要依靠模態(tài)辨識(shí)技術(shù)來獲取其真實(shí)飛行狀態(tài)下的時(shí)變動(dòng)力學(xué)特性。因此，為更好地滿足航天工程高可靠、低成本和快速響應(yīng)的發(fā)射需求，從結(jié)構(gòu)動(dòng)力學(xué)反問題入手，開展運(yùn)載火箭時(shí)變模態(tài)參數(shù)辨識(shí)的關(guān)鍵技術(shù)攻關(guān)和方案深化論證，對(duì)解決運(yùn)載火箭的相關(guān)動(dòng)力學(xué)與控制問題意義重大［1?3］。

由于燃料快速消耗和級(jí)間分離，飛行狀態(tài)下的運(yùn)載火箭會(huì)表現(xiàn)出顯著的時(shí)變特征，準(zhǔn)確獲取其時(shí)變模態(tài)對(duì)姿控系統(tǒng)設(shè)計(jì)、動(dòng)載荷設(shè)計(jì)、天地差異研究和系統(tǒng)優(yōu)化設(shè)計(jì)等都具有重要意義［2?6］。以姿控系統(tǒng)設(shè)計(jì)為例，運(yùn)載火箭主要依靠陀螺敏感姿態(tài)進(jìn)行控制以實(shí)現(xiàn)穩(wěn)定飛行。在實(shí)際飛行中，陀螺感受到的姿態(tài)信息不僅包含剛體姿態(tài)，還包含箭體彈性變形所引起的附加姿態(tài)，姿控系統(tǒng)的設(shè)計(jì)必須考慮彈性振動(dòng)的影響，否則可能會(huì)導(dǎo)致姿態(tài)發(fā)散，飛行失?。??3，7］。然而為了保證良好的氣動(dòng)特性與運(yùn)載效率，運(yùn)載火箭常常呈現(xiàn)出長(zhǎng)細(xì)比大、結(jié)構(gòu)質(zhì)量占比低的特點(diǎn)，導(dǎo)致其彈性振動(dòng)頻率接近剛體穿越頻率且具有較寬的變化范圍，在擴(kuò)大陷波濾波器凹陷寬度的同時(shí)，還需要增加相位穩(wěn)定方案。例如，Ares?I運(yùn)載火箭的長(zhǎng)細(xì)比約為17，其一級(jí)火箭結(jié)構(gòu)的一階彎曲模態(tài)低至0.9 Hz，與一級(jí)火箭PID控制器的帶寬非常接近［8］，燃料快速消耗會(huì)導(dǎo)致固有頻率不斷發(fā)生變化，大大增加了姿控系統(tǒng)的設(shè)計(jì)難度［9］。此外，以往飛行試驗(yàn)的遙測(cè)數(shù)據(jù)表明，由于無法準(zhǔn)確模擬飛行狀態(tài)下運(yùn)載火箭的真實(shí)工作環(huán)境，通過有限元分析或地面振動(dòng)試驗(yàn)獲得的模態(tài)參數(shù)往往與飛行狀態(tài)下的真實(shí)值存在一定差異，這進(jìn)一步增加了姿控系統(tǒng)的設(shè)計(jì)難度［10?12］。因此，充分利用天地試驗(yàn)數(shù)據(jù)和模態(tài)辨識(shí)理論發(fā)展運(yùn)載火箭時(shí)變模態(tài)辨識(shí)技術(shù)，將是未來實(shí)現(xiàn)運(yùn)載火箭彈性振動(dòng)智能控制與實(shí)時(shí)決策的必經(jīng)途徑。

近十年來，隨著辨識(shí)理論與測(cè)試技術(shù)的不斷發(fā)展，工作模態(tài)分析在航天領(lǐng)域已取得一定進(jìn)展。James等［13?14］對(duì)運(yùn)載火箭工作模態(tài)分析方法進(jìn)行了綜述，并指出時(shí)變的系統(tǒng)特性、非定常的運(yùn)行環(huán)境及其所導(dǎo)致的非平穩(wěn)振動(dòng)響應(yīng)，是解決運(yùn)載火箭時(shí)變模態(tài)辨識(shí)問題的主要挑戰(zhàn)。De Vivo等［5］基于歐空局Vega火箭的飛行數(shù)據(jù)，采用自然激勵(lì)技術(shù)（NExT）生成相關(guān)函數(shù)，在短時(shí)時(shí)不變假設(shè)下采用最小二乘復(fù)指數(shù)法估計(jì)模態(tài)參數(shù)，與有限元分析結(jié)果吻合良好。張家雄等［15］采用PolyMAX方法對(duì)飛行器低空飛行試驗(yàn)中典型部位的低頻振動(dòng)實(shí)測(cè)數(shù)據(jù)進(jìn)行辨識(shí)，將模態(tài)辨識(shí)結(jié)果與地面振動(dòng)試驗(yàn)進(jìn)行對(duì)比，發(fā)現(xiàn)兩者固有頻率接近，模態(tài)振型基本吻合，但也存在一定偏差。針對(duì)飛行遙測(cè)數(shù)據(jù)信噪比低、辨識(shí)結(jié)果虛假模態(tài)多的問題，王亮等［16?17］將自回歸滑動(dòng)平均（ARMA）方法與NExT相結(jié)合用于飛行模態(tài)辨識(shí)，并基于穩(wěn)定圖對(duì)真實(shí)模態(tài)進(jìn)行了篩選。隨后王亮等［6］又研究了特征時(shí)段選擇對(duì)飛行模態(tài)辨識(shí)結(jié)果的影響，指出在使用遙測(cè)數(shù)據(jù)進(jìn)行模態(tài)辨識(shí)時(shí)，應(yīng)選擇低階模態(tài)響應(yīng)信噪比較高的數(shù)據(jù)段。南宮自軍等［2］以Ares I?X運(yùn)載火箭飛行模態(tài)辨識(shí)［4］為案例，對(duì)飛行模態(tài)辨識(shí)技術(shù)的研究現(xiàn)狀和未來發(fā)展方向進(jìn)行了綜述。董嚴(yán)等［18］采用ARMA方法對(duì)某火箭飛行試驗(yàn)中三個(gè)測(cè)點(diǎn)的實(shí)測(cè)數(shù)據(jù)進(jìn)行辨識(shí)，得到了其橫向模態(tài)參數(shù)隨時(shí)間的變化規(guī)律。馬志賽等［19］對(duì)子空間方法、時(shí)間序列方法等時(shí)域辨識(shí)方法進(jìn)行了系統(tǒng)性總結(jié)，指出在復(fù)雜噪聲環(huán)境下開展在線、自適應(yīng)時(shí)變模態(tài)辨識(shí)及其工程應(yīng)用是未來的發(fā)展重點(diǎn)。Ma等［20?21］基于時(shí)變自回歸滑動(dòng)平均（TARMA）模型提出了一系列時(shí)變模態(tài)僅輸出遞推辨識(shí)方法，并搭建大長(zhǎng)細(xì)比變質(zhì)量充液筒試驗(yàn)平臺(tái)，完成了辨識(shí)方法的實(shí)驗(yàn)驗(yàn)證，為開展運(yùn)載火箭飛行模態(tài)辨識(shí)奠定了基礎(chǔ)。余磊等［22］發(fā)展了一種基于TARMA模型的時(shí)變模態(tài)批量辨識(shí)方法，并利用Ariane?5運(yùn)載火箭芯級(jí)的數(shù)值仿真數(shù)據(jù)對(duì)所提方法進(jìn)行了數(shù)值驗(yàn)證。馬慶港等［3］對(duì)運(yùn)載火箭工作模態(tài)辨識(shí)方法及流程進(jìn)行了綜述，指出準(zhǔn)確獲取模態(tài)參數(shù)對(duì)運(yùn)載火箭有限元模型修正、自適應(yīng)控制、實(shí)時(shí)狀態(tài)監(jiān)測(cè)等均具有重要意義。

隨著運(yùn)載火箭重量、推力和長(zhǎng)細(xì)比的進(jìn)一步增大，其振動(dòng)特性愈加復(fù)雜，固有頻率呈現(xiàn)出低且密集的特點(diǎn)，具有顯著時(shí)變特性的結(jié)構(gòu)系統(tǒng)與控制、氣動(dòng)、動(dòng)力等分系統(tǒng)之間的耦合問題也日漸突出。為獲取谷神星一號(hào)運(yùn)載火箭的全周期模態(tài)信息，本文旨在發(fā)展適用于運(yùn)載火箭的時(shí)變模態(tài)僅輸出遞推辨識(shí)方法，以期在激勵(lì)未知的情況下，僅基于響應(yīng)數(shù)據(jù)辨識(shí)運(yùn)載火箭的飛行模態(tài)。在此基礎(chǔ)上，借助遞推辨識(shí)方法計(jì)算量和所占用內(nèi)存空間較小，以及對(duì)計(jì)算機(jī)的處理能力和存儲(chǔ)空間要求較低的優(yōu)勢(shì)，開展谷神星一號(hào)運(yùn)載火箭飛行遙測(cè)數(shù)據(jù)的快速處理，準(zhǔn)確辨識(shí)其主要低階飛行模態(tài)，服務(wù)相關(guān)系統(tǒng)設(shè)計(jì)需求。

1 谷神星一號(hào)運(yùn)載火箭

1.1 發(fā)射記錄

谷神星一號(hào)運(yùn)載火箭是中國(guó)民營(yíng)航天企業(yè)北京星河動(dòng)力裝備科技有限公司自主研發(fā)的一款四級(jí)小型固體商業(yè)運(yùn)載火箭，一、二、三級(jí)采用固體發(fā)動(dòng)機(jī)，四級(jí)為液體上面級(jí)。該運(yùn)載火箭直徑1.4 m，全長(zhǎng)約20 m，起飛重量約33 t，500 km太陽(yáng)同步軌道運(yùn)載能力為300 kg，可滿足微小型衛(wèi)星的專屬、共享、搭載等定制化發(fā)射需求。如圖1所示，谷神星一號(hào)運(yùn)載火箭的發(fā)射記錄已有5次。

2020年11月7日，谷神星一號(hào)（遙一）運(yùn)載火箭在酒泉衛(wèi)星發(fā)射中心發(fā)射，成功將天啟星座十一星精確送入預(yù)定軌道。這是星河動(dòng)力實(shí)施的首次發(fā)射任務(wù)，也是中國(guó)民營(yíng)商業(yè)運(yùn)載火箭首次進(jìn)入500 km太陽(yáng)同步軌道，是中國(guó)商業(yè)航天的一次重大突破。

2021年12月7日，谷神星一號(hào)（遙二）運(yùn)載火箭成功發(fā)射，將搭載的5顆小衛(wèi)星送入預(yù)定軌道，首次實(shí)現(xiàn)“一箭五星”。2022年8月和11月，谷神星一號(hào)（遙三、四）運(yùn)載火箭接連發(fā)射兩次，將多顆衛(wèi)星順利送入預(yù)定軌道。2023年1月9日，谷神星一號(hào)（遙五）運(yùn)載火箭在酒泉衛(wèi)星發(fā)射中心發(fā)射，再次將搭載的5顆衛(wèi)星順利送入預(yù)定軌道。截至目前，該型火箭保持著100%的發(fā)射成功率，不斷刷新紀(jì)錄，正式進(jìn)入快響應(yīng)、高密度的規(guī)模化發(fā)射階段。

1.2 時(shí)變模態(tài)辨識(shí)問題

谷神星一號(hào)運(yùn)載火箭為三級(jí)固體發(fā)動(dòng)機(jī)加先進(jìn)液體上面級(jí)，如圖2所示。其中，芯一級(jí)配備單臺(tái)固體發(fā)動(dòng)機(jī)，海平面推力為60 t，燃燒時(shí)間為74 s；芯二級(jí)配備單臺(tái)固體發(fā)動(dòng)機(jī)，推力為28 t，燃燒時(shí)間為70 s；芯三級(jí)配備單臺(tái)固體發(fā)動(dòng)機(jī)，推力為8.8 t，燃燒時(shí)間為69 s。一、二級(jí)級(jí)間采用熱分離方式，分離時(shí)間短，二級(jí)飛行初始穩(wěn)定性好。

發(fā)動(dòng)機(jī)燃料快速消耗和級(jí)間分離會(huì)引起運(yùn)載火箭自身質(zhì)量分布的顯著變化，進(jìn)而導(dǎo)致飛行狀態(tài)下運(yùn)載火箭的模態(tài)參數(shù)（固有頻率、阻尼比、模態(tài)振型等）不斷變化。顯然，姿控等分系統(tǒng)的設(shè)計(jì)均需要根據(jù)運(yùn)載火箭實(shí)時(shí)的模態(tài)參數(shù)進(jìn)行調(diào)整，因此，能否準(zhǔn)確快速獲取飛行狀態(tài)下運(yùn)載火箭的時(shí)變模態(tài)，對(duì)各分系統(tǒng)設(shè)計(jì)至關(guān)重要。

隨著測(cè)試技術(shù)（數(shù)據(jù)采集與傳輸、傳感器布局優(yōu)化等）水平、新型傳感器性能以及計(jì)算機(jī)處理能力的不斷提升，開展運(yùn)載火箭飛行模態(tài)辨識(shí)所需的硬件條件已基本成熟，制約該技術(shù)走向?qū)嵱没闹饕系K在于欠缺高效魯棒的辨識(shí)方法。除響應(yīng)信號(hào)的非平穩(wěn)性、辨識(shí)信息的低冗余度、短數(shù)據(jù)甚至實(shí)測(cè)數(shù)據(jù)不完整等一般性問題外，開展運(yùn)載火箭飛行模態(tài)辨識(shí)的主要挑戰(zhàn)表現(xiàn)在以下三個(gè)方面［19］：

（1） 僅輸出。運(yùn)載火箭在飛行狀態(tài)下才會(huì)表現(xiàn)出時(shí)變特性，而飛行狀態(tài)下作用在結(jié)構(gòu)上的激勵(lì)難以觀測(cè)（如氣動(dòng)載荷、發(fā)動(dòng)機(jī)內(nèi)部擾動(dòng)等），施加可控可測(cè)的人工激勵(lì)難度較高。因此，解決運(yùn)載火箭飛行模態(tài)辨識(shí)問題需要首先解決僅輸出辨識(shí)問題，即在激勵(lì)未知的情況下，僅基于響應(yīng)數(shù)據(jù)獲取運(yùn)載火箭的時(shí)變模態(tài)。

（2） 遞推。為實(shí)現(xiàn)在線監(jiān)測(cè)與實(shí)時(shí)控制，需要對(duì)運(yùn)載火箭的時(shí)變特性進(jìn)行快速獲取，這就要求辨識(shí)方法具有可遞推性，即在采集得到當(dāng)前時(shí)刻新數(shù)據(jù)后，只基于新數(shù)據(jù)所提供的信息去修改原來估計(jì)出的模態(tài)參數(shù)，而不是基于所有數(shù)據(jù)重新進(jìn)行一輪模態(tài)參數(shù)估計(jì)。遞推辨識(shí)過程中每一步的計(jì)算量和所占用的內(nèi)存空間較小，對(duì)計(jì)算機(jī)的處理能力和存儲(chǔ)空間要求較低，可保證數(shù)據(jù)處理過程跟得上數(shù)據(jù)采集過程。

（3） 復(fù)雜噪聲。由于激勵(lì)未知，從處理實(shí)測(cè)響應(yīng)數(shù)據(jù)的難易程度來看，運(yùn)載火箭真實(shí)的發(fā)射與運(yùn)行工況遠(yuǎn)比地面振動(dòng)試驗(yàn)更為惡劣，需要經(jīng)歷復(fù)雜的噪聲環(huán)境。例如，當(dāng)環(huán)境激勵(lì)較弱時(shí)，箭體結(jié)構(gòu)彈性振動(dòng)幅值較小，可能導(dǎo)致實(shí)測(cè)響應(yīng)信號(hào)的信噪比較低，甚至淹沒在環(huán)境噪聲中。動(dòng)量輪等姿控器件、發(fā)動(dòng)機(jī)內(nèi)部燃燒、液體晃動(dòng)、旋轉(zhuǎn)部件等因素會(huì)產(chǎn)生諧波激勵(lì)，導(dǎo)致作用在箭體上的激勵(lì)為非白噪聲，基于有色噪聲激勵(lì)下的響應(yīng)信號(hào)進(jìn)行時(shí)變模態(tài)辨識(shí)常常會(huì)出現(xiàn)虛假模態(tài)。此外，當(dāng)實(shí)測(cè)響應(yīng)信號(hào)中含有脈沖噪聲等不滿足高斯分布的測(cè)量噪聲時(shí)，也會(huì)嚴(yán)重影響數(shù)據(jù)質(zhì)量，甚至導(dǎo)致辨識(shí)結(jié)果失真。

綜上所述，發(fā)展運(yùn)載火箭時(shí)變模態(tài)僅輸出遞推辨識(shí)方法，在同時(shí)滿足“僅輸出”和“遞推”的基礎(chǔ)上兼顧解決“復(fù)雜噪聲”問題，是運(yùn)載火箭飛行模態(tài)辨識(shí)技術(shù)最終走向?qū)嵱没谋亟?jīng)之路。

2 僅輸出遞推辨識(shí)方法

2.1 時(shí)變自回歸滑動(dòng)平均模型

由于運(yùn)載火箭自身顯著的時(shí)變特性，其響應(yīng)信號(hào)具有非平穩(wěn)性，時(shí)域信號(hào)隨時(shí)間的變化規(guī)律是運(yùn)載火箭結(jié)構(gòu)系統(tǒng)時(shí)變特性的最直接體現(xiàn)。相較于頻域辨識(shí)方法，時(shí)域辨識(shí)方法無需將非平穩(wěn)的響應(yīng)信號(hào)變換至頻域進(jìn)行分析，可有效避免數(shù)據(jù)變換過程中的各類誤差［19］。此外，頻域辨識(shí)方法難以滿足遞推辨識(shí)需求，因此，目前常用的僅輸出遞推辨識(shí)方法主要基于狀態(tài)空間模型或時(shí)間序列模型發(fā)展而來［19，23］。本節(jié)主要介紹時(shí)間序列模型中的時(shí)變自回歸滑動(dòng)平均（TARMA）模型，進(jìn)而給出基于TARMA模型的僅輸出遞推辨識(shí)方法。

作為基于非平穩(wěn)時(shí)間序列建立起來的離散模型，TARMA 模型既能反映時(shí)間序列自身數(shù)據(jù)結(jié)構(gòu)及其變化規(guī)律，也能反映產(chǎn)生該時(shí)間序列的待辨識(shí)系統(tǒng)的固有特性。令na 和nc 分別表示自回歸和滑動(dòng)平均階數(shù)，則TARMA 模型可寫為如下形式［19，23］：

式中 x [ t ] 表示維度為k 的離散非平穩(wěn)響應(yīng)信號(hào)；e [ t ] 表示均值為零、協(xié)方差為Σ [ t ] 的不相關(guān)殘差序列，即滿足正態(tài)獨(dú)立分布N （ 0，Σ [ t ] ）；Ai [ t ] 和Ci [ t ] 分別表示TARMA 模型的自回歸和滑動(dòng)平均系數(shù)矩陣。

基于“時(shí)間凍結(jié)”假設(shè)，可得每一時(shí)刻TARMA模型的功率譜密度（Power Spectral Density，PSD）函數(shù)為［23?25］：

式中 ω 表示頻率，單位為rad s；Δt 為采樣間隔；j 為虛數(shù)單位；（ ? ） H 表示矩陣的Hermitian 轉(zhuǎn)置。

每一時(shí)刻TARMA模型的模態(tài)參數(shù)可通過求解如下特征值問題獲得［23?24］：

（ pi [ t ] I - A [ t ]） vi [ t ]= 0， i= 1，2，…，kna? ?（3）

式中 pi [ t ]和vi [ t ]=[ pi [ t ]-na LTi，…，pi [ t ]-1 LTi]T分別表示伴隨矩陣A [ t ] 的特征值和特征向量，其中Li 為模態(tài)振型向量；I 為單位矩陣。伴隨矩陣A [ t ]由自回歸系數(shù)矩陣構(gòu)成，具體如下［23?24，26］：

更進(jìn)一步，可得TARMA模型的固有頻率和阻尼比分別如下：

時(shí)間序列x [ t ] 的取值大小與先后順序反映了待辨識(shí)系統(tǒng)的固有特性，因此，TARMA 模型參數(shù)Ai [ t ]，Ci [ t ] 和Σ [ t ] 中也蘊(yùn)含了該系統(tǒng)固有特性的信息，這一過程可歸納為信息的凝聚性，即大量數(shù)據(jù)所蘊(yùn)含的信息被凝聚為少數(shù)幾個(gè)模型參數(shù)［23］。顯然，在估計(jì)得到TARMA 模型參數(shù)后，即可基于上述過程獲取待辨識(shí)系統(tǒng)的模態(tài)參數(shù)。

在TARMA模型參數(shù)的估計(jì)過程中，自回歸和滑動(dòng)平均階數(shù)會(huì)直接影響估計(jì)精度及其計(jì)算復(fù)雜度。對(duì)于TARMA模型參數(shù)的遞推估計(jì)過程，可根據(jù)殘差特性來選擇模型階數(shù)，即通過計(jì)算殘差平方和與序列平方和的比值來刻畫TARMA模型對(duì)信號(hào)的擬合誤差［23，27］。一般地，較高的模型階數(shù)能夠更好地對(duì)信號(hào)進(jìn)行擬合，有助于減小擬合誤差，而較低的模型階數(shù)則意味著較低的計(jì)算復(fù)雜度和較好的計(jì)算效率。因此，模型階數(shù)的選擇過程往往需要在擬合精度與計(jì)算效率之間做折中。

2.2 辨識(shí)方法

將式（1）中的TARMA模型改寫為偽線性回歸形式［23?25］，如下式所示：

x [ t ]= w [ t ]TΨ [ t ]+ e [ t ] （6）

其中參數(shù)矩陣w [ t ] 和回歸向量Ψ [ t ] 分別具有如下形式：

在TARMA 模型參數(shù)的估計(jì)過程中，需要對(duì)舊數(shù)據(jù)進(jìn)行遺忘，降低舊數(shù)據(jù)對(duì)當(dāng)前估計(jì)結(jié)果的影響，以有效跟蹤待辨識(shí)系統(tǒng)的時(shí)變特性。借助指數(shù)加權(quán)遺忘機(jī)制，引入遺忘因子λ（0 < λ ≤ 1）對(duì)時(shí)刻τ 下的數(shù)據(jù)進(jìn)行加權(quán)（權(quán)值為λt - τ），以實(shí)現(xiàn)對(duì)舊數(shù)據(jù)的遺忘。定義指數(shù)加權(quán)最小二乘費(fèi)用函數(shù)如下［23］：

式中表示Euclidean范數(shù)。

式（9）的解為：

式中 diag{ ? }表示以{ ? }中元素為對(duì)角元素的對(duì)角矩陣。

定義協(xié)方差矩陣P [ t ]如下：

則進(jìn)一步可得：

P [ t ]-1 = λP [ t - 1 ]-1 + Ψ [ t ] Ψ [ t ]T （13）

根據(jù)矩陣逆定理，式（13）可改寫為：

將式（14）代入式（10）可得：

形式：

則式（15）中t 時(shí)刻參數(shù)矩陣的估計(jì)值w? [ t ] 可改寫為：

綜合式（14）～（19），可將指數(shù)加權(quán)遞推偽線性回歸TARMA模型參數(shù)估計(jì)方法［23］總結(jié)如下：

在估計(jì)得到參數(shù)矩陣后，基于式（2）可計(jì)算時(shí)變PSD 函數(shù)，基于式（3）～（5）可計(jì)算時(shí)變模態(tài)參數(shù)。需要說明的是，該方法的初始化條件為w? [ 0 ]T = 0和P [ 0 ]= αΙ，其中系數(shù)α 一般選取較大的正數(shù)，例如α= 104［ 23，27］，旨在保證初始狀態(tài)時(shí)協(xié)方差矩陣足夠大。由于指數(shù)加權(quán)遺忘機(jī)制的引入，α 對(duì)估計(jì)結(jié)果的影響會(huì)隨著數(shù)據(jù)長(zhǎng)度的增加而逐漸消失。一般地，對(duì)漸變系統(tǒng)進(jìn)行辨識(shí)時(shí)通常選擇遺忘因子λ 為常數(shù)，對(duì)于突變系統(tǒng)則需選擇自適應(yīng)的遺忘因子。

3 時(shí)變模態(tài)辨識(shí)結(jié)果

3.1 發(fā)射準(zhǔn)備階段

本節(jié)基于谷神星一號(hào)運(yùn)載火箭的慣性器件遙測(cè)數(shù)據(jù)進(jìn)行時(shí)變模態(tài)辨識(shí)。選用運(yùn)載火箭橫向和法向速度增量的遙測(cè)數(shù)據(jù)，開展基于多分量實(shí)測(cè)非平穩(wěn)振動(dòng)信號(hào)的數(shù)據(jù)快速處理與時(shí)變模態(tài)辨識(shí)。在點(diǎn)火發(fā)射前，運(yùn)載火箭豎立于發(fā)射臺(tái)上，可近似視為懸臂梁邊界條件，由于發(fā)動(dòng)機(jī)尚未點(diǎn)火，此時(shí)其模態(tài)參數(shù)應(yīng)保持不變。為驗(yàn)證該階段運(yùn)載火箭的模態(tài)參數(shù)是否發(fā)生變化，在選擇TARMA 模型階數(shù)na = 20 和nc = 2 的情況下，進(jìn)一步選取遺忘因子λ = 0.984，基于點(diǎn)火發(fā)射前90 s 內(nèi)的遙測(cè)數(shù)據(jù)對(duì)運(yùn)載火箭的固有頻率和PSD 函數(shù)進(jìn)行辨識(shí)，獲取地面風(fēng)載作用下運(yùn)載火箭模態(tài)參數(shù)隨時(shí)間的變化規(guī)律。

發(fā)射準(zhǔn)備階段固有頻率的辨識(shí)結(jié)果如圖3所示，其中，紅色圓點(diǎn)為固有頻率的辨識(shí)結(jié)果，綠色虛線為有限元模型計(jì)算得到的固有頻率參考值。由圖3可知，在發(fā)射前90 s內(nèi)，運(yùn)載火箭的前三階固有頻率基本保持不變，這與實(shí)際工況是一致的。發(fā)射準(zhǔn)備階段PSD函數(shù)的辨識(shí)結(jié)果如圖4所示，顯然圖中脊線均呈現(xiàn)水平狀態(tài)，且與圖3中固有頻率的辨識(shí)結(jié)果完全對(duì)應(yīng)，這也進(jìn)一步表明該階段運(yùn)載火箭的模態(tài)參數(shù)基本保持不變。需要說明的是，從發(fā)射前60 s開始，運(yùn)載火箭的電爆閥、隔離閥等相繼開始工作，相應(yīng)的起爆和開啟動(dòng)作會(huì)引起較大的沖擊振動(dòng)，因此圖4中相繼出現(xiàn)了多條豎直方向的高亮譜線，也導(dǎo)致了圖3中相對(duì)應(yīng)時(shí)刻的固有頻率辨識(shí)結(jié)果較差。

將圖3中的固有頻率辨識(shí)結(jié)果與有限元計(jì)算結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)前兩階固有頻率吻合良好，但有限元計(jì)算出的第三階固有頻率值較辨識(shí)值偏高。前期相關(guān)研究已表明，基于Ares I?X和Vega火箭的有限元模型計(jì)算得到的固有頻率常常與其辨識(shí)值存在偏差［3?5］，這也進(jìn)一步說明了有限元模型與真實(shí)工作狀態(tài)下的運(yùn)載火箭仍存在差異。因此，為更好地對(duì)谷神星一號(hào)運(yùn)載火箭的結(jié)構(gòu)動(dòng)力學(xué)特性進(jìn)行表征，尤其是對(duì)發(fā)射準(zhǔn)備階段第三階模態(tài)進(jìn)行表征，仍需要對(duì)其有限元模型進(jìn)行修正，以提高模型精度。

3.2 飛行階段

如圖2所示，谷神星一號(hào)運(yùn)載火箭芯一級(jí)所配備的固體發(fā)動(dòng)機(jī)的燃燒時(shí)間為74 s。因此，基于運(yùn)載火箭點(diǎn)火發(fā)射后0～74 s的遙測(cè)數(shù)據(jù)可對(duì)一級(jí)飛行階段的模態(tài)參數(shù)進(jìn)行辨識(shí)。固有頻率及其對(duì)應(yīng)PSD函數(shù)的辨識(shí)結(jié)果如圖5所示，其中，紅色圓點(diǎn)為固有頻率的辨識(shí)結(jié)果，背景色為PSD函數(shù)的辨識(shí)結(jié)果。由圖5可知，固有頻率的辨識(shí)結(jié)果與PSD函數(shù)的脊線吻合良好，清晰表明前四階固有頻率均隨燃料消耗而不斷升高，驗(yàn)證了本文辨識(shí)方法對(duì)運(yùn)載火箭時(shí)變模態(tài)進(jìn)行跟蹤的有效性。需要說明的是，由于發(fā)射40 s之后第二階模態(tài)未被較好地激發(fā)出來，導(dǎo)致第二階固有頻率的辨識(shí)精度較差，但PSD函數(shù)的脊線依舊能夠刻畫出第二階固有頻率的大概變化規(guī)律。

此外，圖5中也給出了芯一級(jí)固體發(fā)動(dòng)機(jī)滿載和空載情況下基于有限元模型計(jì)算得到的固有頻率參考值，分別對(duì)應(yīng)實(shí)際飛行的0時(shí)刻和74 s時(shí)刻，在圖5中用黑色菱形表示。顯然，飛行模態(tài)的辨識(shí)結(jié)果和有限元計(jì)算結(jié)果吻合良好，說明本文辨識(shí)方法能夠?qū)\(yùn)載火箭的時(shí)變模態(tài)進(jìn)行準(zhǔn)確辨識(shí)。需要指出的是，在火箭發(fā)動(dòng)機(jī)滿載和空載情況下進(jìn)行有限元分析相對(duì)容易，但在有限元模型中考慮燃料消耗所導(dǎo)致的質(zhì)量連續(xù)時(shí)變特性比較困難，難以獲取每一時(shí)刻下運(yùn)載火箭的模態(tài)信息。因此，開展飛行模態(tài)辨識(shí)可有效彌補(bǔ)有限元分析的這一局限性，能夠連續(xù)獲取每一時(shí)刻下運(yùn)載火箭的模態(tài)參數(shù)，為姿控系統(tǒng)設(shè)計(jì)和天地一致性研究提供重要支撐。

如圖2所示，谷神星一號(hào)運(yùn)載火箭芯二級(jí)所配備的固體發(fā)動(dòng)機(jī)的燃燒時(shí)間為70 s。因此，基于運(yùn)載火箭點(diǎn)火發(fā)射后75～145 s的遙測(cè)數(shù)據(jù)可對(duì)二級(jí)飛行階段的模態(tài)參數(shù)進(jìn)行辨識(shí)。固有頻率及其對(duì)應(yīng)PSD函數(shù)的辨識(shí)結(jié)果如圖6所示。由圖6可知，固有頻率的辨識(shí)結(jié)果與PSD函數(shù)的脊線吻合良好，表明了前兩階固有頻率均隨燃料消耗而不斷升高的規(guī)律。與圖5中的辨識(shí)結(jié)果相比，圖6中模態(tài)辨識(shí)結(jié)果的精度相對(duì)較差，主要是由于二級(jí)飛行階段運(yùn)載火箭的飛行速度更大，所處噪聲環(huán)境更為復(fù)雜，一定程度上增加了僅基于響應(yīng)數(shù)據(jù)進(jìn)行模態(tài)辨識(shí)的難度。此外，圖6中25～30 Hz頻帶內(nèi)雖然模糊出現(xiàn)了PSD函數(shù)的脊線，但在辨識(shí)過程中通過篩除阻尼比過大的極點(diǎn)，該頻帶內(nèi)并未出現(xiàn)大范圍的虛假模態(tài)。

類似地，圖6中也給出了芯二級(jí)固體發(fā)動(dòng)機(jī)滿載和空載情況下基于有限元模型計(jì)算得到的固有頻率參考值，分別對(duì)應(yīng)實(shí)際飛行的75 s時(shí)刻和145 s時(shí)刻。由圖6可知，飛行模態(tài)的辨識(shí)結(jié)果和有限元計(jì)算結(jié)果基本吻合，說明本文辨識(shí)方法依舊能夠?qū)\(yùn)載火箭二級(jí)飛行階段的時(shí)變模態(tài)進(jìn)行跟蹤。

目前，通過有限元分析和地面振動(dòng)試驗(yàn)，一般在運(yùn)載火箭發(fā)射之前就給定了其整個(gè)工作周期的模態(tài)特性，但兩者均難以連續(xù)獲取運(yùn)載火箭的時(shí)變模態(tài)參數(shù)。有限元方法難以對(duì)連續(xù)時(shí)變結(jié)構(gòu)系統(tǒng)進(jìn)行仿真分析，只能對(duì)特定時(shí)刻（或特定燃料消耗情況）下的結(jié)構(gòu)系統(tǒng)進(jìn)行模擬。地面振動(dòng)試驗(yàn)無法避免進(jìn)行多次重復(fù)試驗(yàn)的繁瑣流程，往往采用秒時(shí)刻凍結(jié)進(jìn)行多次時(shí)不變結(jié)構(gòu)系統(tǒng)的模態(tài)試驗(yàn)，再利用插值法獲取全過程的模態(tài)信息，費(fèi)時(shí)費(fèi)力。相較而言，借助時(shí)變模態(tài)辨識(shí)技術(shù)，充分利用天地試驗(yàn)數(shù)據(jù)，可準(zhǔn)確快速地獲取真實(shí)工作狀態(tài)下運(yùn)載火箭的時(shí)變動(dòng)力學(xué)特性，且滿足天地一致性要求，能夠有效服務(wù)運(yùn)載火箭的有限元模型修正和姿控系統(tǒng)設(shè)計(jì)。

4 結(jié)? 論

本文系統(tǒng)性介紹了運(yùn)載火箭時(shí)變模態(tài)辨識(shí)的研究意義和發(fā)展現(xiàn)狀，闡述了現(xiàn)階段開展時(shí)變模態(tài)辨識(shí)所面臨的難題與挑戰(zhàn)，發(fā)展了指數(shù)加權(quán)遞推偽線性回歸TARMA模型參數(shù)估計(jì)方法，能夠在激勵(lì)未知的情況下，僅基于響應(yīng)數(shù)據(jù)快速辨識(shí)運(yùn)載火箭的時(shí)變模態(tài)。在工程應(yīng)用方面，圍繞谷神星一號(hào)運(yùn)載火箭的時(shí)變模態(tài)辨識(shí)問題，基于上述辨識(shí)方法完成了飛行遙測(cè)數(shù)據(jù)的快速處理，準(zhǔn)確辨識(shí)了發(fā)射準(zhǔn)備階段和飛行階段的關(guān)鍵低階時(shí)變模態(tài)，有效服務(wù)了運(yùn)載火箭的有限元模型修正和姿控系統(tǒng)設(shè)計(jì)。

參考文獻(xiàn)：

［1］????? 孟光， 周徐斌， 苗軍. 航天重大工程中的力學(xué)問題［J］. 力學(xué)進(jìn)展， 2016， 46： 267-322.

MENG Guang， ZHOU Xubin， MIAO Jun. Mechanical problems in momentous projects of aerospace engineering［J］. Advances in Mechanics， 2016， 46： 267-322.

［2］????? 南宮自軍， 戴新進(jìn)， 王亮， 等. 航天飛行器飛行模態(tài)辨識(shí)技術(shù)應(yīng)用及展望［J］. 強(qiáng)度與環(huán)境， 2017， 44（1）： 20-26.

NANGONG Zijun， DAI Xinjin， WANG Liang， et al. Application and prospective outlook of the aerospace vehicles operational mode identification［J］. Structure & Environment Engineering， 2017， 44（1）： 20-26.

［3］????? 馬慶港， 王紫揚(yáng)， 康杰， 等. 運(yùn)載火箭結(jié)構(gòu)工作模態(tài)參數(shù)辨識(shí)研究綜述［J］. 彈箭與制導(dǎo)學(xué)報(bào)， 2022， 42（1）： 60-70.

MA Qinggang， WANG Ziyang， KANG Jie， et al. Overview of operational modal parameter identification for structures of rocket launchers［J］. Journal of Projectiles， Rockets， Missiles and Guidance， 2022， 42（1）： 60-70.

［4］????? Bartkowicz T J， James G H. Ares I-X in-flight modal identification［C］. 52nd AIAA/ASME/ASCE/AHS/ASC Structare， Structural Dynamics and Materials Conference. Denver， Colorado， 2011.

［5］????? De Vivo A， Brutti C， Leofanti J L. Vega in-flight modal identification with the operational modal analysis technique［J］. Journal of Spacecraft and Rockets， 2014， 51（5）： 1464-1473.

［6］????? 王亮， 張妍， 蔡毅鵬， 等. 特征時(shí)段選擇對(duì)飛行模態(tài)辨識(shí)結(jié)果的影響［J］. 噪聲與振動(dòng)控制， 2017， 37（1）： 128-131.

WANG Liang， ZHANG Yan， CAI Yipeng， et al. Influence of the measure feature phase selection on operational modal identification of aircrafts［J］. Noise and Vibration Control， 2017， 37（1）： 128-131.

［7］????? Khoshnood A M， Roshanian J， Jafari A A， et al. Simultaneous estimation of two bending vibration frequencies for attitude control of a flexible launch vehicle［J］. Journal of Systems and Control Engineering， 2009， 223： 721-726.

［8］????? Jang J W， Hall R， Bedrossian N， et al. Ares-I bending filter design using a constrained optimization approach［C］//AIAA Guidance， Navigation and Control Conference and Exhibit. Honolulu， Hawaii， 2008： 6289.

［9］????? Tobbe P A， Matras A L， Wilson H E. Modeling and simulation of variable mass， flexible structures［C］//AIAA Modeling and Simulation Technologies Conference. Chicago， Illinois， 2009： 6023.

［10］??? Fransen S， Rixen D， Henriksen T， et al. On the operational modal analysis of solid rocket motors［C］//The 28th IMAC， International Modal Analysis Conference. Jacksonville， Florida， 2010： 453-463.

［11］??? Goursat M， D?hler M， Mevel L， et al. Crystal clear SSI for operational modal analysis of aerospace vehicles［C］//The 28th IMAC， International Modal Analysis Conference. Jacksonville， Florida， 2010： 1421-1430.

［12］??? Manzato S， Peeters B， Debille J. Tracking the evolution of modal properties of a solid propellant launcher during static firing test［C］//The 30th IMAC， International Modal Analysis Conference. Jacksonville， Florida， 2012： 559-570.

［13］??? James G H. Development of operational modal analysis techniques for launch data［J］. Advances in the Astronautical Sciences， 2013， 147： 209-230.

［14］??? James G， Kaouk M， Cao T. Progress in operational analysis of launch vehicles in nonstationary flight［C］//The 31st IMAC， International Modal Analysis Conference. Garden Grove， California， 2013： 59-75.

［15］??? 張家雄， 何詠梅， 張華山， 等. 基于PolyMAX法的飛行器工作模態(tài)分析技術(shù)與應(yīng)用［J］. 航天器環(huán)境工程， 2015， 32（1）： 28-33.

ZHANG Jiaxiong， HE Yongmei， ZHANG Huashan， et al. Operational modal analysis technology based on PolyMAX method and its applications for flying vehicles［J］. Spacecraft Environment Engineering， 2015， 32（1）： 28-33.

［16］??? 王亮， 張妍， 周曉麗， 等. 基于ARMA-NExT和穩(wěn)定圖方法的飛行器工作模態(tài)指示研究［J］. 動(dòng)力學(xué)與控制學(xué)報(bào)， 2016， 14（3）： 258-262.

WANG Liang， ZHANG Yan， ZHOU Xiaoli， et al. Study on the indication of the aircraft operational mode based on ARMA-NExT and stabilization diagram［J］. Journal of Dynamics and Control， 2016， 14（3）： 258-262.

［17］??? 王亮， 蔡毅鵬， 朱辰， 等. 基于ARMA-NExT的飛行器工作模態(tài)辨識(shí)技術(shù)研究［J］. 導(dǎo)彈與航天運(yùn)載技術(shù)， 2017（1）： 18-21.

WANG Liang， CAI Yipeng， ZHU Chen， et al. Operational mode identification of the aircraft based on ARMA-NExT［J］. Missiles and Space Vehicles， 2017（1）： 18-21.

［18］??? 董嚴(yán)， 付小燕， 丁志偉. 基于多測(cè)點(diǎn)數(shù)據(jù)的火箭飛行模態(tài)參數(shù)識(shí)別方法［J］. 固體火箭技術(shù)， 2018， 41（4）： 520-523.

DONG Yan， FU Xiaoyan， DING Zhiwei. Rocket flight modal identification method based on data of multi-measure points［J］. Journal of Solid Rocket Technology， 2018， 41（4）： 520-523.

［19］??? 馬志賽， 丁千， 劉莉， 等. 線性時(shí)變結(jié)構(gòu)模態(tài)參數(shù)時(shí)域辨識(shí)方法的研究進(jìn)展［J］. 機(jī)械工程學(xué)報(bào)， 2018， 54（23）： 137-159.

MA Zhisai， DING Qian， LIU Li， et al. Research progress on time-domain modal parameter estimation methods for linear time-varying structures［J］. Journal of Mechanical Engineering， 2018， 54（23）： 137-159.

［20］??? Ma Z S， Ding Q， Tang Y. Operational modal analysis of a liquid-filled cylindrical structure with decreasing filling mass by multivariate stochastic parameter evolution methods［J］. International Journal of Mechanical Sciences， 2020， 172： 105420.

［21］??? Ma Z S， Li L Q， Ding Q. Multivariate recursive bayesian linear regression and its applications to output-only identification of time-varying mechanical systems［J］. Journal of Vibration and Control， 2021， 27（11-12）： 1395-1406.

［22］??? 余磊， 劉莉， 崔穎， 等. 一種運(yùn)載火箭時(shí)變結(jié)構(gòu)模態(tài)參數(shù)辨識(shí)的確定性演化方法［J］. 宇航學(xué)報(bào)， 2020， 41（4）： 379-388.

YU Lei， LIU Li， CUI Ying， et al. A time-varying structure modal parameter estimation method in deterministic evolution for launch vehicle［J］. Journal of Astronautics， 2020， 41（4）： 379-388.

［23］??? 馬志賽. 線性時(shí)變結(jié)構(gòu)模態(tài)參數(shù)僅輸出遞推辨識(shí)方法研究［D］. 北京： 北京理工大學(xué)， 2017.

MA Zhisai. Output-only modal parameter recursive estimation for linear time-varying structures［D］. Beijing： Beijing Institute of Technology， 2017.

［24］??? Ma Z S， Liu L， Zhou S D， et al. Parametric output-only identification of time-varying structures using a kernel recursive extended least squares TARMA approach［J］. Mechanical Systems and Signal Processing， 2018， 98： 684-701.

［25］??? Ma Z S， Ding Q， Zhou S D. Novel adaptive methods for output-only recursive identification of time-varying systems subject to gross errors［J］. Journal of Vibration and Control， 2020， 26（5-6）： 306-317.

［26］??? Bertha M， Golinval J C. Identification of non-stationary dynamical systems using multivariate ARMA models［J］. Mechanical Systems and Signal Processing， 2017， 88： 166-179.

［27］??? Spiridonakos M D， Fassois S D. Parametric identification of a time-varying structure based on vector vibration response measurements［J］. Mechanical Systems and Signal Processing， 2009， 23（6）： 2029-2048.

Time-varying modal parameter estimation of the CERES-1 launch vehicle

Abstract： It is usually difficult to establish the dynamic model of a launch vehicle that accurately describes its time-varying characteristics. Therefore modal identification techniques are particularly necessary to obtain the time-varying dynamic characteristics of launch vehicles under flight conditions. Aiming at the problem of in-flight modal identification of launch vehicles， an output-only recursive identification method based on the time-dependent autoregressive moving average model is developed by using exponentially weighted mechanisms to track the time-varying characteristics. Without measuring the natural excitation forces， the proposed method can accurately and quickly identify the time-varying modal parameters of launch vehicles by exclusively using the measured response signals. Taking the CERES-1 launch vehicle as an example， time-varying modal parameters before liftoff and during the flight phase are accurately estimated by processing the flight telemetry data. Identification results are consistent with the variation of the finite element analysis results， demonstrating the high achievable accuracy of the proposed method. The proposed in-flight modal identification method can obtain the full-cycle modal information of launch vehicles， which meets the engineering requirements for the finite element model updating and attitude control system design.

Key words： launch vehicle； time-varying modal parameter； CERES-1； recursive identification