Two-Stage Calibration Scheme for Magnetic Measurement System on Guided Munition

In order to calibrate the magnetic measurement system used in guided munition on site, a two-stage calibration (TSC) scheme without reference is proposed in this paper. Analyzing the interfering magnetic field in the projectile and misalignment angles between the projectile coordinate system and measurement coordinate system establishes a proper mathematical equivalent model and derives a calibration method. The first stage is ellipsoid fitting to obtain the equivalent zero-offset, equivalent sensitivity and equivalent non-orthogonal angles of the sensor; the second stage is to calibrate the misalignment angles between the projectile coordinate system and the measurement coordinate system with the three-position calibration (TPC) method. Complete calibration is convenient to operate and does not need an additional reference, which has wide applicability. The simulation results show that the deviation in the measured value after compensation is within 100 nT. The experiment proves that the error of compensated magnetic value is about 150 nT, which meets the accuracy of requirements in guided munitions.


Introduction
In the field of navigation technology, the triaxial magnetic measurement system, with the AMR magnetometer as the core sensor which measures the earth's magnetic field to obtain precise directional information [1], has been widely applied to vehicle detection [2], resource exploration [3], magnetic anomaly detection and navigation [4][5][6], due to the advantages of using a magnetometer such as small space occupation, fast response speed, good anti-overload performance, low cost, etc. [4,7]; as well, it is used in guided munitions. Generally, the magnetic measurement system must be calibrated in advance [8]; many methods need a high-precision turntable or an additional reference, however, in the field of guided munition, due to the long-term manufacturing, storage, transportation and assembling, the interfering magnetism of the projectile will change, meaning that the parameters are inaccurate which are calibrated before, such as zero-offset, sensitivity and non-orthogonal angles [9]. Besides, there are misalignment angles between the projectile coordinate system and measurement coordinate system [10]. Theoretically, a one-degree misalignment angle results in thousands of nT errors of measured value [11], which need to rotate the sensor coordinate system to coincide with the projectile coordinate system.
Aiming towards the calibration for the magnetometer and the magnetic measurement system, some scholars has proposed some effective methods. Valérie Renaudin proposed a complete calibration in magnetic domain in [11], they analyzed the reasons of errors and fitted the measured value of magnetometer as a sphere; however, the equation represents a second-order surface: it may be a hyperboloid, a cone surface or an ellipsoid, there is no specific restriction, and it cannot judge the validity of the calibration after the process is completed. Yanxia Liu proposed a two-stage calibration based on particle swarm optimization in [12]. Particle swarms belong to a stochastic global optimal algorithm, and can calibrate all parameters of sensors and improve the accuracy of triaxial sensor. Hard magnetism and soft magnetism were considered when they established the error model; however, some necessary components are ignored in the analysis, which cannot be estimated and extra deviations still exist, and the particle swarm optimization algorithm has its disadvantages: the more complex the surrounding environment is, the more samples are needed to describe the posterior probability distribution, and the higher the complexity of the algorithm, the less it is suitable for calibrating magnetic measurement systems on guided munitions on site. Supeng Li proposed a method to calibrate the misalignment angle in [13]; the method is based on the principle that the trajectory of the measurement data analyzes the graphical feature to obtain the constrained relationship among the components of the coordinate system. A three-step calibration method for tri-axial sensors was proposed in [14]; the method aims to calibrate the magnetometer and accelerometer simultaneously, one of the two kinds of sensors is calibrated first to define an internal reference, which is not suitable for measurement systems with one sensor. Donghui Liu proposed a simplified ellipsoid fitting to calibrate magnetometers [15]; however, they could not solve the problem of misalignment angles. Lee Jung analyzed the misalignment and soft-iron distortion to calibrate magnetometers, however this is the internal misalignment error of the sensor [16]. Tandon P. proposed a method based on ellipsoid fitting [17]; ellipsoid fitting is a classical method to calibrate field vector sensors, however, projectiles cannot be installed on the calibration equipment. The methods above can improve the accuracy of sensors or systems effectively, however, some of them need additional high-precision equipment or an external reference, and cannot calibrate the magnetic measurement system on the projectile separately. Some scholars proposed some methods with EKF [18][19][20][21] or UKF [22][23][24] to estimate measured value, however, the nonlinear model has errors in linearization, and magnetic measurement systems are affected by the hard magnetism and soft magnetism; they are not applicable for magnetic measurement systems, and the errors need to be eliminated effectively in advance. Yanke Wang proposed a calibration method of magnetometer based on a BP neural network [25]; to improve the training speed and convergence of neural networks, the Levenberg Marquardt backpropagation training algorithm is designed to calibrate the magnetometer. This method can effectively reduce the error caused by the change in magnetometer parameters; the measurement error of a magnetometer can be less than 10 nT, however, the influence of misalignment errors is not considered in their model. Further, the prediction ability of the neural network is related to the selected training model; when predicting different measurement systems, the surrounding environment is different, and the required training models are also different. In addition, the BP neural network method has an over-fitting phenomenon, and the limit position of over-fitting is also difficult to determine to ensure the optimal prediction. Therefore, aiming at the various errors of measurement systems and inconvenient conditions, this paper proposes a complete calibration scheme without additional reference and uses ellipsoid fitting based on the Least Square method and TPC method to obtain equivalent parameters and misalignment angles to finish fast calibration for measurement systems on site.

