摘要 首先推导角速度公式和角加速度公式，并举一个例子说明角速度公式的细节，然后从加速度公式出发推导转动惯量矩阵 J 和姿态动力学方程 Jw’=w×Jw，并解释为什么需要进行一次坐标变换。最后总结从力矩到姿态四元数的完整过程，并对无力矩输入时的进动现象进行仿真。

旋转坐标系下的速度和加速度

  本体系绕世界系以角速度 $\vec{\omega }$ 旋转，设世界系下位置向量 $\vec{r}$ 的一阶导和二阶导分别为 ${}^w\dot{\vec{r}}$ 和 ${}^w\ddot{\vec{r}}$，本体系下位置向量的一阶导和二阶导分别为 ${}^b\dot{\vec{r}}$ 和 ${}^b\ddot{\vec{r}}$，满足
$$\begin{array}{rl}{}^w\vec{r}=& R{}^b\vec{r}\\ {}^w\dot{\vec{r}}=& R{}^b\dot{\vec{r}}+\vec{\omega }\times {}^w\vec{r}\\ {}^w\ddot{\vec{r}}=& R{}^b\ddot{\vec{r}}+2\vec{\omega }\times R{}^b\dot{\vec{r}}+\vec{\omega }\times (\vec{\omega }\times {}^w\vec{r})+\dot{\vec{\omega }}\times {}^w\vec{r}\end{array}\text{\hspace{1em}(1)}$$
式中所有的 $R$ 都是从本体系到世界系的旋转矩阵，完整应写作 ${}_b^{w}R$；所有的 $\omega$ 都在世界系下表示，完整应写作 ${}^w\omega$：
$$\begin{array}{rl}{}^w\vec{r}=& {}_b^{w}R{}^b\vec{r}\\ {}^w\dot{\vec{r}}=& {}_b^{w}R{}^b\dot{\vec{r}}+{}^w\vec{\omega }\times {}^w\vec{r}\\ {}^w\ddot{\vec{r}}=& {}_b^{w}R{}^b\ddot{\vec{r}}+2{}^w\vec{\omega }\times {}_b^{w}R{}^b\dot{\vec{r}}+{}^w\vec{\omega }\times ({}^w\vec{\omega }\times {}^w\vec{r})+{}^w\dot{\vec{\omega }}\times {}^w\vec{r}\end{array}$$
证明：
$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}{(}^w\vec{r})=& \frac{\text{d}}{\text{d}t}(R{}^b\vec{r})\\ {}^w\dot{\vec{r}}=& \dot{R}{}^b\vec{r}+R{}^b\dot{\vec{r}}\\ =& R{}^b\dot{\vec{r}}+{}^w\vec{\omega }\times R{}^b\vec{r}\\ =& R{}^b\dot{\vec{r}}+{}^w\vec{\omega }\times {}^w\vec{r}\end{array}$$
此处以及后面的 ${}^w\vec{\omega }$ 都是在世界系下的坐标，下面简记作 $\vec{\omega }$。
  此处补充一些个人理解。$\dot{R}=\vec{\omega }\times R$ 是因为 $R$ 可以简单地看作三个向量绕旋转轴 $\vec{\omega }$ 转动，不涉及坐标变换；而 ${}^w\dot{\vec{r}}\ne \vec{\omega }\times {}^w\vec{r}$ 是因为此时向量 $\vec{r}$ 不是绕旋转轴 $\vec{\omega }$ 转动，而是个复合运动，$\vec{\omega }$ 是本体系（旋转坐标系）相对于世界系的旋转角速度，$\vec{r}$ 在本体系下还有 ${\vec{r}}_2$ 的角速度，当 $\vec{r}$ 在本体系下的角速度为0时两种情况等价，可以理解成有个本体系固定在 $\vec{r}$ 上。
  对式(1)的第二行再次求导
