-
摘要: 为了降低采样点水平和高程误差对数字高程模型（digital elevation model，DEM）建模精度的影响，受总体最小二乘算法启发，以较高精度的多面函数（multiquadric function，MQ）为基函数，发展了整体最小二乘MQ算法（MQ-T），并分别借助数值实验和实例分析验证模型计算精度。数值实验中，以高斯合成曲面为研究对象，设计了受不同误差分量影响的采样数据，借助MQ-T曲面建模，并将计算结果与传统MQ进行比较。结果表明，当采样点仅受高程误差分量影响时，MQ-T计算结果精度与MQ相当；当采样数据受水平误差分量影响时，MQ-T计算结果中误差小于MQ中误差。实例分析中，以全站仪获取的采样数据为研究对象，借助MQ-T构建测区DEM，并将计算结果与传统插值算法进行比较，如反距离加权（inverse distance weighted，IDW）法、克里金（Kriging）法和澳大利亚国立大学DEM专用插值软件（（Australian National University DEM，ANUDEM）法。精度分析表明，随着采样点密度降低，各种插值算法精度逐步降低；不管采样密度多少，MQ-T计算精度始终高于传统插值算法；对山体阴影图分析表明，MQ-T相比Kriging法有一定峰值削平现象。Abstract: Motivated by the idea of total least squared method, a total error-based multiquadric method (MQ-T) has been developed to decrease the effect of both horizontal and vertical errors inherent in sample points on surface modeling. Two examples including a numerical test and a real-world example were, respectively, employed to test the robustness of MQ-T to sample errors. The numerical test indicates that when sample points are only subject to vertical errors, MQ-T has a similar performance to MQ. When sample points are subject to horizontal errors, MQ-T is more accurate than MQ. In the real-world example, MQ-T was used to construct DEMs with sample points collected by a total station instrument, and its accuracy was compared with those of classical interpolation methods including inverse distance weighting, ordinary Kriging (Kriging) and ANUDEM. Results indicate that with the decrease of sample density, the interpolation accuracies of all methods become lower. Regardless of sample density, MQ-T is always more accurate than the other methods. Yet, compared with Kriging, MQ-T has a peak-cutting problem.
-
Keywords:
- DEM /
- positional error /
- accuracy /
- interpolation /
- multiquadric method
-
数字高程模型(digital elevation model, DEM)是对地球表面地形地貌的离散化数字表达，其主要借助插值算法对不规则采样数据插值完成[1-2]。受仪器噪声、人为操作不当以及外界环境不理想等因素影响，DEM源数据不可避免地含有采样误差，即误差是空间采样数据的基本属性[3-6]。插值过程中，采样误差能够被传播和放大，进而显著影响DEM及其派生品的精度。为此，众多研究人员建议采用光滑插值算法代替准确插值用于曲面建模。例如，文献[7]使用克里金(Kriging)法值移除地球物理数据中的噪声，实现光滑地学曲面建模。文献[8]认为薄板样条算法具有较强光滑性，可以有效去除采样数据中误差。鉴于不规则三角网(triangulated irregular network, TIN)线性插值的理论缺陷，文献[9]借助二元二次多项式曲面加权平均实现了DEM光滑插值。考虑到航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission, SRTM)数据受噪声影响，文献[9]利用Kriging值完成SRTM重采样。为了抑制采样数据中噪声的影响，文献[10]发展了高精度曲面建模的光滑插值算法。理论上，光滑插值算法仅降低采样数据中的高程误差分量对曲面建模的影响，而忽略了水平误差分量。尽管仅移除高程误差分量可实现DEM光滑建模，但水平误差分量会导致DEM及其派生品(如地形和水系特征线)产生平面位置误差[11]，进而降低DEM在洪水风险分析[12]、山体滑坡[13]、遥感影像配准[14]等领域数字模拟结果的精度。因此，在水平位置精度要求较高的应用中，DEM建模时需同时考虑采样数据中的水平和高程误差。
总体最小二乘算法(total least squares, TLS)是近30年发展起来的一种能同时顾及系数矩阵和观测向量误差的数据处理方法[15]。在曲线(面)拟合中，TLS通过极小化所有采样点到待定曲线(面)的正交距离之和，来降低采样数据水平和高程误差分量对模拟结果的影响。大量实例分析表明，当采样点受水平和高程误差分量影响时，TLS模拟结果精度高于最小二乘算法[16-17]。实际应用中，TLS自身并不能作为插值算法直接用于DEM建模。传统插值算法中，多面函数(digital elevation model, MQ)以其高计算精度已成为DEM建模的主要技术手段[18-19]。因此，本文以MQ为基函数，借助TLS处理误差思路，发展总体最小二乘MQ算法(MQ-T)，以降低采样点水平和高程误差分量对模拟结果的影响，进而实现DEM高精度建模。
1 MQ-T原理
1.1 传统MQ
在曲面建模中，MQ基于数值逼近理论，假设任意平滑的规则或者不规则曲面均可以由简单曲面相互叠加进而逼近到任意精度等级[20]。MQ求算未知点(xi, yi)函数值的解析表达式为[21]：
$$ f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^n {{w_j}\varphi \left( {{r_{ij}}} \right)} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_i},{y_i}} \right)} $$ (1) 式中，m表示多项式阶数，本文取m=1；n表示采样点数；pj(xi, yi)和bj分别表示多项式的项和对应权重系数；rij表示待求点到第j个采样点的距离；φ(r)和w分别表示映射函数和对应权重。因此，基于式(1)求算未知点函数值f(xi, yi)的前提是已知b和w。
因采样点总数小于未知参数个数(n < n+m)，直接利用最小二乘算法无法解算未知参数值。因此，可借助岭回归推导MQ基本方程组。岭回归实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、降低精度为代价获得回归系数更为符合实际、更可靠的回归方法[22]。岭回归目标函数是在最小二乘方法目标函数基础上加一正则项。基于岭回归获得MQ目标函数为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i},i = 1 \cdots n \end{array} \right. $$ (2) 式中，fi=f(xi, yi)；w=[w1…wn]T，b=[b0…bm]T；e=[e1…en]T，表示采样点高程误差；φi=[φ(ri1)]…φ(rin)]; pi=[p0(xi, yi)…pm(xi, yi)]；c表示光滑参数。
为了求式(2)极小值，由拉格朗日乘数法构建辅助函数为：
$$ \begin{array}{l} L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{e}} + \\ \sum\limits_{j = 1}^n {{\alpha _j}\left( {{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} - {e_i}} \right)} \end{array} $$ (3) 式中，α表示拉格朗日算子。
对L(w, e, b, α)求关于w、e、b和α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{e}}_i}}} = 0 \Rightarrow c{\mathit{\boldsymbol{e}}_i} = {\mathit{\boldsymbol{\alpha }}_i},i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \\ \;\;\;\;\;\;{\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {e_i} = {f_i},i = 1 \ldots n \end{array} \right. $$ (4) 式(4)中消去w和e得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (5) 式中，P=(pij)，pij=pj(xi, yi)；α=[α1…αn]T；f=[f1…fn]T；K=(kij)为MQ径向基函数，由核方法推导而来，即kij=k(rij)=φiφjT。文献[23]分析认为，三次曲面作为核函数其精度优于Hardy所采用的双曲面法，故本文采用k(r)=1+r3。
将$ \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\alpha _i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} $代入式(1)得：
$$ \begin{array}{l} f\left( {{x_i},{y_i}} \right) = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \\ \;\;\;\;\; = \sum\limits_{j = 1}^m {{\alpha _j}{\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{\varphi }}_j^{\rm{T}}} + \sum\limits_{j = 0}^m {{b_j}{p_j}\left( {{x_j},{y_i}} \right)} \end{array} $$ (6) 其矩阵形式为：
$$ \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{K\alpha + Pb}} $$ (7) 因此，基于MQ插值计算时必须首先解算方程组(5)获得α和b的值，然后基于式(6)或式(7)求算未知点函数值。由以上推导可见，MQ计算中仅考虑采样点高程误差对模拟结果的影响，而忽略了水平误差分量。当采样点同时含有水平和高程误差时，此方案势必使得MQ插值结果产生平面位置误差。
1.2 MQ-T
理论上，式(1)可看作为原始特征空间中的超平面方程，其法向量为(w －1)。因此，基于整体最小二乘算法思想，点的平面位置误差为点到超平面的垂直距离。对于MQ而言，采样点平面位置误差可表达为：
$$ {\varepsilon _i} = \frac{{{f_i} - {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} - {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}}}}{{\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} }} $$ (8) 因此，MQ-T的目标函数表达为：
$$ \left\{ \begin{array}{l} \mathop {\min }\limits_{\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}}} F\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }}} \right) = \frac{1}{2}{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}} + \frac{c}{2}{\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}\mathit{\boldsymbol{\varepsilon }}\\ {\rm{s}}.\;{\rm{t}}.\;{f_i} = {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + {\mathit{\boldsymbol{p}}_i}\mathit{\boldsymbol{b}} + {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n \end{array} \right. $$ (9) 类似MQ推导过程，对式(9)由拉格朗日乘数法构建辅助函数，求其关于w、ε、b、α的导数并等于零得：
$$ \left\{ \begin{array}{l} \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{w}}}} = 0 \Rightarrow \mathit{\boldsymbol{w}} = \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{\alpha }}_i}\mathit{\boldsymbol{\varphi }}_i^{\rm{T}}} \\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{\varepsilon }},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\varepsilon _i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} {\varepsilon _i}\sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} ,i = 1 \cdots n\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial \mathit{\boldsymbol{b}}}} = 0 \Rightarrow \sum\limits_{i = 1}^n {{\mathit{\boldsymbol{p}}_i}{\mathit{\boldsymbol{\alpha }}_i}} = 0\\ \frac{{\partial L\left( {\mathit{\boldsymbol{w}},\mathit{\boldsymbol{e}},\mathit{\boldsymbol{b}},\mathit{\boldsymbol{\alpha }}} \right)}}{{\partial {\mathit{\boldsymbol{\alpha }}_i}}} = 0 \Rightarrow {\mathit{\boldsymbol{\varphi }}_i}\mathit{\boldsymbol{w}} + \sqrt {1 + {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}} = {f_i},i = 1 \cdots n \end{array} \right. $$ (10) 对式(10)消去w和ε得：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{1}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} \mathit{\boldsymbol{\alpha }}\\ \mathit{\boldsymbol{b}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (11) 因wTw=αTKα，则式(11)迭代表达式为：
$$ \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} + {\rm{diag}}\left( {\frac{{1 + {{\left( {{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}} \right)}^{\left( i \right)}}K{\alpha ^{\left( i \right)}}}}{c}} \right)}&\mathit{\boldsymbol{P}}\\ {{\mathit{\boldsymbol{P}}^{\rm{T}}}}&\mathit{\boldsymbol{0}} \end{array}} \right]\left[ \begin{array}{l} {\mathit{\boldsymbol{\alpha }}^{\left( {i + 1} \right)}}\\ {\mathit{\boldsymbol{b}}^{\left( {i + 1} \right)}} \end{array} \right] = \left[ \begin{array}{l} \mathit{\boldsymbol{f}}\\ \mathit{\boldsymbol{0}} \end{array} \right] $$ (12) 式中，i为迭代次数。
比较式(5)和式(12)得出，MQ为MQ-T特殊情况，即当αTKα=0时的情况。因此，式(12)计算时，应先借助式(5)求算α，然后基于式(12)迭代计算，直至α和b收敛。
2 DEM建模算法精度分析
本文分别以数值实验和实例分析验证MQ-T计算精度。数值实验中，以高斯合成曲面为研究对象，分别设计受不同误差分量污染的采样数据，基于MQ-T曲面建模，并与传统MQ方法计算结果比较，验证模型的去噪性。实例分析中，以全站仪获取的采样数据为研究对象，借助MQ-T构建测区DEM，并将其与传统插值算法(如IDW，Kriging和ANUDEM)计算结果比较，验证模型的实用性。
2.1 数值实验
高斯合成曲面表达式为：
$$ \begin{array}{l} z\left( {x,y} \right) = 3{\left( {1 - x} \right)^2}{{\rm{e}}^{ - x2 - {{\left( {y + 1} \right)}^2}}} - \\ 10\left( {\frac{x}{5} - {x^3} - {y^5}} \right){{\rm{e}}^{ - {x^2} - {y^2}}} - \frac{1}{3}{{\rm{e}}^{ - {{\left( {x + 1} \right)}^2} - {y^2}}} \end{array} $$ 计算区域为[－3, 3]×[－3, 3]，模拟区域网格数为101 × 101。从高斯合成曲面中随机采集31×31个采样点用于曲面建模，采样点分别受如下误差影响。
例1：仅受高程误差分量影响。其中，误差来源于正态分布N (0, 0.02)。
例2：仅受水平误差分量影响。其中，x和y坐标误差均来源于N (0, 0.01)。因此，例1中高程误差分量方差和例2中水平误差分量方差相等。此方案可用于比较水平和高程误差对模拟结果精度影响的程度。
例3：同时受水平和高程误差分量影响。其中，水平和高程误差分别来源于正态分布N (0, 0.01)和N (0, 0.02)。
实例中分别采用中误差(root mean squared error，RMSE)和平均误差(mean error, ME)作为精度指标比较MQ-T和MQ模拟结果精度。中误差和平均误差计算的公式为：
$$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^k {{{\left( {{z_i} - {{\hat z}_i}} \right)}^2}} }}{k}} $$ (13) $$ {\rm{ME}} = \frac{{\sum\limits_{i = 1}^k {\left( {{z_i} - {{\hat z}_i}} \right)} }}{k} $$ (14) 式中，zi和${{\hat z}_i} $分别表示第i个检核点的真实值和模拟值；k表示检核点个数。
由表 1可见，当采样点仅受高程误差分量影响时(例1)，两种方法模拟结果相当，表明MQ-T对高程误差具有较好的光滑性。当采样点仅受水平误差影响时(例2)，两种方法模拟结果精度均低于例1，且MQ-T误差小于MQ误差，表明MQ-T能较好消除采样点水平误差分量对曲面模拟的影响。当采样点水平和高程坐标分量均受误差影响时(例3)，两种方法模拟结果精度均低于例1和例2，表明采样点精度越低，模拟结果精度越低；MQ-T中误差明显小于MQ，但前者平均误差略大于后者。综上所述，当采样点含有水平误差分量时，以中误差为误差指标，MQ-T模拟结果精度高于MQ。
表 1 3种误差分布条件下，MQ和MQ-T模拟结果精度比较Table 1. Accuracy Comparison Between the Classical MQ and MQ-T Under Different Error Distributions算例 误差分布 方法 RMSE ME 水平分量 高程分量 例1 - N (0, 0.02) MQMQ-T 0.062 70.062 1 0.0070.007 例2 N (0, 0.01) - MQMQ-T 0.147 80.139 3 -0.003-0.002 例3 N (0, 0.01) N (0, 0.02) MQMQ-T 0.158 70.142 7 -0.014-0.015 2.2 实例分析
测区位于山东省济南市，测区面积为3.36 km2，最小和最大高程分别为109.6 m和364.5 m。利用全站仪获取该测区5 394个采样点，平均点间距(sample interval, SI)为25 m。基于全站仪测角和测距精度指标以及数据采集时各种观测条件，采样点平面和高程精度分别为5 cm和10 cm。为了验证采样点密度对DEM精度的影响，将原始采样数据依次抽稀，每次抽稀后数据点数和对应样点间距如表 2所示。除了MQ-T和MQ外，还采用传统插值算法如IDW、Kriging以及ANUDEM进行10 m分辨率DEM建模。为了验证DEM精度，借助实时动态定位软件RTK获取了该测区500个检核点(见图 1)。各种插值算法构建DEM精度如表 3所示。
表 2 采样点与采样点平均点间距对应关系Table 2. Relationship Between Number of Sample Points and Sample Interva采样点数 点间距/m 2 100 40 2 742 35 3 721 30 5 394 25 表 3 各种插值算法模拟结果中误差/mTable 3. RMSE of Different Interpolation Methods/m点间距 IDW Kriging ANUDEM MQ MQ-T 40 3.93 2.43 2.68 2.38 2.26 35 3.41 2.37 2.55 2.24 2.11 30 3.17 2.22 2.43 2.09 2.06 25 2.69 1.84 2.28 1.73 1.73 平均 3.30 2.22 2.49 2.11 2.04 由表 3可见，随着平均点间距增大，各种插值算法精度均有所下降。具体而言，IDW、Kriging、ANUDEM和MQ降幅分别为46%、32%、15%和38%；MQ-T精度降幅为31%。结果表明，ANUDEM模拟结果受采样点间距影响最小，而IDW影响最大。不管采用何种采样点间距，MQ-T计算精度均高于其他插值方法。MQ-T平均计算精度分别是IDW的1.61倍、Kriging法的1.09倍、ANUDEM法的1.22倍以及MQ的1.03倍。
各种插值方法基于30 m平均点间距采样点构建的10 m分辨率DEM山体阴影图如图 2所示。由图 2可见，IDW结果最粗糙，这主要是因为IDW为准确插值算法，无法光滑去噪。ANUDEM法结果最光滑，去噪效果明显，但导致很多微地形特征信息丢失(如图 2(c)中圆圈标注的区域)。Kriging法、MQ法和MQ-T法结果相似，但前两种方法模拟结果有零散坑洼和凸起(如图 2(b)和2(d)中方框标注区域)。相比较而言，MQ-T法模拟结果稍好于其他方法。但相比Kriging法，MQ-T法结果左上角的山脊有峰值削平现象(如图 2(e)中圆圈标记区域)。
![]()
图 2 各种插值算法构建DEM的山体阴影图Figure 2. Hillshades of IDW, Kriging, ANUDEM, MQ and MQ-T Algorithms3 结语
为了降低采样数据中的水平和高程误差分量对DEM建模精度的影响，本文以较高插值精度的多面函数(MQ)为基函数，借助整体最小二乘算法(TLS)思路，发展了整体最小二乘MQ算法(MQ-T)。数值实验分析表明：(1)当采样点仅受高程误差分量影响时，MQ-T法计算结果精度与MQ法相当；(2)当采样点受水平误差分量影响，以中误差为误差指标，MQ-T计算结果精度高于MQ法。实例分析中，分别基于MQ法和MQ-T法构建测区DEM，并将计算结果与传统插值算法(如IDW、Kriging、ANUDEM)比较。结果分析表明：(1)各种插值算法结果精度均随采样点密度的降低而逐步降低，其中，ANUDEM精度变化最小，IDW精度变化最大；(2)不管采用何种采样点间隔，MQ-T精度始终高于其他插值算法。
-
![]()
图 2 各种插值算法构建DEM的山体阴影图
Figure 2. Hillshades of IDW, Kriging, ANUDEM, MQ and MQ-T Algorithms
表 1 3种误差分布条件下，MQ和MQ-T模拟结果精度比较
Table 1 Accuracy Comparison Between the Classical MQ and MQ-T Under Different Error Distributions
算例 误差分布 方法 RMSE ME 水平分量 高程分量 例1 - N (0, 0.02) MQMQ-T 0.062 70.062 1 0.0070.007 例2 N (0, 0.01) - MQMQ-T 0.147 80.139 3 -0.003-0.002 例3 N (0, 0.01) N (0, 0.02) MQMQ-T 0.158 70.142 7 -0.014-0.015 
下载: 导出CSV
表 2 采样点与采样点平均点间距对应关系
Table 2 Relationship Between Number of Sample Points and Sample Interva
采样点数 点间距/m 2 100 40 2 742 35 3 721 30 5 394 25
下载: 导出CSV
表 3 各种插值算法模拟结果中误差/m
Table 3 RMSE of Different Interpolation Methods/m
点间距 IDW Kriging ANUDEM MQ MQ-T 40 3.93 2.43 2.68 2.38 2.26 35 3.41 2.37 2.55 2.24 2.11 30 3.17 2.22 2.43 2.09 2.06 25 2.69 1.84 2.28 1.73 1.73 平均 3.30 2.22 2.49 2.11 2.04
下载: 导出CSV
-
[1] Chen C F, Yue T X. A Method of DEM Construction and Related Error Analysis[J]. Computers & Geosciences, 2010, 36(6):717-725 http://www.isprs.org/proceedings/XXXVII/congress/5_pdf/58.pdf
[2] Hutchinson M F, Gallant J C. Digital elevation models and representation of terrain shape[M]//Wilson J P, Gallant J C. Terrain Analysis: Principles and Applications. New York: Wiley, 2000
[3] Aguilar F J, Aguilar M A, Agüera F, et al. The Accuracy of Grid Digital Elevation Models Linearly Constructed from Scattered Sample Data[J]. International Journal of Geographical Information Science, 2006, 20(2):169-192 doi: 10.1080/13658810500399670
[4] Fisher P F, Tate N J. Causes and Consequences of Error in Digital Elevation Models[J]. Progress in Physical Geography, 2006, 30(4):467-489 doi: 10.1191/0309133306pp492ra
[5] Reinoso J F. A Priori Horizontal Displacement (HD) Estimation of Hydrological Features when Versioned DEMs are Used[J]. Journal of Hydrology, 2010, 384(1/2):130-141 https://www.sciencedirect.com/science/article/pii/S0022169410000405
[6] Vosselman G. Automated Planimetric Quality Control in High Accuracy Airborne Laser Scanning Surveys[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 74:90-100 doi: 10.1016/j.isprsjprs.2012.09.002
[7] Billings S D, Newsam G N, Beatson R K. Smooth Fitting of Geophysical Data Using Continuous Global Surfaces[J]. Geophysics, 2002, 67(6):1823-1834 doi: 10.1190/1.1527082
[8] Hutchinson M F, Gessler P E. Splines:More than Just a Smooth Interpolator[J]. Geoderma, 1994, 62(1-3):45-67 doi: 10.1016/0016-7061(94)90027-2
[9] 范家明, 陶大欣, 朱建青, 等.一种基于TIN的DEM表面C1光滑插值模型[J].测绘科学技术学报, 2007, 24(3):179-181 http://www.cqvip.com/qk/98086B/200703/24804748.html Fan Jiaming, Tao Daxin, Zhu Jianqing, et al. Surface Smooth Interpolation Model with C1 Continuities of DEM Based on TIN[J]. Journal of Zhengzhou Institute of Surverying and Mapping, 2007, 24(3):179-181 http://www.cqvip.com/qk/98086B/200703/24804748.html
[10] Grohmann C H, Steiner S S. SRTM Resample with Short Distance-Low Nugget Kriging[J]. International Journal of Geographical Information Science, 2008, 22(8):895-906 doi: 10.1080/13658810701730152
[11] Chen C F, Yue T X, Dai H L, et al. The Smoothness of HASM[J]. International Journal of Geographical Information Science, 2013, 27(8):1651-1667 doi: 10.1080/13658816.2013.787146
[12] Fan L, Smethurst J, Atkinson P, et al. Propagation of Vertical and Horizontal Source Data Errors into a TIN with Linear Interpolation[J]. International Journal of Geographical Information Science, 2014(1):1-23
[13] Reinoso J, León C, Mataix J. Estimating Horizontal Displacement between DEMs by Means of Particle Image Velocimetry Techniques[J]. Remote Sensing, 2016, 8(1):14-20
[14] Delacourt C, Allemand P, Casson B, et al. Velocity Field of the "La Clapière" Landslide Measured by the Correlation of Aerial and QuickBird Satellite Images[J]. Geophysical Research Letters, 2004, 31(15):1-6
[15] van Niel T G, McVicar T R, Li L, et al. The Impact of Misregistration on SRTM and DEM Image Differences[J]. Remote Sensing of Environment, 2008, 112(5):2430-2442 doi: 10.1016/j.rse.2007.11.003
[16] Golub G H, van Loan C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6):883-893 doi: 10.1137/0717073
[17] 姚宜斌, 黄书华, 孔建, 等.空间直线拟合的整体最小二乘算法[J].武汉大学学报·信息科学版, 2014, 39(5):571-574 http://ch.whu.edu.cn/CN/abstract/abstract2985.shtml Yao Yibin, Huang Shuhua, Kong Jian, et al. Total Least Algorithm for Straight Line Fitting[J]. Geomatics and Information Science of Wuhan University, 2014, 39(5):571-574 http://ch.whu.edu.cn/CN/abstract/abstract2985.shtml
[18] Michael K, Mathias A. A Weighted Total Least-Squares Algorithm for Fitting a Straight Line[J]. Measurement Science and Technology, 2007, 18(11):3438-3446 doi: 10.1088/0957-0233/18/11/025
[19] Franke R. Scattered Data Interpolation:Tests of Some Method[J]. Mathematics of Computation, 1982, 38(4):181-200 https://www.researchgate.net/publication/238834111_Scattered_Data_Interpolation_Tests_of_Some_Method
[20] Chen C F, Li Y Y, Yan C Q, et al. A Robust Algorithm of Multiquadric Method Based on an Improved Huber Loss Function for Interpolating Remote-Sensing-Derived Elevation Data Sets[J]. Remote Sensing, 2015, 7(3):3347-3371 doi: 10.3390/rs70303347
[21] Hardy R L. Multiquadric Equations of Topography and Other Irregular Surfaces[J]. Journal of Geophysical Research, 1971, 76(8):1905-1915 doi: 10.1029/JB076i008p01905
[22] Chen C F, Yan C Q, Cao X W, et al. A Greedy-Based Multiquadric Method for LiDAR-Derived Ground Data Reduction[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 102:110-121 doi: 10.1016/j.isprsjprs.2015.01.012
[23] Hoerl A E, Kannard R W, Baldwin K F. Ridge Regression:Some Simulations[J]. Communications in Statistics, 1975, 4(2):105-123 doi: 10.1080/03610927508827232
-
期刊类型引用(2)
1. 高原，朱娅男，陈传法，胡占占，胡保健. 高精度DEM建模的加权径向基函数插值方法. 武汉大学学报(信息科学版). 2023(08): 1373-1379 . 
百度学术
2. 李振洪，李鹏，丁咚，王厚杰. 全球高分辨率数字高程模型研究进展与展望. 武汉大学学报(信息科学版). 2018(12): 1927-1942 . 百度学术
其他类型引用(4)
计量
- 文章访问数: 1798
- HTML全文浏览量: 163
- PDF下载量: 266
- 被引次数: 6
