王延存,俞家勇,田茂义,周茂伦,李国玉
(1.山东科技大学 测绘科学与工程学院,山东 青岛 266590;2.青岛秀山移动测量有限公司,山东 青岛 266590)
车载移动测量系统[1]由三维激光扫描仪、惯性测量单元(inertial measurement unit,IMU)、全球定位系统(global positioning system,GPS)等传感器组成,该系统能够快速获取高精度点云数据,可广泛应用于三维建模、大比例尺地形图、城市部件普查等方面,可以有效降低测量人员的劳动强度,提高工作效率。车载移动测量系统凭借自身的优势,成为当前研究的热点。然而点云数据的高程为大地高H,常规工程测量的高程是正常高h,为了能将车载激光点云数据运用于一般工程测量中,必须把大地高H转化为正常高h;大地高与正常高之间的差值为高程异常ζ,即ζ=H-h。高程异常如图1所示。因此如何获取高精度的高程异常改正值,成为点云数据处理的一个研究重点。目前确定高程异常的方法有地球重力场模型法[2-4]和GPS水准法。重力场模型法是通过EGM2008地球重力场模型计算待定点的高程异常值,该方法可以获取高精度高程异常值,但是重力数据获取比较困难且成本较高,因此在工程中应用比较少。GPS水准法是在测区内布设一定数量且分布均匀的控制点,进行GPS和水准联测,得到控制点的高程异常,进而通过曲面拟合获得该区域的高程异常数学表达式,再内插得到该区域其他点的高程异常值。常用的曲面拟合方法包括多项式曲面拟合法、多面函数法、移动曲面法、最小二乘配置法、薄板样条函数法(TPS)以及BP神经网络法等[5-15]。地形变化较大区域,文献[16-17]通过结合不同的拟合方法来适应不同的地区;针对小区域范围,文献[18-19]进行了研究,利用云模型对高程异常进行定性定量分析,形成云形图,直观获得该小区域的高程异常值。以上研究都取得了较好的效果。
图1 高程异常示意图
在实际的工程应用中,由于外界条件的影响,组合导航定位精度较差区域会导致点云高程误差大,此时利用文献[5-19]所提拟合方法求得的模型参数的可靠性降低,不能满足高精度测图要求。对于高程异常值含有粗差的情况,文献[20-24]进行了研究,但是局限于通过降低粗差的权值来提高模型参数的稳定性,而没有剔除粗差。且文献[20-24]都是针对大区域含粗差情况下的高程异常拟合,小区域范围内的稳健高程异常拟合相关研究较少。针对这一情况,本文提出了一种稳健高程异常改正方法,该方法以薄板样条函数为基础,结合RANSAC算法,并以狄克松判别法为准则,获取高程异常值。通过实验证明该方法可以有效的剔除粗差,获得高精度的高程异常改正值。
薄板样条函数插值[25-26]是一种二维的插值方法,在具有3个以上非共线点的情况下,薄板样条函数具有唯一解。假设在空间区域R2内有n个非共线点(t1,t2,…,tn),每个点的坐标为(xi,yi,zi)(i=1,2,…,n),其中zi=z(xi,yi)。
(1)
式中:δi和bj为待定系数;‖·‖为欧氏范数,在二维空间即为两点间的距离。ω是TPS的核函数,ω(r)=r2ln r/16π。当满足条件ωT=0时,式(1)就是空间区域R2上的自然薄板样条函数。
在小区域范围内,高程异常可以看成光滑连续曲面,非常适合通过自然薄板样条函数进行拟合。但当拟合数据中偶然误差较大或含有粗差时,直接利用自然薄板样条函数拟合得到的参数δi和bj不是最优解,进行插值时,距离粗差点越近的区域受到的影响越大。在这种情况下,需要引进一种稳健的拟合方法。本文在自然薄板样条函数的基础上,结合RANSAC算法,通过狄克松判别法寻找并剔除粗差点,进而得到该区域的高程异常值。
RANSAC算法是最有效的稳健估计算法之一,当数据中的粗差数据超过50%时,使用该算法仍然可以得到理想的估计参数。RANSAC算法首先设定一个准则,根据这个准则将数据分为局内点(可以用数学模型表达的数据)和局外点(不可以用数学模型表达的数据),最后利用保留下来的局内点进行参数计算,从而得到最优参数估值。最小迭代数k通过公式(2)计算得到。
(2)
式中:m为计算模型需要的最小数据量;ε为数据错误率。通常情况下,ε、P和k是根据实际情况动态确定的。
狄克松判别法采用极差比的思想,可以快速有效检测粗差的存在。RANSAC算法中,判定条件的选取直接决定最终结果的精度,因此在本文中,利用狄克松判别法作为局内点和局外点的判定条件。
通过RANSAC算法随机选取待拟合点,利用自然薄板样条函数进行拟合,得到每个待拟合点的高程异常残差Δξi,将残差Δξi按照从小到大的顺序排列。
Δξ1≤Δξ2≤…≤Δξn
(3)
(4)
利用自然薄板样条函数进行曲面拟合至少需要4个非共线点,根据RANCAC算法的思想,即需要在待拟合点中随机选取4个点作为种子点,根据式(2)计算最小迭代次数k,并计算自然薄板样条函数模型参数δi和bj的初始值,利用该初始值计算其余待拟合点的高程异常值,并与已知高程异常值做差得到高程异常值差值Δξi,利用狄克松判别法判定寻找该组拟合点对应的粗差值。将粗差值判定为局外点,其余的作为局内点,并统计局内点的数量。迭代统计最大局内点数量,将该组数据作为最终的待拟合点。该稳健拟合方法的流程如图2所示。
①利用ε、P、m根据式(2)计算出最小迭代次数k,本文中m=4。
②随机选取m个点,计算出模型参数δi和bj的初始值,并利用自然薄板样条函数的特性检核参数是否正确。
③根据②得到的参数计算Δξi,根据公式(3)、公式(4),进行自残序列两端逐点判别,如果小于临界值,将其判定为局内点,反之为局外点。
④重复②至③k次,统计每次迭代的局内点数量及局内点点号。
⑤选取最大局内点,根据公式(1)计算出参数δi和bj的最优解。
⑥通过最优解计算检核点残差。
图2 稳健拟合方法流程图
本文采用的数据为南方某地的车载移动测量系统获取的点云数据,该地区为丘陵地区。高程变化比较平缓。在测区内通过随机抽取,共选择了26个均匀分布的地面点,根据某勘测院1∶500水准精化模型获得对应的正常高,在该区域选择7个高程不同的地面点作为检核点,检核点高程分布在32~43 m范围内。高程异常拟合点点位分布如图3所示。
图3 待拟合点点位分布
针对移动测量作业过程中,组合导航系统定位精度受环境影响较大,采集前需要对作业区域进行规划,确保系统精度能够满足地籍测量要求。但在实际作业过程中,依然会存在局部小范围区域GPS定位精度较差情况,一般此区域高程误差在5~10 cm。一般误差较大区域占总的采集范围大约在5%~10%。在选择地面观测数据过程中,为保证拟合点均匀分布,采样拟合点难免会有少数部分在误差较大区域,需要在拟合方法中考虑粗差剔除,因此粗差取值范围在5~10 cm。本文考虑粗差点在拟合点数据中所占比例较小(考虑极端条件10%,拟合点数为26),故分别取1个、3个粗差点进行分析,阐述本文方法的稳健性。
目前最常用的高程异常拟合方法是多项式曲面法,本文采用三次曲面法进行拟合作为对比实验。三次曲面法是把高程异常值构成的曲面看作标准的三次曲面,而在地形变化较大的地区,高程异常值构成的曲面比较复杂,三次曲面法不能真实地描述该曲面,所得参数误差较大。因此该方法适用于高程异常变化比较平缓,且变化趋势近似于三次曲面的平原、低矮丘陵地区。使用该方法时,拟合点位要分布均匀且不含有粗差,否则局部地区可能会有较大误差。利用三次曲面法进行高程异常拟合,只需求解10个待定系数,计算效率较高,且方程中最高次项为三阶,有效避免了“龙格”现象。本文实验数据所在区域为丘陵地区,高程异常值变化比较平缓,适合使用三次曲面法进行拟合。
在实验中,分别利用三次曲面法、自然薄板样条函数法、稳健自然薄板样条函数法进行高程异常拟合,统计3种方法的检核点高程异常值的残差Δξi,结果如图4所示。
图4 3种方法各检核点残差比较
根据图4有以下分析结果:
①从图4(a)可以看出,在数据不包含粗差的情况下,其中三次曲面拟合的残差最大值为3.4 cm,而通过薄板样条函数进行拟合时,残差在2.5 cm以下,说明在该区域薄板样条函数法的拟合精度优于三次曲面法。实验表明,2种方法在拟合点不含有粗差的情况下,其精度都可以满足1∶500地形图要求。
②从图4(b)、图4(c)可以看出,在数据中通过数值模拟加入一个粗差点时,用薄板样条函数法拟合得到参数,求得的第3检核点的残差较大。加入3个粗差点时,薄板样条函数法对应的第3、4、6检核点的残差较大,最大值达到了5.2 cm,且只有一个检核点超过了1:500地形图的高程精度要求;利用三次曲面进行拟合得到的检核点的残差较大,最大值达到了6.8 cm,且在粗差较多的情况下,50%以上的检核点的高程精度都不能满足1∶500地形图的要求。因此,三次曲面法在该区域进行高程异常拟合的抗差性要低于薄板样条函数法。
③利用本文提出的方法,虽然随着待拟合点中粗差的增多,各检核点的残差增大,但是各检核点的残差都控制在3 cm以下,即在模拟的极端条件下(含有10%粗差),本文提出的方法在该区域依然满足1∶500地形图的精度要求。导致这一现象的原因是,传统的自然薄板样条函数法是将所有的数据不加区分地全部用于拟合中,这导致求解得到的模型参数不能反映真实的高程异常值,根据自然薄板样条函数的特性可知,距离粗差点越近的检核点,其残差值越大;而稳健薄板样条函数法对应的残差较小且稳定,是因为稳健法并使用RANSAC思想,并根据狄克松准则将粗差点剔除,寻找到最大局内点,最后利用保留下来的待拟合点(局内点)进行拟合,体现出该方法具有较强的稳健性。
在点云数据含有粗差的情况下,利用传统的拟合方法得到的模型参数是不可靠的。针对这一情况,本文在自然薄板样条函数的基础上,提出了一种稳健高程异常改正方法。该方法的核心思想是结合RANSAC算法,并根据狄克松判别法剔除粗差点,从而达到稳健的效果。实验结果表明,本文提出的方法可以在保证效率的前提下,有效地消除粗差点带来的影响,提高高程异常改正精度,解决了传统方法的不足。同时排除的粗差点还可以作为测图误差较大区域的标志,有助于车载激光点云数据质量把控,提高整体作图精度;该方法不仅可以应用于点云数据的高程异常拟合,还可以应用于其他的包含粗差的曲面拟合和制作DEM。需要指出的是,由于自然薄板样条函数拟合得到的曲面具有连续、光滑的特性,因此本方法在分区域拟合中也具有一定的优势。