任意旋轉(zhuǎn)角三維坐標(biāo)轉(zhuǎn)換的微分羅德里格矩陣算法

羅翔，劉志平

(1.中國礦業(yè)大學(xué) 江蘇省資源環(huán)境信息工程重點實驗室，江蘇 徐州 221116;2.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116)

0 引 言

三維坐標(biāo)轉(zhuǎn)換是空間大地測量應(yīng)用中坐標(biāo)成果換算的基礎(chǔ)工作，其常用模型包括Bursa-Wolf模型、Molodensky模型和武測模型。很多學(xué)者對上述七參數(shù)模型進(jìn)行，研究和改進(jìn)，如針對線性化系數(shù)矩陣病態(tài)問題的重心化方法[1]；解決參數(shù)相關(guān)性的兩步解法[2]；考慮參數(shù)先驗信息的序貫平差法[3]等。這些算法在一定程度上提高了模型適用性，但均將坐標(biāo)旋轉(zhuǎn)角視為微小值，僅適用于小旋轉(zhuǎn)角的坐標(biāo)轉(zhuǎn)換[4-7]。為解決任意旋轉(zhuǎn)角坐標(biāo)轉(zhuǎn)換問題，國內(nèi)外學(xué)者提出多種三維坐標(biāo)轉(zhuǎn)換方法：文獻(xiàn)[8-10]基于四元數(shù)理論[11-12]將旋轉(zhuǎn)矩陣用單位四元數(shù)表示，通過四元數(shù)的迭代更新獲得旋轉(zhuǎn)參數(shù)，適用于任意旋轉(zhuǎn)角，然而迭代存在繁瑣的矩陣求導(dǎo)過程；文獻(xiàn)[13-15]利用9個方向余弦參數(shù)表示旋轉(zhuǎn)矩陣，雖適用于任意旋轉(zhuǎn)角，但模型參數(shù)達(dá)到13個，迭代收斂慢；基于羅德里格矩陣?yán)碚?，文獻(xiàn)[16]采用非迭代解法，忽略了隨機(jī)模型對參數(shù)估計的影響；文獻(xiàn)[17-19]采用參數(shù)偏導(dǎo)數(shù)求解雅可比矩陣，形式復(fù)雜且收斂性差。本文針對上述空間坐標(biāo)轉(zhuǎn)換參數(shù)求解中矩陣求導(dǎo)復(fù)雜、迭代收斂慢、隨機(jī)模型不嚴(yán)密等問題，利用羅德里格矩陣、反對稱矩陣與向量積運算性質(zhì)，通過羅德里格矩陣微分與參數(shù)微分變換導(dǎo)出任意旋轉(zhuǎn)角三維坐標(biāo)轉(zhuǎn)換函數(shù)模型、隨機(jī)模型同步收斂的新方法，簡稱微分羅德里格矩陣算法。最后通過文獻(xiàn)算例和模擬算例將此方法與四元數(shù)、方向余弦、羅德里格矩陣算法進(jìn)行坐標(biāo)轉(zhuǎn)換對比研究，以評估該算法正確性及實用性。

1 三維坐標(biāo)轉(zhuǎn)換模型

1.1 基于羅德里格矩陣的轉(zhuǎn)換模型

設(shè)xi和Xi分別為空間一點在A,B空間直角坐標(biāo)系中的坐標(biāo)，則他們存在如式(1)的轉(zhuǎn)換關(guān)系：

xi=k·R·Xi+T-εi，Dε,

(1)

目前，任意角旋轉(zhuǎn)矩陣R可由旋轉(zhuǎn)角余弦、羅德里格矩陣、單位四元數(shù)等方式構(gòu)建。本文采用羅德里格矩陣表示的三維坐標(biāo)轉(zhuǎn)換模型，其任意角旋轉(zhuǎn)矩陣R表示形式為

(2)

1.2 微分羅德里格矩陣解法

將式(2)代入式(1)，即可得基于羅德里格矩陣的三維坐標(biāo)轉(zhuǎn)換模型，若按矩陣求導(dǎo)方式對該模型進(jìn)行參數(shù)線性化展開，其步驟繁瑣不易實現(xiàn)。下文利用羅德里格矩陣和反對陣矩陣的性質(zhì)構(gòu)造微分羅德里格矩陣算法，避免了復(fù)雜矩陣求導(dǎo)計算和化簡，推導(dǎo)過程如下。

由式(2)可以推導(dǎo)出羅德里格矩陣R和反對稱矩陣S滿足如下性質(zhì)，

(3)

將式(2)代入式(1)可得

式(4)兩端乘(I3-S)，可得

(I3-S)xi+(I3-S)εi=

k(I3+S)Xi+(I3-S)T，

(5)

將反對稱矩陣S，尺度參數(shù)k和平移向量T用初值項和微分項之和表示，即令S=S0+Sd，T=T0+dT，k=k0+dk，代入式(5)進(jìn)行整理化簡，并保留至微分項一階項，可得

(I3-S0)xi-k0(I3+S0)Xi-(I3-S0)T0+

(I3-S)εi=dk·(I3+S0)Xi-

Sxi+k0Xi-T0·d+(I3-S0)·dT,

(6)

式中，反對稱矩陣與向量積運算性質(zhì)滿足Sd·(*)=-S*·dα，*表示任意列向量。

進(jìn)一步，式(6)兩端乘(I3-Sα0)-1進(jìn)行還原，并結(jié)合式(3)化簡，可得

(7)

(8)

利用最小二乘法，可得模型參數(shù)改正數(shù)的迭代估計式，即

(9)

通過式(9)迭代更新參數(shù)估值，繼而更新系數(shù)矩陣A和觀測值L，迭代計算直至滿足迭代停止條件。分析上述推導(dǎo)過程可知，該方法充分利用了羅德里格矩陣、反對稱矩陣以及向量積運算的性質(zhì)，將羅德里格矩陣微分轉(zhuǎn)換為參數(shù)向量微分形式，避免了三維坐標(biāo)轉(zhuǎn)換模型中的矩陣求導(dǎo)，實現(xiàn)了三維坐標(biāo)轉(zhuǎn)換函數(shù)模型和隨機(jī)模型同步迭代逼近的微分羅德里格矩陣解法。

1.3 模型參數(shù)初值估計

(10)

根據(jù)式(1)將重心化坐標(biāo)模型表示為

基于式(11)，利用反對稱矩陣與向量積運算性質(zhì)S·(*)=-S*·α，可以導(dǎo)出反對稱矩陣參數(shù)初值的估值

(12)

進(jìn)一步，平移參數(shù)初值可按式(13)計算

(13)

2 計算結(jié)果與分析

2.1 現(xiàn)有文獻(xiàn)算例

為驗證本文所提的微分羅德里格矩陣算法的正確性和有效性，利用該方法對現(xiàn)有文獻(xiàn)算例進(jìn)行實驗分析。使用文獻(xiàn)[7]的數(shù)據(jù)，選取9對控制點進(jìn)行實驗?？刂泣c轉(zhuǎn)換前后坐標(biāo)系分別為WGS84坐標(biāo)系和西安80坐標(biāo)系，覆蓋范圍較大，但兩組坐標(biāo)旋轉(zhuǎn)角度相對較小。

表1 文獻(xiàn)算例轉(zhuǎn)換坐標(biāo)及其較差

表1中ΔX,ΔY,ΔZ為通過轉(zhuǎn)換參數(shù)求得的轉(zhuǎn)換坐標(biāo)與原始坐標(biāo)的較差。由表1可見，利用本文提出的微分羅德里格矩陣算法反算出來的坐標(biāo)，點位較差最大為36.201 mm，與文獻(xiàn)[7]中所計算36.2mm的最大點位較差一致。說明本文方法能夠準(zhǔn)確地求出小旋轉(zhuǎn)角的轉(zhuǎn)換參數(shù)，驗證了算法正確性和小旋轉(zhuǎn)角有效性。

2.2 模擬實驗算例

為更全面地驗證所提算法的正確性和大旋轉(zhuǎn)角有效性，選取經(jīng)度117.0°～117.8°，緯度33.0°～33.8°，以0.1°為步長形成9×9=81個格網(wǎng)點坐標(biāo)作為源坐標(biāo)，令大地高均為0，利用CGCS2000橢球參數(shù)計算各點空間直角坐標(biāo)，再將格網(wǎng)點坐標(biāo)基于模型式(1)按模擬轉(zhuǎn)換參數(shù)計算目標(biāo)坐標(biāo)。模擬參數(shù)取值方案為3個旋轉(zhuǎn)角取值范圍為[-90，90][-180，180][-180，180]，以15°為步長得到13×25×25=8 125組；平移參數(shù)固定為(x，y，z)=(1.0，1.5，2.0)×103；尺度參數(shù)固定為k=(1+10-5)。考慮到實際中存在的觀測誤差，在經(jīng)過上述參數(shù)理論值轉(zhuǎn)換后的坐標(biāo)三維分量中加入0.005 m的隨機(jī)誤差。

