描述姿态的三种方法:从理论到c语言代码

122 阅读10分钟

前言

在描述物体的姿态时,常用的方法有三种:欧拉角、四元数、旋转矩阵。它们有这各自的特点且可以互相转换。

描述姿态

欧拉角

定义:用三个绕特定坐标轴的连续旋转角度来描述姿态。最常见的顺序是**偏航(Yaw)- 俯仰(Pitch)- 滚动(Roll)。

  • 偏航(ψ):绕Z轴旋转,指向前后方向的变化。
  • 俯仰(θ):绕Y轴旋转,指向上下方向的变化。
  • 滚动(φ):绕X轴旋转,绕自身前后轴线的旋转。 优点:
  • 非常直观:飞行员和工程师可以直接理解“机头抬升20度,右倾10度”这样的描述。 缺点:
  • 万向节死锁:这是欧拉角最致命的缺陷。当俯仰角为±90度时,偏航和滚转的旋转轴会重合,丢失一个自由度,导致系统出现奇异性。
  • 计算复杂:进行多次旋转叠加时,计算效率低。

旋转矩阵

旋转矩阵是一个3x3的正交矩阵

定义:它的每一列代表新坐标系的X、Y、Z轴单位向量在原坐标系下的坐标。可以直接用来旋转三维向量或进行坐标系变换。

优点:

  • 无奇异性:没有万向节死锁问题。
  • 旋转向量直接:与向量乘法即可完成旋转,非常方便。
  • 组合简单:多个旋转可以直接通过矩阵乘法连接。 缺点:
  • 冗余:9个数字只表示3个自由度,存在6个约束(正交且行列式为1)。
  • 效率:存储和计算相比四元数更耗费资源。
  • 数值误差:多次运算后可能不再严格正交,需要重新正交化。

四元数

四元数是一种扩展的复数,由一个实部和三个虚部组成,通常表示为 q=[w,x,y,z]q = [w, x, y, z]q=w+xi+yj+zkq = w + xi + yj + zk

定义: 用于表示旋转的单位四元数满足 w² + x² + y² + z² = 1。它可以被解释为绕一个轴 [x, y, z] 旋转一个角度 θ,其中:

  • w=cos(θ/2)w = cos(θ/2)
  • (x,y,z)=sin(θ/2)(axisx,axisy,axisz)(x, y, z) = sin(θ/2) * (axis_x, axis_y, axis_z) 优点:    - 紧凑:只有4个元素,比旋转矩阵更节省存储空间。    - 无奇异性:没有万向节死锁问题。    - 计算高效:组合旋转和旋转向量所需的运算量比矩阵少。
    • 数值稳定:比旋转矩阵更容易保持单位性质。 缺点:
    • 不直观:人类很难直接理解一个四元数的几何意义。

转换方法

欧拉角与旋转矩阵

欧拉角 → 旋转矩阵
基本信息:
以z y x为顺序
偏航ψ, 俯仰θ, 横滚φ

以z y x为顺序
偏航ψ, 俯仰θ, 横滚φ
R=[
  [cosψ*cosθ,  cosψ*sinθ*sinφ - sinψ*cosφ,  cosψ*sinθ*cosφ + sinψ*sinφ],
  [sinψ*cosθ,  sinψ*sinθ*sinφ + cosψ*cosφ,  sinψ*sinθ*cosφ - cosψ*sinφ],
  [-sinθ,      cosθ*sinφ,                   cosθ*cosφ                 ]
]

旋转矩阵 → 欧拉角
基本信息:
以z y x为顺序
偏航ψ, 俯仰θ, 横滚φ

z y x为顺序
偏航ψ, 俯仰θ, 横滚φ
θ = -arcsin(R[2][0])         // 俯仰角
ψ = atan2(R[1][0], R[0][0])  // 偏航角 
φ = atan2(R[2][1], R[2][2])  // 横滚角 

四元数与旋转矩阵

四元数 → 旋转矩阵

归一化四元数q = [w, x, y, z]
旋转矩阵R = [
  [1 - 2y² - 2z²,     2xy - 2wz,       2xz + 2wy],
  [2xy + 2wz,         1 - 2x² - 2z²,   2yz - 2wx],
  [2xz - 2wy,         2yz + 2wx,       1 - 2x² - 2y²]
]

旋转矩阵 → 四元数

正交矩阵旋转矩阵 R = [
[m00, m01, m02],
[m10, m11, m12],
[m20, m21, m22]
]

