吳永峰，葉 鳴，姚浩祺，韋才敏

（汕頭大學(xué)數(shù)學(xué)系，廣東 汕頭 515063）

嫦娥三號(hào)軟著陸預(yù)測(cè)及優(yōu)化模型

吳永峰，葉 鳴，姚浩祺，韋才敏

（汕頭大學(xué)數(shù)學(xué)系，廣東 汕頭 515063）

本文是對(duì)2014年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A題的解答.本文利用動(dòng)力學(xué)模型，基于變推力發(fā)動(dòng)機(jī)結(jié)合螺旋搜索等算法解決了嫦娥三號(hào)軟著陸預(yù)測(cè)及優(yōu)化問(wèn)題.在問(wèn)題一中，首先建立物理學(xué)模型，解出嫦娥三號(hào)著陸準(zhǔn)備軌道近、遠(yuǎn)月點(diǎn)速度大?。蝗缓笤谥鳒p速階段，建立了主減速變推力動(dòng)力學(xué)模型，對(duì)主減速階段嫦娥三號(hào)所經(jīng)過(guò)的極角θ進(jìn)行求解，加以考慮月球自轉(zhuǎn)帶來(lái)的影響，確定著陸準(zhǔn)備軌道近月點(diǎn)和遠(yuǎn)月點(diǎn)的位置以及嫦娥三號(hào)在近、遠(yuǎn)月點(diǎn)的速度方向.在問(wèn)題二中，建立變推力發(fā)動(dòng)機(jī)燃耗優(yōu)化模型對(duì)于問(wèn)題一中的主減速階段進(jìn)行優(yōu)化；以安全半徑與速度增量為指標(biāo)建立綜合評(píng)價(jià)模型，求出粗避障式安全的降落范圍；精避障階段以保證安全為前提，建立評(píng)價(jià)指標(biāo)為不平坦程度與平均坡度的綜合評(píng)價(jià)模型確定精確降落區(qū)域.最后根據(jù)對(duì)每個(gè)階段建立的力學(xué)模型，進(jìn)行求解.最后綜合起來(lái)確定嫦娥三號(hào)的著陸軌道和在6個(gè)階段的最優(yōu)控制策略.在問(wèn)題三中，建立誤差模型對(duì)設(shè)計(jì)的著陸軌道和控制策略做相應(yīng)的誤差分析，并且利用單因素敏感性分析對(duì)主減速階段與剩余4階段進(jìn)行自變量與因變量敏感度計(jì)算.

主減速動(dòng)力學(xué)模型；變推力；螺旋搜索算法

0 引 言

現(xiàn)代科學(xué)技術(shù)中，航空航天技術(shù)仍扮演著不可或缺的重要角色，月球作為離地球最近的天體，自古以來(lái)一直吸引著人們.要是想對(duì)月球進(jìn)行近距離乃至實(shí)地勘測(cè)，著陸器及相關(guān)軟著陸軌道設(shè)計(jì)是不可缺少的.之前成功實(shí)現(xiàn)月球軟著陸的國(guó)家也僅有美國(guó)及前蘇聯(lián)，其他國(guó)家在此方面的技術(shù)仍不算成熟.

此次由我國(guó)研制的首次裝有變推力發(fā)動(dòng)機(jī)的嫦娥三號(hào)也于2013年12月2日1時(shí)30分成功發(fā)射，12月6日抵達(dá)月球軌道.然而參考國(guó)內(nèi)外相關(guān)文獻(xiàn)，對(duì)于恒推力發(fā)動(dòng)機(jī)的月球登陸器的軌道設(shè)計(jì)的文獻(xiàn)仍占大多數(shù)，而對(duì)于變推力發(fā)動(dòng)機(jī)的相關(guān)文獻(xiàn)較少.

本文依據(jù)2014年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A題所給出的嫦娥三號(hào)的參數(shù)及所設(shè)計(jì)著陸月球前的著陸軌道的相應(yīng)基本要求，根據(jù)文獻(xiàn)[1]和[2]提出一種針對(duì)裝有變推力發(fā)動(dòng)機(jī)的月球登陸器如何設(shè)計(jì)軟著陸登陸軌道的算法，建立數(shù)學(xué)模型嘗試解決以下問(wèn)題：

（1）確定著陸準(zhǔn)備軌道近月點(diǎn)和遠(yuǎn)月點(diǎn)的位置，以及嫦娥三號(hào)相應(yīng)速度的大小與方向.

（2）確定嫦娥三號(hào)的著陸軌道和在6個(gè)階段的最優(yōu)控制策略.

（3）對(duì)于設(shè)計(jì)的著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析.

注：題目和相關(guān)附件可到全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽官方網(wǎng)站http://www.mcm.edu.cn下載

