插值条件下格网DEM坡度计算模型的噪声误差分析

卢华兴，刘学军，王永君，田梓文

1.东南大学 交通学院，江苏 南京 210018；2.南京师范大学 地理科学学院，江苏 南京 210097；3.国家海洋局 第一海洋研究所，山东 青岛 266100

1 引 言

坡度（slope）是表达地形曲面结构和形态的基本参数之一，也是地表过程模拟、土壤侵蚀模型、土地利用分类等环境模型的基本变量，在地貌、地理、土壤、资源与环境等领域有着广泛的应用［1］。目前在基于格网DEM数据的数字地形分析中，格网DEM的坡度计算误差主要来源于DEM数据误差、DEM 数据结构和分析算法［2］，国内外学者从不同的角度对这一问题进行了探讨和分析，如DEM分辨率对坡度计算精度的影响［3－4］、DEM 误差 在 地 形 分 析 中 的 传 递 分 析［5－6］、不同的地形分析算法的精度评价［7－8］等，然而这些研究的不足之处在于并没有考虑DEM误差的空间自相关特性。事实上，DEM作为地形曲面的离散化表达，由于插值技术的存在，DEM误差在空间上是相关的，而DEM误差的空间自相关性已经被证明存在并极大影响DEM及数字地形分析的不确定性，如文献［9］给出了坡度中误差、DEM中误差及DEM误差相关系数之间的计算公式，并从试验上证明当DEM误差相互独立时，坡度误差可达33%，远大于通过计算值和实测值对坡度误差的估计；文献［10］以随机模拟的方式，通过改变DEM误差的空间自相关大小，发现等高线密度（与面积百分比）随着莫兰指数I增大而降低，而河网密度（百分比）却在I＝0.4附近达到了最小；文献［11］采用指数自相关函数模型，从试验角度证实了DEM误差空间自相关性对坡度精度的影响；文献［2］采用随机模拟并结合试验验证的方式，通过选择两种空间自相关函数模型（线形模型、指数模型）论述在两种空间自相关函数模型下的窗口分析的协方差矩阵坡度误差的影响。然而，以上模拟或试验方法在数字地形分析精度评价方面多数采用经验空间自相关模型，经验模型存在如下缺陷：①空间自相关模型选择及其参数的确定上具有盲目性，无法忠实地描述DEM误差空间自相关性规律；②空间自相关模型在整个数字地形分析精度评价过程中一成不变，无法描述DEM格网的误差空间分异特征。

本文首先从DEM统一插值模型入手，推导插值条件下的格网DEM邻近格网之间噪声误差的局部自相关性模型，从理论上描述了格网DEM噪声误差对坡度精度影响，最后从试验角度具体分析3种典型插值条件下的坡度计算的局部空间自相关性对数字地形分析精度的影响，并评价了不同地形分析算法对DEM噪声误差的抗差能力。

2 格网DEM插值的噪声误差分析

DEM插值是指由采样点的高程信息按照特定的空间自相关性规则去估计待测点高程的数学方法，其目的是缺值估计和散点数据的格网化［12］，目前格网DEM多数由散点或地形特征数据的插值获得，DEM插值过程其实是已知点权重分配的过程，即DEM插值的最终结果是已知点高程向量的一种线性组合［13］，设已知点高程向量z：z1，z2，…，zn（n为已知点个数），则待插点P（xP，yP）的高程zP可表示为

根据公式（1），格网P和Q的高程估值为

上式的k和l表示已知点权重向量，对式（2）写成矩阵形式如下

图1 散点插值Fig.1 Elevation points interpolation

插值条件下，格网DEM的噪声误差是指已知点上的噪声误差在插值函数作用下传递到格网点上的噪声误差，格网DEM上的噪声误差的中误差为格网DEM噪声精度。根据方差协方差传播率，并顾及已知点是独立观测数据，于是格网P、Q的方差协方差阵为

上式为2×2阶矩阵，主对角线元素分别表示格网P、Q的指已知点噪声误差经插值的传递误差，表达了格网P、Q插值精度

而非对角线元素、分别格网P、Q之间的协方差，表达了临近格网P、Q之间的误差的相关性

