基于羅德里格矩陣的抗差迭代坐標(biāo)轉(zhuǎn)換算法研究

李國(guó)琴，田林亞，郭英起，張 洋，畢繼鑫

(1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098；2.黑龍江工程學(xué)院 測(cè)繪工程學(xué)院，黑龍江 哈爾濱 150050；3.江蘇蘇州地質(zhì)工程勘察院，江蘇 蘇州 215129)

對(duì)于不同基準(zhǔn)下的空間直角坐標(biāo)的轉(zhuǎn)換，Bursa-Wolf模型、Molodensky模型和武測(cè)模型已得到廣泛應(yīng)用[1]，但是這三種模型主要適用于旋轉(zhuǎn)角為小角度的特殊前提。在工業(yè)安裝測(cè)量、盾構(gòu)機(jī)導(dǎo)向測(cè)量、海底沉管測(cè)量等許多實(shí)際測(cè)量中，坐標(biāo)轉(zhuǎn)換時(shí)經(jīng)常遇到大旋轉(zhuǎn)角的情況，如果仍簡(jiǎn)單套用已有的坐標(biāo)轉(zhuǎn)換模型和算法，就可能得到錯(cuò)誤的計(jì)算結(jié)果和結(jié)論。近些年，一些學(xué)者對(duì)大旋轉(zhuǎn)角情況下的坐標(biāo)轉(zhuǎn)換問(wèn)題做了較多的研究。同濟(jì)大學(xué)陳義教授提出一種顧及旋轉(zhuǎn)矩陣元素之間的相關(guān)性，利用附有限制條件的間接平差求取13個(gè)參數(shù)進(jìn)行坐標(biāo)轉(zhuǎn)換的方法[2]；武漢大學(xué)曾文憲教授提出三維坐標(biāo)轉(zhuǎn)換的非線性模型，實(shí)現(xiàn)了旋轉(zhuǎn)角在50°以內(nèi)的坐標(biāo)轉(zhuǎn)換[3]；文獻(xiàn)[4]～文獻(xiàn)[5]研究了利用單位四元數(shù)法構(gòu)造旋轉(zhuǎn)矩陣進(jìn)行坐標(biāo)轉(zhuǎn)換的方法；文獻(xiàn)[6]～文獻(xiàn)[7]將羅德里格矩陣引入到坐標(biāo)轉(zhuǎn)換中，無(wú)需進(jìn)行繁瑣的三角計(jì)算，但并未顧及公共點(diǎn)自身精度對(duì)坐標(biāo)轉(zhuǎn)換結(jié)果的影響；文獻(xiàn)[8]將穩(wěn)健抗差估計(jì)理論應(yīng)用到坐標(biāo)轉(zhuǎn)換中，但該方法只適用于旋轉(zhuǎn)角為小角度的坐標(biāo)轉(zhuǎn)換。本文提出一種基于羅德里格矩陣的抗差迭代算法，其主要思想：首先計(jì)算基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換參數(shù)的初值，然后通過(guò)線性化模型計(jì)算轉(zhuǎn)換參數(shù)的改正數(shù)，同時(shí)應(yīng)用穩(wěn)健抗差估計(jì)理論中的選權(quán)迭代法進(jìn)行迭代計(jì)算，剔除含有粗差的公共點(diǎn)以及減小誤差過(guò)大的點(diǎn)對(duì)轉(zhuǎn)換結(jié)果的影響。該算法不僅同時(shí)適用小旋轉(zhuǎn)角和大旋轉(zhuǎn)角情況的坐標(biāo)轉(zhuǎn)換，而且可以有效抵抗可能存在的數(shù)據(jù)粗差對(duì)坐標(biāo)轉(zhuǎn)換結(jié)果的影響，計(jì)算模型簡(jiǎn)單，便于程序?qū)崿F(xiàn)，收斂速度快，計(jì)算精度高，可在實(shí)際工作中推廣應(yīng)用。

1 基于羅德里格矩陣的抗差迭代坐標(biāo)轉(zhuǎn)換算法

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

文獻(xiàn)[9]敘述了空間直角坐標(biāo)轉(zhuǎn)換模型，將羅德里格矩陣R代入到該轉(zhuǎn)換模型中，可得到含有3個(gè)平移參數(shù)、3個(gè)反對(duì)稱矩陣參數(shù)、1個(gè)尺度比參數(shù)的羅德里格矩陣坐標(biāo)轉(zhuǎn)換七參數(shù)模型，即

(1)

式中：下標(biāo)q表示欲轉(zhuǎn)換的目的坐標(biāo)系，下標(biāo)p表示原坐標(biāo)系；λ為尺度比參數(shù)；ΔX，ΔY，ΔZ為三個(gè)平移參數(shù)；R為羅德里格矩陣，滿足R=(I+S)(I-S)-1=(I-S)-1(I+S)，I為三階單位陣，反對(duì)稱矩陣S是由三個(gè)獨(dú)立的未知數(shù)參數(shù)a，b，c組成，即

(2)

先對(duì)基于羅德里格矩陣的坐標(biāo)轉(zhuǎn)換模型進(jìn)行線性化，即

(3)

其中，λ0為尺度比參數(shù)初值，R0為羅德里格矩陣初值，ΔX0，ΔY0，ΔZ0為三個(gè)平移參數(shù)初值，得到誤差方程

(4)

可解得轉(zhuǎn)換參數(shù)改正數(shù)

(5)

其中，

B11=B22=B33=1,

B12=B13=B21=B23=B31=B32=0,

B14=κ[(4ab2+4ac2)Xpi+(2a2b+4ac-2b-

2b3-2bc2)Ypi+(2c+2b2c+2c3-2a2c+4ab)Zpi],

B15=κ[(-4b-4a2b)Xpi+(2ab2+4bc-2a-

2a3-2ac2)Ypi+(2b2-2a2-2c2-4abc-2)Zpi],

B16=κ[(-4c-4a2c)Xpi+(4abc+2c2-2-

2a2-2b2)Ypi+(2a3+2a+2ab2-2ac2+4bc)Zpi],

B24=κ[(2a2b-2b-2b3-2bc2-4ac)Xpi-

(4a+4ab2)Ypi+(2a2-2b2-2c2+4abc-2)Zpi],

B25=κ[(2ab2-2ac2-4bc-2a3-2a)Xpi+

(4a2b+4bc2)Ypi+(2b2c-2a2c+4ab-2c3-2c)Zpi],

B26=κ[(2a2+2b2-2c2+4abc+2)Xpi+

(4b2c-4c)Ypi+(2bc2+4ac-2a2b-2b3-2b)Zpi],

B34=κ[(2b2c+2c3-2a2c-4ab+2c)Xpi+

(2b2+2c2-2a2+4abc+2)Ypi+(-4a-4ac2)Zpi],

B35=κ[(2a2-2b2+2c2-4abc+2)Xpi+

(2b2c-2a2c-4ab-2c3-2c)Ypi+(-4b-4bc2)Zpi],

B36=κ[(2ab2-2ac2-4bc+2a3+2a)Xpi+

(2bc2-2a2b-2b3-4ac-2b)Ypi+(4a2c+4b2c)Zpi],

B17=κ[(1+a2-b2-c2)Xpi-(2ab+2c)Ypi+

(-2b+2ac)Zpi],

B27=κ[(2c-2ab)Xpi+(1-a2+b2-

c2)Ypi-(2a+2bc)Zpi],

B37=κ[(2ac+2b)Xpi+(2a-2bc)Ypi+

(1-a2-b2+c2)Zpi].

為了實(shí)現(xiàn)基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換，需要模型參數(shù)初值計(jì)算。首先，計(jì)算尺度比參數(shù)λ的初值?？梢酝ㄟ^(guò)位于兩套坐標(biāo)系的兩個(gè)公共點(diǎn)的坐標(biāo)進(jìn)行計(jì)算，即

(6)