1 基本假設(shè)

1）主減速階段，嫦娥三號(hào)按著預(yù)定著陸點(diǎn)方向減速.

2）所用的數(shù)據(jù)均能客觀有效地反映嫦娥三號(hào)的狀態(tài).

2 模型的建立與求解

2.1 問(wèn)題一的分析

問(wèn)題一要確定著陸準(zhǔn)備軌道近月點(diǎn)和遠(yuǎn)月點(diǎn)的位置，以及嫦娥三號(hào)相應(yīng)的速度與方向，主要由以下步驟組成:

步驟一：計(jì)算嫦娥三號(hào)在軌道近、遠(yuǎn)月點(diǎn)的速度；

步驟二：首先確定著陸準(zhǔn)備軌道的位置，建立主減速動(dòng)力學(xué)模型，通過(guò)該模型，判斷近、遠(yuǎn)月點(diǎn)的位置，以月球經(jīng)緯度進(jìn)行表示.

2.1.1 計(jì)算嫦娥三號(hào)在軌道近、遠(yuǎn)月點(diǎn)的速度大小[3]

由于嫦娥三號(hào)在軌道運(yùn)動(dòng)中只受到月球萬(wàn)有引力作用，所以此時(shí)遵循機(jī)械能守恒定律，并且根據(jù)開(kāi)普勒第二定律的角動(dòng)量守恒定律可得出以下方程：

已知嫦娥三號(hào)的著陸準(zhǔn)備軌道為橢圓軌道，設(shè)為：

所以，ν1=1.692 2×103m/s，ν2=1.613 9×103m/s，即近日點(diǎn)速度為1.692 2×103m/s，遠(yuǎn)日點(diǎn)速度為1.613 9×103m/s.

2.1.2 確定近、遠(yuǎn)月點(diǎn)的位置

（1）著陸準(zhǔn)備軌道的確定

根據(jù)相關(guān)文獻(xiàn)[4]，了解到月球只有一面始終對(duì)準(zhǔn)地球，為保證探測(cè)器在月面著陸后可以連續(xù)與地面通信，著陸區(qū)應(yīng)該選擇在月球正面，且是由近月點(diǎn)位置下降到達(dá)著陸區(qū)，此時(shí)可以確定，嫦娥三號(hào)著陸準(zhǔn)備軌道的近月點(diǎn)一定位于月球與地球之間，其著陸準(zhǔn)備軌道的大致圖像為圖2.且嫦娥三號(hào)的著陸準(zhǔn)備軌道為與月球赤道傾角為90度的極軌道，可以經(jīng)過(guò)預(yù)定著陸點(diǎn)為（19.51W，44.12N）.說(shuō)明嫦娥三號(hào)著陸準(zhǔn)備軌道的大概位置是：位于月球上一條緯度位置為19.51W的經(jīng)線上方.環(huán)月方向?yàn)椋何靼肭蜃阅舷虮?，東半球自北向南.

圖1 著陸準(zhǔn)備軌道坐標(biāo)圖

圖2 著陸準(zhǔn)備軌道與地、月之間的位置

（2）建立主減速動(dòng)力學(xué)模型

Step1.建立力學(xué)方程

因?yàn)椋憸?zhǔn)備軌道近月點(diǎn)的高度遠(yuǎn)小于月球平均半徑，此時(shí)可以看做近月點(diǎn)均位于月球表面，由萬(wàn)有引力公式可得月球?qū)τ阪隙鹑?hào)位于近月點(diǎn)的引力：

同理，由萬(wàn)有引力公式可以得出，當(dāng)嫦娥三號(hào)位于著陸準(zhǔn)備軌道時(shí)，地球與其之間的引力最大為8N，遠(yuǎn)小于此時(shí)月球?qū)ζ涞囊Γ虼?，可以無(wú)需在模型中考慮地球重力對(duì)于嫦娥三號(hào)的影響.

根據(jù)假設(shè)，主減速階段嫦娥三號(hào)在向著預(yù)定著陸區(qū)域進(jìn)行減速，從而到達(dá)預(yù)定著陸點(diǎn)上空，此時(shí)嫦娥三號(hào)運(yùn)動(dòng)軌跡如圖3所示.

對(duì)于本文討論的變推力重力轉(zhuǎn)彎軟著陸的情況，有α-γ＝180°，且要求推力加速度a=F/m保持不變，于是可得到如下動(dòng)力方程：

其中，α，γ分別為推力矢量F和下降速度矢量ν與當(dāng)?shù)厮椒较虻膴A角，從νθ軸開(kāi)始演逆時(shí)針?lè)较蛐D(zhuǎn)為正.由圖2可知γ≤0，h0，hm分別重力轉(zhuǎn)彎軟著陸階段的嫦娥三號(hào)距離月面的初始和終端高度.

同時(shí)，在速度方向上有：

根據(jù)軟著陸下降的飛行路徑s和矢徑r之間的關(guān)系，如圖4所示，可以得出：

其中，V代表速度矢量的模.

圖3 嫦娥三號(hào)軟著陸下降軌道示意圖

圖4 飛行路徑與矢徑的關(guān)系

Step2.求解反推加速度a

在飛行過(guò)程中，嫦娥三號(hào)的瞬時(shí)機(jī)械能為：

結(jié)合（2）式，可以得到單位時(shí)間內(nèi)機(jī)械能的變化量

且由著陸器路程sf的表達(dá)式：

并根據(jù)重力轉(zhuǎn)彎段的終端條件：s=sf，r=rf=rm+h1，可以得到如下方程：

此時(shí)解出方程便可以得到重力拐彎階段所需要的推力加速度a.

Step3. 求解下降飛行的極角θ

下降飛行的極角θ是指嫦娥三號(hào)從初始下降點(diǎn)（近月點(diǎn)）與月心的連線開(kāi)始到著陸點(diǎn)上空與月心的連線所劃過(guò)的角度（見(jiàn)圖3），求得這一角度，不僅可以利用公式得到下降飛行時(shí)的月面距離，還可根據(jù)著陸點(diǎn)的經(jīng)緯度，對(duì)嫦娥三號(hào)初始下降點(diǎn)經(jīng)緯度進(jìn)行推算.

由圖3可知，嫦娥三號(hào)的瞬時(shí)動(dòng)量矩H=rVcos γ＝rνθ，結(jié)合（1）式可得到單位時(shí)間的動(dòng)量矩變化量：

考慮初始條件s=0，H=H0：同時(shí)帶入重力拐彎段下降的終端條件：s=sf，θ=θf(wàn)，即可得到下降飛行的極角θ：

通過(guò)式（9）就可以求解出嫦娥三號(hào)下降飛行所劃過(guò)的極角θ.

（3）主減速動(dòng)力學(xué)模型求解

考慮到預(yù)定著陸點(diǎn)的海拔為-2 641 m，所以模型中的h1=3 000-2 641=359 m.我們令γ0∈[-89°，-1°]，以1°為步長(zhǎng)，搜索關(guān)于極角θ的解，約束條件有：

利用MATLAB軟件對(duì)該主減速動(dòng)力模型求解，求出滿足解的結(jié)果如下表1：

表1 主減速動(dòng)力學(xué)模型滿足解結(jié)果

表格中的結(jié)果均滿足式（7）有解的情況，從結(jié)果可以得出，所消耗的燃料相差無(wú)幾，可以近似地看成所有的策略均滿足最低能耗.從實(shí)際考慮，嫦娥三號(hào)在主加速階段所用時(shí)間應(yīng)較短，則此時(shí)γ0=-6°時(shí)，時(shí)間為567.959 7 s，此時(shí)所走路程4.805 5×105m，同樣也為全局最小.則可以確定γ0=-6°為主減速動(dòng)力模型最優(yōu)解.此時(shí)，變推力加速度a=3.008 6 m/s2，極角θ＝0.272 4 rad.

（4）嫦娥三號(hào)著陸準(zhǔn)備軌道近、遠(yuǎn)月點(diǎn)位置確定

由已知得，預(yù)定著陸點(diǎn)為（19.51W，44.12N），嫦娥三號(hào)著陸準(zhǔn)備軌道位于月球上一條緯度位置為19.51W的經(jīng)線上方.主減速階段時(shí)極角θ＝0.2724·180/π≈15.61°，即此時(shí)近月點(diǎn)（減速下降點(diǎn)）位置為月球上坐標(biāo)（19.51W，59.73N）上空15 km處或者處于坐標(biāo)（19.51W，28.51N）上空15 km處.

由于嫦娥三號(hào)在軌道上環(huán)月方向?yàn)槲靼肭蜃阅舷虮?，所以?9.51W，61.60N）坐標(biāo)不符合實(shí)際情況，為錯(cuò)解.另外由于月球進(jìn)行自轉(zhuǎn)，自轉(zhuǎn)周期為27.32天，自轉(zhuǎn)方向自西向東，所以在主減速階段621.194 5s，月球已經(jīng)轉(zhuǎn)過(guò)經(jīng)度0.09°.

綜上所述，近月點(diǎn)（動(dòng)力下降點(diǎn)）位置處于月球經(jīng)緯度坐標(biāo)（19.42W，28.51N）上方15 km處，近月點(diǎn)速度方向垂直于嫦娥三號(hào)所受重力方向.遠(yuǎn)月點(diǎn)位置處于月球坐標(biāo)（160.58E，28.51S）上方100 km，遠(yuǎn)月點(diǎn)速度方向垂直于嫦娥三號(hào)所受重力方向.