由上式可知，造成临近格网的噪声误差存在相关性的根本原因是插值过程中共用了已知点，若临近格网在插值过程中没有共用已知点，其协方差为0，表示两者不考虑或不存在相关性。

3 格网DEM坡度计算模型的误差分析

格网DEM坡度计算通常采用邻域分析方法（neighborhood analysis，NA），即以固定尺寸如3×3、5×5或7×7等二维滑动窗口逐行或逐列对格网DEM数据进行扫描，每次滑动过程中，对局部窗口里的格网高程值按照特定差分方式估计偏导数，用于估计中心格网的坡度值。一般的，格网DEM窗口分析算法可以写成局部窗口格网高程值的函数，以3×3窗口分析为例（如图2所示），设9个格网对应的高程值为z0、z1、…、z8，窗口分析模型可以写成如下形式

图2 DEM 3×3邻域窗口Fig.2 3×3DEM local moving window

根据误差传播定理中函数误差的传播公式［14］，格网DEM邻域分析的中误差为

式中，右边括号内的变量表示窗口分析函数对局部格网z＝ ［z0z1…z8］的偏导数，Dz为z＝［z0z1…z8］的方差协方差阵，其主对角线为格网的方差，而非对角线上的元素为格网之间的协方差。本文以坡度算法为研究对象，其他局部窗口分析算法的误差分析与此类同。在坡度计算过程中，坡度提取的噪声误差是格网DEM上的噪声误差在插值函数作用下传递到格网点上的噪声误差，即已知点上的随机误差经两次传递得到的噪声误差，格网上坡度的噪声误差的中误差为坡度噪声精度。根据文献［2］的推导，坡度中误差为

式中，Dz为z＝［z0z1z2z3z4z5z6z7z8］T的协方差阵（公式（10）），p＝APz，q＝Aqz，而Ap和Aq是向量z的差分系数，不同的坡度差分算法有不同的系数，为方便计算，在文献［2］的基础上，这里对4种坡度计算模型中的差分系数作进一步整理，得到如表1的差分系数。

表1 坡度计算模型的差分系数Tab.1 Difference coefficient of slope algorithms

需要指出的是，在坡度中误差的计算过程中（公式（10）），对于Dz的求解，文献［2］采用的是经验协方差函数，本文所采用的是理论协方差模型（公式（4）），在计算邻域窗口任意两临近格网的误差空间相关性时，需要在插值过程中记录插值点的点号和权重，并且在协方差阵计算过程中搜索共用点的点号和权重。

4 试验分析

为了对上述理论进行分析论证，试验选择了机载三维激光测量（Lidar）散点数据，散点的垂直中误差（RMSE）为0.3～0.5m，这里取0.4m作为散点高程精度的试验值。考虑到源数据密度太大，对样区随机抽30%左右的散点，得到2512个点作为试验数据（图3（b））。在插值方法上选择DEM地表重建过程中3种典型的插值方法：不规则三角网（triangular irregular network，TIN）的线性插值（linear interpolation for TIN，LIT）、反距离权（inversed distance weighted，IDW），径向基函数插值（racial basic function，RBF），搜索点数为13个，IDW选取的核函数为一次反距离函数，RBF核函数选择二次曲面函数，并建立TIN（图3（c）），作为LIT插值的源数据。

图3 试验数据Fig.3 Test data

具体的试验过程如下：

（1）数据预处理，散点构建TIN，获得TIN数据，获取散点数据的先验方差D，确定插值参数（搜索点数和插值核函数）。

（2）插值处理，逐行逐列进行插值，每插值一次记录一次散点的权重k、l，并记录散点的ID号。

（3）第1次误差传递，按照公式（6），由权重向量k和先验方差D计算DEM精度，获得DEM精度数据。

（4）局部窗口处理，对格网DEM逐行逐列进行局部窗口处理，包括：① 公共点搜索，搜索窗口范围内任意两格网对在插值过程中共用已知点的ID号，并获取共用点权重l；② 计算协方差矩阵，根据公式7，由k、l计算局部窗口9格网的协方差矩阵Dz。

（5）第2次误差传递，按照公式（10），计算坡度精度，获得坡度噪声精度数据，并进行统计分析。具体计算方法的流程图如图4所示。

