A co‐occurrence matrix‐based matching area selection algorithm for underwater gravity‐aided inertial navigation

Key Technologies Research and Development Program, Grant/Award Number: 2016YFB0501700; National Natural Science Foundation of China, Grant/Award Number: 61673060 Abstract The matching area selection algorithm is one of the key technologies for underwater gravity‐aided inertial navigation system, which directly affects the positioning accuracy and matching rate of underwater navigation. The traditional matching area selection algorithms usually use the statistical characteristic parameters of gravity field. However, the traditional algorithms are difficult to reflect the spatial relation characteristic of gravity field, which always miss some latent matching areas with obvious change of gravity field. In order to solve this problem, the matching area selection algorithm based on co‐ occurrence matrix is proposed. The proposed algorithm establishes gravity anomaly co‐ occurrence matrix and extracts spatial relation characteristic parameters to reflect the gravity field. The comprehensive spatial characteristic parameter is built by entropy and is used to select the matching area by maximization of inter‐class variance. The experimental results show that the proposed algorithm can select more effective matching areas than the traditional algorithms.


| INTRODUCTION
Inertial navigation system (INS) is widely used for underwater vehicles [1][2][3]. However, the INS errors diverge with time accumulation, which extremely affect the navigation accuracy. Therefore, INS cannot be applied to long-endurance underwater navigation independently [4,5]. INS errors need to be corrected by other navigation methods [6][7][8][9][10][11]. The gravity field is a stable geophysical field. INS errors can be corrected by gravity matching. In recent years, with the development of satellite altimetry and gravity measurement technology, the resolution and accuracy of underwater gravity field data have been improved [12][13][14]. Thus, the gravity-aided INS is suitable for long-endurance and high-precision autonomous underwater navigation, which meets the requirements of energy economy, detection continuity and safety of underwater vehicles [15,16].
The matching area selection algorithm is one of the key technologies for underwater gravity-aided INS [17]. The gravity matching of underwater vehicles can be only completed in the matching area. In gravity-aided INS, real-time gravity anomaly collected by gravity sensor and gravity anomaly from gravity reference database are compared by the matching algorithm. The precise matching position information can be used to correct the position errors of INS. In the process of gravity matching, the matching accuracy is closely related to the characteristic change of gravity field. The more obvious characteristic change of gravity field is, the higher the matching accuracy is. In the area with no obvious characteristic change of gravity field, the characteristics of gravity field around the INS indication position are similar, which easily lead to mismatch. Therefore, the matching area selection directly affects the matching accuracy and matching rate of gravity matching.
The traditional matching area selection algorithms mainly use statistical characteristic parameters of gravity field. The data structure of gravity field is grid. Each grid has longitude, latitude and gravity anomaly. The traditional statistical characteristic parameters include standard deviation, the latitude and longitude correlation, variation coefficient, kurtosis coefficient, skew coefficient etc. The traditional algorithms select the matching area by comparing and analysing the traditional statistical characteristic parameters [18][19][20]. However, the shortcomings of the traditional algorithms are uncertainty of parameter threshold and ill consideration in spatial characteristics of gravity field. Due to the traditional data structure of gravity field, there is no statistical characteristic parameter that can measure the spatial characteristics of gravity field in the traditional algorithms. Therefore, the traditional algorithms are easy to miss the potential matching area so as to affect the accuracy of gravity matching. Recently, the matching area selection algorithm based on isotropic coefficient has been proposed [17]. It combines the traditional statistical characteristic parameters and isotropic coefficient to make the selection more accurate. Isotropic coefficient is used to evaluate the directional characteristic of gravity field. However, the isotropic coefficient algorithm is still limited. It is also difficult to analyse more spatial characteristics of gravity field. Therefore, we should find more spatial characteristics to select more potential matching areas.
In order to address the limitations of the existing algorithms, a matching area selection algorithm based on co-occurrence matrix is proposed. The gravity anomaly cooccurrence matrix (GACM) is built, which is a second-order function. GACM is used to measure spatial relationship of sampling points of gravity field, which can better reflect the characteristic change of gravity field. Based on GACM, the spatial relation feature parameters are calculated. They are used to constitute the comprehensive spatial characteristic parameter of gravity field by entropy method. The comprehensive spatial characteristic parameter can better analyse the spatial characteristic of gravity field and is used to automatically determine the matching area criteria by maximization of interclass variance method. The experimental results show that the matching area selection by the proposed algorithm is more valid than traditional algorithm and the isotropic coefficient algorithm.
The rest of this paper is organized as follows. Section 2 introduces the principle of gravity-aided INS and the grey level co-occurrence matrix. In Section 3, the matching area selection algorithm based on co-occurrence matrix is presented. Section 4 analyses and compares the experimental results of different algorithms in simulation experiments and sea experiments. Finally, Section 5 summarizes this work and offers the future work on this research topic.

