程 興，高 晨，賈大玲，李 君，王紫揚

（北京宇航系統(tǒng)工程研究所，北京，100076）

0 引 言

近年來，空間發(fā)射任務(wù)逐漸多樣化、復(fù)雜化，因此對運載火箭提出了更高的要求[1～4]。彈道、姿控和制導(dǎo)等專業(yè)的設(shè)計和集成對于運載火箭非常重要，吸引了眾多學(xué)者的關(guān)注并且取得了許多成果[5～7]。

運載火箭的工程化設(shè)計一般采用“獨立設(shè)計+聯(lián)合分析仿真”的模式，即彈道、姿控、制導(dǎo)、載荷、結(jié)構(gòu)、動力等專業(yè)根據(jù)標準彈道各自獨立設(shè)計，各個專業(yè)根據(jù)自身需求進行偏差條件下的包絡(luò)設(shè)計。這樣的方法可以有效保證設(shè)計效率，但是由于各個專業(yè)進行設(shè)計時使用的動力學(xué)模型不同，參數(shù)偏差的取值和考慮因素各有側(cè)重，難以準確分析多項偏差組合對運載火箭姿態(tài)、軌跡、載荷等參數(shù)的影響。因此，如何準確、高效進行總體回路的集成建模，量化分析各種因素對運載火箭質(zhì)量特性以及飛行動力學(xué)的影響，是現(xiàn)階段運載火箭設(shè)計中急需解決的問題。

另一方面，運載火箭是復(fù)雜的系統(tǒng)工程，影響因素多且相互間耦合，尤其是一級飛行段（含助推飛行段）經(jīng)歷惡劣氣動環(huán)境，缺乏對其飛行力學(xué)特征進行量化分析的手段，難以開展動力學(xué)關(guān)鍵參數(shù)的天地一致性確認及非故障狀態(tài)的偏差辨識。從飛行結(jié)果來看，盡管飛行過程中姿態(tài)穩(wěn)定并高精度地將有效載荷送入預(yù)定軌道，但各子級的姿態(tài)及軌道參數(shù)乃至飛行時間均較標稱值存在偏差，尤其是一級飛行段的參數(shù)偏差更為突出。若能有效辨識飛行偏差背后的物理因素，據(jù)此開展設(shè)計改進與優(yōu)化，提升飛行與設(shè)計的一致性，將能夠提高設(shè)計精細化水平，從而有利于提高飛行可靠性及潛力挖掘，需采用有效的辨識算法進行參數(shù)辨識和一致性分析[8～10]。

目前，常用的辨識算法包括卡爾曼濾波、極大似然估計、分割算法等。但是卡爾曼濾波算法對模型精確程度要求較高；基于極大似然估計的辨識方法一般要求解雅克比矩陣，算法計算量較大，工程實現(xiàn)較為困難，而且算法在系統(tǒng)非線性較強時可能失效，需要克服矩陣奇異的問題。運載火箭飛行環(huán)境干擾因素復(fù)雜繁多，因此所設(shè)計算法的魯棒性尤為重要。粒子群算法、遺傳算法等算法不依賴于初始值，具有較好的尋優(yōu)能力，而且計算量較小，具有較好的工程應(yīng)用價值[11～14]。

針對速度偏差、姿態(tài)偏差及關(guān)機時間偏差幅值較大的狀況，以實現(xiàn)火箭動力學(xué)天地一致性分析，本文首先提出源變量動力學(xué)建模技術(shù)，并據(jù)此實現(xiàn)總體回路集成仿真；在此基礎(chǔ)上，利用粒子群算法具有較快的收斂速度，遺傳算法具有較強全局搜索能力的特點，采用分層結(jié)構(gòu)設(shè)計混合優(yōu)化算法進行參數(shù)辨識，保證系統(tǒng)具有較快的搜索速度和較強的全局搜索能力，從而實現(xiàn)運載火箭動力學(xué)參數(shù)的天地一致性分析。最后通過仿真驗證了所提方法的有效性。

1 基于源變量的總體回路集成仿真技術(shù)

液體運載火箭普遍采用彈道、姿控、制導(dǎo)、原始數(shù)據(jù)計算、分離等多專業(yè)獨立設(shè)計的模式[15,16]。該設(shè)計模式下各專業(yè)以標準彈道為起點各自獨立設(shè)計并按專業(yè)準則進行偏差包絡(luò)設(shè)計，既保證了設(shè)計效率，也保證了系統(tǒng)適應(yīng)有限偏差的能力。該模式下，各專業(yè)使用的動力學(xué)數(shù)學(xué)模型各有差異、參數(shù)偏差取值的考慮因素各有側(cè)重，這使得針對特定偏差的分析具有周期長、匹配性低的特點，難以量化分析單項或多項偏差對飛行姿態(tài)、飛行軌跡、飛行載荷的綜合影響，譬如推進劑流量下降后將影響推進劑質(zhì)量及晃動特性，影響全箭質(zhì)量分布及質(zhì)心、轉(zhuǎn)動慣量，過載及速度等。另外，彈道設(shè)計采用瞬時平衡準則，忽略姿態(tài)控制的動態(tài)響應(yīng)過程，姿態(tài)動力學(xué)建模以標準彈道為基礎(chǔ)，通常采用質(zhì)量瞬時凝固法等，即專業(yè)間的模型缺乏一致性，以提高各專業(yè)間參數(shù)的匹配性、量化分析參數(shù)偏差對系統(tǒng)影響為目標，提出了基于源變量的動力學(xué)集成建模技術(shù)。

為了精確量化分析包括發(fā)動機推力異常、推進劑流量異常等偏差或故障對全箭質(zhì)量特性及飛行動力學(xué)的影響，提出源變量建模思路：以發(fā)動機安裝角、方位角、擺角、推進劑秒耗量、加注量、全箭分布質(zhì)量等參數(shù)為基本變量（即源變量），通過在線計算推進劑流量及剩余量、全箭質(zhì)量特性（質(zhì)量、質(zhì)心、轉(zhuǎn)動慣量）、發(fā)動機推力及推進劑晃動特性、結(jié)構(gòu)彈性振動及載荷分布，時時迭代計算全箭飛行軌跡和姿態(tài)動力學(xué)特性，實現(xiàn)總體-彈道-姿控-制導(dǎo)-動力-推進劑晃動-載荷計算的集成建模，并以此推動一體化仿真。

相對文獻[15]和文獻[16]中提出的運載火箭六自由度仿真方法不同，本仿真模型具有3個特點：

a）發(fā)動機推力P為安裝角、安裝象限、擺角的矢量化模型。

式中ckε為第k臺發(fā)動機的安裝角；ckμ為第k臺發(fā)動機的象限角；ckδ為第k臺發(fā)動機的擺動指令角。

該建模方式下，通過模型變換便能簡潔、高效地消除剛-晃-彈動力學(xué)模型中的代數(shù)環(huán)問題而顯著提高仿真效率及精度。這里不作動力學(xué)模型的逐項詳細分析，而僅列出耦合模型形式。動力學(xué)的完整表達形式及其物理含義詳見文獻[1]。

即有：

進一步變換有：

通過式（5）消除了動力學(xué)模型中的代數(shù)環(huán)問題。

相對于將質(zhì)心動力學(xué)方程建立在發(fā)射慣性系或發(fā)射系的傳統(tǒng)建模方法，基于箭體坐標系的質(zhì)心方程（2）能大大簡化變換矩陣T，并顯著提升計算效率及精度。

c）總體回路專業(yè)深度耦合。

相對基于事先計算推進劑晃動特性、彈性振動特性、全量質(zhì)量特性的傳統(tǒng)運載火箭六自由度仿真，本仿真通過離線-在線組合、參數(shù)間耦合機理的再梳理等方法研究，實現(xiàn)基于分布質(zhì)量的全箭質(zhì)量特性參數(shù)（質(zhì)心、轉(zhuǎn)動慣量、液位高度、增壓壓力）在線計算，推進劑晃動特性參數(shù)的在線計算（晃動質(zhì)量、晃動頻率、晃動質(zhì)心高度），以及飛行軌跡、飛行姿態(tài)、晃動響應(yīng)、彈性響應(yīng)等響應(yīng)參數(shù)的在線計算。

源變量建模及仿真方法優(yōu)化的詳細過程可參見文獻[17]。

