CN105865432A - 一种针对陀螺仪多源噪声的混合滤波方法与测试平台 - Google Patents

Abstract

一种针对陀螺仪多源噪声的混合滤波方法与测试平台，所述的测试平台包括实时仿真目标机、三轴速率转台、三轴陀螺仪、星敏感器、星仿真器、姿态确定模块、姿态控制模块以及执行机构；所述针对陀螺仪多源噪声的混合滤波方法，对陀螺仪的器件噪声采用H₂优化进行抑制，对滤波系统的建模误差采用H_∞优化进行抑制，结合H₂优化与H_∞优化设计混合滤波器使得估计误差满足性能指标；所设计的针对陀螺仪多源噪声的混合滤波方法有效提高姿态确定精度，该方法实时运行于测试平台的姿态确定模块；本发明测试平台验证针对陀螺仪多源噪声的混合滤波方法的有效性与工程实用性，适用于航空航天领域的地面仿真验证。

Description

一种针对陀螺仪多源噪声的混合滤波方法与测试平台

技术领域

本发明涉及一种针对陀螺仪多源噪声的混合滤波方法与测试平台，该测试平台可用于验证包括针对陀螺仪多源噪声的混合滤波方法在内的多种滤波方法；针对陀螺仪多源噪声提出了一种混合滤波方法对其抑制，该方法能显著提高陀螺仪的滤波精度，进一步改善姿控系统的控制性能。

背景技术

陀螺仪作为敏感角运动的一种装置，它不仅是运载体航行驾驶的重要仪器仪表，而且是运载体控制系统、惯性导航系统和惯性制导系统的核心元件，因而成为惯性技术的一个重要组成部分。陀螺仪含有多源噪声，可将多源噪声分为两大类，一是噪声统计特性可知的器件噪声，二是能量有界的滤波系统建模误差。实际应用中，在补偿陀螺仪部分噪声(如常值漂移)后，陀螺仪器件噪声还包括高斯白噪声、互相关随机漂移、角度随机游走、随机速率斜漂等，这类噪声由白噪声驱动并且统计特性可知，应用H₂优化可有效抑制；在飞行器尤其是航天器中，陀螺仪需要与其它敏感器，如太阳敏感器、星敏感器、地球敏感器及加速度计等组合使用才能完成姿态信息获取，在组合测量过程中存在滤波系统的建模误差，这也是影响滤波精度的一种干扰，可被视为能量有界的噪声，因此采用H_∞优化可有效抑制滤波系统的建模误差；在飞行器姿控系统中陀螺仪选型确定的情况下，设计滤波方法进一步提高测量精度，进而改善控制性能是一种重要思路。因此，研究针对陀螺仪多源噪声的混合滤波方法具有重要意义。

姿控系统测试平台已经在飞行器研制过程中得到了广泛应用，国内航天五院502所、航天八院812所，还有哈工大、清华、北航等高校都曾研制出姿控系统测试平台。但现有的这些姿控系统测试平台通常针对特定型号的飞行器而研制，未见专用于研究陀螺仪对姿控系统性能影响的测试平台。中国专利CN102530267B和中国专利CN103197669B都涉及航天器的测试平台，但两者都未针对陀螺仪多源噪声的滤波方法展开研究，前者提出一种卫星公用平台，后者则关注卫星的姿态控制模式。目前的姿控系统测试平台无法用于研究针对陀螺仪多源噪声的滤波方法，多种针对陀螺仪噪声的滤波方法也无法在现有的测试平台上完成分析、测试与评估。

发明内容

本发明的技术解决问题是：所设计的一种针对陀螺仪多源噪声的混合滤波方法与测试平台，克服现有姿控系统测试平台对陀螺仪多源噪声研究的缺失，该测试平台研究飞行器飞行过程中陀螺仪噪声对姿控系统控制精度的影响；针对陀螺仪多源噪声提出一种混合滤波方法加以抑制，该方法提高滤波精度，改善控制性能。

