All-parameter calibration method suitable for MEMS-IMU array temperature compensation

Abstract

The invention provides a full-parameter calibration method suitable for MEMS-IMU array temperature compensation, and the method comprises the steps: building an inertial measurement reference through an inner and outer frame rotating shaft, and enabling each MEMS-IMU in a full-temperature range to have a unified measurement scale; secondly, a more comprehensive error model is established for the low-precision MEMS-IMU, besides conventional gyroscope / accelerometer scale factors, mounting angles and zero parameters, an error model considering gyroscope / accelerometer scale factor asymmetry, accelerometer lever arm and time asynchronous errors is established, and the error model is established; the effectiveness of the temperature compensation work under the dynamic condition and the comprehensive navigation performance of the MEMS-IMU array system are ensured; and finally, zero angular velocity measurement is introduced, angular velocities in three directions of an inertial measurement coordinate system are all zero under a static condition, and when the inner shaft or the outer shaft rotates, the angular velocity in the vertical direction is zero, so that higher stability of MEMS-IMU output at different temperatures is ensured, and when a certain IMU in the array fails, the stability of the MEMS-IMU output is ensured. And the whole inertial measurement system works stably and reliably.

Description

Technical Field

[0001] The invention belongs to the technical field of MEMS inertial navigation, and in particular relates to a full-parameter calibration method suitable for MEMS-IMU array temperature compensation. Background Art

[0002] In order to adapt to high and low temperature environments and ensure the working accuracy of the system, the inertial navigation system needs to undergo a temperature compensation test before leaving the factory to compensate for the error of the MEMS-IMU with temperature changes. The simplest temperature compensation method for the inertial navigation system is to collect static data of the inertial instrument at different temperature points to establish the relationship between the output of the inertial instrument and the temperature, and then compensate according to the real-time temperature. However, this method mixes the temperature drift errors of the zero bias, scale factor and installation angle into one, and its temperature compensation effect and system accuracy are often unsatisfactory under dynamic conditions. Subsequent inertial technology engineers have developed a temperature compensation method that performs system-level calibration at various temperature points and then performs temperature curve fitting on all key parameters to ensure dynamic navigation performance under different ambient temperatures. Generally, this system-level calibration method establishes an inertial measurement benchmark with a gyroscope sensitive axis with better full-temperature stability, thereby completing the compensation of the installation angle temperature drift error. However, for low-cost MEMS-IMUs, the sensitive axis of its gyroscope will change significantly from low temperature to high temperature, and it is no longer suitable as an inertial measurement benchmark that requires high stability.

[0003] On the other hand, the array-type MEMS inertial navigation system improves its performance by fusing multiple instrument signals on the same axis. If each IMU establishes its own inertial measurement reference with its own inertial sensor sensitive axis, there will be differences in installation deflection angle and scale factor between each inertial sensor, and there will also be installation lever arm deviations between MEMS accelerometers. If the MEMS inertial sensor is still modeled as sensitive to the same true angular rate, data fusion errors are bound to occur. Of course, directly taking the arithmetic average of the raw data of the inertial sensor to virtualize it as an inertial measurement unit and then performing system-level calibration can also avoid the occurrence of this problem, but when an IMU in the array fails and needs to be isolated, this method may cause the entire inertial measurement system to fail. Therefore, it is particularly necessary to take a method to unify the measurement scale of each IMU in the array. Summary of the invention

[0004] To solve the above problems, the present invention provides a full-parameter calibration method suitable for temperature compensation of MEMS-IMU arrays, which optimizes the traditional method of establishing a coordinate system reference based on the sensitive axis of the inertial sensor to establishing an inertial measurement reference based on the rotation axes of the inner and outer frames, so that all MEMS-IMUs have a unified measurement scale within the full temperature range, which not only ensures that the MEMS-IMU output has higher stability at different temperatures, but also ensures that the entire inertial measurement system works stably and reliably when a certain IMU in the array fails.

[0005] A full parameter calibration method suitable for temperature compensation of a MEMS-IMU array, wherein the MEMS-IMU array is placed inside a dual-axis turntable for calibration, and the method comprises the following steps:

[0006] S1: Construct an inertial measurement coordinate system p, where z of the inertial measurement coordinate system p is p The y axis is the inner frame rotation axis of the dual-axis turntable. p The axis is located in the plane formed by the inner frame rotation axis and the outer frame rotation axis and is perpendicular to the inner frame rotation axis. p With y p 、z p Form a right-handed rectangular coordinate system;

[0007] S2: Construct accelerometer error model and gyroscope error model in inertial measurement coordinate system p;

[0008] S3: construct the state equation and measurement equation of the Kalman filter according to the accelerometer error model and gyroscope error model in the inertial measurement coordinate system p and the navigation system attitude error angle and navigation coordinate system velocity error in the navigation coordinate system n;

[0009] S4: Iteratively solve the state equation and measurement equation of the Kalman filter to obtain the accelerometer error estimate, the gyroscope error estimate, the navigation system attitude error angle estimate, and the navigation coordinate system velocity error estimate;

[0010] S5: Correct the corresponding errors according to the four types of estimated values ​​to complete the calibration of all parameters.

[0011] Furthermore, the accelerometer error model corresponding to any IMU in the MEMS-IMU array is constructed in the inertial measurement coordinate system p as follows:

[0012] δf p =H A X A

[0013] Among them, δf p is the accelerometer specific force measurement error, X A is the accelerometer error state, H A is the observation matrix of the accelerometer error state, and the accelerometer error state X A for:

[0014]

[0015] Where T represents transpose, δK a (:,i), i=1,2,3 represents δK a The i-th column of aDenotes the accelerometer scale factor matrix K a The estimation error of diag(·) represents the construction of a three-dimensional diagonal matrix from the three-dimensional vectors in brackets. The accelerometer scale factor vector K a1 =[k a1x k a1y k a1z ] T , k a1x , k a1y and k a1z They represent the scale factors of the three-axis accelerometers in the IMU, δK is the transformation matrix from the oblique coordinate system formed by the three accelerometer sensitive axes in the IMU to the inertial measurement coordinate system p; as represents the estimated error of the accelerometer asymmetry; δL a Indicates the inner arm error of the inner frame of the dual-axis turntable; represents the accelerometer zero position estimation error equivalently converted to the inertial measurement coordinate system p; δt GA Indicates the measurement delay error of the accelerometer relative to the gyroscope;

[0016] The observation matrix H of the accelerometer error state A for:

[0017]

[0018] Among them, N a =[N ax N ay N az ] T Represents the raw output pulse of the three-axis accelerometer in the IMU; M a represents the accelerometer asymmetric error observation matrix; R represents the transformation matrix that transforms the inertial measurement angular velocity and angular acceleration from the oblique coordinate system to the inertial measurement coordinate system p; I 3×3 represents the three-dimensional identity matrix; represents the projection of the rotation angular rate of the inertial measurement coordinate system p relative to the navigation coordinate system n in the inertial measurement coordinate system p; f p represents the specific force of the inertial measurement frame p.

[0019] Furthermore, the gyro error model corresponding to any IMU in the MEMS-IMU array is constructed in the inertial measurement coordinate system p as follows:

[0020] δω p =H G X G

[0021] Among them, δω p Indicates the gyro output angular velocity measurement error, XG Indicates the gyro error state, H G The observation matrix represents the gyro error state, and the gyro error state X G for:

[0022]

[0023] Where T represents transpose, δK g (:,1),i=1,2,3 represents δK g The i-th column of g Represents the gyro scale factor matrix K g The estimation error of diag(·) indicates that a three-dimensional diagonal matrix is ​​constructed from the three-dimensional vectors in the brackets. The gyro scale factor vector K g1 =[k g1x k g1y k g1z ] T , k g1x , k g1y and k g1z Respectively represent the scale factors of the three-axis gyroscope in the IMU, δK is the transformation matrix from the oblique coordinate system formed by the three gyro sensitive axes in the IMU to the inertial measurement coordinate system p; gs represents the estimated error of the gyro asymmetry; δε p represents the gyro zero position estimation error equivalently converted to the inertial measurement coordinate system p;

[0024] The observation matrix H of the gyro error state G for:

[0025] H G =[diag(N gx ) diag(N gy ) diag(N gz ) M g -I]

[0026] Among them, N g =[N gx N gy N gz ] T Represents the original output pulse of the three-axis gyroscope in the IMU; M g represents the gyro asymmetry error observation matrix; I represents the unit matrix.

[0027] Furthermore, the state equation and measurement equation of the Kalman filter are as follows:

[0028]

[0029] ZCali =H Cali X Cali

