鹿珂珂，劉陵順，唐大全

（海軍航空大學(xué) 航空作戰(zhàn)勤務(wù)學(xué)院，山東，煙臺 264001）

1843 年，Hamilton 首先提出了四元數(shù)的數(shù)學(xué)概念[1]，用于解決空間矢量類似平面問題的復(fù)數(shù)方法.隨著科學(xué)技術(shù)的發(fā)展，四元數(shù)在機器人和航空航天領(lǐng)域得到廣泛應(yīng)用[2]，成為分析研究剛體運動的有效工具.此外，四元數(shù)也被大量用于電腦繪圖及相關(guān)的圖像分析上，用于表示三維物體的旋轉(zhuǎn)及姿態(tài)，還見于控制論、信號處理、姿態(tài)控制、物理、軌道力學(xué)和生物信息學(xué)等學(xué)科領(lǐng)域[3].

1990 年，SHUSTER 等[4]首次提出四元數(shù)乘法與旋轉(zhuǎn)矩陣乘法的順序相反，并提出了與Hamilton 四元數(shù)乘法定義不同的乘法表達(dá)約定.基于該文獻，美國國家航空航天局（national aeronautics and space administration，NASA）的噴氣推進實驗室(jet propulsion laboratory，JPL)提出了另外一種四元數(shù)的表達(dá)約定[5]，與Hamilton 約定對應(yīng)，本文將這種表達(dá)稱為JPL 約定.

此后，在學(xué)術(shù)界使用四元數(shù)的過程中，開始出現(xiàn)混亂和不嚴(yán)謹(jǐn)，也帶來了困惑和爭論[3,6?8].文獻[3]整理了不同領(lǐng)域的部分論文，目前大多數(shù)學(xué)科領(lǐng)域和數(shù)學(xué)分析使用Hamilton 約定，包括大多數(shù)開源軟件包，比如，Eigen[9]、Google Ceres[10]等.目前存在混淆的主要是航空航天和機器人領(lǐng)域.在航空航天領(lǐng)域，國外JPL 約定主要出現(xiàn)在與NASA 相關(guān)的文獻中[4?5,11?12].在機器人領(lǐng)域，則出現(xiàn)了很多使用JPL 約定的經(jīng)典論文[13?14]，受其影響，該領(lǐng)域四元數(shù)混淆情況尤其嚴(yán)重.國內(nèi)由于受經(jīng)典著作[15]的影響，四元數(shù)通常默認(rèn)采用Hamilton 約定，然而，也有一些文獻[16?18]在使用四元數(shù)時，沒有顯式說明采用哪種表達(dá)約定情況下，使用了JPL 約定.四元數(shù)的兩種不同表達(dá)約定造成了大量的混淆，尤其是將四元數(shù)作為工具的文獻[19?21]，不加說明地進行了使用，增加了公式推導(dǎo)和實踐的難度，并成為實踐中的錯誤來源.

文獻[8]首先提出了存在兩種四元數(shù)表達(dá)帶來的問題，重點對數(shù)學(xué)表達(dá)形式上的不同進行了分析，文獻[3]和[7]則從主動旋轉(zhuǎn)和被動旋轉(zhuǎn)的角度，理論上分析論證了兩種四元數(shù)的起源和不同，文獻[6]基于Hamilton 四元數(shù)進行了詳細(xì)的數(shù)學(xué)推導(dǎo)，并提到了未選擇JPL 的原因是由于Hamilton 四元數(shù)應(yīng)用的廣泛性，然而上述研究均缺少對兩種表達(dá)方式的實驗驗證分析，并且沒有清晰地給出文獻閱讀和軟件工具中判別和轉(zhuǎn)換四元數(shù)的方法.因此，本文在相關(guān)研究基礎(chǔ)上，詳細(xì)推導(dǎo)分析了Hamilton 和JPL 約定下旋轉(zhuǎn)四元數(shù)的異同，論證了Hamilton 四元數(shù)在處理連續(xù)旋轉(zhuǎn)時的同形性解釋，推導(dǎo)了兩種四元數(shù)表達(dá)方式的判別、測試和轉(zhuǎn)換方法，并使用實驗數(shù)據(jù)進行了驗證.提議在使用四元數(shù)時，明確指出所采用的表達(dá)約定，普遍采用Hamilton 約定，逐漸舍棄JPL約定.