2.2 問(wèn)題二的模型建立與求解

問(wèn)題二要確定嫦娥三號(hào)的著陸軌道和在6個(gè)階段的最優(yōu)控制策略，主要由以下步驟組成：

步驟一：首先對(duì)6個(gè)階段分別建立模型，求解出各自的最優(yōu)控制策略；

步驟二：將6個(gè)最優(yōu)控制策略進(jìn)行合理性結(jié)合，目的是得出嫦娥三號(hào)的著陸軌道.

2.2.1 對(duì)6個(gè)階段最優(yōu)控制策略模型的建立與求解

（1）著陸準(zhǔn)備軌道

如圖5所示，著陸準(zhǔn)備軌道下的近月點(diǎn)坐標(biāo)（19.51W，28.51N），從近月點(diǎn)開(kāi)始動(dòng)力下降時(shí)，速度與水平方向的夾角大小為6°，遠(yuǎn)月點(diǎn)坐標(biāo)（160.58E，28.51S）.

圖5 著陸準(zhǔn)備軌道示意圖

（2）主減速階段

嫦娥三號(hào)在月球表面實(shí)行軟著陸，由于月球表面不具有大氣，無(wú)法通過(guò)大氣阻力減速，須用發(fā)動(dòng)機(jī)進(jìn)行制動(dòng)，因此燃耗最少就成為了一個(gè)重要的優(yōu)化指標(biāo)，由于主減速階段的耗燃占整個(gè)軟著陸過(guò)程的大部分，基于第一問(wèn)的變推力發(fā)動(dòng)機(jī)的主減速動(dòng)力模型，提出了一種在不考慮月球自轉(zhuǎn)影響下，用變推力發(fā)動(dòng)機(jī)進(jìn)行減速下降的優(yōu)化方法[5].

1.運(yùn)動(dòng)方程

在月心慣性坐標(biāo)系中，采用變推力發(fā)動(dòng)機(jī)，使得嫦娥三號(hào)加速度隨時(shí)間而呈線性變化.其加速度可以表示為：

其中ax，ay為主減速時(shí)加速度，a0，a1，b0，b1為常系數(shù)，t為時(shí)間，主減速軌道任一時(shí)刻的位置和速度分別為：

其中，（x0，y0），（νx0，νy0）分別為初始位置的坐標(biāo)的速度.

式（11）中有方程4個(gè)和變量5個(gè)，控制其中的某個(gè)參數(shù)就可以實(shí)現(xiàn)主減速的軌道優(yōu)化.如若（x0，y0），（νx0，νy0）分別為著陸點(diǎn)的要求坐標(biāo)和主減速結(jié)束時(shí)要求的速度，t1為主減速段總飛行時(shí)間，則加速度系數(shù)為：

控制參數(shù)多，利于實(shí)現(xiàn)燃料的絕對(duì)最省，但不利于計(jì)算的實(shí)現(xiàn).因此，本文主要研究式（10）加速度線性變化.

2.推力加速度和燃耗

因?yàn)槭剑?0）中加速度為推力加速度與引力加速度之和，由此可得嫦娥三號(hào)發(fā)動(dòng)機(jī)需要施加的加速度：

其中：gx，gy為當(dāng)?shù)刂亓铀俣?，隨位置而變，且apx，apy為推力加速度，是時(shí)間的函數(shù)，且：

若嫦娥三號(hào)質(zhì)量為m，發(fā)動(dòng)機(jī)比沖為νe，則需要的發(fā)動(dòng)機(jī)推力：

F的方向可由計(jì)算的推力加速度求得，其對(duì)應(yīng)的質(zhì)量秒流量為：

因?yàn)棣蚭為常數(shù)2 940，結(jié)合前式從而可得出：

由上述推導(dǎo)可知，當(dāng)著陸位置給定時(shí)，已知總飛行時(shí)間，即可求得需要的速度增量，進(jìn)而求得燃料消耗質(zhì)量，此時(shí)我們需要速度增量為t1的函數(shù)；當(dāng)著陸點(diǎn)位置未知時(shí)，需要速度增量為t1和著陸點(diǎn)位置的函數(shù).

3.燃耗優(yōu)化

對(duì)圖6所示的主減速軌跡，設(shè)坐標(biāo)系基準(zhǔn)平面位于初始運(yùn)動(dòng)平面內(nèi)，主減速制動(dòng)時(shí)，因?yàn)橹扅c(diǎn)位于初始運(yùn)動(dòng)平面內(nèi)而無(wú)需改變運(yùn)動(dòng)平面而最省能量.因此，對(duì)不定點(diǎn)著陸優(yōu)化，著陸點(diǎn)位置用以參數(shù)表示即可，即：