Model of Measurement System
Although micro-mini magnetic sensors have different sensitivity mechanisms to magnetic fields and different sensing characteristics, manufacturing errors are related to various factors. For the convenience of analysis, considering the relationship of the input-output of the magnetic sensor, the errors on the measurement system mainly include two types: single-axis sensing error and three-axis coupling error [26].

Error of Single-Axis Sensor
The mathematical model of the single-axis magnetic sensor can be expressed as: where, k n,i (i = x, y, z) is standard sensitivity, δk i (i = x, y, z) is deviation of sensitivity, k i (i = x, y, z) is real sensitivity and B s 0,i (i = x, y, z) is zero offset.

Error of Triaxial Sensor
The errors of the triaxial magnetic sensor are mainly the mismatch error of zero-offset among the three axes, the static sensitivity mismatches error and the non-orthogonal error among the three axes [27]. The mathematical model of zero-offset and sensitivity can be expressed as: where k M,n is the theoretical sensitivity of the sensor and k M,x , k M,y and k M,z are real sensitivities of three axes respectively. Ideally, k M,n = 1.
The non-orthogonal error among the axes is the measurement error caused by the fabrication that the three sensitive axes cannot be orthogonal completely to each other during manufacturing [28]. For the convenience of analysis, determine the ideal coordinate system as 's ' system, determine the real coordinate system as 's' system; Ox s' coincides with Ox s , illustrated in Figure 1. Considering that the non-orthogonal angles δ M,1 , δ M,2 and δ M,3 are all small angles, the non-orthogonal angle matrix C M,nonorth can be simplified as: The mechanical misalignment angles of the magnetic sensor mainly consist of welding alignment errors and assembling alignment errors.

Misalignment Angles
After orthogonal calibration, the measurement coordinate system and the projectile coordinate system are both the orthogonal coordinate system. The conversion relationship between the two can be regarded as the rotation transformation between the two orthogonal coordinate systems, described by Euler angles. The rotation sequence of the two coordinate systems is defined as that the projectile coordinate system first rotates α y around the Y-axis, then rotates α z around the Z-axis and finally rotates α x around the X-axis, then the projectile coordinate system coincides with the orthogonal magnetic sensitive coordinate system, as shown in Figure 2. The conversion matrix is defined as C b s' .
cos α y cos α z sin α z − sin α y cos α z sin α x sin α y − cos α x cos α y sin α z cos α x cos α z sin α x cos α y + cos α x sin α y sin α z cos α x sin α y + sin α x cos α y sin α z − sin α x cos α z cos α x cos α y − sin α x sin α y sin α z   (4)

