基于RANSAC的坐标系转换抗差算法研究

刘猛奎,赵明金,石波,王云鹏,张倩

(1.山东科技大学 测绘科学与工程学院,山东 青岛 266590;2. 山东省国土测绘院,山东 济南 250102)

0 引 言

在大地测量、工程测量、摄影测量等领域中,坐标系之间的转换是必不可少的,比如WGS-84坐标系、西安80坐标系、北京54坐标系以及地方坐标系之间的相互转换．但是由于测量过程中多方面因素的影响,获取的公共点中含有较多粗差点,势必会影响坐标转换参数的求算精度．因此国内外学者针对如何处理坐标转换中含有的粗差点问题进行了深入的研究．

目前,处理公共点粗差数据的方法有两种:一是将粗差归入函数模型；二是将粗差归入随机模型,后者又称抗差估计[1]．抗差估计最早由Huber提出,随后,Kubik[2]等人将抗差估计理论引入测绘界,提出了著名的“Danish法”;周江文[3]提出了等价权概念,将M估计最小二乘化,并给出了IGG1和IGG2两种有效的等价权函数;杨元喜等[4]应用相关等价权原理构造出了抗差估计解式并提出了IGG3方案．以上方法的思想均是基于选权迭代,其本质是选择不同的权函数,通过迭代更新,使含有粗差的观测值的权越来越小,从而实现改正粗差提高精度的目的[5]．在空间坐标转换中,普遍采用基于选权迭代思想的抗差估计算法．常见的等价权函数有Huber法、Tukey法、Danish法、IGG1方案、IGG3方案等,它们普遍应用在坐标系转换中以降低粗差点带来的影响．吴祖海等[6]结合假设检验与降权迭代两种方法来减弱粗差点的影响;徐波等[7]通过IGG1函数进行定权,降低了含粗差观测值的权值;张洋等[8]结合工程实例采用IGG3方案对误差过大的公共点重新定权;倪飞[9]使用了不同的抗差权函数(Tukey、IGG1、IGG3)验证其抵抗粗差的有效性;Ge等[10]针对GPS中的坐标转换问题,分别采用了IGG、IGG3等方案进行抗差估计;Gong等[11]提出了一种附有拉格朗日乘数的加权整体最小二乘方法,同时利用IGG权函数对公共点中混有的粗差点进行降权．上述几种方法的粗差点所占比例均在20%左右,允许的粗差点个数较少,当公共点中的粗差点个数不断增加时,会出现粗差点参与解算的现象,从而导致求解出来的坐标转换参数偏离真值．尤其当粗点个数差超过公共点个数一半时,传统的抗差估计算法将不再适用[12]．

本文运用随机抽样一致性(RANSAC)算法的思想对坐标系抗差算法展开研究分析．该方法最大的优点是当粗差点的个数超过50%时,仍然能够有效地处理粗差点[13]．RANSAC算法的思想不同于传统的滤波方法,在尽可能地使用少量的点预计数学模型参数基础上,利用余下的点来验证模型．对于一组包含异常数据的样本数据集,利用RANSAC算法,可以解算出数据的数学模型参数,并选取出有效的样本数据．RANSAC算法作为最有效的鲁棒估计算法之一,它在计算机视觉领域中应用广泛．目前国内外学者主要将其应用在基础矩阵[14]、特征匹配[15]、图像拼接[16]等方面,而在坐标系转换中粗差点问题解决方面应用较少．本文将RANSAC算法应用到坐标系转换中以解决粗差点个数较多的问题,探讨了一种利用RANSAC算法理论解算坐标转换参数的方法,同时与常用的抗差算法(基于IGG3方案的最小二乘抗差估计算法)进行比较,并利用仿真数据进行验证,算例结果表明,当公共点中含有大量的粗差点时,基于RANSAC的抗差算法进行抗差是可行的,且能克服传统抗差估计算法存在的局限性．

1 基于罗德里格矩阵的坐标系转换模型

为了同时适用于小角度和大角度的空间坐标转换,本文采用了基于罗德里格矩阵求解坐标转换参数的方法[17]．如图1所示,设有两个空间直角坐标系分别为O-XYZ和s-xyz．

图1 两个不同的空间直角坐标系

同一点M在不同坐标系下的坐标分别为(X,Y,Z)和(x,y,z)．它们之间的坐标转换关系为:

(1)

式中:λ为尺度比例因子,假设其初值是1;R为3×3的旋转矩阵．用罗德里格矩阵中的三个独立参数a,b,c表示旋转矩阵R．设反对称矩阵S为

(2)

式中,a,b,c相互独立．R可由S构成罗德里格矩阵

R=(I+S)(I-S)-1.

(3)

式中,I为三阶单位矩阵．

任意两组点分别按照式(1)构造方程组,做差得

(4)

联系式(2)、式(3),可得

(5)

式中:x21=x2-x1;y21=y2-y1;z21=z2-z1;X21=X2-X1;Y21=Y2-Y1;Z21=Z2-Z1．

该方程组中只有两个方程独立,要求解3个未知数需要构建另外一个方程组,整理可得

