-
目前，遥感技术以其全天时、全天候、高分辨率等优势，已成为数字高程模型(digital elevation model，DEM)建模的主要数据源[1]。但在数据采集中，受仪器噪声、恶劣天气、多路径反射以及人为操作不当等因素影响，采样数据中不可避免含有异常值(粗差)[2-4]，进而严重影响DEM建模精度。
根据采样点对建模后DEM逼近的程度，传统插值算法可分为准确插值和光滑插值[5, 6]。其中，准确插值适用于采样点精度较高的情况，而当采样点含有偶然误差且精度较低时，光滑插值可以在移除偶然误差的同时，较好地保持采样数据中的信号特征[7, 8]。 地学领域中，常用的插值算法包括反距离加权(inverse distance weighting，IDW)、普通克里金(ordinary Kriging，OK)和专业DEM插值软件(Australian National University DEM，ANUDEM)等。IDW通过反距离函数模拟采样点的邻域,忽视了空间结构信息和领域以外的信息联系；OK是一种广义的线性回归方法,丢弃了非线性信息和非随机规律; ANUDEM以曲面连续光滑为前提假设，但该假设在实际地面中很难实现。理论上，这些方法均无法抑制粗差对DEM建模的影响。为了克服采样数据中异常值的影响，相关研究人员发展了一系列抗差插值算法。例如，张菊清等提出了加权平均抗差算法[9]。该算法认为，如果采样点与待求点高程值差值比较大，则认为该采样点含有粗差，进而降低其插值时权重。但当异常点与待求点距离非常近时，此方法失效。为了抑制多波束测深数据中的异常值对海底地形的影响，Debese等提出了基于M估计的多项式抗差插值算法[10]；陆丹等提出了基于截断最小二乘(least trimmed squares,LTS)多项式抗差算法[11]。这两种算法均以多项式表达式作为基函数，分别借助M估计和LTS抑制异常值对曲面拟合影响。但当海底地形较复杂时，由于多项式曲面拟合精度较低，该方法无法准确区分实际地形数据和异常值。鉴于此，Chen和Li提出了基于最小绝对偏差(least absolute deviation,LAD)的多面函数(multi-quadric,MQ)抗差算法(MQ-L)[12]。该算法首先借助空间填充算法选择部分采样点作为MQ中心节点，然后借助LAD解算MQ误差方程。与传统算法相比，该算法在一定程度上抑制了DEM构建时采样数据中的粗差影响，但MQ节点数的选择需要充分利用先验知识，致使该算法主观性较强，进而限制其推广使用。
为了克服上述抗差算法的缺陷，本文以MQ为基函数，由Huber损失函数和权重惩罚项组成目标函数，发展了MQ抗差插值算法(MQ-H)，并用数值试验和实例数据验证了本文算法的抗差性。
-
曲面建模旨在借助采样点z(xi,yi)(i=1,…,n)获得曲面函数f(x,y)的准确解析表达式，进而估计未知点函数值。曲面建模的一般模型可表达为：
(1) 式中，f(xi,yi)表示函数模拟值；r(xi,yi)表示采样或模拟误差。对MQ而言，f(x,y)解析表达式为：
(2) 式中，p(x,y)和b分别表示MQ多项式核函数及对应系数；di表示第i个采样点与待求点距离；q(d)和β分别为MQ径向核函数及其权重。龚健雅分析认为，三次曲面作为核函数精度优于Hardy所采用的双曲面法[13]，故本文采用q(d)=1+d3。
式(2)的矩阵表达式：
(3) 式中，f=[f1 … fn]T；b=[b1 … bm]T；β=[β1 … βn]T；
为了求算b和β，基于岭回归[14, 15]构建MQ目标函数：
(4) 式中，L(r)和βTQβ分别表示目标函数的损失函数和权重惩罚项；λ为光滑参数，其最优值可以借助交叉验证求得。当L(r)=r2/2时，式(4)即为传统MQ目标函数，b和β可由如下求得。
对F求关于β的导数得：
(5) (6) 对F求关于b的导数得：
(7) 将式(6)代入式(7)得：
(8) 由式(6)和式(8)得：
(9) 借助矩阵分块求逆公式，式(9)可转换为：
(10) 式中，M=Q+λI；N=PTM-1P。
-
由以上分析得出，传统的MQ目标函数因含有采样值残差的二次平方项，其建模过程不具有抗差性。相比而言，Huber损失函数[16]能够较好地抑制异常值对计算结果的影响，其表达式为：
(11) 式中，σ为残差标准差；t为二次平方项和线性项的界点(调制因子)。
当异常值远偏离常规采样点时，其对计算结果仍有较大影响，因此该部分异常点仍需剔除。式(11)表明，传统的Huber损失函数并不能完全抑制异常值对参数计算影响，本文对Huber损失函数进行了改进，其表达式为：
(12) 式(11)和(12)的比较结果表明，后者可以完全抑制残差较大异常点对计算结果的影响。本文取t1=2.5，t2=3。因此，以式(12)作为MQ-H目标函数中损失函数可较好抑制异常点对建模影响。MQ-H对b和β解算过程如下所示。
令
(13) 当ri∈Aj (j=1,2,3)时，分别对F求关于β和b的导数得：
(14) 式中,m=I2( Pb-f )-t1e1+ t1e3。令：
(15) 则Hessian矩阵可表达为：
(16) 所以：
(17) 借助牛顿迭代法，β和b可表达为：
(18) 式中，k为迭代次数。
借助矩阵分块求逆公式得：
(19) (20) 式中，
基于式(19)和(20)，式(18)可转换为：
(21) 比较式(10)和(21)可得，传统MQ是MQ-H的特殊情况，即I2是单位对角矩阵；e1=0，e3=0。由以上推导过程可见，MQ-H计算过程中，采样点即为中心节点。与MQ-L[12]相比，MQ-H避免了从采样点中随机选择中心节点的盲目性，进而克服了模拟结果的不确定性。
-
选择高斯合成曲面为研究对象，分析MQ-H的抗差性，并与传统MQ算法及MQ-L计算结果进行精度比较。含有误差的高斯合成曲面表达式为：
(22) 式中，r为误差项。本文试验中的误差主要来源于如下4个分布：标准正态分布N(0,1)，被污染的正态分布(1-θ)N(0,1)+θN(0,52)(θ为污染率，分别为10%、20%、30%)，拉普拉斯分布L(0,1)以及柯西分布C(0,1)。
设模拟曲面网格数为101×101，从模拟曲面中随机采集含有误差的51×51个点对曲面进行模拟，并将模拟结果与对应真实值比较。其中，任意一点真实函数值z(x,y)由(x,y)代入去除误差项r的式(22)计算而来。模拟结果精度指标分别为中误差、最大和最小误差。不同采样误差条件下，三种方法模拟结果如表 1所示。其中，MQ-H1表示用式(11)作为MQ-H损失函数，MQ-H2表示用式(12)作为MQ-H损失函数。
表 1 不同误差分布条件下，各种方法模拟结果比较
Table 1. Accuracy Comparison Under Different Error Distributions
误差分布 方法 最优λ 中误差 最大误差 最小误差 N(0,1) MQ-H2 0.600 9 0.216 2 0.702 6 -0.680 2 MQ-H1 0.888 7 0.211 6 0.675 1 -0.698 1 MQ 0.782 9 0.211 1 0.680 6 -0.675 4 MQ-L - 0.399 2 2.155 8 -1.372 9 (1-θ)N(0,1)+ θN(0,52) ,θ=0.1 MQ-H2 0.626 4 0.222 7 1.011 0 -1.006 6 MQ-H1 0.703 0 0.257 9 1.061 0 -0.957 5 MQ 3.271 2 0.461 0 1.961 2 -1.305 1 MQ-L - 0.449 4 2.675 7 -2.406 2 (1-θ)N(0,1)+ θN(0,52) ,θ=0.2 MQ-H2 0.536 9 0.254 1 0.764 3 -0.945 9 MQ-H1 1.781 9 0.382 8 1.824 8 -1.001 4 MQ 6.097 6 0.700 9 2.533 0 -1.743 4 MQ-L - 0.569 0 2.436 2 -3.563 3 (1-θ)N(0,1)+ θN(0,52) ,θ=0.3 MQ-H2 0.602 1 0.354 3 1.101 5 -1.377 3 MQ-H1 0.665 2 0.618 0 2.079 8 -2.570 2 MQ 0.741 2 1.078 8 3.384 4 -6.177 8 MQ-L - 0.685 8 3.087 8 -4.763 2 L(0,1) MQ-H2 0.602 1 0.220 5 1.252 4 -0.729 0 MQ-H1 0.665 2 0.223 2 1.514 4 -0.783 4 MQ 0.741 2 0.232 2 1.543 9 -0.809 1 MQ-L - 0.233 9 1.087 8 -0.763 2 C (0,1) MQ-H2 1.076 0 0.369 8 0.889 4 -1.018 7 MQ-H1 1.704 4 0.391 6 1.063 8 -1.495 2 MQ 5.125 6 1.459 1 5.034 1 -4.724 6 MQ-L - 0.609 9 3.123 6 -3.542 5 由表 1可见，当采样误差服从正态分布时，传统MQ以及抗差MQ(MQ-H1及MQ-H2)精度基本一致，而MQ-L精度低于其他三种方法。当采样误差服从拉普拉斯分布时，MQ-L精度同样低于其他三种方法，这是由于MQ-L核函数最优节点数无法准确确定，导致其精度有损失[12]。当采样误差服从正态分布时，MQ计算精度达到最优，因此其计算结果精度(中误差和最小误差)略高于抗差MQ。以上分析表明，抗差MQ有较高的计算精度，且改进Huber损失函数和未改进函数计算精度相当。
当采样值被粗差污染时，随着污染率的提高，各种方法的模拟精度均有所降低，但抗差MQ(MQ-H1及MQ-H2)精度均高于其他两种方法。当污染率从10%提高到30%时，传统MQ精度降低最明显，中误差从0.461 0增大到1.078 8；MQ-L精度降低次之，中误差从0.449 4增大到0.685 8。两个抗差MQ模型中，MQ-H1受异常值影响也非常明显，中误差从0.257 9增大到0.618 0；而MQ-H2抗差性最好，中误差从0.222 7增大到0.354 3。此结果表明，改进后的Huber损失函数可以有效地抑制异常值对MQ计算结果的影响。
当采样误差服从柯西分布时，传统MQ计算结果完全失真，其各种误差指标值均明显大于其他方法指标值。MQ-L计算结果精度高于传统MQ，但显著低于抗差MQ。这都说明改进Huber损失函数具有较好的抗差性，且MQ-H2精度略高于MQ-H1。
综上所述，相比传统MQ以及MQ-L，抗差MQ具有较高的计算精度，且改进后的Huber损失函数可明显提高抗差MQ计算精度。
-
试验测区位于山东省东部，测区面积为2.86 km2，平均高程39.3 m，高程标准差为12.4 m，测区地形如图 1所示。地面高程采样点数据来源于153幅无人遥测飞艇立体像对，点数为124 615个，平均点间距为4.8 m。分别基于传统MQ及抗差MQ构建测区5 m分辨率DEM，并将计算结果与IDW、OK及ANUDEM的精度比较。各种插值算法生成的DEM等高线及对应山体阴影如图 2~图 6所示。
图 1 测区正射影像
Figure 1. Orthogonal Image
图 2 IDW等高线及山体阴影
Figure 2. Contour and Hillshade of IDW
图 3 OK等高线及山体阴影
Figure 3. Contour and Hillshade of OK
图 4 ANUDEM等高线及山体阴影
Figure 4. Contour and Hillshade of ANUDEM
图 5 MQ等高线及山体阴影
Figure 5. Contour and Hillshade of MQ
图 6 MQ-H2等高线及山体阴影
Figure 6. Contour and Hillshade of MQ-H2
如图 2和图 3所示，在测区左上角以及右侧区域，IDW及OK等高线较为密集，表明采样数据中含有异常值，较多零散分布闭合等高线表明采样数据受偶然误差影响。由于IDW为准确插值，在插值过程中忠实反应原始数据精度、没有去除原始数据中的偶然误差以及粗差，因此其计算结果精度较低。考虑OK具有一定光滑性[17]，其山体阴影图比IDW略光滑。ANUDEM以及传统MQ算法均在一定程度上降低了采样数据中的偶然误差及异常值对计算结果影响，但它们是以损失地形细节信息为代价的，例如，原始地形中的山谷被明显填平，梯田状地形也随之消失(图 4和图 5)。相比较而言，抗差MQ在保持地形细节信息的同时，较好地抑制了采样点中异常值对计算结果的影响(图 6)，例如，零星分布的等高线大部分已消失，而原始地形中梯田以及山谷均得到了较好的保留。
-
为提高DEM建模精度，本文以较高精度的多面函数MQ为基函数，由Huber损失函数及权重惩罚项组成目标函数，发展了MQ抗差插值算法(MQ-H)。数值试验表明，当采样误差服从正态分布时，MQ-H计算精度与传统的MQ相当；当采样误差被粗差污染时，MQ-H能较好地抑制粗差对建模精度影响，其计算误差远小于传统MQ及MQ-L。实例分析表明，相比传统的插值算法(如IDW、OK及ANUMDE)，MQ-H不仅能抑制粗差及偶然误差对DEM建模的影响，而且可较好地捕捉到地形细节信息。
-
摘要: 为了抑制采样点中粗差对数字高程模型(digital elevation model,DEM)建模的影响,以较高精度的多面函数(multi-quadric,MQ)为基函数,由改进Huber损失函数和权重惩罚项组成目标函数,发展了MQ抗差插值算法(MQ-H)。通过优化MQ-H目标函数,采样点权重计算最终转换为方程组求解。以数学曲面为研究对象,将MQ-H计算结果与传统MQ及最小绝对偏差MQ(MQ-L)进行比较,结果表明:当采样误差服从正态分布时,MQ-H计算精度与传统MQ相当,而远高于MQ-L;当采样误差服从拉普拉斯分布时,MQ-H计算精度略高于MQ-L及传统MQ;当采样点被粗差污染时,MQ-H计算精度远高于传统MQ及MQ-L。在实例分析中,以无人遥测飞艇立体像对获取的地面离散高程点为基础数据,基于MQ-H构建测区DEM,并将计算结果与传统插值算法,如反距离加权(inverse distance weighting,IDW)、普通克里金(ordinary Kriging,OK)和专业DEM插值软件ANUDEM(Australian National University DEM)进行比较,结果表明,传统插值方法在不同程度上受采样点中异常值或偶然误差影响,而MQ-H受异常值影响较小,且能准确捕捉到地形细节信息。Abstract: In this paper, we propose a robust multi-quadric method (MQ-H) based on Huber loss function to conduct interpolations of contaminated spatial points, especially those derived from remote-sensing techniques. The objective function of the MQ-H has two main parts ; an improved Huber loss function and a regularized penalty term used to improve robustness and avoid overfitting, respectively. A mathematical surface, subject to model error with different distributions, was employed to comparatively analyze the robustness of the MQ-H, the classical MQ, and a least absolute deviation based MQ (MQ-L). The results indicated that when sample errors follow a normal distribution or a Laplacian distribution, the performance of MQ-H is comparatively better than those of MQ, and more accurate than MQ-L. For sample errors with a contaminated normal distribution and Cauchy distribution, MQ-H is more robust than MQ-L and MQ. Moreover, MQ with the improved Huber loss function is superior to MQ with the classical Huber loss function. A real-world example of DEM construction with stereo-image-derived elevation points indicates that compared to the classical interpolation methods including IDW (inverse distance weighting), OK (ordinary Kriging) and ANUDEM (Australian National University DEM), MQ-H has a better ability to reduce the impact of outliers while maintaining subtle terrain features suitable for qualitative analysis.
-
Key words:
- robust /
- multi-quadric /
- accuracy /
- DEM /
- Huber loss function
-
表 1 不同误差分布条件下，各种方法模拟结果比较
Table 1. Accuracy Comparison Under Different Error Distributions
误差分布 方法 最优λ 中误差 最大误差 最小误差 N(0,1) MQ-H2 0.600 9 0.216 2 0.702 6 -0.680 2 MQ-H1 0.888 7 0.211 6 0.675 1 -0.698 1 MQ 0.782 9 0.211 1 0.680 6 -0.675 4 MQ-L - 0.399 2 2.155 8 -1.372 9 (1-θ)N(0,1)+ θN(0,52) ,θ=0.1 MQ-H2 0.626 4 0.222 7 1.011 0 -1.006 6 MQ-H1 0.703 0 0.257 9 1.061 0 -0.957 5 MQ 3.271 2 0.461 0 1.961 2 -1.305 1 MQ-L - 0.449 4 2.675 7 -2.406 2 (1-θ)N(0,1)+ θN(0,52) ,θ=0.2 MQ-H2 0.536 9 0.254 1 0.764 3 -0.945 9 MQ-H1 1.781 9 0.382 8 1.824 8 -1.001 4 MQ 6.097 6 0.700 9 2.533 0 -1.743 4 MQ-L - 0.569 0 2.436 2 -3.563 3 (1-θ)N(0,1)+ θN(0,52) ,θ=0.3 MQ-H2 0.602 1 0.354 3 1.101 5 -1.377 3 MQ-H1 0.665 2 0.618 0 2.079 8 -2.570 2 MQ 0.741 2 1.078 8 3.384 4 -6.177 8 MQ-L - 0.685 8 3.087 8 -4.763 2 L(0,1) MQ-H2 0.602 1 0.220 5 1.252 4 -0.729 0 MQ-H1 0.665 2 0.223 2 1.514 4 -0.783 4 MQ 0.741 2 0.232 2 1.543 9 -0.809 1 MQ-L - 0.233 9 1.087 8 -0.763 2 C (0,1) MQ-H2 1.076 0 0.369 8 0.889 4 -1.018 7 MQ-H1 1.704 4 0.391 6 1.063 8 -1.495 2 MQ 5.125 6 1.459 1 5.034 1 -4.724 6 MQ-L - 0.609 9 3.123 6 -3.542 5 
下载: 导出CSV
-
[1] Robinson N, Regetz J, Guralnick R P. Earth Env-DEM90:A Nearly-Global, Void-Free, Multi-scale Smoothed, 90 m Digital Elevation Model from Fused ASTER and SRTM Data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 87(1):57-67 [2] Gonga-Saholiariliva N, Gunnell Y, Petit C, et al. Techniques for Quantifying the Accuracy of Gridded Elevation Models and for Mapping Uncertainty in Digital Terrain Analysis[J]. Progress in Physical Geography, 2011, 35(6):739-764 [3] Liu H, Jezek K C, O'Kelly M E. Detecting Outliers in Irregularly Distributed Spatial Data Sets by Locally Adaptive and Robust Statistical Analysis and GIS[J]. International Journal of Geographical Information Science, 2001, 15(8):721-741 [4] 陈传法,岳天祥.基于条件模拟的DEM误差曲面实现研究[J].武汉大学学报·信息科学版,2010,35(2):197-200 Chen Chuanfa, Yue Tianxiang. DEM Error Surface Relization Based on Conditional Simulation[J]. Geomatics and Information Science of Wuhan University, 2010, 35(2):197-200 [5] Billings S D, Beatson R K, Newsam G N. Interpolation of Geophysical Data Using Continuous Global Surfaces[J]. Geophysics, 2002, 67(6):1810-1822 [6] 张锦明,游雄,万刚.径向基函数算法中插值参数对DEM精度的影响[J].武汉大学学报·信息科学版,2013,38(5):608-611 Zhang Jinming, You Xiong, Wan Gang. Effects of Interpolation Parameters in Multi-log Radial Basis Function on DEM Accuracy[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5):608-611 [7] 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 [8] Hutchinson M F, Gessler P E. Splines-More than Just a Smooth Interpolator[J]. Geoderma, 1994, 62(1-3):45-67 [9] 张菊清,陈再辉,魏建忠.DEM空间数据抗差内插模型及其分析[J].测绘科学,2007,32(6):33-34+204 Zhang Juqing, Chen Zaihui, Wei Jianzhong. Robust Interpolation and Analyses for DEM[J]. Science of Surveying and Mapping, 2007, 32(6):33-34+204 [10] Debese N, Moitié R, Seube N. Multibeam Echosounder Data Cleaning through a Hierarchic Adaptive and Robust Local Surfacing[J]. Computers&Geosciences, 2012, 46:330-339 [11] 陆丹,李海森,魏玉阔,等. 基于截断最小二乘估计的多波束异常测深值剔除方法[J].大地测量与地球动力学, 2012, 32(1):89-93 Lu Dan, Li Haisen, Wei Yukuo, et al. A Method of Multi-beam Bathmetriy Outliers Eliminatin Based on Trimmed Least Square Estimation[J]. Journal of Geodesy and Geodynamics, 2012, 32(1):89-93 [12] Chen C F, Li Y Y. A Robust Multiquadric Method for Digital Elevation Model Construction[J]. Mathematical Geosciences, 2013, 45(3):297-319 [13] 龚健雅.关于DTM中多面函数内插法几个问题的研究[J]. 测绘通报,1985(5):12-16 Gong Jianya. Some Researches on Multiquadric Interpolation Problems[J]. Bulliten of Surveying and Mapping, 1985(5):12-16 [14] Hoerl A E, Kannard R W, Baldwin K F. Ridge Regression:Some Simulations[J]. Communications in Statistics, 1975, 4(2):105-123 [15] Kan B, Alpu Ö, Yazc B. Robust Ridge and Robust Liu Estimator for Regression Based on the LTS Estimator[J]. Journal of Applied Statistics, 2012, 40(3):644-655 [16] Huber P J. Robust Statistics[M]. New York:Wiley Blackwell, 1986 [17] Billings S D, Newsam G N, Beatson R K. SmoothFitting of Geophysical Data Using Continuous Global Surfaces[J]. Geophysics, 2002, 67(6):1823-1834 -
点击查看大图
图(6) / 表(1)
计量
- 文章访问数: 1350
- HTML全文浏览量: 38
- PDF下载量: 309
- 被引次数: 0
