Compensation of Optical Pump Magnetometer Using the Improved Mayﬂy Optimization Algorithm

: In order to solve the problem that the cesium optical pump magnetometer is disturbed by the carrier’s interference magnetic ﬁeld during magnetic ﬁeld anomaly detection, an interference magnetic ﬁeld compensation method based on an improved mayﬂy optimization algorithm (IMOA) was proposed in this paper. First, by combining the measurement results of the attitude sensor with the geomagnetic inclination and magnetic declination in the locality, the measurement results of the optical pump magnetometer can be decomposed into the component values under the three axes of the carrier coordinate system. A compensation model including the carrier interference magnetic ﬁeld was established. Then, considering the poor global search performance that existed in the mayﬂy optimization algorithm (MOA), an elite chaotic reverse learning strategy and Levy mutation strategy were introduced to improve the MOA. The compensation performance of the IMOA was estimated with a series of ﬁeld experiments and compared with the stretching particle swarm optimization algorithm. The experiment results indicated that these two methods can effectively compensate the magnetometer’s measurement values, and that the IMOA method more easily jumps out of the local optimum, and has higher compensation accuracy.


Introduction
The geomagnetic field is a passive signal and has the characteristics of one-to-one correspondence with the spatial position.However, the ferromagnetic material on the earth's surface will generate an induced magnetic field due to the action of the geomagnetic field, which will cause the local geomagnetic field to be distorted.The geomagnetic field distortion also provides a theoretical basis for the detection of underwater magnetic targets, that is, the detection of magnetic anomalies [1,2].The magnetic anomaly detection system based on the underwater unmanned vehicle (UUV) in the deep-sea environment can be close to the magnetic target and obtain anomalous magnetic signals with high signal-tonoise ratio.Therefore, as a highly maneuverable and flexible measurement method, UUV equipped with magnetometers began to be used for a series of tasks such as measuring the magnetic characteristics of magnetic targets on the water surface [3,4], locating and tracking underwater magnetic substances [5], and detecting the magnetic field of the seabed [6].
However, the actual measurement result of the magnetometer is superimposed by the earth's magnetic field, the magnetic field of the magnetic target, and the interfering magnetic field of the UUV.Therefore, when the local standard geomagnetic field is known, the target magnetic field may also be submerged in the interference magnetic field of the UUV, resulting in the failure of the magnetic anomaly detection mission.Therefore, the UUV's interference magnetic field compensation method has become one of the key factors restricting the development of underwater magnetic anomaly detection technology.In order to improve the measurement accuracy of the abnormal magnetic fields, it is necessary to conduct in-depth research on the removal of the UUV's interference magnetic fields.
For the compensation research of the magnetometer, Tolles [7] and Leliak [8] decomposed the aircraft's disturbance magnetic field into a constant magnetic field related to hard magnetic materials, an induced magnetic field caused by the change of carrier attitude, and an eddy current magnetic field.The corresponding mathematical model was established with the Tolles-Lawson equation (T-L equation), so that the compensation problem of the carrier interference magnetic field was transformed into the linear regression problem of the compensation parameters.Gebre-Egziabher et al. [9] proposed a two-step calibration method based on the principle that the theoretical magnetic field strength after calibration is equal to the reference magnetic field strength.This method transformed the calibration compensation problem into a linear regression problem by introducing intermediate parameters.The linear least squares method (LSM) was first used to solve the intermediate parameters, and then the calibration compensation parameters were solved by the algebraic solution method.
Aiming at the problem of errors in the measurement data in the two-step method, Wu et al. [10] and Zhang et al. [11] proposed to use the overall least squares method and the truncated overall least squares method to solve the intermediate parameters, respectively.Based on the principle that the component information of the magnetic field after compensation is equal to the reference magnetic field, Zhang et al. [12] introduced intermediate parameters to linearize the compensation parameters and obtained the compensation parameters by the least squares method.Pang et al. [13] proposed a calibration compensation method based on a two-step LSM.After obtaining the intermediate parameters, the LSM was also used to solve the compensation parameters, which further increased the solution accuracy of the compensation parameters.
Yu et al. [14] and Vasconcelos et al. [15] fit the ellipsoid equation based on the magnetometer measurement value, and then used the LSM and the maximum likelihood estimation method to iteratively solve the compensation parameters.In order to fit the ellipsoid, this method requires the magnetometer to traverse various spatial attitude angles as much as possible, and the variation of the roll angle and pitch angle during the actual UUV navigation process is very limited, so this compensation method is not suitable for the compensation of UUV's interference magnetic field.
There are also many scholars who have transformed the compensation problem of the carrier interference magnetic field into the nonlinear parameter optimization problem, and research has been made on the optimal estimation method of compensation parameters.Pang et al. [16] established the calibration model of the fluxgate magnetometer, and the calibration parameters were iteratively solved using the "Fsolve" function in MATLAB.Li et al. [17] and Li et al. [18] used the UKF method and the trust region method to estimate the compensation parameters, respectively.These two methods can obtain better estimation results when the selection of the iteration initial value is appropriate.The above compensation parameter estimation results are greatly affected by the initial values of the parameters, and inappropriate initial values are also likely to cause a divergence of results.
Research using meta-heuristic algorithms has determined that this kind of algorithm has strong applicability to high-dimensional nonlinear optimization search problems, and has the advantages of insensitivity to the initial value, a simple structure and no gradient mechanism [19].Some scholars have begun to introduce heuristic algorithms into the estimation of compensation parameters, and have proved the superiority of such algorithms.The particle swarm optimization (PSO) algorithm has the advantages of insensitivity to the initial value and a fast convergence speed.Accordingly, many scholars have improved this method and applied it to the compensation of the carrier interference magnetic field.For example, Wu et al. [20] and Wu et al. [21] calibrated and compensated the measurement results of a fluxgate magnetometer based on a PSO and a stretching PSO algorithm.Li et al. [22] improved the inertia weight in the PSO algorithm and proposed a damping PSO algorithm, and the results showed that the effect is better than the two-step method.Huang et al. [23] proposed a magnetometer measurement error compensation method based on an immune adaptive particle swarm optimization algorithm.Zhang et al. [24] used the differential evolution (DE) algorithm to estimate the compensation parameters and the local magnetic field.Gao et al. [25] introduced the cuckoo optimization search (CS) algorithm into the compensation of the magnetometer.The experimental results showed that the compensation accuracy of this method is better than that of the UKF method.
The above compensation method is mainly based on the measurement results of the fluxgate magnetometer.Due to the limitations of production, processing and assembly, the fluxgate magnetometer has non-orthogonal errors and zero-scale factor errors [26].It needs to be specially calibrated, and the accumulation of errors will inevitably reduce its compensation accuracy.As a scalar magnetometer, the cesium optical pump probe has high geomagnetic measurement accuracy, no zero drift and is not affected by temperature.It does not need accurate orientation during operation, and is suitable for high-precision, rapid, and continuous measurement under moving conditions, such as airborne magnetic surveys and marine magnetic surveys.
At present, the compensation research for the cesium optical pump magnetometer's measurement error mainly focuses on the aeromagnetic compensation method based on the T-L model.This compensation method assumes that the absolute value of the disturbance magnetic field at the scalar magnetometer is the projection of the vector disturbance magnetic field in the direction of the geomagnetic field [27].This assumption ignores the angle between the direction of the geomagnetic field and the direction of the interference magnetic field, so it is suitable for the case where the interference magnetic field is extremely small.In addition, the direction cosine is mostly provided by the fluxgate magnetometer in the aeromagnetic compensation, but the data collected by the fluxgate are not accurate, and the loading position of the fluxgate cannot be exactly the same as that of the optical pump magnetometer [28].Therefore, there is a certain deviation in the direction cosine value obtained by the fluxgate measurement result, which may lead to poor estimation accuracy of the compensation parameters.Wu et al. [29] decomposed the total field data of the optical pump in the three-axis direction and used the genetic algorithm (GA) to estimate the compensation parameters, but this method used the three-axis magnetometer to calculate the three direction cosines containing the interference external magnetic field, which ignores the error of the fluxgate magnetometer itself, and thus may reduce the compensation accuracy.
In order to solve the problem of interference magnetic field when the remotely operated vehicle (ROV) is equipped with a cesium optical pump magnetometer for magnetic anomaly detection, this paper proposed an improved mayfly optimization algorithm to estimate the compensation parameters.The inertial navigation sensor was introduced to help decompose the total magnetic field measurement result of the optical pump magnetometer.Then, a compensation model based on the component value of the measurement results was constructed.The mayfly optimization algorithm (MOA) is a swarm intelligent optimization algorithm proposed by Zervoudakis and Tsafarakis [30] in 2020.Compared with other intelligent optimization algorithms, this algorithm has certain advantages in local search performance and population diversity [31] and is suitable for high-dimensional nonlinear parameter estimation.However, the algorithm also easily falls into local optimum, so this paper introduced the elite chaotic reverse learning strategy and Levy mutation strategy to improve the mayfly algorithm, and the improved mayfly optimization algorithm (IMOA) was used to estimate the compensation parameters.Finally, the practicability of the proposed compensation method was verified with a series of field experiments.
To summarize, the main contributions of this paper are as follows: • An ROV's interference magnetic field removal method based on the measurement results of a cesium optical pump magnetometer was proposed in this paper, which used the inertial navigation sensor to help decompose the magnetic field measurement results; • This paper proposed an improved mayfly optimization algorithm, which improved the application ability of this algorithm in the field of estimation of the compensation parameters;