| Gravity-aided inertial navigation system
Gravity-aided inertial navigation is a passive autonomous navigation technology. As shown in Figure 1, the gravity-aided INS consists of four parts: INS, gravity anomaly reference database, gravity sensor and the matching algorithm.
The gravity matching can only be completed in the matching area. When underwater vehicles sail to the matching area, the initial position provided by INS is helpful to search for gravity anomaly information from the gravity reference database. Meanwhile, the gravity sensor offers a real-time gravity anomaly measurement information. Then both the gravity anomaly information is sent to the matching algorithm.
The precise matching position information is obtained by matching algorithm to correct the INS position.
In the process of gravity matching, the matching accuracy is closely related to the characteristic changes of gravity field. In the area where the gravity field does not change much, the characteristic of gravity field around the INS indication position is similar to those measured by gravity sensor. In this case, it is difficult to gravity matching. Therefore, matching area selection is the first step of gravity matching.

| Grey-level co-occurrence matrix
Grey-level co-occurrence matrix (GLCM) is also known as the grey-level spatial dependence matrix, which is a statistical method of examining texture [21]. GLCM is based on the estimation of the second-order joint conditional probability density functions. GLCM describes pixel spatial relationship and texture feature of image by calculating how many times pairs of pixels with specific values and in a specified spatial relationship occur in an image.
As shown in Figure 2, GLCM counts the times of two pixels occurred in direction 0°and distance 1 in an image. For example, the times of two pixels (1, 1) occurred in direction 0°a nd distance 1 is 1, and the times of two pixels (1, 2) occurred in direction 0°and distance 1 is 2. The spatial relation feature parameters of the image are calculated from the GLCM.
The data structure of gravity field is similar to the raster image. The gravity field data are stored in a grid. Each grid has longitude, latitude and gravity anomaly. Therefore, as a method of image feature analysis and classification, GLCM can be extended to analyse the characteristic of gravity field and select the matching area.

| THE MATCHING AREA SELECTION ALGORITHM BASED ON CO-OCCURRENCE MATRIX
The flow chart of the matching area selection algorithm based on co-occurrence matrix is shown in Figure 3.
Firstly, the gravity field data are transformed to a certain format. The latitude and longitude coordinates of gravity field are transformed into the grid coordinates. The gravity anomaly is transformed into the gravity anomaly levels. Then, the GACMs are created. The spatial relation feature parameters are calculated from GACMs in different distances and directions. The comprehensive spatial characteristic parameter of the gravity field is built by the entropy method. Finally, the comprehensive spatial characteristic parameter is used to select the matching area by the maximization of inter-class variance method.