本发明的技术解决方案是：一种针对陀螺仪多源噪声的混合滤波方法与测试平台，其特征在于：测试平台包括实时仿真目标机(1)、三轴速率转台(2)、三轴陀螺仪(3)、星敏感器(4)、星仿真器(5)、姿态确定模块(6)、姿态控制模块(7)以及执行机构(8)；所述实时仿真目标机(1)包括飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)，分别实时运算飞行器的姿态动力学与运动学模型；所述三轴速率转台(2)是飞行器的运动仿真器，提供与飞行器相似或等价的运动，用于模拟飞行器的姿态运动；所述三轴陀螺仪(3)是角速度姿态信息敏感器，敏感三轴速率转台(2)的角速度；所述星敏感器(4)通过敏感恒星方位而获得卫星姿态信息，用来修正三轴陀螺仪(3)的漂移与噪声；所述星仿真器(5)为星敏感器(4)提供模拟恒星方位；所述姿态确定模块(6)实时运行针对陀螺仪多源噪声的混合滤波方法，对三轴陀螺仪(3)与星敏感器(4)实测的姿态信息滤波，包括了数字运算单元(61)、RS422/232通信单元(62)、电源单元(63)、复位单元(64)、扩展备用I/O口(65)、存储器(66)以及时钟单元(67)，其中RS422/232通信单元(62)、电源单元(63)、复位单元(64)、扩展备用I/O口(65)、存储器(66)、时钟单元(67)都与数字运算单元(61)相连，数字运算单元(61)是姿态确定模块(6)的运算核心，RS422/232通信单元(62)是姿态确定模块(6)的输入输出端，接收三轴陀螺仪(3)与星敏感器(4)实测姿态信息，输出滤波后的姿态信息，电源单元(63)负责供电，复位单元(64)可恢复数字运算单元(61)的初始状态，扩展备用I/O口(65)提供可扩展的输入输出端口，存储器(66)存储数字运算单元(61)的程序与数据，时钟单元(67)为数字运算单元(61)提供时钟频率；所述姿态控制模块(7)实时运算姿态控制算法，为执行机构(8)提供力矩控制指令；所述执行机构(8)包括了反作用飞轮(81)和磁力矩器(82)，在接收力矩控制指令后，将输出力矩信号传给实时仿真目标机(1)；测试平台的数据流回路如下：飞行器的期望姿态信号通过姿态控制模块(7)的RS232串行数据通信接口导入测试平台，该期望姿态信号通过与姿态确定模块(6)滤波后的姿态信息进行比较，得到偏差信号；将此偏差信号传递给姿态控制模块(7)，该模块中的姿态控制算法解算出指令力矩信号并对指令力矩信号进行分配，分配后的信号经RS422/232串行数据通信接口作用于执行机构(8)；执行机构(8)中反作用飞轮(81)和磁力矩器(82)接收到分配后的指令力矩信号，联合作用输出力矩信号，该力矩信号通过RS422/232串行数据通信接口传入实时仿真目标机(1)中的飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)；飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)分别运算飞行器姿态动力学与运动学模块，得到飞行器转动角速度，该飞行器转动角速度经RS422/232串行数据通信接口传送至三轴速率转台(2)；三轴速率转台(2)受飞行器转动角速度驱动而转动，装配在三轴速率转台(2)上的三轴陀螺仪(3)以及星敏感器(4)将敏感到的姿态信息经RS422/232串行数据通信接口传递给姿态确定模块(6)；姿态确定模块(6)中针对陀螺仪多源噪声的混合滤波方法解算得到滤波后的姿态信息，该姿态信息与期望姿态作比较后得到新的偏差信号，形成了测试平台的数据流回路。

针对陀螺仪多源噪声的混合滤波方法与测试平台，其特征在于：所述的混合滤波方法是一种递推算法，该方法的滤波周期包括以下步骤：首先构建含有陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程；其次设计一个线性滤波器；然后构造增广的状态估计方程与待估计信号估计误差方程；最后是利用MATLAB中的LMI工具箱求解最优化参数；

第一步：设计含陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程；

状态方程含陀螺仪器件噪声与滤波系统建模误差，其状态变量是用四元数形式表示的飞行器姿态角；星敏感器的量测方程为状态方程提供修正信息，状态方程、量测方程以及待估计信号方程为Σ₁：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A\left(t\right)x\left(t\right)+B\left(t\right)n\left(t\right)+w_1\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Dn\left(t\right)+w_1\left(t\right)\\ z\left(t\right)=Lx\left(t\right)\end{array}\right.$

其中，t表示时间，x(t)＝[q₀(t) q₁(t) q₂(t) q₃(t)]^T为四元数形式表示的飞行器姿态角；为x(t)的一阶微分；n(t)＝[n_x(t) n_y(t) n_z(t)]^T为陀螺仪滤波系统的建模误差；w₁(t)＝[δ₀(t) δ₁(t) δ₂(t) δ₃(t)]^T为四元数形式表示的陀螺仪器件噪声；y(t)＝[y₀(t) y₁(t) y₂(t) y₃(t)]^T为四元数形式表示的星敏感器测得的恒星方位；z(t)＝[z₀(t) z₁(t) z₂(t) z₃(t)]^T为系统待估计的信号；式中，矩阵C是量测矩阵，矩阵D是量测方程建模误差系数矩阵，矩阵L是为待估计信号系数矩阵，矩阵C、矩阵D、矩阵L都为四阶单位阵，系统状态矩阵A(t)为状态方程建模误差系数矩阵B(t)为w_x(t)、w_y(t)、w_z(t)分别为飞行器相对轨道坐标系的滚动轴角速率、俯仰轴角速率及偏航轴角速率；

将含有陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程离散化，可以得到如下模型Σ₂：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}x(k+1)=\stackrel{\‾}{A}\left(k\right)x\left(k\right)+\stackrel{\‾}{B}\left(k\right)n\left(k\right)+w_1\left(k\right)\\ y\left(k\right)=Cx\left(k\right)+Dn\left(k\right)+w_1\left(k\right)\\ z\left(k\right)=Lx\left(k\right)\end{array}\right.$

其中x(k+1)为k+1时刻以四元数形式表示的状态变量，x(k)为k时刻以四元数形式表示的状态变量，n(k)为k时刻的陀螺仪滤波系统的建模误差，w₁(k)为k时刻以四元数形式表示的陀螺仪器件噪声，y(k)为k时刻以四元数形式表示的星敏感器测得的恒星方位，z(k)为k时刻待估计的信号；A(k)为k时刻的系统状态矩阵，B(k)为k时刻的状态方程建模误差系数矩阵，I为单位阵，T为采样时间；

第二步：针对陀螺仪器件噪声和滤波系统建模误差设计线性滤波器；

$\left\{\begin{array}{c}\hat{x}(k+1)=\hat{A}\hat{x}\left(k\right)+\hat{B}y\left(k\right)\\ \hat{z}\left(k\right)=\hat{C}\hat{x}\left(k\right)\end{array}\right.$

式中为k+1时刻的状态变量估计值，为k时刻的状态变量估计值，为k时刻的待估计信号估计值，为线性滤波器的系统矩阵，为线性滤波器的输入矩阵，为线性滤波器的输出矩阵；

第三步：构造增广的状态估计方程与待估计信号估计误差方程；

目的是使得增广的状态估计方程与待估计信号估计误差方程

$\left\{\begin{array}{c}\tilde{x}(k+1)=\tilde{A}\tilde{x}\left(k\right)+\tilde{B}n\left(k\right)+{\tilde{B}}_1{w}_1\left(k\right)\\ \tilde{z}\left(k\right)=\tilde{C}\tilde{x}\left(k\right)\end{array}\right.$

式中第一个表达式为增广的状态估计方程，第二个方程是待估计信号的估计误差方程，这两个方程是渐近稳定的，且对应于通道的滤波误差方差的一个上界最小化，对应于通道的滤波误差向量满足其中,

$\tilde{x}\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ \hat{x}\left(k\right)\end{array}\right]\tilde{A}=\left[\begin{array}{cc}\stackrel{\‾}{A}\left(k\right)& 0\\ \hat{B}C& \hat{A}\end{array}\right],\tilde{B}=\left[\begin{array}{c}\stackrel{\‾}{B}\left(k\right)\\ \hat{B}D\end{array}\right],{\tilde{B}}_1=\left[\begin{array}{c}I\\ \hat{B}\end{array}\right]$

$\tilde{C}=\left[\begin{array}{cc}L& -\hat{C}\end{array}\right],\tilde{z}\left(k\right)=z\left(k\right)-\hat{z}\left(k\right)$

第四步：利用MATLAB中的LMI工具箱求解最优化参数；

对给定的常数γ>0，系统Σ₂存在一个混合滤波器可以归纳为以下的优化问题：

$\underset{R,P,Q,S,U,H}{\min }Trace\left(H\right)$

$s.t.\left(i\right)\left[\begin{array}{cccccc}R& *& *& *& *& *\\ R& P& *& *& *& *\\ 0& 0& {\mathrm{\γ}}^{2}I& *& *& *\\ R\stackrel{\‾}{A}\left(k\right)& R\stackrel{\‾}{A}\left(k\right)& R& R& *& *\\ P\stackrel{\‾}{A}\left(k\right)+QC+S& P\stackrel{\‾}{A}\left(k\right)+QC& P+Q& R& P& *\\ L-U& L& 0& 0& 0& I\end{array}\right]>0$

$\left(ii\right)\left[\begin{array}{ccc}H& *& *\\ R\stackrel{\‾}{B}\left(k\right)& R& *\\ P\stackrel{\‾}{B}\left(k\right)+QD& R& P\end{array}\right]\≥0$

其中*处的子块可由矩阵的对称性得到；进而求解出优化问题的最优解R、P、S、Q、U、H，其中R、P、S、Q、U、H为求解最优化问题定义的辅助变量，Trace表示矩阵的迹，即为矩阵对角元素之和，则：

$\hat{A}={(R-P)}^{-1}S,\hat{B}={(R-P)}^{-1}Q,\hat{C}=U$

是所要求的混合滤波器的系数矩阵，相应的

由于优化问题是一个具有线性矩阵不等式约束和线性目标函数的凸优化问题，因此可以应用MATLAB中LMI工具箱的求解器最优化最小值mincx来求解该问题；

第五步：判断终止条件，若满足，则滤波方法运行结束，若不满足，则循环调用；

终止条件为滤波方法的总运行时间，当该方法运行时间没有达到总运行时间，则更新时间，以下一个离散时刻k+1代入第一步，循环调用该滤波算法；若满足终止条件，则退出该方法。

本发明与现有技术相比的优点在于：

(1)设计了一种飞行器测试平台，验证陀螺仪多源噪声对姿控系统控制精度的影响，填补现有测试平台对陀螺仪多源噪声研究的空白；所设计的飞行器测试平台相比传统测试平台，可完成多种滤波方法的测试、分析与评估。

(2)针对陀螺仪多源噪声，设计一种针对陀螺仪多源噪声的混合滤波方法，该方法对陀螺仪器件噪声采用H₂优化抑制，对滤波系统建模误差采用H_∞优化抑制，结合H₂优化与H_∞优化设计混合滤波器，保证了估计误差满足给定的优化指标，提高姿态信息的滤波精度，改善控制性能。

附图说明

图1为本发明的测试平台工作回路；

图2为本发明混合滤波方法设计流程图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步详细说明。

以微纳三轴稳定卫星的姿控系统半物理仿真试验为例，来说明测试平台以及针对陀螺仪多源噪声的混合滤波方法。

如图1所示，给出了本发明的测试平台，其组成部分包括了实时仿真目标机1、三轴速率转台2、光纤陀螺仪3、星敏感器4、星仿真器5、姿态确定模块6、姿态控制模块7以及执行机构8；所述实时仿真目标机1接收执行机构8输出的力矩信号，同步运算微纳卫星的动力学与运动学模型，解算出微纳卫星滚动、俯仰以及偏航轴的角速度与角度信息；所述三轴速率转台2作为卫星的运动仿真器，可以提供与微纳卫星在轨时类似或等价的运动，它受实时仿真目标机1的控制，模拟了微纳卫星的在轨运动；所述光纤陀螺仪3与星敏感器4固定安装在三轴速率转台2的台面上，光纤陀螺仪3敏感三轴速率转台2的角速度，星敏感器4敏感恒星方位，为光纤陀螺仪3提供量测信息；所述星仿真器5是在试验条件下对恒星特征进行仿真，星敏感器提供模拟恒星方位；所述姿态确定模块6实时运行针对光纤陀螺仪多源噪声的混合滤波方法，实现对光纤陀螺仪3与星敏感器4实测的姿态信息滤波，包括了数字运算单元61、RS422/232通信单元62、电源单元63、复位单元64、扩展备用I/O口65、存储器66以及时钟单元67，其中RS422/232通信单元62、电源单元63、复位单元64、扩展备用I/O口65、存储器66、时钟单元67都与数字运算单元61相连，数字运算单元61是姿态确定模块6的运算核心，RS422/232通信单元62是姿态确定模块6的输入输出端，接收三轴陀螺仪3与星敏感器4实测姿态信息，输出滤波后的姿态信息，电源单元63负责供电，复位单元64可恢复数字运算单元61的初始状态，扩展备用I/O口65提供可扩展的输入输出端口，存储器66存储数字运算单元61的程序与数据，时钟单元67为数字运算单元61提供时钟频率；所述姿态控制模块7实时运算姿态控制方法，为执行机构8提供力矩控制指令；所述执行机构8包括了反作用飞轮81和磁力矩器82，在接收力矩控制指令后，将输出力矩信号传给实时仿真目标机1。

该测试平台的数据流回路如下：微纳卫星的期望姿态信号通过姿态控制模块7的RS232串行数据通信接口导入姿控系统，该期望姿态信号通过与姿态确定模块6滤波后的姿态信息进行比较，得到偏差信号；将此偏差信号传递给姿态控制模块7，该模块中的姿态控制算法解算出指令力矩信号并对指令力矩信号进行分配，分配后的信号经RS422/232串行数据通信接口作用于执行机构8；执行机构8中反作用飞轮81和磁力矩器82接收到分配后的指令力矩信号，联合作用输出力矩信号，该力矩信号通过RS422/232串行数据通信接口传入实时仿真目标机1中的飞行器姿态动力学仿真机11与飞行器姿态运动学仿真机12；飞行器姿态动力学仿真机11与飞行器姿态运动学仿真机12分别运算微纳卫星姿态动力学与运动学模块，得到微纳卫星转动角速度，该转动角速度经RS422/232串行数据通信接口作用于三轴速率转台2；三轴速率转台2受微纳卫星转动角速度驱动而转动，装配在三轴速率转台2上的光纤陀螺仪3以及星敏感器4将敏感到的姿态信息经RS422/232串行数据通信接口传递给姿态确定模块6；姿态确定模块6中针对光纤陀螺仪多源噪声的混合滤波方法解算得到滤波后的姿态信息，该姿态信息与期望姿态作比较后得到新的偏差信号，形成了测试平台的数据流回路。

如图2所示，给出的针对光纤陀螺仪多源噪声的混合滤波方法是一种递推算法，该方法的滤波周期包括以下步骤：首先构建含有光纤陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程；其次设计一个线性滤波器；然后构造增广的状态估计方程与待估计信号估计误差方程；最后是利用MATLAB中的LMI工具箱求解最优化参数；

第一步：设计含光纤陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程；

状态方程含光纤陀螺仪器件噪声与滤波系统建模误差，其状态变量是用四元数形式表示的微纳卫星姿态角；星敏感器的量测方程为状态方程提供修正信息，状态方程、量测方程以及待估计信号方程为Σ₁：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A\left(t\right)x\left(t\right)+B\left(t\right)n\left(t\right)+w_1\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Dn\left(t\right)+w_1\left(t\right)\\ z\left(t\right)=Lx\left(t\right)\end{array}\right.$

其中，t表示时间，x(t)＝[q₀(t) q₁(t) q₂(t) q₃(t)]^T为四元数形式表示的飞行器姿态角；为x(t)的一阶微分；n(t)＝[n_x(t) n_y(t) n_z(t)]^T为陀螺仪滤波系统的建模误差；w₁(t)＝[δ₀(t) δ₁(t) δ₂(t) δ₃(t)]^T为四元数形式表示的陀螺仪器件噪声；y(t)＝[y₀(t) y₁(t) y₂(t) y₃(t)]^T为四元数形式表示的星敏感器测得的恒星方位；z(t)＝[z₀(t) z₁(t) z₂(t) z₃(t)]^T为系统待估计的信号；式中，矩阵C是量测矩阵，矩阵D是量测方程建模误差系数矩阵，矩阵L是为待估计信号系数矩阵，矩阵C、矩阵D、矩阵L都为四阶单位阵，系统状态矩阵A(t)为状态方程建模误差系数矩阵B(t)为w_x(t)、w_y(t)、w_z(t)分别为微纳卫星相对轨道坐标系的滚动轴角速率、俯仰轴角速率及偏航轴角速率；

将含有光纤陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程离散化，可以得到如下模型Σ₂：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}x(k+1)=\stackrel{\‾}{A}\left(k\right)x\left(k\right)+\stackrel{\‾}{B}\left(k\right)n\left(k\right)+w_1\left(k\right)\\ y\left(k\right)=Cx\left(k\right)+Dn\left(k\right)+w_1\left(k\right)\\ z\left(k\right)=Lx\left(k\right)\end{array}\right.$

