微分转动与角速度：三维空间中的矩阵向量形式及其Python实现

引言

在刚体动力学和三维运动学中，微分转动与角速度构成了描述旋转运动的核心框架。这两个概念不仅在理论物理中具有基础性地位，更在机器人学、航空航天控制和计算机图形学等领域有着广泛应用。本文将从矩阵向量形式出发，深入探讨微分转动与角速度的数学本质，并通过Python的SymPy和NumPy库展示其计算与实现过程。

理论基础

旋转矩阵与微分转动

旋转矩阵$ \mathbf{R} $是描述刚体方向的正交变换矩阵，满足$ \mathbf{R}^T \mathbf{R} = \mathbf{I} $且$ \det(\mathbf{R}) = 1 $。当考虑时间$ dt $内的微小旋转时，我们得到\ast \ast 微分转动\ast \ast$ d\mathbf{R} $。由正交性约束导出的关键关系为：

$$(\mathbf{R}+d\mathbf{R}{)}^T(\mathbf{R}+d\mathbf{R})=\mathbf{I}$$

展开并忽略高阶项后，得到：

$${\mathbf{R}}^{T}d\mathbf{R}+(d\mathbf{R}{)}^T\mathbf{R}=0$$

这表明$ \mathbf{R}^T d\mathbf{R} $是一个\ast \ast 反对称矩阵\ast \ast ，记作$ \boldsymbol{\Omega} $：

$$\mathbf{\Omega }\triangleq {\mathbf{R}}^{T}d\mathbf{R}$$

满足$ \boldsymbol{\Omega} + \boldsymbol{\Omega}^T = \mathbf{0} $。

角速度的反对称表示

三维角速度向量$ \boldsymbol{\omega} = [\omega_x, \omega_y, \omega_z]^T $可表示为反对称矩阵：

$$[\mathbf{\omega }{]}_{\times }=\left[\begin{array}{ccc}0& -{\omega }_z& {\omega }_y\\ {\omega }_z& 0& -{\omega }_x\\ -{\omega }_y& {\omega }_x& 0\end{array}\right]$$

微分转动与角速度的基本关系为：

$$\mathbf{\Omega }=[{\mathbf{\omega }}_b{]}_{\times }dt$$

其中$ \boldsymbol{\omega}_b $是体坐标系下的角速度。

旋转矩阵微分方程

在体坐标系下，角速度与旋转矩阵的变化率满足：

$$d\mathbf{R}=\mathbf{R}[{\mathbf{\omega }}_b{]}_{\times }dt$$

或等价地：

$$\frac{d\mathbf{R}}{dt}=\mathbf{R}[{\mathbf{\omega }}_b{]}_{\times }$$

在空间坐标系（固定参考系）下，角速度$ \boldsymbol{\omega}_s $满足：

$$\frac{d\mathbf{R}}{dt}=[{\mathbf{\omega }}_s{]}_{\times }\mathbf{R}$$

两种坐标系间的转换关系为：

$${\mathbf{\omega }}_s=\mathbf{R}{\mathbf{\omega }}_b$$
$$[{\mathbf{\omega }}_s{]}_{\times }=\mathbf{R}[{\mathbf{\omega }}_b{]}_{\times }{\mathbf{R}}^T$$

线速度与角速度

对固定在刚体上的点$ \mathbf{p} $（体坐标系中为常向量），其空间坐标为$ \mathbf{r} = \mathbf{R} \mathbf{p} $。线速度$ \mathbf{v} $为：

$$\mathbf{v}=\frac{d\mathbf{r}}{dt}=\frac{d\mathbf{R}}{dt}\mathbf{p}$$

代入角速度表达式得：

$$\mathbf{v}=\mathbf{R}({\mathbf{\omega }}_b\times \mathbf{p})={\mathbf{\omega }}_s\times \mathbf{r}$$

符号推导（SymPy实现）

SymPy库提供了强大的符号计算能力，可严格验证上述关系：

import sympy as sp
from sympy import Matrix, symbols, eye, sin, cos, diff, simplify# 定义符号变量
t = symbols('t')
theta = symbols('theta', cls=sp.Function)(t)
omega_x, omega_y, omega_z = symbols('omega_x omega_y omega_z')# 角速度向量和反对称矩阵
omega = Matrix([omega_x, omega_y, omega_z])
omega_skew = Matrix([[0, -omega_z, omega_y],[omega_z, 0, -omega_x],[-omega_y, omega_x, 0]
])# 绕x轴旋转矩阵
R_x = Matrix([[1, 0, 0],[0, cos(theta), -sin(theta)],[0, sin(theta), cos(theta)]
])# 验证微分方程
dR_dt = diff(R_x, t)
omega_actual = diff(theta, t)
omega_b = Matrix([omega_actual, 0, 0])
omega_b_skew = Matrix([[0, 0, 0],[0, 0, -omega_actual],[0, omega_actual, 0]
])
dR_dt_calculated = R_x @ omega_b_skew
print("微分方程验证:", simplify(dR_dt - dR_dt_calculated) == Matrix.zeros(3, 3))# 验证坐标系转换
R = Matrix([[sp.cos(theta), -sp.sin(theta), 0],[sp.sin(theta), sp.cos(theta), 0],[0, 0, 1]
])
dR_dt = diff(R, t)
omega_s_skew = dR_dt @ R.T
omega_s = Matrix([omega_s_skew[2,1], omega_s_skew[0,2], omega_s_skew[1,0]])
omega_b = Matrix([0, 0, diff(theta, t)])
print("坐标转换验证:", simplify(omega_s - R @ omega_b) == Matrix.zeros(3, 1))