1 三維剛體的旋轉(zhuǎn)表達(dá)

式中Rwb為 旋轉(zhuǎn)矩陣，表示坐標(biāo)系 Fw到 Fb的旋轉(zhuǎn)，滿足如下約束

圖1 三維剛體固聯(lián)坐標(biāo)系旋轉(zhuǎn)示意圖Fig.1 Schematic diagram of 3D rigid body fixed coordinate system rotation

三維剛體旋轉(zhuǎn)除了旋轉(zhuǎn)矩陣的表達(dá)，還可以使用旋轉(zhuǎn)向量、歐拉角和四元數(shù)[22?23].在個別文獻[3]中，為了區(qū)分Hamilton 和JPL 約定下的乘法，采用了兩種四元數(shù)乘法符號，符號 ?用于JPL 約定，Hamilton約定則采用 ⊙，而文獻中常見的Hamilton 表達(dá)一般使用類似實數(shù)變量的乘法[1,4,6,13,21].為了表述方便，本文選擇在Hamilton 和JPL 四元數(shù)約定中采用相同的變量和運算符號，四元數(shù)乘法符號均采用 ?.

式中CH(·)為旋轉(zhuǎn)四元數(shù)到旋轉(zhuǎn)矩陣的映射.對應(yīng)式（3）的關(guān)系式為

2 Hamilton 與JPL 約定的比較

2.1 與旋轉(zhuǎn)矩陣乘法的同形性

選擇在Hamilton 約定基礎(chǔ)上，提出另外的四元數(shù)表達(dá)約定，原因在于部分學(xué)者認(rèn)為旋轉(zhuǎn)四元數(shù)與旋轉(zhuǎn)矩陣的映射關(guān)系滿足[13]

據(jù)此，根據(jù)式（6）可以進一步推導(dǎo)得到

也就是

由式（11）可知，在式（9）的假設(shè)下，處理連續(xù)旋轉(zhuǎn)時，旋轉(zhuǎn)矩陣與Hamilton 四元數(shù)約定不滿足同形性.但是，事實上Hamilton 約定下的四元數(shù)與旋轉(zhuǎn)矩陣的映射關(guān)系由式（5）給出，在此條件下，由式（6）得到

因而

由式（13）可見，Hamilton 四元數(shù)在表示連續(xù)旋轉(zhuǎn)時，按照與旋轉(zhuǎn)矩陣的映射關(guān)系CH(·)，維持了與旋轉(zhuǎn)矩陣相同的鏈?zhǔn)椒▌t，滿足與旋轉(zhuǎn)矩陣乘法的同形性.

在JPL 約定下，滿足

JPL 約定在處理連續(xù)旋轉(zhuǎn)時，形式上較為整潔，但是結(jié)論與Hamilton 約定一致.

2.2 Hamilton 和JPL 約定區(qū)別

Hamilton 和JPL 四元數(shù)表達(dá)約定的區(qū)別體現(xiàn)主要在如下幾個方面：

① 四元數(shù)定義中相互正交單位向量的乘法定義不同.Hamilton 約定中，i?j=?j?i=k，JPL 約定中，j?i=?i?j=k.

②元素排列順序通常不同.設(shè)四元素實部為qw，