其中x(k+1)为k+1时刻以四元数形式表示的状态变量，x(k)为k时刻以四元数形式表示的状态变量，n(k)为k时刻的陀螺仪滤波系统的建模误差，w₁(k)为k时刻以四元数形式表示的陀螺仪器件噪声，y(k)为k时刻以四元数形式表示的星敏感器测得的恒星方位，z(k)为k时刻待估计的信号；A(k)为k时刻的系统状态矩阵，B(k)为k时刻的状态方程建模误差系数矩阵，I为单位阵，T为采样时间；

第二步：针对光纤陀螺仪器件噪声和滤波系统建模误差设计线性滤波器；

$\left\{\begin{array}{c}\hat{x}(k+1)=\hat{A}\hat{x}\left(k\right)+\hat{B}y\left(k\right)\\ \hat{z}\left(k\right)=\hat{C}\hat{x}\left(k\right)\end{array}\right.$

式中为k+1时刻的状态变量估计值，为k时刻的状态变量估计值，为k时刻的待估计信号估计值，为线性滤波器的系统矩阵，为线性滤波器的输入矩阵，为线性滤波器的输出矩阵；

第三步：构造增广的状态估计方程与待估计信号估计误差方程；

目的是使得增广的状态估计方程与待估计信号估计误差方程：

$\left\{\begin{array}{c}\tilde{x}(k+1)=\tilde{A}\tilde{x}\left(k\right)+\tilde{B}n\left(k\right)+{\tilde{B}}_1{w}_1\left(k\right)\\ \tilde{z}\left(k\right)=\tilde{C}\tilde{x}\left(k\right)\end{array}\right.$

式中第一个表达式为增广的状态估计方程，第二个方程是待估计信号的估计误差方程，这两个方程是渐近稳定的，且对应于通道的滤波误差方差的一个上界最小化，对应于通道的滤波误差向量满足其中，

$\tilde{x}\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ \hat{x}\left(k\right)\end{array}\right],\tilde{A}=\left[\begin{array}{cc}\stackrel{\‾}{A}\left(k\right)& 0\\ \hat{B}C& \hat{A}\end{array}\right],\hat{B}=\left[\begin{array}{c}\stackrel{\‾}{B}\left(k\right)\\ \hat{B}D\end{array}\right],{\tilde{B}}_1=\left[\begin{array}{c}I\\ \hat{B}\end{array}\right]$

$\tilde{C}=\left[\begin{array}{cc}L& -\hat{C}\end{array}\right],\tilde{z}\left(k\right)=z\left(k\right)-\hat{z}\left(k\right)$

第四步：利用MATLAB中的LMI工具箱求解最优化参数；

对给定的常数γ>0，系统Σ₂存在一个混合滤波器可以归纳为以下的优化问题：

$\underset{R,P,Q,S,U,H}{\min }Trace\left(H\right)$

$s.t.\left(i\right)\left[\begin{array}{cccccc}R& *& *& *& *& *\\ R& P& *& *& *& *\\ 0& 0& {\mathrm{\γ}}^{2}I& *& *& *\\ R\stackrel{\‾}{A}\left(k\right)& R\stackrel{\‾}{A}\left(k\right)& R& R& *& *\\ P\stackrel{\‾}{A}\left(k\right)+QC+S& P\stackrel{\‾}{A}\left(k\right)+QC& P+Q& R& P& *\\ L-U& L& 0& 0& 0& I\end{array}\right]>0$

$\left(ii\right)\left[\begin{array}{ccc}H& *& *\\ R\stackrel{\‾}{B}\left(k\right)& R& *\\ P\stackrel{\‾}{B}\left(k\right)+QD& R& P\end{array}\right]\≥0$

其中*处的子块可由矩阵的对称性得到；进而求解出优化问题的最优解R、P、S、Q、U、H，其中R、P、S、Q、U、H为求解最优化问题定义的辅助变量，Trace表示矩阵的迹，即为矩阵对角元素之和，则：

$\hat{A}={(R-P)}^{-1}S,\hat{B}={(R-P)}^{-1}Q,\hat{C}=U$

是所要求的混合滤波器的系数矩阵，相应的

由于优化问题是一个具有线性矩阵不等式约束和线性目标函数的凸优化问题，因此可以应用MATLAB中LMI工具箱的最优化极小值求解器mincx来求解该问题。

第五步：判断终止条件，若满足，则滤波方法运行结束，若不满足，则循环调用；

终止条件为滤波方法的总运行时间，当该方法运行时间没有达到总运行时间，则更新时间，以下一个离散时刻k+1代入第一步，循环调用该滤波算法；若满足终止条件，则退出该方法。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (2)

1.一种针对陀螺仪多源噪声的混合滤波方法与测试平台，其特征在于：测试平台包括实时仿真目标机(1)、三轴速率转台(2)、三轴陀螺仪(3)、星敏感器(4)、星仿真器(5)、姿态确定模块(6)、姿态控制模块(7)以及执行机构(8)；所述实时仿真目标机(1)包括飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)，分别实时运算飞行器的姿态动力学与运动学模型；所述三轴速率转台(2)是飞行器的运动仿真器，提供与飞行器相似或等价的运动，用于模拟飞行器的姿态运动；所述三轴陀螺仪(3)是角速度姿态信息敏感器，敏感三轴速率转台(2)的角速度；所述星敏感器(4)通过敏感恒星方位而获得卫星姿态信息，用来修正三轴陀螺仪(3)的漂移与噪声；所述星仿真器(5)为星敏感器(4)提供模拟恒星方位；所述姿态确定模块(6)实时运行针对陀螺仪多源噪声的混合滤波方法，对三轴陀螺仪(3)与星敏感器(4)实测的姿态信息滤波，包括了数字运算单元(61)、RS422/232通信单元(62)、电源单元(63)、复位单元(64)、扩展备用I/O口(65)、存储器(66)以及时钟单元(67)，其中RS422/232通信单元(62)、电源单元(63)、复位单元(64)、扩展备用I/O口(65)、存储器(66)、时钟单元(67)都与数字运算单元(61)相连，数字运算单元(61)是姿态确定模块(6)的运算核心，RS422/232通信单元(62)是姿态确定模块(6)的输入输出端，接收三轴陀螺仪(3)与星敏感器(4)实测姿态信息，输出滤波后的姿态信息，电源单元(63)负责供电，复位单元(64)可恢复数字运算单元(61)的初始状态，扩展备用I/O口(65)提供可扩展的输入输出端口，存储器(66)存储数字运算单元(61)的程序与数据，时钟单元(67)为数字运算单元(61)提供时钟频率；所述姿态控制模块(7)实时运算姿态控制算法，为执行机构(8)提供力矩控制指令；所述执行机构(8)包括了反作用飞轮(81)和磁力矩器(82)，在接收力矩控制指令后，将输出力矩信号传给实时仿真目标机(1)；测试平台的数据流回路如下：飞行器的期望姿态信号通过姿态控制模块(7)的RS232串行数据通信接口导入测试平台，该期望姿态信号通过与姿态确定模块(6)滤波后的姿态信息进行比较，得到偏差信号；将此偏差信号传递给姿态控制模块(7)，该模块中的姿态控制算法解算出指令力矩信号并对指令力矩信号进行分配，分配后的信号经RS422/232串行数据通信接口作用于执行机构(8)；执行机构(8)中反作用飞轮(81)和磁力矩器(82)接收到分配后的指令力矩信号，联合作用输出力矩信号，该力矩信号通过RS422/232串行数据通信接口传入实时仿真目标机(1)中的飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)；飞行器姿态动力学仿真机(11)与飞行器姿态运动学仿真机(12)分别运算飞行器姿态动力学与运动学模块，得到飞行器转动角速度，该飞行器转动角速度经RS422/232串行数据通信接口传送至三轴速率转台(2)；三轴速率转台(2)受飞行器转动角速度驱动而转动，装配在三轴速率转台(2)上的三轴陀螺仪(3)以及星敏感器(4)将敏感到的姿态信息经RS422/232串行数据通信接口传递给姿态确定模块(6)；姿态确定模块(6)中针对陀螺仪多源噪声的混合滤波方法解算得到滤波后的姿态信息，该姿态信息与期望姿态作比较后得到新的偏差信号，形成了测试平台的数据流回路。