其中θ為初始點(diǎn)到著陸點(diǎn)的月心角（即主減速階段的行程角）.

此時(shí)給定θ，t，帶入式（10）－（13）、（17）即可求得最小速度增量Δν.

4.實(shí)際問(wèn)題求解

利用MATLAB軟件，對(duì)第一問(wèn)主減速階段所建立的模型進(jìn)行優(yōu)化（優(yōu)化后燃料消耗為1 053 kg，與第一問(wèn)求解最優(yōu)結(jié)果燃料最小消耗1 058 kg，優(yōu)化效果極小，幾乎可以忽略不計(jì)，所以在主減速階段，我們依然采用問(wèn)題一的求解結(jié)果：近月點(diǎn)（動(dòng)力下降點(diǎn)）位置處于月球經(jīng)緯度坐標(biāo)（19.42W，28.51N）上方15 km處，近月點(diǎn)速度方向垂直于嫦娥三號(hào)所受重力方向.遠(yuǎn)月點(diǎn)位置處于月球坐標(biāo)（160.58E，28.51S）上方100 km，遠(yuǎn)月點(diǎn)速度方向垂直于嫦娥三號(hào)所受重力方向.

（3）快速調(diào)整階段

對(duì)此時(shí)的嫦娥三號(hào)進(jìn)行受力分析，如圖7所示.可以建立力學(xué)方程如下：

其中含有自變量a的關(guān)系式a＝a（φ，θ），輸入φ，θ，即可求出a的值，得出時(shí)間t= νcosφ/acosθ，以及階段結(jié)束后最終速度νf=νsinφ+（g-asinθ），Δν＝at，最后根據(jù)式（16）可得所耗燃料質(zhì)量.

圖6 主減速軌跡圖

圖7 快速調(diào)整階段嫦娥三號(hào)受力分析

（4）粗避障階段

粗避障主要目的是在較大著陸范圍（2 300 m×2 300 m）內(nèi)避開(kāi)明顯危及著陸安全的危險(xiǎn)障礙，為精避障提供較好的安全著陸點(diǎn)選擇區(qū)域，避免出現(xiàn)近距離精避障避無(wú)可避的風(fēng)險(xiǎn)，整體上提高系統(tǒng)安全著陸概率.此外，還要綜合評(píng)估燃料消耗來(lái)選擇最優(yōu)位置.具體算法流程如下：

Step1.圖像分析

首先將距月面2 400 m處的數(shù)字高程圖，2 300 m×2 300 m的范圍分割成每個(gè)碎片具有范圍（100 m×100 m），共有行×列為23×23個(gè)數(shù)的碎片.分別對(duì)每一個(gè)碎片圖像進(jìn)行分析.利用MATLAB讀取每個(gè)碎片圖像，將其每一個(gè)像素點(diǎn)上的高度轉(zhuǎn)換成呈現(xiàn)所有像素點(diǎn)的高度值矩陣.我們找出所有碎片的高度眾數(shù)92，設(shè)立閥值為5，也就是說(shuō)，高度值在[87，97]范圍，我們認(rèn)定其為可以實(shí)行降落的安全高度.高度不在其范圍的像素點(diǎn)即為危險(xiǎn)像素點(diǎn).

Step2. 采用算法確定每一個(gè)碎片的安全半徑.

基于螺旋搜索算法，我們提出了新的改進(jìn)算法：以每一個(gè)碎片的中心為起點(diǎn)，以一個(gè)單元格為長(zhǎng)度，正方形為范圍，逐漸擴(kuò)大正方形的范圍，直至碰上危險(xiǎn)區(qū)域，停止繼續(xù)擴(kuò)大搜索范圍，以之前的范圍確定安全半徑.算法示意圖如圖8：其中，R為安全半徑.

Step3.計(jì)算到達(dá)安全著陸點(diǎn)的速度增量

在粗避障階段，從目標(biāo)上方2 400 m到100 m時(shí)，初始具有垂直速度，水平速度為零.僅考慮水平方向時(shí)，嫦娥三號(hào)到達(dá)安全著陸點(diǎn)時(shí)經(jīng)歷了兩個(gè)階段，在總路程的一半前，做加速運(yùn)動(dòng)，當(dāng)嫦娥三號(hào)運(yùn)動(dòng)路程達(dá)到總路程的一半時(shí)，水平速度達(dá)到最大，此時(shí)改變推力方向，進(jìn)入減速階段，保證在到達(dá)目標(biāo)上方100 m時(shí)，水平速度為零.示意圖如圖9.

圖8 基于螺旋算法的安全半徑搜索算法

圖9 嫦娥三號(hào)粗避障運(yùn)動(dòng)狀態(tài)圖圖

假設(shè)加速階段與減速階段，推力大小相同，推力方向與水平面方向的夾角相同，加減速階段推力方向如圖10所示，即其中F1＝F2，β1＝β2.

因此，我們可以建立力學(xué)方程：

其中x表示嫦娥三號(hào)到安全著陸點(diǎn)的水平平移距離.

此時(shí)ν0為快速調(diào)整姿態(tài)下的最終速度，求出a后，則此時(shí)速度增量可知，從而確定粗避障所耗時(shí)間T=x/Δν，利用式（16）也可確定消耗燃料質(zhì)量.

Step4.建立綜合評(píng)價(jià)模型并求解，根據(jù)安全半徑和速度增量的評(píng)價(jià)值，確定最佳安全著陸點(diǎn)主要步驟如下：

1.評(píng)價(jià)指標(biāo)類型的一致化處理

在已建立的指標(biāo)體系中，指標(biāo)集可能同時(shí)含有“極大型”和“極小型”指標(biāo)，我們分別稱之為優(yōu)質(zhì)因子和劣質(zhì)因子，本文中速度增量為劣質(zhì)因子，安全半徑為優(yōu)質(zhì)因子，因此在評(píng)價(jià)之前必須將評(píng)價(jià)指標(biāo)的類型進(jìn)行一致化處理，即要統(tǒng)一化為極大型指標(biāo).

2.評(píng)價(jià)指標(biāo)的無(wú)量綱化處理

本文采用極差化公式對(duì)指標(biāo)進(jìn)行無(wú)量綱化處理，公式如下：

其中max為數(shù)據(jù)的最大值，min為數(shù)據(jù)的最小值.

3.確定評(píng)價(jià)指標(biāo)權(quán)重

以嫦娥三號(hào)安全著陸為首要任務(wù)，在此基礎(chǔ)上我們考慮令安全半徑和速度增量的權(quán)重為7比3.則此時(shí)綜合評(píng)價(jià)評(píng)分模型為0.7·R*+0.3·Δν*.

4.評(píng)價(jià)模型求解

根據(jù)該綜合評(píng)價(jià)模型利用MATLAB進(jìn)行求解，得到最優(yōu)的著陸點(diǎn)區(qū)域?yàn)閳D像分析時(shí)將原有圖片切割為行×列為23×23時(shí)20行23列的碎片區(qū)域.

（5）精避障階段

精避障主要目的是在粗避障選取的較安全區(qū)域內(nèi)進(jìn)行精確的障礙檢測(cè)，務(wù)必識(shí)別并閃避危及安全的障礙，選擇最佳的著陸點(diǎn)，此時(shí)在這里不再考慮燃料消耗最優(yōu).具體算法流程如下：

Step1. 圖像分析.

首先將在距離月面100 m處懸停拍攝到的數(shù)字高程圖每個(gè)像素點(diǎn)的高度值讀出，將1 000×1 000（0.1 m）的范圍分割成每個(gè)碎片具有范圍40×40（0.1 m），共有行×列為25× 25個(gè)數(shù)的碎片，劃分范圍的依據(jù)為嫦娥三號(hào)著陸器直徑約為4 m，則其降落的安全半徑為20.分別對(duì)每一個(gè)碎片圖像進(jìn)行分析.

Step2. 計(jì)算每個(gè)碎片各自的高度方差diversity（不平坦程度）.

Step3.計(jì)算每個(gè)碎片各自的平均坡度.

坡度slope為兩點(diǎn)的高程差與其水平距離的百分比 其計(jì)算公式如下：

則平均坡度Δslope計(jì)算公式為：

Step4. 建立綜合評(píng)價(jià)模型0.5·Δslope*+0.5·diversity*，根據(jù)高度方差和平均坡度的評(píng)價(jià)值，確定最佳安全著陸點(diǎn).

Step5.評(píng)價(jià)模型求解

根據(jù)該綜合評(píng)價(jià)模型利用MATLAB進(jìn)行求解，我們得到最優(yōu)的著陸點(diǎn)區(qū)域?yàn)閳D像分析時(shí)將原有圖片切割為行×列為25×25范圍為40×40（0.1 m）的6行18列的碎片區(qū)域.

此時(shí)，精避障時(shí)嫦娥三號(hào)運(yùn)動(dòng)狀態(tài)如下：

嫦娥三號(hào)初速度為零，對(duì)其進(jìn)行受力分析，如圖11，建立力學(xué)方程如下