當(dāng)含有多個(gè)公共點(diǎn)時(shí)，可以分別求取各公共點(diǎn)對(duì)應(yīng)的邊長(zhǎng)之比，然后取其平均值計(jì)算尺度比參數(shù)初值。其次，求取3個(gè)反對(duì)稱矩陣參數(shù)的初值，將位于兩套坐標(biāo)系下的公共點(diǎn)1和公共點(diǎn)2的坐標(biāo)分別代入式(1)，做差可消去三個(gè)平移參數(shù)，得

(7)

其中，

Xq21=Xq2-Xq1，Yq21=Yq2-Yq1，

Zq21=Zq2-Zq1,

Xp21=Xp2-Xp1，Yp21=Yp2-Yp1，

Zp21=Zp2-Zp1.

根據(jù)反對(duì)稱矩陣S和羅德里格矩陣R的關(guān)系，可得

(8)

將(I+S)，(I-S)代入，可得

(9)

式(9)中的系數(shù)陣為奇異陣，無(wú)法解出參數(shù)a，b，c初值的唯一解。因此，再將位于兩套坐標(biāo)系下的公共點(diǎn)3代入式(1)，并與公共點(diǎn)1代入式(1)的結(jié)果做差，聯(lián)立式(9)得

(10)

對(duì)式(10)，采用最小二乘法得到反對(duì)稱矩陣參數(shù)a，b，c的初值，進(jìn)而得到羅德里格矩陣的初值R0。最后，將位于兩套坐標(biāo)系下的公共點(diǎn)1的坐標(biāo)代入式(1)，可求出三個(gè)平移參數(shù)ΔX，ΔY，ΔZ的初值(也可以利用多個(gè)公共點(diǎn)通過(guò)式(8)計(jì)算三個(gè)平移參數(shù)，然后取其平均值)，即

(11)

1.2 基于穩(wěn)健抗差估計(jì)理論的選權(quán)迭代

根據(jù)現(xiàn)代平差理論，經(jīng)典最小二乘法不具有抗干擾性和抵抗粗差的能力[10]。如果所選擇使用的公共點(diǎn)中某個(gè)點(diǎn)精度較低或誤差過(guò)大，勢(shì)必會(huì)影響轉(zhuǎn)換參數(shù)的解算精度，進(jìn)而會(huì)影響坐標(biāo)轉(zhuǎn)換的精度[8]。因此，本文將穩(wěn)健抗差估計(jì)理論應(yīng)用到基于羅德里格矩陣的坐標(biāo)轉(zhuǎn)換模型中，以達(dá)到抵抗粗差對(duì)坐標(biāo)轉(zhuǎn)換精度影響的目的。

將參數(shù)初值加上求解的七個(gè)轉(zhuǎn)換參數(shù)改正數(shù)后，坐標(biāo)轉(zhuǎn)換即可進(jìn)行。坐標(biāo)轉(zhuǎn)換后可以獲得公共點(diǎn)的坐標(biāo)差值，坐標(biāo)差值的大小可以反映轉(zhuǎn)換參數(shù)精度的高低。如果在解算坐標(biāo)轉(zhuǎn)換參數(shù)時(shí)公共點(diǎn)中的某個(gè)點(diǎn)精度較低，那么坐標(biāo)轉(zhuǎn)換后該公共點(diǎn)的坐標(biāo)差值會(huì)比其它公共點(diǎn)的坐標(biāo)差值更大，根據(jù)穩(wěn)健抗差估計(jì)理論，可通過(guò)式(12)對(duì)公共點(diǎn)重新定權(quán)，從而降低精度較低的公共點(diǎn)對(duì)降低坐標(biāo)轉(zhuǎn)換精度的影響。

(12)

2 模擬計(jì)算與分析

現(xiàn)模擬一套數(shù)據(jù)作為算例進(jìn)行計(jì)算與分析。設(shè)一個(gè)空間直角坐標(biāo)系下有4個(gè)點(diǎn)p1～p4，現(xiàn)將這些點(diǎn)繞X軸，Y軸，Z軸分別旋轉(zhuǎn)20°，30°，35°，再沿X軸，Y軸，Z軸分別平移230 m，170 m，75 m，設(shè)尺度比λ=1，最終得到在新坐標(biāo)系下4個(gè)點(diǎn)q1～q4的坐標(biāo)，各點(diǎn)的坐標(biāo)數(shù)據(jù)如表1所示。