③相同四元數(shù)對應(yīng)旋轉(zhuǎn)向量相同，但是表示的旋轉(zhuǎn)不同.三維剛體的旋轉(zhuǎn)可以使用旋轉(zhuǎn)向量來表示，旋轉(zhuǎn)向量的方向為轉(zhuǎn)動軸，轉(zhuǎn)過的角度為向量的長度.設(shè)旋轉(zhuǎn)向量的方向為單位向量uwR（上標(biāo)表示向量表示在 Fw中），旋轉(zhuǎn)角度為θ，Hamilton 約定和JPL約定下，旋轉(zhuǎn)四元數(shù)與旋轉(zhuǎn)向量的對應(yīng)關(guān)系均為

但是對于相同的旋轉(zhuǎn)向量，Hamilton 約定采用右手定則，JPL 采用左手定則判斷四元數(shù)表示的旋轉(zhuǎn).相同的四元數(shù)表達(dá)=[qw qv]T，Hamilton 和JPL 約定下旋轉(zhuǎn)矩陣分別為

式中：I3×3為 單位矩陣，下標(biāo)表示矩陣維數(shù)是 3×3；[·]×表示將向量映射成反對稱矩陣.根據(jù)式（15），使用李代數(shù)[25]通過旋轉(zhuǎn)向量來得到旋轉(zhuǎn)矩陣則分別為

④四元數(shù)相乘結(jié)果形式不同.根據(jù)四元數(shù)實部和虛部的排列順序不同，在進行四元數(shù)乘法運算時，左乘四元數(shù)和 右乘四元數(shù)的形式存在不同，詳見表1.

表1 四元數(shù) 對應(yīng)左乘和右乘矩陣形式Tab.1 Left-multiplied and right-multiplied matrix forms corresponding to quaternion

表1 四元數(shù) 對應(yīng)左乘和右乘矩陣形式Tab.1 Left-multiplied and right-multiplied matrix forms corresponding to quaternion