$$\begin{array}{rl}\frac{\text{d}}{\text{d}t}{(}^w\dot{\vec{r}})=& \frac{\text{d}}{\text{d}t}(R{}^b\dot{\vec{r}}+\vec{\omega }\times {}^w\vec{r})\\ {}^w\ddot{\vec{r}}=& \dot{R}{}^b\dot{\vec{r}}+R{}^b\ddot{\vec{r}}+\dot{\vec{\omega }}\times {}^w\vec{r}+\vec{\omega }\times {}^w\dot{\vec{r}}\\ =& R{}^b\ddot{\vec{r}}+\vec{\omega }\times R{}^b\dot{\vec{r}}+\vec{\omega }\times (R{}^b\dot{\vec{r}}+\vec{\omega }\times {}^w\vec{r})+\dot{\vec{\omega }}\times {}^w\vec{r}\\ =& R{}^b\ddot{\vec{r}}+2\vec{\omega }\times R{}^b\dot{\vec{r}}+\vec{\omega }\times (\vec{\omega }\times {}^w\vec{r})+\dot{\vec{\omega }}\times {}^w\vec{r}\end{array}$$
当匀速转动时 $\dot{\vec{\omega }}=0$，最后一项可以省略。

举例

  如图所示，世界系用 $O{x}_w{y}_w{z}_w$ 表示，本体系用 $O{x}_b{y}_b{z}_b$ 表示。一开始时，本体系与世界系重合，向量 $\vec{r}$ 与y轴正向重合，然后本体系沿世界系x轴逆时针旋转，角速度为 ${}^w{\vec{\omega }}_1$；向量 $\vec{r}$ 沿本体系z轴顺时针旋转，角速度为 ${}^b{\vec{\omega }}_2$
$${}^w{\vec{\omega }}_1=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],{}^b{\vec{\omega }}_2=\left[\begin{array}{c}0\\ 0\\ -1\end{array}\right]$$
注意两个角速度向量分别在世界系和本体系下表示。向量 $\vec{r}$ 在本体系下的坐标为
$${}^b\vec{r}=\left[\begin{array}{c}\sin t\\ \cos t\\ 0\end{array}\right]$$
本体系到世界系的旋转矩阵 ${}_b^{w}R$ 为
$${}_b^{w}R=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos t& -\sin t\\ 0& \sin t& \cos t\end{array}\right],{}_w^{b}R=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos t& \sin t\\ 0& -\sin t& \cos t\end{array}\right]$$
向量 $\vec{r}$ 在世界系下的坐标为
$${}^w\vec{r}={}_b^{w}R{}^b\vec{r}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos t& -\sin t\\ 0& \sin t& \cos t\end{array}\right]\left[\begin{array}{c}\sin t\\ \cos t\\ 0\end{array}\right]=\left[\begin{array}{c}\sin t\\ {\cos }^{2}t\\ \sin t\cos t\end{array}\right]$$
两边同时求导得到角速度公式
$$\begin{array}{rl}{}^w\vec{v}={}^w\dot{\vec{r}}=& {}_b^{w}R{}^b\dot{\vec{r}}+{}^w{\vec{\omega }}_1\times {}^w\vec{r}\\ \left[\begin{array}{c}\cos t\\ -\sin 2t\\ \cos 2t\end{array}\right]=& \left[\begin{array}{ccc}1& 0& 0\\ 0& \cos t& -\sin t\\ 0& \sin t& \cos t\end{array}\right]\left[\begin{array}{c}\cos t\\ -\sin t\\ 0\end{array}\right]+\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]\times \left[\begin{array}{c}\sin t\\ {\cos }^{2}t\\ \sin t\cos t\end{array}\right]\\ =& \left[\begin{array}{c}\cos t\\ -\sin t\cos t\\ -{\sin }^{2}t\end{array}\right]+\left[\begin{array}{c}0\\ -\sin t\cos t\\ {\cos }^{2}t\end{array}\right]\end{array}$$
于是可以得到，世界系下的速度向量转换到本体系下并不等于本体系下位置向量的导数。
$$\begin{array}{rl}& {}^b\vec{v}={}_w^{b}R{}^w\vec{v}=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos t& \sin t\\ 0& -\sin t& \cos t\end{array}\right]\left[\begin{array}{c}\cos t\\ -\sin 2t\\ \cos 2t\end{array}\right]=\left[\begin{array}{c}\cos t\\ -\sin t\\ \cos t\end{array}\right]\\ \ne & {}^b\dot{\vec{r}}=\left[\begin{array}{c}\cos t\\ -\sin t\\ 0\end{array}\right]\end{array}$$

姿态动力学方程

$$\vec{M}=J\dot{\vec{\omega }}+\vec{\omega }\times J\vec{\omega }$$
  参考资料[3]的第2章是我觉得所有资料里讲这个方程讲的最清楚的，网上看什么资料都不如看这本英文书，这篇博客就当是我的笔记了，强烈建议有条件还是看原文。我把原文大致翻译了一下:
刚体姿态动力学详细推导与翻译
  角动量定理
$$\frac{\text{d}\vec{H}}{\text{d}t}=\vec{M}$$
其中 $H,M$ 分别为刚体的角动量和作用在刚体上的力矩。

推导1

  下面的推导来自文末参考资料[1]，也是主流推导方法，主要思路是因为转动惯量矩阵在本体系保持不变，利用向量在本体系与世界系中的关系推导。但[1]不够详细，更详细的推导见参考资料[3]。
  刚体角动量是刚体内所有质量微元相对基准点 $S$ 的角动量之和
$$\vec{h}=m\vec{r}\times \vec{v},H=\underset{V}{\iiint }\vec{r}\times \dot{\vec{r}}\text{d}m$$
其中 $\text{d}m$ 为刚体的质量微元。
$$\begin{array}{rl}& \vec{r}\times \vec{v}=\vec{r}\times (\vec{\omega }\times \vec{r})=r^2\vec{\omega }-(\vec{r}\cdot \vec{\omega })\vec{r}\\ =& [(x^2+y^2+z^2)w_x-(x{w}_x+y{w}_y+z{w}_z)x]{\vec{e}}_x\\ +& [(x^2+y^2+z^2)w_y-(x{w}_x+y{w}_y+z{w}_z)y]{\vec{e}}_y\\ +& [(x^2+y^2+z^2)w_z-(x{w}_x+y{w}_y+z{w}_z)z]{\vec{e}}_z\\ =& [(y^2+z^2)w_x-xy{w}_y-xz{w}_z]{\vec{e}}_x\\ +& [-xy{w}_x+(x^2+z^2)w_y-yz{w}_z]{\vec{e}}_y\\ +& [-xz{w}_x-yz{w}_y+(x^2+y^2)w_z]{\vec{e}}_z\\ =& \left[\begin{array}{ccc}{y}^2+z^2& -xy& -xz\\ -xy& x^2+z^2& -yz\\ -xz& -yz& x^2+y^2\end{array}\right]\left[\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right]\end{array}$$
于是
$$H=\underset{V}{\iiint }\vec{r}\times \vec{v}\text{d}m=J\vec{\omega }$$
其中 $J$ 为刚体的转动惯量矩阵，形式为
$$J=\left[\begin{array}{ccc}{J}_x& -J_{xy}& -J_{xz}\\ -J_{xy}& J_y& -J_{yz}\\ -J_{xz}& -J_{yz}& J_z\end{array}\right]$$
其中
$$J_x=\underset{V}{\iiint }(y^2+z^2)\text{d}m$$
其它项类似。可以看出 $J$ 是刚体的属性，在与刚体固连的本体系下不随姿态改变。
  这里需要注意区分世界系和与刚体固连的本体系，也就是为什么需要进行一次坐标变换。在世界系中，刚体上各点的坐标不是定值，所以刚体的转动惯量矩阵在世界系下不是定值，只有在本体系下才是定值，所以这里需要转换坐标系。
  下面推导坐标系变换的最后一步，同时加上坐标系角标来清楚地标明每个变量属于哪个坐标系下的描述。