方案1：選取點(diǎn)1，2，3作為公共點(diǎn)，計(jì)算基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換的初值，然后將這些初值代入坐標(biāo)轉(zhuǎn)換模型(1)中，對(duì)原坐標(biāo)系中的4個(gè)點(diǎn)p1，p2，p3，p4進(jìn)行坐標(biāo)轉(zhuǎn)換，計(jì)算轉(zhuǎn)換后的坐標(biāo)與新坐標(biāo)系下點(diǎn)q1，q2，q3，q4的坐標(biāo)差值，記為|Δ1|，再計(jì)算其點(diǎn)位偏差，記為m1，見(jiàn)表2。完成基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換參數(shù)初值的計(jì)算后，將7個(gè)初值代入到基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換線性化模型(3)，可得到誤差方程式(4)，利用最小二乘法解算7個(gè)轉(zhuǎn)換參數(shù)的改正數(shù)，然后將7個(gè)轉(zhuǎn)換參數(shù)初值加上各自的改正數(shù)后再代入式(1)，對(duì)原坐標(biāo)系中的4個(gè)點(diǎn)p1，p2，p3，p4進(jìn)行坐標(biāo)轉(zhuǎn)換，計(jì)算轉(zhuǎn)換后的坐標(biāo)與新坐標(biāo)系下點(diǎn)q1，q2，q3，q4的坐標(biāo)差值，記為|Δ2|，再計(jì)算其點(diǎn)位偏差m2，見(jiàn)表2。

表1 已知點(diǎn)坐標(biāo)數(shù)據(jù)

方案2：對(duì)原坐標(biāo)系下點(diǎn)p1，分別在X，Y，Z方向加入1 cm，2 cm，2 cm的粗差，然后選取點(diǎn)1，2，3作為公共點(diǎn)，首先通過(guò)傳統(tǒng)的基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換模型對(duì)原坐標(biāo)系中的4個(gè)點(diǎn)p1，p2，p3，p4進(jìn)行坐標(biāo)轉(zhuǎn)換，計(jì)算轉(zhuǎn)換后的坐標(biāo)與新坐標(biāo)系下點(diǎn)q1，q2，q3，q4的坐標(biāo)差值，記為|Δ1|，再計(jì)算其點(diǎn)位偏差，記為m1，見(jiàn)表2。采用本文所述基于羅德里格矩陣的抗差迭代坐標(biāo)轉(zhuǎn)換模型，對(duì)原坐標(biāo)系中4個(gè)點(diǎn)p1，p2，p3，p4進(jìn)行坐標(biāo)轉(zhuǎn)換，計(jì)算轉(zhuǎn)換后的坐標(biāo)與新坐標(biāo)系下點(diǎn)q1，q2，q3，q4的坐標(biāo)差值，記為|Δ2|，再計(jì)算其點(diǎn)位偏差m2，見(jiàn)表2。

1)方案1分析。對(duì)比兩次計(jì)算所得的各點(diǎn)的點(diǎn)位偏差可知，僅利用計(jì)算所得的羅德里格坐標(biāo)轉(zhuǎn)換模型七參數(shù)初值進(jìn)行坐標(biāo)轉(zhuǎn)換，所得結(jié)果精度很低，而通過(guò)七參數(shù)初值和采用線性模型計(jì)算其參數(shù)的改正數(shù)，并將改正數(shù)與初值融合后再進(jìn)行坐標(biāo)轉(zhuǎn)換，精度顯著提高。事實(shí)上，利用最小二乘原理進(jìn)行迭代計(jì)算所得到的轉(zhuǎn)換參數(shù)改正數(shù)的效果會(huì)更好，由于篇幅有限，本文不再陳述。

