CN103852085A - 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法 - Google Patents

一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法 Download PDF

Info

Publication number
CN103852085A
CN103852085A CN201410116682.2A CN201410116682A CN103852085A CN 103852085 A CN103852085 A CN 103852085A CN 201410116682 A CN201410116682 A CN 201410116682A CN 103852085 A CN103852085 A CN 103852085A
Authority
CN
China
Prior art keywords
error
gma
omega
delta
inertial navigation
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
CN201410116682.2A
Other languages
English (en)
Other versions
CN103852085B (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN201410116682.2A priority Critical patent/CN103852085B/zh
Publication of CN103852085A publication Critical patent/CN103852085A/zh
Application granted granted Critical
Publication of CN103852085B publication Critical patent/CN103852085B/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
    • G01C25/00Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C25/00Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass
    • G01C25/005Manufacturing, calibrating, cleaning, or repairing instruments or devices referred to in the other groups of this subclass initial alignment, calibration or starting-up of inertial devices

Abstract

本发明公开了一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,属于惯性技术领域。本发明采用9次翻转路径设计,对系统输出惯性数据应用最小二乘拟合的方法求解得到光纤捷联惯导系统18项误差系数;再通过一次翻转得到零偏误差,进而得到全部21项器件误差参数;本发明提供的技术方案利用六面体或其它相似的可翻转装置即可完成现场标定试验,克服了传统实验室标定的不足,提高了系统实际使用精度。

Description

一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法
技术领域
本发明属于惯性技术领域,涉及一种光纤捷联惯导系统现场标定方法,具体地说,是指一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法。
背景技术
光纤陀螺具有精度高、启动快、动态范围大、抗振动冲击能力强及成本低等优点,是惯性仪表领域的发展趋势。近年来,光纤陀螺技术的迅猛发展推动了光纤捷联惯导系统在陆、海、空、天领域的应用。光纤捷联惯导在使用前必须通过实验室转台标定试验确定出其核心部件即光纤陀螺和加速度计的标度因数和各项误差系数,在后续的导航计算中进行补偿。
但是,通过实验室转台试验标定出的光纤惯导系统各项误差系数并不是固定不变的,包括光纤陀螺零位漂移误差、标度因数、安装误差和加速度计常值误差、标度困数、安装误差等。这些误差参数随着系统使用或存放时间的推移而变化,尤其是光纤陀螺零位漂移和加速度计偏值,每次上电启动都不相同,且时间间隔越长变化越大。光纤陀螺和加速度计误差参数的改变直接导致光纤惯导系统的精度降低,使惯导系统无法使用要求。
因此,通常需要对光纤惯导系统进行半年或三个月一次的定期标定,而且传统的基于精密转台的标定方法复杂、耗时长,这为使用单位增添了巨大的工作量及人力、物力、财力的消耗。因此,在使用现场对光纤惯导的各项误差系数进行标定,不仅可以减少甚至取消定期标定,还可以提高光纤捷联惯导的实际使用精度。但是,在现场没有精密的转台作为测试基准,不能对充纤陀螺捷联惯导进行精确定向,所以传统的基于精密转台的静态多位置标定方法和速率标定方法都无法实施,为光纤捷联惯导的现场标定带来了极大难度。
在没有精确的基准设备和复杂的转位机构的情况下,现场标定就不能采用传统的速率位置法,必须依靠系统级方法来克服基准信息缺乏的困难。系统级方法以光纤惯导的导航误差为观测量,通过惯导系统误差参数与导航误差之间的关系,建立惯导系统误差参数与导航误差之间的方程式,进而求解出惯导系统的误差参数。因此,系统级标定方法可以克服使用现场设备条件不足的缺点,并获得高精度的系统误差标定接过。
参考文献[1](光电工程,刘百奇,房建成.光纤陀螺IMU的六位置旋转现场标定新方法[J].35(1),2008:60-65)公开了一种光纤陀螺IMU的六位置旋转现场标定新方法,本文采用光纤陀螺IMU在六个位置上进行12次旋转,然后根据光纤陀螺IMU的误差模型建立42个非线性输入输出方程,通过旋转积分和对称位置误差相消,消除方程中的非线性项,最终求出陀螺标度因数、陀螺常值漂移、陀螺安装误差和加速度计常值偏置等15个误差系数。但是该方法不能够标定出加速度计通道的标度因数和安装误差。
参考文献[2](测控技术,颜开思,李岁劳,龚柏春,贾继超.[J].30(5),2011:106-109)公开了一种基于平台和正六面体的惯导系统现场标定技术,该文献通过翻转正六面体使对称位置误差相消,并且在对准中获取姿态信息,同时精确标定出陀螺漂移和加速度计零偏。最后对理论分析结果进行了仿真验证,仿真结果表明该方案可以实现外场条件下的陀螺漂移和加速度计零偏的精确标定。但是该方法不能够标定出陀螺、加速度的标度因数和安装误差。
参考文献[3](吴赛成,秦石乔,王省书,胡春生.[J].中国惯性技术学报,19(2),2011:185-189)公开了一种激光陀螺惯性测量单元系统级标定方法,该文献建立了附加约束条件的陀螺和加速度计安装坐标系数学模型,根据陀螺和加速度计的输出误差方程,从惯性导航基本误差方程出发推导了惯性测量单元的系统级误差参数标定Kalman滤波模型,该模型包含了陀螺和加速度计零偏、比例因子、安装误差在内共21维标定误差状态变量,且仅以速度解算误差为观测量。但是该方法标定步骤较多,标定时间过长,缺少实验数据验证。
参考文献[4](专利文献号CN102607594A)公开了一种捷联惯导光纤陀螺系统误差参数现场标定方法,所述现场标定方法通过姿态测量仪器给出载体姿态角,选取姿态作为观测量,标定出光纤捷联惯导系统光纤陀螺各项误差系数。但是该方法需要现场提供姿态测量辅助器件,实时精确测量载体姿态角,并且要与光纤陀螺输出值保持一致。
发明内容
本发明的目的在于提供一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,减少或者甚至取消光纤捷联惯导系统周期性地返实验室校标,提高系统实际使用精度。
本发明采用9次翻转路径设计,包括转动轴、转动顺序和转动角度等;应用最小二乘拟合的方法得到光纤捷联惯导系统全部21项器件误差参数;利用六面体或其它相似的可翻转装置即可完成现场标定试验。具体方法步骤如下:
第一步:将光纤捷联惯导系统通过工装安装在正六面体装置上,锁紧。连接系统、电源和采集计算机之间的线缆,并检查正确。
第二步:将正六面体装置置于水平面上,上电预热使光纤捷联惯导系统达到热平衡状态,并装订光纤捷联惯导系统的初始位置参数,包括初始的经度、纬度和高度。
第三步:采用“静止-转动-静止”进行手动翻转正六面体装置,按照转动路径序列完成9次翻转。转动前后每个位置静止3~5min,并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据。
第四步:对惯性器件误差进行建模。
第五步:对9组惯性器件数据进行处理,采用最小二乘拟合的方法求解,得到除陀螺零偏以外的其他18项误差系数;
第六步:根据标定得到18项误差系数,对保存的9组惯性器件数据进行补偿,得到新的9组惯性器件数据。
第七步:对于第六步中得到的新的9组惯性器件数据,重复第五步和第六步3~5次,即迭代器件18项误差系数值进行重复最小二乘拟合计算,直至收敛。累加每次迭代估计得到器件18项误差系数值,即为最终标定参数值。
第八步:根据上述标定的18项误差系数,对光纤捷联惯导系统进行补偿。
第九步:重复第一步到第二步,按照序列1进行手动手动翻转六面体依次,转动前后分别静止时间30min,保存陀螺仪输出数据,基于该输出数据,计算陀螺零偏误差。
第十步:根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿,采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可,完成了光纤捷联惯导系统21项误差参数的现场标定。
本发明的有益效果在于:
本发明所提出的方法可以在现场完成光纤捷联惯导系统21项误差参数的标定,克服了传统实验室标定的不足,提高了系统实际使用精度。
附图说明
图1为本发明提供的基于最小二乘拟合的光纤捷联惯导系统现场标定方法流程图;
图2A和图2B分别为本发明实施例中静态和摇摆情况下现场标定补偿前后20min导航定位误差对比曲线。
具体实施方式
下面结合附图和实施例对本发明进行详细说明。
本发明提供一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,如图1所示,所述方法包括如下步骤:
第一步:将光纤捷联惯导系统通过工装安装在正六面体装置上,锁紧。连接系统、电源和采集计算机之间的线缆,并检查正确。
第二步:将正六面体装置置于水平面上,上电预热使光纤捷联惯导系统达到热平衡状态,并装订光纤捷联惯导系统的初始位置参数,包括初始的经度、纬度和高度。
第三步:按照表1转动路径序列,表1中转动轴X、Y、Z,光纤捷联惯导系统初始姿态为0时,光纤捷联惯导系统XYZ轴与导航坐标系东北天位置重合。
采用“静止-转动-静止”进行手动翻转正六面体装置,完成9次翻转,转动角允许存在±10°误差。转动前后每个位置静止3~5min,并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据。
表1最小二乘拟合法转动路径序列(转动角单位:度)
Figure BDA0000482609920000041
第四步:对惯性器件误差进行建模,包括光纤陀螺误差模型和加速度计误差模型,分别如下:
δω ib b = δω ibx b δω iby b δω ibz b = gB x gB y gB z + gSF x gMA xy gMA xz gMA yx gSF y gMA yz gMA zx gMA zy gSF z ω ibx b ω iby b ω ibz b
δf ib b = δf ibx b δf iby b δf ibz b = aB x aB y a B z + aSF x 0 0 aMA yx aSF y 0 aMA zx aMA zy aSF z f ibx b f iby b f ibz b
式中δωib b表示陀螺仪的误差输出矢量;δωibx b、δωiby b、δωibz b分别表示由三轴陀螺误差引起的误差角速度。ωibx b、ωiby b、ωibz b分别表示三轴陀螺测量值;gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间安装误差角;gBx、gBy、gBz分别表示三轴陀螺仪零偏误差;
式中
Figure BDA0000482609920000044
表示加速度计的误差输出;δfibx b、δfiby b、δfibz b分别表示由三轴加速度计误差引起的误差加速度;fibx b、fiby b、fibz b分别表示三轴加速度计测量值;aSFx、aSFy、aSFz分别为三轴加速度计的标度因数误差;aBx、aBy、aBz分别为三轴加速度计的零偏;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角;
第五步:对保存的9组惯性器件数据进行处理,采用最小二乘拟合的方法求解,得到除陀螺零偏以外的其他18项误差系数。
所述最小二乘拟合法包含以下几个步骤:
步骤1:建立器件误差与系统速度误差一阶导数变化量的数学模型。
K·X=A
其中K表示转动系数矩阵,X表示误差向量,A表示观测矩阵。误差向量X=[Xa;Xg];
Xa=[aBx aBy aBz aSFx aMAyx aSFy aMAzx aMAzy aSFz]T
Xg=[gSFx gMAxy gMAxz gMAyx gSFy gMAyz gMAzx gMAzy gSFz]T
式中gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间的安装误差角;aSFx、aSFy、aSFz分别为三轴加速度计标度因数误差;aBx、aBy、aBz分别为三轴加速度计零偏;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角;
步骤2:求解观测矩阵A。
观测矩阵
Figure BDA0000482609920000051
表示T2时刻速度误差一阶导数矢量,
Figure BDA0000482609920000052
表示T1时刻速度误差一阶导数矢量。
观测矩阵由三向速度误差一阶导数组成,光纤捷联惯导系统静止时,速度误差即为光纤捷联惯导系统导航输出的速度值。采用标准卡尔曼滤波计算速度误差一阶导数值。
建立状态方程: d dt δV i ( t ) δ V · i ( t ) δ V · · i ( t ) = 0 1 0 0 0 1 0 0 0 δV i ( t ) δ V · i ( t ) δ V · · i ( t ) + 0 0 w i
建立量测方程: Z i = 1 0 0 δV i ( t ) δ V · i ( t ) δ V · · i ( t ) + μ i
其中,Zi表示观测矢量;δVi(t)=Vi(t),i=x,y,z,δVi(t)表示速度误差矢量,Vi(t)表示导航解算速度矢量,wi和μi分别表示状态噪声和观测噪声。采用标准卡尔曼滤波方程进行迭代,滤波完毕即可得到速度误差一阶导数的估计值
Figure BDA0000482609920000055
求取转动前后两次速度误差一阶导数并作差,即可得到系统观测阵A。
步骤3:求解转动系数矩阵K。
转动系数矩阵 K = df x - g · dφ y df y g · dφ x df z 0 1 × 9
式中,df=[dfx dfy dfz]T表示加速度计误差系数阵,dφ=[dφx dφy dφz]T表示陀螺误差系数阵。
d f = C b n ( T 2 ) M 2 - C b n ( T 1 ) M 1 , 其中,
M i = 1 0 0 f ibx b ( T i ) 0 0 0 0 0 0 1 0 0 f ibx b ( T i ) f iby b ( T i ) 0 0 0 0 0 1 0 0 0 f ibx b ( T i ) f iby b ( T i ) f ibz b ( T i ) , 其中i=1,2;
d φ = - ∫ T 1 T 2 C b ( T 1 ) n · C b ( t ) b ( T 1 ) · ω x ( t ) 0 0 0 ω y ( t ) 0 0 0 ω z ( t ) M g · dt
其中,fibx b(Ti)、fiby b(Ti)、fibz b(Ti)分别表示Ti时刻三轴加速度计测量值;
Figure BDA00004826099200000510
为惯组初始时刻T1的姿态矩阵,ωx(t)、ωy(t)、ωz(t)为t时刻的转动角速度,而为t时刻b系到T1时刻b系的转换矩阵,Mg为系数矩阵,这两者与转动轴向有关;所述b系是指载体坐标系,n代表n系,是指导航坐标系,本发明中定义导航坐标系为东北天,即当地地理坐标系。
若光纤捷联惯导系统绕X轴转动,则:
C b ( t ) b ( T 1 ) = 1 0 0 0 cos ( ω x ( t ) ) sin ( ω x ( t ) ) 0 - sin ( ω x ( t ) ) cos ( ω x ( t ) ) , M g = 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0
若光纤捷联惯导系统绕Y轴转动,则:
C b ( t ) b ( T 1 ) | = cos ( ω y ( t ) ) 0 - sin ( ω y ( t ) ) 0 1 0 sin ( ω y ( t ) ) 0 cos ( ω y ( t ) ) , M g = 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
若光纤捷联惯导系统绕Z轴转动,则:
C b ( t ) b ( T 1 ) | = cos ( ω z ( t ) ) sin ( ω z ( t ) ) 0 - sin ( ω z ( t ) ) cos ( ω z ( t ) ) 0 0 0 1 , M g = 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1
步骤4:采用最小二乘方法求解误差向量:X=KT·(KKT)-1·A。由于转动系数矩阵K的秩rank(K)=15,得到的误差系数有15项,包括三轴陀螺仪标度因数误差gSFx、gSFy、gSFz,各轴陀螺仪间安装误差角gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy,三轴加速度计零偏aBx、aBy、aBz;各轴加速度计间安装误差角aMAyx、aMAzx、aMAzy
步骤5:求解三轴加速度计标度因数;
采用第1次转动数据: aSF z = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 1
采用第4次转动数据: aSF y = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 4
采用第7次转动数据: aSF x = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 7
式中,
Figure BDA00004826099200000610
表示转动前从0时刻到T1时刻的时间内天向速度误差一阶导数之和,
Figure BDA00004826099200000611
表示转动后从0时刻到T2时刻的时间内天向速度误差一阶导数之和。
第六步:根据上述步骤4和步骤5中标定得到18项误差系数,对保存的9组惯性器件数据进行补偿,得到新的9组惯性器件数据。其中补偿公式如下:
w x w y w z = ( I + gSF x gMA xy gMA xz gMA yx gSF y gMA yz gMA zx gMA zy gSF z ) ω ibx b ω iby b ω ibz b
f x f y f z = ( I + aSF x 0 0 aMA yx aSF y 0 aMA zx aMA zy aSF z ) f ibx b - aB x f iby b - aB y f ibz b - aB z
式中gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差标定结果;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间安装误差角标定结果;aSFx、aSFy、aSFz分别为三轴加速度计标度因数误差标定结果;aBx、aBy、aBz分别为三轴加速度计零偏标定结果;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角标定结果;ωibx b、ωiby b、ωibz b分别表示三轴陀螺原测量值;fibx b、fiby b、fibz b分别表示三轴加速度计原测量值;I表示三阶单位阵;wx、wy、wz分别表示三轴陀螺补偿后测量值;fx、fy、fz分别表示三轴加速度计补偿后测量值。
第七步:重复第五步和第六步3~5次,即迭代器件18项误差系数值进行重复最小二乘拟合计算,直至收敛。累加每次迭代估计得到器件18项误差系数值,即为最终标定参数值。
第八步:根据上述最终标定的18项误差系数值,按照第六步中补偿公式对光纤捷联惯导系统进行补偿。
第九步:重复第一步到第二步,按照表1的序列1进行手动翻转六面体,转动前后分别静止时间30min,保存陀螺仪输出数据。基于输出数据,计算陀螺零偏误差gB=[gBx gBy gBz]T。计算公式如下:
gB = 1 2 { ( ω ib b ( T 1 ) + ω ib b ( T 2 ) ) - [ C b n ( T 1 ) + C b n ( T 1 ) ] T ( I - φ × ) ω ie n }
其中,φ=[0 0 φz]T,φz表示航向角误差, φ z = Δω ibx b - ΔC b n ( 1,2 ) ω iey n - ΔC b n ( 1,3 ) ω iez n ΔC b n ( 1,1 ) ω iey n ,
Δω ib b = ω ib b ( T 2 ) - ω ib b ( T 1 ) , ΔC b n = C b n ( T 2 ) - C b n ( T 1 )
式中
Figure BDA0000482609920000074
分别表示T1、T2时刻三轴陀螺测量矢量,表示转动前后两时刻陀螺测量误差矢量、
Figure BDA0000482609920000076
表示
Figure BDA0000482609920000077
第一行分量,
Figure BDA0000482609920000078
分别表示T1、T2时刻载体系b系到导航系n系的方向余弦矩阵,表示转动前后两时刻方向余弦矩阵误差,
Figure BDA00004826099200000710
表示第i行第j列元素值,
Figure BDA00004826099200000712
表示导航系下地球速度矢量,
Figure BDA00004826099200000713
分别表示导航系下地球速度矢量的北向分量和天向分量。
第十步:根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿,采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可,完成了光纤捷联惯导系统21项误差参数的现场标定。
实施例
第一步:将光纤捷联惯导系统通过工装安装在正六面体装置上,锁紧。连接系统、电源、采集计算机之间的线缆,并检查正确。
第二步:将正六面体装置置于水平面上,上电预热使光纤捷联惯导系统达到热平衡状态,并装订光纤捷联惯导系统的初始位置参数,包括初始的经度、纬度和高度。
第三步:按照表1转动路径,采用“静止-转动-静止”进行手动翻转正六面体装置,完成9次翻转,转动角允许存在±10°误差。转动前后每个位置静止3~5min,并保存9次转动数据。
第四步:对保存的9组数据进行处理,采用最小二乘拟合的方法求解,并迭代5次,得到除陀螺零偏以外的其他18项误差参数。
第五步:先补偿18项误差参数,重复第一步到第二步,按照表1的序列1进行手动翻转六面体,转动前后分别静止时间30min,保存陀螺输出数据。
第六步:基于第五步数据,计算三轴陀螺零偏误差值,并对光纤捷联惯导系统陀螺输出进行补偿,完成光纤捷联惯导系统21项误差参数的现场标定。
第七步:将系统断电,1天后重新启动光纤捷联惯导系统,首先静态采集23min惯性器件数据,然后再将系统断电3-5h,将系统安装于双轴转台上,先静态3min再摇摆20min,一共采集23min惯性器件数据。
第八步:离线处理两组惯性器件数据,将第五步得到的现场标定结果对数据进行补偿,采用解析式粗对准3min和纯惯性导航,对比补偿前后两组数据的纯惯性导航结果。
结果分析:
(1)采用最小二乘拟合法迭代5次得到的18项误差参数值,和最后补偿得到的三轴陀螺零偏值如表2所示。从表2中可以看出,经过5次迭代,参数估计值渐近收敛,其中三轴加速度计零偏误差值收敛至60μg以内,加速度计标度因数收敛值在5ppm以内,陀螺标度因数收敛值在35ppm以内,各轴间安装误差收敛至14″以内。
表2最小二乘拟合法现场标定实验结果
Figure BDA0000482609920000081
Figure BDA0000482609920000091
(2)对比现场标定补偿前后的数据导航结果如图2A和图2B所示。图2A是20min静态水平定位误差对比曲线,图2B是摇摆情况下水平定位误差对比曲线。从图2A和图2B中可以看出,不管是静态还是动态情况下,系统的导航定位误差减小了1倍以上,因此采用本发明提供的现场标定方法补偿后的数据精度更高。

Claims (6)

1.一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于通过如下步骤实现:
第一步:将光纤捷联惯导系统通过工装安装在正六面体装置上,锁紧;连接系统、电源和采集计算机之间的线缆,并检查正确;
第二步:将正六面体装置置于水平面上,上电预热使光纤捷联惯导系统达到热平衡状态,并装订光纤捷联惯导系统的初始位置参数,包括初始的经度、纬度和高度;
第三步:采用“静止-转动-静止”进行手动翻转正六面体装置,按照转动路径序列完成9次翻转;转动前后每个位置静止3~5min,并保存9次转动过程中光纤捷联惯导系统输出所有惯性器件数据;所述转动路径序列如下:
Figure FDA0000482609910000011
第四步:对惯性器件误差进行建模;
第五步:对保存的9组惯性器件数据进行处理,采用最小二乘拟合的方法求解,得到除陀螺零偏以外的其他18项误差系数;
第六步:根据标定得到18项误差系数,对保存的9组惯性器件数据进行补偿,得到新的9组惯性器件数据;
第七步:对于第六步中得到的新的9组惯性器件数据,重复第五步和第六步3~5次,即迭代器件18项误差系数值进行重复最小二乘拟合计算,直至收敛;累加每次迭代估计得到器件18项误差系数值,即为最终标定参数值;
第八步:根据上述标定的18项误差系数,对光纤捷联惯导系统进行补偿;
第九步:重复第一步到第二步,按照第三步中的序列1进行手动翻转正六面体装置,转动前后分别静止时间30min,保存陀螺仪输出数据,基于该输出数据,计算三轴陀螺零偏误差;
第十步:根据求解得到的三轴陀螺零偏误差对光纤捷联惯导系统误差再次补偿,采用光纤捷联惯导系统输出陀螺值减去标定得到的陀螺零偏值即可,完成了光纤捷联惯导系统21项误差参数的现场标定。
2.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于:第三步中所述转动角允许存在±10°误差。
3.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于:第四步中所述的对惯性器件误差进行建模,包括光纤陀螺误差模型和加速度计误差模型,分别如下:
δω ib b = δω ibx b δω iby b δω ibz b = gB x gB y gB z + gSF x gMA xy gMA xz gMA yx gSF y gMA yz gMA zx gMA zy gSF z ω ibx b ω iby b ω ibz b
δf ib b = δf ibx b δf iby b δf ibz b = aB x aB y a B z + aSF x 0 0 aMA yx aSF y 0 aMA zx aMA zy aSF z f ibx b f iby b f ibz b
式中δωib b表示陀螺仪的误差输出矢量;δωibx b、δωiby b、δωibz b分别表示由三轴陀螺误差引起的误差角速度;ωibx b、ωiby b、ωibz b分别表示三轴陀螺测量值;gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间安装误差角;gBx、gBy、gBz分别表示三轴陀螺仪零偏误差;
式中
Figure FDA0000482609910000025
表示加速度计的误差输出;δfibx b、δfiby b、δfibz b分别表示由三轴加速度计误差引起的误差加速度;fibx b、fiby b、fibz b分别表示三轴加速度计测量值;aSFx、aSFy、aSFz分别为三轴加速度计的标度因数误差;aBx、aBy、aBz分别为三轴加速度计的零偏;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角。
4.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于:所述最小二乘拟合法包含以下几个步骤:
步骤1:建立器件误差与系统速度误差一阶导数变化量的数学模型:
K·X=A
其中K表示转动系数矩阵,X表示误差向量,A表示观测矩阵;误差向量X=[Xa;Xg];
Xa=[aBx aBy aBz aSFx aMAyx aSFy aMAzx aMAzy aSFz]T
Xg=[gSFx gMAxy gMAxz gMAyx gSFy gMAyz gMAzx gMAzy gSFz]T
式中gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间的安装误差角;aSFx、aSFy、aSFz分别为三轴加速度计标度因数误差;aBx、aBy、aBz分别为三轴加速度计零偏;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角;
步骤2:求解观测矩阵A;
观测矩阵
Figure FDA0000482609910000023
表示T2时刻速度误差一阶导数矢量,
Figure FDA0000482609910000024
表示T1时刻速度误差一阶导数矢量;
观测矩阵由三向速度误差一阶导数组成,光纤捷联惯导系统静止时,速度误差即为光纤捷联惯导系统导航输出的速度值;采用标准卡尔曼滤波计算速度误差一阶导数值;
建立状态方程: d dt δV i ( t ) δ V · i ( t ) δ V · · i ( t ) = 0 1 0 0 0 1 0 0 0 δV i ( t ) δ V · i ( t ) δ V · · i ( t ) + 0 0 w i
建立量测方程: Z i = 1 0 0 δV i ( t ) δ V · i ( t ) δ V · · i ( t ) + μ i
其中,Zi表示观测矢量;δVi(t)=Vi(t),i=x,y,z,δVi(t)表示速度误差矢量,Vi(t)表示导航解算速度矢量,wi和μi分别表示状态噪声和观测噪声;采用标准卡尔曼滤波方程进行迭代,滤波完毕即得到速度误差一阶导数的估计值
Figure FDA0000482609910000033
求取转动前后两次速度误差一阶导数并作差,即得到系统观测阵A;
步骤3:求解转动系数矩阵K;
转动系数矩阵 K = df x - g · dφ y df y g · dφ x df z 0 1 × 9
式中,df=[dfx dfy dfz]T表示加速度计误差系数阵,dφ=[dφx dφy dφz]T表示陀螺误差系数阵;
d f = C b n ( T 2 ) M 2 - C b n ( T 1 ) M 1 , 其中,
M i = 1 0 0 f ibx b ( T i ) 0 0 0 0 0 0 1 0 0 f ibx b ( T i ) f iby b ( T i ) 0 0 0 0 0 1 0 0 0 f ibx b ( T i ) f iby b ( T i ) f ibz b ( T i ) , 其中i=1,2;
d φ = - ∫ T 1 T 2 C b ( T 1 ) n · C b ( t ) b ( T 1 ) · ω x ( t ) 0 0 0 ω y ( t ) 0 0 0 ω z ( t ) M g · dt
其中,fibx b(Ti)、fiby b(Ti)、fibz b(Ti)分别表示Ti时刻三轴加速度计测量值;
Figure FDA0000482609910000038
为惯组初始时刻T1的姿态矩阵,ωx(t)、ωy(t)、ωz(t)为t时刻的转动角速度,而
Figure FDA0000482609910000039
为t时刻b系到T1时刻b系的转换矩阵,Mg为系数矩阵,这两者与转动轴向有关;
若光纤捷联惯导系统绕X轴转动,则:
C b ( t ) b ( T 1 ) = 1 0 0 0 cos ( ω x ( t ) ) sin ( ω x ( t ) ) 0 - sin ( ω x ( t ) ) cos ( ω x ( t ) ) , M g = 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0
若光纤捷联惯导系统绕Y轴转动,则:
C b ( t ) b ( T 1 ) | = cos ( ω y ( t ) ) 0 - sin ( ω y ( t ) ) 0 1 0 sin ( ω y ( t ) ) 0 cos ( ω y ( t ) ) , M g = 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0
若光纤捷联惯导系统绕Z轴转动,则:
C b ( t ) b ( T 1 ) | = cos ( ω z ( t ) ) sin ( ω z ( t ) ) 0 - sin ( ω z ( t ) ) cos ( ω z ( t ) ) 0 0 0 1 , M g = 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1
步骤4:采用最小二乘方法求解误差向量:X=KT·(KKT)-1·A;由于转动系数矩阵K的秩rank(K)=15,得到的误差系数有15项,包括三轴陀螺仪标度因数误差gSFx、gSFy、gSFz,各轴陀螺仪间安装误差角gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy,三轴加速度计零偏aBx、aBy、aBz;各轴加速度计间安装误差角aMAyx、aMAzx、aMAzy
步骤5:求解三轴加速度计标度因数;
采用第1次转动数据: aSF z = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 1
采用第4次转动数据: aSF y = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 4
采用第7次转动数据: aSF x = 1 2 ( T 1 + T 2 ) g ( Σ t = 0 t = T 2 δ V · z ( t ) + Σ t = 0 t = T 1 δ V · z ( t ) ) 7
式中,
Figure FDA0000482609910000046
表示转动前从0时刻到T1时刻的时间内天向速度误差一阶导数之和,
Figure FDA0000482609910000047
表示转动后从0时刻到T2时刻的时间内天向速度误差一阶导数之和。
5.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于:第六步中所述的对保存的9组惯性器件数据进行补偿,得到新的9组惯性器件数据,其中补偿公式如下:
w x w y w z = ( I + gSF x gMA xy gMA xz gMA yx gSF y gMA yz gMA zx gMA zy gSF z ) ω ibx b ω iby b ω ibz b
f x f y f z = ( I + aSF x 0 0 aMA yx aSF y 0 aMA zx aMA zy aSF z ) f ibx b - aB x f iby b - aB y f ibz b - aB z
式中gSFx、gSFy、gSFz分别表示三轴陀螺仪标度因数误差标定结果;gMAxy、gMAxz、gMAyx、gMAyz、gMAzx、gMAzy分别表示各轴陀螺仪间安装误差角标定结果;aSFx、aSFy、aSFz分别为三轴加速度计标度因数误差标定结果;aBx、aBy、aBz分别为三轴加速度计零偏标定结果;aMAyx、aMAzx、aMAzy分别表示各轴加速度计间安装误差角标定结果;ωibx b、ωiby b、ωibz b分别表示三轴陀螺原测量值;fibx b、fiby b、fibz b分别表示三轴加速度计原测量值;I表示三阶单位阵;wx、wy、wz分别表示三轴陀螺补偿后测量值;fx、fy、fz分别表示三轴加速度计补偿后测量值。
6.根据权利要求1所述的一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法,其特征在于:计算陀螺零偏误差计算公式如下:
gB = 1 2 { ( ω ib b ( T 1 ) + ω ib b ( T 2 ) ) - [ C b n ( T 1 ) + C b n ( T 1 ) ] T ( I - φ × ) ω ie n }
其中,φ=[0 0 φz]T,φz表示航向角误差, φ z = Δω ibx b - ΔC b n ( 1,2 ) ω iey n - ΔC b n ( 1,3 ) ω iez n ΔC b n ( 1,1 ) ω iey n , Δω ib b = ω ib b ( T 2 ) - ω ib b ( T 1 ) , ΔC b n = C b n ( T 2 ) - C b n ( T 1 ) ;
式中分别表示T1、T2时刻三轴陀螺测量矢量,表示转动前后两时刻陀螺测量误差矢量、表示第一行分量,
Figure FDA0000482609910000058
分别表示T1、T2时刻载体系b系到导航系n系的方向余弦矩阵,
Figure FDA0000482609910000059
表示转动前后两时刻方向余弦矩阵误差,表示
Figure FDA00004826099100000511
第i行第j列元素值,
Figure FDA00004826099100000512
表示导航系下地球速度矢量,
Figure FDA00004826099100000513
分别表示导航系下地球速度矢量的北向分量和天向分量。
CN201410116682.2A 2014-03-26 2014-03-26 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法 Active CN103852085B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410116682.2A CN103852085B (zh) 2014-03-26 2014-03-26 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410116682.2A CN103852085B (zh) 2014-03-26 2014-03-26 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法

Publications (2)

Publication Number Publication Date
CN103852085A true CN103852085A (zh) 2014-06-11
CN103852085B CN103852085B (zh) 2016-09-21

Family

ID=50860026

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410116682.2A Active CN103852085B (zh) 2014-03-26 2014-03-26 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法

Country Status (1)

Country Link
CN (1) CN103852085B (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122412A (zh) * 2014-07-29 2014-10-29 北京机械设备研究所 一种基于北斗二代速度信息的加速度计标定方法
CN104677381A (zh) * 2015-01-29 2015-06-03 中国空空导弹研究院 一种微惯性测量单元的测试系统
CN104897178A (zh) * 2015-07-06 2015-09-09 中国人民解放军国防科学技术大学 一种双惯导联合旋转调制导航与在线相对性能评估方法
CN106884645A (zh) * 2015-12-16 2017-06-23 航天科工惯性技术有限公司 陀螺测斜仪的标定方法
CN107024674A (zh) * 2017-05-26 2017-08-08 北京航空航天大学 一种基于递推最小二乘法的磁强计现场快速标定方法
CN108132060A (zh) * 2017-11-17 2018-06-08 北京计算机技术及应用研究所 一种捷联惯导系统无基准的系统级标定方法
CN108398576A (zh) * 2018-03-06 2018-08-14 中国人民解放军火箭军工程大学 一种静态误差标定系统及方法
CN109883392A (zh) * 2019-03-08 2019-06-14 哈尔滨工程大学 一种基于相位补偿的捷联惯导升沉测量方法
CN110186479A (zh) * 2019-05-30 2019-08-30 北京航天控制仪器研究所 一种惯性器件误差系数确定方法
CN113551688A (zh) * 2021-05-27 2021-10-26 北京航天发射技术研究所 车载定位定向导航设备无依托快速免拆卸标定方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101021546A (zh) * 2007-03-26 2007-08-22 北京航空航天大学 一种光纤陀螺惯性测量单元的现场标定方法
CN101706287A (zh) * 2009-11-20 2010-05-12 哈尔滨工程大学 一种基于数字高通滤波的旋转捷联系统现场标定方法
WO2012021199A2 (en) * 2010-06-17 2012-02-16 The Aerospace Corporation High-frequency, hexapod six degree-of-freedom shaker
CN103245358A (zh) * 2012-06-05 2013-08-14 北京航空航天大学 一种光纤陀螺标度因数非对称性误差的系统级标定方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101021546A (zh) * 2007-03-26 2007-08-22 北京航空航天大学 一种光纤陀螺惯性测量单元的现场标定方法
CN101706287A (zh) * 2009-11-20 2010-05-12 哈尔滨工程大学 一种基于数字高通滤波的旋转捷联系统现场标定方法
WO2012021199A2 (en) * 2010-06-17 2012-02-16 The Aerospace Corporation High-frequency, hexapod six degree-of-freedom shaker
CN103245358A (zh) * 2012-06-05 2013-08-14 北京航空航天大学 一种光纤陀螺标度因数非对称性误差的系统级标定方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
刘百奇等: "光纤陀螺IMU的六位置旋转现场标定新方法", 《光电工程》 *

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104122412A (zh) * 2014-07-29 2014-10-29 北京机械设备研究所 一种基于北斗二代速度信息的加速度计标定方法
CN104677381A (zh) * 2015-01-29 2015-06-03 中国空空导弹研究院 一种微惯性测量单元的测试系统
CN104897178A (zh) * 2015-07-06 2015-09-09 中国人民解放军国防科学技术大学 一种双惯导联合旋转调制导航与在线相对性能评估方法
CN104897178B (zh) * 2015-07-06 2017-07-07 中国人民解放军国防科学技术大学 一种双惯导联合旋转调制导航与在线相对性能评估方法
CN106884645A (zh) * 2015-12-16 2017-06-23 航天科工惯性技术有限公司 陀螺测斜仪的标定方法
CN107024674B (zh) * 2017-05-26 2019-04-26 北京航空航天大学 一种基于递推最小二乘法的磁强计现场快速标定方法
CN107024674A (zh) * 2017-05-26 2017-08-08 北京航空航天大学 一种基于递推最小二乘法的磁强计现场快速标定方法
CN108132060A (zh) * 2017-11-17 2018-06-08 北京计算机技术及应用研究所 一种捷联惯导系统无基准的系统级标定方法
CN108398576A (zh) * 2018-03-06 2018-08-14 中国人民解放军火箭军工程大学 一种静态误差标定系统及方法
CN108398576B (zh) * 2018-03-06 2020-02-07 中国人民解放军火箭军工程大学 一种静态误差标定系统及方法
CN109883392A (zh) * 2019-03-08 2019-06-14 哈尔滨工程大学 一种基于相位补偿的捷联惯导升沉测量方法
CN110186479A (zh) * 2019-05-30 2019-08-30 北京航天控制仪器研究所 一种惯性器件误差系数确定方法
CN113551688A (zh) * 2021-05-27 2021-10-26 北京航天发射技术研究所 车载定位定向导航设备无依托快速免拆卸标定方法及装置

Also Published As

Publication number Publication date
CN103852085B (zh) 2016-09-21

Similar Documents

Publication Publication Date Title
CN103852085B (zh) 一种基于最小二乘拟合的光纤捷联惯导系统现场标定方法
CN103852086B (zh) 一种基于卡尔曼滤波的光纤捷联惯导系统现场标定方法
CN105371844B (zh) 一种基于惯性/天文互助的惯性导航系统初始化方法
CN105180968B (zh) 一种imu/磁强计安装失准角在线滤波标定方法
CN104344836B (zh) 一种基于姿态观测的冗余惯导系统光纤陀螺系统级标定方法
CN102538792B (zh) 一种位置姿态系统的滤波方法
CN103323026B (zh) 星敏感器和有效载荷的姿态基准偏差估计与修正方法
CN101706284B (zh) 提高船用光纤陀螺捷联惯导系统定位精度的方法
CN103245359B (zh) 一种惯性导航系统中惯性传感器固定误差实时标定方法
CN103076025B (zh) 一种基于双解算程序的光纤陀螺常值误差标定方法
CN104374388B (zh) 一种基于偏振光传感器的航姿测定方法
CN102564455B (zh) 星敏感器安装误差四位置标定与补偿方法
CN106017507A (zh) 一种用于中低精度的光纤惯组快速标定方法
CN103575299A (zh) 利用外观测信息的双轴旋转惯导系统对准及误差修正方法
CN106969783A (zh) 一种基于光纤陀螺惯性导航的单轴旋转快速标定技术
CN101246023A (zh) 微机械陀螺惯性测量组件的闭环标定方法
CN101290326A (zh) 石英挠性加速度计测量组件的参数辨识标定方法
CN104344837A (zh) 一种基于速度观测的冗余惯导系统加速度计系统级标定方法
CN102393210A (zh) 一种激光陀螺惯性测量单元的温度标定方法
CN103217159A (zh) 一种sins/gps/偏振光组合导航系统建模及动基座初始对准方法
CN104374401A (zh) 一种捷联惯导初始对准中重力扰动的补偿方法
CN106940193A (zh) 一种基于Kalman滤波的船舶自适应摇摆标定方法
CN104121927A (zh) 一种适用于低精度无方位基准单轴转位设备的惯性测量单元标定方法
CN103759729A (zh) 采用捷联惯导的月球软着陆地面试验用初始姿态获取方法
CN102679999A (zh) 星敏感器安装误差四位置标定与补偿方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant