M-LZD算法中表面变形的可探测性分析

张同刚，岑敏仪，任自珍，金国清

(1. 西南交通大学地球科学与环境工程学院，四川 成都 610031； 2. 西南交通大学高速铁路运营安全空间信息技术国家地方联合工程实验室，四川 成都 610031； 3. 中铁第五勘测设计院集团有限公司，北京 100000)



M-LZD算法中表面变形的可探测性分析

张同刚1，2，岑敏仪1，2，任自珍1，2，金国清3

(1. 西南交通大学地球科学与环境工程学院，四川 成都 610031； 2. 西南交通大学高速铁路运营安全空间信息技术国家地方联合工程实验室，四川 成都 610031； 3. 中铁第五勘测设计院集团有限公司，北京 100000)

摘要：现有的基于LZD方法的DEM变形探测方法都是将表面变形作为粗差，采用不同的稳健估计方法来探测。粗差具有发现和定位能力是粗差探测的前提，而现有方法均没有进行粗差发现和定位能力分析。本文首先采用判断矩阵对其进行了理论分析，然后对结果在不同类型地面的DEM进行了试验。根据试验结果，任意2个判断向量线性无关，判断矩阵中没有零元素，相关数等于多余观测个数。试验结果表明，DEM表面上无论多少粗差均具有发现和定位能力；基于LZD方法的DEM表面变形探测方法在理论上是完备的，可以正确地发现和定位粗差。

关键词：粗差发现；粗差定位；判断矩阵；相关数；LZD；DEM

无控制点DEM匹配差异探测是当前地学领域内的热点和难点问题。随着地面和机载激光扫描技术的发展，可以定期或根据需要对兴趣地点进行重复扫描，获取高分辨率的地面三维模型[1]。提取的多时相的DEM包含了丰富的地面变化信息，这是将其应用于滑坡等固体地球自然灾害监测、防治的关键问题[2-3]。不借助地面控制点完成不同时期DEM的配准来探测地面变形的技术为该问题提供了自动化的解决途径。

当前领域内最有代表性的DEM配准算法是最小高差算法(least Z-difference，LZD)，该方法采用最小二乘技术，能够在不借助控制点的情况下完成2期DEM的配准。在此方法的基础上，为了在有变形的情况下精确完成配准并确定地面变形的范围和幅度，Pilgrim利用稳健估计中的M估计来替代最小二乘，提出M-LZD算法，能够在不超过16%的情况下完成匹配[4]。更复杂的LMS估计和LTS估计也被引入到该领域，如LMS-LZD[5]和LTS-LZD[6]算法；为了进一步提升算法性能，表面局部不变特征[7]和变形量的空间关系[8]也被引入到该领域。随着算法性能的提升，算法也越来越复杂。

利用稳健估计替代最小二乘算法，在匹配过程

中将表面变形作为粗差来对待，是这类方法共同的基本思想。从这个角度来看，探测表面变形的过程实际上就是匹配过程中观测量粗差的探测过程。根据测量误差处理理论，观测量是否具备粗差发现和定位能力是任何粗差探测方法的前提[9-12]，对于不具备粗差定位能力的情况，任何方法准确定位其位置的概率不会超过50%[9]。这是关系到这一类基于稳健估计的LZD算法探测到的表面变形可靠性的基础性问题。文献[13]采用判断矩阵对DEM中任意一个粗差能否发现和定位问题进行了一些探索，但DEM变形是一个区域，与之对应的粗差肯定不止一个。

考虑到LMS估计和LTS估计算法模型非常复杂，本文采用判断矩阵中的相关数[14]，以这类算法中具有代表性的M-LZD算法为依托，对多个粗差发现和定位能力进行理论分析和模拟试验研究。

一、M-LZD算法模型

LZD算法采用刚体转换模型，采用七参数来描述表面间的转换关系[15]，即表面转换关系由3个旋转参数(θx,θy,θz)、3个平移参数(Δx，Δy，Δz)和缩放系数组成，其算法模型可以表示如下

V=AX+dzw

(1)

(2)

式中，A为系数阵，并且列满秩；V是改正数向量；X为待求转换参数；dz为Z坐标差；Fxi,Fyi表示第i点的X和Y方向的坡度，xi,yi表示i点的平面坐标；w为对应观测量的权，由M估计确定。

为方便描述起见，采用单变量来代替A中的元素，这样A可以简写成

(3)

二、判断矩阵和相关数

根据LZD算法的系数阵A，通过简单的矩阵行变换可以得到判断矩阵J[10]

(4)

(5)

判断矩阵中6个判断向量(J1-J6)可以求得6个粗差发现和定位相关数(以下简称“相关数”)(η1～η6)。根据文献[9]

1. 判断向量的线性相关性分析

判断矩阵可以表达成判断向量的形式

(6)

再令

(7)

将式(5)中各项代入J1—J7，并合并同类项系数，得

(8)

(9)

(10)

(11)

(12)

(13)

从式(8)—式(13)可以看出，判断向量J1—J6的形式非常类似，因此讨论6个判断向量两两之间的线性相关性只需要讨论任意2个判断向量之间的线性相关性即可。以J1、J2为例进行讨论。

假设

sJ1+tJ2=0(s、t为常数)

(14)

将式(8)、式(9)代入式(14)，并简化系数得

(15)

根据式(6)和式(7)将上式展开，得

(16)

两两作差，并代入各单变量的具体形式，得

(17)

式(17)成立的条件可以根据系数全部为0或不全部为0分成2种情况来分析。显然当系数全部为0时，判断向量J1、J2线性无关；其他情况下，J1、J2线性相关。下面着重讨论系数不全为0的情况，由于式(17)是和的形式，其成立的条件非常复杂，很难进行系统完善的分析。但显然式(18)是式(17)成立的条件中的一个。

(18)

根据式中符号的物理含义，Fx和Fy是DEM表面上X、Y方向上的坡度。式(18)中前2条比较简单，物理意义很明确，分别是DEM表面为X、Y方向的坡度分别相等，从DEM表面特征来看，也就是该表面在X、Y方向上是规则曲面。如果DEM确实为规则曲面，那么LZD数学模型中系数阵A(式(2))列秩亏，式(1)没有唯一解，这种情况与LZD方法的要求矛盾。式(18)中的其他3个条件的物理含义不太明确，不便于分析，但是前2个条件不成立，式(18)肯定不成立。

对于式(17)中在系数不全为0时，其他成立条件很难通过严格的数学公式来进行分析。结合实际DEM特点对此进行讨论。式(17)中包括DEM的平面坐标和表面坡度。对于实际的规则格网DEM而言，一方面DEM表面坡度与地形相关，没有固定的变化规律；另一方面，平面坐标X、Y具有很强的规律。再考虑到式(17)中各个条件相互独立，这样式(17)在系数不全为零时通常不成立。

对于实际DEM而言，其包含的数据量一般非常大，这样式(17)中包括的条件就非常多，由于这些条件相互独立，因此要使得它们同时成立的情况非常罕见。因此，对于实际DEM匹配而言，其判断向量J1、J2线性无关。类似的，可以证明J2和J3，J3和J4，J4和J5，J5和J6，J1和J6也线性无关。

2. 粗差发现和定位相关数分析

根据式(8)—式(13)，判断矩阵的6个判断向量(J1-J6)形式一致，对于其对应的6个相关数(η1-η6)分析可以任取其中的一个来进行分析。

(18)

三、试验

为了对以上理论进行验证，并兼顾试验的代表性，选择各种类型地形的DEM来进行模拟试验(如图1所示)。这些DEM是实际地形数据，为规则格网数据。格网间距为10m，大小为120×100像素，在试验中作为基准DEM。待匹配DEM通过对基准DEM旋转和平移后，并添加随机误差和变形得到。根据M-LZD算法的探测能力，添加变形为15%。试验中添加的随机误差和粗差均服从正态分布N(0,σ)。

图1　试验模型

试验中在DEM上均匀选取6个点作为必要观测量，分别进行计算相应的判断矩阵和相关数。限于篇幅，表1仅给出了其中一个判断矩阵及其相关数。所有的试验结果中，任意2个判断向量均线性无关，这表明DEM表面存在任意一个粗差都具备发现和定位能力；在判断矩阵中没有出现零元素的情况，相关数均等于多余观测量的个数，这表明最多可以发现粗差的个数为多余观测数的一半，可以定位的粗差比可以发现粗差个数少一个。

表1　判断矩阵和相关数

根据上面的试验结果，在实际应用中只要DEM表面变形的比例不超过50%，均可以发现和定位。本文6组数据的测试结果也表明，添加的15%的变形均被准确定位，因此M-LZD算法中DEM表面变形具备可探测性，算法确定的变形是可靠的。

