坡度變化率建模的低分辨率坡度補償方法

尹 力，葉樂佳，邸凱昌，劉 斌，孫小珠，王長煥，薄 正

（1.中國科學院 空天信息創(chuàng)新研究院 遙感科學國家重點實驗室，北京 100101；2.中國科學院大學，北京 100049；3.中國科學院 比較行星學卓越創(chuàng)新中心，合肥 230026；4.上海宇航系統(tǒng)工程研究所，上海 201109）

引 言

地表坡度作為一項重要的地形參數(shù)，在深空探測中廣泛應用于地形地貌分析[1]、路徑規(guī)劃[2-3]、著陸點選擇[4-5]等科研和工程任務中。坡度圖一般由數(shù)字高程模型（Digital Elevation Model，DEM）提取得到，因此，坡度圖的分辨率和精度受到DEM的限制[6]。各種科研或工程任務均對高質(zhì)量坡度數(shù)據(jù)提出了需求，特別是在著陸器和行星車任務中，如著陸點地形特征分析、行星車導航規(guī)劃等方面都需要高分辨率大范圍的坡度圖。然而，DEM通常由立體影像或激光高度計數(shù)據(jù)生成，以月球和火星為例，激光高度計采樣間隔較大，通常為百米量級，高分辨率立體影像則僅在月球和火星少數(shù)區(qū)域存在[7-8]，導致全球覆蓋的月球和火星DEM都是低分辨率的（格網(wǎng)間距幾十到上百m），高分辨率（米級格網(wǎng)間距）DEM產(chǎn)品僅在局部小范圍內(nèi)可用[9-10]。低分辨率DEM生成的坡度相對平滑，也就是說，在同一區(qū)域內(nèi)，低分辨率DEM生成的坡度要低于高分辨率DEM生成的坡度,即出現(xiàn)了坡度低估現(xiàn)象[11]。因此，現(xiàn)有DEM數(shù)據(jù)生成的坡度產(chǎn)品，往往難以滿足應用要求，尤其是行星著陸探測器可達性分析和巡視導航設計的安全要求。

針對坡度轉(zhuǎn)換及坡度低估問題，國內(nèi)外學者做了眾多研究及探索，研究場景包括地球以及行星。湯國安等[12]分析了在不同地形復雜度條件下DEM的不確定性，建立了多尺度的地形分析模型以及尺度轉(zhuǎn)換模型。王英等[13]以30 mDEM為基礎，利用分形結(jié)合半方差函數(shù)的方法分析低分辨率DEM提取的坡度與高分辨率DEM提取的坡度之間的轉(zhuǎn)換規(guī)律。Zhang等[14]用分形理論，通過變異函數(shù)法定義可提供坡度與空間分辨率之間的關系信息的分形參數(shù)，并討論了分形參數(shù)在不同尺度上的變化，建立了一種低分辨率坡度補償?shù)哪Ｐ?。以西班牙南部不同分辨率的DEM對該模型進行驗證，結(jié)果表明該方法相對于直接從低分辨率DEM導出的坡度在準確性方面有了顯著提升。陳燕等[15]通過坡度轉(zhuǎn)換圖譜的方式，選擇黃土高原典型地貌類型的試驗區(qū)域，對該區(qū)域低分辨率DEM提取出的坡度統(tǒng)計值進行誤差糾正。

針對行星也有相關的坡度補償工作，Wang等[16]以HiRISE影像作為試驗數(shù)據(jù)，提出了一種應用于火星坡度糾正的模型，實驗表明該模型可應用于不同分辨率差異的場景中。Wu等[17]還利用該模型所提出的坡度補償方法對月球上的低分辨率坡度圖進行了補償及分析，驗證了該方法應用于行星場景的可行性與有效性。

分析發(fā)現(xiàn)，該補償方法在坡度突變區(qū)域表現(xiàn)欠佳，大誤差普遍聚集在坡度突變處。本文提出一種改進的低分辨率坡度補償模型，利用拉普拉斯算子提取的待補償坡度的變化率信息作為自變量加入補償函數(shù)中，從而提升坡度補償精度。選取覆蓋多種地形的月球和火星數(shù)據(jù)進行實驗及誤差分析，并與傳統(tǒng)方法做出對比，驗證本文提出方法的有效性與適用性?；诒疚奶岢龅姆椒ㄩ_展應用研究，擬合生成了適用于全月的補償模型，另有坡度分級補償模型作為全月補償模型的補充；對覆蓋“天問一號”的“祝融號”火星車著陸點50 km×50 km范圍內(nèi)的低分辨率坡度進行補償。

1 引入坡度變化率的補償方法

本文提出的低分辨率坡度補償方法流程如圖1所示。首先，選取覆蓋同一區(qū)域的高分辨率DEM（如LROC NAC DEM、HiRISE DEM）和低分辨率DEM（如HRSC-MOLA融合DEM）。從高分辨率DEM中提取高分辨率坡度作為目標值，從低分辨率DEM中提取低分辨率坡度作為待補償值。此外，為了與低分辨率坡度圖逐像素對應，需要對高分辨率坡度進行降采樣，即找到同一低分辨率坡度值與多個高分辨率坡度值之間的對應關系。然后，使用拉普拉斯算子對需要補償?shù)牡头直媛势露葓D進行坡度變化率的計算。最后，擬合目標坡度值與待補償坡度值以及坡度變化率之間的二元線性回歸函數(shù)。

圖1 所提出坡度補償方法的流程圖Fig.1 Flowchart of the proposed slope compensation method

1.1 坡度預處理

設地面上某點在x方向上的高程變化率為fx，在y方向上的變化率為fy，則該地面點坡度為β ，β的計算如式（1）所示

不同的坡度提取方法主要體現(xiàn)在fx與fy的計算方式不同，本文采用三階反距離平方權差分（Horn算法）[18]，其計算方法為

圖2 坡度計算方法示意圖Fig.2 Slope calculation method

1.2 補償模型

文獻[16]所采用的補償函數(shù)如式（3）所示

其中：x為待補償坡度值；y為坡度放大倍數(shù)；t為高分、低分DEM分辨率比值；a,b,c均為擬合參數(shù)。在DEM分辨率比值一定的情況下，即式（3）中的t值為常量時，可寫成

為了更簡潔清晰地表示初始坡度與目標坡度之間的關系，可將式（4）的等號左右同時乘x，將方程轉(zhuǎn)化為式（5）的表達形式，本質(zhì)上是一個一元線性函數(shù)

其中：Z為補償后的低分辨率坡度值；X表示待補償?shù)牡头直媛势露戎?；a,b表示線性擬合公式中的參數(shù)，其中a為斜率，b為截距。

在該模型中，補償后的坡度只與待補償?shù)钠露染€性相關，即默認對應的坡度會隨著DEM分辨率的降低而線性減小。然而，行星表面有一些重要的地形特征，如撞擊坑等，這類地形起伏變化很大，坡度值很大。隨著DEM分辨率的降低，發(fā)生坡度的突變。此時仍使用一元線性補償模型，無法還原原始地形。通過分析文獻[16]中的補償結(jié)果，發(fā)現(xiàn)利用線性補償后的結(jié)果與目標值相減，其殘差與補償前坡度的變化率相關?；谏鲜龇治?，本文提出在已有線性補償函數(shù)中加入坡度變化率項。記為X′，作為一個自變量引入到補償模型中參與擬合。所提模型數(shù)學表達式如式（6）所示

其中：Z為補償后的低分辨率坡度值；X表示待補償?shù)牡头直媛势露葓D；X′表示待補償坡度的變化率；a,b,c表示擬合參數(shù)，可基于最小二乘原理擬合得到。

本文中，坡度的變化率用拉普拉斯算子提取獲得，拉普拉斯算子具備各向同性，即旋轉(zhuǎn)不變性的性質(zhì)，并具有計算簡單等優(yōu)勢[20]。定義二維坡度圖像f(x,y)的拉普拉斯算子為

圖像是由離散像素組成的，其離散形式為

上述公式的右邊實際上是坡度圖像某像素和它周圍的8個像素與圖3所示的模板的乘積。模板左上角的像素坐標為(x?1,y?1)，右下角的像素坐標為(x+1,y+1)。如果使用這個模板滑過圖像并計算每個像素的拉普拉斯算子，這個過程就是使用拉普拉斯算子計算坡度圖像的邊緣信息。

圖3 拉普拉斯算計卷積模板Fig.3 Laplacian convolution template

2 實驗驗證

2.1 實驗區(qū)域與數(shù)據(jù)

針對月球?qū)嶒?，選擇了全月范圍內(nèi)的14個高分辨率LROC NAC DEM產(chǎn)品，可在華盛頓大學（University of Washington）的PDS Geosciences Node網(wǎng)站（https://ode.rsl.wustl.edu/moon/indexProductSearch.aspx）下載得到。從14個DEM中提取的14幅坡度圖及其對應的地形類型和編號如圖4所示，它們分布在整個月球表面，地形特征包括高地、月海、撞擊坑、盆地和過渡帶等。直接下載的DEM的原始分辨率從2、3到5 m不等，本實驗為了統(tǒng)一分析方便，首先將下載的原始數(shù)據(jù)均重采樣為5 m。

圖4 月球坡度實驗數(shù)據(jù)Fig.4 Experimental data of lunar slope

針對火星實驗，選取了中國2020年發(fā)射的“天問一號”火星任務搭載的“祝融號”火星車著陸區(qū)烏托邦平原的6幅1 m分辨率的HiRISE DEM數(shù)據(jù)，可從亞利桑那大學（University of Arizona）的月球與行星實驗室網(wǎng)站（https://www.uahirise.org/dtm/）下載得到。這6個區(qū)域的地形非常接近，提取的坡度見圖5。

圖5 火星坡度實驗數(shù)據(jù)Fig.5 Experimental data of Martian slope

2.2 實驗結(jié)果與分析

對于同一區(qū)域，隨著DEM分辨率的降低，坡度低估現(xiàn)象逐漸明顯，圖6展示了上述月球區(qū)域5的一個示例。在圖6第1行中，將2 m分辨率的LROC NAC DEM依次降采樣至20、100和200 m，在第2行，展示了從這些DEM提取出的坡度圖?？梢钥闯?，隨著DEM的降采樣，坡度出現(xiàn)了明顯的減小趨勢。

圖6 坡度低估現(xiàn)象Fig.6 Phenomenon of slope reduction

為了證明所提出的考慮坡度變化率的補償方法的效果比傳統(tǒng)線性方法有所提高，分別采用式（5）中的線性補償方法和式（6）中考慮坡度變化率的方法對14個月球區(qū)域的坡度數(shù)據(jù)進行參數(shù)擬合，每組數(shù)據(jù)包含一個高分辨率的DEM和一個低分辨率的DEM，其中低分辨率DEM是高分辨率DEM降采樣的模擬DEM。第1步，分別從高分辨率和低分辨率DEM中獲取坡度圖；第2步，將高分辨率DEM提取的坡度圖降采樣到與低分辨率坡度圖相同的分辨率；第3步，計算低分辨率坡度的拉普拉斯濾波圖；最后，使用70%比例的數(shù)據(jù)進行回歸，30%比例的數(shù)據(jù)進行檢查。擬合結(jié)果見表1。

表1 月球數(shù)據(jù)擬合結(jié)果Table 1 Fitting results of lunar data

同理，對6個火星區(qū)域的坡度進行參數(shù)擬合，擬合結(jié)果見表2。

表2 火星數(shù)據(jù)擬合結(jié)果Table 2 Fitting results of Martian data

本文使用平均絕對誤差（Mean Absolute Error,MAE）和均方根誤差（Root Mean Square Error,RMSE）作為評價指標。MAE（表示為EMAE)和RMSE（表示為ERMSE）的計算方法如下

表3 月球補償實驗精度評價Table 3 Evaluation of lunar data experiments

表4 火星實驗精度評價Table 4 Evaluation of Martian data experiments

3 應用研究

3.1 月球多地形整體補償模型及分級補償模型

上文圖4所示的這14幅數(shù)據(jù)覆蓋了不同的地形特征，可以作為月球多種典型地形坡度的代表數(shù)據(jù)。從這14個區(qū)域內(nèi)獲得有效數(shù)據(jù)共45 578 385組，隨機選取有效數(shù)據(jù)的70%作為模型訓練數(shù)據(jù)，30%作為測試數(shù)據(jù)，采用引入坡度變化率的坡度補償方法，擬合得到月球多地形整體坡度補償模型。該補償模型設定高分尺度為5 m，低分尺度為20 m。表達式見式（11），評價結(jié)果如表5所示。

表5 月球多地形整體補償模型精度評價Table 5 Evaluation of compensation model for global moon

將全月坡度分為0°～3°，3°～6°，6°～9°，9°～12°，12°～15°，15°～20°，20°～30°，30°以上，共8級。確定分級標準后，將月球上的低分尺度坡度數(shù)據(jù)逐一劃入對應的分級范圍內(nèi)，同樣，在各個分級內(nèi)，隨機抽取70%的數(shù)據(jù)用于模型擬合，其余30%的數(shù)據(jù)用于模型檢驗。對每個坡度范圍內(nèi)的數(shù)據(jù)分別利用引入坡度變化率的補償方式進行模型擬合，即可得到每個坡度范圍內(nèi)的模型。理論上應相對于整體數(shù)據(jù)所得到的模型更適用于對應坡度范圍的補償工作,因此可在整體補償模型表現(xiàn)不佳的坡度范圍內(nèi)，使用對應的分級補償模型進行補償。模型擬合及模型檢驗結(jié)果如表6所示，結(jié)果表明：大部分分級范圍內(nèi)，整體模型與分級模型的效果接近，可根據(jù)實際情況選擇二者之一，但在坡度大于20°時，整體補償模型效果較差，此時可輔以分級模型以得到更優(yōu)的補償結(jié)果。

表6 分級坡度擬合和檢驗結(jié)果Table 6 Fitting and validating results of graded lunar slopes

3.2 “天問一號”著陸區(qū)坡度分析

