月地高速再入返回器彈道重建與氣動力參數(shù)辨識

李 齊，魏昊功，張正峰，陳春亮，彭 兢，楊孟飛

(1. 北京空間飛行器總體設(shè)計部，北京 100094；2. 中國空間技術(shù)研究院，北京 100094)

0 引 言

2020年12月17日，嫦娥五號返回器攜帶月球樣品，以接近11 km/s的速度再入地球大氣層，以半彈道跳躍式再入飛行后，降落在內(nèi)蒙古四子王旗著陸場，圓滿完成了著陸回收任務(wù)。返回器以近第二宇宙速度半彈道跳躍式再入大氣層，在我國航天工程中尚屬首次，與我國近地軌道再入返回器如神舟飛船和返回式衛(wèi)星相比，再入過程中氣動效應(yīng)更為復(fù)雜，對氣動數(shù)據(jù)的精準(zhǔn)度要求更為苛刻[1]。

返回器再入彈道重建與氣動力參數(shù)辨識有著非常重要的意義。通過獲取真實環(huán)境下的飛行數(shù)據(jù)，利用系統(tǒng)辨識理論和技術(shù)，通過測量數(shù)據(jù)辨識獲得真實飛行的氣動力參數(shù)，對于分析高速再入彈道特性、氣動特性與控制規(guī)律，驗證氣動設(shè)計與理論假設(shè)，驗證地面試驗和數(shù)值計算結(jié)果，改進(jìn)氣動設(shè)計方法和工具，為飛行器控制系統(tǒng)提供正確的氣動數(shù)學(xué)模型都具有重要價值[2]。Ivchenko等[3]通過建立全局最小二乘優(yōu)化問題對探空火箭碎片的GPS數(shù)據(jù)進(jìn)行彈道重建。Karlgaard等[4]采用迭代廣義卡爾曼濾波方法分別完成對SFDT-1和SFDT-2飛行試驗器的彈道重建、大氣重建和氣動力重建，所獲得的飛行試驗數(shù)據(jù)用于系列型號研制。

月地高速再入返回器飛行測量數(shù)據(jù)具有以下特點(diǎn)： 1) 跳躍式再入使得飛行高度變化大，由此導(dǎo)致氣動參數(shù)變化范圍大，同時由于跳躍高度是在80 km以上的稀薄流區(qū)，氣動參數(shù)不確定度大； 2) 由于高速再入大氣層時形成的等離子鞘的影響，雷達(dá)測量存在黑障區(qū)，導(dǎo)致外測數(shù)據(jù)不完整。

對于氣動力參數(shù)不確定性問題，王貴東等[5-6]利用克拉馬-羅界，結(jié)合置信度得到測量噪聲產(chǎn)生的數(shù)據(jù)偏差區(qū)間，并采用蒙特卡洛分析方法得到了時序氣動參數(shù)的偏差估計。隨著計算機(jī)技術(shù)不斷提高，將神經(jīng)網(wǎng)絡(luò)等機(jī)器學(xué)習(xí)算法應(yīng)用于氣動參數(shù)辨識成為可能，但采用學(xué)習(xí)算法存在樣本數(shù)量少的問題。王超等[7]采用反向傳播神經(jīng)網(wǎng)絡(luò)建立了高速飛行器全彈道的非線性氣動力模型，并驗證建模方法的可行性。浦甲倫等[8]以神經(jīng)網(wǎng)絡(luò)為基礎(chǔ)，利用支持向量機(jī)擴(kuò)充樣本，設(shè)計氣動參數(shù)智能辨識方法。風(fēng)洞自由飛試驗同樣存在數(shù)據(jù)量少的問題。對此，張?zhí)戽萚9]基于極大似然法對風(fēng)洞自由飛試驗開展氣動參數(shù)辨識，利用有限測量數(shù)據(jù)獲取飛行器的氣動特性。宋威等[10]開展了大鈍頭返回艙跨聲速風(fēng)洞自由飛試驗，考察動態(tài)穩(wěn)定性問題，采用線性和非線性氣動參數(shù)辨識獲得動態(tài)和靜態(tài)穩(wěn)定性參數(shù)。

本文針對月地高速再入返回器跳躍式再入彈道氣動參數(shù)變化范圍大、氣動參數(shù)不確定度大和數(shù)據(jù)不完整等特點(diǎn)，建立了適用于返回器外形的氣動力參數(shù)模型，完成了彈道重建與氣動力參數(shù)辨識，開展了不同大氣模型對氣動參數(shù)辨識的影響分析，完成了飛行辨識數(shù)據(jù)修正。通過與考慮復(fù)雜氣動效應(yīng)的地面預(yù)測氣動數(shù)據(jù)對比，獲得了相互校驗結(jié)論。

1 彈道重建方法

返回器將內(nèi)、外測數(shù)據(jù)作為輸入，通過數(shù)據(jù)預(yù)處理和彈道重建，獲得沿再入彈道的高度、速度、加速度、角速率、姿態(tài)角以及動壓、馬赫數(shù)等彈道與飛行參數(shù)。將這些輸出參數(shù)再作為輸入，聯(lián)合大氣密度模型和返回器的物理幾何參數(shù)，建立運(yùn)動方程組進(jìn)行反解，最終獲得返回器沿彈道真實飛行狀態(tài)下的氣動力參數(shù)和發(fā)動機(jī)開機(jī)時的力矩參數(shù)。彈道重建與氣動力辨識過程的相互關(guān)系如圖1所示。

圖1 彈道重建與氣動力參數(shù)辨識結(jié)構(gòu)框圖Fig.1 Trajectory reconstruction and aerodynamic parameter identification structure

用于彈道重建和氣動力參數(shù)辨識的返回器飛行數(shù)據(jù)分為內(nèi)測和外測數(shù)據(jù)。其中，內(nèi)測數(shù)據(jù)包括來自慣性測量單元的加速度、四元數(shù)以及來自光纖陀螺的姿態(tài)角速率；外測數(shù)據(jù)包括來自地基/船上雷達(dá)的彈道坐標(biāo)、經(jīng)度、緯度、高度以及來自風(fēng)場探空氣球的溫度、氣壓、風(fēng)速、風(fēng)向等。

1.1 數(shù)據(jù)預(yù)處理

由于返回器本身及飛行環(huán)境條件的復(fù)雜性，傳感器和數(shù)據(jù)采集系統(tǒng)非理想性等眾多因素的影響，不可避免地使飛行實測數(shù)據(jù)中含有確定性誤差和隨機(jī)誤差。因此在將飛行測量數(shù)據(jù)用于氣動參數(shù)辨識之前，應(yīng)首先對試驗數(shù)據(jù)進(jìn)行各種預(yù)處理，以盡可能地消除各種確定性誤差和隨機(jī)誤差，提高氣動辨識結(jié)果的準(zhǔn)確度。具體內(nèi)容包括：

1)時間的統(tǒng)一：將外測數(shù)據(jù)與內(nèi)測數(shù)據(jù)時間原點(diǎn)進(jìn)行統(tǒng)一。

2)內(nèi)測四元數(shù)換算姿態(tài)角

利用內(nèi)測四元數(shù)可以得到返回器在再入坐標(biāo)系(或某地面固連坐標(biāo)系)中的姿態(tài)角。

從再入坐標(biāo)系到體軸系的轉(zhuǎn)換矩陣為:

T=T3·T2·T1

于是，返回器相對于再入坐標(biāo)系中的姿態(tài)角可由下式求得：

(1)

式中:γ為滾轉(zhuǎn)角，ψ為偏航角，?為俯仰角,T為轉(zhuǎn)換矩陣。

如果把上述算法中的再入點(diǎn)經(jīng)緯度用當(dāng)?shù)亟?jīng)緯度代替，得到的就是返回器的當(dāng)?shù)刈藨B(tài)角。

3)外測數(shù)據(jù)坐標(biāo)轉(zhuǎn)換及平滑

對雷達(dá)測量的矢徑、仰角、方位角向再入坐標(biāo)系坐標(biāo)的換算，以及地心固連系中的坐標(biāo)向再入坐標(biāo)系坐標(biāo)的換算。

4)返回器坐標(biāo)系統(tǒng)一

通過坐標(biāo)變換，返回器所受地球重力加速度在再入坐標(biāo)系中的投影可表示為：

(2)

式中：gx,gy,gz為重力加速度；grωx,grωy,grωz為重力加速度徑向分量;aex,aey,aez為離心加速度;acx,acy,acz為科氏加速度。

5)實測密度及風(fēng)的修正

根據(jù)再入前落區(qū)附近0～25 km高度范圍內(nèi)的風(fēng)場實測數(shù)據(jù)對辨識所采用的大氣密度模型進(jìn)行修正。由于實測的一般是水平風(fēng)，不只是風(fēng)速的大小，還有風(fēng)的方向，因此需要通過姿態(tài)角把風(fēng)投影到返回器體軸系中才能進(jìn)行風(fēng)的修正。由于再入坐標(biāo)系原點(diǎn)與落區(qū)相差往往超過20°的地心角，因此進(jìn)行風(fēng)的修正時必須考慮地球的曲率，即考慮再入點(diǎn)與落區(qū)的經(jīng)緯度差別。

6)跳點(diǎn)的識別、剔除與補(bǔ)正

飛行中外界干擾或儀器偶然跳動，可能使測量結(jié)果出現(xiàn)跳點(diǎn)。如何準(zhǔn)確區(qū)分跳點(diǎn)與數(shù)據(jù)分散度，需要建立正確的判據(jù)。通常跳點(diǎn)剔除所采用的判據(jù)，如Chauvenet[11]、Grubbs[12]等提出的判據(jù)，無法直接用于氣動參數(shù)辨識所用的角速率、過載這類振蕩型試驗數(shù)據(jù)；基于狀態(tài)方程采用卡爾曼濾波、分割辨識算法或最大似然法也可用于識別和剔除跳點(diǎn)，但對數(shù)據(jù)預(yù)處理來說過于復(fù)雜。本文采用簡單的滑弧算法。

7)數(shù)據(jù)平滑

數(shù)據(jù)預(yù)處理的另一項重要工作就是采用數(shù)據(jù)平滑器將高頻噪聲去除，留下氣動參數(shù)辨識所需的低頻信息。工程上常采用的平滑算法很多，如多項式平滑、樣條函數(shù)平滑、最優(yōu)Fourier平滑等。本文研究采用了多項式平滑方法。

8)時間延遲修正

濾波器時間延遲、采樣時刻差異及傳導(dǎo)時刻差，通常在測量系統(tǒng)設(shè)計時就是已知的，只需用下式進(jìn)行修正：

z(t)=zm(t+Δt)

(3)

式中:z是測量向量,下標(biāo)m表示測量結(jié)果。

傳感器對外作用的動態(tài)響應(yīng)的時間延遲也能夠以一個常值時間延遲Δt來近似。

9)傳感器安裝位置校正

陀螺對返回器姿態(tài)角速率的測量值與它在飛行器上的安裝位置無關(guān)，但加速度的測量值與其安裝位置有關(guān)，傳感器不僅感受質(zhì)心處的值，同時還感受由角速率引起的附加值。因此，在對加速度數(shù)據(jù)進(jìn)行預(yù)處理時，應(yīng)將實測值換算成質(zhì)心處的值。

1.2 彈道重建算法

輸出誤差模型的一般形式為：

(4)

式中:x∈Rn,u∈Rr,y∈Rm分別是狀態(tài)向量、輸入向量、觀測向量；θ∈Rp為待辨識參數(shù)向量；z∈Rm是測量向量；N是采樣數(shù)據(jù)長度；f,g分別是相應(yīng)維的非線性實值函數(shù)；F是時不變過程噪聲矩陣；w(t)與v(k)分別是過程噪聲和觀測噪聲，假定為獨(dú)立無關(guān)的零均值高斯白噪聲，則極大似然估計的指標(biāo)函數(shù)為:

(5)

式中:v(k)為輸出誤差。

(6)

(7)

(8)

式中:

(9)

迭代求解的一般過程：首先給定辨識參數(shù)初值及狀態(tài)初值，根據(jù)式(1)的狀態(tài)方程及觀測方程計算狀態(tài)變量及觀測量，同時計算輸出誤差；然后根據(jù)式(7)計算協(xié)方差矩陣R；通過求解靈敏度方程計算靈敏度矩陣，代入式(8)進(jìn)行辨識參數(shù)更新，反復(fù)迭代直至最終收斂。

式(4)中狀態(tài)方程的具體形式為:

(10)

其中:

(11)

θ=(Δp,ΔV,Δq,Δa,Δω)

(12)

式(4)中第三行測量方程的具體形式為:

(13)

式中:下標(biāo)m表示測量結(jié)果，εi(i=1,2,3)為高斯白噪聲。

2 攻角和側(cè)滑角計算方法

基于重建后的彈道數(shù)據(jù)計算攻角和側(cè)滑角。由于在較低高度下的飛行速度較低，風(fēng)對飛行器攻角和側(cè)滑角影響較大，因此必須考慮風(fēng)速的影響。體軸系速度:

(14)

(15)

式中:α為攻角，β為側(cè)滑角，V為合速度。

3 氣動力參數(shù)辨識方法

3.1 線性回歸法

將系統(tǒng)模型方程以及回歸方程推廣為包含多個變量的情形：

(16)

式中:z=[z(1),z(2), …,z(N)]T為N×1維測量向量；θ=[θ0,θ1,θ2, …,θn]T為np×1維未知參數(shù)向量，np=n+1；X=[1,ξ1, …,ξn]由回歸量構(gòu)成N×np維向量，回歸量ξj，j=1, 2, …,n是由參數(shù)構(gòu)成的函數(shù)；v=[v(1),v(2), …,v(N)]T由測量誤差構(gòu)成的N×1維向量，隨機(jī)誤差項vm包含未知參數(shù)擾動量以及因變量的隨機(jī)測量誤差。

方程誤差辨識模型中，誤差項可假定具有以下性質(zhì)：

E(vvT)=σ2I

(17)

給定z，參數(shù)θ的最佳估計值來自于最小化測量輸出和模型輸出的方差和，即方程誤差辨識的目標(biāo)函數(shù)為：

(18)

極小化J(θ)，得：

(19)

(20)

其中:np×np維矩陣XTX是對稱矩陣。

3.2 氣動力模型

對于氣動參數(shù)辨識，無量綱的氣動力系數(shù)屬于因變量。通過最小化每一個運(yùn)動方程的誤差方差和，可以得到相應(yīng)的氣動力系數(shù)。氣動力系數(shù)值無法通過飛行直接測量，必須通過其它測量值和運(yùn)動方程解算。方程簡化為：

(21)

(22)

回歸量一般由飛行器狀態(tài)變量和控制變量的直接測量值解算得到?；貧w量基于經(jīng)驗給出備選集合，再從備選回歸量集中選擇與輸出量相關(guān)性高的項。氣動力初始模型如下：

(23)

4 辨識結(jié)果

4.1 彈道重建結(jié)果

重建彈道與外測數(shù)據(jù)對比如圖2所示。本次飛行內(nèi)測數(shù)據(jù)與外測數(shù)據(jù)一致性較好，最大偏差出現(xiàn)在第一次再入跳出點(diǎn)附近，最大高度偏差約為1.8 km。經(jīng)過彈道重建修正后，重建彈道數(shù)據(jù)與外測數(shù)據(jù)最大高度偏差均小于1 km。彈道重建結(jié)果見表1。

圖2 重建彈道高度Fig.2 Reconstructed flight altitude

表1 彈道重建結(jié)果Table 1 Trajectory reconstruction results

4.2 攻角與側(cè)滑角辨識結(jié)果

采用著陸區(qū)風(fēng)場數(shù)據(jù)修正前后返回器再入飛行攻角隨再入時間的變化曲線如圖3所示。由圖可見，大部分時間攻角變化平緩，振幅較?。划?dāng)高度小于25 km后，返回器飛行穩(wěn)定性減弱，姿態(tài)振蕩增大。經(jīng)過著陸場風(fēng)修正后，25 km以下攻角變化規(guī)律與修正前相反，表征了返回器亞跨聲速段風(fēng)對姿態(tài)的影響十分顯著。

圖3 飛行攻角辨識值沿彈道變化Fig.3 Reconstructed angle of attack

采用著陸區(qū)風(fēng)場數(shù)據(jù)修正前后返回器側(cè)滑角沿彈道的變化曲線如圖4所示。由圖可見，飛行側(cè)滑角基本在0°附近振蕩。其中，30～80 km側(cè)滑角基本在±2°范圍內(nèi)變化。30 km以下和80 km以上側(cè)滑角振蕩較大，這主要是由于返回器在亞跨聲速范圍動不穩(wěn)定和高空動穩(wěn)定性減弱所致。風(fēng)修正后，側(cè)滑角均值隨高度和速度的降低逐漸增大。