源變量仿真（六自由度全量模型）與標準彈道（三自由度質(zhì)心模型）的速度、位置偏差對比見圖1。

圖1 源變量仿真與彈道計算參數(shù)的對比Fig.1 Flight Parameter between Svarsim and Trajectory Design

圖1 的對比結(jié)果表明考慮控制過程的動態(tài)響應(yīng)與僅考慮瞬時平衡的飛行速度略有差異，該偏差在系統(tǒng)設(shè)計允許的誤差范圍內(nèi)，再次驗證彈道計算采用瞬時平衡處理的工程有效性，也驗證了基于源變量的多專業(yè)耦合仿真模型的正確性。

源變量集成仿真實現(xiàn)了多專業(yè)間的參數(shù)緊耦合，能量化分析一個或多個參數(shù)偏差對包括飛行軌跡、飛行姿態(tài)在內(nèi)的飛行品質(zhì)參數(shù)的影響。因此以該模型為數(shù)學(xué)基礎(chǔ)，開展運載火箭動力學(xué)的天地差異性辨識。

2 基于源變量仿真的動力學(xué)天地差異性辨識

下面以全程復(fù)現(xiàn)實際飛行中的軌跡、姿態(tài)、推進劑液位高度等關(guān)鍵飛行參數(shù)為目標，研究運載火箭參數(shù)的天地差異性。

影響飛行姿態(tài)、飛行軌跡及關(guān)機時間的因素多，其中包括各發(fā)動機流量偏差、比沖偏差，推力線橫移，各推進劑貯箱的加注量偏差、結(jié)構(gòu)質(zhì)量偏差，氣動法向力系數(shù)偏差、阻力系數(shù)偏差、壓心偏差，全箭質(zhì)心橫移，發(fā)動機零位及伺服機構(gòu)零位漂移等。除加注量偏差、結(jié)構(gòu)質(zhì)量偏差外，其余參數(shù)偏差均隨時間變化，總計216個待辨識量。

為適應(yīng)待辨識偏差參數(shù)多的狀況，采用粒子群算法-遺傳算法組合的新型搜索算法。

目標函數(shù)為仿真與飛行之間的速度偏差、姿態(tài)角偏差、角速度偏差、控制擺角偏差及助推器質(zhì)量（液位高度）偏差均方和最小，具體為

式中 ki=0,…, 4 分別為速度偏差、姿態(tài)偏差、角速度偏差、控制擺角及液位對應(yīng)質(zhì)量偏差的權(quán)重系數(shù)；為仿真中的發(fā)動機關(guān)機時間；為對應(yīng)飛行值；，分別為各助推器氧化劑、燃燒劑質(zhì)量仿真質(zhì)量與飛行質(zhì)量之間的偏差。

姿態(tài)角偏差、角速度及擺角偏差遠小于速度偏差及質(zhì)量偏差，因此通過 ki來調(diào)節(jié)實現(xiàn)參數(shù)對結(jié)果的均衡化，本次分析中具體取值為 k1=10，k2=200，k3=400，

2.1 遺傳算法通過模擬生物進化論實現(xiàn)參數(shù)的辨識和尋優(yōu)

遺傳算法包括 4個基本元素，即編碼、選擇、交叉和變異，具有較強的全局搜索能力。粒子群算法源于對鳥類捕食行為的研究，具有較快的收斂速度。本文基于遺傳算法和粒子群算法，采用分層設(shè)計的思路，設(shè)計參數(shù)辨識混合算法。算法的底層采用遺傳算法，對系統(tǒng)參數(shù)進行全局搜索和辨識；頂層設(shè)計采用粒子群算法，對遺傳算法得到的最優(yōu)解作為初始種群，從而進行快速的局部搜索，得到參數(shù)辨識的結(jié)果。同時，為了提高粒子群算法的快速性，本文通過設(shè)計動態(tài)慣性權(quán)值對算法進行設(shè)計。定義待辨識系統(tǒng)參數(shù)為 xi，中間代個體定義為 yi，最大迭代次數(shù)定義為 G1，變異概率記為Pm。

參數(shù)辨識底層遺傳算法的流程為：

a）步驟1：通過確定待辨識參數(shù)的可能取值范圍，采用二進制的方法對各個待辨識參數(shù)進行編碼，即：

