基于傾側(cè)角反饋控制的預(yù)測校正再入制導(dǎo)方法

趙江，周銳

(北京航空航天大學(xué)自動化科學(xué)與電氣工程學(xué)院，北京100191)

0 引言

高超聲速飛行器再入滑翔制導(dǎo)性能受到了大氣模型參數(shù)的不確定性和動力學(xué)系統(tǒng)的強(qiáng)擾動性等諸多因素的制約，已經(jīng)成為各國航空航天領(lǐng)域研究的重點(diǎn)問題和難點(diǎn)問題［1-3］。再入制導(dǎo)方法通常分為標(biāo)準(zhǔn)軌跡制導(dǎo)和預(yù)測校正制導(dǎo)兩大類。標(biāo)準(zhǔn)軌跡制導(dǎo)方法通過設(shè)計導(dǎo)引律以實(shí)時跟蹤事先確定的標(biāo)準(zhǔn)再入軌跡，對彈載計算機(jī)的性能要求比較低，但是較大的初始軌跡偏差可能會導(dǎo)致制導(dǎo)性能嚴(yán)重退化。預(yù)測校正制導(dǎo)方法則是通過預(yù)測再入飛行器的落點(diǎn)位置偏差來對制導(dǎo)指令進(jìn)行實(shí)時修正，能夠顯著降低初始散布誤差對制導(dǎo)性能的影響，提高制導(dǎo)精度。因此，預(yù)測校正制導(dǎo)成為高超聲速飛行器再入滑翔制導(dǎo)方法的發(fā)展趨勢［4-5］。

文獻(xiàn)［6］利用準(zhǔn)平衡滑翔特性設(shè)計了一種簡單可行的預(yù)測校正制導(dǎo)方法，該方法將再入滑翔飛行中的過程約束轉(zhuǎn)化為傾側(cè)角的邊界約束，通過對運(yùn)動方程進(jìn)行數(shù)值積分來預(yù)測飛行器的待飛航程誤差，實(shí)時校正傾側(cè)角控制指令。在此基礎(chǔ)上，文獻(xiàn)［7 -8］提出了分段的預(yù)測校正制導(dǎo)方法，進(jìn)一步提高了落點(diǎn)預(yù)測精度和迭代計算效率。文獻(xiàn)［9 -10］則以能量為自變量建立三自由度運(yùn)動學(xué)方程，通過求解待飛航程誤差與能量的近似關(guān)系，分別設(shè)計縱向制導(dǎo)律和側(cè)向制導(dǎo)律，增強(qiáng)了預(yù)測校正制導(dǎo)方法的靈活性。文獻(xiàn)［11 -13］分析了大氣模型的不確定性和氣動參數(shù)偏差對落點(diǎn)預(yù)測的影響，在預(yù)測校正制導(dǎo)中引入了參數(shù)估計和修正環(huán)節(jié)，提高了再入制導(dǎo)算法的適應(yīng)性。文獻(xiàn)［14 -15］在標(biāo)準(zhǔn)軌跡制導(dǎo)和預(yù)測校正制導(dǎo)的基礎(chǔ)上提出了混合再入制導(dǎo)方法，通過設(shè)計誤差補(bǔ)償策略和算法切換機(jī)制，將再入制導(dǎo)歸納為最優(yōu)化求解問題，增強(qiáng)了再入制導(dǎo)的抗干擾能力。

為了抑制再入滑翔軌跡的周期性震蕩現(xiàn)象，上述預(yù)測校正制導(dǎo)方法都必須利用準(zhǔn)平衡滑翔條件(QEGC)對制導(dǎo)指令進(jìn)行修正。然而，在飛行航跡角增大或軌道再入速度較低等特殊情況下，QEGC有可能失效［16］。因此，不能完全依賴于該條件進(jìn)行再入制導(dǎo)律設(shè)計。

本文在預(yù)測校正制導(dǎo)方法的基礎(chǔ)上，對升力式高超聲速飛行器的再入制導(dǎo)問題進(jìn)行了研究，提出了一種能夠有效抑制周期性軌跡震蕩的滑翔制導(dǎo)律。

1 再入制導(dǎo)問題

1.1 三自由度運(yùn)動學(xué)方程

高超聲速飛行器的控制問題通常采用六自由度模型來描述再入飛行特性，其結(jié)構(gòu)較為復(fù)雜且主要用于姿態(tài)回路的控制律設(shè)計。本文主要研究再入滑翔制導(dǎo)律設(shè)計，即外環(huán)回路制導(dǎo)律設(shè)計問題。因此，忽略再入飛行過程中風(fēng)速和向心力影響的前提下［17］，建立如下高超聲速再入飛行器的三自由度無量綱運(yùn)動方程［18］

式中:r 為無量綱地心距;V 為飛行器相對地球的無量綱速度;θ 和φ 分別為經(jīng)度和緯度;γ 和ψ 分別為飛行航跡角和航向角;Ω 為地球的自轉(zhuǎn)角速度;σ 為飛行器的傾側(cè)角;L 和D 分別為無量綱的升力加速度和阻力加速度，具體計算公式為

式中:KC= 0.5R0Sref/ m;ρ 為大氣密度;CL和CD分別為升力系數(shù)和阻力系數(shù);R0為地球半徑;Sref為飛行器參考面積;m 為飛行器質(zhì)量。

1.2 再入過程約束

高超滑翔飛行器的再入飛行過程主要考慮熱流密度約束Q、動壓約束q 和過載約束n 的影響，具體計算公式如下所示［19］:

式中:KQ為與飛行器相關(guān)的常值參數(shù);g0為海平面地球引力加速度;Qmax、qmax和nmax分別為熱流密度約束、動壓約束和過載約束的最大幅值，由高超聲速飛行器本身的結(jié)構(gòu)和材料所決定。約束(3)式～(5)式為再入滑翔飛行過程中必須滿足的硬約束條件。

1.3 再入終端約束

不失一般性，引入如下能量形式的e 作為自變量:

本文考慮的再入飛行終端約束主要包括終端高度約束、終端速度約束和終端經(jīng)度與緯度約束，表達(dá)形式為

式中:rf、Vf、θf和φf分別為再入滑翔飛行器地心距、飛行速度、經(jīng)度和緯度的給定終端狀態(tài)。

2 預(yù)測校正制導(dǎo)

2.1 縱向制導(dǎo)律設(shè)計

考慮到再入飛行器初始下降段的熱保護(hù)要求，縱向制導(dǎo)設(shè)計采用升力式滑翔制導(dǎo)常用的二次函數(shù)分段標(biāo)稱攻角剖面，具體計算公式為

式中:α0為再入攻角初值:VT為飛行器攻角幅值開始下降時的臨界速度:K 為常值參數(shù)。

通過對熱流密度約束、動壓約束和過載約束進(jìn)行轉(zhuǎn)化，在高度-速度平面建立再入飛行走廊的下邊界，可以得到飛行器控制傾側(cè)角的幅值約束，表達(dá)形式［20］為

即過程約束下傾側(cè)角的最大幅值為

高超聲速飛行器的再入高度較高，導(dǎo)致初始下降段的氣動力較小，因此該段主要采用開環(huán)制導(dǎo)方式。本文選擇常值傾側(cè)角σ0作為初始下降段的控制變量，其具體計算公式為

式中:ψT為再入初始狀態(tài)的目標(biāo)視線角;Δψ 為航向角誤差。

再入飛行器下降飛行過程中，高度軌跡震蕩的最大幅值通常出現(xiàn)在第一個波谷之后。因此，為了有效地抑制周期性軌跡震蕩，選取飛行器高度變化率首次為0 的狀態(tài)作為初始下降段與滑翔段的交班時刻，即

進(jìn)入滑翔段之后，在每個制導(dǎo)周期內(nèi)，給定初始傾側(cè)角σini的幅值，從當(dāng)前狀態(tài)對運(yùn)動方程進(jìn)行積分(傾側(cè)角符號由側(cè)向制導(dǎo)邏輯決定)，則能夠得到滑翔飛行末端時刻的待飛航程偏差:

式中待飛航程s 選取星下點(diǎn)圓弧長度進(jìn)行近似計算:

則使待飛航程誤差f 為0 的傾側(cè)角幅值可通過割線法快速求解，具體計算公式［6］為

結(jié)合(13)式～(15)式能求得當(dāng)前制導(dǎo)周期內(nèi)傾側(cè)角指令σcmd的幅值(符號由側(cè)向制導(dǎo)邏輯決定)。在不依賴于QEGC 約束的情況下，該制導(dǎo)指令無法導(dǎo)引飛行器完成平穩(wěn)再入飛行。為此，本文在文獻(xiàn)［21］的基礎(chǔ)上，設(shè)計了一種簡易可行的傾側(cè)角反饋控制律以抑制再入滑翔軌跡的周期性震蕩，具體計算公式為

縱向制導(dǎo)律的設(shè)計思想是在傾側(cè)角外環(huán)控制回路添加高度變化率反饋，以抑制再入滑翔過程中產(chǎn)生的高度軌跡震蕩。通常，高超聲速飛行器再入大氣層后會在首次拉起過程中產(chǎn)生最明顯的軌跡震蕩，之后，震蕩幅值隨著速度的下降而逐漸衰弱。因此，在滑翔飛行前期，傾側(cè)角反饋控制律(16)式應(yīng)側(cè)重于抑制再入軌跡震蕩，選取較大的增益;在滑翔飛行后期，反饋控制律(16)式應(yīng)側(cè)重于待飛航程的校正，選取較小的增益，故采用分段遞減函數(shù)表示增益系數(shù):

式中:k1和k2根據(jù)實(shí)際再入初始狀態(tài)及終端約束來設(shè)定。本文中，k1取20 ～30，k2取10 ～15. 通過Monte-Carlo 仿真實(shí)驗及分析可以證明，分段選取增益參數(shù)k 的方法具有可行性。

圖1 參考高度變化率示意圖Fig.1 Example of reference altitude rate

2.2 側(cè)向制導(dǎo)律設(shè)計

側(cè)向制導(dǎo)律的主要任務(wù)是確定傾側(cè)角控制指令的符號。為了實(shí)現(xiàn)再入滑翔飛行的側(cè)向軌跡控制，必須尋找合適的傾側(cè)角反轉(zhuǎn)邏輯。傳統(tǒng)的側(cè)向制導(dǎo)律往往簡單設(shè)計了橫程誤差走廊(或航向角誤差走廊)來實(shí)現(xiàn)傾側(cè)角的反轉(zhuǎn)控制。實(shí)際上，同時考慮橫程誤差和航向角誤差對傾側(cè)角指令符號的影響，才能更為合理和準(zhǔn)確地實(shí)現(xiàn)飛行器的側(cè)向運(yùn)動控制。因此，本文引入歸一化誤差:

式中:Δχ 為橫程誤差;p 為權(quán)值系數(shù)，0≤p≤1. p 的大小決定了歸一化誤差δ 側(cè)重于體現(xiàn)航向角誤差還是橫程誤差。通常情況下p 取0.4 ～0.6 即能滿足側(cè)向制導(dǎo)要求。

以歸一化誤差δ 為基礎(chǔ)設(shè)計側(cè)向誤差走廊可以兼顧橫程誤差和航向角誤差對傾側(cè)角反轉(zhuǎn)指令的影響，其形式如圖2所示。

圖2 側(cè)向誤差走廊設(shè)計Fig.2 Design of lateral error corridor

側(cè)向制導(dǎo)邏輯為:當(dāng)歸一化誤差δ 超出誤差走廊的上邊界時，傾側(cè)角指令符號為負(fù);當(dāng)歸一化誤差δ 超出走廊的下邊界時，傾側(cè)角指令為正;當(dāng)歸一化誤差δ 位于誤差走廊內(nèi)時，傾側(cè)角符號保持不變。該制導(dǎo)邏輯的具體計算公式為

3 仿真分析

本文以升力式高超聲速再入飛行器CAV-H 為仿真對象，其主要特征參數(shù)為:飛行器質(zhì)量m =907.2 kg，氣動參考面積Sref=0.483 9 m2，最大熱流密度約束Qmax=1.0 MW/m2，最大動壓約束qmax=80 kPa，最大過載約束nmax=4.0 g.

飛行器初始再入狀態(tài)如表1所示。終端狀態(tài)約束為:終端目標(biāo)位置(E90°，N30°)，終端高度20 km，終端速度1 800 m/s. 標(biāo)稱攻角剖面參數(shù):K=0.11，α0=22°，VT的Ma=16. 傾側(cè)角指令的取值范圍為-85° ～85°. 傾側(cè)角反饋控制律參數(shù)k1= 25，k2=10.歸一化誤差系數(shù)p =0.5. 仿真計算機(jī)操作系統(tǒng)為Windows XP，主頻3.0 GHz，內(nèi)存2.0 GB. 制導(dǎo)程序采用標(biāo)準(zhǔn)C/C + + 語言編寫，制導(dǎo)周期為1 s.

