CN103630146B

CN103630146B - The laser gyro IMU scaling method that a kind of discrete parsing is combined with Kalman filter

Info

Publication number: CN103630146B
Application number: CN201310419209.7A
Authority: CN
Inventors: 李建利; 焦峰; 王文建; 房建成; 钟麦英
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2013-09-15
Filing date: 2013-09-15
Publication date: 2016-07-20
Anticipated expiration: 2033-09-15
Also published as: CN103630146A

Abstract

The present invention provides a laser gyro IMU calibration method combining discrete analysis and Kalman filtering. The method combines the advantages of discrete analysis and system-level filtering. Firstly, according to the angular velocity and acceleration channel error model of the laser gyro IMU, four The azimuth rotation rate discrete analysis method preliminarily calibrates all 24 calibration coefficients of the angular velocity and acceleration channels; then uses the error value of the calibration coefficient as a state variable, and uses the Kalman filter method with lever arm effect compensation to estimate the error value of the calibration coefficient. And amended to obtain accurate calibration coefficients. This method eliminates the errors caused by random errors such as the deformation of the damping device and the drift of the gyro and accelerometer during the calibration process, improves the accuracy of the calibration coefficients, and reduces the error of the navigation solution.

Description

A Laser Gyro IMU Calibration Method Combining Discrete Analysis and Kalman Filter

技术领域technical field

本发明涉及一种离散解析与Kalman滤波结合的激光陀螺IMU标定方法，可以应用于精确标定位置姿态测量系统(PositionandOrientationSystem,POS)和惯性导航系统(InertialNavigationSystem,INS)，以及INS/GPS(GlobalPositionSystem,全球定位系统)组合导航系统的标定系数，提高导航解算的精度。The invention relates to a laser gyroscope IMU calibration method combining discrete analysis and Kalman filtering, which can be applied to accurately calibrate a position and orientation measurement system (Position and Orientation System, POS) and an inertial navigation system (Inertial Navigation System, INS), as well as INS/GPS (Global Position System, global Positioning system) combined with the calibration coefficient of the navigation system to improve the accuracy of the navigation solution.

背景技术Background technique

惯性测量单元(InertialMeasurementUnit,IMU)是POS、INS、INS/GPS等惯性测量系统的核心部件，惯性测量系统是以牛顿力学定律为基础，利用IMU测量载体的线运动和角运动参数，在给定初始条件下，由计算机推算得到载体的位置速度姿态参数。它是一种自主式的位置速度姿态测量单元，完全依靠机载设备自主地完成测量任务，和外界不发生任何光、电联系，具有隐蔽性好，工作不受气象条件限制的优点，因此在航空、航天和航海领域得到广泛的使用。IMU主要由陀螺仪和加速度计等惯性测量组件构成，陀螺仪的精度很大程度上决定了IMU的测量精度。由于激光陀螺具有高动态、低成本、高可靠、性能价格比高的特点，在军用、民用方面被广泛应用。但在实际应用中，激光陀螺的闭锁效应严重影响了其测量精度，通常采用抖动偏频的方法消除闭锁效应。由于存在机械抖动部件，利用激光陀螺构成的IMU中通常采用减振装置，将IMU内部的高频抖动和外界隔离。Inertial Measurement Unit (IMU) is the core component of inertial measurement systems such as POS, INS, INS/GPS, etc. The inertial measurement system is based on Newton's laws of mechanics, and uses IMU to measure the linear motion and angular motion parameters of the carrier. Under initial conditions, the position, velocity, and attitude parameters of the carrier are calculated by computer. It is an autonomous position, velocity and attitude measurement unit, which completely relies on the airborne equipment to complete the measurement task independently, and does not have any optical or electrical contact with the outside world. It has the advantages of good concealment, and its work is not restricted by weather conditions. It is widely used in the fields of aviation, spaceflight and navigation. The IMU is mainly composed of inertial measurement components such as gyroscopes and accelerometers. The accuracy of the gyroscope largely determines the measurement accuracy of the IMU. Due to the characteristics of high dynamics, low cost, high reliability, and high cost performance, laser gyroscopes are widely used in military and civilian applications. However, in practical applications, the blocking effect of the laser gyroscope seriously affects its measurement accuracy, and the method of dithering the bias frequency is usually used to eliminate the blocking effect. Due to the existence of mechanical jitter components, vibration damping devices are usually used in IMUs composed of laser gyro to isolate the high-frequency jitter inside the IMU from the outside world.

IMU的标定系数是影响POS、INS、INS/GPS等导航系统精度的关键因素。但是在激光陀螺IMU中，由于减振装置采用橡胶材料，在IMU标定翻转的过程中减振装置会因重力和外部输入角速度、加速度等而引起变形，从而导致标定精度的下降。此外，标定过程中陀螺和加速度计的随机漂移等随机误差也同样会引入标定误差。The calibration coefficient of IMU is a key factor affecting the accuracy of navigation systems such as POS, INS, and INS/GPS. However, in the laser gyro IMU, since the vibration damping device is made of rubber material, the vibration damping device will be deformed due to gravity and external input angular velocity, acceleration, etc. during the IMU calibration flipping process, resulting in a decrease in calibration accuracy. In addition, random errors such as random drift of the gyroscope and accelerometer during the calibration process will also introduce calibration errors.

发明内容Contents of the invention

本发明的技术解决问题是：提高带减振装置的激光陀螺IMU标定精度，提供一种离散解析与Kalman滤波结合的激光陀螺IMU标定方法，该方法修正了因减振装置变形、陀螺和加速度计漂移等随机误差而引起的标定误差，提高了标定系数的精度，减小了导航解算的误差。The technical solution of the present invention is to improve the calibration accuracy of the laser gyro IMU with a vibration damping device, and provide a laser gyro IMU calibration method combining discrete analysis and Kalman filtering. The calibration error caused by random errors such as drift improves the accuracy of the calibration coefficient and reduces the error of the navigation solution.

本发明的技术解决方案为：一种离散解析与Kalman滤波结合的激光陀螺IMU标定方法，其特点在于包括下列步骤：The technical solution of the present invention is: a kind of laser gyroscope IMU calibration method that discrete analysis combines with Kalman filter, it is characterized in that comprising the following steps:

步骤(1)、根据激光陀螺及加速度计的特性，建立激光陀螺IMU角速度和加速度通道的确定性误差模型，模型参数包括：x、y、z三轴陀螺输入角速度ω_x、ω_y、ω_z，三轴输出角速度G_x、G_y、G_z，三轴陀螺零偏G_bx、G_by、G_bz，三轴陀螺标度因数S_x、S_y、S_z，三轴陀螺随机误差δG_x、δG_y、δG_z，角速度通道的六个安装误差系数E_xy、E_xz、E_yx、E_yz、E_zx、E_zy；x、y、z三轴加速度计的输入加速度a_x、a_y、a_z，三轴加速度计的输出加速度A_x、A_y、A_z，三轴加速度计零偏A_bx、A_by、A_bz，三轴加速度计标度因数K_x、K_y、K_z，三轴加速度计随机误差δA_x、δA_y、δA_z，和加速度通道的六个安装误差系数M_xy、M_xz、M_yx、M_yz、M_zx、M_zy。Step (1), according to the characteristics of the laser gyro and accelerometer, establish a deterministic error model of the laser gyro IMU angular velocity and acceleration channel, the model parameters include: x, y, z three-axis gyro input angular velocity ω _x , ω _y , ω _z , three-axis output angular velocity G _x , G _y , G _z , three-axis gyro zero bias G _bx , G _by , G _bz , three-axis gyro scale factor S _x , S _y , S _z , three-axis gyro random error δG _x , δG _y , δG _z , the six installation error coefficients E _xy , E _xz , E _yx , E _yz , E _zx , E _zy of the angular velocity channel; the input acceleration a _x , a _y of the x, y, z three-axis accelerometer , a _z , the output acceleration A _x , A _y , A _z of the three-axis accelerometer, the zero bias of the three-axis accelerometer A _bx , A _by , A _bz , the scale factor of the three-axis accelerometer K _x , K _y , K _z , the random errors of the three-axis accelerometer δA _x , δA _y , δA _z , and the six installation error coefficients M _xy , M _xz , M _yx , M _yz , M _zx , M _zy of the acceleration channel.

