四元數(shù)與歐拉角剛體動(dòng)力學(xué)數(shù)值積分算法及其比較

徐小明　鐘萬(wàn)勰

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對(duì)歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動(dòng)力學(xué)的數(shù)值仿真中更具優(yōu)勢(shì).

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動(dòng)力學(xué)； 保辛； 重陀螺

中圖分類(lèi)號(hào)： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個(gè)參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對(duì)于兩者的研究也卷帙浩繁，但對(duì)兩者優(yōu)劣的評(píng)價(jià)卻褒貶不一.

描述剛體在三維空間中的運(yùn)動(dòng)姿態(tài)可采用2類(lèi)12種歐拉角系統(tǒng)，分別對(duì)應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問(wèn)題[2]，即：無(wú)論采用哪種歐拉角系統(tǒng)，都不可避免地會(huì)含有奇異點(diǎn)，使得在該點(diǎn)附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對(duì)于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開(kāi)奇異點(diǎn)；而在大角度旋轉(zhuǎn)時(shí)，若想避開(kāi)奇異點(diǎn)，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計(jì)算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計(jì)算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問(wèn)題，但是計(jì)算過(guò)程較復(fù)雜.

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

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

1剛體旋轉(zhuǎn)及其運(yùn)動(dòng)學(xué)表示

剛體運(yùn)動(dòng)由質(zhì)心平動(dòng)和繞質(zhì)心轉(zhuǎn)動(dòng)組成.如果剛體上有一固定點(diǎn)，則稱(chēng)為剛體定點(diǎn)轉(zhuǎn)動(dòng)問(wèn)題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點(diǎn)，亦為固定點(diǎn).Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點(diǎn)有關(guān)的主軸，則剛體定點(diǎn)轉(zhuǎn)動(dòng)可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時(shí)間變化情況見(jiàn)圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗(yàn)證對(duì)狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢(shì)所在.四元數(shù)保辛積分的約束誤差見(jiàn)圖5，表明在時(shí)間積分過(guò)程中單位長(zhǎng)度約束條件滿(mǎn)足得很好.

4結(jié)束語(yǔ)

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對(duì)2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

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

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類(lèi)方法的思想，可以很好地避免約束違約，仿真效果令人滿(mǎn)意，可作為陀螺等仿真分析的有力工具.

本文僅對(duì)以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對(duì)其他形式并未涉及，對(duì)其積分效果不佳的成因亦未研究，還有很多工作有待展開(kāi).

參考文獻(xiàn)：

[1]張樹(shù)俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國(guó)防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書(shū)， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運(yùn)動(dòng)奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

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

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報(bào)， 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]徐小明， 鐘萬(wàn)勰. 剛體動(dòng)力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬(wàn)勰， 高強(qiáng). 約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動(dòng)力學(xué)與控制學(xué)報(bào)， 2006， 4（3）： 193200.

[10]鐘萬(wàn)勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 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.

（編輯于杰）

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對(duì)歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動(dòng)力學(xué)的數(shù)值仿真中更具優(yōu)勢(shì).

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動(dòng)力學(xué)； 保辛； 重陀螺

中圖分類(lèi)號(hào)： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個(gè)參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對(duì)于兩者的研究也卷帙浩繁，但對(duì)兩者優(yōu)劣的評(píng)價(jià)卻褒貶不一.

描述剛體在三維空間中的運(yùn)動(dòng)姿態(tài)可采用2類(lèi)12種歐拉角系統(tǒng)，分別對(duì)應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問(wèn)題[2]，即：無(wú)論采用哪種歐拉角系統(tǒng)，都不可避免地會(huì)含有奇異點(diǎn)，使得在該點(diǎn)附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對(duì)于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開(kāi)奇異點(diǎn)；而在大角度旋轉(zhuǎn)時(shí)，若想避開(kāi)奇異點(diǎn)，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計(jì)算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計(jì)算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問(wèn)題，但是計(jì)算過(guò)程較復(fù)雜.

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

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

1剛體旋轉(zhuǎn)及其運(yùn)動(dòng)學(xué)表示

