直方圖匹配算法進行坡度變換的精度評價

胡云華，賀秀斌，畢景芝

（1．中國科學院 水利部 成都山地災害與環境研究所，成都610041；2．中國科學院大學，北京100049；3．中國地質大學（北京），北京100083）

由于在地形分析中的重要作用，DEM數據和基于DEM的坡度、坡長、坡向等地形因子的提取算法及其精度問題多年以來一直是研究的熱點［1－3］。湯國安和劉學軍［4－7］等的研究表明，DEM 網格大小是影響坡度提取結果的重要原因，研究發現，隨著DEM分辨率的下降，地形具有坡度衰減、坡長擴張的效應［8－10］，然而在實際應用中，由于條件和經費的限制，在進行大尺度研究區地形因子提取時，很難獲得高分辨率的DEM數據，因此采用合理的方法進行校正是一種有效提高數據精度的方法。楊勤科等［11］以黃土高原為研究區，提出了一種用數字圖像處理中的直方圖匹配算法進行坡度變換的方法，經驗證，該方法能很好地擬合數據的直方圖，從而提高數據的質量，并推廣到其它地區。但該方法存在的一個問題是只從圖像統計的角度對坡度數據進行校正，在校正過程中丟失了像元的空間位置信息，在對像元值進行處理時，形成一刀切的狀況，處理結果在空間域上的精度還有待驗證，因此在使用該算法前，最好能進行空間域和頻率域的雙重驗證。但很多人在應用時沒有注意到這個問題，郭偉玲等［12］用該算法對黃土丘陵溝壑區進行坡長變換，結果表明，該算法可以很好地修正因DEM分辨率下降引起的坡長擴張效應，而且圖像空間格局基本保持不變，但根據其統計的圖像標準差，校正后的坡長誤差較大。劉前進［13］、史云飛［14］、史彩寧［15］、泮雪芹［16］等都使用了該算法對坡度數據進行處理，但都沒有考慮精度的問題。

本文以山東丘陵區北巖子口溝小流域作為研究區，使用直方圖匹配算法將研究區低分辨率DEM獲得的坡度頻率累積直方圖匹配到高分辨率的坡度頻率累積直方圖上，并對匹配后的圖像分別從統計特征和空間特征方面進行精度評價，為算法的應用和改進提供參考。

1 材料與方法

1．1 研究區概況

山東丘陵分隔為魯東和魯中丘陵兩部，研究區位于山東省魯東丘陵中部北巖子口溝小流域，地處棲霞市西北郊（北緯37°19′59″—37°22′30″，東經120°45′00″—120°48′47″），平均海拔175m，流域面積25．87 km2，該區域為低山丘陵地貌，大小侵蝕溝發育，溝壑切割深度1～20m，侵蝕溝呈明顯的“U”型或“V”型，溝兩側大都有人工修筑的梯田。土地利用類型為林地、草地、梯田、果園等。該地區為中緯度西風帶季風氣候，年平均氣溫12．6℃，四季比較分明，春季多風少雨，空氣比較干燥；夏季氣溫較高并多雨，空氣比較濕熱；秋季較為涼爽少風，雨量較適中；冬季較寒冷降雪量適中，多為偏北風。由于地表土質較為疏松，夏季多暴雨，導致區內季節性土壤侵蝕和水土流失嚴重，溝蝕地貌廣泛發育。

1．2 數據來源和預處理

研究數據來源于國土局1∶1萬DLG數字化地形圖數據，包含等高線、高程點等基本信息。在Arc－GIS 9．3中利用研究區的等高線數據，建立TIN數字高程模型，然后用TIN插值分別獲得地面分辨率為5 m和20m的DEM。利用數字化地形圖中的1 223個已知高程點對插值生成的DEM進行驗證，插值數據的中誤差分別為2．85m和3．82m，均低于國家山地標準5m中誤差要求，說明插值結果是可信的。已有研究表明，1∶1萬地形圖適合制作空間分辨率為5m的DEM圖［17］，本研究過程中，以5m的DEM所獲得的坡度值作為標準數據，20m的DEM獲得的坡度值作為校正數據。為了更好地對數據進行對比分析，利用ArcGIS 9．3的重采樣工具將5m和20 m的兩幅DEM圖像都采樣到10m，分別以SLP 5m和SLP 20m來表示，而校正后的坡度圖以RSLP 20 m表示。