| Gravity anomaly co-occurrence matrix
GACM represents the frequency of a sampling point with the gravity anomaly level i and another sampling point with the WANG ET AL.
-251 gravity anomaly level j occurring in a specific spatial relationship. In the proposed algorithm, the gravity anomaly levels N g are defined as 1-16, which determines the size of GACM. As shown in Figure 4, the sampling point relationships of direction and distance are defined. Four directions are horizontal 0°, diagonal 45°, vertical 90°and diagonal 135°. The distance is 1. In this case, the gravity field is represented by 16 GACMs. The GACM is established as follows:

| Great GACM
A single GACM is not enough to describe the spatial characteristic of gravity field. Therefore, multiple GACMs in different directions and distances are created. The equation of GACM is where i and j are gravity anomaly levels of two elements, d is the distance between two elements, and θ is the direction between two elements. where R is normalized parameter.

| Calculate statistical properties of GACM
GACM can reveal the spatial relation characteristic of gravity field. We calculate the spatial relation feature parameters in four directions by the GACMs in four directions, respectively, and then take the average value.

| Spatial relation feature of gravity field
The spatial relation features of gravity field are helpful to the classification of gravity field and the selection of the matching area. The spatial relation feature parameters are calculated from GACMs. The four parameters energy, contrast, correlation and homogeneity are unrelated, which can be used to analyse and reflect the spatial characteristic of gravity field.

| Energy
Energy represents the sum of squared elements in GACM. Energy is also known as uniformity or the angular second moment. The smaller the energy, the more obvious the change of gravity field. The equation of energy is where M is the number of rows and N is the number of columns.

| Contrast
Contrast measures the intensity contrast between a sampling point and its adjacent sampling points over the whole gravity field, and reflects the local variations in GACM. The larger the contrast, the more obvious the change of gravity field. The equation of contrast is where N g is the number of the gravity anomaly levels, M is the number of rows and N is the number of columns.

| Correlation
Correlation measures the joint probability occurrence of the specified sampling point pairs, and reflects how correlated a sampling point is to its adjacent sampling points over the whole gravity field. The equation of correlation is where μx and μ y are means of p x and p y , respectively; σ x and σ y are standard deviation of p x and p y , respectively; M is the number of rows and N is the number of columns.

| Homogeneity
Homogeneity measures the closeness of the distribution of elements in GACM to the GACM diagonal. The smaller the homogeneity is, the more obvious the change of gravity field is. The equation of homogeneity is where M is the number of rows and N is the number of columns.

F I G U R E 4
The spatial relationships in gravity anomaly co-occurrence matrix WANG ET AL. -253

| Comprehensive spatial characteristic parameter based on entropy
The above four spatial relation feature parameters are used to construct the comprehensive spatial characteristic parameter. The construction process by the entropy method is as follows: The evaluation matrix A is composed of four evaluation indices which include energy, contrast, correlation and homogeneity. The equation is where m = 4 is the number of evaluation indices, n ¼ M � N is the number of evaluation objects and the size of the feature matrix is M � N.
Because each spatial relation feature parameter is not in the same quantity level, the evaluation matrix A is needed to normalize. Energy and homogeneity are negative indices, so the equation is Contrast and correlation are positive indices; therefore, the equation is The normalized evaluation matrix is

| Calculation of entropy weight
The entropy value of the evaluation index i in A ′ is calculated as follows: Then, the entropy weight of each spatial relation feature parameter is calculated as follows: 3.3.3 | Comprehensive spatial characteristic parameter The comprehensive spatial characteristic parameter T is built by calculating the weighted summation of the four spatial relation feature parameters:

| Matching area selection criteria based on maximization of inter-class variance
The matching area selection criterion is determined by maximization of inter-class variance method. The assumptions are as follows: the grid size of gravity field is M � N; the grid number of the matching area and the non-matching area are N 0 and N 1 ; the proportion of the matching area and the non-matching area are k 0 and k 1 ; the average value of the comprehensive spatial characteristic parameter in the matching area and the nonmatching area are t 0 and t 1 . The average value of the comprehensive spatial characteristic parameter is t, which is given as follows: where 8 > > < > > : The between-cluster variance of the matching area and the non-matching area is σ: The maximization of inter-class variance was obtained by the traversal method. The adaptive optimal threshold t * maximizes the between-cluster variance σ. The equation is Therefore, the comprehensive spatial characteristic parameter is used to reflect the spatial characteristic of gravity field. The gravity field is divided into two types: matching area and non-matching area. The matching area selection criterion is T ≥ t * .