其中νf，T分別為精避障階段嫦娥三號(hào)最終速度及階段所耗時(shí)間，其中x表示嫦娥三號(hào)到安全著陸點(diǎn)的水平平移距離.此時(shí)輸入θ值，即可得出νf，T的結(jié)果.此時(shí)可得速度增量結(jié)果Δν＝at，利用式（16）算出所耗燃料質(zhì)量.

圖10 兩個(gè)階段嫦娥三號(hào)推力方向圖

圖11 精避障階段嫦娥三號(hào)受力分析圖

（6）緩速下降階段

在緩速下降階段，此時(shí)推力方向與重力方向相反，此時(shí)通過(guò)方程：根據(jù)精避障最終速度，可以求出緩速下降階段所耗時(shí)間和速度增量，同樣利用式（16）可算出所耗燃料質(zhì)量.

3.2.2 確定嫦娥三號(hào)著陸軌道

在著陸準(zhǔn)備軌道上，嫦娥三號(hào)由近月點(diǎn)（19.42W，28.51N），主減速階段初速度方向與水平方向所夾銳角為6°.第二階段主減速階段結(jié)果為第一問(wèn)結(jié)果，根據(jù)3.2.1對(duì)快速調(diào)整、粗避障、精避障、緩速下降四個(gè)階段所建立的力學(xué)方程（18）-（21），以快速調(diào)整時(shí)輸入的θ1，φ及精避障時(shí)輸入的θ2為變量，考慮推力大小范圍 [1 500 N，7 500 N]，采用遍歷算法，輸出經(jīng)歷后四個(gè)階段段所用時(shí)間及所消耗燃料.

此時(shí)六個(gè)階段所耗總時(shí)間最?。篢min=661.199 7 s，此時(shí)θ1=0，φ=85.94°，θ2=0.

所耗六個(gè)階段總?cè)剂献钌伲篢min消耗燃料=1 187.914 3 kg，此時(shí)θ1=0，φ=34.37°，θ2=0.

2.3 問(wèn)題三的分析

問(wèn)題三要對(duì)著陸軌道和控制策略做相應(yīng)的誤差分析和敏感性分析，主要由以下步驟組成：

步驟一：建立誤差分析模型，目的對(duì)著陸軌道和控制策略進(jìn)行誤差分析；

步驟二：建立敏感性分析模型，對(duì)著陸軌道和控制策略進(jìn)行敏感性分析.

2.3.1 誤差分析模型的建立與求解

（1） 模型建立[6]

在此我們針對(duì)主減速階段引入8個(gè)常規(guī)慣性測(cè)量單元典型誤差值（數(shù)值結(jié)果表2），利用這些誤差值計(jì)算出主減速段初始測(cè)量量對(duì)應(yīng)的實(shí)際量范圍.定義初始實(shí)際量計(jì)算公式為非初始實(shí)際量計(jì)算公式為其中為實(shí)際量，為測(cè)量量，xi為初始狀態(tài)偏差（包括初始位置及速度偏差），為傳感器測(cè)量偏差（包括位置、速度和加速度測(cè)量偏差），為傳感器刻度因素誤差系數(shù)（包括位置、速度、加速度刻度因素誤差系數(shù)）.此外，當(dāng)測(cè)量量為零時(shí)，不考慮誤差.

表2 常規(guī)慣性測(cè)量單元典型誤差值

（2）實(shí)際求解：

由公式解出主減速段初始實(shí)際量范圍為：[0，15 000±100，1 700±10，0，0].

主減速段終止實(shí)際量范圍為：[469 900±1 040，359±100.718，0，57±10.57，3.008 6± 0.023 2].

2.3.2 敏感性分析模型建立與求解

（1）模型建立

由推斷可以得出主減速階段為單獨(dú)一個(gè)階段，其自變量為初速度入射角γ，因變量為燃料消耗與時(shí)間；剩下的快速調(diào)整姿態(tài)、粗避障、精避障和緩速下降階段為另一一個(gè)整體階段，其自變量為調(diào)整姿態(tài)時(shí)初速度與水平方向的夾角φ（因θ1，θ2最優(yōu)值為零，所以不做敏感性分析）.

針對(duì)每個(gè)階段的一個(gè)指標(biāo)所改變的百分?jǐn)?shù)，保持其余指標(biāo)不變，求出各因變量所改變的百分?jǐn)?shù)，然后求出所有指標(biāo)的測(cè)試結(jié)果，即可得出敏感性分析結(jié)果.

（2）模型求解

由圖12可知，當(dāng)γ增加相對(duì)自己的百分?jǐn)?shù)值時(shí)，燃料損耗敏感性低，基本沒(méi)有影響，但時(shí)間的敏感性高，時(shí)間改變相對(duì)自己的百分?jǐn)?shù)值大概與γ一樣.由圖13可知，φ的百分?jǐn)?shù)改變對(duì)時(shí)間和燃料消耗的敏感性高，呈負(fù)相關(guān)關(guān)系.