圖4 飛行側(cè)滑角辨識值沿彈道變化Fig.4 Reconstructed sideslip angle

4.3 氣動力參數(shù)辨識結(jié)果

基于辨識算法，首先采用美國1976大氣模型辨識出氣動升力系數(shù)和氣動阻力系數(shù)。圖5為升力系數(shù)、阻力系數(shù)和升阻比沿彈道變化。通過與地面設(shè)計配平值對比可知，辨識得到的升力系數(shù)在高度30 km以上低于地面設(shè)計值，辨識得到的阻力系數(shù)在高度30～90 km之間低于地面設(shè)計值，結(jié)合升力系數(shù)對比結(jié)果，可推測辨識采用的大氣密度值高于實際飛行時同高度大氣密度。由于中高空大氣模型對氣動力參數(shù)(升力系數(shù)和阻力系數(shù))辨識結(jié)果影響很大，為獲得更為準(zhǔn)確的辨識數(shù)據(jù)，需要開展大氣模型分析，盡量減小辨識時的大氣密度輸入誤差，此工作將在下一部分討論。

圖5 升、阻力系數(shù)辨識值與地面設(shè)計值沿彈道分布對比Fig.5 Reconstructed lift, drag coefficients compared with design value

此外，由于跳出大氣段傳感器修正誤差數(shù)值與實際加速度數(shù)值在同一量級，因此跳出大氣段升阻系數(shù)辨識結(jié)果存在較大誤差。因此本次90 km以上的氣動力參數(shù)辨識值無法作為地面預(yù)測數(shù)據(jù)的校驗依據(jù)。

為進(jìn)一步分析大氣模型偏差對氣動力辨識的影響，本文另取兩種大氣模型：中國標(biāo)準(zhǔn)大氣模型[13-14]以及NRLMSISE-00大氣模型[15]，分別開展氣動力參數(shù)辨識。圖6顯示了兩種大氣模型下軸向力系數(shù)辨識結(jié)果與美國1976大氣模型辨識值在各高度下的相對偏差分布。由圖6可見，60 km以上區(qū)域不同大氣模型之間氣動力辨識差異較大，最大可達(dá)到約35%；60 km以下美國海軍大氣模型與國家軍用標(biāo)準(zhǔn)大氣模型辨識所得氣動力系數(shù)基本一致。

圖6 不同大氣模型對應(yīng)軸向力系數(shù)相對偏差沿高度分布Fig.6 Deviation of axial force coefficients with respect to altitude

將辨識采用的大氣密度模型按照三種大氣模型的均值進(jìn)行擬合后，代入返回器氣動參數(shù)辨識，將所得配平攻角與配平升阻比曲線與考慮真實氣體效應(yīng)和燒蝕變形影響的地面預(yù)測曲線進(jìn)行對比，如圖7所示。由圖7可見，除90 km以上的高空和30 km以下的亞跨聲速段外，配平攻角與配平升阻比的辨識值與地面預(yù)測值均較為接近，同時刻同高度下配平攻角差值在1°內(nèi)，配平升阻比偏差在5%以內(nèi)。本次氣動參數(shù)辨識數(shù)據(jù)與地面預(yù)測數(shù)據(jù)完成了相互校驗。

圖7 配平氣動特性地面預(yù)測曲線與辨識修正曲線對比Fig.7 Reconstructed trim angle of attackand lift drag ratio compared to design value

5 結(jié) 論

本文針對月地高速再入返回器跳躍式再入彈道氣動參數(shù)變化范圍大、氣動參數(shù)不確定度大和數(shù)據(jù)不完整等特點(diǎn)，建立了適用于近第二宇宙速度半彈道跳躍式再入的氣動力參數(shù)辨識模型，完成了彈道重建與氣動力參數(shù)辨識。辨識結(jié)果表明，90 km以上高空的大氣模型對辨識影響較大，除此以外，超聲速和高超聲速氣動力辨識結(jié)果與設(shè)計值均較為接近，校驗了相應(yīng)速域的氣動設(shè)計方法。