| The schema of the proposed algorithm
The scheme of the proposed algorithm is as follows: Step 1. Transform gravity field data. The latitude and longitude coordinates of gravity field are transformed to the grid coordinates. The gravity anomaly value is transformed to the gravity anomaly levels.
Step 2. Establish the GACM. The transformed gravity field data is used to establish the GACMs in different directions and distances.
Step 3. Calculate comprehensive spatial characteristic parameter. The four spatial relation feature parameters are calculated from GACMs. The four spatial relation feature parameters are used to build the comprehensive spatial characteristic parameter by the entropy method.
Step 4. Select the matching area. The comprehensive spatial characteristic parameter describes the spatial characteristic of gravity field. The matching area selection criteria is determined by the maximization of inter-class variance method. Then the gravity field is divided into the matching area and the non-matching area.

| The advantages of the proposed algorithm
Compared with the traditional algorithm and the isotropic coefficient algorithm, the main contributions of the proposed algorithm are as follows: First, the traditional algorithm and the isotropic coefficient algorithm are first-order functions which measure the statistical characteristics of gravity field. However, they do not consider spatial relationship of sampling points of gravity field. The proposed algorithm is second-order function, which considers the spatial relationship of sampling points of gravity field.
Second, GACM is a function of the spatial variation in gravity anomaly intensities. If the characteristic of gravity field is more spatial characteristic than intensity characteristic, spatial relation analysis is helpful, but the traditional thresholding method is not effective. Spatial relation analysis attempts to describe gravity field by intuitive quantification, such as rough, smooth, silky or bumpy. The traditional algorithms cannot analyse the spatial relationships of sampling points of gravity field.
Finally, the proposed algorithm applies the entropy method and the maximization of inter-class variance method to determine the matching area selection criterion, which reduces the manual intervention and artificial experience judgement. The adaptive optimal threshold is automatically found by the proposed algorithm.

| Simulation experiments
The semi-physical simulation experiment of the gravity-aided INS is used to demonstrate the effectiveness of the proposed algorithm. Real INS data and real gravity field data are used in simulation experiments. The semi-physical simulation platform is shown in Figure 5.

| Simulation experiment conditions
The gravity field data is as follows: the range is northern latitude 25°-30°and east longitude 125°-130°; the resolution is 1´�1´; the grid of gravity field is 300�300; the maximum gravity anomaly is 91.09 mGal; the minimum gravity anomaly is −75.97 mGal. The real-time gravity field measurement data is obtained by adding white Gaussian noise whose variance is 2 mGal to the gravity field data. The grey image of the gravity field is shown in Figure 6. According to the contour distribution of gravity field, it can be seen that there are obvious areas, gentle areas and fluctuating areas of gravity anomaly change. Therefore, the gravity field data with the typical characteristics is suitable for simulation experiment to verify the validity of the proposed algorithm.
The dual-axis rotational fibre-optic INS is used to output INS trajectory. The INS errors are shown in Table 1.

| Simulation of matching area selection algorithms
In order to make comparison among the traditional algorithm, the isotropic coefficient algorithm and the proposed algorithm, the same gravity field is used in simulation.
The principle and process of the traditional algorithm and the isotropic coefficient algorithm are detailed in [14]. The matching areas selected by the two algorithms are shown in Figures 7 and 8, respectively. The yellow areas are the matching areas. The matching areas 1-3 are selected by the two algorithms.
The comprehensive spatial characteristic parameter describes the gravity field. According to the adaptive optimal threshold of the comprehensive spatial characteristic parameter, the gravity field is divided into matching area and non-matching area. The matching areas selected by the proposed algorithm are shown in Figure 9. The yellow areas are the matching area and the other areas are the non-matching area. If the yellow line is not large enough and dense, it is not suitable to gravity matching for underwater vehicle. If the yellow line is large enough and dense, it is matching area. The matching areas 1-4 are selected by the proposed algorithm.

| Comparison results of different algorithms
The simulation results show that the matching areas 1-3 are common areas selected by the tradition algorithm, the isotropic coefficient algorithm and the proposed algorithm at the same time. The scope of matching areas 1-3 selected by the proposed algorithm is larger, and spatial continuity is better. The reason is that GACM considers the spatial distribution relationship of gravity field in the proposed algorithm, but the traditional algorithm and the isotropic coefficient algorithm only consider the discrete statistics characteristic. The WANG ET AL.
-255 matching area 4 is only selected by the proposed algorithm. The matching area 4 with great characteristic change of gravity field is missed by the traditional algorithm and the isotropic coefficient algorithm. The reason is that the spatial relation characteristic parameters of the proposed algorithm can find more spatial characteristic of gravity field than the statistical parameters. There are no parameters which can measure the spatial relationship characteristics of gravity field in the traditional algorithm and the isotropic coefficient algorithm. In addition, the threshold is very critical to select the matching area. The threshold selection needs to experience judgement and repeat in the traditional algorithm and the isotropic coefficient algorithm. Different experimental areas have different thresholds. However, the threshold can be selected dynamically and autonomously by the maximization of inter-class variance method in the proposed algorithm. Therefore, the traditional algorithm and the isotropic coefficient algorithm may miss some useful matching area. By the contrast, the matching area selected by the proposed algorithm has more continuity and the coverage is larger. Simulation results show that the proposed algorithm is more effective than the traditional algorithm and the isotropic coefficient algorithm.

| Simulation experiment of gravity matching
The simulation of gravity matching is applied to demonstrate the effectiveness of the matching area selected by the proposed algorithm. As shown in Figure 9, the matching area 4 is only selected by the proposed algorithm, while it is omitted by the traditional algorithm and the isotropic coefficient algorithm. Gravity matching navigation experiment is carried out in the matching area 4.
In the gravity matching simulation experiment, the navigation trajectories are suitable for the shape and range of the matching area. The situations of gravity matching navigation are as followed: the navigation distance is 101.82 km; the directions of navigation are four directions (0°, 45°, 90°and 135°); the maximum positioning error is about 10.9 km; the matching algorithms are the mean absolute difference algorithm and the mean squared difference algorithm; matching points are 41; matching times are 32. In the gravity field with resolution 1´�1′, if the distance between the matching trajectory and the real trajectory does not exceed one 1´�1′ grid in longitudinal and latitudinal directions, it indicates that the gravity matching is successful. Namely, the positioning accuracy of navigation system is about 2.42 km. The matching rate is the proportion of the number of successful matches in the total number of matches.
The results of matching error and matching rate are shown in Table 2. The navigation trajectories are suitable for the shape and range of the matching area. The experimental results show that the gravity matching rate in the matching area 4 is above 95%, which indicates that the matching area only selected by the proposed algorithm is effective. The traditional algorithm and the isotropic coefficient algorithm omit the matching area 4. Hence, the proposed algorithm is more valid than the traditional algorithm and the isotropic coefficient algorithm.

| Sea experiment
In order to further prove the effectiveness of the proposed algorithm, the sea experiment is carried out on scientific survey ship in a certain sea area. The sea experiment is based on the high-precision laser INS, the gravimeter, the gravity matching device and the global navigation satellite system (GNSS). The differential positioning trajectory of GNSS represents the real trajectory of the ship.

