四元數與歐拉角剛體動力學數值積分算法及其比較

徐小明　鐘萬勰

摘要：為推廣四元數保辛積分在工程中的應用，對歐拉角表示的狀態方程數值積分與四元數的保辛積分進行比較.重陀螺的數值仿真結果表明四元數保辛積分的數值結果明顯優于歐拉角狀態方程積分.與歐拉角狀態方程積分相比，四元數保辛積分在剛體動力學的數值仿真中更具優勢.

關鍵詞：四元數； 歐拉角； 剛體動力學； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻標志碼： A

0引言

四元數、歐拉角和方向余弦[1]是描述剛體旋轉最主要的3種坐標形式.方向余弦法需要9個參量，應用較少；而另外2種則應用廣泛，如航天器姿態控制和捷聯式慣性導航系統[1]等，對于兩者的研究也卷帙浩繁，但對兩者優劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態可采用2類12種歐拉角系統，分別對應于不同的旋轉軸先后次序.目前公認的用歐拉角描述旋轉的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統，都不可避免地會含有奇異點，使得在該點附近區域進行的數值積分精度不高.對于小角度旋轉，只要采用適當的歐拉角系統便可避開奇異點；而在大角度旋轉時，若想避開奇異點，必須提供2套歐拉角系統交替進行計算.據此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學者[4]在繼續研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復雜.

四元數用于描述剛體旋轉，沒有奇異性問題，可很好地描述剛體的全角度旋轉.然而，四元數需要滿足長度等于1的單位約束，這成為制約其應用的限制.在實際應用中，經常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數表示的剛體動力學方程導入微分代數方程范疇，近年來逐漸有學者[57]展開相關問題的研究.徐小明等[8]提出一種基于四元數理論描述剛體旋轉的保辛數值積分方法.該方法先將問題導入微分代數方程系統，然后利用分析結構力學理論[9]進行逐步積分，該積分保辛.該方法在積分點上嚴格滿足四元數長度等于1的約束條件，而在區段內部則由插值近似，屬于祖沖之類方法[10].數值算例表明仿真效果優異.

本文簡要介紹四元數和歐拉角的基本理論，以重陀螺為例對2種表示形式的數值積分進行比較.對于歐拉角表述，應用比較普遍的狀態方程表述.研究表明，以歐拉角和角速度為狀態變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關注；而采用文獻[8]給出的保辛格式，四元數單位長度約束條件得到滿足，仿真結果優異，能量也達到守恒.

1剛體旋轉及其運動學表示

剛體運動由質心平動和繞質心轉動組成.如果剛體上有一固定點，則稱為剛體定點轉動問題.假設Oxyz為系統的慣性坐標系，O為原點，亦為固定點.Ox′y′z′為隨體坐標系，固定于剛體上.若將剛體的隨體坐標軸取為與固定點有關的主軸，則剛體定點轉動可由式（1）描述.

2種數值積分的系統能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進行數值積分，但是能量卻保持得不好，這驗證對狀態方程應用辛歐拉格式并不能保辛.與之相反，四元數的數值積分保辛，其能量保持得很好，這也是保辛積分的優勢所在.四元數保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結束語

介紹2種剛體旋轉的數值積分，一種基于歐拉角表示，另一種基于四元數表示.以重陀螺的高速旋轉為例，對2種數值積分進行比較發現：以歐拉角、角速度組成狀態變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數值積分存在缺陷；與之相

反，采用四元數表示，根據分析結構力學的保辛積分方法，并結合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態變量的數值積分進行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻：

[1]張樹俠， 孫靜. 捷聯式慣性導航系統[M]. 北京： 國防工業出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數表示下姿態的2階運動奇異性[J]. 科學技術與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態角解算的全角度雙歐法[J]. 北京航空航天大學學報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學的四元數表示及保辛積分[J]. 應用數學和力學， 2014， 35（1）： 111.

[9]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數在剛體姿態仿真中的應用研究[J].飛行力學， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）

摘要：為推廣四元數保辛積分在工程中的應用，對歐拉角表示的狀態方程數值積分與四元數的保辛積分進行比較.重陀螺的數值仿真結果表明四元數保辛積分的數值結果明顯優于歐拉角狀態方程積分.與歐拉角狀態方程積分相比，四元數保辛積分在剛體動力學的數值仿真中更具優勢.

關鍵詞：四元數； 歐拉角； 剛體動力學； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻標志碼： A

0引言

四元數、歐拉角和方向余弦[1]是描述剛體旋轉最主要的3種坐標形式.方向余弦法需要9個參量，應用較少；而另外2種則應用廣泛，如航天器姿態控制和捷聯式慣性導航系統[1]等，對于兩者的研究也卷帙浩繁，但對兩者優劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態可采用2類12種歐拉角系統，分別對應于不同的旋轉軸先后次序.目前公認的用歐拉角描述旋轉的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統，都不可避免地會含有奇異點，使得在該點附近區域進行的數值積分精度不高.對于小角度旋轉，只要采用適當的歐拉角系統便可避開奇異點；而在大角度旋轉時，若想避開奇異點，必須提供2套歐拉角系統交替進行計算.據此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學者[4]在繼續研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復雜.

四元數用于描述剛體旋轉，沒有奇異性問題，可很好地描述剛體的全角度旋轉.然而，四元數需要滿足長度等于1的單位約束，這成為制約其應用的限制.在實際應用中，經常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數表示的剛體動力學方程導入微分代數方程范疇，近年來逐漸有學者[57]展開相關問題的研究.徐小明等[8]提出一種基于四元數理論描述剛體旋轉的保辛數值積分方法.該方法先將問題導入微分代數方程系統，然后利用分析結構力學理論[9]進行逐步積分，該積分保辛.該方法在積分點上嚴格滿足四元數長度等于1的約束條件，而在區段內部則由插值近似，屬于祖沖之類方法[10].數值算例表明仿真效果優異.

本文簡要介紹四元數和歐拉角的基本理論，以重陀螺為例對2種表示形式的數值積分進行比較.對于歐拉角表述，應用比較普遍的狀態方程表述.研究表明，以歐拉角和角速度為狀態變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關注；而采用文獻[8]給出的保辛格式，四元數單位長度約束條件得到滿足，仿真結果優異，能量也達到守恒.

1剛體旋轉及其運動學表示

剛體運動由質心平動和繞質心轉動組成.如果剛體上有一固定點，則稱為剛體定點轉動問題.假設Oxyz為系統的慣性坐標系，O為原點，亦為固定點.Ox′y′z′為隨體坐標系，固定于剛體上.若將剛體的隨體坐標軸取為與固定點有關的主軸，則剛體定點轉動可由式（1）描述.

2種數值積分的系統能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進行數值積分，但是能量卻保持得不好，這驗證對狀態方程應用辛歐拉格式并不能保辛.與之相反，四元數的數值積分保辛，其能量保持得很好，這也是保辛積分的優勢所在.四元數保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結束語

介紹2種剛體旋轉的數值積分，一種基于歐拉角表示，另一種基于四元數表示.以重陀螺的高速旋轉為例，對2種數值積分進行比較發現：以歐拉角、角速度組成狀態變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數值積分存在缺陷；與之相

反，采用四元數表示，根據分析結構力學的保辛積分方法，并結合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態變量的數值積分進行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻：

[1]張樹俠， 孫靜. 捷聯式慣性導航系統[M]. 北京： 國防工業出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數表示下姿態的2階運動奇異性[J]. 科學技術與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態角解算的全角度雙歐法[J]. 北京航空航天大學學報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學的四元數表示及保辛積分[J]. 應用數學和力學， 2014， 35（1）： 111.

[9]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數在剛體姿態仿真中的應用研究[J].飛行力學， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）

摘要：為推廣四元數保辛積分在工程中的應用，對歐拉角表示的狀態方程數值積分與四元數的保辛積分進行比較.重陀螺的數值仿真結果表明四元數保辛積分的數值結果明顯優于歐拉角狀態方程積分.與歐拉角狀態方程積分相比，四元數保辛積分在剛體動力學的數值仿真中更具優勢.

關鍵詞：四元數； 歐拉角； 剛體動力學； 保辛； 重陀螺

中圖分類號： TP391.9； O313.3

文獻標志碼： A

0引言

四元數、歐拉角和方向余弦[1]是描述剛體旋轉最主要的3種坐標形式.方向余弦法需要9個參量，應用較少；而另外2種則應用廣泛，如航天器姿態控制和捷聯式慣性導航系統[1]等，對于兩者的研究也卷帙浩繁，但對兩者優劣的評價卻褒貶不一.

描述剛體在三維空間中的運動姿態可采用2類12種歐拉角系統，分別對應于不同的旋轉軸先后次序.目前公認的用歐拉角描述旋轉的固有缺陷是奇異性問題[2]，即：無論采用哪種歐拉角系統，都不可避免地會含有奇異點，使得在該點附近區域進行的數值積分精度不高.對于小角度旋轉，只要采用適當的歐拉角系統便可避開奇異點；而在大角度旋轉時，若想避開奇異點，必須提供2套歐拉角系統交替進行計算.據此，黃雪樵[3]提出一種“雙歐法”的計算方法；近幾年仍有學者[4]在繼續研究該方法.雙歐法雖然解決奇異性問題，但是計算過程較復雜.

