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

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

引 言

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

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

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

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

1 引入坡度變化率的補(bǔ)償方法

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

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

1.1 坡度預(yù)處理

設(shè)地面上某點(diǎn)在x方向上的高程變化率為fx，在y方向上的變化率為fy，則該地面點(diǎn)坡度為β ，β的計(jì)算如式（1）所示

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

圖2 坡度計(jì)算方法示意圖Fig.2 Slope calculation method

1.2 補(bǔ)償模型

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)驗(yàn)證

2.1 實(shí)驗(yàn)區(qū)域與數(shù)據(jù)

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

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

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

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

2.2 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

本文使用平均絕對(duì)誤差（Mean Absolute Error,MAE）和均方根誤差（Root Mean Square Error,RMSE）作為評(píng)價(jià)指標(biāo)。MAE（表示為EMAE)和RMSE（表示為ERMSE）的計(jì)算方法如下

表3 月球補(bǔ)償實(shí)驗(yàn)精度評(píng)價(jià)Table 3 Evaluation of lunar data experiments

表4 火星實(shí)驗(yàn)精度評(píng)價(jià)Table 4 Evaluation of Martian data experiments

3 應(yīng)用研究

3.1 月球多地形整體補(bǔ)償模型及分級(jí)補(bǔ)償模型

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

表5 月球多地形整體補(bǔ)償模型精度評(píng)價(jià)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級(jí)。確定分級(jí)標(biāo)準(zhǔn)后，將月球上的低分尺度坡度數(shù)據(jù)逐一劃入對(duì)應(yīng)的分級(jí)范圍內(nèi)，同樣，在各個(gè)分級(jí)內(nèi)，隨機(jī)抽取70%的數(shù)據(jù)用于模型擬合，其余30%的數(shù)據(jù)用于模型檢驗(yàn)。對(duì)每個(gè)坡度范圍內(nèi)的數(shù)據(jù)分別利用引入坡度變化率的補(bǔ)償方式進(jìn)行模型擬合，即可得到每個(gè)坡度范圍內(nèi)的模型。理論上應(yīng)相對(duì)于整體數(shù)據(jù)所得到的模型更適用于對(duì)應(yīng)坡度范圍的補(bǔ)償工作,因此可在整體補(bǔ)償模型表現(xiàn)不佳的坡度范圍內(nèi)，使用對(duì)應(yīng)的分級(jí)補(bǔ)償模型進(jìn)行補(bǔ)償。模型擬合及模型檢驗(yàn)結(jié)果如表6所示，結(jié)果表明：大部分分級(jí)范圍內(nèi)，整體模型與分級(jí)模型的效果接近，可根據(jù)實(shí)際情況選擇二者之一，但在坡度大于20°時(shí)，整體補(bǔ)償模型效果較差，此時(shí)可輔以分級(jí)模型以得到更優(yōu)的補(bǔ)償結(jié)果。

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

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

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

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

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

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

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

表7 補(bǔ)償前后對(duì)比分析Table 7 Comparison analysis before and after compensation

圖8 HRSC-MOLA補(bǔ)償前后坡度對(duì)比Fig.8 Slope map before and after compensation

圖9 HiRISE坡度以及HRSC-MOLA坡度補(bǔ)償前后對(duì)比Fig.9 HiRISE slope and HRSC-MOLA slope before and after compensation

4 結(jié) 論

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

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

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

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