Classification of Gravity Matching Areas Using PSO-BP Neural Networks based on PCA and Satellite Altimetry Data over the Western Pacific

For inertial navigation systems (INS), as one of the major methods for underwater navigation, errors diverge over time. With the development of geophysical navigation technology, gravity navigation has become an effective method of navigation. Significant changes in the gravity characteristic of the matching region ensure that gravity matching navigation works effectively. In this paper, we combine artificial intelligence algorithms and statistical metrics to classify gravity-matching navigation regions. Firstly, this paper analyzes and extracts gravity anomaly data from a matching region in different ways. Then, a particle swarm optimization (PSO) algorithm is used to optimize the network weights of a back propagation (BP) NN. Finally, based on principal component analysis (PCA) theory and PSO-BP NN, this paper proposes the PPBA method to classify the matching area. Moreover, the Terrain Contour Matching (TERCOM) matching algorithm and gravity anomaly data from the Western Pacific are used to verify the classification performance of the PPBA method. The experiments prove that the PPBA method has a high classification accuracy, and the classification results are consistent with the matching navigation experimental results. This work can provide a reference for designing navigation regions and navigation routes for submarines.


Introduction
INS [1,2] is one of the major methods for underwater navigation. However, INS errors diverge over time. INS cannot provide accurate navigation information over a long period of time; therefore, it must be corrected through other navigation methods to suppress the growth of position errors. Underwater autonomous navigation is the key technology for the realization of effective marine underwater military strategy [3].
Geophysical navigation [4,5] has become effective in underwater autonomous navigation. At present, geophysical navigation can be aided by terrain, gravity, and geomagnetism. However, underwater terrain matching auxiliary navigation needs to transmit sounding signals, which makes it difficult to ensure the concealment of submarines [6,7]. Moreover, the geomagnetic field is a weak magnetic field [8], which makes it susceptible to interference. At the same time, it is difficult to create a high-resolution geomagnetic reference map for a large region. By contrast, the characteristics of the ocean gravity field are stable, and the matching and positioning accuracy is high [9]. The use of Earth's gravity field to realize the navigation of underwater vehicles is an important research direction for underwater autonomous navigation.
Gravity matching navigation can work effectively when the gravity characteristics of the matching navigation region change obviously, that is, the region to be matched has good matching characteristics. Therefore, the characteristic parameters of the underwater gravity field are of great significance, as they reflect the characteristics of the gravity field of the matching region. It is urgent to discern the characteristic parameters of the underwater gravity field, select the matching metric criteria, and select the methods of regional matching.
Traditional classification methods have used one of the characteristic gravity parameters such as standard deviation, latitude and longitude roughness, anomaly variation coefficient, kurtosis coefficient, or skew coefficient for analysis and comparison; the disadvantage is that this brings a certain randomness and singleness. Considering the above two shortcomings, the singleness of the parameter and no consideration of directionality, the comprehensive characteristic parameter method has been proposed by various experts and scholars. The analytic hierarchy process was proposed as a new selection criterion of gravity matching [10,11]. Li [12] proposed a new principal component-weighted average normalization method to obtain the overall characteristic parameter, which is used to divide the excellent suitability region, general suitability region, and non-suitability region; Wang studied the influence of these data on a gravity-aided INS from the perspective of a gravity anomaly reference map to evaluate the feasibility and performance of gravity navigation in different navigation regions in advance [13]. Zitova found, in research on scene matching, that a scene region suitable for the matching algorithm must be selected before realizing scene matching [14]. The correlation radius and correlation function are used to estimate the availability of the matching region effectively [15]. Fractal theory (FD) [16][17][18][19] was proposed as a comprehensive characteristic parameter method which considered the directivity of the matching region to analyze a gravity anomaly sequence and assess the gravity field characteristics for the first time. This method proved that the matching region selected by the improved method not only had larger coverage but also better directivity. However, when using a navigation sensor and a gravity anomaly reference map under current technical conditions, how to quantitatively evaluate and classify the gravity matching in different marine environments remains a problem in need of a solution.
Artificial neural networks (ANN) play an important role in the domain of intelligent control, and they have allowed for great progress in many fields, such as neuroscience, mathematics, statistics, and computer science [20][21][22]. BP NN have become one of the most widespread applications due to the fact that the BP NN algorithm is simple and plastic. The trainers in this group are mostly mathematical optimization techniques that aim to find the maximum performance (minimum error). However, there exist two disadvantages: falling into local minima easily and having a slow convergence [23,24]. Local optima refer to the sub-optimal solutions in a search space that might mistakenly be assumed as the global optimum. The challenge here is the unknown number of runs that a trainer requires to be restarted with different initial solutions to find the global optimum. To address these problems, scholars have proposed many solutions, including the Conjugate gradient [25,26], Newton [27], Gauss-Newton [28], and Levenberg-Marquardt methods. In addition, some of the most popular multi-solution trainers in the literature are: Genetic Algorithm (GA) [29,30], PSO [31][32][33], Ant Colony Optimization (ACO) [34], Artificial Bee Colony (ABC) [35,36], and Differential Evolution (DE) [37]. The reason as to why such optimization techniques have been employed as training algorithms is their high performance in terms of approximating the global optimum. At the same time, there are some recent real-life applications of recent ANN. Wang proposed an approach for source localization in underwater ocean waveguides based on machine learning that achieves a good localization performance with high SNRs [38]. Seyedali employed the recently proposed Grey Wolf Optimizer (GWO) for training Multi-Layer Perceptron (MLP), which proved that the GWO algorithm is able to provide very competitive results in terms of improved local optima avoidance and also demonstrates a high level of accuracy in the classification and approximation of the proposed trainer [39]. Abbas used a fuzzy model of control parameters of the whale optimization algorithm which was designed to train an MLP-NN and proposed fuzzy logic to adjust the Chimp optimization algorithm's control parameters to classify marine mammals. As the simulation results show, Abba's algorithms can identify the global optimal, avoid local optimization, and have better performance for classifying [40,41]. Wang proposed a new meta-heuristic method called the sine-cosine algorithm (SCA), which was used to train RBF NNs and which provided better results in terms of convergence rate and classification accuracy compared to the standard algorithms [42]. This all shows that ANN has been widely used in classification learning.
PSO is a swarm intelligence optimization method based on the swarm behavior of biota. The algorithm has the advantages of fewer parameters, rapid convergence, and strong global search ability, and it is widely used in pattern recognition, signal processing, intelligent control, and other fields. According to this idea, this paper adopts PSO to optimize a BP NN.
In this paper, we introduce an efficient PSO algorithm to optimize a BP NN and combine these with the PCA method to propose the PPBA method, which is used to classify and identify the gravity matching region in the entire research area. The PPBA method starts with a large number of characteristic gravity data and ends with the labels "suitable" or "unsuitable" by using the trained PPBA method. The experimental results for the use of TERCOM algorithms show that the navigation errors in good matching regions are much better than those in non-matching regions. The correctness of the PPBA classification results is consistent with the experimental matching navigation results. Thus, the classification evaluates the feasibility and performance of gravity matching navigation in different navigation regions in advance. This provides a reference for designing navigation regions and navigation routes for submarines.
The rest of this paper is organized as follows. In Section 2, the PPBA method is constructed. In Section 3, the experimental materials are introduced. Section 4 presents the experimental results. Section 5 is a discussion, and Section 6 draws conclusions.

