无人机姿态解算代码讲解

无人机姿态解算代码讲解

姿态解算,对于无人机飞控算法来说极其重要,无人机的自稳飞行以及后面的一系列的高度数据融合,GPS数据融合、光流数据融合等都会有影响,所以理解好姿态解算这个算法还是很有必要的。
关于太多的理论的知识已经有很多的网上的解说,这次我写的是我这几年对姿态结算的代码的理解,怎么去调参数,哪些代码起着什么样的作用。
先上个姿态解算的结构框图: 在这里插入图片描述
姿态解算输入数据为三轴加速计,三轴陀螺仪。输出数据为欧拉角。

void IMUupdate(float gxi,float gyi,float gzi,float axi,float ayi,float azi,__Euler *euler)//姿态解算输入数据为三轴加速计,三轴陀螺仪。
{
float recipNormmag=0.0f;
float halfvx=0.0f, halfvy=0.0f, halfvz=0.0f;
float halfex=0.0f, halfey=0.0f, halfez=0.0f;
 float qa=0.0f, qb=0.0f, qc=0.0f;
 float q0q0 = q0*q0;
 float q0q1 = q0*q1;
 float q0q2 = q0*q2;
 float q0q3 = q0*q3;
 float q1q1 = q1*q1;
 float q1q2 = q1*q2;
 float q1q3 = q1*q3;
 float q2q2 = q2*q2;   
 float q2q3 = q2*q3;
 float q3q3 = q3*q3;
 static int last_time=0;int now_time=0;
 //前面都是算法的变量定义
 if(euler==NULL)
  return;
  now_time=SystemGetMirco();
  if(last_time!=0)
 {
  sampleFreq=1000000.0f/(float)(SystemGetMircoTim());
  //定义姿态解算时间周期,算法运行时间跟实际时间必须一致,否则算法会出错
 }
 last_time=now_time;
 ax=axi;ay=ayi;az=azi;gx=gxi;gy=gyi;gz=gzi;
 // Compute feedback only if accelerometer measurement valid (avoids NaN in accelerometer normalisation)
 if(!((ax == 0.0f) && (ay == 0.0f) && (az == 0.0f))) 
 //
 {
 // Normalise accelerometer measurement
 //加速计做归一化,确保计算的数据比例一样
  recipNormmag = Fylib_Qrsqrt(ax * ax + ay * ay + az * az);
  ax *= recipNormmag;
  ay *= recipNormmag;
  az *= recipNormmag;
// Estimated direction of gravity and vector perpendicular to magnetic flux
  halfvx = 2*(q1 * q3 - q0 * q2);
  halfvy = 2*(q0 * q1 + q2 * q3);
  halfvz = q0*q0 - q1*q1 - q2*q2 + q3*q3;
  //计算各轴XYZ在重力方向上的分量,后面互补滤波需要用到这个数据
  // Error is sum of cross product between estimated and measured direction of gravity
  halfex = (ay * halfvz - az * halfvy);// + Mag_ex * _w_mag;
  halfey = (az * halfvx - ax * halfvz);// + Mag_ey * _w_mag;
  halfez = (ax * halfvy - ay * halfvx);// + Mag_ez * _w_mag;//(Mag_ez - corrrect_values1) * _w_mag;
//各轴分量与加速计在重力下各轴分量之差,这个是互补滤波的核心,halfvx 是估计分量误差,acc x y z是加速计分量,不一样的角度加速计在重力作用下的分量是固定的,所以互补滤波以此作为绝对参数数据
	integralFBx += twoKi * halfex * (1.0f / sampleFreq); // integral error scaled by Ki
    integralFBy += twoKi * halfey * (1.0f / sampleFreq);
    integralFBz += twoKi * halfez * (1.0f / sampleFreq);
    //互补滤波KI分量,KI参数可调,会影响更相信加速计或陀螺仪
    gx += integralFBx; // apply integral feedback
   gy += integralFBy;
   gz += integralFBz;
   //KI分量加入欧拉角补偿
   // Apply proportional feedback
  
  gx += twoKp * halfex;
  gy += twoKp * halfey;
  gz += twoKp * halfez;
  互补滤波KP分量加入补偿
  qa = q0;
 qb = q1;
 qc = q2;
 q0 += (-qb * gx - qc * gy - q3 * gz);
 q1 += (qa * gx + qc * gz - q3 * gy);
 q2 += (qa * gy - qb * gz + q3 * gx);
 q3 += (qa * gz + qb * gy - qc * gx); 
 四元素更新
 // Normalise quaternion
 recipNormmag = Fylib_Qrsqrt(q0 * q0 + q1 * q1 + q2 * q2 + q3 * q3);
 q0 *= recipNormmag;
 q1 *= recipNormmag;
 q2 *= recipNormmag;
 q3 *= recipNormmag;
 //四元素归一化
 euler->pitch_x = atan2(2*q0*q1 + 2*q2*q3, 1 -2*q1*q1 - 2*q2*q2)* MV_RAD_DEG;//pitch
 euler->roll_y = asin(2*q0*q2 - 2*q1*q3)*MV_RAD_DEG;//roll
 euler->yaw_z = atan2(2*q0*q3+ 2*q1*q2 , 1 -2*q2*q2 - 2*q3*q3)*MV_RAD_DEG + 270;// - Gps_declination;//yaw
 //用四元素来代表欧拉角,这个公式很多论文都有推导方法,看懂其推倒也是需要一定的数学功底
 
}

以上为传统的四元素互补滤波解算姿态角,用的过程主要调试KP KI两个分量。另外加速计跟陀螺仪在时间上会有一定的反应时间差,所以加速计陀螺仪理论上是要加时间延迟缓冲的。熟悉使用这个姿态解算,不仅在无人机可用,其他很多需要平衡类的产品均可用。关注我,我会持续更新有关无人机各个模块的知识,带你完整的熟悉无人机各个模块。技术交流也可以加我扣扣1836703877.

### 基于 MATLAB 的无人机姿态实现 #### 扩展卡尔曼滤波器(EKF)用于无人机姿态 扩展卡尔曼滤波器是一种常用的非线性系统状态估计工具,在无人机姿态中有广泛应用。它通过对 IMU 数据和其他传感器数据的融合,能够有效提高姿态估计精度[^1]。 以下是基于 EKF 的无人机姿态的核心步骤及其 MATLAB 实现: 1. **定义状态变量** 定义无人机的状态向量 \( \mathbf{x} = [\phi, \theta, \psi]^T \),其中 \( \phi \), \( \theta \), 和 \( \psi \) 分别表示滚转角、俯仰角和偏航角。 2. **建立运动学模型** 使用四元数或者欧拉角描述无人机姿态变化过程,并构建离散时间的动力学方程。假设加速度计和陀螺仪分别提供线性和角速率测量,则可以得到如下形式的状态转移矩阵: ```matlab function x_next = state_transition(x_current, w, dt) % 输入参数说明: % x_current: 当前时刻的状态向量 [q0; q1; q2; q3] % w: 陀螺仪测得的角度增量 [wx; wy; wz] % dt: 时间步长 Omega = [0 -w(1)*dt -w(2)*dt -w(3)*dt; w(1)*dt 0 w(3)*dt -w(2)*dt; w(2)*dt -w(3)*dt 0 w(1)*dt; w(3)*dt w(2)*dt -w(1)*dt 0]; Q = expm(Omega/2); % 四元数更新公式 x_next = Q * x_current; % 更新后的四元数 end ``` 3. **设计观测模型** 利用加速度计的数据来校正预测值中的偏差。通常情况下,可以通过重力矢量的方向重建角度信息作为观测量。 ```matlab function z_hat = observation_model(q_estimated) g_vector_body_frame = quatrotate(conj(q_estimated'), [0; 0; 9.81]'); z_hat = atan2(-g_vector_body_frame(1), sqrt(g_vector_body_frame(2)^2 + g_vector_body_frame(3)^2)); end ``` 4. **执行 EKF 迭代** 结合上述两部分完成完整的 EKF 循环操作流程,具体代码片段展示如下所示: ```matlab P = eye(size(Q)); % 初始化协方差矩阵 R = diag([sigma_acc^2 sigma_gyro^2]); % 测量噪声协方差 for k=1:length(data)-1 u_k = data(k).gyro; % 获取当前时刻输入信号 y_k = data(k).accelero; % 获取当前时刻输出信号 % 预测阶段 x_pred = f(x_prev,u_k); F = jacobian_f(); % 计雅各比矩阵Jf P_pred = F*P*F'+Qk; % 更正阶段 H = jacobian_h(); K = P_pred*H'/(H*P_pred*H'+R); e_z = h(x_pred)-y_k; x_updated = x_pred-K*e_z; P = (eye(length(K))-K*H)*P_pred; save_results(k)=x_updated; end ``` #### 错误状态卡尔曼滤波器(ESKF) 另一种常用的方法是错误状态卡尔曼滤波器,其主要特点是围绕名义轨迹展开泰勒级数近似处理非线性函数关系[^2]。这种方法特别适合多源异构传感环境下的高动态场景应用场合。 --- ### 总结与建议 无论是采用传统的 EKF 抑或是更先进的 ESKF 方法,都需要针对实际硬件平台特性调整相应配置参数并验证性能表现。此外,为了进一步提升鲁棒性还可以引入额外辅助设备比如 GPS 或者视觉 SLAM 提供外部约束条件从而形成闭环控制系统架构[^3]。
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值