李國琴，田林亞，郭英起，張 洋，畢繼鑫

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

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

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

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

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

(1)

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

(2)

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

(3)

其中，λ0為尺度比參數(shù)初值，R0為羅德里格矩陣初值，ΔX0，ΔY0，ΔZ0為三個平移參數(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].

為了實現(xiàn)基于羅德里格矩陣坐標(biāo)轉(zhuǎn)換，需要模型參數(shù)初值計算。首先，計算尺度比參數(shù)λ的初值。可以通過位于兩套坐標(biāo)系的兩個公共點的坐標(biāo)進(jìn)行計算，即

(6)

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

(7)

其中，

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

Zq21=Zq2-Zq1,

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

Zp21=Zp2-Zp1.

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

(8)

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

(9)

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

(10)

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

(11)

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

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

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

(12)

2 模擬計算與分析

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

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

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

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

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

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

3 結(jié)束語

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

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

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

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

[3] 曾文憲，陶本藻. 三維坐標(biāo)轉(zhuǎn)換的非線性模型[J]. 武漢大學(xué)學(xué)報(信息科學(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]. 礦山測量，2012(5)：28-30.

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

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

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

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

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