步骤(2)、根据步骤(1)中建立的激光陀螺IMU角速度和加速度通道确定性误差模型，采用四方位旋转速率离散解析方法，初步标定出三轴陀螺零偏G_bx、G_by、G_bz，三轴陀螺标度因数S_x、S_y、S_z，角速度通道的六个安装误差系数E_xy、E_xz、E_yx、E_yz、E_zx、E_zy，三轴加速度计零偏A_bx、A_by、A_bz，三轴加速度计标度因数K_x、K_y、K_z，和加速度通道六个安装误差系数M_xy、M_xz、M_yx、M_yz、M_zx、M_zy共24个标定系数。Step (2), according to the laser gyro IMU angular velocity and acceleration channel deterministic error model established in step (1), using the four-axis rotation rate discrete analysis method, preliminarily calibrate the three-axis gyro zero bias G _bx , G _by , G _bz , three-axis gyro scale factors S _x , S _y , S _z , six installation error coefficients E _xy , E _xz , E _yx , E _yz , E _zx , E _zy of the angular velocity channel, three-axis accelerometer zero bias A _bx _. _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} calibration factor.

步骤(3)、建立激光陀螺IMU角速度和加速度通道的随机误差模型，模型参数包括：三轴陀螺随机误差δG_x、δG_y、δG_z，三轴陀螺随机游走ε_x、ε_y、ε_z、三轴陀螺随机游走驱动噪声三轴陀螺标度因数误差δS_x、δS_y、δS_z，角速度通道六个安装误差系数偏差δE_xy、δE_xz、δE_yx、δE_yz、δE_zx、δE_zy；三轴加速度计随机误差δA_x、δA_y、δA_z，三轴加速度计随机游走三轴加速度计随机游走驱动噪声三轴加速度计标度因数误差δK_x、δK_y、δK_z，加速度通道安装误差系数偏差δM_xy、δM_xz、δM_yx、δM_yz、δM_zx、δM_zy。根据随机误差模型和惯性导航解算的基本原理，建立Kalman滤波的状态方程。Step (3), establish the random error model of the laser gyro IMU angular velocity and acceleration channel, the model parameters include: three-axis gyro random error δG _x , δG _y , δG _z , three-axis gyro random walk ε _x , ε _y , ε _z , three-axis gyroscope random walk driving noise Three-axis gyro scale factor error δS _x , δS _y , δS _z , six installation error coefficient deviations of angular velocity channel δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy ; random error of three-axis accelerometer δA _x , δA _y , δA _z , three-axis accelerometer random walk Triaxial Accelerometer Random Walk Drive Noise Three-axis accelerometer scale factor error δK _x , δK _y , δK _z , acceleration channel installation error coefficient deviation δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy . According to the random error model and the basic principle of inertial navigation solution, the state equation of Kalman filter is established.

步骤(4)、选取东、北、天向的速度误差δV_E、δV_N、δV_U作为观测量，建立带杆臂效应补偿的Kalman滤波量测方程。Step (4), select the velocity errors δV _E , δV _N , and δV _U in the east, north, and sky directions as observations, and establish a Kalman filter measurement equation with lever-arm effect compensation.

步骤(5)、用Kalman滤波估计步骤(3)中24个随机误差模型参数ε_x、ε_y、ε_z、δS_x、δS_y、δS_z、δE_xy、δE_xz、δE_yx、δE_yz、δE_zx、δE_zy、δK_x、δK_y、δK_z、δM_xy、δM_xz、δM_yx、δM_yz、δM_zx、δM_zy，并与步骤(2)中的24个标定系数对应相加进行修正。重复此步骤，直到所有的标定系数都收敛为止，收敛条件为连续十分钟内24个标定系数在修正前后相对变化小于0.05％。此时得到精确的标定系数，用于补偿激光陀螺IMU数据。Step (5), use Kalman filter to estimate 24 random error model parameters in step (3) ε _x , ε _y , ε _z , δS _x , δS _y , δS _z , δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy , δK _x , δK _y , δK _z , δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy , and correspondingly added to the 24 calibration coefficients in step (2) for correction. Repeat this step until all the calibration coefficients converge, and the convergence condition is that the relative change of 24 calibration coefficients before and after correction within ten minutes is less than 0.05%. At this time, accurate calibration coefficients are obtained, which are used to compensate the laser gyro IMU data.

步骤(2)所述的四方位旋转速率离散解析方法中，四方位是将激光陀螺IMU分别调整到x、y、z、-x轴分别指天，其余两轴水平；当x、y、z轴分别指天时，转台的外框轴逆时针方向以5°/s—100°/s之间的任一恒定角速率匀速旋转两周；当-x指天时，转台的外框轴顺时针方向以5°/s—100°/s之间的任一恒定角速率匀速旋转两周。In the four-direction rotation rate discrete analysis method described in step (2), the four-direction is to adjust the laser gyro IMU to the x, y, z, and -x axes respectively pointing to the sky, and the remaining two axes are horizontal; when x, y, z When the axes point to the sky, the outer frame axis of the turntable rotates counterclockwise for two rounds at a constant angular rate between 5°/s and 100°/s; when -x points to the sky, the outer frame axis of the turntable rotates clockwise Rotate at a constant speed for two rounds at any constant angular rate between 5°/s and 100°/s.

步骤(3)所述的Kalman滤波状态方程中，状态变量共有30维，包括东、北、天向的速度误差δV_E、δV_N、δV_U，俯仰、横滚、航向角误差φ_E、φ_N、φ_U，以及三轴陀螺随机游走ε_x、ε_y、ε_z，陀螺标度因数误差δS_x、δS_y、δS_z，角速度通道六个安装误差系数偏差δE_xy、δE_xz、δE_yx、δE_yz、δE_zx、δE_zy，三轴加速度计随机游走三轴加速度计标度因数误差δK_x、δK_y、δK_z，加速度通道安装误差系数偏差δM_xy、δM_xz、δM_yx、δM_yz、δM_zx、δM_zy。In the Kalman filter state equation described in step (3), the state variables have a total of 30 dimensions, including velocity errors δV _E , δV _N , and δV _U in the east, north, and sky directions, pitch, roll, and heading angle errors φ _E , φ _N , φ _U , and three-axis gyro random walk ε _x , ε _y , ε _z , gyro scaling factor errors δS _x , δS _y , δS _z , six installation error coefficient deviations of angular velocity channels δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy , three-axis accelerometer random walk Three-axis accelerometer scale factor error δK _x , δK _y , δK _z , acceleration channel installation error coefficient deviation δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy .

步骤(5)所述的带杆臂效应补偿的Kalman滤波量测方程中，量测量Z是经杆臂效应补偿后的激光陀螺IMU东、北、天向的速度误差δV_E、δV_N、δV_U。In the Kalman filter measurement equation with lever-arm effect compensation described in step (5), the measurement Z is the speed error δV _E , δV _N , δV of the laser gyro IMU east, north, and celestial directions after the lever-arm effect compensation _U.

$Z Z = = [\begin{matrix} δ δ {V V}_{E E.} \\ {δV δV}_{N N} \\ {δV δV}_{U u} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} \\ {V V}_{N N} \\ {V V}_{U u} \end{matrix}] - - [\begin{matrix} d d {V V}_{l l E E.}^{n no} \\ {dV dV}_{l l N N}^{n no} \\ {dV dV}_{l l U u}^{n no} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} - - d d {V V}_{l l E E.}^{n no} \\ {V V}_{N N} - - {dV dV}_{l l N N}^{n no} \\ {V V}_{U u} - - {dV dV}_{l l U u}^{n no} \end{matrix}]$

其中，V_E、V_N、V_U分别是激光陀螺IMU东、北、天向的速度，分别为杆臂效应在东、北、天向附加的速度误差。Among them, V _E , V _N , and V _U are the speeds of the laser gyro IMU in the east, north, and sky directions, respectively, are the additional velocity errors of the lever-arm effect in the east, north, and sky directions, respectively.