Gravity Characteristic Selection Method
Significant changes in the gravity characteristic of the matching region ensure that gravity matching navigation works effectively. Better matching navigation results may be gained when gravity changes are much rougher in the matching area. Conversely, smooth gravity changes may lead to large errors and even invalid positioning results. Drawing on the definition of the statistical parameters of terrain features, the standard deviation of gravity anomaly, the standard deviation of slope, roughness, difference entropy of gravity anomaly, FD, average gravity difference, kurtosis coefficient, and relevant coefficient comprise the characteristic parameters of the gravity field.
At present, underwater gravity field data are usually stored in the form of grid data. It is assumed that the gravity data grid of a certain gravity field region is M × N, and the gravity anomaly at Grid(i,j) of the gravity field is G(i,j). In order to analyze the adaptation characteristics of the local gravity field, a local window with a size of m × n is designed to calculate the various statistical parameters of the local gravity field. Finally, the matching characteristics of each local gravity field are combined to obtain the gravity matching characteristics of the whole region. The main characteristic parameters of the gravity field are analyzed and deduced as follows.

1.
Standard deviation of gravity anomaly (σ): where g = 1 mn ∑ m i=1 ∑ n j=1 , g(i, j) is the mean of the gravity anomaly, and σ is the standard deviation of the gravity anomaly, which is used to describe the overall variation characteristics and dispersion degree of the gravity anomaly data. That is, when the gravity field features are obvious, the σ is large. On the contrary, when the σ is small, it means that the fluctuation of the gravity anomaly map is relatively flat.

2.
Standard deviation of slope (σ S ): where is the standard deviation of the slope, S is the slope of the gravity field, S ∅ and S λ are the slope of gravity field in the latitude and longitude directions, respectively, ∆Grid is the length of the grid side of the gravity field, and S is the mean of the gravity anomaly in the local calculated window. The slope of the gravity field reflects the rate of gravity change in the longitude and latitude directions, and the σ S reflects the characteristics of the rate of the slope change. The dispersion of gravity is positively correlated with the σ S , that is, the faster the slope changes, the more obvious the change in the gravity field.

