Three-Dimensional Underwater Path Planning of Submarine Considering the Real Marine Environment

The submarine is an underwater ship that can perform a variety of combat missions in a complex marine environment. Path planning of submarines has always been the focus of military marine engineering research. In most practical applications, there are numerous marine physical phenomena in the marine environment, such as pycnocline, density fronts, mesoscale eddy, etc., which have an important impact on the navigation of submarines. First, the artificial potential field heuristic factor is introduced into the ant colony algorithm to improve its convergence speed, and the artificial potential field ant colony optimization (APF-ACO) is obtained. In addition, this article uses the unit-body to reflect the regional physical elements and quantifies the physical marine phenomenon in the form of the cost function, which is used to solve the problem of submarine path planning in the complex marine environment. In this article, the algorithm is tested in a real marine data environment. The experimental results show that the algorithm can realize the utilization of the ocean sound speed environment, ocean density environment and ocean current environment, and obtain a path more suitable for submarine underwater navigation.


I. INTRODUCTION
The submarine is a ship that can fight independently underwater. It appeared in the First World War and was widely used in various military operations afterwards. It has been made an important part of modern naval operations, and mainly undertakes the important strategic role of combat readiness, cruise and deterrence against the enemy.
With the continuous changes in the international situation, underwater navigation of submarines has become a hot research topic. Unlike surface ships, submarines face more unknown threats underwater, mainly including: seabed topography and bottom quality, marine physical phenomena, enemy sonar detection, etc [1].
In order to ensure safe navigation, submarines need to choose the path that best meets the need of combat missions among the many navigation paths. Therefore, path planning is a vital step for submarine underwater operations. Under The associate editor coordinating the review of this manuscript and approving it for publication was Heng Wang . normal circumstances, three factors need to be considered in submarine path planning: 1) Concealment factors. The combat effectiveness and threat of submarines are mainly focused on their concealment characteristics. Therefore, when performing combat missions, submarines should ensure their concealment as much as possible. Mainly consider the characteristics of marine acoustic velocity distribution and other characteristics to plan the navigation path. 2) Safety factors. The submarine's underwater navigation environment includes the enemy's sonar detection range, dangerous physical marine phenomena, seabed topography and submarine reefs and other factors that threaten the submarine's survival. Therefore, these areas should be prevented in submarine path planning to achieve the goal of safe navigation. 3) Economic factors. The submarine is difficult to replenish underwater, so it is necessary to consider the submarine's endurance when performing tasks. The usual method is to plan the shortest path for the submarine from the departure to the destination and try to use marine currents and other marine elements to save fuel [2]. Due to the complexity and nonlinearity of the underwater 3D environment, the 2D land-based algorithms cannot be directly used to solve the underwater path planning problems. In response to this issue, many scholars have focused on the development of underwater path planning algorithms. The most commonly used path planning algorithms are as follows: Geometric model search algorithms: The most representative algorithms are Dijsktra algorithm and A * algorithm. Since Dijkstra algorithm traverses all nodes and obtains the shortest path, the obtained shortest path has a high success rate and satisfactory robustness [3], [4]. However, when applied to large-scale complex path topology networks, node traversal and low efficiency are fatal drawbacks [5]. Arinaga et al. applied Dijkstra algorithm to a global path search for AUVs in an underwater environment. The algorithm can avoid a set of obstacles and reach the end but does not consider the impact of the marine environment [3]. The A * algorithm is based on Dijkstra algorithm to add the estimated cost of the target point to the current node [6]. However, the A * algorithm gradually determines the next path grid by comparing the heuristic function values of the adjacent grids of the current path. The A * algorithm is not guaranteed to be optimal when there are multiple minima [7]. Garau et al. implemented a path search with A * algorithm and considered the influence of marine environmental factors [8].
Sample planning-based approaches: The most representative algorithms are PRM algorithm and RRT * algorithm. At the heart of the PRM algorithm is sampling and building a roadmap. In [9] and [10], the algorithm was verified to be applicable to underwater 3D path planning under complex constraints. RRT * algorithm has a powerful spatial search ability [11]. For single query problems, the RRT * algorithm has a faster convergence speed than the PRM algorithm [12]. Carreras et al. employed the RRT * algorithm to perform 2D AUV path planning, and the 3D results show that the adaptability of this method in the real complex environment is Satisfactory [11] Artificial potential field (APF)-based approaches: Solari et al. introduced the artificial potential field algorithm into underwater path planning and found that it is feasible in a static environment, but it cannot perform accurate dynamic obstacle avoidance [13]. To solve this problem, Cheng et al. added the speed synthesis algorithm of AUV to the APF algorithm, which improved the convergence speed of the algorithm while avoiding obstacles accurately [14]. Jantapremjit et al. not only realize automatic obstacle avoidance by applying the APF algorithm, but also introduced the state-dependent Riccati equation method to optimize the optimal high-order sliding mode control, which improved the robustness of the AUV motion [15].
Evolutionary algorithms (EAs)-based approaches: EDA is a stochastic optimization algorithm based on statistical principles. Compared with other evolutionary algorithms, EDA has stronger global search ability and faster convergence speed. But the algorithm is prone to fall into premature convergence [16]. The PSO algorithm performs path planning operations by simulating the individual cooperation mechanism in the group [17]. It has the advantages of easy implementation and rapid convergence, but it is easy to fall into the local optimal solution [18]. Liu et al. used the PSO algorithm for AUV path planning. Simulation experiments show that the algorithm is simple easy to implement, not sensitive to the population size, and has a faster convergence speed [17]. GA is a computational model that simulates genetic selection and natural elimination. Similar to PSO algorithm, GA also finds the optimal solution through iteration of random solutions [19]. Compared with the PSO algorithm, GA can find the global optimal solution according to the adaptive evaluation, but the calculation time is too long [20]. The difference between DE algorithm and GA is that DE algorithm uses the difference vector between individuals in the population to achieve individual variation. The robustness of DE algorithm is better than that of genetic algorithm [21].
Heuristic algorithms (HAs)-based approaches: The most representative algorithms are ACO and SA. ACO is a probabilistic algorithm with the characteristics of distributed computing, positive feedback and heuristic search [22]. The advantage of ACO is that it can be applied to the underwater 3D path search problem, and its parameters are relatively small and do not need to be adjusted manually [23]. This is why this paper chooses the ant colony algorithm as the basic algorithm of the optimization algorithm. However, the convergence speed of ACO is slow, and it is easy to fall into the local optimum, so ACO needs to be optimized. Ma et al. introduced alarm pheromone in the ACO algorithm (AP-ACO) for path planning of underwater vehicles. The experimental results show that compared with the ordinary ACO algorithm, AP-ACO has faster convergence speed and stability [24]. However, the above algorithms only consider the impact of marine currents on underwater path planning, which is obviously inconsistent with the real marine environment. SA performs path planning by simulating the annealing process of solid matter. SA has the characteristics of flexible use and less initial conditions [25]. [26] used SA for underwater path planning. The results show that the algorithm can complete the large turning radius path planning of AUV, but it has the shortcomings of slow convergence speed and poor randomness.
In addition to the above algorithms, in recent years, many emerging hybrid algorithms have been proposed and proven to be effective. For example, Orozco-Rosas et al. integrated membrane computing, the pseudobacterial genetic algorithm and artificial potential field algorithm to obtain a hybrid path planning algorithm based on membrane pseudobacterial potential field (MemPBPF). The algorithm not only realizes the effectiveness of the path planning of the mobile robot, but also can realize the parallel operation on the multi-core computer, which improves the efficiency [27], [28].
However, some of the above algorithms are not suitable for underwater 3D path planning, and some only consider the impact of ocean currents on path planning, which is obviously inconsistent with the real marine environment.
Since various marine phenomena have important effects on underwater navigation in the real marine environment, the impact of the marine environment on the planning results should be taken into account when planning the path. But the research in this area is not sufficient. Yu et al. proposed a fast-stepping method in a large three-dimensional marine battlefield environment to realize underwater vehicle route planning. This method not only considers obstacles, enemy detection probability, marine depth and route length, but also considers performance constraints such as the safe depth of the aircraft and the turning radius. This method can better realize path planning in a complex marine environment, but the planning speed is slow [29]. Zhu et al. introduced marine factors such as seabed topography, marine currents, marine climate, sea depth, and sea water transparency into the submarine's path planning, which improved the path's adaptability to the complex marine environment. However, this method does not consider physical marine phenomena, all of which still have limitations [2]. Liu et al. introduced the maneuvering rules of underwater submarine vehicles into path planning. This algorithm can find a suitable search range according to different state distributions, and has good adaptability. The results of simulation experiments prove that the algorithm can obtain a safe dynamic path with optimized performance. However, the algorithm still simplifies the complex marine environment [30]. Therefore, it can be seen that there is still a lot of work to be done in path planning in a complex marine environment. This article hopes to design a submarine underwater path planning algorithm that can adapt to the real marine environment by introducing physical marine phenomena into path planning.
This article mainly considers introducing the following marine phenomena into underwater path planning. 1) Marine acoustic environment. The acoustic environment of the marine mainly includes the acoustic velocity and the acoustic velocity gradient. The acoustic velocity affects the range of the anti-submarine sonar. Therefore, the lower the acoustic velocity, the closer the submarine can get to the target position, so as to achieve the combat purpose of concealed attack. The greater the absolute value of the acoustic velocity gradient, the more severe the refraction of acoustic rays, which is more conducive to the concealment of the submarine. Moreover, the acoustic velocity gradient can represent the size of the sonic cline to a certain extent. In order to avoid the detection of surface anti-submarine equipment, the submarine should navigate below the sonic cline as much as possible. 2) Pycnocline. Pycnocline is a physical phenomenon that occurs between two layers of fluids with different densities, and it has an important impact on the safe navigation of submarines. Pycnocline is divided into positive and negative. The positive pycnocline is also called ''liquid seafloor'', which can carry submarines to sail above it to achieve the purpose of concealed approach to the target and saving fuel. The negative pycnocline is called the ''submarine cliff''. Owing to its physical characteristics, the submarine will ''fall deep'' when encountering it and face danger. Therefore, the submarine should take different maneuvering routes when encountering different types of the pycnocline. 3) Density front. The density front mainly refers to the narrow water body in the horizontal direction between two bodies of water with different densities. When the submarine navigates through the front with a sudden change in density, the hull may appear large inclination and heave. At this time, if the maneuvering measures are not appropriate, it may dive to a large depth or suddenly emerge from the water, endangering the safety and concealment of the submarine [31]. 4) Mesoscale eddy. Mesoscale eddy is a complex and large-scale marine vortex phenomenon. Submarine encounters a strong mesoscale eddy and will get out of control. What's more, it will greatly weaken the performance of underwater communication and detection and pose a great threat to the safety of underwater targets such as submarines.
Based on the improved artificial potential field ant colony algorithm, this article quantifies the above-mentioned marine physical phenomena through the cost function, and introduces it into the path planning algorithm.
The rest of this article is organized as follows. In the second section, the environment model will be introduced. Section III will introduce the Improved Artificial Potential Field Ant Colony Optimization (APF-ACO) used in this article. In Section 4, a cost function that takes into account marine environmental impacts will be presented. In Section V, a series of experiments will be implemented to verify the reliability of the newly proposed algorithm and the performance of the marine-related cost function in a real marine data environment. The conclusion will be given in Section VI.

