◎李佳彤 于曉航 宋陽

數(shù)學(xué)建模在數(shù)據(jù)處理與仿真的應(yīng)用

◎李佳彤 于曉航 宋陽

本文主要針對(duì)嫦娥三號(hào)著落過程中粗、精避障問題，基于已有算法完成了著陸位置的選取任務(wù)。進(jìn)而優(yōu)化了算法，突破性的預(yù)先排除深坑、大坑等不利地形，避免為月球車進(jìn)一步開展工作造成障礙；同時(shí)開創(chuàng)性的更新了算法流程，上百倍的縮小了程序運(yùn)算時(shí)間。我們采用基于高程圖數(shù)據(jù)的危險(xiǎn)區(qū)掃描識(shí)別方法，通過計(jì)算掃描區(qū)域的粗糙度和坡度來判別地貌、識(shí)別危險(xiǎn)區(qū)，然后掃描計(jì)算每一選定區(qū)域的粗糙度值、坡度值和相對(duì)目標(biāo)著陸點(diǎn)的距離等因素來加權(quán)生成該點(diǎn)著陸成本圖，最后根據(jù)成本圖優(yōu)選著陸點(diǎn)，仿真計(jì)算結(jié)果表明，該方法能夠快速識(shí)別最優(yōu)著陸點(diǎn)。

模型建立

記嫦娥三號(hào)著陸器矩形邊長(zhǎng)為s，為了提高著陸安全性，考慮控制誤差，應(yīng)使著陸點(diǎn)附近控制誤差范圍內(nèi)的地形也符合安全著陸要求，因此，確定最小計(jì)算單位為3s個(gè)網(wǎng)格。按柵格化掃描模式開始掃描，對(duì)每一個(gè)掃描區(qū)：計(jì)算3s個(gè)網(wǎng)格內(nèi)的所有規(guī)則樣點(diǎn)的方差δ，作為該網(wǎng)格的粗糙度。然后將3s個(gè)網(wǎng)格平均劃分成九個(gè)小的單位s*s，取出每一個(gè)小單位的角位置的四個(gè)點(diǎn)，任取三個(gè)點(diǎn)可以組成一個(gè)平面，計(jì)算這個(gè)平面與月球的水平面的夾角α?？梢缘玫剿膫€(gè)坡度值 1α ， 2α，3α，4α，在每一個(gè)小單位中都可以得到相應(yīng)的1iα，2iα ，3iα ，4iα ，求對(duì)應(yīng)的的α均值，作為該掃描區(qū)域的坡度值。以逐行掃描的方式，每次移動(dòng)三個(gè)網(wǎng)格，進(jìn)行相同的計(jì)算過程，當(dāng)計(jì)算完一列后，向右側(cè)移動(dòng)一列，完成整個(gè)區(qū)域的掃描計(jì)算。然后，計(jì)算每一個(gè)網(wǎng)格到預(yù)定著陸點(diǎn)的水平距離，生成一個(gè)距離矩陣R。

在選擇最優(yōu)著陸點(diǎn)之前，將得到的粗糙度圖、四個(gè)坡度圖和距離圖歸一化。根據(jù)給定的著陸器安全著陸的最大粗糙度maxδ、坡度maxα和著陸器最大水平機(jī)動(dòng)距離maxR，若δ（i，j） >δm ax 或α（i，j）>αm ax 或R （i，j） >Rmax ，則將該網(wǎng)格內(nèi)的值設(shè)置成1，并且，以該點(diǎn)為中心的s個(gè)網(wǎng)格內(nèi)的值均設(shè)置成1，以保證著陸器在著陸過程中，不會(huì)與危險(xiǎn)區(qū)中的障礙物相撞。從而，生成歸一化的粗糙度圖、坡度圖和距離圖。

根據(jù)對(duì)預(yù)定的著陸區(qū)域附近地形的初步分析（任務(wù)開始前已由地面地形分析給出），可以設(shè)定粗糙度、坡度和距離在最優(yōu)著陸點(diǎn)選擇中的權(quán)值a，b，c。由公式