本发明的原理是：四方位旋转速率离散解析方法是根据激光陀螺IMU的角速度和加速度通道确定性误差模型，建立包含24个标定系数的多元线性方程组，利用求解多元线性方程组解的原理计算激光陀螺IMU的所有24个标定系数；而Kalman滤波方法则是利用递推最小方差估计的原理，从量测信息中消除高斯白噪声，提取出被估计参数的信息，进而修正四方位旋转速率离散解析方法初步计算出的标定系数，最终提高标定系数精度。The principle of the present invention is: the four-direction rotation rate discrete analysis method is based on the deterministic error model of the angular velocity and acceleration channel of the laser gyro IMU, establishes a multivariate linear equation group containing 24 calibration coefficients, and uses the principle of solving the multivariate linear equation group solution to calculate All 24 calibration coefficients of the laser gyro IMU; and the Kalman filter method uses the principle of recursive minimum variance estimation to eliminate Gaussian white noise from the measurement information, extract the information of the estimated parameters, and then correct the dispersion of the four-direction rotation rate The calibration coefficients initially calculated by the analytical method can finally improve the accuracy of the calibration coefficients.

本发明与现有技术相比的优点在于：本发明在离散解析标定方法的基础上，考虑了标定过程中的减振装置变形、陀螺和加速度计漂移等随机误差，并采用Kalman滤波方法估计并消除了这些随机误差，实现了带减振装置的激光陀螺IMU高精度标定，提高了标定系数的精度，减小了导航解算的误差。Compared with the prior art, the present invention has the advantages that: on the basis of the discrete analytical calibration method, the present invention considers random errors such as vibration damping device deformation, gyroscope and accelerometer drift in the calibration process, and uses the Kalman filter method to estimate and These random errors are eliminated, and the high-precision calibration of the laser gyro IMU with a vibration reduction device is realized, the accuracy of the calibration coefficient is improved, and the error of the navigation solution is reduced.

附图说明Description of drawings

图1为本发明的原理框图。Fig. 1 is a functional block diagram of the present invention.

图2为本发明的四方位旋转速率离散解析标定方案。Fig. 2 is a four-direction rotation rate discrete analytical calibration scheme of the present invention.

图3为本发明的步骤(5)的Kalman滤波及修正标定系数流程图。Fig. 3 is the flow chart of Kalman filtering and correcting calibration coefficient in step (5) of the present invention.

图4为本发明所采用的Kalman滤波算法原理图。FIG. 4 is a schematic diagram of the Kalman filtering algorithm adopted in the present invention.

具体实施方式detailed description

如图1所示，本发明的具体方法如下：As shown in Figure 1, the concrete method of the present invention is as follows:

(1)根据激光陀螺及加速度计的特性，建立激光陀螺IMU角速度和加速度通道的误差模型。(1) According to the characteristics of the laser gyroscope and accelerometer, the error model of the angular velocity and acceleration channels of the laser gyroscope IMU is established.

角速度通道误差模型为：The angular velocity channel error model is:

$\begin{matrix} {G G}_{x x} = = {G G}_{b b x x} + + {S S}_{x x} {ω ω}_{x x} + + {E E.}_{x x y the y} {ω ω}_{y the y} + + {E E.}_{x x z z} {ω ω}_{z z} + + {δG δG}_{x x} \\ {G G}_{y the y} = = {G G}_{b b y the y} + + {E E.}_{y the y x x} {ω ω}_{x x} + + {S S}_{y the y} {ω ω}_{y the y} + + {E E.}_{y the y z z} {ω ω}_{z z} + + {δG δ G}_{y the y} \\ {G G}_{z z} = = {G G}_{b b z z} + + {E E.}_{z z x x} {ω ω}_{x x} + + {E E.}_{z z y the y} {ω ω}_{y the y} + + {S S}_{z z} {ω ω}_{z z} + + {δG δ G}_{z z} \end{matrix} - - - - - - ((11))$

其中，ω_x、ω_y、ω_z为x、y、z三轴陀螺输入角速度，G_x、G_y、G_z为三轴输出角速度，G_bx、G_by、G_bz为三轴陀螺零偏，S_x、S_y、S_z为陀螺标度因数，δG_x、δG_y、δG_z为三轴陀螺随机误差，E_xy、E_xz、E_yx、E_yz、E_zx、E_zy为角速度通道的六个安装误差系数。Among them, ω _x , ω _y , ω _z are the input angular velocity of the x, y, z three-axis gyroscope, G _x , G _y , G _z are the three-axis output angular velocity, G _bx , G _by , G _bz are the zero bias of the three-axis gyroscope , S _x , S _y , S _z are gyro scaling factors, δG _x , δG _y , δG _z are three-axis gyro random errors, E _xy , E _xz , E _yx , E _yz , E _zx , E _zy are angular velocity channels The six installation error coefficients.

加速度通道误差模型为：The acceleration channel error model is:

$\begin{matrix} {A A}_{x x} = = {A A}_{b b x x} + + {K K}_{x x} {a a}_{x x} + + {M m}_{x x y the y} {a a}_{y the y} + + {M m}_{x x z z} {a a}_{z z} + + {δA δA}_{x x} \\ {A A}_{y the y} = = {A A}_{b b x x} + + {M m}_{y the y x x} {a a}_{x x} + + {K K}_{y the y} {a a}_{y the y} + + {M m}_{y the y z z} {a a}_{z z} + + {δA δA}_{y the y} \\ {A A}_{z z} = = {A A}_{b b z z} + + {M m}_{z z x x} {a a}_{x x} + + {M m}_{z z y the y} {a a}_{y the y} + + {K K}_{z z} {a a}_{z z} + + {δA δA}_{z z} \end{matrix} - - - - - - ((22))$

其中，a_x、a_y、a_z为x、y、z三轴输入加速度，A_x、A_y、A_z为三轴输出加速度，A_bx、A_by、A_bz为三轴加速度计零偏，K_x、K_y、K_z为标度因数，δA_x、δA_y、δA_z为三轴加速度计随机误差，M_xy、M_xz、M_yx、M_yz、M_zx、M_zy为加速度通道的六个安装误差系数。Among them, a _x , a _y , a _z are three-axis input accelerations of x, y, z, A _x , A _y , A _z are three-axis output accelerations, A _bx , A _by , A _bz are three-axis accelerometer zero bias , K _x , K _y , K _z are scaling factors, δA _x , δA _y , δA _z are random errors of the three-axis accelerometer, M _xy , M _xz , M _yx , M _yz , M _zx , M _zy are acceleration channels The six installation error coefficients.

(2)根据步骤(1)中建立的激光陀螺IMU角速度和加速度通道误差模型，采用四方位旋转速率离散解析方法，初步标定出三轴陀螺零偏G_bx、G_by、G_bz，陀螺标度因数S_x、S_y、S_z，角速度通道的六个安装误差系数E_xy、E_xz、E_yx、E_yz、E_zx、E_zy，三支加速度计零偏A_bx、A_by、A_bz，三支加速度计标度因数K_x、K_y、K_z，和加速度通道六个安装误差系数M_xy、M_xz、M_yx、M_yz、M_zx、M_zy共24个标定系数。具体步骤如下：(2) According to the error model of the laser gyro IMU angular velocity and acceleration channel established in step (1), the four-axis rotation rate discrete analysis method is used to preliminarily calibrate the three-axis gyro zero bias G _bx , G _by , G _bz , and the gyro scale Factors S _x , S _y , S _z , six installation error coefficients E _xy , E _xz , E _yx , E _yz , E _zx , E _zy of the angular velocity channel, three accelerometer zero biases A _bx , A _by , A _bz , three accelerometer scale factors K _x , K _y , K _z , and six installation error coefficients M _xy , M _xz , M _yx , M _yz , M _zx , M _zy of the acceleration channels, a total of 24 calibration coefficients. Specific steps are as follows:

(a)将激光陀螺IMU通过一个专用的过渡板，与三轴转台的内框轴固连。如图2中4个方位所示，将激光陀螺IMU分别调整到x、-x、y、z轴分别指天，与转台旋转轴z_p重合，其余两轴水平。当如方位1所示x轴指天、如方位3所示y轴指天以及如方位4所示z轴指天时，转台的外框轴逆时针方向以恒定角速度Ω(5°/s＜Ω＜100°/s，通常取10°/s)匀速旋转两周；当如图2中方位2所示-x指天时，转台的外框轴顺时针方向以恒定角速度Ω(5°/s＜Ω＜100°/s，通常取10°/s)匀速旋转两周。(a) Connect the laser gyro IMU to the inner frame shaft of the three-axis turntable through a special transition plate. As shown in the four orientations in Figure 2, adjust the laser gyro IMU to point to the sky with the x, -x, y, and z axes respectively, coincide with the turntable rotation axis z _p , and the other two axes are horizontal. When the x-axis points to the sky as shown in position 1, the y-axis points to the sky as shown in position 3, and the z-axis points to the sky as shown in position 4, the outer frame axis of the turntable moves counterclockwise at a constant angular velocity Ω (5°/s<Ω <100°/s, usually take 10°/s) to rotate at a constant speed for two weeks; when -x refers to the sky as shown in the position 2 in Figure 2, the outer frame axis of the turntable is clockwise at a constant angular velocity Ω (5°/s<Ω<100°/s, usually 10°/s) rotates at a constant speed for two weeks.

(b)根据上述方案可建立如下方程：(b) According to the above scheme, the following equation can be established:

角速度通道方程:Angular velocity channel equation:

$[\begin{matrix} {G G}_{x x 11} & {G G}_{y the y 11} & {G G}_{z z 11} \\ {G G}_{x x 22} & {G G}_{y the y 22} & {G G}_{z z 22} \\ {G G}_{x x 33} & {G G}_{y the y 33} & {G G}_{z z 33} \\ {G G}_{x x 44} & {G G}_{y the y 44} & {G G}_{z z 44} \end{matrix}] = = [\begin{matrix} 11 & \overset{&OverBar; &OverBar;}{ω ω} & 00 & 00 \\ 11 & - - \overset{&OverBar; &OverBar;}{ω ω} & 00 & 00 \\ 11 & 00 & \overset{&OverBar; &OverBar;}{ω ω} & 00 \\ 11 & 00 & 00 & \overset{&OverBar; &OverBar;}{ω ω} \end{matrix}] [\begin{matrix} {G G}_{b b x x} & {G G}_{b b y the y} & {G G}_{b b z z} \\ {S S}_{x x} & {E E.}_{x x y the y} & {E E.}_{x x z z} \\ {E E.}_{y the y x x} & {S S}_{y the y} & {E E.}_{y the y z z} \\ {E E.}_{z z x x} & {E E.}_{z z y the y} & {S S}_{z z} \end{matrix}] - - - - - - ((33))$

其中,G_x1、G_y1、G_z1为第1个方位下惯性系统x、y、z三轴陀螺的输出，G_x2、G_y2、G_z2为第2个方位下三轴陀螺的输出，G_x3、G_y3、G_z3为第3个方位下三轴陀螺的输出，G_x4、G_y4、G_z4为第4个方位下三轴陀螺的输出，Ω为转台输入角速率，ω_ie为地球自转角速率(°/h)，φ为当地纬度。in, G _x1 , G _y1 , G _z1 are the outputs of the inertial system x, y, z three-axis gyroscope in the first azimuth, G _x2 , G _y2 , G _z2 are the outputs of the three-axis gyroscope in the second azimuth, G _x3 , G _y3 , G _z3 are the output of the three-axis gyro in the third azimuth, G _x4 , G _y4 , G _z4 are the output of the three-axis gyro in the fourth azimuth, Ω is the input angular rate of the turntable, ω _ie is the rotation angle of the earth Speed (°/h), φ is the local latitude.

加速度通道方程:Acceleration channel equation:

$[\begin{matrix} {A A}_{x x 11} & {A A}_{y the y 11} & {A A}_{z z 11} \\ {A A}_{x x 22} & {A A}_{y the y 22} & {A A}_{z z 22} \\ {A A}_{x x 33} & {A A}_{y the y 33} & {A A}_{z z 33} \\ {A A}_{x x 44} & {A A}_{y the y 44} & {A A}_{z z 44} \end{matrix}] = = [\begin{matrix} 11 & g g & 00 & 00 \\ 11 & - - g g & 00 & 00 \\ 11 & 00 & g g & 00 \\ 11 & 00 & 00 & g g \end{matrix}] [\begin{matrix} {A A}_{b b x x} & {A A}_{b b y the y} & {A A}_{b b z z} \\ {K K}_{x x} & {M m}_{y the y x x} & {M m}_{z z x x} \\ {M m}_{x x y the y} & {K K}_{y the y} & {M m}_{z z y the y} \\ {M m}_{x x z z} & {M m}_{y the y z z} & {K K}_{z z} \end{matrix}] - - - - - - ((44))$

其中，A_x1、A_y1、A_z1为第1个方位下惯性系统x、y、z三轴加速度计的输出，A_x2、A_y2、A_z2为第2个方位下三轴加速度计的输出，A_x3、A_y3、A_z3为第3个方位下三轴加速度计的输出，A_x4、A_y4、A_z4为第4个方位下三轴加速度计的输出，g为当地重力加速度。Among them, A _x1 , A _y1 , A _z1 are the output of the inertial system x, y, z three-axis accelerometer in the first orientation, A _x2 , A _y2 , A _z2 are the output of the three-axis accelerometer in the second orientation , A _x3 , A _y3 , A _z3 are the outputs of the triaxial accelerometer in the third orientation, A _x4 , A _y4 , A _z4 are the outputs of the triaxial accelerometer in the fourth orientation, and g is the local gravity acceleration.

(c)求解角速度和加速度通道方程，可得标定系数。(c) Solve the angular velocity and acceleration channel equations to obtain the calibration coefficients.

角速度通道标定系数：Angular velocity channel calibration coefficient:

$[\begin{matrix} {G G}_{b b x x} & {S S}_{x x} \\ {G G}_{b b y the y} & {E E.}_{y the y x x} \\ {G G}_{b b z z} & {E E.}_{z z x x} \end{matrix}] = = \frac{11}{22} [\begin{matrix} {G G}_{x x 22} + + {G G}_{x x 11} & (({G G}_{x x 11} - - {G G}_{x x 22})) / / \overset{&OverBar; &OverBar;}{ω ω} \\ {G G}_{y the y 22} + + {G G}_{y the y 11} & (({G G}_{y the y 11} - - {G G}_{y the y 22})) / / \overset{&OverBar; &OverBar;}{ω ω} \\ {G G}_{z z 22} + + {G G}_{z z 11} & (({G G}_{z z 11} - - {G G}_{z z 22})) / / \overset{&OverBar; &OverBar;}{ω ω} \end{matrix}] - - - - - - ((55))$

$[\begin{matrix} {E E.}_{x x y the y} & {E E.}_{x x z z} \\ {S S}_{y the y} & {E E.}_{y the y z z} \\ {E E.}_{z z y the y} & {S S}_{z z} \end{matrix}] = = \frac{11}{\overset{&OverBar; &OverBar;}{ω ω}} [\begin{matrix} {G G}_{x x 33} - - {G G}_{b b x x} & {G G}_{x x 44} - - {G G}_{b b x x} \\ {G G}_{y the y 33} - - {G G}_{b b y the y} & {G G}_{y the y 44} - - {G G}_{b b y the y} \\ {G G}_{z z 33} - - {G G}_{b b z z} & {G G}_{z z 44} - - {G G}_{b b z z} \end{matrix}] - - - - - - ((66))$

加速度通道标定系数：Acceleration channel calibration factor:

$[\begin{matrix} {A A}_{b b x x} & {A A}_{b b y the y} & {A A}_{b b z z} \\ {K K}_{x x} & {M m}_{y the y x x} & {M m}_{z z x x} \\ {M m}_{x x y the y} & {K K}_{y the y} & {M m}_{z z y the y} \\ {M m}_{x x z z} & {M m}_{y the y z z} & {K K}_{z z} \end{matrix}] = = \frac{11}{22 g g} [\begin{matrix} g g (({A A}_{x x 11} + + {A A}_{x x 22})) & g g (({A A}_{y the y 11} + + {A A}_{y the y 22})) & g g (({A A}_{z z 11} + + {A A}_{z z 22})) \\ {A A}_{x x 11} - - {A A}_{x x 22} & {A A}_{y the y 11} - - {A A}_{y the y 22} & {A A}_{z z 11} - - {A A}_{z z 22} \\ 22 (({A A}_{x x 33} - - {A A}_{b b x x})) & 22 (({A A}_{y the y 33} - - {A A}_{b b y the y})) & 22 (({A A}_{z z 33} - - {A A}_{b b z z})) \\ 22 (({A A}_{x x 44} - - {A A}_{b b x x})) & 22 (({A A}_{y the y 44} - - {A A}_{b b y the y})) & 22 (({A A}_{z z 44} - - {A A}_{b b z z})) \end{matrix}] - - - - - - ((77))$

