CN114877858A - 一种高动态和磁干扰环境下的姿态估计算法 - Google Patents

一种高动态和磁干扰环境下的姿态估计算法 Download PDF

Info

Publication number
CN114877858A
CN114877858A CN202210486754.7A CN202210486754A CN114877858A CN 114877858 A CN114877858 A CN 114877858A CN 202210486754 A CN202210486754 A CN 202210486754A CN 114877858 A CN114877858 A CN 114877858A
Authority
CN
China
Prior art keywords
state
formula
magnetic interference
equation
expressed
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210486754.7A
Other languages
English (en)
Other versions
CN114877858B (zh
Inventor
孙伟
韩冬
陈龙
连杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Xian Research Institute Co Ltd of CCTEG
Original Assignee
Xidian University
Xian Research Institute Co Ltd of CCTEG
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University, Xian Research Institute Co Ltd of CCTEG filed Critical Xidian University
Priority to CN202210486754.7A priority Critical patent/CN114877858B/zh
Publication of CN114877858A publication Critical patent/CN114877858A/zh
Application granted granted Critical
Publication of CN114877858B publication Critical patent/CN114877858B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C1/00Measuring angles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/18Stabilised platforms, e.g. by gyroscope
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/16Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/18Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Automation & Control Theory (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computing Systems (AREA)
  • Navigation (AREA)

Abstract

本发明提供一种高动态和磁干扰环境下的姿态估计算法,该姿态估计算法包括第一层姿态估计算法、第二层姿态估计算法,其中第一层算法解算俯仰、滚转角,选用方向余弦矩阵的最后一行作为俯仰、滚转角的状态变量,选用该状态变量的误差状态在EKF框架下进行更新与校正;第二层算法解算航向角,选用方向余弦矩阵的第一列的前两个元素作为航向角的状态变量,同时为了避免磁干扰对姿态解算的影响,使用磁力计模长乘以模块磁场倾角增量再取对数函数的方法表征磁干扰的强度,并将磁干扰的强度和观测噪声相关联。解决了姿态估计算法在高动态环境下精度低、易发散,且姿态角受磁干扰影响的问题,有效提高了姿态估计的精度和鲁棒性。

Description

一种高动态和磁干扰环境下的姿态估计算法
【技术领域】
本发明涉及一种姿态估计算法,具体涉及一种高动态和磁干扰环境下的姿态估计算法。
【背景技术】
卡尔曼滤波是姿态估计中应用最为广泛的一种方法,传统卡尔曼滤波只使用于线性系统,而三维旋转运动由于其转动空间是流形,无论用哪种其数学描述方法都是在非线性空间中工作。为了适应非线性系统,提出使用扩展卡尔曼滤波(EKF)进行姿态估计,其核心思想就是利用泰勒级数,将非线性系统在工作点附近进行一阶泰勒展开并忽略高阶项。但是在高动态环境中,运动方程的模型噪声和陀螺仪噪声增大,观测方程的线性化点开始远离工作点,最终导致EKF的估计结果产生偏差或不一致,甚至会导致滤波器发散。
由于加速度计并未包含航向信息,要想获得无累计误差的航向角则需要引入磁力计,但载体在运动过程中,不可避免的会运动到有外部磁场干扰的环境,当距离干扰源比较近时干扰磁场的强度会远远大于地磁场强度,使得磁力计数据失去参考价值。
【发明内容】
为了解决上述问题,本发明针对传统的姿态估计算法在高动态环境下估计值容易发散,姿态角受磁干扰影响的问题,提出一种高动态和磁干扰环境下的姿态估计算法,其将俯仰角、滚转角和航向角分层解算,解决了姿态估计算法在高动态环境下精度低、易发散,且姿态角受磁干扰影响的问题,有效提高了姿态估计的精度和鲁棒性。
本发明是通过以下技术方案实现的,提供一种高动态和磁干扰环境下的姿态估计算法,包括以下步骤:
S1第一层姿态估计算法:基于线性化程度更高的误差状态解算俯仰角、滚转角,并选用方向余弦矩阵的最后一行作为俯仰角、滚转角状态变量,选用该状态变量的误差状态在EKF框架下进行更新与校正;
S2第二层姿态估计算法:选用方向余弦矩阵的第一列的前两个元素作为航向角的状态变量,并选用磁力计校正航向信息,使用磁力计模值和模块磁场倾角增量相乘再取对数函数的方法度量磁干扰强度。
特别的,所述S1具体按照以下步骤实施:
S11名义状态预测
选用从载体坐标系b系到导航坐标系n系的姿态旋转矩阵
Figure BDA0003630294470000021
的最后一行r作为表达俯仰θ、滚转角γ的状态变量,并将陀螺仪的零偏bg当成状态变量估计,第一层的名义状态向量x1按如下公式表示:
x1=[r bg] (1),
公式(1)中,r按如下公式计算:
r=[-sinθ sinγcosθ cosγcosθ] (2),
状态向量x1的运动方程按如下公式表示:
Figure BDA0003630294470000031
公式(3)中,bΦk为[tk-1,tk]时间段内的等效旋转矢量,(bΦk)^bΦk的反对称矩阵;
S12误差状态预测和校正
误差状态的运动方程为:
Figure BDA0003630294470000032
Figure BDA0003630294470000033
公式(4)中,δφ为r的误差状态,δbg为bg的误差状态,将公式(4)按矩阵形式表示为:
Figure BDA0003630294470000034
公式(5)中,Fk-1为误差状态δx1的状态转移矩阵,其按照如下公式计算:
Figure BDA0003630294470000035
误差状态的先验协方差
Figure BDA0003630294470000036
的预测方程如下:
Figure BDA0003630294470000037
于公式(7)中,Q为预测方程的过程噪声矩阵;
引入加速度计对预测结果进行校正,观测方程按如下公式表示:
z1,k=h1(x1,k,vk) (8),
公式(8)中,z1,k=am为加速度计测量值,vk为加速度计测量噪声,由于状态变量是选用姿态旋转矩阵的最后一行,所以h1(·)是一个线性函数,
根据链式求导法则,观测方程关于δx1的雅克比矩阵如下:
Figure BDA0003630294470000041
公式(9)中,Hx1=[I3 03],
Figure BDA0003630294470000042
因此,误差状态的后验校正过程按如下公式表示:
Figure BDA0003630294470000043
Figure BDA0003630294470000044
Figure BDA0003630294470000045
公式(10)中,R1为观测方程的噪声矩阵,I6为6*6的单位矩阵;
S13误差状态的注入和复位
将经过加速度计校正后的误差状态按如下公式注入到名义状态中:
Figure BDA0003630294470000046
Figure BDA0003630294470000047
公式(11)中,Exp(δφ)是误差状态到上SO(3)的指数映射,其按照如下公式表示:
Exp(δφ)=I3+u^sinφ+(1-cosφ)(u^)2 (12),
将误差状态注入到名义状态之后,需要将误差状态清零,具体按如下公式操作:
Figure BDA0003630294470000051
Figure BDA0003630294470000052
特别的,所述S11中,需引入单子样+前一周期的等效旋转矢量补偿不可交换误差,具体按如下公式表示:
Figure BDA0003630294470000053
公式(14)中,Φk为等效旋转矢量,ΔΘk为角增量,ΔΘk=ΔT·(bωm,k-bg,k-1)^
特别的,所述S2具体按照以下步骤实施:
S21航向信息解算
选用从载体坐标系b系到导航坐标系n系的姿态旋转矩阵
Figure BDA0003630294470000054
的第一列的前两个元素
Figure BDA0003630294470000055
作为表达航向角
Figure BDA0003630294470000056
的状态变量,其中第二层的状态向量x2为:
Figure BDA0003630294470000057
运动方程按如下公式表示:
Figure BDA0003630294470000061
公式(16)中,
Figure BDA0003630294470000062
Figure BDA0003630294470000063
Figure BDA0003630294470000064
计算所得,
状态x2的先验协方差的预测方程按如下公式表示:
Figure BDA0003630294470000065
S22航向信息校正
由于磁力计数据并不能直接校正状态x2,因此需将其转换到和状态x2同一量纲,状态x2的后验校正过程按如下公式表示:
Figure BDA0003630294470000066
Figure BDA0003630294470000067
Figure BDA0003630294470000068
公式(18)中,
Figure BDA0003630294470000069
S23磁干扰的检测
使用模块磁场倾角增量和磁力计模值来评价磁干扰强度,具体评价函数按如下公式表示:
mlen=log(||bmm||sin(|dip-dipref|)) (19),
公式(19)中,dipref为无磁干扰情况下的模块磁场倾角,dip为加速度计输出矢量和磁力计输出矢量的夹角,dip按如下公式进行计算:
Figure BDA0003630294470000071
公式(20)中,am为加速度计的测量值,mm为磁力计的测量值。
特别的,所述S2中,为了减小磁干扰对航向信息的影响,按如下公式将磁干扰强度和观测矩阵R2关联:
Figure BDA0003630294470000072
公式(21)中,σm为磁力计噪声的标准差,σd为磁干扰的基础标准差。
本发明提供一种高动态和磁干扰环境下的姿态估计算法,其将俯仰、滚转角的解算基于线性化程度更高的误差状态进行,并选用方向余弦矩阵的最后一行作为其状态变量,减小了系统非线性;同时,磁力计只用来校正航向信息,彻底避免了磁干扰对俯仰、滚转角的影响;使用模块磁场倾角增量和磁场模值相乘再取对数函数的方法度量磁干扰强度,并用其动态调节观测噪声的大小,使得航向角在短时间的磁干扰环境中也不受影响。
【附图说明】
图1为本发明一种高动态和磁干扰环境下的姿态估计算法的计算框图;
图2为磁干扰检测图;
图3为本发明与EKF算法、ESKF算法在磁干扰环境下的对比结果图;
图4为本发明与EKF算法、ESKF算法在高动态环境下的对比结果图。
【具体实施方式】
为了使本发明的目的、技术方案及优点更加清楚明白,以下对本发明进一步详细说明。
请参阅图1,本发明针对传统的姿态估计算法在高动态环境下估计值容易发散,姿态角受磁干扰影响的问题,提出一种基于误差状态的双层解耦姿态估计算法,即一种高动态和磁干扰环境下的姿态估计算法,具体按照如下步骤实施:
S1第一层姿态估计算法
于本发明中,第一层姿态估计算法只对俯仰角θ、滚转角γ进行估计,且只用加速度计对预测值做校正,将磁力计和俯仰角、滚转角彻底解耦,从根本上避免了磁干扰对俯仰角、滚转角的影。为了使本发明提供的估计算法适应高动态环境,使用线性化程度更高的误差状态做更新和校正。将系统真实状态分解成名义状态(nominal-state)和误差状态(error-state),名义状态是大的状态量继承了系统的非线性,但不考虑噪声,如此就不包含不确定度,误差状态是真实状态和名义状态之差,包含了传感器噪声和模型噪声。
S11名义状态预测
选用从载体坐标系(b系)到导航坐标系(n系)的姿态旋转矩阵
Figure BDA0003630294470000081
的最后一行r作为表达俯仰、滚转角的状态变量,并将陀螺仪的零偏bg也当成状态变量估计。
其中,第一层的名义状态向量x1按如下公式表示:
x1=[r bg] (1),
公式(1)中,r按如下公式进行计算:
r=[-sinθ sinγcosθ cosγcosθ] (2);
状态向量x1的运动方程按如下公式表示:
Figure BDA0003630294470000091
公式(3)中,bΦk为[tk-1,tk]时间段内的等效旋转矢量,(bΦk)^bΦk的反对称矩阵,在此处使用等效旋转矢量的目的是:在高动态的环境下载体不满足定轴转动的假设,若直接假设角增量ΔΘk≈ΔT·[bωk×],则会引入不可交换误差降低姿态角的解算精度,因此,按如下公式引入单子样+前一周期的等效旋转矢量补偿不可交换误差:
Figure BDA0003630294470000092
公式(14)中,Φk为等效旋转矢量,ΔΘk为角增量,ΔΘk=ΔT·(bωm,k-bg,k-1)^
S12误差状态预测和校正
名义状态r和bg的误差状态分别定义为δφ、δbg,其中,误差状态的运动方程如下:
Figure BDA0003630294470000101
Figure BDA0003630294470000102
将公式(4)按矩阵形式表示为:
Figure BDA0003630294470000103
公式(5)中,Fk-1是误差状态δx1的状态转移矩阵,按如下公式表示:
Figure BDA0003630294470000104
误差状态的先验协方差
Figure BDA0003630294470000105
的预测方程按如下公式表示:
Figure BDA0003630294470000106
公式(7)中,Q为预测方程的过程噪声矩阵;
由于陀螺仪中存在零偏,只进行预测会造成姿态漂移,因此此处引入加速度计对预测结果进行校正,观测方程按如下公式表示:
z1,k=h1(x1,k,vk) (8),
公式(8)中,z1,k=am为加速度计测量值,vk为加速度计测量噪声。由于状态变量是选用姿态旋转矩阵的最后一行,所以h1(·)是一个线性函数;
根据链式求导法则,观测方程关于δx1的雅克比矩阵如下:
Figure BDA0003630294470000111
公式(9)中,Hx1=[I3 03],
Figure BDA0003630294470000112
因此,误差状态的后验校正过程按如下公式表示:
Figure BDA0003630294470000113
Figure BDA0003630294470000114
Figure BDA0003630294470000115
公式(10)中,R1为观测方程的噪声矩阵,I6为6*6的单位矩阵。
S13误差状态的注入和复位
将经过加速度计校正后的误差状态按如下公式注入到名义状态中:
Figure BDA0003630294470000116
Figure BDA0003630294470000117
公式(11)中,Exp(δφ)是误差状态到上SO(3)的指数映射,其计算公式根据罗德里格斯公式给出,如下:
Exp(δφ)=I3+u^sinφ+(1-cosφ)(u^)2 (12),
将误差状态注入到名义状态之后,需要将误差状态清零,具体按如下同时操作:
Figure BDA0003630294470000121
Figure BDA0003630294470000122
S2第二层姿态估计算法
第二层姿态估计算法使用Kalman滤波来估计航向角
Figure BDA0003630294470000123
使用磁力计信息校正预测值。
S21航向信息解算
选用从b系到n系的姿态旋转矩阵
Figure BDA0003630294470000124
的第一列的前两个元素
Figure BDA0003630294470000125
作为表达航向角的状态变量。其中第二层的状态向量x2为:
Figure BDA0003630294470000126
运动方程按如下公式表示:
Figure BDA0003630294470000127
公式(16)中,
Figure BDA0003630294470000128
Figure BDA0003630294470000129
Figure BDA00036302944700001210
计算出来;
状态x2的先验协方差的预测方程如下:
Figure BDA00036302944700001211
S22航向信息校正
校正过程需要引入磁力计对航向信息进行校正,而磁力计数据并不能直接校正状态x2,因此需将其转换到和状态x2同一量纲,状态x2的后验校正过程按如下公式表示:
Figure BDA0003630294470000131
Figure BDA0003630294470000132
Figure BDA0003630294470000133
公式(18)中,
Figure BDA0003630294470000134
S22磁干扰的检测
模块磁场倾角是模块检测到的磁场通量线与地球表面之间的夹角,当磁干扰发生时模块磁场倾角会发生改变,但是其角度变化方向不确定,可能增大,也可能减小。磁力计模值可用来度量磁干扰的强度,其变化方向也不能完全确定,一般情况下,磁干扰越强,模值越大。
鉴于磁干扰发生的情况非常复杂,为了更精准的度量磁干扰强度,本发明使用模块磁场倾角增量和磁力计模值来评价磁干扰强度,具体评价函数按如下公式表示:
mlen=log(||bmm||sin(|dip-dipref|)) (19),
公式(19)中,dipref为无磁干扰情况下的模块磁场倾角,dip为加速度计输出矢量和磁力计输出矢量的夹角,按如下公式计算:
Figure BDA0003630294470000135
公式(20)中,am为加速度计的测量值,mm为磁力计的测量值;
为了减小磁干扰对航向信息的影响,本发明将磁干扰强度和观测矩阵R2按如下公式关联:
Figure BDA0003630294470000141
公式(21)中,σm为磁力计噪声的标准差,σd为磁干扰的基础标准差。
为了验证本发明提供的高动态和磁干扰环境下的姿态估计算法的有效性,选用XSENS公司的MTI-G-710模块对其进行验证,实验过程中记录传感器模块中加速度计、陀螺仪和磁力计的原始信号,传感器模块采样频率为100HZ,之后再将传感器原始信号导入到matlab进行分析和实验。以下通过磁干扰实验、高动态环境实验对本发明提供的姿态估计算法进行验证。
1、磁干扰实验
首先选用该实验可对本发明提到的使用模块磁场倾角增量和磁场模值相乘再取对数函数度量磁干扰强度的方法进行验证。具体实验方法如下:
将模块放置在无磁干扰的环境中并保持静止,然后将一块磁铁由远及近靠近模块再逐渐远离。分别使用模块磁场倾角、磁力计模值、本发明提到的方法去度量磁干扰程度,对比结果如图2所示。
由图2可知,在磁干扰实验中,模块磁场倾角没有及时反映磁干扰,而磁场强度出现了先下降而后上升的趋势,只有本发明提供的算法能及时并正确的反映磁干扰,因此本发明提供的方法明显由于前两者。
其次验证磁干扰对姿态角的影响,此处先将模块拿到空中做小幅度运动,然后将模块保持静止,再慢慢用磁力靠近,分别运行本文算法、EKF算法、ESKF算法,观察姿态角受影响程度,对比结果如图3所示。
由图3可知,EKF算法对磁干扰没有抵抗力,不仅航向角受到影响,俯仰、横滚角也受到了影响,而本发明提供的算法和ESKF算法则几乎不受磁干扰的影响,且由图3(3)可知,本发明提供的算法在磁铁靠近时抵抗力更优于ESKF算法。
2、高动态环境实验
将模块放置在三轴转台上,控制转台进行高动态运动,将转台输出的平台角度作为真值。使用本发明所提供的算法、EKF算法、ESKF算法分别进行姿态估计,并将解算出的姿态角和真值做差。对比结果如图4所示:
由图4可知,EKF算法在高动态环境下表现很差,在运动的初期就表现出了巨大的误差;ESKF算法和本发明在俯仰、滚转角上误差小于±2°,两者相差不大,在航向角上ESKF算法虽然在运动初期可以保持良好的精度,但随着高动态环境的持续,航向角出现了10°的漂移。
综上所述,本发明针对在高动态场景下由于非线性误差过大导致姿态估计算法容易发散和姿态角容易受到磁干扰影响的问题,提供了一种高动态和磁干扰环境下的姿态估计算法,其使用误差状态做更新和校正,利用误差状态是小量的特性,保证线性化一直有效;并使用旋转矩阵的最后一行作为状态变量,进一步降低了系统的非线性,使得本发明提供的算法可以适应高动态环境;通过分层计算,俯仰角、滚转角不使用磁力计校正,彻底规避了磁力计干扰对这两个角度的影响,同时通过综合考虑磁力计模值和模块磁场倾角增量来度量磁干扰程度,并将其和第二层姿态估计算法的观测矩阵关联,使得航向角也能在短时间的磁干扰环境中保持一定的精度。

Claims (5)

1.一种高动态和磁干扰环境下的姿态估计算法,其特征在于,包括以下步骤:
S1第一层姿态估计算法:基于线性化程度更高的误差状态解算俯仰角、滚转角,并选用方向余弦矩阵的最后一行作为俯仰角、滚转角状态变量,选用该状态变量的误差状态在EKF框架下进行更新与校正;
S2第二层姿态估计算法:选用方向余弦矩阵的第一列的前两个元素作为航向角的状态变量,并选用磁力计校正航向信息,使用磁力计模值和模块磁场倾角增量相乘再取对数函数的方法度量磁干扰强度。
2.根据权利要求1所述一种高动态和磁干扰环境下的姿态估计算法,其特征在于,所述S1具体按照以下步骤实施:
S11名义状态预测
选用从载体坐标系b系到导航坐标系n系的姿态旋转矩阵
Figure FDA0003630294460000011
的最后一行r作为表达俯仰θ、滚转角γ的状态变量,并将陀螺仪的零偏bg当成状态变量估计,第一层的名义状态向量x1按如下公式表示:
x1=[r bg] (1),
公式(1)中,r按如下公式计算:
r=[-sinθ sinγcosθ cosγcosθ] (2),
状态向量x1的运动方程按如下公式表示:
Figure FDA0003630294460000021
公式(3)中,bΦk为[tk-1,tk]时间段内的等效旋转矢量,(bΦk)^为bΦk的反对称矩阵;
S12误差状态预测和校正
误差状态的运动方程为:
Figure FDA0003630294460000022
Figure FDA0003630294460000023
公式(4)中,δφ为r的误差状态,δbg为bg的误差状态,将公式(4)按矩阵形式表示为:
Figure FDA0003630294460000024
公式(5)中,Fk-1为误差状态δx1的状态转移矩阵,其按照如下公式计算:
Figure FDA0003630294460000025
误差状态的先验协方差
Figure FDA0003630294460000026
的预测方程如下:
Figure FDA0003630294460000027
于公式(7)中,Q为预测方程的过程噪声矩阵;
引入加速度计对预测结果进行校正,观测方程按如下公式表示:
z1,k=h1(x1,k,vk) (8),
公式(8)中,z1,k=am为加速度计测量值,vk为加速度计测量噪声,由于状态变量是选用姿态旋转矩阵的最后一行,所以h1(·)是一个线性函数,
根据链式求导法则,观测方程关于δx1的雅克比矩阵如下:
Figure FDA0003630294460000031
公式(9)中,Hx1=[I3 03],
Figure FDA0003630294460000032
因此,误差状态的后验校正过程按如下公式表示:
Figure FDA0003630294460000033
Figure FDA0003630294460000034
Figure FDA0003630294460000035
公式(10)中,R1为观测方程的噪声矩阵,I6为6*6的单位矩阵;
S13误差状态的注入和复位
将经过加速度计校正后的误差状态按如下公式注入到名义状态中:
Figure FDA0003630294460000036
Figure FDA0003630294460000037
公式(11)中,Exp(δφ)是误差状态到上SO(3)的指数映射,其按照如下公式表示:
Exp(δφ)=I3+u^sinφ+(1-cosφ)(u^)2 (12),
将误差状态注入到名义状态之后,将误差状态清零,具体按如下公式操作:
Figure FDA0003630294460000041
Figure FDA0003630294460000042
3.根据权利要求2所述一种高动态和磁干扰环境下的姿态估计算法,其特征在于,所述S11中,需引入单子样+前一周期的等效旋转矢量补偿不可交换误差,具体按如下公式表示:
Figure FDA0003630294460000043
公式(14)中,Φk为等效旋转矢量,ΔΘk为角增量,ΔΘk=ΔT·(bωm,k-bg,k-1)^
4.根据权利要求3所述一种高动态和磁干扰环境下的姿态估计算法,其特征在于,所述S2具体按照以下步骤实施:
S21航向信息解算
选用从载体坐标系b系到导航坐标系n系的姿态旋转矩阵
Figure FDA0003630294460000044
的第一列的前两个元素
Figure FDA0003630294460000045
作为表达航向角
Figure FDA0003630294460000046
的状态变量,其中第二层的状态向量x2为:
Figure FDA0003630294460000051
运动方程按如下公式表示:
Figure FDA0003630294460000052
公式(16)中,
Figure FDA0003630294460000053
Figure FDA0003630294460000054
Figure FDA0003630294460000055
计算所得,
状态x2的先验协方差的预测方程按如下公式表示:
Figure FDA0003630294460000056
S22航向信息校正
由于磁力计数据并不能直接校正状态x2,因此需将其转换到和状态x2同一量纲,状态x2的后验校正过程按如下公式表示:
Figure FDA0003630294460000057
Figure FDA0003630294460000058
Figure FDA0003630294460000059
公式(18)中,
Figure FDA00036302944600000510
S23磁干扰的检测
使用模块磁场倾角增量和磁力计模值来评价磁干扰强度,具体评价函数按如下公式表示:
mlen=log(||bmm||sin(|dip-dipref|)) (19),
公式(19)中,dipref为无磁干扰情况下的模块磁场倾角,dip为加速度计输出矢量和磁力计输出矢量的夹角,dip按如下公式进行计算:
Figure FDA0003630294460000061
公式(20)中,am为加速度计的测量值,mm为磁力计的测量值。
5.根据权利要求4所述一种高动态和磁干扰环境下的姿态估计算法,其特征在于,所述S2中,为了减小磁干扰对航向信息的影响,按如下公式将磁干扰强度和观测矩阵R2关联:
Figure FDA0003630294460000062
公式(21)中,σm为磁力计噪声的标准差,σd为磁干扰的基础标准差。
CN202210486754.7A 2022-05-06 2022-05-06 一种高动态和磁干扰环境下的姿态估计算法 Active CN114877858B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210486754.7A CN114877858B (zh) 2022-05-06 2022-05-06 一种高动态和磁干扰环境下的姿态估计算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210486754.7A CN114877858B (zh) 2022-05-06 2022-05-06 一种高动态和磁干扰环境下的姿态估计算法

Publications (2)

Publication Number Publication Date
CN114877858A true CN114877858A (zh) 2022-08-09
CN114877858B CN114877858B (zh) 2023-04-14

Family

ID=82673636

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210486754.7A Active CN114877858B (zh) 2022-05-06 2022-05-06 一种高动态和磁干扰环境下的姿态估计算法

Country Status (1)

Country Link
CN (1) CN114877858B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003065793A (ja) * 2001-08-22 2003-03-05 Japan Aviation Electronics Industry Ltd 慣性装置
CN102607561A (zh) * 2012-02-28 2012-07-25 西安费斯达自动化工程有限公司 基于加速度计的飞行器欧拉角修正模型
CN102830622A (zh) * 2012-09-05 2012-12-19 北京理工大学 一种四旋翼飞行器自抗扰自动飞行控制方法
CN103822633A (zh) * 2014-02-11 2014-05-28 哈尔滨工程大学 一种基于二阶量测更新的低成本姿态估计方法
CN103954303A (zh) * 2014-05-15 2014-07-30 东南大学 一种用于磁力计导航系统航向角动态计算及校正方法
CN106643737A (zh) * 2017-02-07 2017-05-10 大连大学 风力干扰环境下四旋翼飞行器姿态解算方法
US20190169979A1 (en) * 2017-12-04 2019-06-06 Hrl Laboratories, Llc Continuous Trajectory Calculation for Directional Drilling
CN110672078A (zh) * 2019-10-12 2020-01-10 南京理工大学 一种基于地磁信息的高旋弹丸姿态估计方法
CN111982100A (zh) * 2020-07-07 2020-11-24 广东工业大学 一种无人机的航向角解算算法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2003065793A (ja) * 2001-08-22 2003-03-05 Japan Aviation Electronics Industry Ltd 慣性装置
CN102607561A (zh) * 2012-02-28 2012-07-25 西安费斯达自动化工程有限公司 基于加速度计的飞行器欧拉角修正模型
CN102830622A (zh) * 2012-09-05 2012-12-19 北京理工大学 一种四旋翼飞行器自抗扰自动飞行控制方法
CN103822633A (zh) * 2014-02-11 2014-05-28 哈尔滨工程大学 一种基于二阶量测更新的低成本姿态估计方法
CN103954303A (zh) * 2014-05-15 2014-07-30 东南大学 一种用于磁力计导航系统航向角动态计算及校正方法
CN106643737A (zh) * 2017-02-07 2017-05-10 大连大学 风力干扰环境下四旋翼飞行器姿态解算方法
US20190169979A1 (en) * 2017-12-04 2019-06-06 Hrl Laboratories, Llc Continuous Trajectory Calculation for Directional Drilling
CN110672078A (zh) * 2019-10-12 2020-01-10 南京理工大学 一种基于地磁信息的高旋弹丸姿态估计方法
CN111982100A (zh) * 2020-07-07 2020-11-24 广东工业大学 一种无人机的航向角解算算法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
陈雨荻: "多旋翼飞行器的姿态与导航信息融合算法研究", 《中国优秀硕士学位论文全文数据库(电子期刊)》 *
韩冬: "随钻姿态测量重力加速度自适应提取算法", 《仪器仪表学报》 *

Also Published As

Publication number Publication date
CN114877858B (zh) 2023-04-14

Similar Documents

Publication Publication Date Title
CN110398257B (zh) Gps辅助的sins系统快速动基座初始对准方法
CN103630137B (zh) 一种用于导航系统的姿态及航向角的校正方法
Pylvänäinen Automatic and adaptive calibration of 3D field sensors
EP3933166B1 (en) Attitude measurement method
CN111238535B (zh) 一种基于因子图的imu误差在线标定方法
US7620514B2 (en) Method and arrangement for correcting an angle-measuring and/or distance-measuring sensor system
US20050065728A1 (en) Method and apparatus for compensating attitude of inertial navigation system and method and apparatus for calculating position of inertial navigation system using the same
EP2135033A2 (en) Auto-calibration of orientation sensing system
CN112798021B (zh) 基于激光多普勒测速仪的惯导系统行进间初始对准方法
CN107270891B (zh) 基于抗差估计的惯性地磁匹配定位方法
Wang et al. A new magnetic compass calibration algorithm using neural networks
CN109059907A (zh) 轨迹数据处理方法、装置、计算机设备和存储介质
CN116817896B (zh) 一种基于扩展卡尔曼滤波的姿态解算方法
CN115046539A (zh) Mems电子罗盘动态校准方法
CN105937911A (zh) 一种磁传感器姿态解算方法
CN106200377A (zh) 一种飞行器控制参数的估计方法
US9863867B2 (en) Automatically updating hard iron and soft iron coefficients of a magnetic sensor
CN111707292B (zh) 一种自适应滤波的快速传递对准方法
CN110672078A (zh) 一种基于地磁信息的高旋弹丸姿态估计方法
CN114877858B (zh) 一种高动态和磁干扰环境下的姿态估计算法
CN107036576B (zh) 基于差商法磁测旋转飞行器滚转角的实时解算方法
CN115839726A (zh) 磁传感器和角速度传感器联合标定的方法、系统及介质
Hemanth et al. Calibration of 3-axis magnetometers
Łuczak Novel algorithm for tilt measurements using MEMS accelerometers
CN109470267A (zh) 一种卫星姿态滤波方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant