數學建模在數據處理與仿真的應用

◎李佳彤 于曉航 宋陽

數學建模在數據處理與仿真的應用

◎李佳彤 于曉航 宋陽

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

模型建立

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

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

根據對預定的著陸區域附近地形的初步分析（任務開始前已由地面地形分析給出），可以設定粗糙度、坡度和距離在最優著陸點選擇中的權值a，b，c。由公式

得到著陸區域的著陸成本圖數據，最低點即為選出的著陸點。

模型分析

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

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

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

模型優化

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

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

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

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

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

硬算的話需要時間復雜度

可以發現當j2=j+s時，相同，所以存在重復計算，因i<=p

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

對于粗避障階段：

對于精避障階段：

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

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

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

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

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

權重因子確定

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

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

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

Fig3　粗避障選擇的粗略著陸

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

此階段坡度影響應大于粗糙度，大于距離，所以再次調整三個權重，使得qD=1e-13，qDis=1e-7，qSlope=1，再次運行程序，此時結果為：

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）

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

粗避障階段仿真結果

完成以上改進，我們做了完整的粗避障模擬。最終由計算機判斷著陸成本Q取最小值初步得到粗略著陸位置x=1875，y=185。成本Q值為0.6875491.051560，權重系數 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。直觀地從圖中也可看出，在紅旗著陸區起伏高度（顏色）是最均勻的，不僅避開了所有的“大坑”而且顯然易于降落，也利于降落后進一步開展勘測工作。

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

精避障階段仿真結果

精避障階段的避障策略采用的算法與粗避障階段基本保持一致，但是需要修改的變量有：掃描單位30s=，掃描分為是3s，降落區域大小為100m*100m，權重因子也相應改變即粗糙度、坡度和距離的權重分別為0.45、0.45、0.1。最終由計算機判斷著陸成本Q取最小值初步得到粗略著陸位置x=191，y=0。降落安全范圍為191≤x≤281，0≤y≤90大小是9平方米的區域。成本Q值為2.152388，權重系數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值最小的區域降落。在圖中由小紅旗標出降落地點。

（作者單位：北京師范大學）

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