图4 坡度精度计算过程Fig.4 Flow chart of DEM accuracy calculation

试验按照上述步骤分别计算3种插值模型下的4种差分算法的坡度精度，并考虑顾及和不顾及噪声误差的邻域空间自相关性两种方式。经过计算，获得了24幅坡度精度数据，如表2所示，表格横向表示顾及和不顾及邻域格网的误差相关性下4种差分算法，纵向表示不同的插值模型。

对表中24幅坡度中误差数据进行统计分析，获取每种插值方法和差分方式下坡度噪声精度的最小值（min）、最大值（max）、均值（mean），结果如表3所示。

考察表1中坡度中误差的均值（图5），发现顾及DEM误差的空间自相关性能明显提高坡度计算模型的精度，而不同的插值方法和差分算法也不同程度地影响坡度噪声精度，具体表现在：在格网DEM插值方法上，IDW方法生成的格网DEM在提取坡度时具有较高的噪声精度，其次依次是RBF和LIT，通过对计算程序的跟踪分析，造成这种现象的原因是IDW采用归一化权重，而其他两种插值方式存在权重大于1或负数的情况；在差分方式上，基于二阶差分的坡度噪声精度最高，其次依次是三阶不带权差分、三阶反距离权差分和三阶反距离平方权差分。这与文献［2］的结论有所不同，其根本原因在于：

表2 格网DEM坡度中误差Tab.2 RMSE of slope extracted from grid DEM

表3 各试验样区坡度中误差统计Tab.3 Statistics of root mean square error for slope （°）

图5 坡度噪声精度统计图Fig.5 Statistics chart of slope noise accuracy

（1）协方差计算模型不同。DEM误差自相关原因在于内插时公共点的存在，与内插邻域大小、公共点数等有关。在文献［2］中，协方差模型采用指数模型和线性模型，这是一种经验模型，其前提是DEM误差自相关应符合指数或者线性变化趋势。而本文采用的是理论模型，其协方差值是通过内插计算模型和原始数据误差计算获得的。通过对试验中的协方差数据进一步分析，发现DEM噪声误差空间自相关分布并不符合线性或者指数模型（如图6所示）。

图6 DEM 3×3窗口噪声误差协方差分布Fig.6 Spatial distribution of noise error covariance at local windows 3×3

（2）DEM噪声误差自相关模型是各向异性的。文献［2］中的指数模型或者线性模型，其关键变量是距离而与方向无关，即DEM自相关程度是各项同性的，本文通过计算获取的协方差具有各项异性（如图6所示），即不同方向的自相关性变异程度是不同的。

（3）在坡度精度评价上，坡度计算精度与参与计算格网点数有关。二阶差分采用周围4个格网点，而三阶差分系列采用的是周围8个格网点（如表1所示），前者相关系数为＝10，而后者为＝36，根据公式（10），参与计算的格网点的数量会影响坡度噪声精度。

（4）通过研究，揭示了DEM噪声误差不是均一的，DEM噪声误差描述需要从全局考虑，即给出DEM误差的空间分布，而不应该只给出中误差这一全局性精度指标。

5 结 论

DEM建立和地形分析过程中邻域操作，使得DEM误差在空间上具有特定自相关性，顾及DEM噪声误差的空间自相关性这不但提高DEM数据本身的质量描述，同时也改善了基于DEM的地形分析的精度。DEM噪声误差的空间自相关性源自于DEM插值过程中共用的散点，在DEM数字地形分析特别是邻域分析不能忽略误差的空间自相关性。本文通过理论推导和试验分析，以坡度为研究对象，在顾及和不顾及DEM噪声误差空间自相关性两种情况下对4种常见坡度计算模型的噪声精度进行分析，本文分析方法的主要特点有：① 顾及了DEM噪声误差的空间自相关性；②采用DEM噪声误差理论空间自相关模型。

经过理论推导和试验分析，有如下初步结论：①顾及DEM误差空间自相关性能够从整体上减少坡度中误差，在数字地形分析精度分析中不能忽略噪声误差的空间自相关性，以提高精度预测准确性；② 坡度计算模型的中误差与格网DEM插值方法和坡度计算过程中的差分算法有关，在差分方式上，二阶差分的精度最高，其次依次是三阶不带权差分、三阶反距离权差分和三阶反距离平方权差分，在格网DEM插值方法上，IDW方法生成的格网DEM在提取坡度时具有较高的坡度，其次依次是RBF和LIT。