1. 先算四个分量的平方值:
实部平方 = (1 + m00 + m11 + m22) / 4
虚部x平方 = (1 + m00 - m11 - m22) / 4
虚部y平方 = (1 - m00 + m11 - m22) / 4
虚部z平方 = (1 - m00 - m11 + m22) / 4

2. 找最大的平方值,分四种情况:
情况1:实部平方最大
实部 w = 开方(实部平方)
除数 = 1 / (4 × w)
虚部 x = (m21 - m12) × 除数
虚部 y = (m02 - m20) × 除数
虚部 z = (m10 - m01) × 除数

情况2:虚部x平方最大
虚部 x = 开方(虚部x平方)
除数 = 1 / (4 × x)
实部 w = (m21 - m12) × 除数
虚部 y = (m01 + m10) × 除数
虚部 z = (m02 + m20) × 除数

情况3:虚部y平方最大
虚部 y = 开方(虚部y平方)
除数 = 1 / (4 × y)
实部 w = (m02 - m20) × 除数
虚部 x = (m01 + m10) × 除数
虚部 z = (m12 + m21) × 除数

情况4:虚部z平方最大
虚部 z = 开方(虚部z平方)
除数 = 1 / (4 × z)
实部 w = (m10 - m01) × 除数
虚部 x = (m02 + m20) × 除数
虚部 y = (m12 + m21) × 除数

欧拉角与四元数

欧拉角 → 四元数
以z y x为顺序 
偏航ψ, 俯仰θ, 横滚φ

以z y x为顺序 
偏航ψ, 俯仰θ, 横滚φ
四元数_x = [cos(φ/2), sin(φ/2), 0, 0]
四元数_y = [cos(θ/2), 0, sin(θ/2), 0]
四元数_z = [cos(ψ/2), 0, 0, sin(ψ/2)]
四元数 = 四元数_z × 四元数_y × 四元数_x
/*********************************************/
实部 = cos(ψ/2)cos(θ/2)cos(φ/2) + sin(ψ/2)sin(θ/2)sin(φ/2)
虚部x = cos(ψ/2)cos(θ/2)sin(φ/2) - sin(ψ/2)sin(θ/2)cos(φ/2)
虚部y = cos(ψ/2)sin(θ/2)cos(φ/2) + sin(ψ/2)cos(θ/2)sin(φ/2)
虚部z = sin(ψ/2)cos(θ/2)cos(φ/2) - cos(ψ/2)sin(θ/2)sin(φ/2)

四元数 → 欧拉角
以z y x为顺序 
偏航ψ, 俯仰θ, 横滚φ

以z y x为顺序 
偏航ψ, 俯仰θ, 横滚φ
先归一化四元数 q = [w, x, y, z]
欧拉角 [ψ, θ, φ](偏航、俯仰、横滚)
θ = asin(2(wy - xz))
φ = atan2(2(wx + yz), 1 - 2(x² + y²))
ψ = atan2(2(wz + xy), 1 - 2(y² + z²))
// 更准确的:
θ = asin(2(wy - xz))
ψ = atan2(2(wz + xy), 1 - 2(y² + z²))  // 偏航角
φ = atan2(2(wx + yz), 1 - 2(x² + y²))  // 横滚角
/************************************************/
四元数 = [实部, 虚部x, 虚部y, 虚部z]
欧拉角 [偏航角, 俯仰角, 横滚角]
俯仰角 = 反正弦(2(实部×虚部y - 虚部x×虚部z))
横滚角 = 反正切2(2(实部×虚部x + 虚部y×虚部z), 1 - 2(虚部x² + 虚部y²))
偏航角 = 反正切2(2(实部×虚部z + 虚部x×虚部y), 1 - 2(虚部y² + 虚部z²))

示例代码

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846

typedef struct {
    double yaw;    // 偏航角 ψ (Z轴)
    double pitch;  // 俯仰角 θ (Y轴) 
    double roll;   // 横滚角 φ (X轴)
} EulerAngles;

typedef struct {
    double m[3][3];
} Matrix3x3;

typedef struct {
    double w, x, y, z;
} Quaternion;

// 安全数学函数
double safe_asin(double x) {
    if (x >= 1.0) return PI/2;
    if (x <= -1.0) return -PI/2;
    return asin(x);
}

// 弧度与度数转换函数
double rad_to_deg(double radians) {
    return radians * 180.0 / PI;
}

double deg_to_rad(double degrees) {
    return degrees * PI / 180.0;
}

