基于羅德里格矩陣的抗差迭代坐標轉換算法研究

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

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

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

1 基于羅德里格矩陣的抗差迭代坐標轉換算法

1.1基于羅德里格矩陣坐標轉換模型

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

(1)

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

(2)

先對基于羅德里格矩陣的坐標轉換模型進行線性化，即

(3)

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

(4)

可解得轉換參數改正數

(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].

為了實現基于羅德里格矩陣坐標轉換，需要模型參數初值計算。首先，計算尺度比參數λ的初值。可以通過位于兩套坐標系的兩個公共點的坐標進行計算，即

(6)

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

(7)

其中，

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

Zq21=Zq2-Zq1,

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

Zp21=Zp2-Zp1.

根據反對稱矩陣S和羅德里格矩陣R的關系，可得

(8)

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

(9)

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

(10)

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

(11)

1.2 基于穩健抗差估計理論的選權迭代

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

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

(12)

2 模擬計算與分析

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

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

表1 已知點坐標數據

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

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

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

3 結束語

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

表2 各點轉換后坐標差值及其點位偏差 cm

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

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

[3] 曾文憲，陶本藻. 三維坐標轉換的非線性模型[J]. 武漢大學學報(信息科學版)，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] 姜柱，劉慶元，周曼. 單位四元數在坐標轉換中的應用[J]. 礦山測量，2012(5)：28-30.

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

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

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

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

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