3.
Roughness of gravity anomaly (R): where R is the roughness of the gravity field [43], and R ∅ and R λ are the slope of the gravity field in the latitude and longitude directions, respectively. The R is used to reflect the smoothness of the trend surface of the gravity field, depicting the fluctuation of the gravity field in this region. The anomalous divergence of gravity is positively correlated with R, that is, the greater the R, the more obvious the change in the gravity field. Difference entropy of gravity anomaly (H): where , H is the difference entropy of the gravity anomaly [44], P ij is the probability of difference entropy, and D ij is the difference in the gravity anomaly. Entropy can be used to represent the amount of information contained in an event. To a certain extent, entropy reflects the change in gravity data, and the dispersion of gravity is negatively correlated with entropy, that is, where the change in gravity data is small, the entropy value is large, and where the change is large, the entropy value is small.

4.
Relevant coefficient (RC): where RC is the relevant coefficient of the gravity field, and RC ∅ and RC λ are the relevant coefficients of the gravity field in the latitude and longitude directions, respectively. The RC of the gravity field is used to reflect the degree of the gravity anomaly, where the gravity anomaly is negatively correlated with the relevant coefficient. The larger the relevant coefficient, the higher the degree of relevance between adjacent grid points. At this time, the change of the gravity anomaly in the region is small.

6.
Fractal dimension (FD): where A(ε) is the surface area of the fractal dimension, ε is the scale to use when measuring C 1 = 2 − D, which is the slope of the fitted line, and C 0 is a constant of the fractal dimension. Peleg et al. adopted Mandelbrot's idea and extended it to surface region calculations [46]. In this extension, the map can be viewed as a gravity field surface whose gravity anomaly from the normal ground is proportional to the image's gray value. Then, all points at distance ε from the surface on both sides create a blanket of thickness, 2ε. The estimated surface region is the volume of the blanket divided by 2ε. For different values of ε, the blanket region can be calculated by A(ε) = Fε 2−D . The FD is used to reflect the complex irregularity of the shape of the local gravity field. The dispersion of gravity is positively correlated with the FD, that is, the larger the FD, the more obvious the local gravity field changes are.

7.
Kurtosis coefficient (PK): where K c1 and K c2 are the coefficients of PK. The term kurtosis was first coined by statistician; it is a measure of how flat or peaked a distribution of data is. The PK is used to describe the sharpness or flatness of the gravity anomaly distribution, which is achieved by comparison with the PK of the standard normal distribution. In this paper, the above parameters are used to describe the statistical characteristics of the gravity field from different dimensions to divide the gravity matching regions, as shown in Table 1. However, a single parameter has limitations in describing the characteristics of a gravity field, as it can only reflect the characteristics of the gravity field from a certain aspect. It is completely insufficient and inaccurate to express the matchable of a gravity matching region by any one of the above statistical parameters. Therefore, parameters of multiple dimensions are usually used for comparison and analysis in the process of actually analyzing the matching of the navigation region. At the same time, methods that combine parameters are used to determine whether a navigation region is good for matching.

The Proposed Principal-Component-Analysis-Based Method
In situations where a wide variety of gravity field statistical characteristic parameters reflect different aspects of gravity field information, how to ensure that the loss of the original gravity field information is limited to within a small range, and how to establish the selection criteria for the matching region by comprehensively using the above different gravity field characteristic parameters have become key issues.
The PCA method is one of the most widely-used data dimensionality reduction algorithms, which is used to analyze the inherent structure of data. The method helps to reduce data dimensionality by rotating coordinate axes. The proposed PCA involved an eigenvalue decomposition to produce eigenvalues and eigenvectors for representing the variation of the origin data. A set of correlated high-dimensional sensed data were transformed into a set of uncorrelated lower-dimensional data. In order to compare and analyze the size and features of the characteristic parameters, the normalization method was used to standardize the origin data to eliminate the influence of data dimensionality. This method uses calculations based on the mean and standard deviation of the original data. The original data that are processed conform to a standard normal distribution, that is, the mean is 0, and the standard deviation is 1. The steps of the PCA method can be summarized by the following: (1) Start; (2) Step (1) Select characteristic parameter X 1 , X 2 , X 3 , · · · , X k ; (3) Step (2) Standardize X and compute its covariance matrix; (4) Step (3) Calculate the eigenvalues and eigenvectors of the covariance matrix; (5) Step (4) Calculate the PCA expression and the new principal component index; (6) End.

BP NN
Since Rulmhart and the parallel distributed processing group put forward the BP algorithm in 1986 [47], ANN researchers began to pay attention to BP NN. BP networks can learn and store a large number of input-output pattern mapping relationships without revealing the mathematical equations that describe these mapping relationships in advance. Analogous to the least squares linear regression problem, when solving a data fitting line, a method is used to make the "deviation" between the predicted value and the actual value as small as possible. The BP NN updates the parameters of the network once through one forward propagation and one back propagation. The first layer is always called the input layer, whereas the last layer is called the output layer. Other layers between the input and output layers are called hidden layers. The process of backpropagation uses the principle of the gradient descent method to continuously adjust the weights and thresholds of the network to minimize the sum of the squared errors of the network and, finally, to make the network approach the real relationship. However, the flaw of the gradient descent method is that, if the location of the initial point is not selected properly, the local optimal solution is found instead of the global optimal solution. At the same time, slow learning efficiency and slow convergence are also disadvantages of BP networks.