2.根据权利要求1所述的针对陀螺仪多源噪声的混合滤波方法与测试平台，其特征在于：所述的混合滤波方法是一种递推算法，该方法的滤波周期包括以下步骤：首先构建含有陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程；其次设计一个线性滤波器；然后构造增广的状态估计方程与待估计信号估计误差方程；最后是利用MATLAB中的LMI工具箱求解最优化参数；

第一步：设计含陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程状态方程含陀螺仪器件噪声与滤波系统建模误差，其状态变量是用四元数形式表示的飞行器姿态角；星敏感器的量测方程为状态方程提供修正信息，状态方程、量测方程以及待估计信号方程为Σ₁：

${\mathrm{\Σ}}_1:\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=A\left(t\right)x\left(t\right)+B\left(t\right)n\left(t\right)+w_1\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Dn\left(t\right)+w_1\left(t\right)\\ z\left(t\right)=Lx\left(t\right)\end{array}\right.$

其中，t表示时间，x(t)＝[q₀(t) q₁(t) q₂(t) q₃(t)]^T为四元数形式表示的飞行器姿态角；为x(t)的一阶微分；n(t)＝[n_x(t) n_y(t) n_z(t)]^T为陀螺仪滤波系统的建模误差；w₁(t)＝[δ₀(t) δ₁(t) δ₂(t) δ₃(t)]^T为四元数形式表示的陀螺仪器件噪声；y(t)＝[y₀(t) y₁(t) y₂(t) y₃(t)]^T为四元数形式表示的星敏感器测得的恒星方位；z(t)＝[z₀(t) z₁(t) z₂(t) z₃(t)]^T为系统待估计的信号；式中，矩阵C是量测矩阵，矩阵D是量测方程建模误差系数矩阵，矩阵L是为待估计信号系数矩阵，它们矩阵C、矩阵D、矩阵L都为四阶单位阵，系统状态矩阵A(t)为状态方程建模误差系数矩阵B(t)为