// 欧拉角度数转换为弧度
EulerAngles euler_deg_to_rad(EulerAngles euler_deg) {
    EulerAngles euler_rad;
    euler_rad.yaw = deg_to_rad(euler_deg.yaw);
    euler_rad.pitch = deg_to_rad(euler_deg.pitch);
    euler_rad.roll = deg_to_rad(euler_deg.roll);
    return euler_rad;
}

// 欧拉角弧度转换为度数
EulerAngles euler_rad_to_deg(EulerAngles euler_rad) {
    EulerAngles euler_deg;
    euler_deg.yaw = rad_to_deg(euler_rad.yaw);
    euler_deg.pitch = rad_to_deg(euler_rad.pitch);
    euler_deg.roll = rad_to_deg(euler_rad.roll);
    return euler_deg;
}

// 归一化四元数
Quaternion normalize_quaternion(Quaternion q) {
    double norm = sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
    if (norm > 1e-6) {
        q.w /= norm;
        q.x /= norm;
        q.y /= norm;
        q.z /= norm;
    }
    return q;
}

// ZYX顺序:偏航ψ→俯仰θ→横滚φ
// 欧拉角 → 旋转矩阵
Matrix3x3 euler_to_matrix(EulerAngles euler) {
    double cos_psi = cos(euler.yaw);
    double sin_psi = sin(euler.yaw);
    double cos_theta = cos(euler.pitch);
    double sin_theta = sin(euler.pitch);
    double cos_phi = cos(euler.roll);
    double sin_phi = sin(euler.roll);
    
    Matrix3x3 R;
    
    // ZYX顺序的旋转矩阵:R = R_z(ψ) * R_y(θ) * R_x(φ)
    R.m[0][0] = cos_psi * cos_theta;
    R.m[0][1] = cos_psi * sin_theta * sin_phi - sin_psi * cos_phi;
    R.m[0][2] = cos_psi * sin_theta * cos_phi + sin_psi * sin_phi;
    
    R.m[1][0] = sin_psi * cos_theta;
    R.m[1][1] = sin_psi * sin_theta * sin_phi + cos_psi * cos_phi;
    R.m[1][2] = sin_psi * sin_theta * cos_phi - cos_psi * sin_phi;
    
    R.m[2][0] = -sin_theta;
    R.m[2][1] = cos_theta * sin_phi;
    R.m[2][2] = cos_theta * cos_phi;
    
    return R;
}

// 旋转矩阵 → 欧拉角 (ZYX顺序)
EulerAngles matrix_to_euler(Matrix3x3 R) {
    EulerAngles euler;
    
    // 俯仰角 θ = -arcsin(R[2][0])
    double sin_pitch = -R.m[2][0];
    euler.pitch = safe_asin(sin_pitch);
    
    double cos_theta = cos(euler.pitch);
    if (fabs(cos_theta) > 1e-6) {
        // 修正:不需要除以cosθ
        euler.yaw = atan2(R.m[1][0], R.m[0][0]);    // 偏航角
        euler.roll = atan2(R.m[2][1], R.m[2][2]);   // 横滚角
    } else {
        // 万向锁处理保持不变
        if (euler.pitch > 0) {
            euler.roll = 0.0;
            euler.yaw = atan2(R.m[0][1], R.m[1][1]);
        } else {
            euler.roll = 0.0;
            euler.yaw = atan2(-R.m[0][1], -R.m[1][1]);
        }
    }
    
    return euler;
}

// 四元数 → 旋转矩阵
Matrix3x3 quaternion_to_matrix(Quaternion q) {
    q = normalize_quaternion(q);
    
    Matrix3x3 R;
    double xx = q.x * q.x, yy = q.y * q.y, zz = q.z * q.z;
    double xy = q.x * q.y, xz = q.x * q.z, yz = q.y * q.z;
    double wx = q.w * q.x, wy = q.w * q.y, wz = q.w * q.z;
    
    R.m[0][0] = 1 - 2 * (yy + zz);
    R.m[0][1] = 2 * (xy - wz);
    R.m[0][2] = 2 * (xz + wy);
    
    R.m[1][0] = 2 * (xy + wz);
    R.m[1][1] = 1 - 2 * (xx + zz);
    R.m[1][2] = 2 * (yz - wx);
    
    R.m[2][0] = 2 * (xz - wy);
    R.m[2][1] = 2 * (yz + wx);
    R.m[2][2] = 1 - 2 * (xx + yy);
    
    return R;
}