$$\begin{array}{rl}{}^w\vec{M}=& {}^w\dot{\vec{h}}\\ {}^w\vec{M}=& R{}^b\dot{\vec{h}}+{}^w\vec{\omega }\times {}^w\vec{h}\\ R^{-1}{}^w\vec{M}=& {}^b\dot{\vec{h}}+R^{-1}{}^w\vec{\omega }\times R^{-1}{}^w\vec{h}\\ {}^b\vec{M}=& {}^b\dot{\vec{h}}+{}^b\vec{\omega }\times {}^b\vec{h}\\ {}^b\vec{M}=& J{}^b\dot{\vec{w}}+{}^b\vec{\omega }\times J{}^b\vec{\omega }\end{array}$$
根据上面结果可以看出，姿态动力学方程 $\vec{M}=J\dot{\vec{w}}+\vec{w}\times J\vec{w}$ 中的所有量，包括角速度向量，都在本体系下描述。

推导2

下面的推导来自文末参考资料[2]，主要思路是所有变量都在世界系中处理。但我觉得这种方法不对。
  推导1中得到转动惯量矩阵 $J$ 后，$J$ 在本体系下保持不变，但在世界系下是变量，可以记作
$$J(t)=\sum _i\sum _j{J}_{ij}{\vec{a}}_i(t){\vec{a}}_j(t{)}^{\text{T}}$$
其中 $J_{ij}$ 是本体系下保持不变的转动惯量矩阵元素，${\vec{a}}_i,{\vec{a}}_j$ 表示本体系的坐标轴，例如当两个系重合时
$${\vec{a}}_2{\vec{a}}_3^{\text{T}}=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]\left[\begin{array}{ccc}0& 0& 1\end{array}\right]=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right]$$
但 ${\vec{a}}_i,{\vec{a}}_j$ 在旋转时是个变量，所以
$$\begin{array}{rl}M=& \frac{\text{d}H}{\text{d}t}=\frac{\text{d}}{\text{d}t}(J(t)\vec{w}(t))\\ =& \dot{J}(t)\vec{w}(t)+J(t)\dot{\vec{w}}(t)\\ =& \sum _i\sum _j{J}_{ij}({\dot{\vec{a}}}_i{\vec{a}}_j^{\text{T}}+{\vec{a}}_i{\dot{\vec{a}}}_j^{\text{T}})\vec{w}+J\dot{w}\\ =& \sum _i\sum _j{J}_{ij}[(\vec{w}\times {\vec{a}}_i){\vec{a}}_j^{\text{T}}\vec{w}+{\vec{a}}_i(\vec{w}\times {\vec{a}}_j{)}^{\text{T}}\vec{w}]+J\dot{w}\\ =& \sum _i\sum _j(\vec{w}\times J_{ij}{\vec{a}}_i){\vec{a}}_j^{\text{T}}\vec{w}+J\dot{w}\\ =& \vec{w}\times \left(\sum _i\sum _j{J}_{ij}{\vec{a}}_i{\vec{a}}_j^{\text{T}}\right)\vec{w}+J\dot{w}\\ =& \vec{w}(t)\times J(t)\vec{w}(t)+J(t)\dot{w}(t)\end{array}$$
其中 $(\vec{w}\times {\vec{a}}_j{)}^{\text{T}}\vec{w}=0$。

举例

  这个例子描述了没有外部力矩输入时角速度却一直在变化的现象，这一现象在参考资料[3]的第13章中有详细描述。
  设转动惯量矩阵 $J$、角速度初始值 ${\vec{w}}_0$ 分别为