•
The proposed compensation method can achieve high compensation accuracy on actual magnetic field data and has feasibility and validity for different experiment data.

The Coordinate System
When the ROV is equipped with a magnetometer to measure the magnetic field information in the water, the posture of the ROV will inevitably change.In order to illustrate the measurement results of the magnetometer more conveniently, it is necessary to introduce the following coordinate systems: the first is the geographic coordinate system (also known as the north-east-down coordinate system).The direction of its x g axis points to the north, the direction of the y g axis points to the east, and the direction of the z g axis is perpendicular to the local horizontal plane downward.Then, a carrier coordinate system is established, in which the x c -axis points forward along the ROV's sailing direction, the y c -axis points to the right-hand direction perpendicular to the x c -axis, and the z c -axis points downwards perpendicular to the x c y c plane.
According to the three-dimensional space coordinate system conversion principle, the coordinate conversion relationship between the geographic coordinate system and the carrier coordinate system is shown in Equation (1) [25].Among them, H c e and H g e represents the projection of the geomagnetic field in the carrier coordinate system and the geographic coordinate system, respectively.K c g represents the coordinate transformation matrix from the geographic coordinate system to the carrier coordinate system.α, β, γ are the heading angle, pitch angle, and roll angle measured by the inertial navigation sensor, respectively. (1)

Interference Magnetic Field Mathematical Model of the ROV
Due to the slow sailing speed of ROV, the ROV's interference magnetic field mainly includes the permanent interference magnetic field generated by the hard magnetic material and the induced magnetic field generated by the soft magnetic material.
The permanent interference magnetic field can be regarded as the superposition of the magnetic fields generated by all the magnetic domains in the fixed hard magnetic material at a certain point in space.The magnetic moment of these magnetic domains will not change for a long period of time.Therefore, the magnitude of the permanent interfering magnetic field at a point near the ROV depends on the distance between the point and the ROV, regardless of the external environment.The permanent interference magnetic field is represented by H p , and its expression in the carrier coordinate system is shown in Equation (3).
The induced interference magnetic field is mainly caused by the magnetization effect of the external magnetic field on the soft magnetic materials constituting the ROV.The magnitude and direction of the induced interference magnetic field are mainly determined by the external magnetic field.Therefore, the magnitude and direction of the induced disturbance magnetic field at a point near the ROV changes with the variation of the ROV's attitude.The induced interference magnetic interference is represented by H i , and its expression in the carrier coordinate system is shown in Equation ( 4) [18].The K i means the coefficient matrix of the induced interference magnetic field, and the value of each element of the coefficient matrix depends on the relative permeability of the soft magnetic material and the distance between the measurement point with the ROV.
Based on the above analysis, the ROV's interference magnetic field at point P in space can be expressed by Equation (5), and the schematic diagram is shown in Figure 1.
The induced interference magnetic field is mainly caused by the magnetization effect of the external magnetic field on the soft magnetic materials constituting the ROV.The magnitude and direction of the induced interference magnetic field are mainly determined by the external magnetic field.Therefore, the magnitude and direction of the induced disturbance magnetic field at a point near the ROV changes with the variation of the ROV's attitude.The induced interference magnetic interference is represented by Hi, and its expression in the carrier coordinate system is shown in Equation ( 4) [18].The Ki means the coefficient matrix of the induced interference magnetic field, and the value of each element of the coefficient matrix depends on the relative permeability of the soft magnetic material and the distance between the measurement point with the ROV.
Based on the above analysis, the ROV's interference magnetic field at point P in space can be expressed by Equation ( 5), and the schematic diagram is shown in Figure 1.According to the literature [32], when the ROV's thruster rotates, the thruster will generate a leakage magnetic field with a sinusoidal periodic distribution.The static interference magnetic field of the ROV described above and the geomagnetic field change very slowly compared with the leakage magnetic field.Therefore, the data collector uses low-pass filtering to remove the high-frequency leakage magnetic field in the measured magnetic field signal in this experiment.Then, the filtered magnetic field measurement results are saved and used in the subsequent compensation parameter evaluation process.

The Compensation Model
At the optical pump magnetometer, there is a vector relationship between the ideal geomagnetic field, the ROV's interference magnetic field, and the total magnetic field measurement results, as shown in Figure 2.Among them, the Hm is the vector expression of the optical pump magnetometer's measurement result, He represents the ideal geomagnetic field.The scalar measurement result of the optical pump magnetometer can be regarded as the modulus value of the vector sum of the earth's magnetic field and the ROV's interference magnetic field at the magnetometer, from which Equation ( 6) can be obtained.According to the literature [32], when the ROV's thruster rotates, the thruster will generate a leakage magnetic field with a sinusoidal periodic distribution.The static interference magnetic field of the ROV described above and the geomagnetic field change very slowly compared with the leakage magnetic field.Therefore, the data collector uses low-pass filtering to remove the high-frequency leakage magnetic field in the measured magnetic field signal in this experiment.Then, the filtered magnetic field measurement results are saved and used in the subsequent compensation parameter evaluation process.

The Compensation Model
At the optical pump magnetometer, there is a vector relationship between the ideal geomagnetic field, the ROV's interference magnetic field, and the total magnetic field measurement results, as shown in Figure 2.Among them, the H m is the vector expression of the optical pump magnetometer's measurement result, H e represents the ideal geomagnetic field.The scalar measurement result of the optical pump magnetometer can be regarded as the modulus value of the vector sum of the earth's magnetic field and the ROV's interference magnetic field at the magnetometer, from which Equation ( 6) can be obtained.In actual measurement, the local magnetic inclination angle I and magnetic declination angle D can be obtained according to the local longitude and latitude information.Thereby, the scalar measurement value of the optical pump magnetometer can be decom- In actual measurement, the local magnetic inclination angle I and magnetic declination angle D can be obtained according to the local longitude and latitude information.Thereby, the scalar measurement value of the optical pump magnetometer can be decomposed into three components in the geomagnetic coordinate system using the Equation (7), namely H g mx , H g my and H g mz .Combined with the attitude information measured by the inertial navigation sensor, the measurement result of the magnetic field component containing the interfering magnetic field in the carrier coordinate system can be obtained.The measurement results in the carrier coordinate system include the earth's magnetic field and the ROV's interference magnetic field, as shown in Equation (8).
Combining the Equations ( 1), ( 5) and ( 8) and performing mathematical transformation, Equation ( 9) can be obtained, where H c e is the geomagnetic field value in the carrier coordinate system after compensation.It can be seen from ( 9) that the key to compensate ROV's interference magnetic field is to solve two sets of compensation parameters, that is, the coefficient matrix K i and the vector H p .Therefore, it is necessary to construct a multi-objective optimization function as shown in Equation (10).According to Equation (10), at least three sets of equations are required to obtain the compensation parameters of the magnetometer.Therefore, in order to solve this optimization problem, it is necessary to select an area where the magnetic field gradient changes very little.The geomagnetic field intensity in this area can be regarded as a fixed value in a short time.The carrier is controlled to conduct several attitude transformations in this area to obtain multiple sets of magnetic field measurement values.When the reference value of the background geomagnetic field, magnetic inclination, and magnetic declination are known, the optimal estimation result of the compensation parameters can be obtained with a certain optimization algorithm.Then, the ROV's interference magnetic field can be removed from the measurement results using the Equation (9).