| Sea experiments conditions
The gravity field data is obtained by fusion processing of multisource data, including satellite altimetry inversion gravity, shipborne gravity and airborne gravity. The gravity field data is as follows: the range is northern latitude 20°-21°and east longitude 117.5°-118.5°; the resolution is 1´�1´; the accuracy is ±5 mGal; the maximum value is 68.755 mGal; the minimum value is −16.362 mGal. The mean value is 7.3857 mGal. The gravity field contour map is shown in Figure 10. The real-time measurement gravity anomaly is obtained by gravimeter. The measurement accuracy of gravimeter is ±3 mGal.

| Result of matching area selection algorithms
The matching areas selected by the traditional algorithm, the isotropic coefficient algorithm and the proposed algorithm are shown in Figures 11-13, respectively. The yellow areas are the area with strong matching. -257 The yellow areas selected by the traditional algorithm and the isotropic coefficient algorithm are small and discontinuous, which are not suitable for the matching area. The matching areas 1 and 2 selected by the proposed algorithm are large and continuous, which are suitable for gravity matching navigation experiment.
Through the analysis and comparison of the results, it is found that the traditional algorithm and the isotropic coefficient algorithm did not select the matching area. This is because these two algorithms use the traditional statistical characteristic parameters of gravity field which cannot reflect the matching characteristic of this area. However, the proposed algorithm selected the matching areas 1 and 2. This is because the proposed algorithm uses the spatial relation characteristic parameters of gravity field which can reflect the matching characteristic of this area. Therefore, the proposed algorithm is more effective than the traditional algorithm and the isotropic coefficient algorithm.

| Experiment of gravity matching navigation
The experiment of gravity matching navigation is applied to demonstrate the effectiveness of the matching area selected by the proposed algorithm.
We have planned a route for gravity matching navigation experiment in matching area 1. As shown in Figure 14, the ship sailed in the direction of the red arrow. The black line is the real trajectory of the ship, the blue line is the INS trajectory and the red line is the matching trajectory. The matching trajectory is closer to real trajectory. It indicates that the gravity matching is successful in the matching area 1.
The error statistical results of INS trajectory and the matching trajectory are shown in Table 3. The longitude error, F I G U R E 1 0 Gravity field contour map F I G U R E 1 1 Matching areas selected by the traditional algorithm F I G U R E 1 2 Matching areas selected by the isotropic coefficient algorithm F I G U R E 1 3 Matching areas selected by the proposed algorithm 258latitude error and position error of the matching trajectory are obviously lower than that of INS trajectory. Compared with INS errors, the mean of absolute error in latitude direction is reduced by 43.54%; the standard deviation of error in latitude direction is reduced by 49.51%; the mean of absolute error in longitude direction is reduced by 35.21%; the standard deviation of error in longitude direction is reduced by 37.76%; the mean position error is reduced by 41.11%; the standard deviation position error is reduced by 17.61%. It can be seen from the error statistics that the matching errors are reduced and the matching accuracy is improved in the matching area 1.
Through simulation experiment and sea experiment, it has been proved that the matching areas selected by the proposed algorithm is effective. Therefore, the proposed algorithm is more effective and superior than the traditional algorithm and the isotropic coefficient algorithm.

| CONCLUSION
A matching area selection algorithm based on co-occurrence matrix is proposed in this study. In the proposed algorithm, two major improvements are as followed: (1) the spatial relation characteristics are used in the matching area selection, which compensate the shortage of the other algorithms to describe the singularity of gravity field characteristic. (2) The entropy method and maximization of inter-class variance method are used to determine the matching area selection criteria, which reduce the manual intervention and artificial experience judgement. The experimental results show that the proposed algorithm is more effective than the traditional algorithm and the isotropic coefficient algorithm.
The matching area selection algorithm aims to improve the matching rate and the matching accuracy of gravity matching. Therefore, in the future work, the spatial relation characteristics of gravity field will be applied to gravity matching algorithm and to further improve the accuracy of matching algorithm. -259