[0030] Among them, X Cali is the state quantity of the Kalman filter, Represents X Cali The differential of Φ Cali is the system matrix of the Kalman filter, W is the system noise matrix; Z Cali is the measurement of the Kalman filter, H Cali is the measurement matrix of the Kalman filter; where the state quantity X Cali According to the accelerometer error model and gyroscope error model in the inertial measurement coordinate system p, and the navigation system attitude error angle and navigation coordinate system velocity error in the navigation coordinate system n, the details are as follows:

[0031] X Cali =[(φ n ) T (δv n ) T (X G ) T (X A ) T ] T

[0032] Among them, φ n is the navigation system attitude error angle in navigation coordinate system n, δv n is the velocity error in the navigation coordinate system n, X A is the accelerometer error state in the inertial measurement coordinate system p, X G is the gyro error state in the inertial measurement coordinate system p, and T represents the transposition.

[0033] Furthermore, the Kalman filter measurement Z Cali and the measurement matrix Η Cali It is expressed as follows:

[0034]

[0035] in, express The error, represents the projection of the rotational angular velocity of the inertial measurement coordinate system p relative to the body coordinate system b in the inertial measurement coordinate system p; 3×3 , O 3×15 , O 3×25 Respectively represent matrices of order 3×3, 3×15, and 3×25 whose elements are all zero; I 3×3 represents the three-dimensional identity matrix; A G represents the angular velocity measurement coefficient matrix; H GThe observation matrix representing the gyro error state; Represents the transformation matrix from the navigation coordinate system n to the inertial measurement coordinate system p; represents the projection of the Earth's rotation angular velocity in the navigation coordinate system n; (·×) represents the antisymmetric operation operator.

[0036] Furthermore, when the dual-axis turntable is not rotating, So there is A G =I 3×3 ; When the outer frame axis rotates, the x perpendicular to the inner and outer frame axis plane p The angular rate of the axis sensitivity is 0, at this time A G =diag([1 0 0]); when the inner frame rotates, the x value perpendicular to the inner frame axis p Axis and y p The angular rates of the axis sensitivity are all 0. At this time, A G =diag(

[110] ).

[0037] Beneficial effects:

[0038] The present invention provides a full-parameter calibration method suitable for temperature compensation of MEMS-IMU arrays. Firstly, a stable and unified inertial measurement reference is established. Considering the non-orthogonality of the rotation axes of the inner and outer frames, the traditional method of establishing the inertial measurement reference with the sensitive axis of the inertial sensor is optimized to establishing the inertial measurement reference with the rotation axes of the inner and outer frames, so that each MEMS-IMU in the full temperature range has a unified measurement scale; secondly, a more comprehensive error model is established for the low-precision MEMS-IMU, considering the maximum error model parameters that can be stimulated in a conventional gravity environment, a linearized error model is constructed for the indirect Kalman filter, and the error parameters to be estimated and the corresponding error parameters are established. The system equation ensures the effectiveness of temperature compensation under dynamic conditions and the comprehensive navigation performance of the MEMS-IMU array system. Finally, in addition to the traditional zero-speed measurement, the zero angular velocity measurement is introduced to update the Kalman filter speed / angular velocity measurement algorithm. Under static conditions, the rotation angular velocities in the three directions are all zero. When the inner or outer axis rotates, the angular velocity in the vertical direction is zero. In this way, the measurement scale and benchmark of the inertial sensors in the array are unified and the accuracy is improved, which not only ensures that the MEMS-IMU output has higher stability under different temperatures, but also ensures that when an IMU in the array fails, the entire inertial measurement system can work stably and reliably. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 A flow chart of a full parameter calibration method suitable for MEMS-IMU array temperature compensation provided by the present invention;

[0040] Figure 2 Definition of the inertial measurement coordinate system in the MEMS-IMU array provided by the present invention;

[0041] Figure 3 The present invention provides a full-parameter self-calibration calculation process for the MEMS array. DETAILED DESCRIPTION

[0042] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0043] The present invention provides a full parameter calibration method suitable for MEMS-IMU array temperature compensation, wherein the MEMS-IMU array is placed inside a dual-axis turntable for calibration, such as Figure 1 As shown, the method comprises the following steps:

[0044] S1: Construct an inertial measurement coordinate system p, where z of the inertial measurement coordinate system p is p The y axis is the inner frame rotation axis of the dual-axis turntable. p The axis is located in the plane formed by the inner frame rotation axis and the outer frame rotation axis and is perpendicular to the inner frame rotation axis. p With y p 、z p Form a right-handed rectangular coordinate system;