// 旋转矩阵 → 四元数
Quaternion matrix_to_quaternion(Matrix3x3 R) {
    double m00 = R.m[0][0], m01 = R.m[0][1], m02 = R.m[0][2];
    double m10 = R.m[1][0], m11 = R.m[1][1], m12 = R.m[1][2];
    double m20 = R.m[2][0], m21 = R.m[2][1], m22 = R.m[2][2];
    
    Quaternion q;
    double tr = m00 + m11 + m22;
    
    if (tr > 0) {
        double S = sqrt(tr + 1.0) * 2;
        q.w = 0.25 * S;
        q.x = (m21 - m12) / S;
        q.y = (m02 - m20) / S;
        q.z = (m10 - m01) / S;
    } else if ((m00 > m11) && (m00 > m22)) {
        double S = sqrt(1.0 + m00 - m11 - m22) * 2;
        q.w = (m21 - m12) / S;
        q.x = 0.25 * S;
        q.y = (m01 + m10) / S;
        q.z = (m02 + m20) / S;
    } else if (m11 > m22) {
        double S = sqrt(1.0 + m11 - m00 - m22) * 2;
        q.w = (m02 - m20) / S;
        q.x = (m01 + m10) / S;
        q.y = 0.25 * S;
        q.z = (m12 + m21) / S;
    } else {
        double S = sqrt(1.0 + m22 - m00 - m11) * 2;
        q.w = (m10 - m01) / S;
        q.x = (m02 + m20) / S;
        q.y = (m12 + m21) / S;
        q.z = 0.25 * S;
    }
    
    return normalize_quaternion(q);
}

// 欧拉角 → 四元数 (ZYX顺序)
Quaternion euler_to_quaternion(EulerAngles euler) {
    double half_psi = euler.yaw * 0.5;
    double half_theta = euler.pitch * 0.5;
    double half_phi = euler.roll * 0.5;
    
    double cos_psi = cos(half_psi);
    double sin_psi = sin(half_psi);
    double cos_theta = cos(half_theta);
    double sin_theta = sin(half_theta);
    double cos_phi = cos(half_phi);
    double sin_phi = sin(half_phi);
    
    Quaternion q;
    
    // ZYX顺序:q = q_z(ψ) * q_y(θ) * q_x(φ)
    q.w = cos_psi * cos_theta * cos_phi + sin_psi * sin_theta * sin_phi;
    q.x = cos_psi * cos_theta * sin_phi - sin_psi * sin_theta * cos_phi;
    q.y = cos_psi * sin_theta * cos_phi + sin_psi * cos_theta * sin_phi;
    q.z = sin_psi * cos_theta * cos_phi - cos_psi * sin_theta * sin_phi;
    
    return normalize_quaternion(q);
}

// 四元数 → 欧拉角 (ZYX顺序) - 修正版本
EulerAngles quaternion_to_euler(Quaternion q) {
    q = normalize_quaternion(q);
    
    EulerAngles euler;
    
    // 俯仰角 θ = asin(2*(w*y - x*z))
    double sin_pitch = 2.0 * (q.w * q.y - q.x * q.z);
    euler.pitch = safe_asin(sin_pitch);
    
    // 检查万向锁 (俯仰角接近 ±90°)
    double cos_theta = cos(euler.pitch);
    if (fabs(cos_theta) > 1e-6) {
        // 非万向锁情况
        // 偏航角 ψ = atan2(2*(w*z + x*y), 1 - 2*(y*y + z*z))
        euler.yaw = atan2(2.0 * (q.w * q.z + q.x * q.y), 
                          1.0 - 2.0 * (q.y*q.y + q.z*q.z));
        
        // 横滚角 φ = atan2(2*(w*x + y*z), 1 - 2*(x*x + y*y))
        euler.roll = atan2(2.0 * (q.w * q.x + q.y * q.z),
                           1.0 - 2.0 * (q.x*q.x + q.y*q.y));
    } else {
        // 万向锁情况:俯仰角接近 ±90°
        // 此时只能确定 ψ - φ 或 ψ + φ
        if (euler.pitch > 0) { // θ = +90°
            euler.roll = 0.0;
            euler.yaw = 2.0 * atan2(q.x, q.w);
        } else { // θ = -90°
            euler.roll = 0.0;
            euler.yaw = -2.0 * atan2(q.x, q.w);
        }
    }
    
    return euler;
}