結(jié)合所有待辨識參數(shù)的可能取值范圍，通過二進制的方法對參數(shù)進行編碼。例如，可以將ix的范圍通過相應(yīng)位數(shù)的二進制串進行描述：

b）步驟2：如果迭代至最大次數(shù)，則算法終止；反之，執(zhí)行步驟3。

c）步驟3：通過式（9）、式（10）選出一定數(shù)目的中間個體，并且進行交叉運算。

d）步驟4：對待優(yōu)化的所有個體以給定概率進行變異操作。

e）步驟5：計算辨識目標函數(shù)。

f）步驟 6：將當代最優(yōu)個體與歷史最優(yōu)個體進行比較并更新。

g）步驟7：判斷算法得到的最優(yōu)解是否滿足停止條件，若滿足，則系統(tǒng)輸出最優(yōu)解，反之則繼續(xù)進行尋優(yōu)。

當采用遺傳算法迭代至代數(shù)G1后，將系統(tǒng)的最優(yōu)解取出，作為粒子群算法的初始值。采用粒子群優(yōu)化算法對遺傳算法得到的結(jié)果進行優(yōu)化。采用粒子群算法迭代至代數(shù)G2之后，若滿足停止條件，則將粒子群優(yōu)化后的結(jié)果作為系統(tǒng)辨識的最優(yōu)結(jié)果輸出；若不滿足停止條件，則用粒子群種群中的 q個粒子隨機替換掉遺傳算法子群中的 q個粒子，對遺傳算法子群再一次進行優(yōu)化求解，不斷循環(huán)，系統(tǒng)滿足停止條件，從而輸出最優(yōu)解。

2.2 采用粒子群算法對系統(tǒng)參數(shù)進行辨識

采用粒子群算法對系統(tǒng)參數(shù)進行辨識時，慣性權(quán)值和參數(shù)收斂速度密切相關(guān)，較大的慣性權(quán)值有利于全局搜索；較小的慣性權(quán)值有利于算法的局部開發(fā)，加速算法的收斂。本文通過設(shè)計動態(tài)慣性權(quán)值來提高粒子收斂速度，通過式（11）、式（12）來更新粒子的位置和速度：

式中 r1，r2為服從均勻分布且相互獨立的隨機數(shù)；g1，g2為給定常數(shù)；vi(j)為待辨識參數(shù)粒子i在第j代的速度；為粒子i在第j代的位置； Pi( j)為粒子i在第j代的最優(yōu)位置為第j代全局最優(yōu)位置。設(shè)動態(tài)慣性權(quán)值?為

可以看出，隨著采用粒子群算法迭代次數(shù)的增多，系統(tǒng)慣性權(quán)值減小，而較小的慣性權(quán)值可以加快算法的收斂速度。

粒子群算法的流程為：

a）步驟1：對粒子群算法的種群進行初始化，即采用遺傳算法的最優(yōu)解作為粒子群算法種群的初始位置，并隨機產(chǎn)生每個粒子的速度。

b）步驟2：判斷系統(tǒng)是否到達最大迭代數(shù)，若到達，則算法終止，反之則執(zhí)行步驟3。

c）步驟3：計算粒子群算法種群每個粒子的適應(yīng)度函數(shù)值。

d）步驟4：更新粒子群算法種群的個體極值和全局極值。

e）步驟5：更新粒子群算法中每個粒子的速度和位置。

f）步驟 6：判斷優(yōu)化得到的最優(yōu)解是否滿足停止準則，若滿足，則輸出最優(yōu)解；反之，跳轉(zhuǎn)到步驟 2繼續(xù)進行參數(shù)優(yōu)化。

2.3 提高辨識效率和有效性處理方式

a）包括增壓壓力在內(nèi)的部分參數(shù)直接采用飛行參數(shù)而非設(shè)計參數(shù)；用于氣動插值的飛行速度、高度參數(shù)也直接用飛行參數(shù)而非在線計算參數(shù)；

b）氣動數(shù)據(jù)偏差包括法向力系數(shù)偏差、壓心系數(shù)偏差，采用乘法形式的修正系數(shù)，為時變參數(shù)，時間及修正值均通過仿真確定，為提高效率，約束時間最小間隔；

