數(shù)學(xué)建模嫦娥三號(hào)運(yùn)行軌跡及著陸點(diǎn)分析

2014高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編號(hào)專用頁賽區(qū)評(píng)閱編號(hào)（由賽區(qū)組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：賽區(qū)評(píng)閱記錄（可供賽區(qū)評(píng)閱時(shí)使用）：評(píng)閱人評(píng)分備注全國統(tǒng)一編號(hào)（由賽區(qū)組委會(huì)送交全國前編號(hào)）：全國評(píng)閱編號(hào)（由全國組委會(huì)評(píng)閱前進(jìn)行編號(hào)）：

嫦娥三號(hào)軟著陸軌道設(shè)計(jì)與控制策略摘要：根據(jù)題目附錄和文獻(xiàn)[4]中提供的嫦娥三號(hào)的運(yùn)行參數(shù)，利用角動(dòng)量守恒及向量幾何的方法，分別確定了近日點(diǎn)、遠(yuǎn)日點(diǎn)的位置向量和速度向量。與文獻(xiàn)[4]的真實(shí)數(shù)據(jù)比較發(fā)現(xiàn)吻合良好。本文重點(diǎn)關(guān)注優(yōu)化減速控制與著陸點(diǎn)避障兩方面：前者燃耗最大，后者決定著陸成敗。首先，在多重坐標(biāo)變換基礎(chǔ)上，建立了飛行器制動(dòng)的動(dòng)力學(xué)方程。并以燃耗為最優(yōu)化性能指標(biāo)、近月點(diǎn)狀態(tài)為初始條件、著陸點(diǎn)狀態(tài)為終端條件，利用極值原理求解飛行器的著陸軌跡，及其最優(yōu)控制參數(shù)。其次，對(duì)避障階段采集的高程圖采取水平剖分、比較高程方差的方法，解出最優(yōu)降落點(diǎn)。關(guān)鍵詞：軟著陸；最優(yōu)軌道；避障問題重述嫦娥三號(hào)于2013年12月2日1時(shí)30分成功發(fā)射，12月6日抵達(dá)月球軌道，于北京時(shí)間12月14號(hào)在月球表面實(shí)施軟著陸。嫦娥三號(hào)在著陸準(zhǔn)備軌道上的運(yùn)行質(zhì)量為2.4t，安裝在其下部的主減速發(fā)動(dòng)機(jī)能夠產(chǎn)生1500N到7500N的可調(diào)節(jié)推力，其比沖（即單位質(zhì)量的推進(jìn)劑產(chǎn)生的推力）為2940m/s，可以滿足調(diào)整速度的控制要求。嫦娥三號(hào)四周安裝了姿態(tài)調(diào)整的發(fā)動(dòng)機(jī)，在給定主減速發(fā)動(dòng)機(jī)的推力方向后，能夠自動(dòng)通過多個(gè)發(fā)動(dòng)機(jī)的脈沖組合實(shí)現(xiàn)各種姿態(tài)的調(diào)整控制。嫦娥三號(hào)的預(yù)定著陸點(diǎn)為19.51W，44.12N，海拔為-2641m。嫦娥三號(hào)在高速飛行的情況下，為了保證嫦娥三號(hào)能準(zhǔn)確地在月球預(yù)定區(qū)域內(nèi)實(shí)現(xiàn)軟著陸，關(guān)鍵的問題是著陸軌道與控制策略的設(shè)計(jì)。其著陸軌道設(shè)計(jì)的基本要求如下：著陸準(zhǔn)備軌道為近月點(diǎn)15km，遠(yuǎn)月點(diǎn)100km的橢圓形軌道；著陸軌道為從近月點(diǎn)至著陸點(diǎn)，其軟著陸過程共分為6個(gè)階段，要求滿足每個(gè)階段在關(guān)鍵點(diǎn)所處的狀態(tài)；盡量減少軟著陸過程的燃料消耗。

根據(jù)上述的基本要求，建立數(shù)學(xué)模型解決下面的問題：

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

（2）確定嫦娥三號(hào)的著陸軌道和在6個(gè)階段的最優(yōu)控制策略。（3）對(duì)于設(shè)計(jì)的著陸軌道和控制策略進(jìn)行相應(yīng)的誤差分析和敏感性分析。2、問題分析2.1技術(shù)背景月球軟著陸的方式有兩種，分別是：1）從月球轉(zhuǎn)移軌道直至月球軌道，從而實(shí)現(xiàn)軟著陸；2）從月球的停泊軌道變軌到近月點(diǎn)，然后實(shí)現(xiàn)軟著陸。[2]嫦娥三號(hào)要求準(zhǔn)確的在月球預(yù)定區(qū)域內(nèi)實(shí)現(xiàn)軟著陸，它降落的過程為如下幾個(gè)階段：嫦娥三號(hào)先進(jìn)行霍曼變軌，由高度約100km的停泊軌道進(jìn)入到近月點(diǎn)高度約15km的橢圓軌道；到達(dá)近月點(diǎn)時(shí)，探測(cè)器的制動(dòng)發(fā)動(dòng)機(jī)點(diǎn)火，進(jìn)入動(dòng)力下降階段，即上圖中的DA段；在距離月面約2km時(shí)，水平速度減為0，即上圖中的AF段，調(diào)整姿態(tài)后，探測(cè)器以自由落體的方式降落到月面，即上圖中的FL段。2.2著陸過程的主要考慮因素著陸期間主要考慮兩大因素：一、主減速段是用時(shí)最長、推進(jìn)劑消耗最多，因此最小化推進(jìn)劑消耗是該段制導(dǎo)律的主要設(shè)計(jì)目標(biāo)。二、避障階段決定飛行器能否順利著陸，因此統(tǒng)計(jì)分析高程數(shù)據(jù)并采取相應(yīng)橫向推力最關(guān)重要。2.3力學(xué)分析如上圖示，以月心為坐標(biāo)原點(diǎn)，月球自傳軸為OZ軸，OX軸過零經(jīng)度線，利用右手準(zhǔn)則確定OY軸，建立空間直角坐標(biāo)系。嫦娥三號(hào)探測(cè)器繞月運(yùn)動(dòng)中，受到指向月心的引力F引力，以及相背于速度方向的由發(fā)動(dòng)機(jī)產(chǎn)生的推力F引力是在月心坐標(biāo)下考慮。F其參數(shù)為：G為萬有引力常量且G=6.67×10M為月球的質(zhì)量，M=7.349×10m為航天器的質(zhì)量r=（發(fā)動(dòng)機(jī)推力是在飛行器軌道坐標(biāo)系下考慮的：FVe為發(fā)動(dòng)機(jī)的沖力，是由發(fā)動(dòng)機(jī)構(gòu)造、性質(zhì)、燃料共同決定。飛行過程中是不變的。mdmdtF是推力方向，由飛行器的飛行姿態(tài)決定。F因此運(yùn)動(dòng)的加速度為：d2.4坐標(biāo)系2.4.1直角坐標(biāo)系的旋轉(zhuǎn)我們假設(shè)月球是一個(gè)不旋球體，不考慮側(cè)向運(yùn)動(dòng)，我們建立的月球探測(cè)器的數(shù)學(xué)模型，利用燃耗最優(yōu)原則，依據(jù)龐特里亞金最大值原理設(shè)計(jì)了最優(yōu)軌道。首先，我們建立了以O(shè)XYZ為原點(diǎn)的月心慣性坐標(biāo)系，OX軸是月球的0經(jīng)度線；OZ是月球的自轉(zhuǎn)方向；OY是根據(jù)右手坐標(biāo)來確立的。要在月心慣性坐標(biāo)下計(jì)算推力，需進(jìn)行坐標(biāo)旋轉(zhuǎn)變換。a根據(jù)歐拉姿態(tài)角的變化過程，OXYZ經(jīng)過兩次有順序的坐標(biāo)軸旋轉(zhuǎn)，得到兩個(gè)矩陣：先將軌道坐標(biāo)OX1Y1ZM再將之沿繞OX方向旋轉(zhuǎn)角度-βM因此得到總的旋轉(zhuǎn)矩陣：TAX1Y1Z1為原點(diǎn)在探測(cè)器的軌道坐標(biāo)系,AZ1指向從月心到著陸器的延伸線方向，AY1垂直指向運(yùn)動(dòng)方向,AX1按右手坐標(biāo)系確定。制動(dòng)發(fā)動(dòng)機(jī)推力F的方向與探測(cè)器縱軸重合，θ為在月心慣性坐標(biāo)中，推力為：

