劉猛奎,趙明金,石波,王云鵬,張倩

(1.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院,山東 青島 266590;2. 山東省國土測繪院,山東 濟(jì)南 250102)

0 引 言

在大地測量、工程測量、攝影測量等領(lǐng)域中,坐標(biāo)系之間的轉(zhuǎn)換是必不可少的,比如WGS-84坐標(biāo)系、西安80坐標(biāo)系、北京54坐標(biāo)系以及地方坐標(biāo)系之間的相互轉(zhuǎn)換．但是由于測量過程中多方面因素的影響,獲取的公共點(diǎn)中含有較多粗差點(diǎn),勢必會(huì)影響坐標(biāo)轉(zhuǎn)換參數(shù)的求算精度．因此國內(nèi)外學(xué)者針對如何處理坐標(biāo)轉(zhuǎn)換中含有的粗差點(diǎn)問題進(jìn)行了深入的研究．

目前,處理公共點(diǎn)粗差數(shù)據(jù)的方法有兩種:一是將粗差歸入函數(shù)模型；二是將粗差歸入隨機(jī)模型,后者又稱抗差估計(jì)[1]．抗差估計(jì)最早由Huber提出,隨后,Kubik[2]等人將抗差估計(jì)理論引入測繪界,提出了著名的“Danish法”;周江文[3]提出了等價(jià)權(quán)概念,將M估計(jì)最小二乘化,并給出了IGG1和IGG2兩種有效的等價(jià)權(quán)函數(shù);楊元喜等[4]應(yīng)用相關(guān)等價(jià)權(quán)原理構(gòu)造出了抗差估計(jì)解式并提出了IGG3方案．以上方法的思想均是基于選權(quán)迭代,其本質(zhì)是選擇不同的權(quán)函數(shù),通過迭代更新,使含有粗差的觀測值的權(quán)越來越小,從而實(shí)現(xiàn)改正粗差提高精度的目的[5]．在空間坐標(biāo)轉(zhuǎn)換中,普遍采用基于選權(quán)迭代思想的抗差估計(jì)算法．常見的等價(jià)權(quán)函數(shù)有Huber法、Tukey法、Danish法、IGG1方案、IGG3方案等,它們普遍應(yīng)用在坐標(biāo)系轉(zhuǎn)換中以降低粗差點(diǎn)帶來的影響．吳祖海等[6]結(jié)合假設(shè)檢驗(yàn)與降權(quán)迭代兩種方法來減弱粗差點(diǎn)的影響;徐波等[7]通過IGG1函數(shù)進(jìn)行定權(quán),降低了含粗差觀測值的權(quán)值;張洋等[8]結(jié)合工程實(shí)例采用IGG3方案對誤差過大的公共點(diǎn)重新定權(quán);倪飛[9]使用了不同的抗差權(quán)函數(shù)(Tukey、IGG1、IGG3)驗(yàn)證其抵抗粗差的有效性;Ge等[10]針對GPS中的坐標(biāo)轉(zhuǎn)換問題,分別采用了IGG、IGG3等方案進(jìn)行抗差估計(jì);Gong等[11]提出了一種附有拉格朗日乘數(shù)的加權(quán)整體最小二乘方法,同時(shí)利用IGG權(quán)函數(shù)對公共點(diǎn)中混有的粗差點(diǎn)進(jìn)行降權(quán)．上述幾種方法的粗差點(diǎn)所占比例均在20%左右,允許的粗差點(diǎn)個(gè)數(shù)較少,當(dāng)公共點(diǎn)中的粗差點(diǎn)個(gè)數(shù)不斷增加時(shí),會(huì)出現(xiàn)粗差點(diǎn)參與解算的現(xiàn)象,從而導(dǎo)致求解出來的坐標(biāo)轉(zhuǎn)換參數(shù)偏離真值．尤其當(dāng)粗點(diǎn)個(gè)數(shù)差超過公共點(diǎn)個(gè)數(shù)一半時(shí),傳統(tǒng)的抗差估計(jì)算法將不再適用[12]．