以上仅探讨了DEM噪声误差的坡度精度评价，本文进一步的研究工作是DEM逼近误差对地形参数的精度影响分析。

［1］ ZHOU Qiming，LIU Xuejun.Digital Terrain Analysis［M］.Beijing：Science Press，2006.（周启鸣，刘学军.数字地形分析［M］.北京：科学出版社，2006.）

［2］ LIU Xuejun，BIAN Lu，LU Huaxing，et al.The Accuracy Assessment on Slope Algorithms with DEM Error Spatial Autocorrelation［J］.Acta Geodaetica et Catographica Sinica，2008，37（2）：200－206.（刘学军，卞璐，卢华兴，等.顾及DEM误差自相关的坡度计算模型精度分析［J］.测绘学报，2008，37（2）：200－206.）

［3］ STENFAN K.The Effect of DEM Raster Resolution on First Order，Second Order and Compound Terrain Derivatives［J］.Transaction in GIS，2004，8（1）：83－111.

［4］ THIERRY L，STEPHANE J，CHRISTOPHE F.Very High Resolution Digital Elevation Models：Do They Improve Models of Plant Species Distribution？ ［J］.Ecological Modeling，2006，198：139－153.

［5］ OKSANEN J，SARJAKOSHI T.Error Propagation Analysis of DEM－based Surface Derivatives ［J］.Computers and Geosciences，2005，31（2）：1015－1027.

［6］ FLORINSKY I.Accuracy of Local Topographic Variables Derived from Digital Elevation Models［J］.International Journal of Geographical Information Sciences，1998，12（1）：47－61.

［7］ JONES K.A Comparison of Algorithms Used to Compute Hill Slope as a Property of the DEM ［J］.Computer and Geosciences，1998，24（4）：315－323.

［8］ LIU Xuejun.On the Accuracy of the Algorithm for Interpretation Grid－based Digital Elevation Model［D］.Wuhan：Wuhan University，2002.（刘学军.基于规则格网数字高程模型解译算法误差分析与评价［D］.武汉：武汉大学，2002.）

［9］ GOODECHILD M.The Spatial Data Infrastructure of Environmental Modeling［M］.London：Longman，1991.

［10］ LEE J.Digital Elevation Models：Issues of Data Accuracy and Application［EB／OL］.［2010－08－12］.http：∥www.esri.com／library／userconf／proc96.

［11］ HEUVELINK B.Error Propagation in Environmental Modeling with GIS［M］.Boca Raton：CRC Press，1998.

［12］ LI Xin，CHENG Guodong，LU Ling.Comparison of Spatial Interpolation Methods［J］.Advance in Earth Science，2000，15（3）：260－265.（李新，程国栋，卢玲.空间插值方法比较［J］.地球科学进展，2000，15（3）：260－265.）

［13］ LU Huaxing.Research on DEM Error Model［D］.Nanjing：Nanjing Normal University，2008.（卢华兴.DEM误差模型研究［D］.南京：南京师范大学，2008.）

［14］ TAO Benzao，YU Zongchou.Basic Principle of Survey Adjustment［M］.Beijing：Publishing House of Surveying and Mapping，1996.（陶本藻，於宗俦.测量平差基础［M］.北京：测绘出版社，1996.）

［15］ FLEMING M，HOFFER R.Machine Processing of Landsat MSS Data and DMA Topographic Data for Forest Cover Type Mapping［R］.West Lafayette：Purdue University Laboratory for Applications of Remote Sensing，1979.

［16］ SHARPNACK D，AKIN G.An Algorithm for Computing Slope and Aspect from Elevations［J］.Photogrammetric Survey，1969，35：247－248.

［17］ HORN B.Hill Shading and the Reflectance Map［C］∥Proceedings of IEEE，1981，69（1）：14－47.

［18］ UNWIN.Introductory Spatial Analysis［M］.London and New York：Methuen，1981.