II. MODELS
The environmental model represents the real marine environmental information with an abstract physical model, and its representation method directly affects the efficiency of the path planning. This section mainly introduces the environment model used in this article

A. SPACE COORDINATE SYSTEM
This article uses the method of spatially equally dividing grids to divide the path planning environment. As shown in Figure 1, the minimum grid interval corresponds to the minimum longitude, latitude, and depth interval of the database. The plane with the depth value of 0 is regarded as the isosurface with the Z value of 0. The vertical downward along the depth is the positive direction of the Z axis. The direction in which the latitude increases along the origin O is the positive direction of the X axis. The increasing direction along the longitude is the positive direction of the Y axis. Therefore, the underwater planning space is constructed as a rectangular parallelepiped equally divided along the longitude, latitude, and depth. We stipulate that the space is equally divided into (l + 1) parallel planes perpendicular to the Y axis and the Z axis along the X axis. Then this article divides each plane into (m + 1) × (n + 1) grids along the Y-axis and Z-axis. In this way, the planning space is preliminarily divided into (l + 1) × (m + 1) × (n + 1) subspaces.
In the planning space, each point in the subspace can be spatially positioned with coordinates (x, y, z), as shown in Figure 2. In this way, the environmental modeling of the planned space is initially completed.

B. UNIT-BODY
Considering that most of the physical phenomena in the complex marine environment cannot be characterized by a simple data point, it is necessary to calculate the physical quantities within a certain range to get the true performance of marine parameters, such as mesoscale eddy. Therefore, on the basis of the above environment modeling, this article has proposed a data cube modeling method.
A data cube is composed of 8 adjacent data points, as shown in Figure 3, where adjacent data cubes share 4 data points on adjacent faces. In this way, a new subspace distribution form is formed, in which the characteristic value of each data cube is calculated by corresponding 8 data in it, such as average value, curl, gradient, etc. The article mainly constructs such a data cube as a subspace of the environment space based on the following reasons.
1) Compared with the scale of the submarine itself, the scale of the marine environment is much larger and the temporal and spatial changes are not drastic. Therefore, the method of reorganizing data points can reduce the amounts of useless calculations and improve efficiency.
2) The reconstruction of data can give each data cube more information, enabling it to realize the expression of marine physical phenomena. For the convenience of explanation, this data cube will be referred to as the unit − body in the following.

III. IMPROVED ARTIFICIAL POTENTIAL FIELD ANT COLONY ALGORITHM (APF-ACO)
The potential field ant colony algorithm proposed in this paper is mainly aimed at the specific application scenarios of submarine underwater navigation. Compared with ordinary underwater vehicles, submarine has more stringent requirements on path planning algorithms due to its special military purposes, among which convergence speed and global optimality are the two most important aspects.
Due to its ergodicity and positive feedback, the common ant colony algorithm has a slow convergence speed and is easy to fall into local optimum, which is not in line with the background of submarine combat applications. To solve this problem, this paper proposes an improved ant colony optimization (APF-ACO) that introduces artificial potential field factors. In the initial stage of the algorithm, the introduction of artificial potential field factors can increase the difference between the grids to be selected, thereby speeding up the speed of ants choosing the correct path and improving the convergence speed of the algorithm. In the middle and late stages of the algorithm, the artificial potential field factor can input the information of the target point into the operation to induce the ants to always choose the route closest to the target point and avoid the result falling into local optimum.
The specific explanation of APF-ACO is as follows.

A. IMPROVED APF ALGORITHM
The original artificial potential field function is prone to deadlock when there are a large number of obstacles around the target point, which makes the algorithm unable to effectively converge. In order to solve this problem, this article VOLUME 10, 2022 improves the repulsion function of the artificial potential field algorithm. The improved repulsion function is shown in (1).
where ξ = (X − X d ) n , representing the n-th power of the Euclidean distance between the current node and the target point. ρ 0 is the value of the influence range of the obstacle, if the distance between the carrier and the obstacle is greater than ρ 0 , the repulsion is 0. n is a positive adjustment constant. The introduction of ξ can reduce the effect of repulsion when the carrier is close to the target point, so that the carrier can reach the target point normally.
For U rep (X ) to find the negative gradient, the repulsion calculation formula of the improved algorithm can be obtained as: where F rep1 (X ) and F rep2 (X ) respectively describe the repulsive force of obstacles pointing to the robot and the attraction of robots pointing to the target point. In this way, the repulsion function expression of the improved algorithm is constructed.
The gravitational expression of the artificial potential field is where k p is the target gain coefficient, and the attractive force of the gravitational potential field to the carrier is the direction of the negative gradient of the gravitational potential energy: where n 1 is the unit vector with direction between the carrier and the target point, ε 1 is the Euclidean distance between the carrier and the target point, ε 1 = X − X d . In this way, the resultant force of the artificial potential field can be obtained:

B. IMPROVED APF HEURISTIC FUNCTION
The heuristic function of the traditional ant colony algorithm only considers the distance between the current position and the next node. When there are a lot of obstacles near the endpoint, the algorithm is easy to fall into the local optimum.
Due to the positive feedback of the ant colony algorithm, the final planned path may not be the globally optimal path.
To solve the problems of the traditional ant colony algorithm, firstly, the heuristic function is improved by the distance and position information of the node and the endpoint searched by the artificial potential field method. Later, in the later stage of the path search, to avoid the potential field ant colony algorithm from falling into the local optimum, it is necessary to reduce the influence of the heuristic function on the path search. Therefore, the decreasing parameter of the heuristic function is introduced.
The improved heuristic function is shown in (7).
where η d is distance heuristic function, η F is potential field heuristic function.
where d Sj is the distance between the starting point S and the next node j, and d jG is the distance between the next node j and the target point G. According to this formula, the farther from the starting point S and the closer to the target point G, the larger d Sj and the smaller d jG , the greater the possibility of selecting node j. The improved distance heuristic function moves the ant away from the initial recognition point and approaches the target when selecting the next node to prevent the algorithm from falling into a deadlock.
η F = a F t ·ζ ·cosθ (9) where F t is the resultant force of the potential field. θ is the angle between the connection direction of the current node and the next node and the direction of the potential field force. The use of the base a coefficient expression is mainly to consider the unity of the order of magnitude of the two multiplication expressions. The existence of artificial potential field makes the ant colony algorithm converge at a faster speed, but at the same time it is more likely to fall into a local optimum in the latter part of the iteration. Therefore, the attenuation coefficient ζ is introduced. With the continuous iteration of the algorithm, the force of the situation gradually decreases. (10) where N k is the current iteration number, N max is the maximum iteration number. The heuristic function of this paper is obtained by multiplying the above formulas. The heuristic function can introduce the influence of distance and potential force into the path planning operation to improve the convergence speed of the algorithm and obtain a better path.

C. IMPROVED PHEROMONE DIFFUSION MODEL
The basic ant colony system partial pheromone update rule is The information speed update rule is uniform and does not consider the influence of pheromone diffusion on the distribution of pheromone, which deviates from the real ant colony system. To more truly restore the high-efficiency characteristics of the ant colony system and improve the convergence speed of the existing ant colony system algorithm, this article proposes a three-dimensional spherical pheromone diffusion model to improve the efficiency of the ant colony's use of pheromone.
We assume that the pheromone concentration obeys the Gaussian distribution, and the pheromone at the current point only diffuses to the adjacent advancing direction grid. The simplified pheromone diffusion model is a circumscribed sphere, as shown in Figure 4. where l ob represents the diffusion radius of the pheromone. We assume that the step length of each movement of the ant is 1 in the y and z directions, the l ob is √ 3 . In (7), calculate the potential field included angle θ for the adjacent grid of the current grid, and the adjacent grid with the smallest θ value is the pheromone diffusion direction. From this, determine the direction of the next grid j, and calculate the pheromone concentration diffused to point j as shown in (12).
where δ is the diffusion coefficient of the pheromone, and d is the Euclidean distance between the current node i and the next candidate node j. Based on the above improvements, calculate the state transition rule of the ant from point i to the next grid point: where α and β are the two adjustable parameters which can affect the performance of the algorithm. τ ij = τ ij + τ ij , τ ij is the global update rule of pheromone, and its formula is as follows. where The flowchart of the algorithm is shown in Figure 5.

IV. MARINE-RELATED COST FUNCTION
In order to introduce the physical phenomena in the marine into the route planning of the submarine, this article extracts the expression of its physical characteristics and uses the cost function to quantify it. By introducing it into the heuristic function of the algorithm, marine physical phenomena can affect the result of path planning.

A. CONCEALMENT COST FUNCTION
The concealment cost function F C proposed in this article for submarine underwater path planning mainly considers two parameters: acoustic velocity and acoustic velocity gradient. In addition, to eliminate the influence of the dimension on the cost function, the cost function is processed with the idea of standardization.

1) ACOUSTIC VELOCITY C
The concealment cost function of submarine path planning is proportional to the acoustic velocity. The formula is where C refers to the average value of the acoustic velocity of 8 data points in the unit-body, C max refers to the largest VOLUME 10, 2022 value of all the acoustic velocity data points obtained by this method, and C min is the smallest value.

2) ACOUSTIC VELOCITY GRADIENT C
The concealment cost function is inversely proportional to the absolute value of the acoustic velocity gradient. Moreover, the acoustic velocity gradient can represent the size of the sonic cline to a certain extent. In order to avoid the detection of surface anti-submarine equipment, the submarine should navigate below the sonic cline as much as possible. The expression is as follows.
where C(i, j, k −1) represents the acoustic velocity gradient value at the point directly above T i , expressed as the difference between the average acoustic velocity of the upper and lower surfaces in the vertical direction in the unit body, the positive direction is downward. The expression is as follows.
In summary, given the weight of each parameter, the concealment cost function can be formed: (19) where α c + β c = 1.

B. SECURITY COST FUNCTION
The safety cost function F s of path planning mainly considers three parameters: pycnocline, density front, and mesoscale eddy.

1) PYCNOCLINE
Submarines should adopt different route options when facing different types of pycnocline. This article considers this point and proposes the cost function shown as follows.
where ρ(T i ) is the density gradient of point T i , and its expression is

2) DENSITY FRONT
The density front poses a threat to the safe navigation of submarines. Considering its physical characteristics, this article uses the horizontal gradient of density to represent it. In addition, since the plane has two directions of latitude and longitude, there is a difference in horizontal density in the two directions. In order to reflect its characteristics more accurately, this article decides to use the larger of the two-direction density difference to characterize the horizontal density difference of the unit-body. The expression is as follows.
where x ρ(T i ) and y ρ(T i ) respectively represent the horizontal density gradients in the two directions at point T i . The expression is as follows.

3) MESOSCALE EDDY
In order to represent mesoscale eddy in a complex marine environment, this article decides to use the curl of the vector field in the unit-body to characterize it. Considering that the vorticities of the mesoscale eddy are all vertically upward, this article uses the calculation of the curl of the vector field composed of four points on the same plane of the unit-body to express the cost function. The expression is as follows. (24) where u is the velocity vector field composed of four points on the top plane of the unit-body. Since the mesoscale eddy is only on the horizontal plane, its curl only contains the quantity in the k direction. The value of its modulus can represent the degree of vortex in the flow field. In summary, the safety cost function formula of underwater path planning is as follows.
where α s , β s and γ s are weighting factors.

C. ECONOMIC COST FUNCTION
The economic cost function mainly considers the impact of marine currents. Let θ be the angle between the marine current and the horizontal, ϕ is the angle between the submarine's heading and the horizontal when passing through the unit body, V is the marine current velocity, and V is the submarine's velocity. Then when there is a marine current area in the unit-body, and the direction of the marine current is opposite to the submarine's heading, the function is as 37022 VOLUME 10, 2022 follows.
When the direction of the marine current is the same as the heading of the submarine: In this way, the cost function expression of submarine underwater path planning can be obtained as: In this way, we can get the final expression of the algorithm heuristic function as shown below.

V. SIMULATION EXPERIMENT AND ANALYSIS
To

A. ALGORITHM COMPARISON ANALYSIS
We compare the path planning results of the APF-ACO and optimized ACO in a variety of obstacle environments. The parameter settings of the algorithm are shown in Table 1. We design experiments to compare the path planning results of the four different implementations described above. Each test was independently executed ten times for each implementation of each test environment. Table 3 shows the arithmetic mean, worst, best and standard deviation of 10 independent tests on each environment (Env1, Env2, . . . , Env5) to account for the randomness of the implementations.
As can be seen from Table 2 and Figure 6, the Optimized ACO and APF-ACO can achieve better results than the Original ACS and the Original APF in each test environment, proving the superiority of the two improved algorithms. In addition to this, compared with the Optimized ACO, APF-ACO performs better on the best planning result and the average planning result. Except for the same optimal path planning results in the third test environment, the optimal path planning results of APF-ACO are better than those of the Optimized ACO in the other four test environments. In addition, the control of the worst experimental results and standard deviation of APF-ACO is also significantly better than that of the Optimized ACO, which makes the algorithm have better stability and meets the needs of submarine underwater path planning in practical situation In this article, five test environments with different obstacle distributions are constructed to study and analyze the performance of each algorithm in different experimental environments.
The establishment of different experimental environments mainly considers the distribution concentration of obstacles near the starting point and the end point and the overall distribution shape of obstacles. For example, the test environment 1 improves the distribution density of obstacles near the end point. test environment 5 improves the distribution density of obstacles near the starting point. The distribution of obstacles in test environment 2 is more concentrated, while the distribution of obstacles in test environments 3 and 4 is more discrete.
It was concluded in the previous part of this article that APF-ACO and Optimized ACO have comparative advantages, so the path planning results of these two algorithms are compared separately. The visualization results are shown in Figure 7 and Figure 8. It can be seen that both of the algorithms can obtain a feasible, stable and satisfactory path. However, it can also be seen from the figure that the optimized ACO algorithm presents more tortuous paths and more turning points when the obstacles are denser. In contrast, the APF-ACO algorithm benefits from the introduction of potential field force parameters, and the constructed path can avoid obstacles earlier to obtain a more stable path. Figure 9 shows the comparison of the convergence speed of the APF-ACO algorithm and the optimized ACO algorithm. As can be seen from the figure, the APF-ACO algorithm proposed in this article can converge faster, and the convergence speed is improved by 65.7% compared with the optimized ACO algorithm. Not only that, the APF-ACO algorithm benefits from the introduction of potential field force parameters, which can obtain better path selection in the early stage of iteration, and can obtain better convergence results after convergence.

B. ENVIRONMENT AND ALGORITHM PREPARATION
In order to verify the application of the algorithm in a real marine environment, this article constructs a virtual underwater environment for simulation experiments, as shown in Figure 10, where the ellipsoid represents obstacles, and the seabed terrain is simulated by mathematical functions. The expression is as follows.
where X , Y and Z represent the 3D space domain, the parameters x 0 and y 0 represent the peak position of the seabed, the parameters σ , a and b control the shape of the mountain peak in the seabed, and h represents the height of the mountain peak in the seabed. The simulation environment is composed of 40 × 20 × 20 unit-bodies, and the passability of each unit-body is calculated by using a grid, and the passable grid is calculated in advance through a cost function to calculate its passing cost. In the path planning, in order to reflect the reliability of the cost function to the greatest extent, the step length of the ant is always set to 1unit of length and the side length of the steering window is set to 1.
In order to confirm the influence of every marine phenomenon on path planning, only one marine-related cost function is considered in each experiment, and the weights of cost functions related to other marine elements are reset to 0.
In addition, the marine environment model involved in this article is set as follows.

1) ACOUSTIC ENVIRONMENT
This article uses real marine acoustic velocity data to construct marine acoustic velocity environment. The data used is the acoustic velocity analysis and forecast data provided by the National Marine Data Center (NMDC) of China, and the time is November 11, 2021. In this experiment, the α c and β c are set to 1, and other weight coefficients are set to 0. The cross-section shows the vertical gradient of acoustic velocity, red represents the positive acoustic velocity gradient, and blue represents the negative acoustic velocity gradient. The vertical sections show the magnitude of the acoustic velocity, red represents the higher acoustic velocity, and blue represents the lower acoustic velocity.

2) PYCNOCLINE
In order to verify the influence of the density cline on the path planning algorithm, this article uses the density field analysis and forecast data provided by NMDC on November 8, 2021, to construct a simulation environment. The blue area in the cross-section represents the negative pycnocline, and the red area represents the positive pycnocline. In this experiment, the α s is set to 1, and other weight coefficients are set to 0.

3) DENSITY FRONT
This article uses the density field analysis and forecast data provided by NMDC on November 8, 2021, to construct a simulation environment containing density fronts. The red area in the two vertical sections represents the density front that is positive in the X direction, and the blue area represents the density front that is negative in the X direction. In this experiment, the β s is set to 1, and other weight coefficients are set to 0.

4) MESOSCALE EDDY
This article uses the National Oceanic and Atmospheric Administration (NOAA) HYCOM (Hybrid Coordinate Ocean Model) analysis ocean current data to build a simulation environment containing mesoscale eddy. The red arrow is the flow vector, and it can be seen that there is an obvious vortex. In this experiment, the γ s is set to 1, and other weight coefficients are set to 0.