$$J=\left[\begin{array}{ccc}2& 0& 0\\ 0& 2& 0\\ 0& 0& 1\end{array}\right],{\vec{w}}_0=\left[\begin{array}{c}1\\ 0\\ 1\end{array}\right]$$
当力矩为0时，由姿态动力学方程得
$$\begin{array}{rl}\dot{\vec{w}}=& -J^{-1}(\vec{w}\times J\vec{w})\\ =& -\left[\begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 0\\ 0& 0& 1\end{array}\right](\left[\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right]\times \left[\begin{array}{c}2w_x\\ 2w_y\\ w_z\end{array}\right])\\ =& -\left[\begin{array}{ccc}0.5& 0& 0\\ 0& 0.5& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}-w_y{w}_z\\ w_x{w}_z\\ 0\end{array}\right]\\ \left[\begin{array}{c}{\dot{w}}_x\\ {\dot{w}}_y\\ {\dot{w}}_z\end{array}\right]=& \left[\begin{array}{c}0.5w_y\\ -0.5w_x\\ 0\end{array}\right]\end{array}$$
解得
$$\begin{array}{rl}{w}_x(t)=& \cos 0.5t\\ w_y(t)=& -\sin 0.5t\\ w_z(t)=& 1\end{array}$$

附录：向量叉乘的乘法分配律

  补充一个在推导过程中遇见的一个公式。当 $A$ 是正交矩阵时，
$$A(\vec{b}\times \vec{c})=(A\vec{b})\times (A\vec{c})$$
$A$ 不是正交矩阵时暂时没发现什么有意思的性质。前面还有一个公式
$$A(A^{-1}\vec{b}\times \vec{b})$$
也没找到能合并$A$ 和 $A^{-1}$ 的方法。
证明：
$$\begin{array}{rl}(A\vec{b})\times (A\vec{c})=& [(b_x{\vec{a}}_1)+(b_y{\vec{a}}_2)+(b_z{\vec{a}}_3)]\times [(c_x{\vec{a}}_1)+(c_y{\vec{a}}_2)+(c_z{\vec{a}}_3)]\\ =& (b_x{\vec{a}}_1)\times (c_y{\vec{a}}_2)+(b_x{\vec{a}}_1)\times (c_z{\vec{a}}_3)+(b_y{\vec{a}}_2)\times (c_x{\vec{a}}_1)\\ & +(b_y{\vec{a}}_2)\times (c_z{\vec{a}}_3)+(b_z{\vec{a}}_3)\times (c_x{\vec{a}}_1)+(b_z{\vec{a}}_3)\times (c_y{\vec{a}}_2)\\ =& (b_x{c}_y-b_y{c}_x){\vec{a}}_3+(b_z{c}_x-b_x{c}_z){\vec{a}}_2+(b_y{c}_z-b_z{c}_y){\vec{a}}_1\\ =& A(\vec{b}\times \vec{c})\end{array}$$

参考

1. 章仁为.卫星轨道姿态动力学与控制[M].北京航空航天大学出版社.1998.
2. 刚体动力学中的欧拉方程是如何推导出来的？-知乎
3. de Ruiter. Spacecraft dynamics and control: an introduction. 2013.

进动现象仿真

  仿真一下前文提到的角速度向量在本体系下绕圈的现象，这一现象使用了本体极锥和空间极锥的概念（不知道是不是这么翻译），在参考资料[3]中有详细介绍。
  仿真结果见视频 https://www.bilibili.com/video/BV1Ck4y1F7Ba/。视频中白色线条代表刚体，黄色和橙黄色分别代表本体系的z轴和xy轴，红色线条表示角速度向量，绿色线条表示角动量向量。视频前半段摄像头在世界系下，后半段跟随刚体转动。
  基于 simucpp 和 raylib 的部分代码如下，完整代码见 https://gitee.com/xd15zhn/atttitudemotion。