中國發(fā)射的“天問一號”火星探測器于2021年5月15日成功著陸于火星的烏托邦平原，著陸器攜帶的“祝融號”火星車已在火星表面開展巡視探測任務。烏托邦平原地表相對平坦，撞擊坑和石塊等分布較少，據(jù)推斷，該地區(qū)近地表含有大量水冰，對該地區(qū)的研究，對于了解水在火星演化中的作用甚至是潛在宜居性具有重要的意義[21]。著陸區(qū)的坡度對火星車的行動能力具有重要影響，研究人員使用MOLA DEM提取坡度數(shù)據(jù)發(fā)現(xiàn)，坡度大的區(qū)域多位于撞擊坑的邊緣和高原與平原的交界處，而烏托邦平原內(nèi)部則坡度較小[22]。大尺度坡度分析可服務于著陸區(qū)地質(zhì)概況研究，然而火星車行進可達性分析則需要更加精細的坡度。

在著陸點附近并沒有大范圍的高分辨率HiRISE DEM，選取覆蓋“天問一號”著陸點50 km×50 km范圍的200 m分辨率的HRSC-MOLA融合DEM，基于該數(shù)據(jù)，提取得到200 m分辨率的坡度圖，將區(qū)域內(nèi)僅有的兩幅小范圍HiRISE DEM數(shù)據(jù)制作生成坡度圖，作為參照值進行補償精度的評價，上述實驗數(shù)據(jù)的覆蓋范圍如圖7所示，圖中大范圍數(shù)據(jù)為HRSC-MOLA融合DEM，黃色框和綠色框中的小范圍數(shù)據(jù)為HiRISE DEM，紅色五角星為“祝融號”火星車著陸點。

圖7 實驗數(shù)據(jù)范圍Fig.7 Experimental data coverage

采用該區(qū)域內(nèi)的兩幅HiRISE DEM及其降采樣的DEM作為擬合數(shù)據(jù)，同2.2節(jié)的步驟，擬合生成的適用于該區(qū)域的200 m分辨率的坡度補償模型為式（12）

采用該式對50 km×50 km范圍的HRSC-MOLA坡度進行補償，以HiRISE數(shù)據(jù)作為標準值，對補償結(jié)果進行評價，評價結(jié)果見表7。補償前后的坡度圖如圖8所示，兩幅HiRISE DEM覆蓋區(qū)域的數(shù)據(jù)對比如圖9所示。由圖8及圖9可見，通過補償，坡度得到了整體放大，提升了地形模型的細節(jié)豐富度。

表7 補償前后對比分析Table 7 Comparison analysis before and after compensation

圖8 HRSC-MOLA補償前后坡度對比Fig.8 Slope map before and after compensation

圖9 HiRISE坡度以及HRSC-MOLA坡度補償前后對比Fig.9 HiRISE slope and HRSC-MOLA slope before and after compensation

4 結(jié) 論

提出了一種改進的低分辨率坡度補償方法，在傳統(tǒng)的一元線性補償方法基礎上，考慮到坡度變化率的影響，使用拉普拉斯算子提取低分坡度變化率信息，并將其作為自變量之一，嵌入至補償模型中。主要得到了以下結(jié)論。

1）以覆蓋全月多種地形的14個區(qū)域以及“天問一號”著陸區(qū)烏托邦平原的6個火星區(qū)域作為實驗區(qū)，利用從高分辨率的LROC NAC DEM和HiRISE DEM提取出的坡度作為參考值，逐區(qū)域利用引入坡度變化率的補償模型與已有的一元線性模型對從降采樣的LROC NAC DEM和HiRISE DEM提取出的坡度進行補償。結(jié)果表明：不論何種地形特征，不論月球場景還是火星場景，利用引入坡度變化率的補償方法所得到的精度均優(yōu)于傳統(tǒng)的線性補償方法。本文方法也可用于其它類型地表的坡度補償，由于數(shù)據(jù)限制，本方法只在月球和火星場景下進行了驗證。

2）利用本文的補償模型，提出了適用于月球多種地形的低分辨率坡度補償方法，并給出20 m分辨率DEM對應的坡度補償?shù)? m量級對應坡度的整體補償模型參數(shù)。結(jié)果表明：應用整體補償模型至分級的坡度補償時，在低坡度區(qū)域，整體補償模型與分級模型效果接近，而在高坡度區(qū)域，整體補償模型表現(xiàn)相對較差，此時可輔以分級模型作為整體模型的補充。

3）對覆蓋“祝融號”火星車著陸點50 km×50 km的HRSC-MOLA低分辨率坡度進行補償，與該范圍內(nèi)的一幅1 m分辨率的HiRISE坡度進行比較，補償比例達到80%以上，補償后的MAE和RMSE均小于1°，與補償前相比，地形細節(jié)豐富度有一定程度的增加，再次證明了該方法的有效性。