本文運(yùn)用隨機(jī)抽樣一致性(RANSAC)算法的思想對坐標(biāo)系抗差算法展開研究分析．該方法最大的優(yōu)點(diǎn)是當(dāng)粗差點(diǎn)的個(gè)數(shù)超過50%時(shí),仍然能夠有效地處理粗差點(diǎn)[13]．RANSAC算法的思想不同于傳統(tǒng)的濾波方法,在盡可能地使用少量的點(diǎn)預(yù)計(jì)數(shù)學(xué)模型參數(shù)基礎(chǔ)上,利用余下的點(diǎn)來驗(yàn)證模型．對于一組包含異常數(shù)據(jù)的樣本數(shù)據(jù)集,利用RANSAC算法,可以解算出數(shù)據(jù)的數(shù)學(xué)模型參數(shù),并選取出有效的樣本數(shù)據(jù)．RANSAC算法作為最有效的魯棒估計(jì)算法之一,它在計(jì)算機(jī)視覺領(lǐng)域中應(yīng)用廣泛．目前國內(nèi)外學(xué)者主要將其應(yīng)用在基礎(chǔ)矩陣[14]、特征匹配[15]、圖像拼接[16]等方面,而在坐標(biāo)系轉(zhuǎn)換中粗差點(diǎn)問題解決方面應(yīng)用較少．本文將RANSAC算法應(yīng)用到坐標(biāo)系轉(zhuǎn)換中以解決粗差點(diǎn)個(gè)數(shù)較多的問題,探討了一種利用RANSAC算法理論解算坐標(biāo)轉(zhuǎn)換參數(shù)的方法,同時(shí)與常用的抗差算法(基于IGG3方案的最小二乘抗差估計(jì)算法)進(jìn)行比較,并利用仿真數(shù)據(jù)進(jìn)行驗(yàn)證,算例結(jié)果表明,當(dāng)公共點(diǎn)中含有大量的粗差點(diǎn)時(shí),基于RANSAC的抗差算法進(jìn)行抗差是可行的,且能克服傳統(tǒng)抗差估計(jì)算法存在的局限性．

1 基于羅德里格矩陣的坐標(biāo)系轉(zhuǎn)換模型

為了同時(shí)適用于小角度和大角度的空間坐標(biāo)轉(zhuǎn)換,本文采用了基于羅德里格矩陣求解坐標(biāo)轉(zhuǎn)換參數(shù)的方法[17]．如圖1所示,設(shè)有兩個(gè)空間直角坐標(biāo)系分別為O-XYZ和s-xyz．

圖1 兩個(gè)不同的空間直角坐標(biāo)系

同一點(diǎn)M在不同坐標(biāo)系下的坐標(biāo)分別為(X,Y,Z)和(x,y,z)．它們之間的坐標(biāo)轉(zhuǎn)換關(guān)系為:

(1)

式中:λ為尺度比例因子,假設(shè)其初值是1;R為3×3的旋轉(zhuǎn)矩陣．用羅德里格矩陣中的三個(gè)獨(dú)立參數(shù)a,b,c表示旋轉(zhuǎn)矩陣R．設(shè)反對稱矩陣S為

(2)

式中,a,b,c相互獨(dú)立．R可由S構(gòu)成羅德里格矩陣

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

(3)

式中,I為三階單位矩陣．

任意兩組點(diǎn)分別按照式(1)構(gòu)造方程組,做差得

(4)

聯(lián)系式(2)、式(3),可得

(5)

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

該方程組中只有兩個(gè)方程獨(dú)立,要求解3個(gè)未知數(shù)需要構(gòu)建另外一個(gè)方程組,整理可得

(6)

求解a,b,c,結(jié)合式(2)、(3)可以得到旋轉(zhuǎn)矩陣的初值．將旋轉(zhuǎn)矩陣帶入式(1)可解得平移矢量初值．

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

(7)

這里的未知參數(shù)共有7個(gè),即λ、a、b、c、ΔX、ΔY、ΔZ,將式(7)進(jìn)行線性化處理得到誤差方程:

v=AΔx+l.

(8)

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

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

(9)

式中,

P為單位矩陣．

2 基于RANSAC的坐標(biāo)轉(zhuǎn)換抗差算法

本章首先介紹RANSAC算法的基本理論,然后詳細(xì)敘述利用該算法剔除粗差點(diǎn),最后對RANSAC算法與基于IGG3方案的最小二乘抗差算法進(jìn)行了比較．

2.1 RANSAC算法

RANSAC即“隨機(jī)抽樣一致性”．它是從一組含有“異常點(diǎn)”的觀測數(shù)據(jù)集合中,通過迭代方式預(yù)計(jì)模型參數(shù)．RANSAC算法的思想不同于傳統(tǒng)的濾波方法,它是盡可能地使用少量的點(diǎn)預(yù)計(jì)數(shù)學(xué)模型參數(shù),再利用余下的點(diǎn)來驗(yàn)證模型．RANSAC算法建立在樣本集中既包括正常數(shù)據(jù)又包括異常數(shù)據(jù)的基礎(chǔ)上．異常數(shù)據(jù)通常是由環(huán)境因素導(dǎo)致的測量誤差或者不合理的假設(shè)造成的．需要明確的是RANSAC是一種隨機(jī)算法,每次計(jì)算得出的結(jié)果可能會(huì)存在差異,但是總會(huì)存在一種合理的方案[18]．

RANSAC算法通過迭代的方式不斷選取數(shù)據(jù)中的隨機(jī)子集來滿足模型參數(shù)求解的需要,被選擇的子集一般被假設(shè)為局內(nèi)點(diǎn)．通常采用如下方法:

1)存在一個(gè)用于解釋觀測數(shù)據(jù)的模型,模型的參數(shù)是由假設(shè)的局內(nèi)點(diǎn)計(jì)算得到．

2)用上一步得到的參數(shù)化模型去檢測剩余的數(shù)據(jù),若存在某個(gè)數(shù)據(jù)點(diǎn)適用于此模型,則把它歸為局內(nèi)點(diǎn)．

3)若模型中包含足夠多的局內(nèi)點(diǎn)時(shí),則認(rèn)為此模型合理．

4)選擇上述模型中的所有局內(nèi)點(diǎn),重新利用最小二乘算法解算出模型參數(shù)．

5)模型的評估,通常用模型的錯(cuò)誤率表示．

對于上述過程,模型的選擇依據(jù)有兩個(gè)標(biāo)準(zhǔn):舍棄因局內(nèi)點(diǎn)太少而無法得出準(zhǔn)確結(jié)果的模型以及選擇比當(dāng)前模型更加適合的模型．

RANSAC算法涉及到三個(gè)不確定的參數(shù):1)模型的閾值t,用來剔除不適合模型的異常點(diǎn)．t可以視為對局內(nèi)點(diǎn)噪聲均方差的假設(shè),對于不同的模型需要預(yù)設(shè)不同的閾值,該參數(shù)對RANSAC性能影響很大,它直接決定模型參數(shù)一致集的大?。?)算法的迭代次數(shù)k,與數(shù)據(jù)點(diǎn)的數(shù)量和粗點(diǎn)所占的比例有關(guān)．該參數(shù)直接影響剩余數(shù)據(jù)點(diǎn)參與模型參數(shù)的檢驗(yàn)次數(shù),從而影響算法的效率．迭代的次數(shù)決定著模型參數(shù)的精度．3)局內(nèi)點(diǎn)的數(shù)量的閾值d,用來表明已經(jīng)找到合適的模型．我們要根據(jù)特定的問題和數(shù)據(jù)集通過實(shí)驗(yàn)來確定參數(shù)t和d,而參數(shù)k可以根據(jù)式(10)進(jìn)行推斷[19]．

(10)

式中:p表示迭代過程中被選為局內(nèi)點(diǎn)的概率;w表示在模型誤差允許的范圍內(nèi)選取一個(gè)局內(nèi)點(diǎn)的概率;q表示適用于模型的最少數(shù)據(jù)個(gè)數(shù)．

2.2 RANSAC算法剔除粗差點(diǎn)

結(jié)合2.1節(jié)內(nèi)容,將RANSAC算法應(yīng)用到坐標(biāo)轉(zhuǎn)換中．基本步驟如下:

2)利用上述解算的結(jié)果去測試剩余的公共點(diǎn)．在剩余的公共點(diǎn)中選擇一個(gè)公共點(diǎn),根據(jù)式(1)中求解的尺度比例因子、羅德里格參數(shù)和平移參數(shù)計(jì)算其在源坐標(biāo)系下的坐標(biāo)經(jīng)轉(zhuǎn)換后與目標(biāo)坐標(biāo)系下的三個(gè)方向的坐標(biāo)差值．將三方向上坐標(biāo)差值的平方和的平方根與設(shè)定的模型閾值t進(jìn)行比較,若小于設(shè)定的閾值,則認(rèn)為此點(diǎn)為局內(nèi)點(diǎn),存儲下來;反之,繼續(xù)在剩余的公共點(diǎn)中選取其它的公共點(diǎn)進(jìn)行測試．直到剩余的公共點(diǎn)全部測試完成．

3)對步驟2)中的局內(nèi)點(diǎn)個(gè)數(shù)與局內(nèi)點(diǎn)數(shù)量閾值d進(jìn)行比較,d的初始大小事先確定．若局內(nèi)點(diǎn)個(gè)數(shù)大于d,則把d更新為當(dāng)前的局內(nèi)點(diǎn)個(gè)數(shù),利用當(dāng)前的局內(nèi)點(diǎn)重新計(jì)算尺度比例因子、羅德里格參數(shù)和平移參數(shù),并存儲下來;反之,d的大小保持不變．

4)重復(fù)步驟1)～3),循環(huán)次數(shù)逐一增加,直到等于迭代次數(shù)k,結(jié)束循環(huán)．此時(shí)d所表示的局內(nèi)點(diǎn)數(shù)量最多,輸出其對應(yīng)的尺度比例因子、羅德里格參數(shù)和平移參數(shù)．

具體算法的流程如圖2所示．

其中,t,k,d表示的含義與2.1節(jié)的含義相同．

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

在坐標(biāo)轉(zhuǎn)換中,本文提出的RANSAC抗差算法的核心在于使用盡可能少(3個(gè))的公共點(diǎn)去計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù),之后利用剩余的公共點(diǎn)去驗(yàn)證直至找到最大的一致集．而基于IGG3方案的最小二乘抗差算法則是對誤差較大的公共點(diǎn)進(jìn)行重新定權(quán),逐漸降低其在求解轉(zhuǎn)換參數(shù)中的作用,具體步驟可參考文獻(xiàn)[8]．本文第三章通過仿真數(shù)據(jù)分別驗(yàn)證了兩種方案．

3 算例分析

圖2 基于RANSAC算法的坐標(biāo)轉(zhuǎn)換流程圖

為了便于研究多粗差點(diǎn)對兩種抗差算法的影響,本文的實(shí)驗(yàn)數(shù)據(jù)來源于數(shù)據(jù)仿真,并針對大角度、小角度情形仿真兩組實(shí)驗(yàn)數(shù)據(jù),仿真數(shù)據(jù)的可靠性已得到證實(shí)[20]．其過程如下:

圖3 公共點(diǎn)在源坐標(biāo)系下和目標(biāo)坐標(biāo)系下的分布圖

3)產(chǎn)生含有粗差點(diǎn)的公共點(diǎn)坐標(biāo)．在20個(gè)公共點(diǎn)中每次隨機(jī)選擇m個(gè)點(diǎn)(m依次取1、2、…、13),將選擇的公共點(diǎn)在目標(biāo)坐標(biāo)系下的三方向坐標(biāo)與粗差值相加．其中,粗差值的產(chǎn)生符合正態(tài)分布N(0.5,0.12)(粗差值個(gè)數(shù)為39個(gè))．

4)按照本文第2章的內(nèi)容,分別使用最小二乘算法(方案一)、基于IGG3方案的最小二乘抗差估計(jì)算法(方案二)、基于RANSAC的抗差估計(jì)算法(方案三)進(jìn)行坐標(biāo)轉(zhuǎn)換參數(shù)的求解．

不同粗差點(diǎn)下的坐標(biāo)轉(zhuǎn)換參數(shù)與真值的差值如表1所示,不同方案下的點(diǎn)位中誤差大小如圖4所示．

表1 不同粗差點(diǎn)個(gè)數(shù)下的坐標(biāo)轉(zhuǎn)換參數(shù)與真值的差值

表1(續(xù)) 不同粗差點(diǎn)個(gè)數(shù)下的坐標(biāo)轉(zhuǎn)換參數(shù)與真值的差值

圖4 不同方案下的點(diǎn)位中誤差大小

分析表1、圖4可知:

1)當(dāng)粗差點(diǎn)個(gè)數(shù)在1～6范圍內(nèi),利用基于IGG3方案的最小二乘抗差估計(jì)算法與RANSAC算法求解的轉(zhuǎn)換參數(shù)與真值的差值均不大,即三個(gè)平移參數(shù)差值均在10 mm左右,尺度比例因子差值在10-7量級,三個(gè)歐拉角參數(shù)差值均在10″以內(nèi)．此外,兩者的點(diǎn)位中誤差大小也均在4 cm左右．

2)對于基于IGG3方案的最小二乘抗差估計(jì)算法而言,當(dāng)粗差點(diǎn)個(gè)數(shù)為7時(shí),求解的轉(zhuǎn)換參數(shù)與真實(shí)值之差突然變大．而點(diǎn)位中誤差也在含有7個(gè)粗差點(diǎn)的情況下驟增,精度大小由厘米級跳到分米級．而此后隨著粗差點(diǎn)個(gè)數(shù)的增加,利用基于IGG3方案的最小二乘估計(jì)抗差算法得出的轉(zhuǎn)換參數(shù)與真值的差值越來越大,且點(diǎn)位中誤差的大小始終保持在分米級．

3)對于RANSAC算法而言,當(dāng)粗差點(diǎn)個(gè)數(shù)在7～12時(shí),其解算出來的轉(zhuǎn)換參數(shù)與真實(shí)值之差仍保持在很小的范圍內(nèi),點(diǎn)位中誤差也基本保持在厘米級．即使當(dāng)粗差點(diǎn)個(gè)數(shù)占公共點(diǎn)個(gè)數(shù)的50%,依舊能夠進(jìn)行坐標(biāo)轉(zhuǎn)換參數(shù)的求解．

為了更加直觀地體現(xiàn)基于IGG3方案的最小二乘抗差估計(jì)算法和RANSAC算法的差異性,特選取6個(gè)粗差點(diǎn)和10個(gè)粗差點(diǎn)兩種情況．不同情況下的20個(gè)點(diǎn)位在不同方案下的坐標(biāo)較差分別如圖5和圖6所示．