PSO
The PSO algorithm was put forward by Kennedy and Eberhar in the United States in 1995 [48], which is a kind of evolutionary algorithm based on intelligence, and which aims to simulate the unpredictable movement of birds. The main idea is to regard a solution for the problem as a particle. The particle's process of searching for its optimal solution refers to the process of searching for food. Similar to a GA, the PSO algorithm is an iterative optimization algorithm, but without operations such as crossover and mutation. The search process is carried out by particles continuously updating their speed and position in the solution space to follow the optimal particle.
The PSO algorithm can be used to optimize the network weights of a BP NN. A traditional BP NN uses error BP to adjust the network connection weights. This method easily falls into the local optimal solution, while the PSO algorithm can search over a larger space to avoid the above problems to a certain degree. The above process is described in Figure 1.
1995 [48], which is a kind of evolutionary algorithm based on in to simulate the unpredictable movement of birds. The main id for the problem as a particle. The particle's process of searchi refers to the process of searching for food. Similar to a GA, the tive optimization algorithm, but without operations such as cro search process is carried out by particles continuously updatin in the solution space to follow the optimal particle.
The PSO algorithm can be used to optimize the network w ditional BP NN uses error BP to adjust the network connection ily falls into the local optimal solution, while the PSO algorith space to avoid the above problems to a certain degree. The abo Figure 1.
Among them, represents the d-th dimensional com tion, which refers to the i-th individual particle.
represents ponents of the optimal value in the particle swarm; Vid mensional components of particle i's speed after iterating for a t the d-th dimensional components of particle i's position after th rithm has the advantages of a smaller calculation and rapid co high-dimensional functions.  Among them, Pbest represents the d-th dimensional components of the best position, which refers to the i-th individual particle. Pgd represents the d-th dimensional components of the optimal value Gbest in the particle swarm; Vid(t) represents the d-th dimensional components of particle i's speed after iterating for a t-th time; Xid(t) represents the d-th dimensional components of particle i's position after the t-th time. The PSO algorithm has the advantages of a smaller calculation and rapid convergence in dealing with highdimensional functions.