四、结论

采用判断矩阵对基于LZD算法的DEM表面粗差发现和定位能力进行理论分析，并进行了试验验证，结果表明DEM表面具有粗差发现和定位能力，因此基于LZD方法将表面变形作为粗差来处理的方法在理论上是完备的，可以正确地发现和定位粗差。

采用各种类型的DEM进行了一系列的模拟试验，不同DEM上的试验结果一致，并与理论分析相符。试验结果表明任意2个判断向量之间均是线性无关的；判断矩阵中没有出现零元素，相关数与多余观测数相等。

参考文献：

[1]武永斌, 卢小平, 陈曦东, 等. 机载LiDAR铁路测绘关键技术及应用[J]. 测绘通报, 2015(9): 64-67.

[2]张同刚, 岑敏仪, 冯义从, 等. 采用截尾最小二乘估计的DEM匹配方法[J]. 测绘学报, 2009, 38(2): 144-151.

[3]汪燕麟, 殷义程, 施昆. 地震灾区中地面三维激光扫描测绘技术的应用方法分析[J]. 测绘通报, 2015(6): 65-68.

[4]PILGRIMLJ.RobustEstimationAppliedtoSurfaceMatching[J].ISPRSJournalofPhotogrammetryandRemoteSensing, 1996(51)：243-257.

[5]LIZ,XUZ,CENM,etal.RobustSurfaceMatchingforAutomatedDetectionofLocalDeformationsUsingLeast-Median-of-SquaresEstimator[J].PhotogrammetricEngineering&RemoteSensing, 2001, 67(11): 1283-1292.

[6]ZHANGT,CENM.RobustDEMCo-registrationMethodforTerrainChangesAssessmentUsingLeastTrimmedSquaresEstimator[J].AdvancesinSpaceResearch, 2008, 41(11): 1827-1835.

[7]SHARPGC,LEESW,WEHEDK.ICPRegistrationUsingInvariantFeatures[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 2002, 24(1): 90-102.

[8]LIUY,WEIB.DevelopingStructuralConstraintsforAccurateRegistrationofOverlappingRangeImages[J].RoboticsandAutonomousSystems, 2004, 47(1): 11-30.

[9]李德仁, 袁修孝. 误差处理与可靠性理论[M]. 武汉: 武汉大学出版社, 2002: 235-292.

[10]岑敏仪, 卓健成, 李志林, 等. 判断观测值粗差能否发现和定位的一种验前方法[J]. 测绘学报, 2003, 32(2): 134-138.

[11]岑敏仪, 顾利亚, 李志林, 等. 基于判断矩阵的观测量粗差发现和定位相关性分析[J]. 测绘学报, 2005, 34(1): 24-29.

[12]戴激光, 王新亮, 郭一洋. 一种稳健的倾斜航空影像粗差探测方法[J]. 测绘通报, 2014(5): 28-31.

[13]ZHANG T, CEN M, REN Z, et al. Ability to Detect and Locate Gross Errors on DEM Matching Algorithm[J]. International Journal of Digital Earth, 2010, 3(1): 72-82.

[14]顾利亚, 岑敏仪, 李志林. 粗差发现和定位能力与相关系数的关系[J]. 武汉大学学报(信息科学版)，2005, 30(7): 621-624.

[15]ROSENHOLM D, TORELEGÅRD K. Three-dimensional Absolute Orientation of Stereo Models Using Digital Elevation Models[J]. Photogrammetric Engineering & Remote Sensing, 1988, 54(10): 1385-1389.

Analysis of Surface Deformation Detectability Based on M-LZD Algorithm

ZHANG Tonggang，CEN Minyi，REN Zizhen，JIN Guoqing

收稿日期：2015-06-23； 修回日期： 2016-01-04

基金项目：中国铁路总公司重大科研项目(2014G004-B)；中国铁路总公司科研试验任务(Z2013-038-5)；长江学者和创新团队发展计划(IRT13092)

作者简介：张同刚(1977—)，男，博士，副教授，研究方向为数字摄影测量与模式识别。E-mail：swjtuztg@gmail.com

中图分类号：P237

文献标识码：B

文章编号：0494-0911(2016)06-0010-04

引文格式: 张同刚，岑敏仪，任自珍，等． M-LZD 算法中表面变形的可探测性分析[J]． 测绘通报，2016 ( 6) : 10-13． DOI: 10． 13474 /j． cnki． 11-2246． 2016． 0179．