1．3 坡度的計算

使用ArcGIS 9．3的slope工具對DEM求取坡度，ArcGIS 9．3采用的是三階反權距離平方差分法提取x和y方向上的變化率fx和fy，其坡度SLP表示為：

式中：SLP——像元的坡度；fx和fy——在x和y方向上的變化率。

1．4 圖像頻率統計

對兩幅坡度圖像的像元值進行統計，分別得到圖像的頻率直方圖和累積頻率直方圖。從坡度頻率圖（圖1）和累積頻率圖（圖2）上可以看出，兩種不同分辨率的DEM獲得的坡度數據差別很大，從統計值來看，5m的DEM所提取的坡度最大值為54．3°，而20 m的DEM提取的坡度的最大值為37．35°，DEM分辨率降低帶來的坡度變緩效應十分明顯。

圖1 坡度頻率直方圖

圖2 坡度累積頻率直方圖

2 直方圖匹配

2．1 匹配的原理

如圖3所示，實線條代表精度較高的數據的累積頻率直方圖，虛線條代表精度較低的數據的累積直方圖，在x軸上讀取一點X1其對應的累積頻率值是Y1，而此時在精度較高的累積直方圖曲線上，和Y1相等的Y2值對應的坡度值為X2，而直方圖匹配的原理就是由X1→Y1→Y2→X2從而推出X1和X2之間的關系，將待校正的坡度值X1帶入函數關系中，求解出其校正值X2。

圖3 直方圖匹配原理［11］

2．2 兩組坡度函數關系的求解

根據直方圖匹配原理，在Matlab軟件中，利用分段線性插值函數interp1對兩條曲線進行分別擬合，即相當于建立了兩條曲線的函數關系式Y1和Y2，再利用等式Y1（X1）＝Y2（X2）得到關于X1和X2的一組對應關系（圖4）。

圖4 原始坡度和校正值的對應關系

從得到的累積頻率直方圖（圖2）曲線可以看出，兩組數據的累積曲線有一個明顯的交點，這意味著對于精度較低的數據而言，不是做單純的銳化或者平滑處理，而是在坡度較小時，要做平滑處理，在坡度較大時要做銳化處理。從X1和X2的對應曲線（圖4）上，也能看出這種特征，不同的坡度范圍對應著不同的曲線斜率，根據曲線特征可以大致分為四段，每一段都接近于直線，因此可以用一階線性方程對其進行擬合，擬合公式及誤差見表1。

表1 變換數據的分段線性擬合公式及誤差

在ArcGIS Workstation中利用AML語言實現上述分段函數變換（主要用了con條件命令，處理代碼見附件），處理完成后分別對標準坡度圖SLP 5m，待校正坡度圖SLP 20m，和校正后的坡度圖RSLP 20m進行直方圖統計，結果如圖5和圖6所示。

從圖5—6可以看出，經過直方圖匹配算法校正后，RSLP 20m的頻率直方圖和累積頻率直方圖和標準圖像直方圖都非常接近，達到了校正的目的。

圖5 坡度累積頻率直方圖對比

圖6 坡度頻率直方圖對比

3 算法精度評價

3．1 頻率域精度評價

算法在頻率域的精度評價實際上是比較圖像的各統計參數同標準數據的各統計參數的接近程度。利用 ArcGIS 9．3的 Band collection statistics工具，可求得各坡度數據的最小值、最大值、平均值、標準差4個參數。圖像的標準差（standard deviation：STD）是一幅圖像數據離散程度的一種度量，圖像標準差越大，像元數值和整幅圖像的平均值之間差異越大，DEM或者坡度圖的標準差值實際代表了地形的粗糙程度，圖像標準差越小，圖像離散程度越小，地形變化越平緩；圖像標準差越大，圖像離散程度越大，地形起伏變化越明顯，對三幅坡度圖像進行統計分析，結果見表2。

表2 統計分析結果 （°）