// 打印函数
void print_euler_rad(EulerAngles euler) {
    printf("偏航(ψ): %.6f rad, 俯仰(θ): %.6f rad, 横滚(φ): %.6f rad\n", 
           euler.yaw, euler.pitch, euler.roll);
}

void print_euler_deg(EulerAngles euler) {
    EulerAngles euler_deg = euler_rad_to_deg(euler);
    printf("偏航(ψ): %.6f°, 俯仰(θ): %.6f°, 横滚(φ): %.6f°\n", 
           euler_deg.yaw, euler_deg.pitch, euler_deg.roll);
}

void print_matrix(Matrix3x3 R) {
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            printf("%.6f ", R.m[i][j]);
        }
        printf("\n");
    }
}

void print_quaternion(Quaternion q) {
    printf("w: %.6f, x: %.6f, y: %.6f, z: %.6f\n", q.w, q.x, q.y, q.z);
}

// 测试函数
int main() {
    printf("=== ZYX顺序转换测试 ===\n");
    
    // 测试1:普通情况
    printf("\n--- 测试1:普通情况 ---\n");
    EulerAngles euler_rad = {0.5, 0.3, 0.2};  // ψ, θ, φ
    
    printf("原始欧拉角: ");
    print_euler_rad(euler_rad);
    print_euler_deg(euler_rad);
    
    // 欧拉角 → 旋转矩阵
    Matrix3x3 R = euler_to_matrix(euler_rad);
    printf("\n旋转矩阵:\n");
    print_matrix(R);
    
    // 旋转矩阵 → 欧拉角
    EulerAngles euler_back = matrix_to_euler(R);
    printf("\n矩阵转换回的欧拉角: ");
    print_euler_rad(euler_back);
    
    // 欧拉角 → 四元数
    Quaternion q = euler_to_quaternion(euler_rad);
    printf("\n四元数: ");
    print_quaternion(q);
    
    // 四元数 → 欧拉角
    EulerAngles euler_from_q = quaternion_to_euler(q);
    printf("四元数转换回的欧拉角: ");
    print_euler_rad(euler_from_q);
    
    // 一致性检查
    printf("\n一致性误差:\n");
    printf("矩阵转换: ψ=%.6f, θ=%.6f, φ=%.6f\n", 
           fabs(euler_rad.yaw - euler_back.yaw),
           fabs(euler_rad.pitch - euler_back.pitch),
           fabs(euler_rad.roll - euler_back.roll));
    printf("四元数转换: ψ=%.6f, θ=%.6f, φ=%.6f\n",
           fabs(euler_rad.yaw - euler_from_q.yaw),
           fabs(euler_rad.pitch - euler_from_q.pitch),
           fabs(euler_rad.roll - euler_from_q.roll));
    
    // 测试2:万向锁情况
    printf("\n--- 测试2:万向锁情况 (俯仰角=90°) ---\n");
    EulerAngles gimbal_lock = {0.8, PI/2, 0.4};  // 精确的90度
    printf("原始欧拉角: ");
    print_euler_deg(gimbal_lock);
    
    Matrix3x3 R_lock = euler_to_matrix(gimbal_lock);
    EulerAngles euler_lock_back = matrix_to_euler(R_lock);
    printf("矩阵转换回的欧拉角: ");
    print_euler_deg(euler_lock_back);
    
    Quaternion q_lock = euler_to_quaternion(gimbal_lock);
    EulerAngles euler_lock_q = quaternion_to_euler(q_lock);
    printf("四元数转换回的欧拉角: ");
    print_euler_deg(euler_lock_q);
    
    // 测试3:接近万向锁情况
    printf("\n--- 测试3:接近万向锁情况 (俯仰角≈90°) ---\n");
    EulerAngles near_gimbal_lock = {0.8, PI/2 - 0.001, 0.4};
    printf("原始欧拉角: ");
    print_euler_deg(near_gimbal_lock);
    
    Matrix3x3 R_near_lock = euler_to_matrix(near_gimbal_lock);
    EulerAngles euler_near_back = matrix_to_euler(R_near_lock);
    printf("矩阵转换回的欧拉角: ");
    print_euler_deg(euler_near_back);
    
    Quaternion q_near_lock = euler_to_quaternion(near_gimbal_lock);
    EulerAngles euler_near_q = quaternion_to_euler(q_near_lock);
    printf("四元数转换回的欧拉角: ");
    print_euler_deg(euler_near_q);
    
    return 0;
}