(3)建立激光陀螺IMU角速度和加速度通道的随机误差模型，并根据随机误差模型和惯性导航解算的基本原理，建立Kalman滤波的状态方程。(3) Establish the random error model of the laser gyro IMU angular velocity and acceleration channel, and establish the state equation of the Kalman filter according to the basic principle of the random error model and inertial navigation solution.

(a)激光陀螺IMU角速度和加速度通道的随机误差模型。(a) Random error model of laser gyro IMU angular velocity and acceleration channels.

角速度通道的随机误差模型：Random error model of angular velocity channel:

$\begin{matrix} {δG δG}_{x x} = = {δS δS}_{x x} {ω ω}_{x x} + + {δE δE}_{x x y the y} {ω ω}_{y the y} + + {δE δE}_{x x z z} {ω ω}_{z z} + + {ϵ ϵ}_{x x} + + {w w}_{{ϵ ϵ}_{x x}} \\ {δG δ G}_{y the y} = = {δE δE}_{y the y x x} {ω ω}_{x x} + + {δS δS}_{y the y} {ω ω}_{y the y} + + {δE δE}_{y the y z z} {ω ω}_{z z} + + {ϵ ϵ}_{y the y} + + {w w}_{{ϵ ϵ}_{y the y}} \\ {δG δ G}_{z z} = = {δE δE}_{z z x x} {ω ω}_{x x} + + {δE δE}_{z z y the y} {ω ω}_{y the y} + + {δS δS}_{z z} {ω ω}_{z z} + + {ϵ ϵ}_{z z} + + {w w}_{{ϵ ϵ}_{z z}} \end{matrix} - - - - - - ((88))$

其中，δG_x、δG_y、δG_z表示三轴陀螺随机误差，ε_x、ε_y、ε_z为三轴陀螺随机游走，w_εx、w_εy、w_εz为随机游走驱动噪声、δS_x、δS_y、δS_z为陀螺标度因数误差，δE_xy、δE_xz、δE_yx、δE_yz、δE_zx、δE_zy角速度通道六个安装误差系数偏差。Among them, δG _x , δG _y , δG _z represent the random error of the three-axis gyroscope, ε _x , ε _y , ε _z are the random walk of the three-axis gyro, w _εx , w _εy , w _εz are the driving noise of the random walk, δS _x , δS _y , δS _z are gyroscope scale factor errors, and δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy angular velocity channels have six installation error coefficient deviations.

加速度通道的随机误差模型：The random error model of the acceleration channel:

$\begin{matrix} {δA δA}_{x x} = = {δK δK}_{x x} {a a}_{x x} + + {δM δ M}_{x x y the y} {a a}_{y the y} + + {δM δ M}_{x x z z} {a a}_{z z} + + {&dtri; &dtri;}_{x x} + + {w w}_{{&dtri; &dtri;}_{x x}} \\ {δA δA}_{y the y} = = {δM δ M}_{y the y x x} {a a}_{x x} + + {δK δK}_{y the y} {a a}_{y the y} + + {δM δ M}_{y the y z z} {a a}_{z z} + + {&dtri; &dtri;}_{y the y} + + {w w}_{{&dtri; &dtri;}_{y the y}} \\ {δA δA}_{z z} = = {δM δ M}_{z z x x} {a a}_{x x} + + {δM δ M}_{z z y the y} {a a}_{y the y} + + {δK δK}_{z z} {a a}_{z z} + + {&dtri; &dtri;}_{z z} + + {w w}_{{&dtri; &dtri;}_{z z}} \end{matrix} - - - - - - ((99))$

其中，δA_x、δA_y、δA_z表示三轴加速度随机误差，为加速度计随机游走，为随机游走驱动噪声，δK_x、δK_y、δK_z为加速度计标度因数误差，δM_xy、δM_xz、δM_yx、δM_yz、δM_zx、δM_zy为加速度通道安装误差系数偏差。Among them, δA _x , δA _y , δA _z represent the random error of the three-axis acceleration, is the accelerometer random walk, is the random walk driving noise, δK _x , δK _y , and δK _z are the accelerometer scale factor errors, and δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , and δM _zy are the deviations of the acceleration channel installation error coefficients.

(b)根据激光陀螺随机误差模型以及惯性导航解算的基本原理，建立Kalman滤波的状态方程。(b) According to the random error model of laser gyroscope and the basic principle of inertial navigation solution, the state equation of Kalman filter is established.

选取东、北、天向的速度误差δV_E、δV_N、δV_U，俯仰、横滚、航向角误差φ_E、φ_N、φ_U，以及步骤(3)中随机误差模型参数ε_x、ε_y、ε_z、δS_x、δS_y、δS_z、δE_xy、δE_xz、δE_yx、δE_yz、δE_zx、δE_zy、δK_x、δK_y、δK_z、δM_xy、δM_xz、δM_yx、δM_yz、δM_zx、δM_zy为状态变量X，X的维数为30维。Select the velocity errors δV _E , δV _N , and δV _U in the east, north, and sky directions, the pitch, roll, and heading angle errors φ _E , φ _N , and φ _U , and the random error model parameters ε _x , ε in step (3) _y , ε _z , δS _x , δS _y , δS _z , δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy , δK _x , δK _y , δK _z , δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , and δM _zy are state variables X, and the dimension of X is 30 dimensions.

状态方程为：The state equation is:

$\overset{\cdot &Center Dot;}{X x} = = F f \cdot &Center Dot; X x + + G G \cdot &Center Dot; w w - - - - - - ((1010))$

其中，X为系统状态变量：Among them, X is the system state variable:

$\begin{matrix} X x = = [[\begin{matrix} {δV δV}_{E E.} & {δV δV}_{N N} & {δV δV}_{U u} & {φ φ}_{E E.} & {φ φ}_{N N} & {φ φ}_{U u} & {&dtri; &dtri;}_{x x} & {&dtri; &dtri;}_{y the y} & {&dtri; &dtri;}_{z z} & {δK δK}_{x x} & {δK δK}_{y the y} & {δK δK}_{z z} & {δM δ M}_{x x y the y} & {δM δ M}_{x x z z} & {δM δ M}_{y the y x x} & {δM δ M}_{y the y z z} & {δM δ M}_{z z x x} & {δM δ M}_{z z y the y} \end{matrix} \\ \begin{matrix} {ϵ ϵ}_{x x} & {ϵ ϵ}_{y the y} & {ϵ ϵ}_{z z} & {δS δS}_{x x} & {δS δS}_{y the y} & {δS δS}_{z z} & {δE δE}_{x x y the y} & {δE δE}_{x x z z} & {δE δE}_{y the y x x} & {δE δE}_{y the y z z} & {δE δE}_{z z x x} & {δE δE}_{z z y the y} {]]}^{T T} \end{matrix} \end{matrix} - - - - - - ((1111))$

w为系统噪声矢量：w is the system noise vector:

w＝[w_axw_ayw_azw_gxw_gyw_gz]^T(12)w＝[w _ax w _ay w _az w _gx w _gy w _gz ] ^T (12)

G为噪声阵：G is the noise matrix:

$G G = = [\begin{matrix} {C C}_{b b}^{n no} & 00_{33 \times \times 33} \\ 00_{33 \times \times 33} & {C C}_{b b}^{n no} \\ 00_{24 twenty four \times \times 33} & 00_{24 twenty four \times \times 33} \end{matrix}] - - - - - - ((1313))$

是从IMU载体系到导航系的姿态转移矩阵，分别用θ、φ、γ表示IMU的航向、俯仰、横滚角，则: is the attitude transfer matrix from the IMU carrier system to the navigation system, and θ, φ, γ represent the heading, pitch, and roll angles of the IMU respectively, then:

${C C}_{b b}^{n no} = = [\begin{matrix} cos cos γ γ cos cos φ φ + + sin sin γ γ sin sin θ θ sin sin φ φ & cos cos θ θ sin sin φ φ & sin sin γ γ cos cos φ φ - - cos cos γ γ sin sin θ θ sin sin φ φ \\ - - cos cos γ γ sin sin φ φ + + sin sin γ γ sin sin θ θ cos cos φ φ & cos cos θ θ cos cos φ φ & - - sin sin γ γ sin sin φ φ - - cos cos γ γ sin sin θ θ cos cos φ φ \\ - - sin sin γ γ cos cos θ θ & sin sin θ θ & cos cos γ γ cos cos θ θ \end{matrix}] = = [\begin{matrix} {T T}_{1111} & {T T}_{1212} & {T T}_{1313} \\ {T T}_{21 twenty one} & {T T}_{22 twenty two} & {T T}_{23 twenty three} \\ {T T}_{3131} & {T T}_{3232} & {T T}_{3333} \end{matrix}] - - - - - - ((1414))$

T_ij(i,j＝1,2,3)表示姿态转移矩阵的9个元素。T _ij (i,j=1,2,3) represents the attitude transfer matrix of 9 elements.

F为系统转移矩阵：F is the system transition matrix:

$F f = = [\begin{matrix} {A A}_{33 \times \times 66} & {C C}_{33 \times \times 1212} & 00_{33 \times \times 1212} \\ {B B}_{33 \times \times 66} & 00_{33 \times \times 1212} & {D D.}_{33 \times \times 1212} \\ 00_{24 twenty four \times \times 66} & 00_{24 twenty four \times \times 1212} & 00_{24 twenty four \times \times 1212} \end{matrix}] - - - - - - ((1515))$

${A A}_{33 \times \times 66} = = [\begin{matrix} \frac{{V V}_{N N}}{R R} tan the tan L L & 22 {ω ω}_{i i e e} sin sin L L + + \frac{{V V}_{E E.} tan the tan L L}{R R} & - - 22 {ω ω}_{i i e e} cos cos L L - - \frac{{V V}_{E E.}}{R R} & 00 & g g & 00 \\ - - 22 {ω ω}_{i i e e} sin sin L L + + \frac{22 {V V}_{E E.} tan the tan L L}{R R} & - - \frac{{V V}_{U u}}{R R} & - - \frac{{V V}_{N N}}{R R} & - - g g & 00 & 00 \\ 22 {ω ω}_{i i e e} cos cos L L + + \frac{22 {V V}_{E E.}}{R R} & \frac{22 {V V}_{N N}}{R R} & 00 & 00 & 00 & 00 \end{matrix}] - - - - - - ((1616))$

其中，V_E、V_N、V_U分别是激光陀螺IMU东、北、天向的速度，L是当地纬度，R为地球半径，ω_ie为地球自转角速率，g为地球重力加速度。Among them, V _E , V _N , and V _U are the velocities of the laser gyro IMU in the east, north, and sky directions respectively, L is the local latitude, R is the radius of the earth, ω _ie is the angular rate of the earth's rotation, and g is the acceleration of gravity of the earth.

${B B}_{33 \times \times 66} = = [\begin{matrix} 00 & - - \frac{11}{R R} & 00 & 00 & {ω ω}_{i i e e} sin sin L L + + \frac{{V V}_{E E.} tan the tan L L}{R R} & - - {ω ω}_{i i e e} cos cos L L - - \frac{{V V}_{E E.}}{R R} \\ \frac{11}{R R} & 00 & 00 & - - {ω ω}_{i i e e} sin sin L L - - \frac{{V V}_{E E.} tan the tan L L}{R R} & 00 & - - \frac{{V V}_{N N}}{R R} \\ \frac{11}{R R} tan the tan L L & 00 & 00 & {ω ω}_{i i e e} cos cos L L + + \frac{{V V}_{E E.}}{R R} & \frac{{V V}_{N N}}{R R} & 00 \end{matrix}] - - - - - - ((1717))$

${C C}_{33 \times \times 1212} = = [\begin{matrix} {T T}_{1111} & {T T}_{1212} & {T T}_{1313} & {T T}_{1111} {a a}_{x x} & {T T}_{1212} {a a}_{y the y} & {T T}_{1313} {a a}_{z z} & {T T}_{1111} {a a}_{y the y} & {T T}_{1111} {a a}_{z z} & {T T}_{1212} {a a}_{x x} & {T T}_{1212} {a a}_{z z} & {T T}_{1313} {a a}_{x x} & {T T}_{1313} {a a}_{y the y} \\ {T T}_{21 twenty one} & {T T}_{22 twenty two} & {T T}_{23 twenty three} & {T T}_{21 twenty one} {a a}_{x x} & {T T}_{22 twenty two} {a a}_{y the y} & {T T}_{23 twenty three} {a a}_{z z} & {T T}_{21 twenty one} {a a}_{y the y} & {T T}_{21 twenty one} {a a}_{z z} & {T T}_{22 twenty two} {a a}_{x x} & {T T}_{22 twenty two} {a a}_{z z} & {T T}_{23 twenty three} {a a}_{x x} & {T T}_{23 twenty three} {a a}_{y the y} \\ {T T}_{3131} & {T T}_{3232} & {T T}_{3333} & {T T}_{3131} {a a}_{x x} & {T T}_{3232} {a a}_{y the y} & {T T}_{3333} {a a}_{z z} & {T T}_{3131} {a a}_{y the y} & {T T}_{3131} {a a}_{z z} & {T T}_{3232} {a a}_{x x} & {T T}_{3232} {a a}_{z z} & {T T}_{3333} {a a}_{x x} & {T T}_{3333} {a a}_{y the y} \end{matrix}] - - - - - - ((1818))$

${D D.}_{33 \times \times 1212} = = [\begin{matrix} {T T}_{1111} & {T T}_{1212} & {T T}_{1313} & {T T}_{1111} {ω ω}_{x x} & {T T}_{1212} {ω ω}_{y the y} & {T T}_{1313} {ω ω}_{z z} & {T T}_{1111} {ω ω}_{y the y} & {T T}_{1111} {ω ω}_{z z} & {T T}_{1212} {ω ω}_{x x} & {T T}_{1212} {ω ω}_{z z} & {T T}_{1313} {ω ω}_{x x} & {T T}_{1313} {ω ω}_{y the y} \\ {T T}_{21 twenty one} & {T T}_{22 twenty two} & {T T}_{23 twenty three} & {T T}_{21 twenty one} {ω ω}_{x x} & {T T}_{22 twenty two} {ω ω}_{y the y} & {T T}_{23 twenty three} {ω ω}_{z z} & {T T}_{21 twenty one} {ω ω}_{y the y} & {T T}_{21 twenty one} {ω ω}_{z z} & {T T}_{22 twenty two} {ω ω}_{x x} & {T T}_{22 twenty two} {ω ω}_{z z} & {T T}_{23 twenty three} {ω ω}_{x x} & {T T}_{23 twenty three} {ω ω}_{y the y} \\ {T T}_{3131} & {T T}_{3232} & {T T}_{3333} & {T T}_{3131} {ω ω}_{x x} & {T T}_{3232} {ω ω}_{y the y} & {T T}_{3333} {ω ω}_{z z} & {T T}_{3131} {ω ω}_{y the y} & {T T}_{3131} {ω ω}_{z z} & {T T}_{3232} {ω ω}_{x x} & {T T}_{3232} {ω ω}_{z z} & {T T}_{3333} {ω ω}_{x x} & {T T}_{3333} {ω ω}_{y the y} \end{matrix}] - - - - - - ((1919))$

(4)选取δV_E、δV_N、δV_U作为观测量Z，建立带杆臂效应补偿的Kalman滤波量测方程。(4) Select δV _E , δV _N , and δV _U as observations Z, and establish a Kalman filter measurement equation with lever-arm effect compensation.