Hard Magnetism
The magnetic field is generated by the hard magnetic materials in the projectile [29]. Hard magnetic materials have high coercivity and remanence value. Once magnetized by an external magnetic field, the remanence can be retained for a long time and will not disappear easily. The projectile body is mainly made of some ferromagnetic materials such as steel alloy; the shells of much of the electromechanical equipment on the projectile are also made of ferromagnetic materials. During the manufacture or storage of the projectile, it is parked at a fixed position for a long time, and it is continuously magnetized by the geomagnetic field in that direction at that location, so that the hard magnetic material in the projectile is magnetized and has strong magnetism. When the sensor on the projectile, the error of hard magnetic materials can be equivalent to a bias, which can be shown as:

Soft Magnetism
Soft magnetism is produced by the soft magnetic materials in the equipment. For the convenience of analysis, the soft magnetic material in the projectile is decomposed into countless, infinitely thin soft magnetic wires parallel to the axial, vertical and transverse directions of the projectile coordinate system equivalently. The axial, vertical and transverse soft magnetic wires are affected by the geomagnetic field, respectively. The axial soft magnetic filament clusters are magnetized by the O x b axial component B b G,x of the geomagnetic field in the equipment coordinate system, and the inductive magnetic field is proportional The inductive magnetic field vector of the axial soft magnetic filament cluster is shown in Figure 3, which is decomposed by projection on the coordinate axis of the projectile coordinate system In the same way, the vertical soft magnetic filament clusters and transverse soft magnetic filament clusters of the projectile are magnetized by the geomagnetic field components B b G,y and B b G,z , respectively, and all the vertical soft magnetic filament clusters and transverse soft magnetic filament clusters on the projectile are correspondingly magnetized. The magnitudes of the induced magnetic field before and after magnetization of the ge- in the projectile coordinate system can be shown as:

Eddy Current Magnetic Field
The eddy current magnetic field is excited by the electric eddy current to generate a corresponding inductive magnetic field around it. The eddy current magnetic field is related to factors such as the rate of change of the environmental magnetic field, the motion parameters of the conductor, the geometry of the conductor, the permeability and the conductivity of the conductor. During the external ballistic flight of the projectile, the eddy current magnetic field of the projectile is mainly caused by a series of closed loops in the metal of the projectile shell when the projectile performs high-speed axial translation, rotation around the axis or maneuvering in the outer ballistic flight. The magnetic flux passing through the enclosed area of the closed loop changes with the changes of the threeaxis components of the geomagnetic field in the projectile coordinate system, and then generates eddy currents in the loop. This electric eddy current generates an eddy current magnetic field to prevent this change [30,31].
According to the theory about electromagnetic field, the influence of the eddy current magnetic field on the projectile magnetic sensor can be described as: The changes in the projectile are mainly caused by the steering gear, the rudder blades and the attitude adjustment engine, etc. [32]. For the convenience of analysis, the eddy current magnetic field of the projectile can be described as a three-dimensional random vector, that is: It meets the following conditions: Therefore, the actual measured magnetic field value of the projectile can be expressed as: Substituting it into the magnetic sensor measurement model, the mathematical model of the measurement information of the projectile can be obtained as: Then, the relationship of input and output of the magnetic sensor can be simplified as: where K M is an equivalent sensitivity matrix, K eddy is an equivalent eddy current sensitivity coefficient matrix, B s offset is an equivalent zero-offset and B s noise is an equivalent noise in output whose mean value is zero.
If the influence of the eddy current magnetic field of the projectile is ignored, it can be approximated as a linear model: where: The parameters (equivalent sensitivity coefficient matrix K M , equivalent zero-offset B s offset ) in the output information of the magnetic sensor are obtained by the first stage calibration method, and the geomagnetic information can be estimated by the following formula:

Calibration Scheme
According to the analysis above, the calibration process can be decomposed into two stages: calibration of sensor and calibration of measurement system: (1) The first stage is to finish the unitization and orthogonalization of the coordinate system of the sensor; (2) The second stage is to calibrate misalignment angles between the measurement system and projectile coordinate system by the TPC method based on the first stage.

First-Stage Calibration
In the calibration area where the geomagnetic field is stable, the geomagnetic information is measured by an ideal orthogonal magnetometer that does not contain various measurement error factors. When rotating the sensor in space, the measured value can be described as a sphere [33]. Then, the measured value of the ideal magnetic sensor in different attitudes satisfies the following equation: According to the mathematical model mentioned before in the article, the measured value of the geomagnetic field is: Therefore, there is: After simplification, it satisfies the quadratic standard equation: Then, the measurement data satisfies a quadratic ellipsoid surface equation, and its geometric meaning is: the three-dimensional coordinates determined by the measurement data are on the three-dimensional ellipsoid surface, and the center of the ellipsoid is equivalent zero-offset. The radius of each semi-axis of the ellipsoid is related to the product of the equivalent sensitivity coefficient of the corresponding axis of the magnetic sensor and the modulus of the geomagnetic field. The angle between each semi-axis of the ellipsoid and the sensitive axis is the equivalent non-orthogonal angle between each axis of the magnetic sensor.
Suppose the equation of the ellipsoid as the following equation: where ξ = [a, b, c, d, e, f , p, q, r, g] T are the coefficients that need to be calculated and T are the measured values of the tri-axial measurement vector. Fitting measured data to obtain the optimal ellipsoid fitting coefficients [a, b, c, d, e, f , p, q, r, g] T , in the process of fitting the ellipsoid surface, the minimum square sum of the algebraic distance between the measurement data and the ellipsoid surface is determined as the criterion as: where:  1, 2, . . . . . . N) is the measurement data of the projectile in the i-th rotation attitude. Calculating with vector (X − X 0 ) T A (X − X 0 ) = 1, then matrix. According to the Cholesky decomposition theorem: when P is a real symmetric positive definite matrix of order n, then P has a unique LDU decomposition, shown as: where D is the diagonal matrix D = diag (d1, d2, . . . , dn) and di > 0 (i = 1, 2, . . . n); L is the unit lower triangular matrix; U is the unit upper triangular matrix.
Due to P T = P, then According to the uniqueness of the matrix Cholesky decomposition, U T = L; L T = U.
The matrix C s s in corresponding model is: Then a complete mathematical model can be established: The parameter can be obtained, Therefore, equivalent sensitivity and equivalent non-orthogonal angles are:

Three-Position Method
Although the geomagnetic field is a global long-term changing magnetic field, it changes slowly and is almost constant. Therefore, there is a constraint relationship of the measured values when the measurement system is in different directions. According to this theory, the projectile is placed on a non-magnetic platform, the X-axis, Y-axis and Z-axis of the projectile coordinate system point to the north, up and east, respectively. Record the measured value at the initial position as B 1 s' , rotate 180 • around the other two axes in the projectile coordinate system and record the measured value as B 2 s' and B 3 s' . There are six rotation sequence schemes as shown in Table 1. Considering the ease of operation, without loss of generality, plan 3 is used in the following description, the sequence of rotation is rotating 180 • around the Z-axis after the initial position, then rotating 180 • around X-axis, illustrated by Figure 4.
Determine B s ij = B s i + B s j /2, i and j are the number 1, 2 and 3 of the positions. Using the Newton-Raphson method, suppose the X = [α x , α y , α z , B c x , B c y , B c z ], define the initial value X k , Taylor expands X i (i = 1, 2, 3 . . .) at X k , omit high-order terms with small values. where: J is the Jacobian matrix of the nonlinear equations: − sin α y cos α z 0 0 sin α x cos α y + cos α x sin α y sin α z 0 0 cos α x cos α y − sin α x sin α y sin α z   (33) Using the least-squares method, ∆X= (N T N −1 N T ∆M can be obtained, and there is an iterative calculation X k+1 = X k + ∆X until the required accuracy threshold. When using the Newton-Raphson method, its convergence requires the initial value to be selected appropriately. Define the initial value:

Error Analysis
Suppose the measurement error conforms to the normal distribution, for which the mean value is zero and the three-axis magnetic measurement values are independent of each other, satisfies ∆H i s'~N (0, σ 2 ), as: Then, Then, According to ∆X= (N T N −1 N T ∆M, the statistics of X satisfy the equation When α x , α y , α z are small values, the quadratic term and other higher-order terms in the Taylor expansion can be ignored, and the covariance matrix of ∆X can be simplified as: where Therefore, it can be seen that the calibration accuracy of the misalignment angles of the magnetic field information is related to the measurement accuracy of the three-axis magnetic sensor, the magnetic field vector and the projection component of H c in the calibration coordinate system. In order to estimate α x , α y and α z with equal accuracy, it is necessary to project the magnetic field vector B c to the calibration coordinate system with equal components. Determine the intensity of the geomagnetic field as ||B||; it needs to meet the following equation: The angle between the magnetic field vector of the calibration site and each axis in the calibration coordinate system is ±arc cos 1/ √ 3 = ±54.73 • , as:

Simulation
In order to verify the feasibility of this scheme, suppose the intensity of geomagnetic field is 60,000 nT, add necessary hard magnetism and soft magnetism; the simulation data is shown in Table 2. Table 2. Simulation data.

Hard Magnetism (nT) Soft Magnetism (nT)
[500, 300, −300] [30,100,100] Fitting the data as an ellipsoid surface is shown in Figure 5. According to the analysis above, the parameters are as shown in Table 3. Compensate the measured value and the fitting diagram as shown in Figure 6. It can be seen from Figure 6 that the fitting graph of the measured value after compensation is close to a sphere, which is a proper graph under ideal conditions, which proves that this method can eliminate the influence of the deviation of parameters of sensor and other magnetic fields.
Suppose the misalignment angles are [−1 • , 2 • , 3 • ], the three-axis components in the calibration coordinate system are [35,468,35,468,35,468] nT, as shown in Table 4, to verify the feasibility of the TPC method. To obtain the optimal results, the available methods are the Las Vegas method and the Monte Carlo method. Because of the limited number of points in the actual sampling, it does not meet the applicated conditions of the Las Vegas method, therefore the distribution of misalignment angles is judged by the Monte Carlo method 10,000 times. The results are shown in Table 5 and Figure 7.  It can be seen from Table 5, Figure 7, in the results of the calibration algorithm based on the TPC method, that the errors of the misalignment angle and the magnetic field component all meet the characteristics of normal distribution, their average values are zero and the three misalignment angles are estimated with equal accuracy when the components are equal. The standard deviation of components is close to the theoretical value 70.7107 nT, it is consistent with the analysis in Section 3 of this paper.

Experiment and Discussion
In order to verify the scheme described in the article, HMC1053 is the magnetometer in this experiment. The magnetoresistive sensor (AMR) of Honeywell of the United States represents the highest level in the industry; the miniature three-axis magnetoresistive sensor HMC1053 has the characteristics of small size, low power consumption and fast dynamic response. It is widely used in the navigation test of unmanned aerial vehicles, unmanned vehicles and various small smart carriers; as well, it has been used in guided munitions of France, and the "CORECT" combined measurement module for rocket ballistic correction jointly developed by Switzerland and Germany. In this experiment, the sensor is installed on the board by welding; it is important to install the measurement system on a piece of equipment made of metal materials as shown in Figure 8, FVM400 fluxgate is used to measure an appropriate area where the geomagnetic field is stable for calibration. The accuracy of the device is 1 nT. Place the system at a certain position and rotate it around the X-axis, Y-axis and Z-axis with multiple attitudes. Without loss of generality, the selected method is to rotate 360 • on the vertical plane when the system is moving 45 • on the horizontal plane. Obtain measured data of multiple postures, obtain the modulus as shown in Figure 9.  Figure 9 shows that, before calibration, the peak-to-peak variation range of the magnetic field modulus is about 44,090 nT. The ellipsoid fitting graph is as shown in Figure 10. Equivalent zero offset, equivalent sensitivity and equivalent non-orthogonal angles are shown in Table 6. The modulus of measured value after compensation in Figure 11.  Figure 11 shows that after the first stage calibration, the peak-to-peak variation range of the magnetic field modulus is about 2200 nT. Compared with the modulus before the first stage, the error reduces by about 20 times.
Then, to calibrate misalignment angles with TPC method, a straight line in the plane of a plastic plank is selected as the X-axis of the calibration coordinate system, the straight line perpendicular to the X-axis in this plane is the Z-axis and the calibration coordinate system is determined according to the right-hand rule. Using FVM400 to measure the magnetic field components and place the plate in the proper position so that the components in the calibration coordinate system are approximately equal, the components are 35,977 nT, 35,871 nT and 36,100 nT.
First, place the projectile on the calibration plane so that the projectile coordinate system is consistent with the calibrated coordinate system, marked as position 1, as shown in Figure 12. Then, rotate 180 • around the Z-axis of the projectile coordinate system, marked as position 2 as shown in Figure 13. Then, rotate 180 • around the X-axis of projectile coordinate system, marked as position 3 as shown in Figure 14. Record the measured value of the three positions in Table 6. It can be seen from Table 7 before compensating for misalignment errors, the error of the single axis is thousands of nT and, due to the existence of the misalignment angles, the measured values of the X-axis and the Z-axis have also changed when rotating around the X-axis and Z-axis; in theory, they should not change. According to the calibration algo-rithm elaborated in Section 3, the misalignment angles are [−2.7903 • , −3.2721 • , 5.0245 • ] and the compensated values are obtained as shown in Table 8. The results show that the TPC method can calibrate misalignment angles, compensate measured value according to analysis in Section 2, the error after calibration is within 150 nT and the maximal error is reduced from 6028 nT to 144 nT. The calibration improves the accuracy of the measured value significantly, which indicates that the misalignment error is eliminated effectively. The results of the experiment show that the proposed scheme has high precision, and the complete process can be finished in a few minutes, which has high operability and applicability, making it suitable for calibrating magnetic measurement systems on guided munition on site.

Conclusions
In this paper, a TSC scheme without additional reference to calibrate tri-axial magnetic measurement system on guide munition is presented. This scheme solves the problem that the projectile coordinate system is not parallel to the coordinate system of the magnetic measurement system, it is suitable for the calibration of the measurement system of various types of projectiles on site without other devices or references.
The simulation and semi-physical experiment verified that the error of the sensor can be eliminated effectively by the first stage and misalignment angles can be calibrated by the second stage. After the complete scheme is applied, the deviation of the measured value has reduced nearly two orders, accuracy of the magnetic measurement system meets the requirements for estimation of the attitude of guided munitions.
Author Contributions: Conceptualization, Y.X. and X.Z.; literature search, Y.X. and X.Z.; methodology, Y.X. and X.Z.; validation, Y.X. and X.Z.; writing-original draft preparation Y.X.; writing-review and editing X.Z.; project administration X.Z. Both authors have read and agreed to the published version of the manuscript. Data Availability Statement: All data, and related codes of the calibration scheme will be sent to the e-mail of the corresponding author upon request, and appropriate reasons will be provided.