2)方案2分析。由于p1點(diǎn)精度低、誤差過(guò)大，相較于方案1中的第二種算法，各點(diǎn)坐標(biāo)轉(zhuǎn)換后的點(diǎn)位偏差明顯變大，尤其是p1點(diǎn)轉(zhuǎn)換后的點(diǎn)位偏差最大。綜合分析表2中方案2的m1和m2，可見(jiàn)采用穩(wěn)健估計(jì)原理進(jìn)行坐標(biāo)轉(zhuǎn)換，除了p1點(diǎn)外，其它各點(diǎn)坐標(biāo)轉(zhuǎn)換后的點(diǎn)位偏差都減小了。對(duì)于點(diǎn)p1，雖然其點(diǎn)位偏差變大了，但其實(shí)它是更接近于該點(diǎn)真實(shí)值的。

3 結(jié)束語(yǔ)

本文將羅德里格矩陣和穩(wěn)健抗差估計(jì)理論綜合應(yīng)用于空間直角坐標(biāo)轉(zhuǎn)換，研究適用于任意旋轉(zhuǎn)角的空間直角坐標(biāo)轉(zhuǎn)換的抗差迭代算法，采用C#語(yǔ)言進(jìn)行了編程實(shí)現(xiàn)，利用仿真數(shù)據(jù)和算法程序進(jìn)行實(shí)驗(yàn)計(jì)算與分析。方案1的兩次計(jì)算結(jié)果表明在無(wú)粗差情況下，利用羅德里格矩陣坐標(biāo)轉(zhuǎn)換線性化模型計(jì)算轉(zhuǎn)換參數(shù)改正數(shù)后再進(jìn)行坐標(biāo)轉(zhuǎn)換，相較于直接利用轉(zhuǎn)換參數(shù)初值進(jìn)行坐標(biāo)轉(zhuǎn)換的精度會(huì)明顯提高。方案2的兩次計(jì)算結(jié)果表明在公共點(diǎn)含粗差情況下，將穩(wěn)健抗差估計(jì)方法應(yīng)用于基于羅德里格矩陣的坐標(biāo)轉(zhuǎn)換是可行的，能夠有效地抵抗粗差對(duì)轉(zhuǎn)換結(jié)果的影響。

表2 各點(diǎn)轉(zhuǎn)換后坐標(biāo)差值及其點(diǎn)位偏差 cm

[1] 劉大杰，施一民，過(guò)靜珺. 全球定位系統(tǒng)(GPS)的原理與數(shù)據(jù)處理[M]. 上海：同濟(jì)大學(xué)出版社，1999.

[2] 陳義，沈云中，劉大杰. 適用于大旋轉(zhuǎn)角的三維基準(zhǔn)轉(zhuǎn)換的一種簡(jiǎn)便模型[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2004，29(12)：1101-1104.

[3] 曾文憲，陶本藻. 三維坐標(biāo)轉(zhuǎn)換的非線性模型[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2003，28(5)：566-568.

[4] SCHENK T. From Point Based to Feature-Based Aerial Triangulation[J]. ISPRS Journal of Photogrammetry & Remote Sensing，2004(58)：315-329.

[5] 姜柱，劉慶元，周曼. 單位四元數(shù)在坐標(biāo)轉(zhuǎn)換中的應(yīng)用[J]. 礦山測(cè)量，2012(5)：28-30.

[6] 韓夢(mèng)澤，李克昭. 基于羅德里格矩陣的空間坐標(biāo)轉(zhuǎn)換[J]. 測(cè)繪工程，2016(4)：25-27.

[7] 楊凡，李廣云，王力，等. 一種基于羅德里格矩陣的最小二乘迭代坐標(biāo)轉(zhuǎn)換方法[J]. 工程勘察，2010(9)：80-84.

[8] 郭英起，唐彬，張秋江，等. 基于空間直角坐標(biāo)系的高精度坐標(biāo)轉(zhuǎn)換方法研究[J]. 大地測(cè)量與地球動(dòng)力學(xué)，2012(2)：125-128.

[9] 潘國(guó)榮，趙鵬飛. 基于空間向量的三維基準(zhǔn)轉(zhuǎn)換模型[J]. 大地測(cè)量與地球動(dòng)力學(xué)，2009(6)：79-82.

[10] 黃維彬. 近代平差理論及其應(yīng)用[M]. 北京:八一出版社，1992.