剛體運(yùn)動(dòng)由質(zhì)心平動(dòng)和繞質(zhì)心轉(zhuǎn)動(dòng)組成.如果剛體上有一固定點(diǎn)，則稱(chēng)為剛體定點(diǎn)轉(zhuǎn)動(dòng)問(wèn)題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點(diǎn)，亦為固定點(diǎn).Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點(diǎn)有關(guān)的主軸，則剛體定點(diǎn)轉(zhuǎn)動(dòng)可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時(shí)間變化情況見(jiàn)圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗(yàn)證對(duì)狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢(shì)所在.四元數(shù)保辛積分的約束誤差見(jiàn)圖5，表明在時(shí)間積分過(guò)程中單位長(zhǎng)度約束條件滿(mǎn)足得很好.

4結(jié)束語(yǔ)

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對(duì)2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

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

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類(lèi)方法的思想，可以很好地避免約束違約，仿真效果令人滿(mǎn)意，可作為陀螺等仿真分析的有力工具.

本文僅對(duì)以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對(duì)其他形式并未涉及，對(duì)其積分效果不佳的成因亦未研究，還有很多工作有待展開(kāi).

參考文獻(xiàn)：

[1]張樹(shù)俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國(guó)防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書(shū)， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運(yùn)動(dòng)奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

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

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報(bào)， 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]徐小明， 鐘萬(wàn)勰. 剛體動(dòng)力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬(wàn)勰， 高強(qiáng). 約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動(dòng)力學(xué)與控制學(xué)報(bào)， 2006， 4（3）： 193200.

[10]鐘萬(wàn)勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 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.

（編輯于杰）

摘要：為推廣四元數(shù)保辛積分在工程中的應(yīng)用，對(duì)歐拉角表示的狀態(tài)方程數(shù)值積分與四元數(shù)的保辛積分進(jìn)行比較.重陀螺的數(shù)值仿真結(jié)果表明四元數(shù)保辛積分的數(shù)值結(jié)果明顯優(yōu)于歐拉角狀態(tài)方程積分.與歐拉角狀態(tài)方程積分相比，四元數(shù)保辛積分在剛體動(dòng)力學(xué)的數(shù)值仿真中更具優(yōu)勢(shì).

關(guān)鍵詞：四元數(shù)； 歐拉角； 剛體動(dòng)力學(xué)； 保辛； 重陀螺

中圖分類(lèi)號(hào)： TP391.9； O313.3

文獻(xiàn)標(biāo)志碼： A

0引言

四元數(shù)、歐拉角和方向余弦[1]是描述剛體旋轉(zhuǎn)最主要的3種坐標(biāo)形式.方向余弦法需要9個(gè)參量，應(yīng)用較少；而另外2種則應(yīng)用廣泛，如航天器姿態(tài)控制和捷聯(lián)式慣性導(dǎo)航系統(tǒng)[1]等，對(duì)于兩者的研究也卷帙浩繁，但對(duì)兩者優(yōu)劣的評(píng)價(jià)卻褒貶不一.

描述剛體在三維空間中的運(yùn)動(dòng)姿態(tài)可采用2類(lèi)12種歐拉角系統(tǒng)，分別對(duì)應(yīng)于不同的旋轉(zhuǎn)軸先后次序.目前公認(rèn)的用歐拉角描述旋轉(zhuǎn)的固有缺陷是奇異性問(wèn)題[2]，即：無(wú)論采用哪種歐拉角系統(tǒng)，都不可避免地會(huì)含有奇異點(diǎn)，使得在該點(diǎn)附近區(qū)域進(jìn)行的數(shù)值積分精度不高.對(duì)于小角度旋轉(zhuǎn)，只要采用適當(dāng)?shù)臍W拉角系統(tǒng)便可避開(kāi)奇異點(diǎn)；而在大角度旋轉(zhuǎn)時(shí)，若想避開(kāi)奇異點(diǎn)，必須提供2套歐拉角系統(tǒng)交替進(jìn)行計(jì)算.據(jù)此，黃雪樵[3]提出一種“雙歐法”的計(jì)算方法；近幾年仍有學(xué)者[4]在繼續(xù)研究該方法.雙歐法雖然解決奇異性問(wèn)題，但是計(jì)算過(guò)程較復(fù)雜.

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

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