c）辨識基于Matlab/Simlink環(huán)境，充分利用高版本支持整數(shù)辨識的功能，而對部分參數(shù)進行整數(shù)規(guī)格化處理，譬如設(shè)定壓心補償系數(shù)=0.9～1.1，遺傳尋優(yōu)中設(shè)定其上限為1100，下限為900，分層值為1，這樣便大幅提升尋優(yōu)效率。

尋優(yōu)辨識之后的飛行速度偏差變化歷程見圖2，姿態(tài)偏差變化歷程見圖3。另外，關(guān)機偏差從飛行相對于標準彈道提前4.73 s降低到當前的0.07 s，即可認為關(guān)機時間復(fù)現(xiàn)飛行狀態(tài)（圖中數(shù)據(jù)已歸一化處理）。

圖2 源變量仿真與飛行之間的速度偏差（辨識結(jié)果）Fig.2 Velocity Bias between Svarsim and Flight (by Identification)

圖3 源變量仿真與飛行之間的姿態(tài)偏差Fig.3 Attitude Bias between Svarsim and Flight (by Identification)

2.4 部分參數(shù)辨識結(jié)果

a）部分助推器流量偏差較大，其中表1為3#助推器的氧化劑、燃燒劑流量偏差結(jié)果（經(jīng)過歸一化處理）。

表1 3#助推器氧化劑、燃燒劑流量偏差Tab.1 3# Booster Flow Deviation of Oxidizer and Fuel

該數(shù)據(jù)表明實際飛行中的流量偏大，對應(yīng)的推力亦偏大，因此提前耗盡；同時關(guān)機前的飛行速度也一直大于彈道設(shè)計值，與實際飛行情況吻合。

b）針對圖3中60 s前后X向速度偏差較大的狀況，開展氣動阻力系數(shù)偏差辨識，結(jié)果表明氣動阻力修正難以消除該項偏差。

c）100 s后的X向速度偏差持續(xù)增大，則可能為推進劑溫升所致，分析中的火箭采用推進劑加溫后自身增壓方案，存在增壓氣體與推進劑換熱的狀況，缺乏推進劑溫度測量數(shù)據(jù)，難以進一步精細化分析。

d）對彈道計算中存在的“底部阻力”項進行了對比性分析，結(jié)果表明，不考慮“底部阻力”項，辨識后的彈道與飛行彈道更接近，否則100 s后的X向、Y向速度偏差將增加，其中X向速度偏差將增加到辨識結(jié)果的兩倍以上。

e）氣動壓心及法向力系數(shù)亦存在偏差，部分時刻的壓心偏差已經(jīng)接近甚至超過設(shè)計要求考慮的±5%偏差范圍，即存在工程設(shè)計不包絡(luò)（已經(jīng)出現(xiàn)的飛行偏差）風險。

3 結(jié)束語

基于源變量的動力學(xué)集成仿真技術(shù)有機耦合多專業(yè)間的參數(shù)，為運載火箭姿態(tài)動力學(xué)天地一致性分析提供了數(shù)學(xué)基礎(chǔ)，本次包括遺傳算法在內(nèi)的新型搜索算法很好地解決了辨識參數(shù)多的障礙（本次辨識中超過210個參數(shù)），辨識出的“不考慮底部阻力更接近飛行”、“部分時刻氣動壓心偏差已經(jīng)接近甚至超過±5%的設(shè)計考慮偏差邊界”等結(jié)果，將對后續(xù)優(yōu)化設(shè)計提供支撐；同時為火箭飛行測量參數(shù)優(yōu)化設(shè)置提供參考。

在本次任務(wù)分析的基礎(chǔ)上，快速開展了多次不同型號的飛行數(shù)據(jù)分析，顯示出本方法的任務(wù)適應(yīng)能力。

分析中還發(fā)現(xiàn)，待辨識的偏差參數(shù)過多使得分析效率偏低（本次任務(wù)中的210個參數(shù)，大約需要自動執(zhí)行3天時間）；也存在阻力修正系數(shù)、壓心系數(shù)等部分參數(shù)容易陷入局部最優(yōu)點的狀況，后續(xù)將探索針對性更強的辨識尋優(yōu)算法。