w_x(t)、w_y(t)、w_z(t)分别为飞行器相对轨道坐标系的滚动轴角速率、俯仰轴角速率及偏航轴角速率；

将含有陀螺仪器件噪声和滤波系统建模误差的状态方程与星敏感器的量测方程离散化，可以得到如下模型Σ₂：

${\mathrm{\Σ}}_2:\left\{\begin{array}{c}x(k+1)=\stackrel{\‾}{A}\left(k\right)x\left(k\right)+\stackrel{\‾}{B}\left(k\right)n\left(k\right)+w_1\left(k\right)\\ y\left(k\right)=Cx\left(k\right)+Dn\left(k\right)+w_1\left(k\right)\\ z\left(k\right)=Lx\left(k\right)\end{array}\right.$

其中x(k+1)为k+1时刻以四元数形式表示的状态变量，x(k)为k时刻以四元数形式表示的状态变量，n(k)为k时刻的陀螺仪滤波系统的建模误差，w₁(k)为k时刻以四元数形式表示的陀螺仪器件噪声，y(k)为k时刻以四元数形式表示的星敏感器测得的恒星方位，z(k)为k时刻待估计的信号；A(k)为k时刻的系统状态矩阵，B(k)为k时刻的状态方程建模误差系数矩阵，I为单位阵，T为采样时间；

第二步：针对陀螺仪器件噪声和滤波系统建模误差设计线性滤波器

$\left\{\begin{array}{c}\hat{x}(k+1)=\hat{A}\hat{x}\left(k\right)+\hat{B}y\left(k\right)\\ \hat{z}\left(k\right)=\hat{C}\hat{x}\left(k\right)\end{array}\right.$

式中为k+1时刻的状态变量估计值，为k时刻的状态变量估计值，为k时刻的待估计信号估计值，为线性滤波器的系统矩阵，为线性滤波器的输入矩阵，为线性滤波器的输出矩阵；

第三步：构造增广的状态估计方程与待估计信号估计误差方程

目的是使得增广的状态估计方程与待估计信号估计误差方程为：

$\left\{\begin{array}{c}\tilde{x}(k+1)=\tilde{A}\tilde{x}\left(k\right)+\tilde{B}n\left(k\right)+{\tilde{B}}_1{w}_1\left(k\right)\\ \tilde{z}\left(k\right)=\tilde{C}\tilde{x}\left(k\right)\end{array}\right.$

式中第一个表达式为增广的状态估计方程，第二个方程是待估计信号的估计误差方程，这两个方程是渐近稳定的，且对应于通道的滤波误差方差的一个上界最小化，对应于通道的滤波误差向量满足其中，

$\tilde{x}\left(k\right)=\left[\begin{array}{c}x\left(k\right)\\ \hat{x}\left(k\right)\end{array}\right],\tilde{A}=\left[\begin{array}{cc}\stackrel{\‾}{A}\left(k\right)& 0\\ \hat{B}C& \hat{A}\end{array}\right],\tilde{B}=\left[\begin{array}{c}\stackrel{\‾}{B}\left(k\right)\\ \hat{B}D\end{array}\right],{\tilde{B}}_1=\left[\begin{array}{c}I\\ \hat{B}\end{array}\right]$

$\tilde{C}=\[\begin{array}{cc}L& -\hat{C}\end{array}\],\tilde{z}\left(k\right)=z\left(k\right)-\hat{z}\left(k\right)$

第四步：利用MATLAB中的LMI工具箱求解最优化参数

对给定的常数γ＞0，系统Σ₂存在一个混合滤波器可以归纳为以下的优化问题：

$\underset{R,P,Q,S,U,H}{min}Trace\left(H\right)$

$\begin{array}{cc}s.t.\left(i\right)& \left[\begin{array}{cccccc}R& *& *& *& *& *\\ R& P& *& *& *& *\\ 0& 0& {\mathrm{\γ}}^{2}I& *& *& *\\ R\stackrel{\‾}{A}\left(k\right)& R\stackrel{\‾}{A}\left(k\right)& R& R& *& *\\ P\stackrel{\‾}{A}\left(k\right)+QC+S& P\stackrel{\‾}{A}\left(k\right)+QC& P+Q& R& P& *\\ L-U& L& 0& 0& 0& I\end{array}\right]>0\end{array}$

$\begin{array}{cc}\left(ii\right)& \left[\begin{array}{ccc}H& *& *\\ R\stackrel{\‾}{B}\left(k\right)& R& *\\ P\stackrel{\‾}{B}\left(k\right)+QD& R& P\end{array}\right]\≥0\end{array}$

其中*处的子块可由矩阵的对称性得到；进而求解出优化问题的最优解R、P、S、Q、U、H，其中R、P、S、Q、U、H为求解最优化问题定义的辅助变量，Trace表示矩阵的迹，即为矩阵对角元素之和，则：

$\hat{A}={(R-P)}^{-1}S,\hat{B}={(R-P)}^{-1}Q,\hat{C}=U$

是所要求的混合滤波器的系数矩阵，相应的

由于优化问题是一个具有线性矩阵不等式约束和线性目标函数的凸优化问题，因此可以应用MATLAB中LMI工具箱的求解器最优化最小值mincx来求解该问题；

第五步：判断终止条件，若满足，则滤波方法运行结束，若不满足，则循环调用；

终止条件为滤波方法的总运行时间，当该方法运行时间没有达到总运行时间，则更新时间，以下一个离散时刻k+1代入第一步，循环调用该滤波算法；若满足终止条件，则退出该方法。