The PCA-PSO-BP Combined Algorithm Method
Thus far, we have analyzed the advantages and the disadvantages of BP NN and the PSO algorithm. Data from these networks are typically multidimensional. To address this problem, this paper combines a PSO-BP NN with the PCA method, which can effectively reduce data dimensionality. Using the gravity field characteristic parameters and matchable labels as input, the mapping relationship between the two is obtained by the PSO-BP NN. Firstly, the gravity field characteristic parameters of the gravity matching region are used as input. Secondly, the corresponding matchable labels are obtained by the PSO-BP NN, which is trained for classification. Finally, the prediction selection of the gravity matching region is completed. The steps of the PPBA method are as follows: Step 1: PCA preprocess the original data. The gravity field characteristic parameters σ, σ S , R, H, RR, AG, FD, and PK are used as the original data of PCA, which are defined as X 1 , X 2 , X 3 , · · · , X k . Firstly, the original data are normalized to obtain the standard matrix Y. Matrix Y consists of the characteristic parameters and can be expressed as: where 2 , y i is the i-th normalized sample column vector, y ij is normalized with k raw samples (rows) on n process variables (columns), x i is the i-th original data, x i ∈ X i , x is the mean of the original data, and s is the standard deviation of the original data. Then, covariance matrix R demonstrates the correlation of the standard matrix, which can be expressed as: After the eigenvalue decomposition of R, that is, |R − λI k | = 0, k eigenvectors λ i can be calculated (i is 1, 2, · · · , p, · · · , k). The rate of information contribution can be calculated as follows: where b i is the rate of information contribution of the i-th principal components. Usually, p(p ≤ k) principal components are taken to replace the m original variables so that the cumulative contribution rate of the principal components is greater than 95%. Finally, through linear combination, k normalized gravity field characteristic parameters are recombined to generate p new parameters, Z 1 , Z 2 , · · · , Z p , which are not related to each other. The linear combination can be expressed as: where a ij is the loading of the j-th eigen parameter on the i-th principal components. After processing via the PCA method, the new gravity field characteristic parameters are sorted according to the contribution rate, where the larger the contribution rate, the more important the principal component. This paper selects the first k principal components Z 1 , Z 2 , · · · , Z k as input for the NN.
Step 2: using a gravity matching algorithm to calculate matching probability (MP). In this paper, the TERCOM algorithm is selected to carry out the gravity matching navigation experiment, which is used to calculate the MP as the index to measure the matchable. The method divides the experimental sea into regions of the same size and conducts a large number of underwater gravity matching navigation experiments to calculate the MP. Based on a large number of simulation experiments and empirical analyses, it is determined that P suit is the MP ≥ 85% in the matching area. According to the accuracy requirements of gravity matching navigation, we set labels with MP > P suit as "suitable" and MP ≤ P suit as "unsuitable." This paper will select the MP as an input for the BP NN.
Step 3: BP NN initializing. The new parameters Z 1 , Z 2 , · · · , Z p are set as the input layer of the BP NN, as determined by λ, which is the matrix's dimensionality. The number N serves as the output layer, as determined by the output properties. The number Q serves as a hidden layer; however, it has no exact mathematical theorem proving what value to take. On the basis of the empirical function Q = √ M + N + a, a = 1 ∼ 10, we combine the classification evaluation indicators to determine the number of layers of the hidden function.
Step 4: the PSO algorithm optimizes the BP NN. Each input of the NN is seen as a particle of PSO. We initialize the speed matrix and the position matrix of each particle. The position matrix will be given a random number between 0 and 1. Then, we calculate the fitness of each particle using the following formula: where K is the number of training samples, Yreal is the true value of the first i samples, and Yi is the predictive value of the first i samples. The optimal solution is the position of minimum fitness for the particle.
Step 5: compare the fitness of each particle. We obtain each particle's extreme value and the global optimal value. I f Present < Pbest, Present = Pbest, Pbest = x; else Pbest the same I f Present < Gbest, Present = Gbest, Gbest = x; else Pbest the same Among them, Present represents the current fitness of the particle, Pbest represents the best position of the current particle itself, and Gbest is the best position found by all particles in the current whole population.
Step 6: update the speed and position of each particle. After each iteration, the particle's speed (V id ) and the X id will be updated by the individual extreme value and the globe extreme value. Particles update their speed and position through the following formulas: where ω is the inertia weight using the linearly decreasing weight strategy, c 1 and c 2 are acceleration constants, and r 1 and r 2 are uniform random numbers in a range from 0 to 1.
Step 7: update the individual extreme value and the global extreme value. Comparing the individual moderate value and the globe moderate value before iteration obtains the best representation of the individual extreme value.
Step 8: stopping iteration. Through the termination condition of the algorithm, we judge whether the algorithm should stop or not. If the number of iterations reaches the maximum or an expectant error, then we return to Step 8, or else we go back to Step 4.
Step 9: use the eigenvalues in Step 1 and the adaptability labels in Step 2 as input parameters to train a reliable PSO-BP NN.
Step 10: divide the experimental region into multiple sub-regions of the same size, repeat Step 1 to obtain several independent eigenvalues, and input them into the trained PSO-BP NN to complete the classification and identification of the gravity matching region. Figure 2 depicts a flowchart of the PPBA classification system for the classification and identification of the gravity matching region.

Materials
As satellite altimetry and satellite gravity technologies continue to advance, the resolution of the marine gravity field has reached 1′ × 1′ based on data processing and further refinements. For example, the Scripps Institution of Oceanography and the University of California San Diego released the 1′ × 1′ global marine gravity anomaly model grav.img.24.1, with an accuracy of 3-8 mGal [31,[48][49]. To show and implement the PPBA classification method in a straightforward manner, we chose a publicly available collection of sea gravity anomalies as our experimental datasets [50]. This paper sourced the gravity anomaly data from the University of California San Diego website (http://topex.ucsd.edu/, accessed on 4 August 2021), with an original resolution of 1′ × 1′ global marine gravity anomaly model grav.img.29.1. The research area for this PPBA experiment that was used for gravity matching navigation was the Western Pacific region, spanning a longitude of 129° E-135° E and a latitude of 24° N-30 °N. Figure 3 is a schematic view of the gravity anomaly reference map in the Western Pacific region (the satellite imagery in figures 3 is from the General bathymetric Chart of the Oceans (GEBCO), https://www.gebco.net/, accessed on 23 June 2022). Currently, the measurement accuracy

Materials
As satellite altimetry and satellite gravity technologies continue to advance, the resolution of the marine gravity field has reached 1 × 1 based on data processing and further refinements. For example, the Scripps Institution of Oceanography and the University of California San Diego released the 1 × 1 global marine gravity anomaly model grav.img.24.1, with an accuracy of 3-8 mGal [31,48,49]. To show and implement the PPBA classification method in a straightforward manner, we chose a publicly available collection of sea gravity anomalies as our experimental datasets [50]. This paper sourced the gravity anomaly data from the University of California San Diego website (http://topex.ucsd.edu/, accessed on 4 August 2021), with an original resolution of 1 × 1 global marine gravity anomaly model grav.img.29.1. The research area for this PPBA experiment that was used for gravity matching navigation was the Western Pacific region, spanning a longitude of 129 • E-135 • E and a latitude of 24 • N-30 • N. Figure 3 is a schematic view of the gravity anomaly reference map in the Western Pacific region (the satellite imagery in Figure 3 is from the General bathymetric Chart of the Oceans (GEBCO), https://www.gebco.net/, accessed on 23 June 2022). Currently, the measurement accuracy of the latest LaCoste & Romberg-type and KSS31-type marine gravimeters reaches ±1.0 mGal under relatively poor sea conditions. Taking into consideration the influences of individual errors (i.e., zero drift correction, Eotvos correction, inaccurate location determination, and spatial correc-tion) in marine measurements, an accuracy of ±1.0 mGal − ±2.0 mGal is expected [51]. These allow for the use of the gravimeter and marine gravity reference map in underwater aided navigation. of the latest LaCoste & Romberg-type and KSS31-type marine gravimeters reaches ±1.0 mGal under relatively poor sea conditions. Taking into consideration the influences of individual errors (i.e., zero drift correction, Eotvos correction, inaccurate location determination, and spatial correction) in marine measurements, an accuracy of ±1.0 mGal − ±2.0 mGal is expected [51]. These allow for the use of the gravimeter and marine gravity reference map in underwater aided navigation. The bilinear interpolation method [52] was employed to transform the gravity anomaly data into 180 m × 180 m grid-resolution gravity data, as shown in Figure 4 (the gravity anomaly data in Figure 4 are from the University of California San Diego website (http://topex.ucsd.edu/, accessed on 4 August 2021). Figure 4 shows the gravity anomaly reference map in the test region with 180 m × 180 m resolution in two dimensions (4a) and three dimensions (4b), respectively. Statistically, the maximum value of the gravity anomaly in this region is +152. 8

Experiments and Results
In order to test the effectiveness and superiority of the PPBA method in the classification and identification of the gravity matching region, we used various characteristic The bilinear interpolation method [52] was employed to transform the gravity anomaly data into 180 m × 180 m grid-resolution gravity data, as shown in Figure 4 (the gravity anomaly data in Figure 4 are from the University of California San Diego website (http://topex.ucsd.edu/, accessed on 4 August 2021). Figure 4 shows the gravity anomaly reference map in the test region with 180 m × 180 m resolution in two dimensions (4a) and three dimensions (4b), respectively. Statistically, the maximum value of the gravity anomaly in this region is +152.8 mGal, while the minimum value is −168.3 mGal, and the average value is −0.6055 mGal. The experimental region was divided into 400 samples; the range of each sample was 15 × 15 Grid 2 , with a size of approximately 34 km × 34 km. All tests were run with Windows 10 (21H2 Version), an Intel(R) Core (TM) i5-7200U CPU @ 2.50 GHz 2.71 GHz, and Matlab2021a. of the latest LaCoste & Romberg-type and KSS31-type marine gravimeters reaches ±1.0 mGal under relatively poor sea conditions. Taking into consideration the influences of individual errors (i.e., zero drift correction, Eotvos correction, inaccurate location determination, and spatial correction) in marine measurements, an accuracy of ±1.0 mGal − ±2.0 mGal is expected [51]. These allow for the use of the gravimeter and marine gravity reference map in underwater aided navigation. The bilinear interpolation method [52] was employed to transform the gravity anomaly data into 180 m × 180 m grid-resolution gravity data, as shown in Figure 4 (the gravity anomaly data in Figure 4 are from the University of California San Diego website (http://topex.ucsd.edu/, accessed on 4 August 2021). Figure 4 shows the gravity anomaly reference map in the test region with 180 m × 180 m resolution in two dimensions (4a) and three dimensions (4b), respectively. Statistically, the maximum value of the gravity anomaly in this region is +152. 8

Experiments and Results
In order to test the effectiveness and superiority of the PPBA method in the classification and identification of the gravity matching region, we used various characteristic

Experiments and Results
In order to test the effectiveness and superiority of the PPBA method in the classification and identification of the gravity matching region, we used various characteristic parameters to qualitatively illustrate the relationship between the characteristics of the marine gravity anomaly reference map and the location accuracy of gravity matching, since different parameters may characterize properties with variable emphases. We used Formulas (1)- (12) to analyze the characteristic parameters of the 400 samples (introduced in Section 3) of the experimental region. Because of the influence of unit and dimension, we normalized the feature data using a regularization method to compare and analyze the size and characteristics of the data. Figure 5 shows a three-dimensional map and the histogram probability statistics of each characteristic parameter, standardized using 400 samples. Figure 5a1,a2-h1,h2 shows a three-dimensional map and the histogram probability statistics of the σ, σ S , R, H, FD, RC, AG, and PK, respectively. According to the correlation with gravity anomaly shown in Figure 5, it can be seen that the larger the σ, the larger the σ S , the larger the R, the larger the FD, the smaller the H, the larger the AG, the smaller the RC, the larger the PK, and the richer the characteristic information contained in the gravity field, which is a benefit for gravity matching.

Characteristic Parameters of the Gravity Field in the Experiment Region
Sensors. 2022, 11, x FOR PEER REVIEW 12 of 25 parameters to qualitatively illustrate the relationship between the characteristics of the marine gravity anomaly reference map and the location accuracy of gravity matching, since different parameters may characterize properties with variable emphases. We used Formulas 1-12 to analyze the characteristic parameters of the 400 samples (introduced in Section 3) of the experimental region. Because of the influence of unit and dimension, we normalized the feature data using a regularization method to compare and analyze the size and characteristics of the data. Figure 5 shows a three-dimensional map and the histogram probability statistics of each characteristic parameter, standardized using 400 samples. Figure 5a1,a2-h1,h2 shows a three-dimensional map and the histogram probability statistics of the σ, σ , R, H, FD, RC, AG, and PK, respectively. According to the correlation with gravity anomaly shown in Figure 5, it can be seen that the larger the σ, the larger the σ , the larger the R, the larger the FD, the smaller the H, the larger the AG, the smaller the RC, the larger the PK, and the richer the characteristic information contained in the gravity field, which is a benefit for gravity matching.  In summary, due to the limited amount of information contained in a single feature quantity, a single gravity field characteristic parameter cannot effectively evaluate matching. We used the PCA method to recalculate and analyze the eight characteristic parameters, using the specific process described in Section 2.2. Among them, the eigenvalues, variance contribution rate, and cumulative contribution rate statistics are given in Table 2. As shown in Table 2, the cumulative contribution rate of the first five principal components reaches 95%, which meets the selection criteria for principal components. The digital data with eight components was reduced to five components. We take the first five principal components as the following research objects. At the same time, we compute the load number of each characteristic parameter for the first five principal components, as shown in Table 3. As shown in Table 3, the σ, σ S , R, FD, and RC occupied the load number with a larger absolute number in the first principal component. At the same time, the H and the AG occupied the load number with a larger absolute number in the second and fourth principal components. In the third principal component, the PK occupied the load number with a larger absolute number. In the fifth principal component, the RC occupied the load number with a larger absolute number. Through the eigenvalues in Table 2 and the load number in Table 3, the principal component linear equation can be calculated as:
statistical histogram of the PCA method in the experimental region. Li [12] normalized the principal component linear combination matrix to one characteristic parameter, which easily and efficiently classifies the matching region. According to Figure 6a, the gravity characteristics of the northwest and southeast regions chang significantly. According to Figure 6b, characteristic parameters greater than 0.3 accounted for approximately 10%, characteristic parameters greater than 0 accounted for 50%, and values between −0.3 and +0.2 accounted for approximately 88%. According to Figure 6, it can be predicted that regions with characteristic parameters of the experimental region greater than 0.1 have good matchability, and the northwest and southeast parts of the experimental region are more suitable for matching navigation.

Results of the PPBA Method and the Classification Evaluation Index
The data in Equation (20) were processed by the PCA are as input data of the NN, and a PPBA method based on the PSO-optimized BP NN was proposed to perform the classification and prediction of the matching region. To test the effectiveness and superiority of the PPBA method, four types of tests were used. Tset1 is a no-feature analysis of the gravity field and a no-parameter optimization for the BP NN. Tset2 uses the PCA for the gravity field and a no-parameter optimization for the BP NN. Tset3 is a no-feature analysis of the gravity field and uses the PSO-BP NN. Tset4 uses the PCA for the gravity field and PSO-BP NN. We used a confusion matrix and average classification accuracy (AA) as the classification evaluation index. The confusion matrix is shown in Table 4. From the confusion matrix [53], more advanced classification indicators can be obtained: Accuracy, Precision-T, Precision-F, Specificity, and Sensitivity. This paper selects three of these and redefines them. The Accuracy (A), Precision-T (PT), Precision-F (PF), and average classification accuracy (AA) are as follows: 1. Accuracy (A): 2. Precision-T (PT):

Results of the PPBA Method and the Classification Evaluation Index
The data in Equation (20) were processed by the PCA are as input data of the NN, and a PPBA method based on the PSO-optimized BP NN was proposed to perform the classification and prediction of the matching region. To test the effectiveness and superiority of the PPBA method, four types of tests were used. Tset1 is a no-feature analysis of the gravity field and a no-parameter optimization for the BP NN. Tset2 uses the PCA for the gravity field and a no-parameter optimization for the BP NN. Tset3 is a no-feature analysis of the gravity field and uses the PSO-BP NN. Tset4 uses the PCA for the gravity field and PSO-BP NN. We used a confusion matrix and average classification accuracy (AA) as the classification evaluation index. The confusion matrix is shown in Table 4. From the confusion matrix [53], more advanced classification indicators can be obtained: Accuracy, Precision-T, Precision-F, Specificity, and Sensitivity. This paper selects three of these and redefines them. The Accuracy (A), Precision-T (PT), Precision-F (PF), and average classification accuracy (AA) are as follows: 1.

4.
Average classification accuracy (AA): According to the four test schemes and classification evaluation indicators, the test was carried out, and the statistical test results are shown in the Table 5. From the above tests, it was shown that the PPBA method we proposed is effective. Therefore, we randomly selected four candidate regions from among the 400 samples as Region-A, Region-B, Region-C, and Region-D for gravity matching navigation simulation analysis and recorded the MP. At the same time, we used the PPBA method to classify the four regions using the labels of "suitable" and "unsuitable." The results of the simulation and classification are shown in Table 6. As shown in Table 6, the results of matching probabilities in Region-A were 92.25% and, in Region-B, were 89.50%, and both of them were greater than 85%, which we set as the probability of success for matching navigation. The results of classification by PPBA for both Region-A and Region-B were "suitable," which is in good agreement with the result of matching probabilities. At the same time, the results of matching probabilities for both Region-C and Region-D were less than 85%, which is in good agreement with the "unsuitable" given by the classification result of PPBA. The above analysis shows that the PPBA method shows promise.

Classification Results by PPBA
To further verify the correctness of the classification results in the four regions, we used the TERCOM algorithm with the continuous optimizing method to conduct INS simulation experiments. The motion mode was uniform linear motion. The simulation parameters were set as follows: an accelerometer of 0.01 • /h with a gyro of 10 −3 m/s 2 , a sailing speed of 10 n mile·h −1 , a sailing direction of northeast at 70 • , an initial position error of 0 m, a velocity error of 0.04 m·s −1 , and a sailing direction error of 0.05 • . In addition, the real-time measured data of the gravimeter were the sum of the sampled value in the gravity anomaly reference map and random noise with a standard deviation of 1 mGal, while the number of sampling points and the sampling period were set as 110 and 2 minutes, respectively. The total sailing time was 3.6 hours. The matching results of the inertial navigation simulation experiment are shown in Figure 7, and the error statistics are shown in the Table 7.
the real-time measured data of the gravimeter were the sum of the sampled value in the gravity anomaly reference map and random noise with a standard deviation of 1 mGal, while the number of sampling points and the sampling period were set as 110 and 2 minutes, respectively. The total sailing time was 3.6 hours. The matching results of the inertial navigation simulation experiment are shown in Figure 7, and the error statistics are shown in the Table 7.    According to Table 7, the latitudinal, longitudinal, and radial errors in regions A and B were much smaller than in regions C and D. Therefore, it is further verified that, compared with regions C and D, regions A and B had good matching performance. Furthermore, the results are consistent with the classification of PPBA. Finally, when the gravity anomaly data of the entire experimental region were input into the PPBA system, we obtained the matching classification map of the experimental region. The gravity matching navigation matchable region in the Western Pacific is shown in Figure 8. jectory (red line) and the matching trajectory (yellow line) had a good match in regions with significant changes in gravity anomalies. At the same time, Figure 7c−d shows that the real trajectory (red line) and the matching trajectory (yellow line) had a large deviation in regions with smooth changes in gravity anomalies. According to Table 7, the latitudinal, longitudinal, and radial errors in regions A and B were much smaller than in regions C and D. Therefore, it is further verified that, compared with regions C and D, regions A and B had good matching performance. Furthermore, the results are consistent with the classification of PPBA. Finally, when the gravity anomaly data of the entire experimental region were input into the PPBA system, we obtained the matching classification map of the experimental region. The gravity matching navigation matchable region in the Western Pacific is shown in Figure 8.

Discussion
As seen in Figure 8, the white portion is the suitability region of the matching region, and the black portion is the unsuitability region. The suitable regions are mainly concentrated in the northwest and southeast, which is consistent with our previous prediction. This provides a reference for planning navigation regions and designing navigation routes for underwater submarines.

Discussion
As seen in Figure 8, the white portion is the suitability region of the matching region, and the black portion is the unsuitability region. The suitable regions are mainly concentrated in the northwest and southeast, which is consistent with our previous prediction. This provides a reference for planning navigation regions and designing navigation routes for underwater submarines.
The gravity characteristics of the matching area change significantly, which causes the gravity matching navigation to work effectively. We proposed the PPBA method to extract the gravity characteristics of the matching area for analysis and then input them into the NN to realize the classification of the area. In order to test the effectiveness and superiority of the PPBA method in the classification and prediction of the matching region, we designed four tests. As seen in Table 5, A and AA of test1 were the lowest at 87.00% and 87.60%, respectively. Compared with test1, test2 and test3 performed better in terms of A and AA, where A increased by 4.25% and 6.50%, and AA increased by 3.14% and 4.38%, respectively. Therefore, using PCA for the gravity anomaly and the BP NN using PSO to optimize the parameter was shown to be effective. At the same time, test4, using the PPBA method, yielded the best performance, where A was 96.75%, AA was 96.22%, PT was 92.86%, and PF was 99.57%. According to other studies [45], many traditional methods use one of