四元数与欧拉角刚体动力学数值积分算法及其比较

徐小明　钟万勰

摘要：为推广四元数保辛积分在工程中的应用，对欧拉角表示的状态方程数值积分与四元数的保辛积分进行比较.重陀螺的数值仿真结果表明四元数保辛积分的数值结果明显优于欧拉角状态方程积分.与欧拉角状态方程积分相比，四元数保辛积分在刚体动力学的数值仿真中更具优势.

关键词：四元数； 欧拉角； 刚体动力学； 保辛； 重陀螺

中图分类号： 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.

（编辑于杰）