從表2可以看出，5m的DEM所提取的坡度值精度明顯高于20m的DEM，20m的DEM對地形進行了平滑，使坡度的最大值從54．30°降到了37．35°，坡度平均值從11．96°降到了10．47°，表征圖像起伏程度的標準差從8．45°下降到6．17°。經過校正后，20mDEM獲得的坡度數據各項指標都有所增加，和標準數據的指標非常相近，有效地減小了由于DEM分辨率下降帶來的地形描述誤差，恢復了地表的粗糙程度。但是這種接近是基于統計值的，只能說在頻率域內算法有很好的校正效果，但算法的精度還需要經過涉及具體空間點的指標的評價。

3．2 空間域精度評價

算法空間域精度評價實際上是針對圖像的每一個像元，比較它變換前后相對于標準數據的接近程度，評價指標有圖像相關性和各種誤差參數等。兩幅圖像之間的相關性大小可以通過協方差和相關系數來表示。兩個變量的協方差表示的是兩個變量總體誤差的方差，若兩個隨機變量X和Y相互獨立，則它們上的協方差為0，如果協方差不為0，說明它們之間存在著一定的關系，而這種關系可以通過相關系數將其進行量化，ArcGIS 9．3中圖層的協方差和相關系數公式分別如下［18］：

式中：Covij——兩個大小相等的圖層i和j的協方差；Z——數據圖層的像元值；μ——圖層數據的平均值；N——像元總數；Corrij——圖層i和j的相關系數。

兩幅圖像數據之間的誤差大小通常還可以用平均誤差ME、中誤差RMSE和誤差標準差STD來進行描述，平均誤差是數據和真實值的誤差的平均值，中誤差是衡量觀測精度的一種數字標準，亦稱“標準差”或“均方根差”，它是相同觀測條件下的一組真誤差平方中數的平方根，誤差標準差體現的是誤差的離散程度，從誤差標準差可以看出數據誤差的波動程度。它們的計算公式分別如下［19］：

式中：ε——實測值與采樣值之間的差值；n——采樣點個數。

對未校正和校正后的20m坡度數據相對于5m坡度數據的協方差、相關系數、平均誤差、中誤差、誤差標準差進行計算，結果見表3。

表3 校正和未校正的坡度數據誤差指標分析結果

從圖像協方差和相關系數來看，雖然經過校正后，坡度的協方差有所提高，但相關系數和原來的未校正的坡度圖大體相當，說明針對空間上某一具體位置點，校正后的坡度值相對于標準數據的精度并沒有得到很好地提高。經過校正后，坡度的平均誤差接近于0，說明校正后的坡度數據從整體上和標準數據非常接近，但是從中誤差和誤差標準差來看，校正數據的中誤差甚至還有所加大，一般認為，兩幅圖像的中誤差是判斷兩幅圖像接近程度的主要指標，經過校正后中誤差沒有有效地降低，說明算法在具體空間點上的校正精度還存在不足。

4 結論和建議

綜上所述，利用不同空間分辨率的DEM進行地面坡度提取時，DEM網格的大小會明顯影響地表的起伏程度，隨著DEM分辨率的減小，所獲得的地形坡度逐漸減小，地形變得平緩。利用數字圖像處理中的直方圖匹配算法，可以有效地對平緩后的坡度數據進行校正，恢復原地表粗糙程度，建立的校正函數關系甚至可以推廣到其它地區。但是應該注意的是，該算法是建立在數據統計特征基礎上的，校正過程拋開了數據的空間位置信息，雖然從整體上恢復了地形的起伏程度，但校正后數據的中誤差和校正前相當，說明算法在具體空間位置上的校正結果還不理想。在實際應用過程中該算法可以用于地形整體粗糙程度的還原，獲取研究區的地形坡度統計信息，用于經驗性的土壤侵蝕模型或水文模型的運算。但不建議將算法用于分布式土壤侵蝕模型的建立，因為土壤侵蝕過程不僅取決于地形還取決于地表植被覆蓋、土壤質地、降雨和水土保持措施等，分布式土壤侵蝕模型是在一個空間位置點同時考慮以上諸多因素進行運算的，該算法處理的結果沒有很好地空間還原精度，將處理后的數據帶入模型所得到的結果不一定可靠。提高算法在每一個像元點校正的精度，是算法需要進一步改進的地方。

附：對SLP 20m數據進行坡度校正的AML代碼如下：

Grid：S1＝con（S＜＝2．29，0，S）

Grid：S2＝con（S＞2．29＆S＜＝5．5，1．715＊S－3．822，S1）

Grid：S3＝con（S＞5．5＆S＜30，1．347＊S－1．675，S2）

Grid：S4＝con（S＞＝30，2．004＊S－21．483，S3）

其中S是待校正的SLP 20m，S4是校正后的RSLP 20m。

［1］ Carter J．The effect of data precision on the calculation of slope and aspect using gridded DEMs［J］．Cartographica，1992，29（1）：22－34．

［2］ Thompson J A，Bell J C，Butler C A．Digital elevation model resolution：Effects on terrain attribute calculation and quantitative soil－landscape modeling［J］．Geoderma，2001，100（1）：67－89．

［3］ 劉學軍，龔健雅，周啟鳴，等．DEM結構特征對坡度坡向的影響分析［J］．地理與地理信息科學，2004，20（6）：1－5．

［4］ 湯國安，龔健雅，陳正江，等．數字高程模型地形描述精度量化模擬研究［J］．測繪學報，2001，3（4）：361－365．

［5］ 湯國安，劉學軍，閭國年．數字高程模型及地學分析的原理與方法［M］．北京：科學出版社，2006：188－189．

［6］ 賈敦新，湯國安，王春．DEM數據誤差與地形描述誤差對坡度精度的影響［J］．地球信息科學學報，2009，11（1）：43－49．

［7］ 王峰，王春梅．地形因子與DEM分辨率關系的初步研究：以蒙陰縣為例［J］．水土保持研究，2009，16（4）：225－229．

［8］ Gao J．Resolution and Accuracy of Terrain Representation by Grid DEMs at a Micro－scale［J］．International Journal of Geographical Information Science，1997，11（2）：199－221．

［9］ David M W，Mccabe G J．Differences in Topographic Characteristics Computed from 100－and 1000－m Resolution Digital Elevation Model Data［J］．Hydrological Processes，2000，14（6）：987－1002．

［10］ 牛亮，楊勤科．DEM尺度變換中直方圖相似度計算與應用［J］．水土保持研究，2010，17（3）：120－125．

［11］ Yang Q K，Jupp D，Li R，et al．Re－scaling lower resolution slope by histogram matching［C］∥Zhou Q M，Lees B G，Tang G A．Advances in Digital Terrain A－nalysis．Berlin：Springer，2008：193－210．

［12］ 郭偉玲，楊勤科，程琳，等．區域土壤侵蝕定量評價中的坡長因子尺度變換方法［J］．中國水土保持科學，2010，8（4）：73－78．

［13］ 劉前進，于興修．北方土石山區土壤侵蝕強度垂直景觀格局：以沂蒙山區為例［J］．地理研究，2010，29（8）：1471－1483．

［14］ 史云飛，張玲玲．魯中南山地丘陵區土壤侵蝕強度景觀格局的動態變化［J］．生態學雜志，2012，31（8）：2059－2065．

［15］ 史彩寧，常慶瑞，王春梅．基于GIS的延安市土壤侵蝕強度等級評價［J］．水土保持研究，2010，17（3）：28－31．

［16］ 泮雪芹，劉占仁，孟曉云，等．云蒙湖流域不同土地利用類型的土壤侵蝕特征分析［J］．水土保持研究，2012，8（4）：6－9．

［17］ 楊勤科，李銳，梁偉．區域水土流失地形因子的地圖學分析［J］．水土保持研究，2006，13（1）：56－58．

［18］ Snedecor G W，Cochran W G ．Statistical Methods．6th ed［M］．Ames，Iowa：The Iowa State University Press，1968．

［19］ 劉飛，范建容，郭芬芬，等．藏北高原區DEM高程與坡度值提取的誤差分析［J］．水土保持通報，2011，31（6）：148－151．