序號 四元數(shù) Hamilton JPL]■■■■■■[pv■■■■■■]× pv?pTv 0■■■■■■ ˉq=■■■■■■?[pv■■■■■■qv qw L(ˉp)=pwI4×4+× pv?pTv 0■■■■■■■■■■■■L(ˉp)=pwI4×4+1ˉp=■■■■■■pv pw]]R(ˉq)=qwI4×4+■■■■■■?[qv × qv?qTv 0■■■■■■R(ˉq)=qwI4×4+■■■■■■[qv × qv?qTv 0■■■■■■■■■■■■ 0 ?pTv pv■■■■■■ ˉq=[pv■■■■■■qw qv■■■■■■■■■■■■L(ˉp)=pwI4×4+]L(ˉp)=pwI4×4+]■■■■■■2■■■■■■ 0 ?pTv pv ?[pv××ˉp=■■■■■■pw pv R(ˉq)=qwI4×4+■■■■■■ 0 ?qTv qv ?[qv]■■■■■■R(ˉq)=qwI4×4+]■■■■■■×■■■■■■ 0 ?qTv qv[qv ×

⑤四元數(shù)的微分形式不同.對于Hamilton 約定，假設(shè)旋轉(zhuǎn)運動過程中t時刻的旋轉(zhuǎn)四元數(shù)為qˉbw(t)，可以得到

根據(jù)式（6），設(shè)角速度在坐標(biāo)系 Fw中的表示為那么

對于JPL 約定，旋轉(zhuǎn)四元數(shù)虛部在前，實部在后，在相同的假設(shè)下，計算坐標(biāo)系 Fb相對于 Fw旋轉(zhuǎn)四元

數(shù)的微分時，有

式中：

以上是Hamilton 和JPL 四元數(shù)約定中存在的一些主要區(qū)別.基于上述結(jié)論，可以推導(dǎo)得到不同約定下不同形式的旋轉(zhuǎn)四元數(shù)的積分、擾動方程以及對擾動的雅克比矩陣求解等，具體可以參考文獻[6].

3 Hamilton 和JPL 約定的區(qū)分與轉(zhuǎn)換

3.1 文獻中四元數(shù)約定區(qū)分

在閱讀和引用機器人和航空航天領(lǐng)域論文的時候，尤其需要注意四元數(shù)采用的約定，在這兩個學(xué)科領(lǐng)域存在兩種約定的廣泛使用，可以按照如下步驟進行四元數(shù)約定的區(qū)分：

①查看文獻中所采用四元數(shù)約定的說明.大部分文獻在使用四元數(shù)過程中，會明確說明采用了Hamilton 四元數(shù)，JPL 約定則很少會顯式說明.

②查看四元數(shù)定義和運算.對于兩種四元數(shù)約定區(qū)分最直觀的方法是查看i?j乘 積結(jié)果.若i?j=k，則采用Hamilton 約定；若i?j=?k，則采用JPL 約定.

③明確四元數(shù)表示旋轉(zhuǎn)的含義，根據(jù)規(guī)范四元數(shù)與旋轉(zhuǎn)矩陣的映射關(guān)系進行確認(rèn).

④根據(jù)四元數(shù)乘法形式確定.根據(jù)四元數(shù)相乘的先后順序，對比表1 中的乘法矩陣形式，確定所采用的四元數(shù)約定.

⑤根據(jù)旋轉(zhuǎn)四元數(shù)的微分方程判斷.通常角速度由慣性測量單元直接得到在當(dāng)前坐標(biāo)系的表示，此時可參照式（19）和（22），根據(jù)具體微分形式確定四元數(shù)約定.

可參照圖2 的具體流程確認(rèn)文獻中所采用的四元數(shù)約定.需要注意的是，要首先明確四元數(shù)的元素排列順序，以及四元數(shù)所表示旋轉(zhuǎn)的含義.如本文第1 節(jié)所述，同一旋轉(zhuǎn)，有兩種表述：一種是從參考坐標(biāo)系 Fw到當(dāng)前坐標(biāo)系 Fb的旋轉(zhuǎn)為與Rwb，另一種是將點坐標(biāo)（向量）從在參考坐標(biāo)系 Fw中的Pw旋轉(zhuǎn)到當(dāng)前坐標(biāo)系 Fb的表示Pb的 旋轉(zhuǎn)為與Rbw.兩種表示的旋轉(zhuǎn)矩陣互為逆矩陣，旋轉(zhuǎn)四元數(shù)互為共軛.不同文獻中四元數(shù)符號形式可能不同，比如，參考文獻[24]中的四元數(shù)與本文Hamilton 約定下的含義相同，對應(yīng)旋轉(zhuǎn)的第一種表述，即坐標(biāo)系bk到坐標(biāo)系bk+1的 旋轉(zhuǎn)，參考文獻[19]中的，采用的是第二種旋轉(zhuǎn)表述，表示將坐標(biāo)系G中的向量轉(zhuǎn)換到坐標(biāo)系L，對應(yīng)本文JPL 約定下的四元數(shù)由于兩種四元數(shù)約定的廣泛應(yīng)用，不排除即使文獻中明確說明所采用的四元數(shù)約定，但由于作者沒有意識到存在兩種四元數(shù)約定，受參考文獻影響，在使用四元數(shù)進行公式推導(dǎo)過程中，仍然出現(xiàn)混淆和錯誤的情況.

圖2 四元數(shù)約定區(qū)分流程圖Fig.2 Flow chart of quaternion convention identification

3.2 軟件工具中四元數(shù)約定測試和轉(zhuǎn)換

對于軟件工具，可以測試i?j乘積結(jié)果，若為i?j=k，則采用Hamilton 約定；若為i?j=?k，則采用JPL 約定.也可以通過測試四元數(shù)到旋轉(zhuǎn)矩陣的轉(zhuǎn)換，得到軟件所采用的四元數(shù)約定.不同約定下四元數(shù)與旋轉(zhuǎn)矩陣的數(shù)量關(guān)系可以通過軟件輸出和手工計算結(jié)果進行對比.

3.3 Hamilton 和JPL 約定的應(yīng)用探討

兩種四元數(shù)的理論都是完善和成熟的，這也意味著在表達(dá)旋轉(zhuǎn)的時候，可以使用其中任意一種表達(dá)約定，并且從上面的推導(dǎo)可以看到，兩者在數(shù)學(xué)運算上存在相似性，可以通過簡單的數(shù)學(xué)運算實現(xiàn)相互轉(zhuǎn)換，并且在計算資源占用上不存在顯著差異.但是，同時存在兩種表達(dá)約定，則導(dǎo)致了研究人員需要額外的成本來區(qū)分和轉(zhuǎn)換使用，而這項工作并不會帶來額外的收益，并且引入潛在的隱患.如果要選擇保留其中一種表達(dá)約定的話，相比之下，Hamilton 四元數(shù)提出距今已有近180 年，歷史時間悠久，受眾和應(yīng)用領(lǐng)域廣泛，根據(jù)參考文獻[3]，兩種四元數(shù)應(yīng)用領(lǐng)域部分文獻初步調(diào)研統(tǒng)計結(jié)果如表2 所示.可見，除了航空和機器人領(lǐng)域，其他諸如數(shù)學(xué)、機械等領(lǐng)域，均采用了Hamilton 約定，大多數(shù)常用軟件和科普網(wǎng)站也均采用了Hamilton 約定.此外，從前面的數(shù)學(xué)推導(dǎo)也可以看到，兩種四元數(shù)在表達(dá)旋轉(zhuǎn)時，均具備同形性，雖然Hamilton 四元數(shù)的連續(xù)旋轉(zhuǎn)與直覺不太一致，但是Hamilton 四元數(shù)的旋轉(zhuǎn)采用右手直角坐標(biāo)系，表達(dá)更符合人們的使用習(xí)慣.因此，建議在學(xué)術(shù)研究領(lǐng)域更為廣泛地應(yīng)用Hamilton 表達(dá)約定，逐漸舍棄JPL 的表達(dá)約定，消除兩種約定同時存在和使用帶來的困惑和使用成本.

表2 兩種四元數(shù)約定的應(yīng)用領(lǐng)域統(tǒng)計[3]Tab.2 Statistics in the application field of the two quaternion conventions[3]

4 仿真驗證

Eigen 庫是用來進行線性代數(shù)、矩陣、向量操作等運算的C++庫，作為跨平臺使用的軟件庫，在機器人的嵌入式應(yīng)用中大量使用，而Matlab 作為商業(yè)數(shù)學(xué)軟件，在科學(xué)研究和工程設(shè)計領(lǐng)域應(yīng)用廣泛，選擇這兩個有代表性的軟件，分別對其中的四元數(shù)約定進行測試.

選取文獻[26]的航姿參考系統(tǒng)（attitude and heading reference system，AHRS）傳感器數(shù)據(jù)，分別利用Hamilton 和JPL 四元數(shù)對Mahony 算法效率和運算結(jié)果進行驗證.其中，傳感器數(shù)據(jù)是從AHRS 設(shè)備中記錄得到，經(jīng)過標(biāo)定校準(zhǔn)的陀螺儀、加速度計和磁力計數(shù)據(jù)，運動過程是順序的繞X,Y,Z軸從0°轉(zhuǎn)動到90°，然后到?90°.傳感器數(shù)據(jù)如圖3 所示.

圖3 傳感器數(shù)據(jù)Fig.3 Sensor data

所采用的實驗環(huán)境為ThinkPad 翼14 筆記本電腦，處理器為i7-10510U，操作系統(tǒng)采用Ubuntu 18.04，使用Ubuntu 版本的Matlab R2018b 64bit.首先分別對兩種四元數(shù)與旋轉(zhuǎn)矩陣、姿態(tài)角和軸角之間的相互轉(zhuǎn)換關(guān)系進行測試，相同轉(zhuǎn)換運算下，兩者的用時分別是JPL 四元數(shù)0.011 038 s，Hamilton 四元數(shù)0.011 022 s，可見在忽略系統(tǒng)資源調(diào)度的情況下，兩者用時幾乎一致，主要源于在不同約定下，轉(zhuǎn)換約定的數(shù)學(xué)運算除去符號和順序等，均保持一致.

Mahony 算法采用Hamilton 四元數(shù)進行推導(dǎo)，使用JPL 四元數(shù)進行改寫時，需要依據(jù)表1 修改四元數(shù)乘法，采用式（17）的坐標(biāo)系轉(zhuǎn)換關(guān)系，以及式（22）的姿態(tài)更新方程.通過多次運行程序，對運算所需時間取平均，得到利用Hamilton 四元數(shù)進行姿態(tài)解算用時0.132 963 s，歐拉姿態(tài)角見圖4（a），四元數(shù)變化見圖4（b）.利用JPL 四元數(shù)進行姿態(tài)解算用時0.133 241 s，歐拉姿態(tài)角見圖5（a），四元數(shù)變化見圖5（b）.從運行時間上來看，兩者由于較大數(shù)據(jù)量的處理，存在細(xì)微的不同，但該時差在系統(tǒng)運行的誤差范圍內(nèi).這符合理論上的推導(dǎo)結(jié)果，由Hamilton 約定改為JPL 約定，只需要同時變更四元數(shù)虛部的符號，并且更改相應(yīng)的運算公式，而運算轉(zhuǎn)換之間本身耗時基本一致，因此，兩種約定在所需計算資源上一致.

圖4 采用Hamilton 四元數(shù)的解算結(jié)果Fig.4 Solving results with Hamilton quaternion

圖5 采用JPL 四元數(shù)的解算結(jié)果Fig.5 Solving results with JPL quaternion

對比圖4 和圖5，可以發(fā)現(xiàn)，兩者姿態(tài)和四元數(shù)輸出結(jié)果也幾乎完全一致，細(xì)微差別源于算法運行過程中數(shù)據(jù)維護的細(xì)微差異，如計算精度和四元數(shù)歸一化等，以及由于萬向鎖原因，在俯仰角 θ接近90?時，歐拉角的奇異性使得滾轉(zhuǎn)角 ?和 偏航角 ψ時出現(xiàn)的區(qū)別，但是這與兩種四元數(shù)表達(dá)方式關(guān)系不大.需要注意的是，圖4（b）和圖5（b）中四元數(shù)數(shù)值變化一致，但是兩者表達(dá)的含義不同，參見式（16）和式（17），而對于四元數(shù)轉(zhuǎn)換得到的歐拉角均為按照常規(guī)定義給出的結(jié)果，因此也得到了相同結(jié)果.

5 結(jié) 論

處理四元數(shù)連乘時，為了滿足與旋轉(zhuǎn)矩陣乘法的同形性，學(xué)者和機構(gòu)提出了四元數(shù)的JPL 表達(dá)約定，本文通過系統(tǒng)分析對比Hamilton 和JPL 四元數(shù)約定下的四元數(shù)定義和推導(dǎo)，指出在四元數(shù)的Hamilton約定下，也同樣具有同形性，Hamilton 和JPL 四元數(shù)具有相同的運算結(jié)果和效率，鑒于Hamilton 四元數(shù)更為悠久的歷史，眾多學(xué)科領(lǐng)中更為廣泛的應(yīng)用，建議在航空航天和機器人領(lǐng)域也一致使用Hamilton 約定，舍棄JPL 四元數(shù)約定，并且在相關(guān)研究中，明確說明所采用的四元數(shù)表達(dá)約定.本文提出的區(qū)分和轉(zhuǎn)換四元數(shù)約定方法，有助于提高研究人員閱讀旋轉(zhuǎn)四元數(shù)相關(guān)文獻時的效率，并避免在理論推導(dǎo)和實踐過程中出現(xiàn)難以發(fā)現(xiàn)的錯誤.