/**************************************************************
// 验证刚体自由转动时的进动现象
simucpp版本：2.1.14
**************************************************************/
int main(void) {
    
    Mat J(3, 3, vecdble{
    2, 0, 0, 0, 2, 0, 0, 0, 1});
    Mat Jinv = J.inv();
    /* 初始化仿真器 */
    Simulator sim1;
    auto intQ = new MStateSpace(&sim1, BusSize(4, 1), true, "intQ");
    auto intW = new MStateSpace(&sim1, BusSize(3, 1), true, "intW");
    auto misoQ = new MFcnMISO(&sim1, BusSize(4, 1), "misoQ");
    auto misoW = new MFcnMISO(&sim1, BusSize(3, 1), "misoW");
    auto muxtau = new Mux(&sim1, BusSize(3, 1), "muxtau");
    UInput **intau = new UInput*[3];
    for (uint32_t i = 0; i < 3; i++) {
    
        intau[i] = new UInput(&sim1, "inw"+to_string(i));
        sim1.connectU(intau[i], muxtau, BusSize(i, 0));
    }
    sim1.connectM(intQ, misoQ);
    sim1.connectM(intW, misoQ);
    sim1.connectM(misoQ, intQ);
    sim1.connectM(intW, misoW);
    sim1.connectM(muxtau, misoW);
    sim1.connectM(misoW, intW);
    sim1.Set_SimStep(0.01);
    intQ->Set_InitialValue(Mat(vecdble{
    1, 0, 0, 0}));
    intW->Set_InitialValue(Mat(vecdble{
    1, 0, 1}));
    intau[0]->Set_Function([](double t){
     return 0; });
    intau[1]->Set_Function([](double t){
     return 0; });
    intau[2]->Set_Function([](double t){
     return 0; });
    misoQ->Set_Function([](Mat *u){
    
        // u[0]为姿态四元数(4x1)，u[1]为本体系下的角速度向量(3x1)
        Mat W(4, 4, vecdble{
      // 角速度向量w1对应的四元数乘法左矩阵(4x4)
            0, -u[1].at(0, 0), -u[1].at(1, 0), -u[1].at(2, 0),
            u[1].at(0, 0), 0, u[1].at(2, 0), -u[1].at(1, 0),
            u[1].at(1, 0), -u[1].at(2, 0), 0, u[1].at(0, 0),
            u[1].at(2, 0), u[1].at(1, 0), -u[1].at(0, 0), 0,
        });
        return 0.5*W*u[0];
    });
    misoW->Set_Function([&J, &Jinv](Mat *u){
    
        Vector3d omega = u[0];  // 本体系下的角速度向量(3x1)
        Vector3d tau = Mat(3, 1, vecdble{
      // 本体系下的力矩向量(3x1)
            u[1].at(0, 0), u[1].at(1, 0), u[1].at(2, 0)
        });
        Vector3d ans = Jinv*(tau - (omega & (J*omega)));
        return Mat(ans);
    });
    sim1.Initialize();
    sim1.Simulate_FirstStep();

    Mat curq, curw, curR, omegaI, hmom;
    /*初始化场景*/
    Camera camera;
    SetWindowMonitor(1);
	SetConfigFlags(FLAG_MSAA_4X_HINT);
	SetConfigFlags(FLAG_FULLSCREEN_MODE);
    SetTargetFPS(50);
    InitGraph(0, 0, "Universe");
	Init_Camera(&camera);

    /*场景循环*/
    while (!WindowShouldClose()) {
    
        curq = intQ->Get_OutValue();  // 姿态四元数
        curR = Quaternion_to_Rotation(curq);  // 旋转矩阵
        curw = intW->Get_OutValue();  // 本体系下的角速度向量
        omegaI = curR * curw;  // 世界系下的角速度向量
        hmom = curR * (J*curw);  // 世界系下的角动量
        if (Shift_View())
            SetReferenceFrame(Mat33_to_Matrix(curR));
        else
            SetReferenceFrame(MatrixIdentity());

		Update_Camera(&camera);
        BeginDrawing();
        ClearBackground(BLACK);
        BeginMode3D(camera);
            DrawGrid(120, 2);
            Draw_Body(curR);  // 画刚体和固连坐标系
            DrawLine3D(Vector3Zero(), Mat_to_Vector3(omegaI), RED);  // 画世界系下的角速度向量
            DrawLine3D(Vector3Zero(), Mat_to_Vector3(hmom), GREEN);  // 画世界系下的角动量向量
        EndMode3D();
        EndDrawing();

        if (IsKeyDown(KEY_Q)) {
    
            sim1.Simulation_Reset();
        }
        if (Pause()) continue;
        for (int i = 0; i < 5; i++)
            sim1.Simulate_OneStep();
    }
    CloseGraph();
    return 0;
}