圖12 主減速階段：自變量γ對(duì)于因變量的影響圖

圖13 整體階段：自變量φ對(duì)于因變量的影響圖

3 結(jié)論

通過(guò)對(duì)以上三個(gè)問(wèn)題的探討，本文研究了如何對(duì)裝有變推力發(fā)動(dòng)機(jī)的月球著陸器進(jìn)行軟著陸軌道設(shè)計(jì).在問(wèn)題一中，通過(guò)對(duì)主減速階段建立動(dòng)力學(xué)模型，判斷著陸器在著陸準(zhǔn)備軌道上所處的狀態(tài)及地點(diǎn).在問(wèn)題二中，對(duì)主減速階段所建立的模型進(jìn)行優(yōu)化，并且通過(guò)著陸軌道各個(gè)階段嫦娥三號(hào)所處狀態(tài)建立相應(yīng)模型確定各個(gè)階段的最優(yōu)控制策略，同時(shí)通過(guò)建立力學(xué)方程，將各個(gè)階段進(jìn)行合理性整合，確訂最后的著陸軌道.在問(wèn)題三中建立誤差分析及敏感性分析模型，使得該軟著陸軌道設(shè)計(jì)模型具有可推廣性，可應(yīng)用于相關(guān)月球著陸器軟著陸軌道設(shè)計(jì)之中.

[1]王鵬基，張熇，曲廣吉.月球軟著陸重力轉(zhuǎn)彎軌道設(shè)計(jì)與分析 [C]//中國(guó)宇航學(xué)會(huì)深空探測(cè)技術(shù)專業(yè)委員會(huì)第二屆學(xué)術(shù)會(huì)議論文集.北京：[出版者不詳]，2005：328-333.

[2]張洪華，梁俊， 黃翔宇，等.嫦娥三號(hào)自主避障軟著陸控制技術(shù) [J]. 中國(guó)科學(xué):技術(shù)科學(xué)，2014，44（6）：559-568.

[3]王建偉，李興.近日點(diǎn)和遠(yuǎn)日點(diǎn)速度的兩種典型解法 [J].物理教師，2013，34（6）：58.

[4]王曉晨，潘晨，五軒.十問(wèn)嫦娥：揭密嫦娥三號(hào)探測(cè)器 [J].中國(guó)航天，2014（1）：15-26.

[5]林勝勇.一種月球軟著陸優(yōu)化算法 [J].上海航天，2008（4）：31-34.

[6]劉浩敏，馮軍華，崔祜濤，等.月球軟著陸制導(dǎo)律設(shè)計(jì)及其誤差分析 [J].系統(tǒng)仿真學(xué)報(bào)，2009，21（4）：936-943.

Prediction and Optimization Model of Chang’e 3 Soft Landing

WU Yongfeng，YE Ming，YAO Haoqi，WEI Caimin
（Department of Mathematics，Shantou University，Shantou 515063，Guangdong，China）

In this paper,the chang’e 3 soft landing predictions and optimization problems are solved by using the dynamic model,based on the variable thrust engine with spiral search algorithm.First,physics model is established in Question 1 to work out the goddess of the moon landing to No.3 track point velocity near and far.Then,in the main reduction stage,main reduction variable thrust dynamics model is established.The polar angle of chang’e 3 was studied,taking into consideration of the effects of rotation of the moon,landing for orbit determination in recent months,the position of the point and far point,and the goddess of the moon on the No.3 in the near and far points speed direction.In Question 2,optimization model of variable thrust engine fuel consumption is established for the problem of main reduction stage optimization.Comprehensive evaluation model for safety radius and velocity increment index is established.The coarse type of obstacle avoidance is the scope of safe landing.Pure phase of obstacle avoidance is a premise to ensure safety.Evaluation index for the uneven degree and average slope of the comprehensive evaluation model are established to determine the precise landing area.The mechanical model for each stage is established.Finally,the goddess of the moon landing orbit and the optimal control strategy in the six stages are determined.In Question 3,the error model is established for the design of landing trajectory and control strategy for the corresponding error analysis.The single factor sensitivity analysis is used for calculation of independent variable and dependent variable sensitivity in main reduction phase and residual phase 4.

main reduction dynamics model；variable thrust;spiral search algorithm

V 476.3

A

1001-4217（2015）03-0018-13

2014-11-24

吳永峰（1993-），男，福建泉州人，汕頭大學(xué)數(shù)學(xué)系2012級(jí)本科生；

韋才敏（1977-），男，博士，教授，研究方向：隨機(jī)模型理論與方法.E-mail：wcm1234460@163.com

汕頭大學(xué)國(guó)家基金培育項(xiàng)目（NFC12002）