数值计算（NumPy实现）

NumPy结合SciPy提供高效的数值计算能力：

import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D# 反对称矩阵函数
def skew(omega):return np.array([[0, -omega[2], omega[1]],[omega[2], 0, -omega[0]],[-omega[1], omega[0], 0]])# 基本旋转计算
theta = np.pi/4
R_z = np.array([[np.cos(theta), -np.sin(theta), 0],[np.sin(theta), np.cos(theta), 0],[0, 0, 1]
])# 微分转动计算
omega_b = np.array([0, 0, 1.0])
dt = 0.01
dR = R_z @ skew(omega_b) * dt
R_new = R_z + dR# 指数映射精确解
R_exact = R_z @ expm(skew(omega_b) * dt)# 运动学计算
p_body = np.array([1, 0, 0])
r_space = R_z @ p_body
v_space = np.cross(omega_b, r_space)  # v = ω × r
omega_s = R_z @ omega_b# 轨迹可视化
times = np.linspace(0, 2*np.pi, 100)
positions = []
for t in times:R_t = np.array([[np.cos(t), -np.sin(t), 0],[np.sin(t), np.cos(t), 0],[0, 0, 1]])positions.append(R_t @ p_body)
positions = np.array(positions)fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(positions[:,0], positions[:,1], positions[:,2], 'b-', linewidth=2)
ax.quiver(0, 0, 0, 0, 0, 1.2, color='r', arrow_length_ratio=0.1)
ax.scatter(positions[::10,0], positions[::10,1], positions[::10,2], c='g', s=50)
ax.set_xlabel('X', fontsize=12)
ax.set_ylabel('Y', fontsize=12)
ax.set_zlabel('Z', fontsize=12)
ax.set_title('刚体点圆周运动轨迹', fontsize=14)
ax.grid(True)
plt.tight_layout()
plt.show()

应用场景

机器人运动学

在机械臂控制中，末端执行器的速度通过雅可比矩阵与关节速度关联：

$${\mathbf{v}}_{ee}=\mathbf{J}\dot{\mathbf{q}}$$

其中旋转部分雅可比直接与角速度相关：

# 六轴机械臂雅可比矩阵示例
J = np.zeros((6, 6)) 
# ... 填充位置和旋转部分雅可比 ...
joint_velocities = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6])
end_effector_velocity = J @ joint_velocities
angular_velocity = end_effector_velocity[3:]

航天器姿态动力学

四元数常用于航天器姿态表示，其微分方程与角速度直接相关：

def quat_derivative(q, omega):q0, q1, q2, q3 = qreturn 0.5 * np.array([[-q1, -q2, -q3],[q0, -q3, q2],[q3, q0, -q1],[-q2, q1, q0]]) @ omega

计算机图形学旋转插值

球面线性插值(SLERP)在动画中提供平滑旋转过渡：

def slerp(q1, q2, t):dot = np.dot(q1, q2)theta = np.arccos(np.clip(dot, -1.0, 1.0))sin_theta = np.sin(theta)return (np.sin((1-t)*theta)*q1 + np.sin(t*theta)*q2) / sin_theta

物理意义与几何解释

角速度$ \boldsymbol{\omega} $本质上是\ast \ast 瞬时旋转轴\ast \ast 方向与\ast \ast 旋转速率\ast \ast 的组合。反对称矩阵$ [\boldsymbol{\omega}]_{\times} $的零空间即为瞬时转轴方向，其秩为2，反映了旋转的自由度特性。

微分转动$ d\boldsymbol{\theta} = \boldsymbol{\omega} dt $仅在$ dt \to 0 $时具有矢量特性。有限旋转的不可交换性导致角速度不能直接积分得到有限转角，而需通过指数映射实现：

$$\mathbf{R}=\exp ([\mathbf{\theta }{]}_{\times })$$

其中$ \boldsymbol{\theta} = \boldsymbol{\omega} t $仅在角速度恒定时成立。这种非线性特性是旋转群$ SO(3) $流形结构的直接体现。

数值稳定性与实现技巧

实际应用中需注意：

1. 旋转矩阵正交化：数值积分可能导致矩阵失去正交性

def normalize_rotation(R):U, _, Vt = np.linalg.svd(R)return U @ Vt

2. 角速度坐标系的明智选择：
   • 体坐标系$ \boldsymbol{\omega}_b $：惯量张量恒定，适合动力学方程
   • 空间坐标系$ \boldsymbol{\omega}_s $：物理直观，适合运动学分析
3. 奇点处理：欧拉角在特定姿态存在奇点，四元数或旋转矩阵更鲁棒

结论

微分转动与角速度构成了三维旋转运动描述的双重基础：微分转动从几何上描述了姿态的无穷小变化，而角速度从物理上量化了旋转的瞬时状态。矩阵向量形式不仅提供了优雅的数学框架，更直接引导出高效的数值实现。

通过SymPy的符号推导能力，我们可严格验证理论关系；借助NumPy的数值计算功能，我们能在实际应用中高效求解。这种从理论到实践的闭环理解，对机器人控制、航天器导航和计算机动画等领域具有根本重要性。

旋转运动描述的本质挑战源于$ SO(3) $群的非线性特性——旋转空间不是平坦的欧几里得空间，而是具有曲率的黎曼流形。这一深层数学结构解释了为何有限旋转不可交换，也指引着我们继续探索更先进的姿态表示方法，如四元数、旋转向量和李群理论，它们都在不同角度上拓展了本文所述的基本框架。