鹿珂珂，劉陵順，唐大全

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

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

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

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

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

1 三維剛體的旋轉表達

式中Rwb為 旋轉矩陣，表示坐標系 Fw到 Fb的旋轉，滿足如下約束

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

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

式中CH(·)為旋轉四元數(shù)到旋轉矩陣的映射.對應式（3）的關系式為

2 Hamilton 與JPL 約定的比較

2.1 與旋轉矩陣乘法的同形性

選擇在Hamilton 約定基礎上，提出另外的四元數(shù)表達約定，原因在于部分學者認為旋轉四元數(shù)與旋轉矩陣的映射關系滿足[13]

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

也就是

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

因而

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

在JPL 約定下，滿足

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

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

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

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

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

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

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

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

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

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

表1 四元數(shù) 對應左乘和右乘矩陣形式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 約定，假設旋轉運動過程中t時刻的旋轉四元數(shù)為qˉbw(t)，可以得到

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

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

數(shù)的微分時，有

式中：

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 Hamilton 和JPL 約定的應用探討

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

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

4 仿真驗證

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

選取文獻[26]的航姿參考系統(tǒng)（attitude and heading reference system，AHRS）傳感器數(shù)據(jù)，分別利用Hamilton 和JPL 四元數(shù)對Mahony 算法效率和運算結果進行驗證.其中，傳感器數(shù)據(jù)是從AHRS 設備中記錄得到，經過標定校準的陀螺儀、加速度計和磁力計數(shù)據(jù)，運動過程是順序的繞X,Y,Z軸從0°轉動到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ù)與旋轉矩陣、姿態(tài)角和軸角之間的相互轉換關系進行測試，相同轉換運算下，兩者的用時分別是JPL 四元數(shù)0.011 038 s，Hamilton 四元數(shù)0.011 022 s，可見在忽略系統(tǒng)資源調度的情況下，兩者用時幾乎一致，主要源于在不同約定下，轉換約定的數(shù)學運算除去符號和順序等，均保持一致.

Mahony 算法采用Hamilton 四元數(shù)進行推導，使用JPL 四元數(shù)進行改寫時，需要依據(jù)表1 修改四元數(shù)乘法，采用式（17）的坐標系轉換關系，以及式（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ù)量的處理，存在細微的不同，但該時差在系統(tǒng)運行的誤差范圍內.這符合理論上的推導結果，由Hamilton 約定改為JPL 約定，只需要同時變更四元數(shù)虛部的符號，并且更改相應的運算公式，而運算轉換之間本身耗時基本一致，因此，兩種約定在所需計算資源上一致.

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

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

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

5 結 論

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