胡云華，賀秀斌，畢景芝

（1．中國(guó)科學(xué)院 水利部 成都山地災(zāi)害與環(huán)境研究所，成都610041；2．中國(guó)科學(xué)院大學(xué)，北京100049；3．中國(guó)地質(zhì)大學(xué)（北京），北京100083）

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

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

1 材料與方法

1．1 研究區(qū)概況

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

1．2 數(shù)據(jù)來(lái)源和預(yù)處理

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

1．3 坡度的計(jì)算

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

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

1．4 圖像頻率統(tǒng)計(jì)

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

圖1 坡度頻率直方圖

圖2 坡度累積頻率直方圖

2 直方圖匹配

2．1 匹配的原理

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

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

2．2 兩組坡度函數(shù)關(guān)系的求解

根據(jù)直方圖匹配原理，在Matlab軟件中，利用分段線性插值函數(shù)interp1對(duì)兩條曲線進(jìn)行分別擬合，即相當(dāng)于建立了兩條曲線的函數(shù)關(guān)系式Y(jié)1和Y2，再利用等式Y(jié)1（X1）＝Y(jié)2（X2）得到關(guān)于X1和X2的一組對(duì)應(yīng)關(guān)系（圖4）。

圖4 原始坡度和校正值的對(duì)應(yīng)關(guān)系

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

表1 變換數(shù)據(jù)的分段線性擬合公式及誤差

在ArcGIS Workstation中利用AML語(yǔ)言實(shí)現(xiàn)上述分段函數(shù)變換（主要用了con條件命令，處理代碼見(jiàn)附件），處理完成后分別對(duì)標(biāo)準(zhǔn)坡度圖SLP 5m，待校正坡度圖SLP 20m，和校正后的坡度圖RSLP 20m進(jìn)行直方圖統(tǒng)計(jì)，結(jié)果如圖5和圖6所示。

從圖5—6可以看出，經(jīng)過(guò)直方圖匹配算法校正后，RSLP 20m的頻率直方圖和累積頻率直方圖和標(biāo)準(zhǔn)圖像直方圖都非常接近，達(dá)到了校正的目的。

圖5 坡度累積頻率直方圖對(duì)比

圖6 坡度頻率直方圖對(duì)比

3 算法精度評(píng)價(jià)

3．1 頻率域精度評(píng)價(jià)

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

表2 統(tǒng)計(jì)分析結(jié)果 （°）

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

3．2 空間域精度評(píng)價(jià)

算法空間域精度評(píng)價(jià)實(shí)際上是針對(duì)圖像的每一個(gè)像元，比較它變換前后相對(duì)于標(biāo)準(zhǔn)數(shù)據(jù)的接近程度，評(píng)價(jià)指標(biāo)有圖像相關(guān)性和各種誤差參數(shù)等。兩幅圖像之間的相關(guān)性大小可以通過(guò)協(xié)方差和相關(guān)系數(shù)來(lái)表示。兩個(gè)變量的協(xié)方差表示的是兩個(gè)變量總體誤差的方差，若兩個(gè)隨機(jī)變量X和Y相互獨(dú)立，則它們上的協(xié)方差為0，如果協(xié)方差不為0，說(shuō)明它們之間存在著一定的關(guān)系，而這種關(guān)系可以通過(guò)相關(guān)系數(shù)將其進(jìn)行量化，ArcGIS 9．3中圖層的協(xié)方差和相關(guān)系數(shù)公式分別如下［18］：

式中：Covij——兩個(gè)大小相等的圖層i和j的協(xié)方差；Z——數(shù)據(jù)圖層的像元值；μ——圖層數(shù)據(jù)的平均值；N——像元總數(shù)；Corrij——圖層i和j的相關(guān)系數(shù)。

兩幅圖像數(shù)據(jù)之間的誤差大小通常還可以用平均誤差ME、中誤差RMSE和誤差標(biāo)準(zhǔn)差STD來(lái)進(jìn)行描述，平均誤差是數(shù)據(jù)和真實(shí)值的誤差的平均值，中誤差是衡量觀測(cè)精度的一種數(shù)字標(biāo)準(zhǔn)，亦稱“標(biāo)準(zhǔn)差”或“均方根差”，它是相同觀測(cè)條件下的一組真誤差平方中數(shù)的平方根，誤差標(biāo)準(zhǔn)差體現(xiàn)的是誤差的離散程度，從誤差標(biāo)準(zhǔn)差可以看出數(shù)據(jù)誤差的波動(dòng)程度。它們的計(jì)算公式分別如下［19］：

式中：ε——實(shí)測(cè)值與采樣值之間的差值；n——采樣點(diǎn)個(gè)數(shù)。

對(duì)未校正和校正后的20m坡度數(shù)據(jù)相對(duì)于5m坡度數(shù)據(jù)的協(xié)方差、相關(guān)系數(shù)、平均誤差、中誤差、誤差標(biāo)準(zhǔn)差進(jìn)行計(jì)算，結(jié)果見(jiàn)表3。

表3 校正和未校正的坡度數(shù)據(jù)誤差指標(biāo)分析結(jié)果

從圖像協(xié)方差和相關(guān)系數(shù)來(lái)看，雖然經(jīng)過(guò)校正后，坡度的協(xié)方差有所提高，但相關(guān)系數(shù)和原來(lái)的未校正的坡度圖大體相當(dāng)，說(shuō)明針對(duì)空間上某一具體位置點(diǎn)，校正后的坡度值相對(duì)于標(biāo)準(zhǔn)數(shù)據(jù)的精度并沒(méi)有得到很好地提高。經(jīng)過(guò)校正后，坡度的平均誤差接近于0，說(shuō)明校正后的坡度數(shù)據(jù)從整體上和標(biāo)準(zhǔn)數(shù)據(jù)非常接近，但是從中誤差和誤差標(biāo)準(zhǔn)差來(lái)看，校正數(shù)據(jù)的中誤差甚至還有所加大，一般認(rèn)為，兩幅圖像的中誤差是判斷兩幅圖像接近程度的主要指標(biāo)，經(jīng)過(guò)校正后中誤差沒(méi)有有效地降低，說(shuō)明算法在具體空間點(diǎn)上的校正精度還存在不足。

4 結(jié)論和建議

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

附：對(duì)SLP 20m數(shù)據(jù)進(jìn)行坡度校正的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］ 劉學(xué)軍，龔健雅，周啟鳴，等．DEM結(jié)構(gòu)特征對(duì)坡度坡向的影響分析［J］．地理與地理信息科學(xué)，2004，20（6）：1－5．

［4］ 湯國(guó)安，龔健雅，陳正江，等．?dāng)?shù)字高程模型地形描述精度量化模擬研究［J］．測(cè)繪學(xué)報(bào)，2001，3（4）：361－365．

［5］ 湯國(guó)安，劉學(xué)軍，閭國(guó)年．?dāng)?shù)字高程模型及地學(xué)分析的原理與方法［M］．北京：科學(xué)出版社，2006：188－189．

［6］ 賈敦新，湯國(guó)安，王春．DEM數(shù)據(jù)誤差與地形描述誤差對(duì)坡度精度的影響［J］．地球信息科學(xué)學(xué)報(bào)，2009，11（1）：43－49．

［7］ 王峰，王春梅．地形因子與DEM分辨率關(guān)系的初步研究：以蒙陰縣為例［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ì)算與應(yīng)用［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］ 郭偉玲，楊勤科，程琳，等．區(qū)域土壤侵蝕定量評(píng)價(jià)中的坡長(zhǎng)因子尺度變換方法［J］．中國(guó)水土保持科學(xué)，2010，8（4）：73－78．

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

［14］ 史云飛，張玲玲．魯中南山地丘陵區(qū)土壤侵蝕強(qiáng)度景觀格局的動(dòng)態(tài)變化［J］．生態(tài)學(xué)雜志，2012，31（8）：2059－2065．

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

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

［17］ 楊勤科，李銳，梁偉．區(qū)域水土流失地形因子的地圖學(xué)分析［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］ 劉飛，范建容，郭芬芬，等．藏北高原區(qū)DEM高程與坡度值提取的誤差分析［J］．水土保持通報(bào)，2011，31（6）：148－151．