基于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算法也能有效剔除粗差點,克服粗差的影響．