在实际的标定过程中，由于激光陀螺IMU与转台的中心不完全重合，即存在杆臂。在转台旋转的时候，激光陀螺IMU会产生一个附加的线速度。因此，需要将激光陀螺IMU在导航坐标系下因杆臂效应附加的线速度进行补偿。ω_b为载体坐标系下的旋转角速度，r_b是IMU敏感中心到转台旋转中心的位置矢量，是从载体坐标系到导航坐标系的转移矩阵。In the actual calibration process, since the center of the laser gyro IMU and the turntable do not completely coincide, there is a lever arm. As the turntable rotates, the laser gyro IMU generates an additional linear velocity. Therefore, it is necessary to compensate the additional linear velocity of the laser gyro IMU due to the lever-arm effect in the navigation coordinate system. ω _b is the rotational angular velocity in the carrier coordinate system, r _b is the position vector from the sensitive center of the IMU to the rotational center of the turntable, is the transfer matrix from the carrier coordinate system to the navigation coordinate system.

r_b＝[r_xr_yr_z]^Tω_b＝[ω_xω_yω_z]^T(20)r _b ＝[r _x r _y r _z ] ^T ω _b ＝[ω _x ω _y ω _z ] ^T (20)

由杆臂效应引起的速度误差在导航坐标系下可表示为：The velocity error caused by the lever-arm effect can be expressed in the navigation coordinate system as:

${dV dV}_{l l}^{n no} = = [\begin{matrix} d d {V V}_{l l E E.}^{n no} \\ {dV dV}_{l l N N}^{n no} \\ {dV dV}_{l l U u}^{n no} \end{matrix}] = = {C C}_{b b}^{n no} (({ω ω}_{b b} \times \times {r r}_{b b})) - - - - - - ((21 twenty one))$

系统的量测方程:The measurement equation of the system:

Z＝HX+η(22)Z=HX+η(22)

其中，Z为观测矢量，H为观测矩阵，η为量测噪声。Among them, Z is the observation vector, H is the observation matrix, and η is the measurement noise.

$Z Z = = [\begin{matrix} δ δ {V V}_{E E.} \\ {δV δV}_{N N} \\ {δV δV}_{U u} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} \\ {V V}_{N N} \\ {V V}_{U u} \end{matrix}] - - [\begin{matrix} d d {V V}_{l l E E.}^{n no} \\ {dV dV}_{l l N N}^{n no} \\ {dV dV}_{l l U u}^{n no} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} - - d d {V V}_{l l E E.}^{n no} \\ {V V}_{N N} - - {dV dV}_{l l N N}^{n no} \\ {V V}_{U u} - - {dV dV}_{l l U u}^{n no} \end{matrix}] - - - - - - ((23 twenty three))$

H＝[I_3×30_3×27](24)H＝[I _3×3 0 _3×27 ](24)

η＝[η_Eη_Nη_U]^T(25)η=[η _E η _N η _U ] ^T (25)

(5)用Kalman滤波估计随机误差模型参数，并进行修正标定系数。(5) Estimate the random error model parameters with Kalman filter, and correct the calibration coefficients.

如图3所示，首先用当前24个标定系数补偿IMU数据，并带入捷联惯导解算求解出IMU在东、北、天向的速度；其次利用步骤(4)中的杆臂效应补偿方法补偿IMU东、北、天向的速度；最后采用Kalman滤波来估计步骤(3)中的24个随机误差模型参数，并与步骤(2)中的24个标定系数对应相加进行修正。As shown in Figure 3, first use the current 24 calibration coefficients to compensate the IMU data, and bring it into the strapdown inertial navigation solution to solve the speed of the IMU in the east, north, and sky directions; secondly, use the lever-arm effect in step (4) The compensation method compensates the speed of the IMU in the east, north, and sky directions; finally, the Kalman filter is used to estimate the 24 random error model parameters in step (3), and they are added to the 24 calibration coefficients in step (2) for correction.

Kalman滤波基本算法编排如图4所示:The basic algorithm arrangement of Kalman filter is shown in Figure 4:

状态一步预测方程:State one-step prediction equation:

${\overset{Λ Λ}{X x}}_{k k / / k k - - 11} = = {φ φ}_{k k,, k k - - 11} {\overset{Λ Λ}{X x}}_{k k - - 11} - - - - - - ((2626))$

状态估值计算方程:State valuation calculation equation:

${\overset{Λ Λ}{X x}}_{k k} = = {\overset{Λ Λ}{X x}}_{k k / / k k - - 11} + + {K K}_{k k} (({Z Z}_{k k} - - {H h}_{k k} {\overset{Λ Λ}{X x}}_{k k / / k k - - 11})) - - - - - - ((2727))$

滤波增量方程:Filter increment equation:

${\overset{Λ Λ}{X x}}_{k k} = = {\overset{Λ Λ}{P P}}_{k k / / k k - - 11} {H h}_{k k}^{T T} {(({H h}_{k k} {P P}_{k k / / k k - - 11} {H h}_{k k}^{T T} + + {R R}_{k k}))}^{- - 11} - - - - - - ((2828))$

一步预测均方误差方程:One-step prediction mean square error equation:

${\overset{Λ Λ}{P P}}_{k k / / k k - - 11} = = {φ φ}_{k k,, k k - - 11} {P P}_{k k - - 11} {φ φ}_{k k,, k k - - 11}^{T T} + + {Γ Γ}_{k k - - 11} {Q Q}_{k k - - 11} {Γ Γ}_{k k - - 11}^{T T} - - - - - - ((2929))$

估计均方误差方程:Estimate the mean square error equation:

${\overset{Λ Λ}{P P}}_{k k} = = ((I I - - {K K}_{k k} {H h}_{k k})) {P P}_{k k / / k k - - 11} {((I I - - {K K}_{k k} {H h}_{k k}))}^{T T} + + {K K}_{k k} {R R}_{k k} {K K}_{k k}^{T T} - - - - - - ((3030))$

Kalman滤波估计出的24个随机误差模型参数按以下方法修正步骤(2)中的24个标定系数：The 24 random error model parameters estimated by the Kalman filter are corrected by the following method to correct the 24 calibration coefficients in step (2):

角速度通道标定系数修正:Angular velocity channel calibration coefficient correction:

$[\begin{matrix} {\overset{&OverBar; &OverBar;}{G G}}_{b b x x} & {\overset{&OverBar; &OverBar;}{G G}}_{b b y the y} & {\overset{&OverBar; &OverBar;}{G G}}_{b b z z} \\ {\overset{&OverBar; &OverBar;}{S S}}_{x x} & {\overset{&OverBar; &OverBar;}{E E.}}_{y the y x x} & {\overset{&OverBar; &OverBar;}{E E.}}_{z z x x} \\ {\overset{&OverBar; &OverBar;}{E E.}}_{x x y the y} & {\overset{&OverBar; &OverBar;}{S S}}_{y the y} & {\overset{&OverBar; &OverBar;}{E E.}}_{z z y the y} \\ {\overset{&OverBar; &OverBar;}{E E.}}_{x x z z} & {\overset{&OverBar; &OverBar;}{E E.}}_{y the y z z} & {\overset{&OverBar; &OverBar;}{S S}}_{z z} \end{matrix}] = = [\begin{matrix} {G G}_{b b x x} & {G G}_{b b y the y} & {G G}_{b b z z} \\ {S S}_{x x} & {E E.}_{y the y x x} & {E E.}_{z z x x} \\ {E E.}_{x x y the y} & {S S}_{y the y} & {E E.}_{z z y the y} \\ {E E.}_{x x z z} & {E E.}_{y the y z z} & {S S}_{z z} \end{matrix}] + + [\begin{matrix} {ϵ ϵ}_{g g x x} & {ϵ ϵ}_{g g y the y} & {ϵ ϵ}_{g g z z} \\ {δS δS}_{x x} & {δE δE}_{y the y x x} & {δE δE}_{z z x x} \\ {δE δE}_{x x y the y} & {δS δS}_{y the y} & {δE δE}_{z z y the y} \\ {δE δE}_{x x z z} & {δE δE}_{y the y z z} & {δS δS}_{z z} \end{matrix}] - - - - - - ((3131))$

其中，为修正后得到的三轴陀螺零偏，为修正后的陀螺标度因数，为修正后的六个角速度通道安装误差系数。in, is the zero bias of the three-axis gyroscope obtained after correction, is the corrected gyro scaling factor, Install error coefficients for the corrected six angular velocity channels.

加速度通道标定系数修正:Acceleration channel calibration coefficient correction:

$[\begin{matrix} {\overset{&OverBar; &OverBar;}{A A}}_{b b x x} & {\overset{&OverBar; &OverBar;}{A A}}_{b b y the y} & {\overset{&OverBar; &OverBar;}{A A}}_{b b z z} \\ {\overset{&OverBar; &OverBar;}{K K}}_{x x} & {\overset{&OverBar; &OverBar;}{M m}}_{y the y x x} & {\overset{&OverBar; &OverBar;}{M m}}_{z z x x} \\ {\overset{&OverBar; &OverBar;}{M m}}_{x x y the y} & {\overset{&OverBar; &OverBar;}{K K}}_{y the y} & {\overset{&OverBar; &OverBar;}{M m}}_{z z y the y} \\ {\overset{&OverBar; &OverBar;}{M m}}_{x x z z} & {\overset{&OverBar; &OverBar;}{M m}}_{y the y z z} & {\overset{&OverBar; &OverBar;}{K K}}_{z z} \end{matrix}] = = [\begin{matrix} {A A}_{b b x x} & {A A}_{b b y the y} & {A A}_{b b z z} \\ {K K}_{x x} & {M m}_{y the y x x} & {M m}_{z z x x} \\ {M m}_{x x y the y} & {K K}_{y the y} & {M m}_{z z y the y} \\ {M m}_{x x z z} & {M m}_{y the y z z} & {K K}_{z z} \end{matrix}] + + [\begin{matrix} {ϵ ϵ}_{a a x x} & {ϵ ϵ}_{a a y the y} & {ϵ ϵ}_{a a z z} \\ {δK δK}_{x x} & {δM δ M}_{y the y x x} & {δM δ M}_{z z x x} \\ {δM δ M}_{x x y the y} & {δK δK}_{y the y} & {δM δ M}_{z z y the y} \\ {δM δ M}_{x x z z} & {δM δ M}_{y the y z z} & {δK δK}_{z z} \end{matrix}] - - - - - - ((3232))$

其中，为修正后得到的三轴加速度计零偏，为修正后的加速度计标度因数，为修正后的六个加速度通道安装误差系数。in, is the zero bias of the three-axis accelerometer obtained after correction, is the corrected accelerometer scale factor, Install error coefficients for the corrected six acceleration channels.

重复步骤(5)，直到所有24个标定系数都收敛为止，收敛条件为连续十分钟内24个标定系数在修正前后相对变化小于0.05％。此时得到精确的标定系数，用于补偿激光陀螺IMU数据。Repeat step (5) until all 24 calibration coefficients converge, and the convergence condition is that the relative changes of the 24 calibration coefficients before and after correction within ten consecutive minutes are less than 0.05%. At this time, accurate calibration coefficients are obtained, which are used to compensate the laser gyro IMU data.

Claims

1. a laser gyroscope IMU calibration method combining discrete analysis and Kalman filter, it is characterized in that comprising the following steps:

Step (1), according to the characteristics of the laser gyro and accelerometer, establish a deterministic error model of the laser gyro IMU angular velocity and acceleration channel, the model parameters include: x, y, z three-axis gyro input angular velocity ω _x , ω _y , ω _z , three-axis gyro output angular velocity G _x , G _y , G _z , three-axis gyro zero bias G _bx , G _by , G _bz , three-axis gyro scale factor S _x , S _y , S _z , three-axis gyro random error δG _x , δG _y , δG _z , six installation error coefficients E _xy , E _xz , E _yx , E _yz , E _zx , E _zy of the angular velocity channel; input acceleration a _x , a of the x, y, z three-axis accelerometer _y , a _z , three-axis accelerometer output acceleration A _x , A _y , A _z , three-axis accelerometer zero bias A _bx , A _by , A _bz , three-axis accelerometer scale factor K _x , K _y , K _z , the random errors of the three-axis accelerometer δA _x , δA _y , δA _z and the six installation error coefficients M _xy , M _xz , M _yx , M _yz , M _zx , M _zy of the acceleration channels;

Step (2), according to the laser gyro IMU angular velocity and acceleration channel deterministic error model established in step (1), using the four-axis rotation rate discrete analysis method, preliminarily calibrate the three-axis gyro zero bias G _bx , G _by , G _bz , three-axis gyro scale factors S _x , S _y , S _z , six installation error coefficients E _xy , E _xz , E _yx , E _yz , E _zx , E _zy of the angular velocity channel, three-axis accelerometer zero bias A _bx _. _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} _{_} Calibration factor;

Step (3), establish the random error model of the laser gyro IMU angular velocity and acceleration channel, the model parameters include: three-axis gyro random error δG _x , δG _y , δG _z , three-axis gyro random walk ε _x , ε _y , ε _z , three-axis gyroscope random walk driving noise Three-axis gyro scale factor error δS _x , δS _y , δS _z , six installation error coefficient deviations of angular velocity channel δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy ; random error of three-axis accelerometer δA _x , δA _y , δA _z , three-axis accelerometer random walk Triaxial Accelerometer Random Walk Drive Noise Three-axis accelerometer scale factor error δK _x , δK _y , δK _z , acceleration channel installation error coefficient deviation δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy , calculated according to random error model and inertial navigation The basic principle of Kalman filter to establish the state equation;

Step (4), select the velocity errors δV _E , δV _N , and δV _U in the east, north, and sky directions as observations, and establish a Kalman filter measurement equation with lever-arm effect compensation;

Step (5), use Kalman filter to estimate 24 random error model parameters in step (3) ε _x , ε _y , ε _z , δS _x , δS _y , δS _z , δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy , δK _x , δK _y , δK _z , δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy , and corrected by adding correspondingly to the 24 calibration coefficients in step (2), repeat this step until Until all the calibration coefficients converge, the convergence condition is that the relative change of 24 calibration coefficients before and after correction within ten minutes is less than 0.05%; in the Kalman filter measurement equation with lever arm effect compensation, the measurement Z is the Velocity errors δV _E , δV _N , and δV _U of the laser gyro IMU in the east, north, and sky directions after arm effect compensation;

Z Z = = [\begin{matrix} δ δ {V V}_{E E.} \\ {δV δV}_{N N} \\ {δV δV}_{U u} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} \\ {V V}_{N N} \\ {V V}_{U u} \end{matrix}] - - [\begin{matrix} d d {V V}_{l l E E.}^{n no} \\ {dV dV}_{l l N N}^{n no} \\ {dV dV}_{l l U u}^{n no} \end{matrix}] = = [\begin{matrix} {V V}_{E E.} - - d d {V V}_{l l E E.}^{n no} \\ {V V}_{N N} - - {dV dV}_{l l N N}^{n no} \\ {V V}_{U u} - - {dV dV}_{l l U u}^{n no} \end{matrix}]

Among them, V _E , V _N , and V _U are the speeds of the laser gyro IMU in the east, north, and sky directions, respectively, are the additional velocity errors of the lever-arm effect in the east, north, and sky directions, respectively.

2. the laser gyroscope IMU calibration method that discrete analysis according to claim 1 combines with Kalman filter, it is characterized in that: in step (2) described four azimuth rotation rate discrete analytic methods, four azimuths are laser gyroscope IMU respectively Adjust the x, y, z, and -x axes to point to the sky, and the other two axes are horizontal; when the x, y, and z axes point to the sky, the outer frame axis of the turntable will rotate at a rate between 5°/s-100°/s in the counterclockwise direction. Rotate at any constant angular rate between 2 rounds at a constant speed; when -x refers to the sky, the outer frame axis of the turntable rotates clockwise at any constant angular speed between 5°/s and 100°/s for two rounds at a constant speed.

3. the laser gyroscope IMU calibration method that discrete analysis according to claim 1 combines with Kalman filter is characterized in that: in the described Kalman filter state equation of step (3), state variable has 30 dimensions, comprises east, north, Velocity errors δV _E , δV _N , δV _U , pitch, roll, and heading angle errors φ _E , φ _N , φ _U , and three-axis gyro random walk ε _x , ε _y , ε _z , three-axis gyro Scale factor error δS _x , δS _y , δS _z , angular velocity channel six installation error coefficient deviation δE _xy , δE _xz , δE _yx , δE _yz , δE _zx , δE _zy , random walk of three-axis accelerometer Three-axis accelerometer scale factor error δK _x , δK _y , δK _z , acceleration channel installation error coefficient deviation δM _xy , δM _xz , δM _yx , δM _yz , δM _zx , δM _zy .