直方图匹配算法进行坡度变换的精度评价

胡云华，贺秀斌，毕景芝

（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．