1剛體旋轉(zhuǎn)及其運(yùn)動(dòng)學(xué)表示

剛體運(yùn)動(dòng)由質(zhì)心平動(dòng)和繞質(zhì)心轉(zhuǎn)動(dòng)組成.如果剛體上有一固定點(diǎn)，則稱(chēng)為剛體定點(diǎn)轉(zhuǎn)動(dòng)問(wèn)題.假設(shè)Oxyz為系統(tǒng)的慣性坐標(biāo)系，O為原點(diǎn)，亦為固定點(diǎn).Ox′y′z′為隨體坐標(biāo)系，固定于剛體上.若將剛體的隨體坐標(biāo)軸取為與固定點(diǎn)有關(guān)的主軸，則剛體定點(diǎn)轉(zhuǎn)動(dòng)可由式（1）描述.

2種數(shù)值積分的系統(tǒng)能量隨時(shí)間變化情況見(jiàn)圖4.在歐拉角表示中，雖然采用辛歐拉格式進(jìn)行數(shù)值積分，但是能量卻保持得不好，這驗(yàn)證對(duì)狀態(tài)方程應(yīng)用辛歐拉格式并不能保辛.與之相反，四元數(shù)的數(shù)值積分保辛，其能量保持得很好，這也是保辛積分的優(yōu)勢(shì)所在.四元數(shù)保辛積分的約束誤差見(jiàn)圖5，表明在時(shí)間積分過(guò)程中單位長(zhǎng)度約束條件滿(mǎn)足得很好.

4結(jié)束語(yǔ)

介紹2種剛體旋轉(zhuǎn)的數(shù)值積分，一種基于歐拉角表示，另一種基于四元數(shù)表示.以重陀螺的高速旋轉(zhuǎn)為例，對(duì)2種數(shù)值積分進(jìn)行比較發(fā)現(xiàn)：以歐拉角、角速度組成狀態(tài)變量，然后直接使用辛歐拉格式不

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

反，采用四元數(shù)表示，根據(jù)分析結(jié)構(gòu)力學(xué)的保辛積分方法，并結(jié)合祖沖之類(lèi)方法的思想，可以很好地避免約束違約，仿真效果令人滿(mǎn)意，可作為陀螺等仿真分析的有力工具.

本文僅對(duì)以歐拉角、角速度組成狀態(tài)變量的數(shù)值積分進(jìn)行研究，對(duì)其他形式并未涉及，對(duì)其積分效果不佳的成因亦未研究，還有很多工作有待展開(kāi).

參考文獻(xiàn)：

[1]張樹(shù)俠， 孫靜. 捷聯(lián)式慣性導(dǎo)航系統(tǒng)[M]. 北京： 國(guó)防工業(yè)出版社， 1992： 4880.

[2]趙曉穎， 溫立書(shū)， 么彩蓮. 歐拉角參數(shù)表示下姿態(tài)的2階運(yùn)動(dòng)奇異性[J]. 科學(xué)技術(shù)與工程， 2012， 12（3）： 634637.

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

[4]李躍軍， 閻超. 飛行器姿態(tài)角解算的全角度雙歐法[J]. 北京航空航天大學(xué)學(xué)報(bào)， 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]徐小明， 鐘萬(wàn)勰. 剛體動(dòng)力學(xué)的四元數(shù)表示及保辛積分[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2014， 35（1）： 111.

[9]鐘萬(wàn)勰， 高強(qiáng). 約束動(dòng)力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動(dòng)力學(xué)與控制學(xué)報(bào)， 2006， 4（3）： 193200.

[10]鐘萬(wàn)勰， 高強(qiáng)， 彭海軍. 經(jīng)典力學(xué)辛講[M]. 大連： 大連理工大學(xué)出版社， 2013： 202241.

[11]周江華， 苗育紅， 李宏， 等. 四元數(shù)在剛體姿態(tài)仿真中的應(yīng)用研究[J].飛行力學(xué)， 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.

（編輯于杰）