[0045] It should be noted that for a multi-IMU array, x i y i z i Represents the body coordinate system of the i-th IMU, and its relationship with the inertial measurement coordinate system is as follows: Figure 2 For any IMU, when the inner frame rotates, its z i The gyro on the axis is sensitive to the main rotation angular rate, so the inner frame rotation axis is used as the z p axis; at the same time, when the outer frame axis rotates, y i The gyro on the y-axis is sensitive to the main angular rate of rotation, so the y p The axis is set to be located in the plane where the inner frame axis and the outer frame axis are located and perpendicular to the inner frame axis. Finally, according to x p With y p 、z p Use a right-handed rectangular coordinate system to determine x p ; When the inner and outer frames rotate, the inertial measurement coordinate system p is fixed to each IMU and rotates with the inner and outer frames.

[0046] That is to say, when the inner frame axis rotates, if the i (i = x, y, z) gyroscope is sensitive to most of the rotation angular rates, then the i (i = x, y, z) axis of the inertial measurement reference coordinate system is established to coincide with the inner frame axis; when the outer frame axis rotates, if the j (j = x, y, z) gyroscope is sensitive to most of the rotation angular rates, then the j (j = x, y, z) axis of the inertial measurement reference coordinate system is established to be located in the plane formed by the inner frame axis and the outer frame axis and perpendicular to the inner frame axis, and the k (k = x, y, z) axis of the inertial measurement reference coordinate system and the i (i = x, y, z) axis and the j (j = x, y, z) axis form a right-handed rectangular coordinate system; it should be noted that constructing an inertial measurement coordinate system with inner and outer rotation axes is the key to the present invention. The definition of this coordinate system breaks the traditional calibration method of establishing an inertial measurement reference with the sensitive axis of an inertial instrument, unifies the measurement scales of all IMUs, and maintains high stability over the entire temperature range. In addition, the present invention also defines the body coordinate system b: that is, the inertial measurement coordinate system when the inner and outer frames are both in the electrical zero position. Navigation coordinate system n: adopts the northeast sky geographic coordinate system, x n and n The axis lies in the local horizontal plane, and x n Point to the east, y n Point to the north, z n Point upward along the vertical line of the ground.

[0047] S2: Construct accelerometer error model and gyroscope error model in inertial measurement coordinate system p;

[0048] S3: construct the state equation and measurement equation of the Kalman filter according to the accelerometer error model and gyroscope error model in the inertial measurement coordinate system p and the navigation system attitude error angle and navigation coordinate system velocity error in the navigation coordinate system n;

[0049] S4: Iteratively solve the state equation and measurement equation of the Kalman filter to obtain the accelerometer error estimate, the gyroscope error estimate, the navigation system attitude error angle estimate, and the navigation coordinate system velocity error estimate;

[0050] S5: Correct the corresponding errors according to the four types of estimated values ​​to complete the calibration of all parameters, that is, use the accelerometer error estimate, gyroscope error estimate, navigation system attitude error angle estimate, and navigation coordinate system velocity error estimate to correct the angular acceleration measured by the accelerometer, the angular velocity measured by the gyroscope, the attitude angle under the navigation system, and the velocity under the navigation coordinate system, respectively, to complete the calibration of the angular acceleration, angular velocity, attitude angle under the navigation system, and velocity under the navigation coordinate system.

[0051] It should be noted that the MEMS array contains multiple IMUs, each of which contains x g ,y g ,z g Three-axis gyro and x a ,ya ,z a Three-axis accelerometer. Taking one of the IMUs as an example, the following is a detailed derivation of the process of constructing the accelerometer error model and the gyroscope error model in the inertial measurement coordinate system p.

[0052] 1.1 Accelerometer error model

[0053] Model the output specific force from the accelerometer pulse:

[0054]

[0055] Among them, f p represents the specific force of the inertial measurement coordinate system, diag(·) indicates that a three-dimensional diagonal matrix is ​​constructed from the three-dimensional vectors in brackets. K a1 =[k a1x k a1y k a1z ] T , k a1x , k a1y and k a1z Represent the scale factors of the three accelerometers, K is the transformation matrix from the oblique coordinate system formed by the three accelerometer sensitive axes to the inertial measurement coordinate system; as =[k sax k say k saz ] T , k sax , k say and k saz Respectively represent the asymmetry of the scale factors of the three accelerometers; N a =[N ax N ay N az ] T Represents the original output pulse of the triaxial accelerometer; S a =[sign(N ax ) sign(N ay ) sign(N az )] T , sign(·) means taking the sign of the content in the brackets; represents the lever arm compensation term in the accelerometer; Represents the accelerometer zero position equivalently converted to the inertial measurement p system; Compensates for the accelerometer measurement delay.