表1 飛行器初始再入狀態(tài)Tab.1 Initial reentry conditions of vehicle

3.1 標(biāo)準(zhǔn)條件下制導(dǎo)方法仿真分析

為了驗證本文提出制導(dǎo)算法的優(yōu)越性和可靠性，設(shè)計以下兩個算例:算例1 采用傳統(tǒng)的預(yù)測校正再入制導(dǎo)方法;算例2 采用本文改進(jìn)的預(yù)測校正制導(dǎo)算法。圖3(a)～圖3(h)給出了標(biāo)準(zhǔn)條件下(不 考慮再入點(diǎn)散布誤差)的制導(dǎo)算法仿真結(jié)果。

圖3 標(biāo)準(zhǔn)條件下制導(dǎo)仿真結(jié)果Fig.3 Simulation results of reentry guidance in standard conditions

通過對比可得出以下結(jié)論:

1)算例1 的落點(diǎn)經(jīng)度、緯度分別為90.02°和30.02°，精度誤差為3.15 km，算例2 的落點(diǎn)經(jīng)度、緯度分別為89.98°和30.05°，精度誤差為5.91 km，兩種算法的經(jīng)度、緯度誤差均小于0.1°，落點(diǎn)距目標(biāo)距離不超過10 km，滿足了制導(dǎo)方法的精度要求。

2)仿真算例主要針對再入滑翔段制導(dǎo)，即高度80 km 降至20 km，不包括高度20 km 降至地面0 km的末制導(dǎo)問題。因此，與給定終端飛行高度20 km和終端速度1 800 m/s 相比，兩種算法仿真結(jié)果的落點(diǎn)高度誤差小于1.5 km，落點(diǎn)速度誤差小于30 m/s.且與算例1 相比，算例2 飛行高度軌跡未出現(xiàn)周期性震蕩，飛行速度曲線較為平滑，驗證了本文提出的傾側(cè)角反饋控制律具有可行性。

3)熱流密度、動壓和過載約束均未超過給定的最大值限度，與算例1 相比，算例2 中上述典型硬約束的曲線更加平滑且無周期性震蕩，熱/壓防護(hù)效果良好。

4)與算例1 相比，算例2 中傾側(cè)角反饋控制律僅在幅值上做出了修正，沒有增加制導(dǎo)指令反轉(zhuǎn)的次數(shù)，體現(xiàn)了制導(dǎo)算法的可靠性。

5)值得注意的是，算例1 和算例2 的仿真結(jié)果主要用于說明本文提出的傾側(cè)角反饋控制律能夠有效地抑制軌跡震蕩，因而算例1 選取了未添加QEGC 的預(yù)測校正制導(dǎo)算法進(jìn)行對比。從理論上講，添加QEGC 與傾側(cè)角反饋控制率都可以有效地抑制再入軌跡震蕩，但傳統(tǒng)的QEGC 僅僅能夠限制高度震蕩的上邊界，且在航跡角增大或軌道再入速度較低等特殊情況下有可能失效，可見本文提出的算法更具有普遍適用性。

3.2 擾動條件下制導(dǎo)方法仿真分析

再入滑翔飛行過程中，大氣模型參數(shù)的不確定性和動力學(xué)系統(tǒng)的強(qiáng)擾動性都可能造成高超聲速飛行器偏離預(yù)定軌跡和目標(biāo)，從而導(dǎo)致任務(wù)失敗。為了驗證本文提出的制導(dǎo)算法在各種參數(shù)偏差條件下都具有魯棒性，針對不同初始擾動情況下的再入滑翔制導(dǎo)進(jìn)行了Monte Carlo 仿真實(shí)驗。再入點(diǎn)初始散布偏差項、分布類型及偏差限度見表2.

表2 再入點(diǎn)散布誤差Tab.2 Dispersions at the reentry interface

圖4(a)～圖4(i)給出了擾動條件下制導(dǎo)算法的仿真曲線和落點(diǎn)散布情況統(tǒng)計。通過分析可得出以下結(jié)論:

1)如圖4(a)所示，從星下點(diǎn)軌跡可看出，在誤差和擾動存在的情況下，仿真飛行過程沒有偏離原始的再入滑翔軌跡，整體上能滿足預(yù)定的飛行航程要求。

2)如圖4(b)所示，在Monte Carlo 數(shù)值仿真算例中，高度曲線未出現(xiàn)周期性的振蕩現(xiàn)象，再入飛行軌跡比較平穩(wěn)，且終端高度和速度約束能夠滿足制導(dǎo)要求，熱流、動壓和過載也未超過給定的最大約束限制，如圖4(f)～圖4(h)所示，驗證了本文提出算法的魯棒性。

3)根據(jù)落點(diǎn)散布誤差可得，如圖4(i)所示，再入落點(diǎn)距給定終端目標(biāo)的距離不超過15 km，且經(jīng)度、緯度誤差均小于0.15°，滿足了預(yù)測校正制導(dǎo)的精度要求。

4)從傾側(cè)角曲線的變化趨勢可知，在再入滑翔飛行的前半段，控制指令主要受到高度反饋?zhàn)兓实挠绊懀鴼w一化誤差走廊更多地作用于再入滑翔飛行的后半段。

5)對比不同擾動條件對預(yù)測校正制導(dǎo)的影響程度能夠發(fā)現(xiàn)，速度偏差和航跡角偏差對制導(dǎo)精度的影響相對較大，而高度偏差對制導(dǎo)精度的影響相對較小。

4 結(jié)論

本文針對大升阻比高超聲速飛行器再入制導(dǎo)問題，提出一種改進(jìn)的傾側(cè)角反饋控制律和側(cè)向誤差走廊，研究分析和仿真結(jié)果表明:

1)考慮軌跡震蕩抑制的傾側(cè)角反饋控制、不改變傳統(tǒng)預(yù)測校正制導(dǎo)律的表達(dá)形式和傾側(cè)角反轉(zhuǎn)的次數(shù)，易于實(shí)現(xiàn)。

2)歸一化的側(cè)向誤差走廊設(shè)計同時兼顧了橫程誤差和航向角誤差對傾側(cè)角反轉(zhuǎn)邏輯的影響，使側(cè)向運(yùn)動控制更加合理和有效。

3)標(biāo)準(zhǔn)條件下的仿真結(jié)果表明，與傳統(tǒng)的預(yù)測校正制導(dǎo)相比，本文提出的制導(dǎo)方法能夠有效地抑制再入軌跡的周期性震蕩。

4)擾動條件下的Monte Carlo 仿真結(jié)果表明，本文提出的制導(dǎo)方法具有良好的魯棒性。

圖4 擾動條件下制導(dǎo)仿真結(jié)果Fig.4 Simulation results of reentry guidance with random disturbance

References)

［1］Sarah N D，Nesrin S K. Survey of planetary entry guidance algorithms［J］. Progress in Aerospace Sciences，2014，68(1):22 -28.

［2］Zhao J，Zhou R. Reentry trajectory optimization for hypersonic vehicle satisfying complex constraints［J］. Chinese Journal of Aeronautics，2013，26(6):1544 -1553.

［3］袁宴波，張科，薛曉東. 基于Radau 偽譜法的制導(dǎo)炸彈最優(yōu)滑翔彈道研究［J］. 兵工學(xué)報，2014，35(8):1179 -1186.YUAN Yan-bo，ZHANG Ke，XUE Xiao-dong. Optimization of glide trajectory of guided bombs using a Radau pseudo-spectral method［J］.Acta Armamentarii，2014，35(8):1179 -1186. (in Chinese)

［4］Zhao J，Zhou R，Jin X. Progress in reentry trajectory planning for hypersonic vehicle［J］. Journal of Systems Engineering and Electronics，2014，25(4):627 -639.

［5］趙江，周銳. 考慮禁飛區(qū)規(guī)避的預(yù)測校正再入制導(dǎo)方法［J］.北京航空航天大學(xué)學(xué)報，2015，41(5):864 -870.ZHAO Jiang，ZHOU Rui. Predictor-corrector reentry guidance with no-fly zone constraints［J］. Journal of Beijing University of Aeronautics and Astronautics，2015，41(5):864 -870. (in Chinese)

［6］Xue S，Lu P. Constrained predictor-corrector entry guidance［J］.Journal of Guidance Control and Dynamics，2010，33(4):1273 -1281.

［7］王智，唐碩，閆曉東. 高超聲速滑翔飛行器約束預(yù)測校正再入制導(dǎo)［J］. 飛行力學(xué)，2012，30(2):175 -180.WANG Zhi，TANG Shuo，YAN Xiao-dong. Constrained predictorcorrector reentry guidance for hypersonic glide vehicle［J］. Flight Dynamics，2012，30(2):175 -180. (in Chinese)

［8］劉冠南，周浩，陳萬春. 高超聲速飛行器再入多段導(dǎo)引方法研究［J］. 飛行力學(xué)，2012，30(4):337 -341.LIU Guan-nan，ZHOU Hao，CHEN Wan-chun. A study on the multiple phrase reentry guidance method of hypersonic vehicle［J］.Flight Dynamics，2012，30(4):337 -341. (in Chinese)

［9］王青，冉茂鵬，趙洋. 基于預(yù)測校正法的高超聲速飛行器再入制導(dǎo)［J］. 北京航空航天大學(xué)學(xué)報，2013，39(12):1564 -1568.WANG Qing，RAN Mao-peng，ZHAO Yang. Reentry guidance for hypersonic vehicle based on predictor-corrector method［J］. Journal of Beijing University of Aeronautics and Astronautics，2013，39(12):1564 -1568. (in Chinese)

［10］王青，莫華東，吳振東，等. 基于能量的高超聲速飛行器再入混合制導(dǎo)方法［J］. 北京航空航天大學(xué)學(xué)報，2014，40(5):580 -585.WANG Qing，MO Hua-dong，WU Zheng-dong，et al. Energybased hybrid reentry guidance for hypersonic vehicles［J］. Journal of Beijing University of Aeronautics and Astronautics，2014，40(5):580 -585. (in Chinese)

［11］曾憲法，王潔瑤，王小虎. 基于能量和解析預(yù)測校正的滑翔制導(dǎo)［J］. 系統(tǒng)工程與電子技術(shù)，2013，35(12):2583 -2589.ZENG Xian-fa，WANG Jie-yao，WANG Xiao-hu. Gliding guidance based on energy and analytical predictor-corrector［J］. Systems Engineering and Electronics，2013，35(12):2583 -2589.(in Chinese)

［12］李強(qiáng)，夏群利，何鏡，等. 基于大氣預(yù)估的再入飛行機(jī)動減速制導(dǎo)方法［J］. 兵工學(xué)報，2013，34(9):1091 -1096.LI Qiang，XIA Qun-li，HE Jing，et al. Maneuvering-deceleration guidance algorithm based on atmosphere estimation for reentry vehicle［J］. Acta Armamentarii，2013，34(9):1091 -1096. (in Chinese)

［13］梁子璇，任章. 基于在線氣動參數(shù)修正的預(yù)測制導(dǎo)方法［J］.北京航空航天大學(xué)學(xué)報，2013，39(7):853 -857.LIANG Zi-xuan，REN Zhang. Predictive reentry guidance with aerodynamic parameter online correction［J］. Journal of Beijing University of Aeronautics and Astronautics，2013，39(7):853 -857. (in Chinese)

［14］Wang J，Qu X，Ren Z. Hybrid reentry guidance based on the online trajectory planning［J］. Journal of Astronautics，2012，33(9):1217 -1224.

［15］王俊波，田源，任章. 基于最優(yōu)化問題的混合再入制導(dǎo)方法［J］.北京航空航天大學(xué)學(xué)報，2010，36(6):736 -740.WANG Jun-bo，TIAN Yuan，REN Zhang. Mixed guidance method for reentry vehicles based on optimization［J］. Journal of Beijing University of Aeronautics and Astronautics，2010，36(6):736 -740. (in Chinese)

［16］Shen Z，Lu P. Onboard generation of three dimensional constrained entry trajectories［J］. Journal of Guidance，Control，and Dynamics，2003，26(1):111 -121.

［17］Leavitt J A，Mease K D. Feasible trajectory generation for atmospheric entry guidance［J］. Journal of Guidance，Control，and Dynamics，2007，30(2):473 -481.

［18］Zhao J，Zhou R，Jin X. Reentry trajectory optimization based on a multistage pseudospectral method［J］. Scientific World Journal，2014(878193):1 -13.

［19］Zhao J，Zhou R，Jin X. Gauss pseudospectral method applied to multi-objective spacecraft trajectory optimization［J］. Journal of Computational and Theoretical Nanoscience，2014，11 (10):2242 -2246.

［20］Lu P，Xue S. Rapid generation of accurate entry landing footprints［J］. Journal of Guidance Control and Dynamics，2010，33(3):756 -767.

［21］Lu P. Entry guidance:a unified method［J］. Journal of Guidance，Control，and Dynamics，2014，37(3):713 -728.