將本文微分羅氏矩陣與單位四元數(shù)法[8]、方向余弦法[14]、常規(guī)羅德里格矩陣[19]進(jìn)行參數(shù)求解試驗對比，每組試驗數(shù)據(jù)重復(fù)100次，迭代準(zhǔn)則不變，統(tǒng)計每組試驗數(shù)據(jù)迭代次數(shù)、單位權(quán)中誤差、參數(shù)真誤差等指標(biāo)評判方法的正確性和有效性。計算結(jié)果如表2所示。

表2 模擬算例計算結(jié)果統(tǒng)計

由表2知，4種算法均獲得了正確的解算參數(shù)，其中：尺度參數(shù)誤差均值量級在10-8量級；旋轉(zhuǎn)角參數(shù)誤差均值在10-3秒量級；平移參數(shù)誤差在10-1m量級，驗證了4種算法正確性和大旋轉(zhuǎn)角適用性。進(jìn)一步綜合比較迭代次數(shù)、單位權(quán)中誤差和參數(shù)精度指標(biāo)，可以得出本文算法有以下特點。

(1)更少的迭代次數(shù)。從迭代次數(shù)統(tǒng)計看，本文所提算法的迭代次數(shù)均值在5次以下，明顯優(yōu)于方向余弦算法(20次)、四元數(shù)算法(14次)和線性羅氏矩陣(19次)，且迭代次數(shù)中誤差小，說明該算法迭代穩(wěn)定、收斂速度快。

(2)更小的單位權(quán)中誤差。從單位權(quán)中誤差統(tǒng)計看，3種算法計算結(jié)果基本一致，均在給定的誤差水平(5 mm)，但從其統(tǒng)計量中誤差可看出，本文算法中誤差最小，表明本文算法比其他兩種方法更具有統(tǒng)計優(yōu)勢。

(3)更高的參數(shù)準(zhǔn)確性。從求解參數(shù)的真誤差統(tǒng)計可以看出：4種算法尺度參數(shù)真誤差幾乎完全一致；角度和平移參數(shù)真誤差中，本文算法誤差最大值、均值、中誤差均小于其他3種算法，這表明本文算法求解參數(shù)的準(zhǔn)確性優(yōu)于其他3種算法。

為了進(jìn)一步直觀表達(dá)新算法的優(yōu)點，統(tǒng)計迭代次數(shù)和參數(shù)真誤差的頻數(shù)直方圖，結(jié)果如圖 1所示。由圖1可看出，本文算法、方向余弦算法、四元數(shù)算法、常規(guī)羅氏矩陣迭代次數(shù)分別集中(累計60%以上)在10次以下、10～30次、10～20次、10～30次；尺度參數(shù)真誤差幾乎完全一致；角度參數(shù)真誤差分別集中在0～2×10-3″，0～4×10-3″，0～4×10-3″，0～4×10-3″間；平移參數(shù)真誤差分別集中在0～2，0～3，0～3，0～3 dm間。綜上，迭代次數(shù)和參數(shù)真誤差中，本文算法均優(yōu)于其他3種算法，表明其收斂速度快、參數(shù)準(zhǔn)確性高，與前文的分析結(jié)論一致。

圖1 迭代次數(shù)和參數(shù)真誤差統(tǒng)計

3 結(jié) 論

本文從常規(guī)任意旋轉(zhuǎn)角空間坐標(biāo)轉(zhuǎn)換模型線性化復(fù)雜、迭代收斂速度慢、精度低等問題出發(fā)，利用羅德里格矩陣和反對稱矩陣的微分形式推導(dǎo)出一種適用于任意旋轉(zhuǎn)角的空間坐標(biāo)轉(zhuǎn)換參數(shù)快速迭代求解算法——微分羅德里格矩陣算法。通過現(xiàn)有文獻(xiàn)算例和模擬算例分析，基于迭代次數(shù)、中誤差、參數(shù)真誤差等指標(biāo)，得出了以下結(jié)論：

(1)本文提出的微分羅德里格矩陣算法適用于任意旋轉(zhuǎn)角的三維坐標(biāo)轉(zhuǎn)換。

(2)通過大量模擬數(shù)據(jù)驗證，相較于經(jīng)典四元數(shù)算法、方向余弦算法、常規(guī)羅德里格矩陣算法，本文算法在收斂速度、參數(shù)精度上具有顯著優(yōu)勢。