得到著陸區(qū)域的著陸成本圖數(shù)據(jù)，最低點(diǎn)即為選出的著陸點(diǎn)。

模型分析

基于上文的理論分析，我們進(jìn)行了仿真模擬。首先考慮到粗避障的精度需求，我們選取了100*100的掃描單位。在附件中給出的“距2400m處的數(shù)字高程圖”中讀取數(shù)據(jù)進(jìn)行掃描，我們最優(yōu)的得到了一塊基于完全平整的區(qū)域，高度值基本為0m，在這塊區(qū)域中粗糙度和坡度基本為零，可以認(rèn)為是完美的落地區(qū)域。但是我們沒有盲目的確認(rèn)落地點(diǎn)而是進(jìn)行了其他方式的測(cè)驗(yàn)，結(jié)果把高程圖重新在MATLAB中繪出三維立體旋轉(zhuǎn)圖后我們發(fā)現(xiàn)，原來我們?cè)谠虑虮砻娴慕德鋮^(qū)域掉到了一個(gè)“大坑”里面，坐標(biāo)為x=1014，y=234（見Fig 1藍(lán)色最深區(qū)域）。這樣的位置即使我們可以確認(rèn)其相當(dāng)“平坦”，但是有兩個(gè)無法避免的弊端：其一月球車一旦降落在其中，很難爬坡出來進(jìn)行勘測(cè)；其二月球車在下降的過程中，很容易撞擊到周圍的“山峰”——紅色環(huán)形。因此我們需要首先在高程地圖上面把這些“大坑”排除。這就是本算法首先值得優(yōu)化也是必須要優(yōu)化的地方。

Fig1　距月2400m處的數(shù)字高程圖在MATLAB下繪制的三維圖像

另外，在運(yùn)算過程中我們進(jìn)一步的優(yōu)化了粗糙度的定義方式以及更新算法流程，將算法運(yùn)行時(shí)間大大減小。

模型優(yōu)化

避開“大坑”。當(dāng)遇到大的隕石坑時(shí)，粗糙度會(huì)變得巨大，可以直接將其剔除，使其不在后來的系數(shù)計(jì)算中被考慮。具體方法是在計(jì)算出各個(gè)位置的粗糙度后，將明顯過大的區(qū)域的D（粗糙度）值設(shè)為無窮大，這樣在后續(xù)的運(yùn)算中就不會(huì)參與計(jì)算，因此剩下參與運(yùn)算的位置就是所有所有不包含大隕石坑的點(diǎn)。這樣的方法將會(huì)良好的避免了之前我們出現(xiàn)的將嫦娥三號(hào)降落在大坑中的情況。

優(yōu)化運(yùn)行時(shí)間。改進(jìn)一：取高程圖中ap,q是從ai,j到中間枚舉的所有元素，每次僅代表一個(gè)元素，范圍在ai,j到ai+3s-1,j+3s-1中變化。在粗避障中，s取基本掃描單位100格，基礎(chǔ)的方差公式為

在公式中由于要每次計(jì)算除以2s，一方面增加運(yùn)算使運(yùn)行速度下降；另一方面，除法產(chǎn)生了小數(shù)，使最終精度會(huì)受到影響。因此我們改進(jìn)了D的定義

這樣成功的降低運(yùn)算速度并且提高精度。

改進(jìn)二：在計(jì)算一個(gè)掃描區(qū)域的粗糙度即方差時(shí)需要用到該區(qū)域的平均值或和值A(chǔ)all，因此需先計(jì)算每個(gè)掃描區(qū)的和值。從Fig 1中我們看到，空心小圓圈表示高程圖中的數(shù)據(jù)點(diǎn)，在“距2400m處的數(shù)字高程圖”中有2300*2300個(gè)數(shù)據(jù)點(diǎn)。紅色小方框代表上一次用來計(jì)算的區(qū)域，黑色小方塊代表下一時(shí)刻計(jì)算和值的區(qū)域，在不斷的平移中我們記錄所有的和值并且按照上式進(jìn)行計(jì)算。在這里發(fā)現(xiàn)其中有重復(fù)的運(yùn)算。為了簡(jiǎn)化，每次我們不重新計(jì)算黑色方框中的數(shù)據(jù)點(diǎn)，而是采用減去前一列再加上后一列的方式，這樣看似簡(jiǎn)單的運(yùn)算實(shí)際上大大簡(jiǎn)化了運(yùn)算復(fù)雜程度。

硬算的話需要時(shí)間復(fù)雜度

可以發(fā)現(xiàn)當(dāng)j2=j+s時(shí)，相同，所以存在重復(fù)計(jì)算，因i<=p

這樣計(jì)算A的總體時(shí)間復(fù)雜度就在預(yù)處理階段，發(fā)現(xiàn)pA（i,j）與pA（i+1,j）只有開頭和末尾的元素不同，存在重復(fù)計(jì)算，因此只需計(jì)算所有pA（1, j），然后

對(duì)于粗避障階段：

對(duì)于精避障階段：

可見均有max（m·s,（n-s）·（m-s））=（n-s）·（m-s）。因此計(jì)算A的總的時(shí)間復(fù)雜度從O（（n-s）·（m-s）·s2）降為了O（（n-s）·（m-s））

改進(jìn)三：受到A的改進(jìn)的啟發(fā)，在接下來計(jì)算粗糙度即方差D時(shí)，希望也能由D（i,j-1）得到D（i,j），優(yōu)化重復(fù)計(jì)算的部分。

為此，將重定義的D展開：

pD（i,j,0）和pD（i,j,1）是預(yù)處理數(shù)組，預(yù)處理復(fù)雜度為O（（n-s）·（m -s）·s ），同改進(jìn)二中對(duì)pE的處理方法，pD可以優(yōu)化到O（（n-s）·（m-s ）），然后可以利用上式在O（（n-s）·（m-s））內(nèi)計(jì)算完所有D，所以整個(gè)復(fù)雜度從O（（n-s）·（m-s）·s2）降為了O（（n-s）·（m-s）·s ）。

通過以上三個(gè)角度的改進(jìn)，算法整體復(fù)雜度從O（（n-s）·（m-s）·s 2）降為了O（（n-s）·（m-s）·s ），進(jìn)而我們進(jìn)行的仿真模擬。原始算法在掃描單位為100*100的條件下需要運(yùn)算1497.431（s），經(jīng)過以上的算法優(yōu)化的方式得到的優(yōu)化后運(yùn)算時(shí)間為60.431（s），時(shí)間減小兩個(gè)量級(jí)，為原來的0.0404。優(yōu)化效果顯著。

權(quán)重因子確定

權(quán)重因子顯然是算法中相當(dāng)重要的一環(huán)。設(shè)置比例合適的粗糙度、坡度、距離三個(gè)因素的權(quán)重因子，計(jì)算每個(gè)不包含隕石坑的掃描范圍的加權(quán)成本值，從中選擇成本最小的位置作為著陸點(diǎn)即可完成粗避障階段的策略選擇。優(yōu)秀的權(quán)重因子的搭配將會(huì)提供降落位置的優(yōu)秀策略，而不合適的權(quán)重因子則會(huì)導(dǎo)致諸如：耗費(fèi)燃料、粗糙度小但陡峭、不陡峭但粗糙等等情況，給降落帶來困難。

關(guān)于權(quán)重因子的確定：為了方便觀察和計(jì)算，三個(gè)加權(quán)值應(yīng)該處于同一個(gè)數(shù)量級(jí)上，所以首先約定三個(gè)影響因素的值與對(duì)應(yīng)權(quán)值的乘積都在0.110量級(jí)上。據(jù)此，初步根據(jù)三個(gè)因素的平均值定義權(quán)值，分別為qD=1e-15，qDis=1e-6，qSlope=1，初次計(jì)算出對(duì)應(yīng)的著陸位置后，三者加權(quán)值為：

