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

◎李佳彤 于曉航 宋陽

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

◎李佳彤 于曉航 宋陽

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

模型建立

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

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

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

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

模型分析

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

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

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

模型優(yōu)化

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

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

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

這樣成功的降低運算速度并且提高精度。

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

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

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

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

對于粗避障階段：

對于精避障階段：

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

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

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

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

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

權(quán)重因子確定

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

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

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

Fig3　粗避障選擇的粗略著陸

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

此階段坡度影響應(yīng)大于粗糙度，大于距離，所以再次調(diào)整三個權(quán)重，使得qD=1e-13，qDis=1e-7，qSlope=1，再次運行程序，此時結(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

著陸點：（1875，185）

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

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

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

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

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

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

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

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