[0056] According to the transformation matrix Its own meaning can also be expressed as:

[0057]

[0058] Among them, μ i,j , i = x, y, z; j = x, y, z represents the sensitive axis of the accelerometer i and the inertial measurement coordinate system j p The angle between the axes.

[0059] In order to facilitate the subsequent error equation expression, K a The specific elements in are shown in formula 3

[0060]

[0061] Taking the intersection of the inner and outer frame rotation axes as the origin of the P system coordinates, the inner rod arm compensation term It can be expressed as:

[0062]

[0063] in

[0064]

[0065] in, Express Perform antisymmetric operations;

[0066] j=1,2,3 means C a p The jth column of aij , i = x, y, z, j = x, y, z means the distance between the sensitive center of accelerometer i and the origin of the inertial measurement coordinate system is j p Projection on the axis, ω p and are the inertial measurement angular velocity and angular acceleration respectively. The previous sampling period t is used in formula (5) k-1 The inertial measurement angle increment Δθ at time p (t k-1 ) and the current t k The inertial measurement angle increment Δθ at time p (t k ) to perform real-time calculations.

[0067]

[0068] Among them, T s Indicates the sampling period;

[0069] Gyro / timer asynchronism compensation It can be expressed as

[0070]

[0071] is the projection of the rotation angular rate of the inertial measurement coordinate system relative to the navigation coordinate system in the inertial measurement coordinate system, t GA is the 1D measurement delay of the accelerometer relative to the gyroscope to be estimated.

[0072] Due to the error in the accelerometer parameters, its actual output specific force is

[0073]

[0074] In formula (8), each parameter to be calibrated is added with a superscript ~, indicating that the current estimated value of the parameter contains errors, and its specific meaning is consistent with the corresponding parameter in formula (1). At this time, the accelerometer specific force measurement error δf p It can be expressed as:

[0075]

[0076] where δK a K a The estimation error of δK as represents the estimated error of the accelerometer asymmetry, δL a Indicates the inner lever arm error; represents the accelerometer zero position estimation error; Indicates the accelerometer delay compensation error; M a Represents the accelerometer asymmetric error observation matrix, which can be specifically expressed as

[0077]

[0078] X A The accelerometer error state obtained after sorting can be expressed as

[0079]

[0080] where δK a (:,i), i=1,2,3 represents δK a The i-th column of GA Represents the measurement delay error of the accelerometer relative to the gyroscope.

[0081] H A The observation matrix representing the accelerometer error state is

[0082]

[0083] Where diag(*) represents a diagonal matrix whose diagonal elements are all *.

[0084] 2.2 Gyro Error Model

[0085] Similarly, starting from the gyro pulse, the gyro output angular velocity model is established, so we have

[0086] ω p =K g [I 3×3 +diag(K gs )diag(S)]N g -ε p (13)

[0087] where ω p represents the gyro output angular velocity of the inertial measurement coordinate system, K g1 =[k g1x k g1y k g1z ] T , k g1x , k g1y and k g1z Represent the scale factors of the three gyroscopes, K is the transformation matrix from the oblique coordinate system formed by the three gyro sensitive axes in the IMU to the inertial measurement coordinate system; gs =[k sgx k sgy k sgz ] T , k sgx , k sgy and k sgz Respectively represent the asymmetry of the three gyro scale factors; N g =[N gx N gy N gz ] T Indicates the gyro output pulse, ε p is the 3D gyro zero position equivalently converted to the p system; S g =[sign(N gx ) sign(N gy ) sign(N gz )] T .

[0088] Its meaning is shown in formula 14

[0089]

[0090] where η i,j , i = x, y, z; j = x, y, z represents the relationship between the i gyro sensitive axis and the j inertial measurement coordinate system p The angle between the axes.

[0091] In order to facilitate the subsequent error equation expression, K g The specific elements in are shown in formula 15.

[0092]

[0093] Due to the error in the parameters to be calibrated, the actual output angular velocity is

[0094]

[0095] In formula (16), each parameter to be calibrated is added with a superscript ~, indicating that the current estimated value of the parameter contains errors, and its connotation is consistent with the corresponding parameter in formula (13).

[0096] Therefore, the gyro output angular velocity error can be expressed as:

[0097]