qD*D=0.002243 qDis*Dis=1.250293 qSlope*Slope=0.739042

Fig3　粗避障選擇的粗略著陸

比較發(fā)現(xiàn)依然相差過大，很容易導(dǎo)致著陸位置不合適，而此時(shí)找到的位置為：（2179，1398），由圖也觀察到此處并不適合著陸。如圖所示在一個(gè)大坑旁邊，不適合降落。

此階段坡度影響應(yīng)大于粗糙度，大于距離，所以再次調(diào)整三個(gè)權(quán)重，使得qD=1e-13，qDis=1e-7，qSlope=1，再次運(yùn)行程序，此時(shí)結(jié)果為：

qD*D=0.220227 qDis*Dis=0.143785 qSlope*Slope=0.687549

min_cost=qD*D+qDis*Dis+qSlope*Slo pe=0.6875491.051560

著陸點(diǎn)：（1875，185）

由圖可以看到，此位置十分合適。證明了本算法的正確性。

粗避障階段仿真結(jié)果

完成以上改進(jìn)，我們做了完整的粗避障模擬。最終由計(jì)算機(jī)判斷著陸成本Q取最小值初步得到粗略著陸位置x=1875，y=185。成本Q值為0.6875491.051560，權(quán)重系數(shù) qD=1e-13，qDis=1e-7，qSlope=1，選擇 位 置 的 相 對(duì) 粗糙 度 為qD*D=0.220227， 相 對(duì) 坡度為 qSlope*Slope=0.6875491.051560，相對(duì)距離qDis*Dis=0.143785，實(shí)際橫向移動(dòng)距離為1199.1038m，運(yùn)行時(shí)間run time：60.856 s。直觀地從圖中也可看出，在紅旗著陸區(qū)起伏高度（顏色）是最均勻的，不僅避開了所有的“大坑”而且顯然易于降落，也利于降落后進(jìn)一步開展勘測(cè)工作。

Fig4　粗避障選擇的粗略著陸位置

精避障階段仿真結(jié)果

精避障階段的避障策略采用的算法與粗避障階段基本保持一致，但是需要修改的變量有：掃描單位30s=，掃描分為是3s，降落區(qū)域大小為100m*100m，權(quán)重因子也相應(yīng)改變即粗糙度、坡度和距離的權(quán)重分別為0.45、0.45、0.1。最終由計(jì)算機(jī)判斷著陸成本Q取最小值初步得到粗略著陸位置x=191，y=0。降落安全范圍為191≤x≤281，0≤y≤90大小是9平方米的區(qū)域。成本Q值為2.152388，權(quán)重系數(shù)qD=1e-13，qDis=1e-7，qSlope=1，選擇位置的相對(duì)粗糙度為qD*D= 0.832249，相對(duì)坡度為qSlope*Slope= 1.292611，相對(duì)距離qDis*Dis= 0.027528，實(shí)際橫向移動(dòng)距離為52.60 m，運(yùn)行時(shí)間run time： 9.750 s。直觀地從圖中也可看出，在紅旗著陸區(qū)起伏高度（顏色）是最均勻的，由于中心區(qū)域有一個(gè)大坑，不適宜降落，加之我們適度調(diào)低了距離的權(quán)重（消耗燃料極少）。因此模擬結(jié)果最終選擇了一個(gè)看似較為偏遠(yuǎn)，但是粗糙度和坡度搭配最為合理，整體降落成本Q值最小的區(qū)域降落。在圖中由小紅旗標(biāo)出降落地點(diǎn)。

（作者單位：北京師范大學(xué)）

Fig5　精避障區(qū)域示意圖和選擇的精確著陸位置