四元數用于描述剛體旋轉，沒有奇異性問題，可很好地描述剛體的全角度旋轉.然而，四元數需要滿足長度等于1的單位約束，這成為制約其應用的限制.在實際應用中，經常采用的正交化修正等方法只能緩解長度的偏移，無法從根本上解決問題；黃雪樵[3]在其雙歐法中也有所討論.目前，對于單位約束最佳的解決方案是將四元數表示的剛體動力學方程導入微分代數方程范疇，近年來逐漸有學者[57]展開相關問題的研究.徐小明等[8]提出一種基于四元數理論描述剛體旋轉的保辛數值積分方法.該方法先將問題導入微分代數方程系統，然后利用分析結構力學理論[9]進行逐步積分，該積分保辛.該方法在積分點上嚴格滿足四元數長度等于1的約束條件，而在區段內部則由插值近似，屬于祖沖之類方法[10].數值算例表明仿真效果優異.

本文簡要介紹四元數和歐拉角的基本理論，以重陀螺為例對2種表示形式的數值積分進行比較.對于歐拉角表述，應用比較普遍的狀態方程表述.研究表明，以歐拉角和角速度為狀態變量形成的1階微分方程在使用差分近似積分時，精度損失很大，能量不守恒；該現象被周江華等[11]稱為“睡陀螺”.這是一個值得注意的問題，卻未得到廣泛關注；而采用文獻[8]給出的保辛格式，四元數單位長度約束條件得到滿足，仿真結果優異，能量也達到守恒.

1剛體旋轉及其運動學表示

剛體運動由質心平動和繞質心轉動組成.如果剛體上有一固定點，則稱為剛體定點轉動問題.假設Oxyz為系統的慣性坐標系，O為原點，亦為固定點.Ox′y′z′為隨體坐標系，固定于剛體上.若將剛體的隨體坐標軸取為與固定點有關的主軸，則剛體定點轉動可由式（1）描述.

2種數值積分的系統能量隨時間變化情況見圖4.在歐拉角表示中，雖然采用辛歐拉格式進行數值積分，但是能量卻保持得不好，這驗證對狀態方程應用辛歐拉格式并不能保辛.與之相反，四元數的數值積分保辛，其能量保持得很好，這也是保辛積分的優勢所在.四元數保辛積分的約束誤差見圖5，表明在時間積分過程中單位長度約束條件滿足得很好.

4結束語

介紹2種剛體旋轉的數值積分，一種基于歐拉角表示，另一種基于四元數表示.以重陀螺的高速旋轉為例，對2種數值積分進行比較發現：以歐拉角、角速度組成狀態變量，然后直接使用辛歐拉格式不

能保辛，能量衰減很快，數值積分存在缺陷；與之相

反，采用四元數表示，根據分析結構力學的保辛積分方法，并結合祖沖之類方法的思想，可以很好地避免約束違約，仿真效果令人滿意，可作為陀螺等仿真分析的有力工具.

本文僅對以歐拉角、角速度組成狀態變量的數值積分進行研究，對其他形式并未涉及，對其積分效果不佳的成因亦未研究，還有很多工作有待展開.

參考文獻：

[1]張樹俠， 孫靜. 捷聯式慣性導航系統[M]. 北京： 國防工業出版社， 1992： 4880.

[2]趙曉穎， 溫立書， 么彩蓮. 歐拉角參數表示下姿態的2階運動奇異性[J]. 科學技術與工程， 2012， 12（3）： 634637.

[3]黃雪樵. 克服歐拉方程奇異性的雙歐法[J]. 飛行力學， 1994， 12（4）： 2837.

[4]李躍軍， 閻超. 飛行器姿態角解算的全角度雙歐法[J]. 北京航空航天大學學報， 2007， 33（5）： 505508.

[5]NIKRAVESH P E， WEHAGE R A， KWON K. Euler parameters in computational kinematics and dynamics： Part 1[J]. J Mechanisms， Transmissions & Automation Des， 1985， 107（3）： 358365.

[6]BETSCH P， SIEBERT R. Rigid body dynamics in terms of quaternions： Hamiltonian formulation and conserving numerical integration[J]. Int J Numer Methods Eng， 2009， 79（4）： 444473.

[7]UDWADIA F E， SCHUTTE A D. An alternative derivation of the quaternion equations of motion for rigidbody rotational dynamics[J]. J Appl Mech， 2010， 77（4）： 44505.

[8]徐小明， 鐘萬勰. 剛體動力學的四元數表示及保辛積分[J]. 應用數學和力學， 2014， 35（1）： 111.

[9]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[10]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數在剛體姿態仿真中的應用研究[J].飛行力學， 2000， 18（4）： 2832.

[12]HAIRER E， LUBICH C， WANNER G. Geometric numerical integration： structurepreserving algorithms for ordinary differential equations[M]. Berlin： Springer， 2006： 189.

（編輯于杰）