由圖5可知,利用基于IGG3方案的最小二乘抗差估計(jì)算法和RANSAC算法得出的坐標(biāo)較差走勢基本一致,且大小在厘米級附近．從而進(jìn)一步證明了在粗差點(diǎn)較少的情況下,基于IGG3方案的最小二乘抗差估計(jì)算法和RANSAC算法均能很好地降低粗差點(diǎn)的影響．而由圖6可知,方案二下20個(gè)點(diǎn)位的坐標(biāo)較差,最大值在分米級上,甚至比不加抗差算法的結(jié)果更大,這就意味著基于IGG3方案的最小二乘抗差估計(jì)算法失效．反之,利用RANSAC算法所得到的X/Y/Z較差則保持在厘米級上,與含有6個(gè)粗差點(diǎn)時(shí)所得出的結(jié)果基本一致．

圖5 不同方案下的X/Y/Z坐標(biāo)較差(粗差點(diǎn)個(gè)數(shù)為6)

圖6 不同方案下的X/Y/Z坐標(biāo)較差(粗差點(diǎn)個(gè)數(shù)為10)

為了避免仿真數(shù)據(jù)具有偶然性,重新仿真一組小角度坐標(biāo)轉(zhuǎn)換的數(shù)據(jù),仿真過程同上,參數(shù)變化如下:源坐標(biāo)系下的公共點(diǎn)坐標(biāo)取值范圍均在-800～800 m,3個(gè)平移參數(shù)分別為250 m、300 m、500 m,3個(gè)歐拉角分別為0.015 rad、0.018 rad、0.02 rad,隨機(jī)誤差值服從N(0,0.0152),粗差值的產(chǎn)生符合正態(tài)分布N(0.4,0.12)．按照本文第2章的內(nèi)容,分別使用最小二乘算法(方案一)、基于IGG3方案的最小二乘抗差估計(jì)算法(方案二)、基于RANSAC的抗差估計(jì)算法(方案三)進(jìn)行坐標(biāo)轉(zhuǎn)換參數(shù)的求解,得到不同方案下的點(diǎn)位中誤差大小,如圖7所示．

圖7 不同方案下的點(diǎn)位中誤差大小

觀察圖7,此組數(shù)據(jù)下的點(diǎn)位中誤差走勢與圖4基本相同,分析同上．

4 結(jié)束語

RANSAC算法常應(yīng)用于影像自動(dòng)匹配,而在坐標(biāo)轉(zhuǎn)換中也會(huì)存在類似問題,例如,地圖矯正．故本文將RANSAC算法應(yīng)用在坐標(biāo)轉(zhuǎn)換中以解決多粗差點(diǎn)的問題,提高坐標(biāo)轉(zhuǎn)換的精度．本文選用適用于大角度、小角度坐標(biāo)轉(zhuǎn)換的羅德里格矩陣模型求解七參數(shù),同時(shí)比較基于IGG3方案的最小二乘抗差算法與RANSAC抗差算法．通過算例分析可得:

1)當(dāng)公共點(diǎn)中的粗差點(diǎn)較少時(shí),利用基于IGG3方案的最小二乘抗差估計(jì)算法和RANSAC算法均能得出高精度的坐標(biāo)轉(zhuǎn)換參數(shù),且二者差異不大．

2)隨著粗差點(diǎn)的增加,當(dāng)達(dá)到一個(gè)臨界值時(shí),基于IGG3方案的最小二乘抗差估計(jì)算法失效,而利用RANSAC算法仍能夠有效剔除粗差點(diǎn),滿足坐標(biāo)轉(zhuǎn)換的精度要求．

本文提出的基于RANSAC坐標(biāo)轉(zhuǎn)換抗差算法主要適用于粗差點(diǎn)個(gè)數(shù)多的情況,彌補(bǔ)了基于選權(quán)迭代思想的傳統(tǒng)抗差算法的不足．即使當(dāng)公共點(diǎn)中的粗差點(diǎn)個(gè)數(shù)占50%時(shí),利用RANSAC算法也能有效剔除粗差點(diǎn),克服粗差的影響．