The Optimal Estimation Algorithm
In this section, we propose a compensation parameters estimation method based on the improved mayfly optimization algorithm (IMOA).The original mayfly optimization algorithm (MOA) is first described.Then, an IMOA based on an elite chaotic reverse learning strategy and Levy flight strategy is proposed.Finally, the compensation parameters estimation steps based on IMOA are given.

The Primitive Mayfly Optimization Algorithm
Like the classic PSO algorithm and GA, the MOA is also a heuristic search algorithm, which was inspired by the social behavior of mayfly groups, including behaviors such as aggregation and mating.The mayfly algorithm combines the advantages of PSO, GA and other algorithms in its structure.It not only has the ability of random search, but also enables the population to evolve according to the law of survival of the fittest.Compared with other algorithms, this algorithm possesses the advantages of strong local search performance, high population diversity, and fast convergence speed.Therefore, it is suitable for solving high-latitude nonlinear optimization problems.
Assuming that in the d-dimensional search space, the number of both the male and female mayflies in the population is m, the position of the male mayfly and female mayfly is expressed as x = (x 1 , x 2 , . . ., x d ) and y = (y 1 , y 2 , . . ., y d ), respectively.The speed of each mayfly is represented by v = (v 1 , v 2 , . . ., v d ).The quality of each mayfly's search result is evaluated by the objective function f (x) defined in the previous section.
After the t-th iteration, the position, velocity, and historical optimal position of the i-th male mayfly in the j-th dimension search space are denoted as x t ij , v t ij , and p ij , respectively, and the population optimal position in the j-th dimension search space is denoted as g j .Since male mayflies generally gather and dance on the water's surface to attract female mayflies, each male mayfly adjusts its flight speed according to its own flying experience and that of the group.Then each male mayfly updates its current speed and position according to Equations ( 11) and ( 12): Among them, a 1 , a 2 are the attraction coefficients of male mayflies, and β is the visibility coefficient.The Euclidean distance between the x i with the p ij is represented as r p .The r g represents the Euclidean distance between the x i and the g j .d represents the dance coefficient, which is used to attract female mayflies, r is a random coefficient, and r ∈ [−1 , 1].The mayfly at the optimal position continuously updates its speed and position by introducing random element d to lead the mayfly group to fly to a better position.
Different from male mayflies, female mayflies will decide whether to fly close to the male mayfly according to the quality of its position.After the t-th iteration, assuming that the position and velocity of the i-th female mayfly in the j-th dimensional search space are y t ij and v t ij , the female mayfly updates its current velocity and position according to Equations ( 13) and ( 14): where a 3 is the attraction coefficient of the female mayfly, r m represents the Euclidean distance between male and female mayflies.The fl is the random walk coefficient, which means that the female mayfly will fly randomly to find a better male mayfly when not attracted.The ultimate goal of the mayfly population is to mate to produce better offspring.In order to avoid falling into the local optimum, the mating process in the algorithm is mainly to randomly select a part of the samples in males and females, respectively.According to the mechanism of survival of the fittest, the male optimal individual and the female optimal individual mate, and the second optimal male individual mates with the female individual.The offspring produced after mating are generated by Equations ( 15) and (16).
Among them, off 1 and off 2 represent the two generated offspring, L is a random number and L ∈ [−1, 1], m represents a male mayfly, and f represents a female mayfly.

The Improved Mayfly Optimization Algorithm
Like most heuristic algorithms, the MOA also has the problem of easily falling into premature convergence in high-dimensional complex nonlinear problems.On the one hand, the MOA uses random initialization to generate the initial population which easily leads to a certain blindness in the iterative search process and may result in long search time, slow convergence speed, and the formation of local optimum.On the other hand, in the later stage of iteration, the inertia weight factor, dance coefficient and random walk coefficient will decrease gradually in order to meet the local search performance.If the algorithm falls into a local optimum in the early stage, it is difficult to jump out of the local optimum in the later stage.
In view of the above defects, we adopted the following measures to improve MOA: (1) we introduced the elite chaotic reverse learning strategy into the population initialization of the MOA.By combining the chaotic solution and the reverse solution, this strategy selects the mayfly population on a global scale to improve the initial population's diversity and quality, avoiding the uncertainty caused by the random initial population, and improving the convergence efficiency; (2) the Levy mutation strategy was used to perturb some mayfly individuals that may fall into the local optimum, so as to improve the population's diversity and anti-stagnation ability, and help the algorithm to jump out of the local optimum.
The elite chaotic reverse learning strategy mainly adopts two stages to initialize the mayfly population.Among the chaotic operators, the cubic chaotic mapping has better uniformity, so it can be considered to use cubic chaotic mapping to initialize the mayfly population.The initialization process is shown in Equations ( 17) and (18).First, the cubic mapping variable is generated by Equation (17), where y i ∈ [−1, 1], and the initial value y 0 is randomly generated.Then, the cubic mapping variables are applied to the mayfly population using Equation (18), where x i is the i-th mayfly after the cubic chaotic mapping, lb and ub are the lower and upper bounds of the mayfly population position, respectively.
On the basis of the cubic chaotic mapping, the lens imaging reverse learning strategy [33] is introduced to reversely solve the above mayfly population.The specific mathematical model is shown in Equation ( 19).Among them, z i is the i-th mayfly after the reverse learning, and η is the lens's zoom factor.
Two mayfly populations are obtained after the cubic chaotic mapping and the lens reverse learning strategy.In order to select the mayfly with better quality, the mayfly population is selected according to the fitness order.The mathematical model of screening is shown in Equation (20).The i-th mayfly individual selected is represented by X i : In the later stage of MOA iteration, mayfly individuals will rapidly assimilate, which may easily lead to the stagnation of the mayfly population at the local optimal position.In order to solve this problem, according to the roulette probability, some mayfly individuals are randomly selected for Levy mutation.Then the positions of the mutated mayfly before and after the mutation are compared according to the objective function.If the mutated position is better, the mutated mayfly individual is chosen to enter the next iteration, otherwise the mutation is invalid.
The Levy mutation strategy comes from Levy flight, which is a non-Gaussian random process proposed by French mathematician Levy.Levy flight is also an ubiquitous phenomenon in nature.Due to its large and small flight steps, both global optimization and local optimization can be taken into account in the optimization problem, which helps the population fly out of the local optimum.The Levy variation is shown below: Among them, β ∈ [0, 2], µ obeys N (0, σ 2 ) distribution, and ν obeys N (0, 1) distribution.

