基于單位四元數矩陣實表示的坐標轉換算法研究

秦 鋒,鹿 松,張振虎

〔1.中國核工業集團三門核電有限公司,浙江 臺州 317112;2.上?？辈煸O計研究院(集團)有限公司,上海 200093〕

0 引 言

空間三維坐標轉換一直是測量領域的重要研究內容。目前很多坐標轉換都涉及大旋轉角的空間坐標轉換,傳統的布爾莎七參數模型只適用于小角度旋轉,難以滿足要求。當前關于大角度坐標轉換模型算法一般有誤差方程線性化算法[1-4]、羅德里格矩陣算法[5-9]、四元數法[10-14]、奇異值分解(SVD)算法等[15-16]。其中,誤差方程線性化算法利用一階微分將誤差方程線性化后迭代求解,該算法存在計算復雜、迭代不收斂或收斂到非極值,以及求解過程中可能會出現矩陣奇異或矩陣條件數過大等問題。羅德里格矩陣算法利用反對稱矩陣構建線性方程進行求解,根據文獻[17]論證其精度略低于七參數迭代模型。奇異值分解法其數學基礎為矩陣分析中子空間的旋轉問題[18],該算法適用性強,但需要進行奇異值分解,難以普遍應用。四元數法利用單位四元數表示旋轉,不用進行迭代,可提供平滑差值和避免萬向鎖問題,其用于坐標轉換參數求解時,適用性強,可以用于任意角度的三維坐標轉換問題求解,但四元數法理論較復雜,其乘積計算不便。本文基于單位四元數矩陣實表示的方式[19]將四元數體上的求解轉化成矩陣上的求解,并給出了利用單位四元數進行坐標轉化的兩種算法的證明過程,公式推導簡便,容易理解,經數據驗證,兩種算法解算結果準確一致。

1 四元數基礎知識

1.1 四元數定義及性質

四元數是簡單的超復數,由1個實部q0和3個虛部q1、q2、q3組成,可表示為:

(1)

式中,i、j、k滿足i2=j2=k2=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j。其單位四元數滿足:

(2)

其共軛四元數為:

(3)

四元數模定義為:

(4)

四元數模滿足以下公式:

(5)

(6)

四元數的逆:

(7)

1.2 四元數運算

四元數加法如下:

(8)

四元數乘法(格拉斯曼積)如下:

(9)

從式(8)和式(9)容易驗證,四元數加法滿足結合律和交換律,乘法滿足結合律和分配律,但不滿足交換律,這導致了四元數之間運算復雜。為了簡化四元數之間的運算,根據參考文獻[19],借助四元數的矩陣實表示,可以將四元數域上的運算轉化為實數域上向量和矩陣的運算。實數域上向量和矩陣的運算比四元數域的計算簡單,且計算過程明了,便于公式推導和實際應用。

1.3 四元數的矩陣實表示

(10)

(11)

(12)

(13)

(14)

(15)

(16)

2 計算原理及公式推導

2.1 空間坐標轉換模型

空間坐標轉換模型可表示為:

(17)

將坐標轉換模型以向量形式表示:

(18)

為便于計算,將坐標進行中心化處理:

(19)

由式(18)和式(19)可以推導出平移參數計算公式:

(20)

而式(18)變為:

(21)

2.2 利用單位四元數進行坐標轉換的第一種方法證明

構造如下優化函數f1:

(22)

因R為正交矩陣,故有‖Rri‖2=‖ri‖2,上式化簡為:

(23)

(24)

注意:上式的“·”表示四元數的點積,其結果為一個標量。

根據引理1和引理2,將式(24)采用矩陣實表示進行表示:

(25)

(26)

(27)

式(27)為關于λ的一元二次方程,尺度系數計算公式如式(28)所示。當λ取式(28)時,f1取最小值。

(28)

求得單位四元數后,可根據下式計算旋轉矩陣R,再結合式(20)及式(28)可解出平移參數。

(29)

2.3 利用單位四元數進行坐標轉換的第二種方法證明

將式(22)用單位四元數表示如下:

(30)

采用矩陣實表示表達四元數:

(31)

(32)

則函數f1可表示為:

(33)

(34)

利用特征分解計算出單位四元數后,將式(30)展開,得:

(35)

式(32)為關于λ的一元二次方程,當λ取式(33)值時,f1取最小值。

(36)

求得單位四元數后,可根據式(29)計算旋轉矩陣R,再結合式(20)及式(36)可解出平移參數。

3 算例分析

算例一:

以激光跟蹤儀測量的某6個公共點轉站前和轉站后的空間三維坐標為例,公共點坐標如表1所示,采用本文兩種算法計算坐標轉換單位四元數、尺度系數和平移參數,并將四元數轉化為旋轉矩陣。另采用奇異值分解法(SVD)計算結果作為對比。

表1 激光跟蹤儀轉站前后公共點測量數據(單位/mm)

在對數據進行中心化處理后,采用本文第一種算法得到的特征分解矩陣M、特征向量矩陣U1、特征值矩陣V1如式(37)-(39)所示。

(37)

(38)

(39)

從上式可以看出,V1中最大特征值對應U1中的第四列特征向量,即為所求單位四元數。

本文第二種計算方法先假定尺度系數λ為1,在對數據進行中心化處理后,可計算得到特征分解矩陣N、特征向量矩陣U2及特征值矩陣V2。V2中最小特征值對應U2中的特征向量即為所求單位四元數。因篇幅所限,不再一一展示矩陣N、U2及V2。

本文兩種方法所求的單位四元數、尺度系數、平移參數如表2所示。兩種方法計算的單位四元數符號相反,但其絕對值一致,不影響旋轉矩陣的計算;尺度系數和平移參數完全一致,說明本文兩種方法是一致的。

表2 本文兩種算法計算的坐標轉換參數

為便于與奇異值分解算法(SVD)計算結果進行對比,采用表2中計算得到的四元數,根據公式(29)計算本文兩種方法得到的旋轉矩陣。因本文兩種方法得到的旋轉矩陣一致,故只展示第一種方法計算得到的旋轉矩陣R。

(40)

奇異值分解算法(SVD)計算出的旋轉矩陣、尺度系數及平移參數如表3所示。

表3 奇異值分解算法(SVD)計算的坐標轉換參數

對比表2和表3可以看出,本文算法和SVD法計算的結果高度一致,說明本文算法是正確的。

算例二:

以文獻[3]中算例為例,計算坐標轉換參數。首先將文獻[3]中設計坐標系轉換成左手系,與測量坐標系一致;根據本文兩種算法,計算單位四元數、尺度系數、平移參數及旋轉矩陣。因本文兩種算法計算結果完全一致,故只展示第一種算法計算結果。計算結果如表4所示。

表4 文獻[3]算例本文算法計算結果

本文算法與文獻[3]和文獻[4]計算的坐標較差數據如表5所示,可發現本文算法與文獻[4]的坐標較差數據完全一致。

表5 坐標較差數據(單位/mm)

4 結 語

本文基于單位四元數矩陣實表示的方式推導并證明了利用單位四元數進行空間三維坐標轉換的兩種算法,兩種算法論證過程清晰,適用性強,不用進行迭代,方便編程,適合大角度的空間坐標轉換。