DEM提取坡向变率中的误差特征与消除方法

2013-08-08谢轶群汤国安

地理与地理信息科学 2013年2期

谢轶群，汤国安，2＊，江岭，2

（1.南京师范大学地理科学学院，江苏南京 210023；2.南京师范大学虚拟地理环境教育部重点实验室，江苏南京 210023）

0 引言

坡向变率（Slope of Aspect，SOA）［1］是反映地形在水平方向变化特征的关键指标，也是地形定量指标体系的重要组成部分，在基于DEM的数字地形分析中具有特殊意义［2－10］。坡向变率的提取一般采用对DEM坡向矩阵提取坡度的方法获得。但由于坡向值的取值区间［0，360）无法以渐变的方式表达由360到0的过渡，尤其会夸大在正北方向上的变化率，对该方法的提取结果若不加修正则会出现大量误差（后文中简称SOA误差），不利于后续分析及应用。因此，研究SOA误差的特征，并在此基础上提出完善的消除方法有重要意义。

汤国安等提出了利用正反DEM结合消除误差的方法［11］，并得到广泛应用［5－10］。但仍存在两个问题：1）误差特征的概括尚不全面，事实上，误差不只存在于正北坡位置，只要提取SOA时坡向矩阵中两坡向值之差的绝对值大于180，则会出现误差；2）由于在提取坡度时需要采用相应算法，而原有的误差消除方法并没有考虑到坡度提取算法对结果造成的影响，因此，无法取得理想的效果。本文在对SOA误差成因研究的基础上，通过分析，揭示了正反DEM方法中存在的缺陷，并提出了基于“六分法”的误差消除方法。

1 误差成因与解决途径

1.1 SOA误差的成因

对于地面任何一个微分单元，坡度为过该点的切面的倾角，而坡向为该切面法线在水平面的投影方向。坡向线为矢量线，所以，邻坡坡向线的夹角应为两向量的夹角。由向量夹角的定义［12］可知，两向量夹角的取值范围为［0°，180°］，因此，相邻两面坡的坡向线的夹角不应大于180°；而由于坡向值的取值范围是［0，360），致使相邻坡向值之差的取值范围为［0，360），超出了正确的范围。至于北坡误差，其实为一种达到最大值的特例。

1.2 正反DEM消除误差方法的缺陷

提取坡度时，三阶反距离平方权差分算法的精度较高［13］，且在诸多GIS平台软件中得到了广泛应用，公式如下：

其中，z代表邻域内的栅格值，i、j为邻域内的行、列号，分别自上至下、自左至右逐渐增大。

当用该算法提取SOA时，栅格值为坡向值，因此，结合式（2）、式（3）及本文对SOA误差特征的概括，可以看出当提取SOA时，若fx或fy分子式中对应的三对坡向差值的绝对值存在大于180的情况，则会造成误差。

正反DEM的误差消除方法［11］，可以理解为在每个栅格处分别根据正反DEM提取SOA，并取最小值作为最终结果。图1所示为正反DEM分别对应的坡向，其值相差180，而由于坡向值不能存在负值和大于等于360的情况，因此，正反DEM坡向矩阵中同一位置坡向值的关系如下：

其中，aspectrev、aspectobv分别为正、反DEM对应栅格的坡向值。

图1 正反地形坡向示意Fig.1 Aspects of obverse and reverse terrains

正反DEM的误差修正方法主要存在以下缺陷：（1）求算fx或fy时，忽略了正反DEM对应的式（2）、式（3）分子中3组差值正负关系的不一致性。图2所示为基于正反DEM在相同位置分别提取的局部坡向矩阵，由正、反DEM提取的SOA值分别为SOA1和SOA2。图2a（左）中坡向矩阵对应的X、Y方向上6对差值中，不存在绝对值大于180的错误情况，而图2a（右）的反DEM坡向矩阵中则出现了这种错误。但是，由于X方向上3对差值间的正负关系发生了变化，使差值结合后反而小于正确情况，导致SOA1大于SOA2。根据其SOA求算公式［6］，SOA2被误作正确结果。

（2）求算fx和fy时，忽略了正反DEM对应的式（2）、式（3）分子中6对差值正误关系的不一致性。由于采用三阶差分算法提取SOA时，会同时使用邻域窗口内X、Y方向上的6对差值，因此，必须保证6对差值都在正确范围内（绝对值小于180），并且相互之间的正负关系正确。虽然借助反DEM，修正了提取SOA时正DEM坡向矩阵中存在的错误差值，但同时可能生成新的错误差值（如南坡变为北坡），因此，若提取SOA时，新生成的错误差值与修正的差值处于同一个邻域窗口中，则无法求算正确的SOA。图2b为基于正反DEM在相同位置分别提取的邻域坡向矩阵，其中，左图在前两行对应的差值存在错误；经过反地形处理后，原有误差得以消除，但在第3行与第2列生成了新误差。因此，由两矩阵提取得到的SOA均为错，无法对误差进行修正。

图2 正反地形DEM邻域坡向矩阵Fig.2 Neighborhood aspects in obverse and reverse DEMs

1.3 SOA误差消除方法——六分法

1.3.1 三阶反距离平方权差分算法的拆分通过正反地形DEM结合消除SOA误差的方法，由于没有考虑到式（2）、式（3）中差值间的关系，无法有效地消除误差。因此，须对式（2）、式（3）分子中的三组差值分别进行修正，既要保证每一组差值在数值上正确，又要保证相应差值间的正负关系正确。综合上述分析，本文提出了一种新的误差消除方法——六分法。为避免误差消除过程中各差值的相互影响，确保误差修正的独立性，在计算fx和fy前，首先需要将两式的分子拆分为6个子式：式（5）—式（10），之后逐个子式进行计算与修正，使进行组合时差值的值与正负关系均正确。

1.3.2 差值修正标准的建立计算获得式（5）—式（10）的6对差值后，需要对其建立相应标准进行修正，修正标准包含两部分：

（1）差值间正负关系的修正标准。在邻域坡向矩阵中，式（5）—式（7）对应的坡向值相减顺序均为“右减左”（X 方向），式（8）—式（10）均为“下减上”（Y方向），具有内部的一致性，且在最后对fx和fy进行结合时，会对其分别取平方。因此，可以对X，Y方向上的差值采用相同的修正标准，保证它们的正负关系正确，且对两方向上的差值进行结合时，不需考虑fx和fy之间的正负关系。

如图3所示，坡向值的取值范围为［0，360），并以正北方向为0，沿顺时针方向增加。由于坡向为矢量，具有方向性，因此，计算SOA时，需采用向量方法计算坡向值的差值，而差值间的正负关系也必须利用向量的特征进行推导。将坡向按图3中的规则以向量的形式表示后，式（5）—式（10）中的差值即为被减向量顺时针转到减向量方向时，需经过的绝对值最小的角度α｜min｜（在正确情况下），该角度可以为负值；反之，在差值错误的情况下，其结果所反映的则是绝对值最大的角度α｜max｜，导致错误的原因则是坡向值差值的取值区间［0，360）与向量夹角的取值区间［0，180］不同，且坡向值不能连续地表示坡向的周期性，即无法连续地由360过渡到0。由于α｜min｜与α｜max｜均表示被减向量沿顺时针方向转过的角度，因此，α｜min｜·α｜max｜＜0，且｜α｜min｜｜＋｜α｜max｜｜＝360。

图3 坡向相减示意Fig.3 Subtraction of aspects

本文以α｜min｜作为提取式（5）—式（10）中差值的标准（也可以逆时针方向建立相应的标准）。图3表示了在北坡的两种情况下，方向向量A1和A2间差值的正负性。其中，图3a中｜A1｜∈（0，90），｜A2｜∈（270，360），因此，｜A2｜－｜A1｜＞180。但因A1沿顺时针方向转到A2所经过的α｜min｜∈（－180，0），对该差值进行修正时，需要变换其值的正负性。同理，图3b中｜A1｜与｜A2｜的差值也需要进行正负性的变换，使得式（5）—式（7）与式（8）—式（10）中差值间的正负关系正确。

（2）差值绝对值的修正标准。根据误差产生的原理，差值的绝对值必须在［0，180］区间内。由于差值只能取α｜min｜与α｜max｜之一，其中｜α｜min｜｜∈［0，180］，｜α｜max｜｜∈（180，360），且｜α｜min｜｜＋｜α｜max｜｜＝360，因此，对于所有绝对值大于180的差值，需要以360－｜α｜max｜｜的形式将其转换为｜α｜min｜｜。

综上，由于α｜min｜·α｜max｜＜0 且｜α｜min｜｜＝360－｜α｜max｜｜，对差值的修正可根据式（11）完成。

其中，diff和diff＿r分别表示修正前、后的差值。

1.3.3 差值的结合——提取SOA值根据式（11）完成对6对差值的修正后，可以得到绝对值与正负关系均正确的结果：f′xi，f′yi（i＝1，2，3）；再由式（2）、式（3）得到式（12）、式（13），将f′xi，f′yi分别带入式（12）、式（13），即可得到结合后的正确的fx，fy。最后将fx，fy带入式（1），即可得到最终的SOA提取结果。图4所示的流程图表示了六分法对SOA进行提取与修正的整个过程。

图4 六分法流程Fig.4 Architecture of the senary division method

2 实验结果与分析

由于采用低空间分辨率DEM提取地形因子时受空间尺度不确定性的影响［14］，很多细部的坡向变化无法表达，得到的仅为概括性的坡向结果，误差已经很大［14－17］，因此，通过该坡向数据难以获得高精度的SOA。据此，实验采用5m分辨率的DEM数据为基本信息源，并基于ArcEngine开发提取算法。

通过对两种方法的结果求差，即可得到正反DEM方法结果的误差分布情况。进一步提取差值图的绝对值后，将其中绝对值大于1、大于5、大于10及大于30的点位提取出来（图5），以直观地显示不同级别差值的具体分布状况（黑色为高于阈值的点）。表1给出了不同阈值d所对应区域的栅格数占总栅格数的百分比P（d）及区域内部差值绝对值的均值E。