The Compensation Parameters Estimation Process Based on the IMOA
The main procedures for estimating compensation parameters using the IMOA are as follows: Step 1: Setting the parameters in the algorithm, mainly including the population size, the maximum number of iterations, the upper and lower bounds of the compensation parameters, the attraction coefficients, the dance coefficient, and the random walk coefficient; Step 2: Initializing the positions of male and female mayflies, respectively, using the elite chaotic reverse learning strategy, setting the initial speed value and calculating the initial fitness function value; Step 3: Performing cyclic calculation on the male mayfly.Updating the speed, position and fitness function value of the male mayfly using the Equations ( 10)-( 12) and then updating the global optimal fitness value; Step 4: Performing cyclic calculation on the female mayfly.Updating the speed, position and fitness function value of the female mayfly using the Equations ( 13), ( 14) and (10) and then updating the global optimal fitness value; Step 5: Selecting partial samples from the updated mayfly populations randomly.Obtaining the progeny mayfly according to the Equations ( 15) and ( 16) and then updating the global optimal fitness value; Step 6: Deciding whether to perform Levy mutation according to the roulette probability and then updating the global optimal fitness value; Step 7: Ending the iteration if the maximum number of iterations is reached, otherwise go to the Step 2 to reiterate.

Field Experimental Verification and Analysis
The field experiment was carried out and a series of magnetic field data were collected in order to verify the performance of the compensation method proposed in this paper.The experiment data collection was divided into two groups, the first group of experiment data was used to estimate the compensation parameters, and the second group of data was used to verify the application ability of the estimated compensation parameters.