(6)

求解a,b,c,结合式(2)、(3)可以得到旋转矩阵的初值．将旋转矩阵带入式(1)可解得平移矢量初值．

由式(1)、式(2)、式(3)可得

(7)

这里的未知参数共有7个,即λ、a、b、c、ΔX、ΔY、ΔZ,将式(7)进行线性化处理得到误差方程:

v=AΔx+l.

(8)

故对式(8),利用最小二乘法可得:

Δx=(ATPA)-1(ATPl).

(9)

式中,

P为单位矩阵．

2 基于RANSAC的坐标转换抗差算法

本章首先介绍RANSAC算法的基本理论,然后详细叙述利用该算法剔除粗差点,最后对RANSAC算法与基于IGG3方案的最小二乘抗差算法进行了比较．

2.1 RANSAC算法

RANSAC即“随机抽样一致性”．它是从一组含有“异常点”的观测数据集合中,通过迭代方式预计模型参数．RANSAC算法的思想不同于传统的滤波方法,它是尽可能地使用少量的点预计数学模型参数,再利用余下的点来验证模型．RANSAC算法建立在样本集中既包括正常数据又包括异常数据的基础上．异常数据通常是由环境因素导致的测量误差或者不合理的假设造成的．需要明确的是RANSAC是一种随机算法,每次计算得出的结果可能会存在差异,但是总会存在一种合理的方案[18]．

RANSAC算法通过迭代的方式不断选取数据中的随机子集来满足模型参数求解的需要,被选择的子集一般被假设为局内点．通常采用如下方法:

1)存在一个用于解释观测数据的模型,模型的参数是由假设的局内点计算得到．

2)用上一步得到的参数化模型去检测剩余的数据,若存在某个数据点适用于此模型,则把它归为局内点．

3)若模型中包含足够多的局内点时,则认为此模型合理．

4)选择上述模型中的所有局内点,重新利用最小二乘算法解算出模型参数．

5)模型的评估,通常用模型的错误率表示．

对于上述过程,模型的选择依据有两个标准:舍弃因局内点太少而无法得出准确结果的模型以及选择比当前模型更加适合的模型．

RANSAC算法涉及到三个不确定的参数:1)模型的阈值t,用来剔除不适合模型的异常点．t可以视为对局内点噪声均方差的假设,对于不同的模型需要预设不同的阈值,该参数对RANSAC性能影响很大,它直接决定模型参数一致集的大小．2)算法的迭代次数k,与数据点的数量和粗点所占的比例有关．该参数直接影响剩余数据点参与模型参数的检验次数,从而影响算法的效率．迭代的次数决定着模型参数的精度．3)局内点的数量的阈值d,用来表明已经找到合适的模型．我们要根据特定的问题和数据集通过实验来确定参数t和d,而参数k可以根据式(10)进行推断[19]．

(10)

式中:p表示迭代过程中被选为局内点的概率;w表示在模型误差允许的范围内选取一个局内点的概率;q表示适用于模型的最少数据个数．

2.2 RANSAC算法剔除粗差点

结合2.1节内容,将RANSAC算法应用到坐标转换中．基本步骤如下:

2)利用上述解算的结果去测试剩余的公共点．在剩余的公共点中选择一个公共点,根据式(1)中求解的尺度比例因子、罗德里格参数和平移参数计算其在源坐标系下的坐标经转换后与目标坐标系下的三个方向的坐标差值．将三方向上坐标差值的平方和的平方根与设定的模型阈值t进行比较,若小于设定的阈值,则认为此点为局内点,存储下来;反之,继续在剩余的公共点中选取其它的公共点进行测试．直到剩余的公共点全部测试完成．

3)对步骤2)中的局内点个数与局内点数量阈值d进行比较,d的初始大小事先确定．若局内点个数大于d,则把d更新为当前的局内点个数,利用当前的局内点重新计算尺度比例因子、罗德里格参数和平移参数,并存储下来;反之,d的大小保持不变．

4)重复步骤1)～3),循环次数逐一增加,直到等于迭代次数k,结束循环．此时d所表示的局内点数量最多,输出其对应的尺度比例因子、罗德里格参数和平移参数．

具体算法的流程如图2所示．

其中,t,k,d表示的含义与2.1节的含义相同．

2.3 RANSAC抗差算法与基于IGG3方案的最小二乘抗差算法比较

在坐标转换中,本文提出的RANSAC抗差算法的核心在于使用尽可能少(3个)的公共点去计算坐标转换参数,之后利用剩余的公共点去验证直至找到最大的一致集．而基于IGG3方案的最小二乘抗差算法则是对误差较大的公共点进行重新定权,逐渐降低其在求解转换参数中的作用,具体步骤可参考文献[8]．本文第三章通过仿真数据分别验证了两种方案．

3 算例分析

图2 基于RANSAC算法的坐标转换流程图

为了便于研究多粗差点对两种抗差算法的影响,本文的实验数据来源于数据仿真,并针对大角度、小角度情形仿真两组实验数据,仿真数据的可靠性已得到证实[20]．其过程如下:

图3 公共点在源坐标系下和目标坐标系下的分布图

3)产生含有粗差点的公共点坐标．在20个公共点中每次随机选择m个点(m依次取1、2、…、13),将选择的公共点在目标坐标系下的三方向坐标与粗差值相加．其中,粗差值的产生符合正态分布N(0.5,0.12)(粗差值个数为39个)．

4)按照本文第2章的内容,分别使用最小二乘算法(方案一)、基于IGG3方案的最小二乘抗差估计算法(方案二)、基于RANSAC的抗差估计算法(方案三)进行坐标转换参数的求解．

不同粗差点下的坐标转换参数与真值的差值如表1所示,不同方案下的点位中误差大小如图4所示．

表1 不同粗差点个数下的坐标转换参数与真值的差值

表1(续) 不同粗差点个数下的坐标转换参数与真值的差值

图4 不同方案下的点位中误差大小

分析表1、图4可知:

1)当粗差点个数在1～6范围内,利用基于IGG3方案的最小二乘抗差估计算法与RANSAC算法求解的转换参数与真值的差值均不大,即三个平移参数差值均在10 mm左右,尺度比例因子差值在10-7量级,三个欧拉角参数差值均在10″以内．此外,两者的点位中误差大小也均在4 cm左右．

2)对于基于IGG3方案的最小二乘抗差估计算法而言,当粗差点个数为7时,求解的转换参数与真实值之差突然变大．而点位中误差也在含有7个粗差点的情况下骤增,精度大小由厘米级跳到分米级．而此后随着粗差点个数的增加,利用基于IGG3方案的最小二乘估计抗差算法得出的转换参数与真值的差值越来越大,且点位中误差的大小始终保持在分米级．

3)对于RANSAC算法而言,当粗差点个数在7～12时,其解算出来的转换参数与真实值之差仍保持在很小的范围内,点位中误差也基本保持在厘米级．即使当粗差点个数占公共点个数的50%,依旧能够进行坐标转换参数的求解．

为了更加直观地体现基于IGG3方案的最小二乘抗差估计算法和RANSAC算法的差异性,特选取6个粗差点和10个粗差点两种情况．不同情况下的20个点位在不同方案下的坐标较差分别如图5和图6所示．

由图5可知,利用基于IGG3方案的最小二乘抗差估计算法和RANSAC算法得出的坐标较差走势基本一致,且大小在厘米级附近．从而进一步证明了在粗差点较少的情况下,基于IGG3方案的最小二乘抗差估计算法和RANSAC算法均能很好地降低粗差点的影响．而由图6可知,方案二下20个点位的坐标较差,最大值在分米级上,甚至比不加抗差算法的结果更大,这就意味着基于IGG3方案的最小二乘抗差估计算法失效．反之,利用RANSAC算法所得到的X/Y/Z较差则保持在厘米级上,与含有6个粗差点时所得出的结果基本一致．

图5 不同方案下的X/Y/Z坐标较差(粗差点个数为6)

图6 不同方案下的X/Y/Z坐标较差(粗差点个数为10)

为了避免仿真数据具有偶然性,重新仿真一组小角度坐标转换的数据,仿真过程同上,参数变化如下:源坐标系下的公共点坐标取值范围均在-800～800 m,3个平移参数分别为250 m、300 m、500 m,3个欧拉角分别为0.015 rad、0.018 rad、0.02 rad,随机误差值服从N(0,0.0152),粗差值的产生符合正态分布N(0.4,0.12)．按照本文第2章的内容,分别使用最小二乘算法(方案一)、基于IGG3方案的最小二乘抗差估计算法(方案二)、基于RANSAC的抗差估计算法(方案三)进行坐标转换参数的求解,得到不同方案下的点位中误差大小,如图7所示．

图7 不同方案下的点位中误差大小

观察图7,此组数据下的点位中误差走势与图4基本相同,分析同上．

4 结束语

RANSAC算法常应用于影像自动匹配,而在坐标转换中也会存在类似问题,例如,地图矫正．故本文将RANSAC算法应用在坐标转换中以解决多粗差点的问题,提高坐标转换的精度．本文选用适用于大角度、小角度坐标转换的罗德里格矩阵模型求解七参数,同时比较基于IGG3方案的最小二乘抗差算法与RANSAC抗差算法．通过算例分析可得:

1)当公共点中的粗差点较少时,利用基于IGG3方案的最小二乘抗差估计算法和RANSAC算法均能得出高精度的坐标转换参数,且二者差异不大．

2)随着粗差点的增加,当达到一个临界值时,基于IGG3方案的最小二乘抗差估计算法失效,而利用RANSAC算法仍能够有效剔除粗差点,满足坐标转换的精度要求．

本文提出的基于RANSAC坐标转换抗差算法主要适用于粗差点个数多的情况,弥补了基于选权迭代思想的传统抗差算法的不足．即使当公共点中的粗差点个数占50%时,利用RANSAC算法也能有效剔除粗差点,克服粗差的影响．