FF是推力方向，如上圖所示。AZ?推力方向與AZφ是推力在AX1YΘ和φ隨飛行器飛行狀態(tài)而變，那么：F同理，在月心慣性坐標(biāo)下，總的加速度為：a2.4.2直角坐標(biāo)系與經(jīng)緯坐標(biāo)系的變換如上圖所示，以月心為坐標(biāo)原點(diǎn)，月球自轉(zhuǎn)軸為豎軸建立空間直角坐標(biāo)系OXYZ，假設(shè)近月點(diǎn)的空間直角坐標(biāo)為(x,y,z)，用μ和σ分別表示近月點(diǎn)位置的經(jīng)度和緯度，可由以下經(jīng)緯度轉(zhuǎn)換成為直角坐標(biāo)系坐標(biāo)的公式：Z=r?sin3.模型假設(shè)不考慮月球的自轉(zhuǎn)。因?yàn)轱w行器的速度相比月球自轉(zhuǎn)的速度快、著陸時(shí)間快得多。不考慮除月球以外其他天體引力對(duì)飛行器的擾動(dòng)。不考慮著陸區(qū)域附近因地形不同而造成的質(zhì)量不同，進(jìn)而影響引力場(chǎng)。即假設(shè)引力場(chǎng)在同樣高度是均勻的。假設(shè)月球引力場(chǎng)是從月心點(diǎn)發(fā)出的，引力指向月心點(diǎn)。假設(shè)飛行器燃料足夠使用。4、符號(hào)及字母描述5.近月點(diǎn)與遠(yuǎn)月點(diǎn)參數(shù)的確定5.1近月點(diǎn)的參數(shù)根據(jù)文獻(xiàn)資料數(shù)據(jù)[4]，由北京工程研究所和空間智能控制技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室的研究員們?cè)谶@篇文獻(xiàn)中對(duì)嫦娥三號(hào)的發(fā)射與著陸的設(shè)計(jì)做出了詳細(xì)的說明,參與嫦娥三號(hào)設(shè)計(jì)與制造的他們?cè)谶@篇文獻(xiàn)中對(duì)嫦娥三號(hào)在登月過程中的各個(gè)階段的各項(xiàng)參數(shù)都做出了詳細(xì)而精確的推到與證明，通過這篇文獻(xiàn)我們可以得到探測(cè)器近月點(diǎn)的位置信息為：經(jīng)度：19.0464°W緯度：28.9989°N高度：15km俯視姿態(tài)角：85°資料中所給的近月點(diǎn)的速度V=1.6957km/s 5.1.1計(jì)算近月點(diǎn)直角坐標(biāo)由2.4.2中構(gòu)建的方程組，設(shè)點(diǎn)N為探測(cè)器在近月點(diǎn)處的位置，且該點(diǎn)所處在0度經(jīng)線的位置上；已知近月點(diǎn)的位置為（19.0464°W，28.9989°N），高度為月球半徑+15km，著陸點(diǎn)的位置為（19.51°W，44.12°N）可得到近月點(diǎn)的直角坐標(biāo)：XYZ5.2遠(yuǎn)月點(diǎn)的參數(shù)5.2.1用近月點(diǎn)計(jì)算遠(yuǎn)月點(diǎn)位置如上圖所示，由于近月點(diǎn)、月心、遠(yuǎn)月點(diǎn)在同一直線上，已知近月點(diǎn)的位置為經(jīng)度：19.0464°W，緯度：28.9989°N。運(yùn)用投影的方法可得到遠(yuǎn)月點(diǎn)的位置，計(jì)算方法如下：經(jīng)度=180°-近月點(diǎn)經(jīng)度=160.9536°W緯度=28.9989°S5.2.2遠(yuǎn)月點(diǎn)的直角坐標(biāo)由2.4.2中的方程組結(jié)合遠(yuǎn)月點(diǎn)的經(jīng)緯度數(shù)據(jù)可以得到遠(yuǎn)月點(diǎn)的直角坐標(biāo)為X1Y15.3飛行器在近月點(diǎn)和遠(yuǎn)月點(diǎn)的速度5.3.1推導(dǎo)、計(jì)算月球的第一宇宙速率根據(jù)萬有引力定律：GMmr2=mv(月球質(zhì)量M=7.349*1022kg，G=6.67*10-11，月球的半徑r=1.737013由此可以得到月球的第一宇宙速度為：V1=1.68km/s由于遠(yuǎn)月點(diǎn)處于嫦娥三號(hào)探測(cè)衛(wèi)星的變軌點(diǎn)也就是變軌后的軌道與環(huán)月運(yùn)動(dòng)軌道的交點(diǎn)處，因此，遠(yuǎn)月點(diǎn)的速度即為探測(cè)器繞月運(yùn)行的速度，即月球的第一宇宙速度，既得遠(yuǎn)月點(diǎn)的速度為V1=1.68km/s。5.3.2利用角動(dòng)量守恒估算近月點(diǎn)速度首先由角動(dòng)量的定義可知：L其中L為角動(dòng)量，r為運(yùn)動(dòng)軌跡向量，p為動(dòng)量，v為速度向量，m為探測(cè)器的重量。然后，剛根據(jù)角動(dòng)量守恒定律：r×v=常量。已知，遠(yuǎn)日點(diǎn)（初次變軌時(shí)的位置）的路徑長度：r1=100km+月球半徑rr代入數(shù)據(jù)可得v2與資料中所給的v2=1.6957km/s相比較，計(jì)算所得的數(shù)據(jù)與資料所給數(shù)據(jù)產(chǎn)生的誤差（3.8%）<5.3.3計(jì)算近月點(diǎn)速度方向由5.1.1中圖所示，設(shè)點(diǎn)N為探測(cè)器在近月點(diǎn)處的位置，且該點(diǎn)所處在0度經(jīng)線的位置上；點(diǎn)L為探測(cè)器著陸點(diǎn)的位置，由此構(gòu)建方向向量OL，ON以及探測(cè)器的速度方向向量V，已知近月點(diǎn)的位置為（19.0464°W，28.9989°N），高度為月球半徑r+15km，著陸點(diǎn)的位置為（19.51°W，44.12°N）根據(jù)向量的原理整理總結(jié)可得如下方程組：ON?將ON，OL數(shù)據(jù)代入結(jié)合5.1.1中所計(jì)算出的近月點(diǎn)N的空進(jìn)直角坐標(biāo)，運(yùn)用matlab進(jìn)行運(yùn)算后整理可得如下線性方程組：1719.5綜上所述及計(jì)算我們可以得到近月點(diǎn)的速度方向向量為：V6、主減速階段的軌道方程6.1動(dòng)力學(xué)方程根據(jù)2.3中牛頓力學(xué)分析；經(jīng)過2.4.1坐標(biāo)方程變換后，飛行器的動(dòng)力學(xué)方程如下：axaya(其中Q=dmd6.2控制方程整理6.1中的方程，將其打開引入簡記符號(hào)后，系統(tǒng)狀態(tài)方程可表示為：dX其中，a=b=c=cosαsinθsinφ-sinαsinθcosφggg以上式子中α，由上式運(yùn)動(dòng)方程可知，決定飛行器狀態(tài)的變量為：sr=（X,Y,Z）6.3燃耗最優(yōu)控制（1）性能指標(biāo)α，β由以上結(jié)論可得出性能指標(biāo)為：J=既降落過程中消耗的燃料最少的情況。（2）初始狀態(tài)有上文數(shù)據(jù)可得，系統(tǒng)的初始狀態(tài)為：v（3）控制變量Q（θ，φ（4）終止?fàn)顟B(tài)v注：動(dòng)力下降過程中推進(jìn)劑消耗約為1.4t。7.調(diào)整與控制階段的控制策略7.1選擇最優(yōu)著陸地點(diǎn)在粗避障階段，將2.4km時(shí)的高程圖用Matlab進(jìn)行剖分，將其剖分100m*100m的方格。根據(jù)高程圖海拔數(shù)據(jù)，計(jì)算以方格為單位的每塊區(qū)域的方差，方差最小的區(qū)域即為最平整的區(qū)域，即探測(cè)器選擇的著陸區(qū)域。同理，在精避障階段將高程圖剖分為5m*5m的方格，根據(jù)海拔數(shù)據(jù)資料，計(jì)算以方格為單位的每塊區(qū)域的方差，方差最小的區(qū)域即為最平整的區(qū)域，即探測(cè)器的著陸位置。數(shù)據(jù)結(jié)果如下：（1）2400米避障定位：圖示區(qū)域：2300m*2300m

剖分塊尺寸：100m*100m，共529個(gè)塊

最優(yōu)剖分塊的高程方差：2.6882

最優(yōu)剖分塊中點(diǎn)的相對(duì)位置：(250m,1950m)（2）100米避障定位：圖示區(qū)域：100m*100m剖分塊尺寸：5m*5m，共400個(gè)塊最優(yōu)剖分塊的高程方差：13.0012最優(yōu)剖分塊中點(diǎn)的相對(duì)位置：(52.5m,27.5m)參考文獻(xiàn)[1]孫剛.測(cè)繪衛(wèi)星的發(fā)展及技術(shù)現(xiàn)狀.測(cè)繪科學(xué)與工程；Vo1．27，No．1Mar．2007[2]軟著陸./view/63739.htm?fr=aladdin#3；Sep.14[3]楊峰，曹麒麟，段海星．基于DNSBlocklist的反垃圾郵件系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J]．計(jì)算機(jī)工程與應(yīng)用；2003(7)：11·12．45．[4]張洪華,關(guān)軼峰,黃翔宇,李驥,趙宇,于萍,張曉文,楊巍,梁俊,王大軼.嫦娥三號(hào)著陸器動(dòng)力下降的制導(dǎo)導(dǎo)航與控制.中國科學(xué)：技術(shù)科學(xué)；Vol.44,377-384,2014[5]葉培建,彭兢.深空探測(cè)與我國深空探測(cè)展望.中國工程科學(xué),2006,8:13–18[6]KlumppAR.ApolloLunarDescentGuidance,R-695.Cambridge:CharlesStarkDraperLaboratoryofMassachusettsInstituteofTechnology,1971.1–4[7]LeeAY,ElyT,SostaricR,etal.Preliminarydesignoftheguidance,navigation,andcontrolsystemofthealtairlunarlander.AIAA2010-7717[8]McHenryRL,BrandTJ,LongAD,etal.SpaceShuttleascentguidance,navigation,andcontrol.JournaloftheAstronauticalS,1979,27:1–38[9]KlumppAR.Apollolunardescentguidance.Automatica,1974,10:133–146[10]SavagePG.Strapdowninertialnavigationintegrationalgorithmdesignpart1:Attitudealgorithms.JGuidControlDynam,1998,21:19–28[11]王寨,李鐵壽,王大軼.探月衛(wèi)星變軌時(shí)的姿態(tài)控制研究.航天控制,2005,23:11–14[13]王大軼，李鐵壽，馬興瑞．月球最優(yōu)軟著陸兩點(diǎn)邊值問題的數(shù)值解法[J]．航天控制，2000，3：44—49．[14]周凈揚(yáng)，周獲.月球探測(cè)器軟著陸精確建模及最優(yōu)軌道設(shè)計(jì).宇航學(xué)報(bào).V01．28No．6.Nov.2007附錄部分MatLab代碼經(jīng)緯坐標(biāo)與直角坐標(biāo)的轉(zhuǎn)換---------------------------functionp=polar2orth(alpha,beta,radius)%alpha：經(jīng)度，正負(fù)號(hào)（E:+,W:-）%beta：緯度，正負(fù)號(hào)（N:+,S:-）p(1)=radius*cos(beta)*cos(alpha);p(2)=radius*cos(beta)*sin(alpha);p(3)=radius*sin(beta);end（2）2400M處高程圖最優(yōu)落點(diǎn)位置算法------------------------------clear;clcM_img100_src=imread('A2.tif');M_img100=double(M_img100_src);N=50;l=length(M_img100)/N;flat_var=zeros(l);flat_mean=zeros(l);%求方格內(nèi)的均值fori=1:lforj=1:lflat_mean(i,j)=mean(mean(M_img100((i-1)*N+1:i*N,(j-1)*N+1:j*N)));endendfori=1:lforj=1:lsum2=0;forp=(i-1)*N+1:i*Nforq=(j-1)*N+1:j*Nsum2=sum2+(M_img100(p,q)-flat_mean(i,j))^2;endendflat_var(i,j)=sum2/N^2;endendmin_var=min(min(flat_var));[x_m,y_m]=find(flat_var==min_var);x_m=x_m(1)y_m=y_m(1)min_varpos_w=x_m*N-N/2pos_h=y_m*N-N/2imshow(M_img100_src);holdon;plot(pos_w,pos_h,'^w','MarkerSize',10,'MarkerFaceColor','w')plot(pos_w,1:length(M_img100),'w')plot(1:length(M_img100),pos_h,'w')plot(length(M_img100)/2,length(M_img100)/2,'+b','MarkerSize',10)（3）100M處高程圖最優(yōu)落點(diǎn)算法-------------------------clear;clcM_img2400_src=imread('A1.tif');M_img2