发明名称 一种基于二阶量测更新的低成本姿态估计方法
摘要 本发明提供的是一种基于二阶量测更新的低成本姿态估计方法。本发明通过对三轴微机械陀螺、三轴微机械加速度计和三轴磁强计的输出数据,进行滤波处理得到载体的姿态信息。针对在室内或磁干扰较强场所,磁强计输出会使横摇和纵摇误差变大,传统方法难以解决的问题。本方法在滤波的量测更新阶段,创新性地采用二阶量测更新,即先进行加速度计量测更新再进行磁强计量测更新。以此修正标准量测更新算法从而使磁强计更新只影响方位角。利用本方法可以使用低成本的微惯性测量单元和磁强计进行姿态估计,并且估计精度高、实时性好、适应强磁干扰环境。本方法适用于车辆、无人机和船舰等载体的姿态估计。
申请公布号 CN103822633B 申请公布日期 2016.12.07
申请号 CN201410047878.0 申请日期 2014.02.11
申请人 哈尔滨工程大学 发明人 周广涛;孙艳涛;姜鑫;梁宏;时贵敏;赵博;郝勤顺;李佳璇;林萌萌;张丽丽
分类号 G01C21/20(2006.01)I 主分类号 G01C21/20(2006.01)I
代理机构 代理人
主权项 一种基于二阶量测更新的低成本姿态估计方法,其特征是:(1)在k=0时刻,即初始时刻由初始对准得到初始姿态φ<sub>0</sub>,所述初始姿态包括横摇、纵摇和航向;(2)在k≥1时刻,利用微机械陀螺输出的载体三轴的角速率数据,计算出姿态角的变化值Δφ,由式<img file="FDA0000970844220000011.GIF" wi="322" he="95" />计算出滤波开始前的粗略姿态角;(3)采集微机械加速度计和磁力计输出的载体三轴线的加速度和磁力数据;(4)设滤波的状态变量为x=[q<sub>e</sub> b<sub>g</sub> b<sub>a</sub>]<sup>T</sup>∈R<sup>9×1</sup>;式中:q<sub>e</sub>表示四元数误差<img file="FDA0000970844220000012.GIF" wi="49" he="55" />矢量部分,b<sub>g</sub>表示陀螺漂移,b<sub>a</sub>表示加速度计零偏;假设四元数误差<img file="FDA0000970844220000013.GIF" wi="259" he="63" />根据四元数相关性质建立状态方程:<maths num="0001"><math><![CDATA[<mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mover><mi>q</mi><mo>&CenterDot;</mo></mover><mi>e</mi></msub></mtd></mtr><mtr><mtd><msub><mover><mi>b</mi><mo>&CenterDot;</mo></mover><mi>g</mi></msub></mtd></mtr><mtr><mtd><msub><mover><mi>b</mi><mo>&CenterDot;</mo></mover><mi>a</mi></msub></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><mo>-</mo><mo>&lsqb;</mo><msub><mi>y</mi><mi>g</mi></msub><mo>&times;</mo><mo>&rsqb;</mo></mrow></mtd><mtd><mrow><mo>-</mo><mn>0.5</mn><mi>I</mi></mrow></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd></mtr></mtable></mfenced><mfenced open = "[" close = "]"><mtable><mtr><mtd><msub><mi>q</mi><mi>e</mi></msub></mtd></mtr><mtr><mtd><msub><mi>b</mi><mi>g</mi></msub></mtd></mtr><mtr><mtd><msub><mi>b</mi><mi>a</mi></msub></mtd></mtr></mtable></mfenced><mo>+</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><mo>-</mo><mn>0.5</mn><msub><mi>v</mi><mi>g</mi></msub></mtd></mtr><mtr><mtd><msub><mi>v</mi><msub><mi>b</mi><mi>g</mi></msub></msub></mtd></mtr><mtr><mtd><msub><mi>v</mi><msub><mi>b</mi><mi>a</mi></msub></msub></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA0000970844220000014.GIF" wi="933" he="294" /></maths>式中:<img file="FDA0000970844220000015.GIF" wi="622" he="231" />下标g代表陀螺、a代表加速度计、m代表磁力计,y<sub>g</sub>表示陀螺实际输出,传感器噪声v<sub>a</sub>、v<sub>g</sub>以及v<sub>m</sub>为零均值高斯白噪声;取观测量为<img file="FDA0000970844220000016.GIF" wi="382" he="143" />由于q<sub>e</sub>为小量,再忽略二阶小项,建立量测方程:<maths num="0002"><math><![CDATA[<mrow><msub><mi>y</mi><mi>a</mi></msub><mo>-</mo><mi>C</mi><mrow><mo>(</mo><mover><mi>q</mi><mo>^</mo></mover><mo>)</mo></mrow><mover><mi>g</mi><mo>~</mo></mover><mo>=</mo><mn>2</mn><mo>&lsqb;</mo><mi>C</mi><mrow><mo>(</mo><mover><mi>q</mi><mo>^</mo></mover><mo>)</mo></mrow><mover><mi>g</mi><mo>~</mo></mover><mo>&times;</mo><mo>&rsqb;</mo><msub><mi>q</mi><mi>e</mi></msub><mo>+</mo><msub><mi>v</mi><mi>a</mi></msub><mo>+</mo><msub><mi>b</mi><mi>a</mi></msub></mrow>]]></math><img file="FDA0000970844220000017.GIF" wi="846" he="71" /></maths><maths num="0003"><math><![CDATA[<mrow><msub><mi>y</mi><mi>m</mi></msub><mo>-</mo><mi>C</mi><mrow><mo>(</mo><mover><mi>q</mi><mo>^</mo></mover><mo>)</mo></mrow><mover><mi>m</mi><mo>~</mo></mover><mo>=</mo><mn>2</mn><mo>&lsqb;</mo><mi>C</mi><mrow><mo>(</mo><mover><mi>q</mi><mo>^</mo></mover><mo>)</mo></mrow><mover><mi>m</mi><mo>~</mo></mover><mo>&times;</mo><mo>&rsqb;</mo><msub><mi>q</mi><mi>e</mi></msub><mo>+</mo><msub><mi>v</mi><mi>m</mi></msub></mrow>]]></math><img file="FDA0000970844220000018.GIF" wi="773" he="71" /></maths>式中:<img file="FDA0000970844220000019.GIF" wi="101" he="55" />表示由估计四元数组成的姿态矩阵,y<sub>a</sub>和y<sub>m</sub>分别为加速度计和磁强计输出,定义<img file="FDA00009708442200000110.GIF" wi="267" he="63" />和<img file="FDA00009708442200000111.GIF" wi="488" he="62" />g为重力加速度,α为磁倾角;(5)利用步骤(4)中建立的状态方程进行滤波的时间更新,由下式<maths num="0004"><math><![CDATA[<mrow><msubsup><mover><mi>x</mi><mo>^</mo></mover><mi>k</mi><mo>-</mo></msubsup><mo>=</mo><msub><mi>&Phi;</mi><mi>k</mi></msub><msub><mover><mi>x</mi><mo>^</mo></mover><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub></mrow>]]></math><img file="FDA00009708442200000112.GIF" wi="318" he="79" /></maths><maths num="0005"><math><![CDATA[<mrow><msubsup><mi>P</mi><mi>k</mi><mo>-</mo></msubsup><mo>=</mo><msub><mi>&Phi;</mi><mi>k</mi></msub><msub><mi>P</mi><mrow><mi>k</mi><mo>-</mo><mn>1</mn></mrow></msub><msup><msub><mi>&Phi;</mi><mi>k</mi></msub><mi>T</mi></msup><mo>+</mo><msub><mi>Q</mi><mi>k</mi></msub></mrow>]]></math><img file="FDA00009708442200000115.GIF" wi="533" he="79" /></maths>分别计算出一步预测量<img file="FDA00009708442200000113.GIF" wi="51" he="56" />和状态误差协方差阵<img file="FDA00009708442200000114.GIF" wi="83" he="55" />式中:Φ<sub>k</sub>为系统状态转移矩阵,Q<sub>k</sub>为系统噪声协方差阵;(6)在滤波的量测更新阶段采用二阶量测更新,先进行加速度计量测更新;取此阶段观测量<img file="FDA0000970844220000021.GIF" wi="339" he="63" />根据建立的量测方程有<img file="FDA0000970844220000022.GIF" wi="481" he="71" />其更新过程如下:先计算滤波增益:<img file="FDA0000970844220000023.GIF" wi="693" he="71" />然后计算状态估计量:<img file="FDA0000970844220000024.GIF" wi="607" he="70" />再计算状态误差协方差阵:<img file="FDA0000970844220000025.GIF" wi="1014" he="71" />式中:下标a表示加速度计量测更新阶段,R<sub>a</sub>为加速度计噪声协方差阵此时<img file="FDA0000970844220000026.GIF" wi="51" he="62" />的9个状态都会更新,利用<img file="FDA0000970844220000027.GIF" wi="54" he="62" />更新<img file="FDA0000970844220000028.GIF" wi="129" he="63" />用<img file="FDA0000970844220000029.GIF" wi="74" he="61" />代表已更新状态,取<img file="FDA00009708442200000210.GIF" wi="275" he="63" />即前三个变量值,利用式<img file="FDA00009708442200000211.GIF" wi="219" he="63" />进行校正,然后进行四元数规范化:<maths num="0006"><math><![CDATA[<mrow><mover><mi>q</mi><mo>^</mo></mover><mo>&DoubleLeftArrow;</mo><mfrac><mover><mi>q</mi><mo>^</mo></mover><mrow><mo>|</mo><mo>|</mo><mover><mi>q</mi><mo>^</mo></mover><mo>|</mo><mo>|</mo></mrow></mfrac></mrow>]]></math><img file="FDA00009708442200000212.GIF" wi="171" he="134" /></maths>加速度计更新结束后,将<img file="FDA00009708442200000213.GIF" wi="70" he="63" />的q<sub>e</sub>部分设置为零;(7)利用步骤(6)中得到的<img file="FDA00009708442200000214.GIF" wi="70" he="63" />进行磁力计量测更新;取此阶段观测量为<img file="FDA00009708442200000215.GIF" wi="434" he="63" />则得<img file="FDA00009708442200000216.GIF" wi="571" he="71" />其更新过程如下:取此时状态误差协方差阵为<maths num="0007"><math><![CDATA[<mrow><msubsup><mi>P</mi><mrow><mi>m</mi><mo>,</mo><mi>k</mi></mrow><mo>-</mo></msubsup><mo>=</mo><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><msub><mi>P</mi><mrow><mi>a</mi><mo>,</mo><mi>k</mi></mrow></msub><mrow><mo>(</mo><mn>1</mn><mo>:</mo><mn>3</mn><mo>,</mo><mn>1</mn><mo>:</mo><mn>3</mn><mo>)</mo></mrow></mrow></mtd><mtd><msub><mn>0</mn><mrow><mn>3</mn><mo>&times;</mo><mn>6</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mn>0</mn><mrow><mn>6</mn><mo>&times;</mo><mn>3</mn></mrow></msub></mtd><mtd><msub><mn>0</mn><mrow><mn>6</mn><mo>&times;</mo><mn>6</mn></mrow></msub></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA00009708442200000217.GIF" wi="590" he="182" /></maths>由此计算滤波增益:<img file="FDA00009708442200000218.GIF" wi="1007" he="143" />然后计算状态估计量:<img file="FDA00009708442200000219.GIF" wi="717" he="87" />再计算状态误差协方差阵:<img file="FDA00009708442200000220.GIF" wi="1070" he="86" />式中:下标m表示磁力计量测更新阶段,P<sub>a,k</sub>(1:3,1:3)表示由P<sub>a,k</sub>前三行前三列所组成的矩阵,0<sub>mn</sub>为m行n列的零矩阵并且<img file="FDA00009708442200000221.GIF" wi="381" he="71" />(8)根据滤波计算出的四元数误差,校正步骤(2)中的粗略姿态角,得到精确的姿态角。
地址 150001 黑龙江省哈尔滨市南岗区南通大街145号哈尔滨工程大学科技处知识产权办公室