[0098] δω p Denotes the gyro output angular velocity measurement error, δK g K g The estimation error, δK gs represents the estimated error of the asymmetry of the gyro scale factor, δε p Indicates the gyro zero position error equivalently converted to the p system. g It represents the asymmetric error measurement matrix of the gyro scale factor, which can be specifically expressed as

[0099]

[0100] X G Indicates the error state of the gyro, H G The observation matrix representing the gyro error state is

[0101]

[0102] H G =[diag(N gx )diag(N gy )diag(N gz )M g -I] (20)

[0103] The Kalman filter system equation can be established based on the strapdown inertial navigation error equation and the gyro error model shown in equations (11) and (17). The accurate calibration of the inertial measurement error parameters can be achieved by using a multiple-iteration Kalman filter estimation algorithm.

[0104] The iterative Kalman filter calibration process is introduced in detail below.

[0105] The calculation of inertial measurement error parameters adopts the "inertial system anti-disturbance self-alignment + backtracking Kalman filter" algorithm with multiple iterations. In the process of error coefficient estimation, according to the general 19-position flip scheme, the static and dynamic flip data of each position are continuously sampled, and the speed error and attitude error formula are combined to iteratively calculate the error parameters. Figure 3 shown.

[0106] Figure 3 The multi-position rolling sequence adopts the universal 19-position calibration sequence, as shown in Table 1.

[0107] Table 1 General 19-position rollover order

[0108]

[0109] The Kalman filter equation can be specifically expressed as:

[0110]

[0111] X Cali is the state quantity of the Kalman filter, Represents X Cali The differential of Φ Cali is the system matrix of the Kalman filter, W is the system noise matrix; Z Cali is the measurement of the Kalman filter, H Cali is the measurement matrix of the Kalman filter.

[0112] X Cali Specifically, it can be expressed as

[0113] X Cali =](φ n ) T (δv n ) T (X G ) T (X A ) T ] T (twenty two)

[0114] In the above formula, φ n is the navigation system attitude error angle vector, δv n Represents the velocity error vector of the navigation coordinate system.

[0115] System matrix Φ Cali It can be obtained by combining the strapdown inertial navigation attitude error differential equation and velocity error differential equation with the accelerometer error model and gyroscope error model shown in equations 9 and 17, as shown in equation 23.

[0116]

[0117] The system noise matrix W~N(0,Q), where Q is the system noise variance matrix.

[0118] O 3×37 Represents a 3×37 matrix whose elements are all zero.

[0119] Projection of the rotational angular velocity of the inertial measurement coordinate system relative to the body coordinate system in the inertial measurement coordinate system It can be expressed as:

[0120]

[0121] in represents the projection of the gyro output angular rate in the inertial measurement coordinate system, Represents the transformation matrix from the navigation coordinate system to the inertial measurement coordinate system, Represents the projection of the Earth's rotation angular velocity in the navigation system.

[0122] Due to the error,

[0123]

[0124] Indicates error Indicates error Indicates error

[0125] Subtracting equation 24 from equation 25, we get

[0126]

[0127] Therefore, the measurement equation in Equation 21 can be expressed as

[0128]

[0129] When the indexing mechanism does not rotate, Therefore, there is Angular velocity measurement coefficient matrix A G =I 3×3 ; When the outer frame axis rotates from the 1st to the 6th order in Table 1, the x perpendicular to the inner and outer frame axis plane p The angular rate of the axis sensitivity is 0, that is At this time A G =diag(

[100] ); When the inner frame rotates in the 7th order and the 14th to 18th order in Table 1, the x perpendicular to the inner frame axis p Axis and y p The angular rates of the axis sensitivity are all 0, that is, At this time A G =diag(

[110] ).

[0130] The above calculation method can be used to complete the iterative Kalman filter calibration calculation. Other IMUs in the array can use the same coordinate system definition and calculation method to complete their own full parameter calibration calculation and unification of inertial measurement benchmarks.

[0131] In summary, the present invention provides a full parameter calibration method suitable for MEMS-IMU array temperature compensation, which can be summarized into the following three steps:

[0132] The first step is to establish an inertial measurement reference coordinate system; based on the inner and outer frame rotation axes, a unified inertial measurement reference is established for all IMUs in the array, the inner frame axis is established as one axis of the inertial measurement reference coordinate system, and the axis perpendicular to the inner frame axis in the plane formed by the inner and outer frame axes is established as the other axis of the inertial measurement reference coordinate system. This takes into account the non-orthogonality of the inner and outer frame rotation axes, and can maintain the spatial stability of the inertial measurement reference in the full temperature range of -40℃ to 60℃.