C. RESULTS AND ANALYSIS OF MARINE-RELATED COST FUNCTION EXPERIMENT
The visualization results of the marine environmental experiments are shown in Figure 11. It can be seen from figure VOLUME 10, 2022 FIGURE 11. Visualization results of the effects of acoustic (a, b, c), pycnocline (d, e, f), density front (g, h, i), and mesoscale eddy (j, k, l) on path planning in three marine environments.  10 that the marine environment cost functions introduced in this article have a significant impact on the results of path planning.
Compared with the control test group that does not consider the influence of the marine environment, the path obtained by the algorithm proposed in this article can maintain the navigation under the position with a large acoustic velocity gradient in the acoustic environment and choose the position with the lower acoustic velocity to pass as far as possible. In the density field environment, the path can avoid the negative pycnocline and density front location seas and navigate above the positive pycnocline location. In the flow field environment, compared with the control group that only considers the flow field direction and size, the path obtained by the algorithm proposed in this article can avoid areas with large vorticity values, further ensuring the navigation safety of the carrier. Table 3 and Figure 12 shows the comparison of the navigation cost of the two paths in all experiments. It can be seen from the results that the algorithm proposed in this article can effectively reduce the submarine navigation cost underwater. In the acoustic environment, the pycnocline environment, the density front environment, and the mesoscale eddy environment, the average navigation cost is reduced by 43.5%, 74.41%, 30.11%, and 25.95%, respectively, which has a significant effect.
In short, the algorithm proposed in this article is sensitive to the parameters of the marine environment, which can truly reflect the influence of the marine environment on the path planning results, and improve the survivability of the submarine underwater.

VI. CONCLUSION
In this article, we propose an improved artificial potential field ant colony Optimization (APF-ACO) considering the elements of the marine environment, which is used to deal with the route planning problem of submarines in the complex marine environment. We design an algorithm comparison experiment between APF-ACO and Optimized ACO, and test the statistical significance of the experimental results. The experimental results show that compared with Optimized ACO, APF-ACO improves the convergence speed by about 65.7%, and has more stable and excellent path planning results.
In addition, this article puts forward the concept of unitbody to express the characteristics of local physical field, and use cost function to quantify marine physical phenomena, which successfully introduce the influence of the marine environment into the calculation of path planning. We used real marine data to construct experimental models, and the performance of the algorithm in the acoustic environment, the pycnocline environment, the density front environment and the mesoscale eddy environment is verified respectively. The experimental results show that compared with the control group, APF-ACO can make the path sensitive to the marine environment, reducing the path cost by 43.5%, 74.41%, 30.11% and 25.95% in the above four marine environments, respectively.
The submarine's underwater path planning also needs to consider its own special mission background and complex dynamic characteristics. Under different mission backgrounds, the concealment and economic requirements of submarines are different, and the indicators of the cost function should be adjusted flexibly; the complex dynamic characteristics of submarines will affect factors such as the range of pitch angle variation and the minimum turning radius. The above contents have a non-negligible impact on the underwater path planning of submarines, and these are also the focus of our follow-up research. VOLUME 10, 2022 BAO LI received the B.Sc. degree in measurement & control technology instrument from the School of Automation, Northwestern Polytechnical University, Xi'an, China, in 2008, and the Ph.D. degree in navigation from the Department of Navigation Engineering, Naval University of Engineering, Wuhan, China, in 2013. He is currently a Lecturer at the Naval University of Engineering. His work focuses on GNSS signal processing, navigation countermeasure, and GNSS precision positioning technology.
ZHIWEN NING was born in Hunan, China, in 1997. He received the bachelor's degree in aircraft manufacturing engineering from the Civil Aviation Flight University of China, in July 2019. He is currently pursuing the master's degree in navigation, guidance, and control with the Naval University of Engineering.
Since 2019, he has published articles at the 2020 International Guidance, Navigation and Control Academic Conference. His research interests include integrated navigation and the application of carrier interference magnetic field compensation.
Mr. Ning was awarded in the China National Scholarship in 2017 and the title of Sichuan Excellent University Graduate in 2019.
YANG CHANG was born in Qinhuangdao, Hebei, China, in 1993. He received the bachelor's degree in electrical engineering and automation from the Ocean University of China, in July 2016. He is currently pursuing the master's degree in navigation, guidance, and control with the Naval University of Engineering.
Since 2019, he has published articles in academic journals, such as the Journal of Sensor Technology, and some international conferences. His research interests include the application of integrated navigation technology and fiber optic gyroscope error compensation.
Mr. Chang was awarded the title of Outstanding Student while studying at the Ocean University on China. VOLUME 10, 2022