Field Experimental System and the Measurement Results
The experimental apparatus mounded on the ROV mainly consisted of an inertial navigation sensor (to collect the ROV's attitude information), a data collector (to collect measurement results and perform filtering), and a cesium optical pump magnetometer.Since the ROV is instrumented with a lot of electronic devices and ferromagnetic substances that can generate magnetic fields, it is a challenging task to compensate for ROV's interference magnetic field during sailing.Therefore, a certain installation distance between the magnetometer and the ROV should be taken into account to isolate some stray magnetic fields.As shown in Figure 3, the cesium optical pump magnetometer is fixed inside a sealed cabin and is installed two meters in front of the ROV through an aluminum alloy measurement rod.The main performance specifications of the magnetometer are listed in Table 1.Except that the propeller can rotate around its shaft, other parts of the ROV are fixed by aluminum alloy screws.Therefore, in addition to the leakage magnetic field generated when the thruster works, the interference magnetic field generated by ROV is mainly a static interference magnetic field.
stances that can generate magnetic fields, it is a challenging task to compensate for ROV's interference magnetic field during sailing.Therefore, a certain installation distance between the magnetometer and the ROV should be taken into account to isolate some stray magnetic fields.As shown in Figure 3, the cesium optical pump magnetometer is fixed inside a sealed cabin and is installed two meters in front of the ROV through an aluminum alloy measurement rod.The main performance specifications of the magnetometer are listed in Table 1.Except that the propeller can rotate around its shaft, other parts of the ROV are fixed by aluminum alloy screws.Therefore, in addition to the leakage magnetic field generated when the thruster works, the interference magnetic field generated by ROV is mainly a static interference magnetic field.

Performance
Parameter Value measurement range 10,000~1,000,000 nT sampling frequency 3~100 Hz sensitivity 0.002 nT dynamic noise 0.1 nT The experimental steps are as follows.First, an experiment site with clean magnetic field was selected and the magnetic field information of this site was obtained.The ROV was controlled to change its pitch angle, roll angle and heading angle many times.At this time, the measurement results of the magnetometer mainly included the geomagnetic field and the interference magnetic field of ROV.Then, the compensation parameters were evaluated using the method in Section 3. Finally, according to the compensation parameters, the interference magnetic field of ROV was calculated and removed from the measurement results.At this time, only the geomagnetic field remained in the measurement results.Of course, if this technology is applied to the actual abnormal magnetic field detection, after removing the interference magnetic field of ROV, the remaining magnetic field is the superposition of the abnormal magnetic field of the magnetic target and the geomagnetic field.
The field experiment site was selected in a water area with a uniform magnetic field far away from artificial constructions.In order to improve the accuracy of estimating the compensation parameters, the ROV was required to navigate on different heading angles and was required to change its pitching and rolling angles as much as possible, as shown in Figure 4.The resolution of the inertial sensor was 0.1°, and the dynamic measurement accuracy was 0.5°.The longitude and latitude measured by the GPS were 114.480°, and  The experimental steps are as follows.First, an experiment site with clean magnetic field was selected and the magnetic field information of this site was obtained.The ROV was controlled to change its pitch angle, roll angle and heading angle many times.At this time, the measurement results of the magnetometer mainly included the geomagnetic field and the interference magnetic field of ROV.Then, the compensation parameters were evaluated using the method in Section 3. Finally, according to the compensation parameters, the interference magnetic field of ROV was calculated and removed from the measurement results.At this time, only the geomagnetic field remained in the measurement results.Of course, if this technology is applied to the actual abnormal magnetic field detection, after removing the interference magnetic field of ROV, the remaining magnetic field is the superposition of the abnormal magnetic field of the magnetic target and the geomagnetic field.
The field experiment site was selected in a water area with a uniform magnetic field far away from artificial constructions.In order to improve the accuracy of estimating the compensation parameters, the ROV was required to navigate on different heading angles and was required to change its pitching and rolling angles as much as possible, as shown in Figure 4.The resolution of the inertial sensor was 0.1 • , and the dynamic measurement accuracy was 0.5 • .The longitude and latitude measured by the GPS were 114.480 • , and 31.099• .The average altitude of the experiment site was 50 m.Thus, the north, east and down component values calculated by the International Geomagnetic Reference Field (IGRF) were 33,475.6nT, −2779 nT and 37,322.6 nT, respectively, and the calculated total magnetic field was 50,212.7 nT.The above magnetic field component values and total magnetic field can be regarded as the reference value for estimating compensation parameters.Before the compensation experiment, the ROV magnetic field measurement system was kept at the shore for several minutes and the ambient magnetic noise was tested.The results showed that the measurement results of the cesium optical pump magnetometer fluctuated within 4.23 nT when the ROV was not working.That is to say, the magnetic noise of the environment was about 4 nT.The compensation parameters estimation process was completed in MATLAB software after obtaining the magnetic field data.magnetic field was 50,212.7 nT.The above magnetic field component values and total magnetic field can be regarded as the reference value for estimating compensation parameters.Before the compensation experiment, the ROV magnetic field measurement system was kept at the shore for several minutes and the ambient magnetic noise was tested.The results showed that the measurement results of the cesium optical pump magnetometer fluctuated within 4.23 nT when the ROV was not working.That is to say, the magnetic noise of the environment was about 4 nT.The compensation parameters estimation process was completed in MATLAB software after obtaining the magnetic field data.The measurement results of the inertial sensor are shown in Figure 5. Figure 6 shows the measurement results of the total magnetic field using the cesium optical pump magnetometer and the decomposed magnetic field component values.During compensation navigation, the thruster on the ROV is in continuous operation.The leakage magnetic field generated by the rotating propeller is reflected in the measurement results of the magnetometer in the form of a high-frequency interference magnetic field, and the different speeds of the propeller will cause the frequency of the leakage magnetic field to change.However, this kind of dynamic interference magnetic field can be removed through the low-pass filter of the data collector.Since the fluctuation amplitude of the environmental magnetic field at the experiment site was within 4.23 nT, it can be concluded that the oscillation of the output value of the optical pump magnetometer in Figure 6 was mainly caused by the attitude change of ROV.The maximum absolute error of the measured total magnetic field relative to the reference magnetic field of the IGRF model was 117.25 nT.According to the literature [34], if the geomagnetic field measurement error is greater than 30 nT, it is easy to cause the failure of underwater magnetic anomaly detection tasks.The measurement results of the inertial sensor are shown in Figure 5. Figure 6 shows the measurement results of the total magnetic field using the cesium optical pump magnetometer and the decomposed magnetic field component values.During compensation navigation, the thruster on the ROV is in continuous operation.The leakage magnetic field generated by the rotating propeller is reflected in the measurement results of the magnetometer in the form of a high-frequency interference magnetic field, and the different speeds of the propeller will cause the frequency of the leakage magnetic field to change.However, this kind of dynamic interference magnetic field can be removed through the low-pass filter of the data collector.Since the fluctuation amplitude of the environmental magnetic field at the experiment site was within 4.23 nT, it can be concluded that the oscillation of the output value of the optical pump magnetometer in Figure 6 was mainly caused by the attitude change of ROV.The maximum absolute error of the measured total magnetic field relative to the reference magnetic field of the IGRF model was 117.25 nT.According to the literature [34], if the geomagnetic field measurement error is greater than 30 nT, it is easy to cause the failure of underwater magnetic anomaly detection tasks.magnetic field was 50,212.7 nT.The above magnetic field component values and total magnetic field can be regarded as the reference value for estimating compensation parameters.Before the compensation experiment, the ROV magnetic field measurement system was kept at the shore for several minutes and the ambient magnetic noise was tested.The results showed that the measurement results of the cesium optical pump magnetometer fluctuated within 4.23 nT when the ROV was not working.That is to say, the magnetic noise of the environment was about 4 nT.The compensation parameters estimation process was completed in MATLAB software after obtaining the magnetic field data.The measurement results of the inertial sensor are shown in Figure 5. Figure 6 shows the measurement results of the total magnetic field using the cesium optical pump magnetometer and the decomposed magnetic field component values.During compensation navigation, the thruster on the ROV is in continuous operation.The leakage magnetic field generated by the rotating propeller is reflected in the measurement results of the magnetometer in the form of a high-frequency interference magnetic field, and the different speeds of the propeller will cause the frequency of the leakage magnetic field to change.However, this kind of dynamic interference magnetic field can be removed through the low-pass filter of the data collector.Since the fluctuation amplitude of the environmental magnetic field at the experiment site was within 4.23 nT, it can be concluded that the oscillation of the output value of the optical pump magnetometer in Figure 6 was mainly caused by the attitude change of ROV.The maximum absolute error of the measured total magnetic field relative to the reference magnetic field of the IGRF model was 117.25 nT.According to the literature [34], if the geomagnetic field measurement error is greater than 30 nT, it is easy to cause the failure of underwater magnetic anomaly detection tasks.

Compensation Parameters Estimation Result and Analysis
The IMOA was used to estimate the compensation parameters based on the measurement data of the optical pump magnetometer.The setting of main parameters in the IMOA are reported in Table 2.In addition, the value range of each element in the coeffi- The measurement value The value of geomagnetic field according to IGRF The three components of magnetic field (nT) The sampling points

Compensation Parameters Estimation Result and Analysis
The IMOA was used to estimate the compensation parameters based on the measurement data of the optical pump magnetometer.The setting of main parameters in the IMOA are reported in Table 2.In addition, the value range of each element in the coefficient matrix K i was [−2, 2], and the value range of each element in the vector H p was [−2000, 2000].The IMOA was compared with the stretching particle swarm optimization algorithm (SPSOA) [21] which has been widely used in the interference magnetic field compensation.In our previous research work (reference [35]), the compensation performance of MOA was compared with that of PSOA, and the results showed that PSOA has better compensation results.It can be seen from the reference [21] that SPSOA has better optimization ability than PSOA.Therefore, this paper mainly compared the compensation performance of IMOA and SPSOA.Figure 7 shows the searching process of compensation parameters as functions of iteration numbers using the IMOA.It can be seen from Figure 7 that in the early stage of the search, each compensation parameter is randomly distributed in its corresponding value space.As the iteration progresses, each mayfly eventually converges to the vicinity of the optimal solution by learning from the flight experience of the group and its own flight experience.This shows that the update and mutation mechanism of the IMOA can promote the mayfly population to fly to the vicinity of the global optimal point.This also fully proves that the algorithm has significant potential in the nonlinear optimization problems.The time costs of IMOA and SPSOA are 131.107s and 43.026 s, respectively, in the process of compensation parameters optimization.It can be seen that although the computation time consumed by IMOA is longer, it is acceptable for the interval time of compensation experiment and magnetic anomaly detection task.Moreover, with the development of computer hardware technology, the computation time will be significantly reduced.After the compensation parameters are obtained, the subsequent calculation time of ROV's interference magnetic field removal is negligible.
Figure 8 shows the descent process of the fitness function value in the iterative search process using the SPSOA and IMOA, respectively.We can determine that the fitness function values of the two optimization algorithms can decrease rapidly in the early stage of the searching process.The SPSOA is the first to reach convergence after about 1000 iterations.The IMOA still shows a steady downward trend after a rapid decline, and finally reaches convergence after 2100 iterations.In terms of convergence speed, although the IMOA is slower than SPSOA, the minimum fitness function value calculated by IMOA is smaller than that obtained by the SPSOA.It can also be seen from the fitness function value decline curve that the SPSOA is more likely to fall into the local optimum than IMOA.
The estimated compensation parameters using the two algorithms are shown in Table 3. Figure 9 compares the compensation results using the estimated compensation parameters in Table 3.As shown in Figure 9, the ROV's interference magnetic field can be eliminated remarkably from the measurement result of the magnetometer after compensation.The residual of the geomagnetic field after compensation of IMOA is reduced to less than 21 nT, in contrast to about 31 nT of the SPSOA.The results indicate that the IMOA exhibits better compensation performance and magnetic anomaly detection capability.Figure 8 shows the descent process of the fitness function value in the iterative search process using the SPSOA and IMOA, respectively.We can determine that the fitness  of the searching process.The SPSOA is the first to reach convergence after about 1000 iterations.The IMOA still shows a steady downward trend after a rapid decline, and finally reaches convergence after 2100 iterations.In terms of convergence speed, although the IMOA is slower than SPSOA, the minimum fitness function value calculated by IMOA is smaller than that obtained by the SPSOA.It can also be seen from the fitness function value decline curve that the SPSOA is more likely to fall into the local optimum than IMOA.The estimated compensation parameters using the two algorithms are shown in Table 3. Figure 9 compares the compensation results using the estimated compensation parameters in Table 3.As shown in Figure 9, the ROV's interference magnetic field can be eliminated remarkably from the measurement result of the magnetometer after compensation.The residual of the geomagnetic field after compensation of IMOA is reduced to less than 21 nT, in contrast to about 31 nT of the SPSOA.The results indicate that the IMOA exhibits better compensation performance and magnetic anomaly detection capability.Three statistic indicators, including mean value (ME), root mean square error (RMSE) and the standard deviation (STD) of the compensation residual are introduced to evaluate the compensation accuracy using different methods more intuitively.Among them, ME can evaluate the degree of clustering of residuals, RMSE can show the compensation accuracy of the algorithms and STD can display the dispersion degree of the residuals of the algorithms.The expressions of the three indicators are shown as follows, where the ri, i = 1, 2, …, n, is the deviation between the compensated geomagnetic value and the geomagnetic reference value.The statistical indicators comparisons among the two compensation methods are listed in Table 4.The statistical data show that the compensation performance using the IMOA is better than SPSOA.Three statistic indicators, including mean value (ME), root mean square error (RMSE) and the standard deviation (STD) of the compensation residual are introduced to evaluate the compensation accuracy using different methods more intuitively.Among them, ME can evaluate the degree of clustering of residuals, RMSE can show the compensation accuracy of the algorithms and STD can display the dispersion degree of the residuals of the algorithms.The expressions of the three indicators are shown as follows, where the r i , i = 1, 2, . . ., n, is the deviation between the compensated geomagnetic value and the geomagnetic reference value.The statistical indicators comparisons among the two compensation methods are listed in Table 4.The statistical data show that the compensation performance using the IMOA is better than SPSOA.

Application Verification of the Compensation Parameters
In order to further prove the effectiveness of the compensation method proposed in this paper, the verification experiment is carried out in the same place.The second group of experiment data is compensated using the estimated compensation parameters in Table 3.The obtained total magnetic field residuals after compensation is shown in Figure 10.The statistic indicators are shown in the Table 5.It can be found from the Figure 10 that the compensation parameters estimated by the two algorithms are applicable to the second group of data.However, from the prospective of statistic indicators, the IMOA method has better compensation performance than that of the SPSO.It demonstrates that the proposed compensation method proposed in this paper has good applicability for the compensation of the ROV's interference magnetic field.According to the literature [34], this paper notes that when the compensated magnetic field error is less than 30 nT, the underwater abnormal magnetic field can be detected.Since the maximum residual of the total geomagnetic field after compensation still reaches 20 nT, the magnetic field measurement data after compensation processing can identify the magnetic target with high probability when the abnormal magnetic field value of the magnetic target is greater than 20 nT.In addition, the eddy current interference magnetic field of ROV was not considered in this paper, so this method is mainly applicable to cases where the vehicle moves at low a speed, or where the vehicle attitude changes slowly.
When the ROV magnetic field measurement system detects the abnormal magnetic field in a concerned water area, the measurement results can be compensated in real time  According to the literature [34], this paper notes that when the compensated magnetic field error is less than 30 nT, the underwater abnormal magnetic field can be detected.Since the maximum residual of the total geomagnetic field after compensation still reaches 20 nT, the magnetic field measurement data after compensation processing can identify the magnetic target with high probability when the abnormal magnetic field value of the magnetic target is greater than 20 nT.In addition, the eddy current interference magnetic field of ROV was not considered in this paper, so this method is mainly applicable to cases where the vehicle moves at low a speed, or where the vehicle attitude changes slowly.
When the ROV magnetic field measurement system detects the abnormal magnetic field in a concerned water area, the measurement results can be compensated in real time according to the compensation parameters estimated in advance, so as to remove the ROV's interference magnetic field and the local geomagnetic field.When the residual magnetic field is greater than 20 nT, the local magnetic target can be identified.According to the literature [4], when the magnetometer moves in a straight line near the magnetic target, the abnormal magnetic field measured by the magnetometer shows a process of gradually increasing and then decreasing.

Figure 1 .
Figure 1.ROV's static interference magnetic field in the point P.

Figure 1 .
Figure 1.ROV's static interference magnetic field in the point P.

Figure 2 .
Figure 2. Vector relationship between total magnetic field, geomagnetic field, and disturbing magnetic field.

Figure 2 .
Figure 2.Vector relationship between total magnetic field, geomagnetic field, and disturbing magnetic field.

Figure 3 .
Figure 3.The structure layout of the ROV magnetic field measurement system.

Figure 3 .
Figure 3.The structure layout of the ROV magnetic field measurement system.

Figure 4 .
Figure 4.The compensation sailing of the ROV magnetic field measurement system.

Figure 5 .Figure 4 .
Figure 5.The measurement results of inertial navigation sensor.

Figure 4 .
Figure 4.The compensation sailing of the ROV magnetic field measurement system.

Figure 6 .
Figure 6.The measurement and decomposition results of the optical pump magnetometer: (a) the measured and true value of the total magnetic field; and (b) the decomposition results of the total magnetic field.

Figure 7 .
Figure 7.The searching process of compensation parameters: (a) the searching process of Hpx; (b) the searching process of Hpy; (c) the searching process of Hpz; (d) the searching process of K11; (e) the searching process of K12; (f) the searching process of K13; (g) the searching process of K21; (h) the searching process of K22; (i) the searching process of K23; (j) the searching process of K31; (k) the searching process of K32; (l) the searching process of K33.

Figure 7 .
Figure 7.The searching process of compensation parameters: (a) the searching process of H px ; (b) the searching process of H py ; (c) the searching process of H pz ; (d) the searching process of K 11 ; (e) the searching process of K 12 ; (f) the searching process of K 13 ; (g) the searching process of K 21 ; (h) the searching process of K 22 ; (i) the searching process of K 23 ; (j) the searching process of K 31 ; (k) the searching process of K 32 ; (l) the searching process of K 33 .

Figure 8 .
Figure 8.The descent process of the fitness function value.

Figure 9 .
Figure 9.Comparison of compensation results using different algorithms: (a) comparison of total magnetic field before and after compensation; and (b) residual comparison after compensation.

Figure 9 .
Figure 9.Comparison of compensation results using different algorithms: (a) comparison of total magnetic field before and after compensation; and (b) residual comparison after compensation.

JFigure 10 .
Figure 10.Compensation residual of the second group of experiment data: (a) comparison of total magnetic field before and after compensation; and (b) residual comparison after compensation.

Figure 10 .
Figure 10.Compensation residual of the second group of experiment data: (a) comparison of total magnetic field before and after compensation; and (b) residual comparison after compensation.

Table 1 .
Main performance specifications of the magnetometer.

Table 1 .
Main performance specifications of the magnetometer.

Table 2 .
Values of main parameters in the IMOA.

Table 3 .
Compensation parameters values estimated using the two methods.The descent process of the fitness function value.

Table 3 .
Compensation parameters values estimated using the two methods.

Table 5 .
Statistic indicators comparisons of the second group of experiment data.

Table 5 .
Statistic indicators comparisons of the second group of experiment data.