[0133] The second step is to refine the full-parameter error model. In addition to the conventional gyro / accelerometer scale factor, installation angle and zero position parameters, the gyro / accelerometer scale factor asymmetry, accelerometer arm and gyro / accelerometer time synchronization errors are considered to build a linearized error model for the indirect Kalman filter. That is, a MEMS-IMU linear calibration error model considering the gyro / accelerometer scale factor asymmetry, accelerometer arm and gyro / accelerometer time synchronization errors is established, and the error parameters to be estimated and the corresponding system equations are established to ensure the effectiveness of temperature compensation under dynamic conditions and the comprehensive navigation performance of the MEMS-IMU array system.

[0134] The third step is the iterative Kalman filter measurement update calculation method. On the basis of the general 19-position system-level calibration test data, zero speed is used as the measurement under static conditions, and the angular velocity of the vertical direction of the rotation axis during rotation is zero as the measurement. The iterative Kalman filter system-level calibration is implemented to realize the calibration of three-dimensional gyro measurement and three-dimensional plus measurement to the inertial measurement coordinate system.

[0135] It can be seen that the present invention optimizes the traditional way of establishing a coordinate system reference with the sensitive axis of the inertial sensor to establish an inertial measurement reference with the rotation axis of the inner and outer frames, so that all MEMS-IMUs have a unified measurement scale within the full temperature range, which not only ensures that the MEMS-IMU output has higher stability at different temperatures, but also ensures that when a certain IMU in the array fails, the entire inertial measurement system can work stably and reliably. At the same time, the traditional MEMS-IMU calibration error model is further refined to ensure the adaptability and accuracy of temperature compensation under dynamic conditions.

[0136] Of course, the present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may certainly make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A full parameter calibration method for temperature compensation of a MEMS-IMU array, wherein the MEMS-IMU array is placed inside a dual-axis turntable for calibration, characterized in that: The method comprises the following steps: S1: Construct an inertial measurement coordinate system p, where z of the inertial measurement coordinate system p is p The y axis is the inner frame rotation axis of the dual-axis turntable. p The axis is located in the plane formed by the inner frame rotation axis and the outer frame rotation axis and is perpendicular to the inner frame rotation axis. p With y p 、z p Form a right-handed rectangular coordinate system; S2: Construct accelerometer error model and gyroscope error model in inertial measurement coordinate system p; S3: construct the state equation and measurement equation of the Kalman filter according to the accelerometer error model and gyroscope error model in the inertial measurement coordinate system p and the navigation system attitude error angle and navigation coordinate system velocity error in the navigation coordinate system n; S4: Iteratively solve the state equation and measurement equation of the Kalman filter to obtain the accelerometer error estimate, the gyroscope error estimate, the navigation system attitude error angle estimate, and the navigation coordinate system velocity error estimate; S5: Correct the corresponding errors according to the four types of estimated values ​​to complete the calibration of all parameters.

2. A full parameter calibration method suitable for MEMS-IMU array temperature compensation as claimed in claim 1, characterized in that: The accelerometer error model corresponding to any IMU in the MEMS-IMU array is constructed in the inertial measurement coordinate system p as follows: δf p =H A X A Among them, δf p is the accelerometer specific force measurement error, X A is the accelerometer error state, H A is the observation matrix of the accelerometer error state, and the accelerometer error state X A for: Where T represents transpose, δK a (:,i), i=1,2,3 represents δK a The i-th column of a Denotes the accelerometer scale factor matrix K a The estimation error of diag(·) represents the construction of a three-dimensional diagonal matrix from the three-dimensional vectors in brackets. The accelerometer scale factor vector K a1 =[k a1x k a1y k a1z ] T , k a1x , k a1y and k a1z They represent the scale factors of the three-axis accelerometers in the IMU, δK is the transformation matrix from the oblique coordinate system formed by the three accelerometer sensitive axes in the IMU to the inertial measurement coordinate system p; as represents the estimated error of the accelerometer asymmetry; δL a Indicates the inner arm error of the inner frame of the dual-axis turntable; δ▽ p represents the accelerometer zero position estimation error equivalently converted to the inertial measurement coordinate system p; δt GA Indicates the measurement delay error of the accelerometer relative to the gyroscope; The observation matrix H of the accelerometer error state A for: Among them, N a =[N ax N ay N az ] T Represents the raw output pulse of the three-axis accelerometer in the IMU; M a represents the accelerometer asymmetric error observation matrix; R represents the transformation matrix that transforms the inertial measurement angular velocity and angular acceleration from the oblique coordinate system to the inertial measurement coordinate system p; I 3×3 represents the three-dimensional identity matrix; represents the projection of the rotation angular rate of the inertial measurement coordinate system p relative to the navigation coordinate system n in the inertial measurement coordinate system p; f p represents the specific force of the inertial measurement frame p.