通过表1可以看出，正反DEM法的结果中仍然存在大量误差。在实验样区内，最大误差值约为70.7，误差区域栅格均值约为5.87。此外，由六分法得出的更加精确的SOA矩阵的平均SOA值约为58.0，相比正反DEM法的57.8，消除了其修正结果中占SOA总值约0.34%的误差。

表1 六分法与正反DEM法的对比Table 1 Comparison between senary division and traditional method

图5 六分法与正反DEM法的差值Fig.5 Differences between the results of senary division method and traditional method

图6中的误差点取的是新旧方法提取结果中差值的绝对值大于1的点位，可以看出，正反DEM法无法消除的误差主要分布在山脊或山谷处，即坡向变率相对较大的位置。其原因主要是在这些位置，邻域坡向矩阵中的差值通常较大，差值间关系更为复杂，正反DEM方法难以消除这些误差。

SOA值较高的区域往往是数字地形分析中更值得关注的部分，在图6中也可看出，正反DEM法的结果中误差较大的点位又趋向分布在这些高值区域。实验样区的统计结果表明，在正反DEM法与六分法提取结果中差值大于5的区域内，80.5%以上的SOA值都大于83且具有13.1以上的平均差值，其栅格数占样区内SOA值大于83的栅格总数的6.8%以上。鉴于较高的SOA值在反映很多地形要素特征时的重要性，表2给出了按不同SOA区间分类的情况下，六分法对正反DEM法结果的改进效果，包括在某区间［a，b）内时，错误栅格数n占区间总栅格数的百分比P（n），区间内的栅格数N占栅格总数的百分比P（N），误差的均值Ew，及区间内误差总值w占全局误差总值的百分比P（w），这里的误差均指误差的绝对值。

表2 六分法在不同SOA区间内的改进效果Table 2 Performances of senary division method in different intervals of SOA

从表2可以看出，正反DEM法提取结果中存在的误差分布特性与图6中一致，集中分布在SOA值较高的区域，特别是当SOA大于80时，其区域内包含的误差值约占误差总值的86.98%，在［60，90）这一区间内的累积值则达到94.82%。此外，在SOA大于80的区域内，有12.7%以上的SOA值存在误差，且其误差绝对值的均值达到5.6以上，因此，六分法在SOA的提取、分析与应用中有重要作用。

3 结语

本文通过对SOA误差的产生原理与消除方法的研究，首先，对SOA误差的成因进行了完善：由于坡向值［0，360）无法以连续的形式表达坡向的周期性，即在360到0时出现数值上的跃变（跳跃间断点），使得在计算坡向值的差值时，会得到大于180°的错误向量夹角［12］。其次，由于正反DEM法更多地关注误差产生的机理，而没有进一步考虑与提取算法结合时会产生的问题，使得该模型难以全面消除误差，且可能生成新的错误，给基于SOA的研究结果带来不确定的影响。本文根据对SOA误差成因及正反DEM法缺陷的深入分析，提出了新的SOA误差修正模型，通过对坡度提取算法的拆分以及相应误差修正标准的建立，实现了对SOA误差更为全面、理想的消除，也使得SOA可以更为广泛、正确的应用于未来的数字地形分析研究中。今后的研究将针对多地形因子对地貌特征表达时的相互影响及关系展开，以进一步完善地形因子的相关研究。

［1］ TANG G A.A Research on the Accuracy of Digital Elevation Models［M］.Beijing，New York：Science Press，2000.

［2］陈楠，林宗坚，汤国安，等.数字高程模型的空间信息不确定性分析［J］.测绘通报，2005（17）：14－17.

［3］王盛萍，张志强，张建军，等.黄土残塬沟壑区流域次生植被物种分布的地形响应［J］.生态学报，2010，30（22）：6102－6112.

［4］周访滨，刘学军.基于栅格DEM自动划分微观地貌形态的研究［J］.武汉理工大学学报（信息与管理工程版），2008，30（2）：173－175.

［5］张勇.黄土高原地面坡谱研究［D］.西北大学，2003.55－75.

［6］何丽丽.基于不同比例尺和分辨率DEM的数字地形分析［D］.西南大学，2007.18－20.

［7］韩海辉.基于SRTM－DEM的青藏高原地貌特征分析［D］.兰州大学，2009.27－29.

［8］周毅，汤国安，王春，等.基于DEM增强黄土典型地貌表达效果的方法研究［J］.测绘通报，2009（11）：34－36.

［9］刘洋，杨晏立，何政伟，等.分水线、汇水线的多分辨率多阈值提取分析［J］.地理空间信息，2010，8（1）：41－44.

［10］陈婷，周汝良，朱大运，等.基于DEM的2种提取地形特征线算法对比研究［J］.林业调查规划，2011，36（6）：1－4.

［11］汤国安，李发源，刘学军.数字高程模型教程［M］.北京，科学出版社，2010.145－149，157－158.

［12］ R.柯朗，F.约翰.微积分和数学分析引论（第一卷，第二分册）［M］.北京：科学出版社，2001.434.