3. A full parameter calibration method suitable for MEMS-IMU array temperature compensation as claimed in claim 1, characterized in that: The gyro error model corresponding to any IMU in the MEMS-IMU array is constructed in the inertial measurement coordinate system p as follows: here p =H G X G Among them, δω p Indicates the gyro output angular velocity measurement error, X G Indicates the gyro error state, H G The observation matrix represents the gyro error state, and the gyro error state X G for: Where T represents transpose, δK g (:,1),i=1,2,3 represents δK g The i-th column of g Represents the gyro scale factor matrix K g The estimation error of diag(·) indicates that a three-dimensional diagonal matrix is ​​constructed from the three-dimensional vectors in the brackets. The gyro scale factor vector K g1 =[k g1x k g1y k g1z ] T , k g1x , k g1y and k g1z Respectively represent the scale factors of the three-axis gyroscope in the IMU, δK is the transformation matrix from the oblique coordinate system formed by the three gyro sensitive axes in the IMU to the inertial measurement coordinate system p; gs represents the estimated error of the gyro asymmetry; δε p represents the gyro zero position estimation error equivalently converted to the inertial measurement coordinate system p; The observation matrix H of the gyro error state G for: H G =[diag(N gx ) diag(N gy ) diag(N gz ) M g -I] Among them, N g =[N gx N gy N gz ] T Represents the original output pulse of the three-axis gyroscope in the IMU; M g represents the gyro asymmetry error observation matrix; I represents the unit matrix.

4. A full parameter calibration method suitable for MEMS-IMU array temperature compensation as claimed in claim 1, characterized in that: The state equation and measurement equation of the Kalman filter are as follows: Z Cali =Η Cali X Cali Among them, X Cali is the state quantity of the Kalman filter, Represents X Cali The differential of Φ Cali is the system matrix of the Kalman filter, W is the system noise matrix; Z Cali is the measurement of the Kalman filter, H Cali is the measurement matrix of the Kalman filter; where the state quantity X Cali According to the accelerometer error model and gyroscope error model in the inertial measurement coordinate system p, and the navigation system attitude error angle and navigation coordinate system velocity error in the navigation coordinate system n, the details are as follows: X Cali =[(φ n ) T (dv n ) T (X G ) T (X A ) T ] T Among them, φ n is the navigation system attitude error angle in navigation coordinate system n, δv n is the velocity error in the navigation coordinate system n, X A is the accelerometer error state in the inertial measurement coordinate system p, X G is the gyro error state in the inertial measurement coordinate system p, and T represents the transposition.

5. A full parameter calibration method suitable for MEMS-IMU array temperature compensation as claimed in claim 4, characterized in that: Kalman filter measurement Z Cali and the measurement matrix Η Cali It is expressed as follows: in, express The error, represents the projection of the rotational angular velocity of the inertial measurement coordinate system p relative to the body coordinate system b in the inertial measurement coordinate system p; 3×3 , O 3×15 , O 3×25 Respectively represent matrices of order 3×3, 3×15, and 3×25 whose elements are all zero; I 3×3 represents the three-dimensional identity matrix; A G represents the angular velocity measurement coefficient matrix; H G The observation matrix representing the gyro error state; Represents the transformation matrix from the navigation coordinate system n to the inertial measurement coordinate system p; represents the projection of the Earth's rotation angular velocity in the navigation coordinate system n; (·×) represents the antisymmetric operation operator.

6. A full parameter calibration method suitable for MEMS-IMU array temperature compensation as claimed in claim 5, characterized in that: When the dual-axis turntable is not rotating, So there is A G =I 3×3 ; When the outer frame axis rotates, the x perpendicular to the inner and outer frame axis plane p The angular rate of the axis sensitivity is 0, at this time A G =diag([1 0 0]); when the inner frame rotates, the x value perpendicular to the inner frame axis p Axis and y p The angular rates of the axis sensitivity are all 0. At this time, A G =diag([1 1 0]).

Images