An Enhanced Virtual Force Algorithm for Diverse k-Coverage Deployment of 3D Underwater Wireless Sensor Networks

The combination of Wireless Sensor Networks (WSNs) and edge computing not only enhances their capabilities, but also motivates a series of new applications. As a typical application, 3D Underwater Wireless Sensor Networks (UWSNs) have become a hot research issue. However, the coverage of underwater sensor networks problem must be solved, for it has a great significance for the network’s capacity for information acquisition and environment perception, as well as its survivability. In this paper, we firstly study the minimal number of sensor nodes needed to build a diverse k-coverage sensor network. We then propose a k-Equivalent Radius enhanced Virtual Force Algorithm (called k-ERVFA) to achieve an uneven regional coverage optimization for different k-coverage requirements. Theoretical analysis and simulation experiments are carried out to demonstrate the effectiveness of our proposed algorithm. The detailed performance comparisons show that k-ERVFA acquires a better coverage rate in high k-coverage sub-regions, thus achieving a desirable diverse k-coverage deployment. Finally, we perform sensitivity analysis of the simulation parameters and extend k-ERVFA to special cases such as sensor-sparse regions and time-variant situations.


Introduction
The research on Underwater Wireless Sensor Networks (UWSNs) has a wide application prospect in marine hydrological data collection, marine pollution detection, water quality monitoring [1][2][3], etc. The coverage optimization of UWSNs has a great significance for the network's capacity for information acquisition and environment perception, as well as its survivability/lifetime [4]. In addition, the coverage effectiveness of UWSNs has a direct bearing on the properties of the network, such as its communication bandwidth and computation capacity, thus determining the quality of service to a certain extent [5].
Existing research on coverage optimization mainly aims at ground sensor networks and can be classified into three categories: the target coverage [6], the regional coverage [7], and the barrier coverage [8][9][10]. The target coverage requires that sensor networks can monitor target nodes. In regional coverage, any point inside the region must be covered by at least one active node. The barrier coverage studies the probability of detecting moving objects when they pass through the monitoring region. In [6,7,[10][11][12], a variety of coverage control methods were proposed, among which the Virtual-Force Algorithm (VFA) by Zou et al. [7] has attracted considerate attention and is regarded • We analyze and derive the minimal number of sensor nodes needed to build a specific diverse k-coverage UWSN.

•
We design an enhanced virtual force algorithm k-ERVFA to solve the non-uniform k-coverage optimization problem of sub-regions with different interest levels in UWSNs.

•
We extend our proposed algorithm to special cases such as sensor sparse regions and time-variant situations.
The rest of this paper is organized as follows: Section 2 summarizes related studies on sensor deployment methods of UWSNs and the k-coverage problem of WSN. In Section 3, we introduce the 3D-UWSNs model with different k-coverage requirements. We then derive the minimal number of nodes needed to meet the specific k-coverage requirement. In Section 4, we give a detailed description of the k-ERVFA algorithm. Simulation results are given in Section 5. In Section 6, we give discussions about simulation parameters and extend k-ERVFA to special cases. Section 7 concludes the whole paper.

Virtual Force Algorithm in Wireless Sensor Network
The VFA has attracted considerable attention and is considered an effective method to solve the coverage problems. For example, Zou et al. [7] achieved the coverage control in two-dimensional sensor networks based on VFA, which moves the node according to the calculated position in each step. In this algorithm, a node will be subject to three resultant forces (obstacles, neighbor nodes, and targets to be covered), including gravity and repulsion. Tao et al. [19] proposed PFCEA (Potential Field-based Coverage-Enhancing Algorithm) based on the rotatable directional sensing model. Furthermore, Huang et al. [20] put forward node re-deployment algorithm PRMCA (Probability model based Rotate or Move along a fixed direction Coverage-enhancing Algorithm) based on the probability model and the virtual force, where the sensor node rotates or moves along the fixed direction in a two-dimensional plane in order to achieve the coverage effect. Heo et al. [21] suggested a coverage control algorithm IDCA (Intelligent Deployment and Clustering Algorithm), and the virtual force between nodes is represented by the result of the expected deployment density divided by the current deployment density. Ma et al. [22] proposed a coverage-control VFA-ACE (Area Coverage Enhance) algorithm based on the virtual potential field, aiming at three-dimensional directional sensor networks. Tan et al. [23] proposed a three-dimensional space deployment algorithm applied to continuous target tracking. Virtual force in this algorithm is generated by the inter-node force, the obstacle repulsive force, the monitored-path attractive force, and the tracking target together.

Existing Studies on k-Coverage of Wireless Sensor Networks
The k-coverage problem of wireless sensor networks has been extensively studied and can be divided into two categories: the two-dimensional cases and the three-dimensional cases. Huang et al. [24] studied two-dimensional k-coverage with geometric analysis and achieved improved coverage. Wan et al. [25] derived the upper and lower bound of the minimal density of nodes in the two-dimensional k-coverage problem using the methods based on the probability model. The studies mentioned above focused on the thick deployment regions, such as a square or a circular area. With finite thin strip regions, Balister et al. [26] calculated the deployment density of sensor nodes needed to maintain both network connectivity and the k-coverage property. Qiu et al. [27] proposed a distributed cooperation scheme based on a local k-order Voronoi diagram in which nodes cooperate in hole detection and recovery. Esnaashari et al. [28] proposed CLA-EDS (a Cellular Learning Automata-based Enhanced Deployment Strategy), a modification of the CLA-DS (a Cellular Learning Automata-based Deployment Strategy) [29], in which a learning algorithm is applied to meet the time-variant and diverse k-coverage requirements in the two-dimensional cases. Yu et al. [30] studied the k-coverage problem in WSNs and proposed protocols based on the coverage contribution area. The proposed protocols achieved low sensor spatial density and prolonged the network lifetime.
For three-dimensional k-coverage problems, Alam and Haas [31] proposed a deployment strategy based on the Voronoi diagram on the tessellation of the truncated octahedral in 3D space. However, k-coverage was not studied in [31], and only one-coverage was guaranteed. Huang et al. [32] developed a polynomial time algorithm to solve the α-coverage problem in three-dimensional sensor networks. Ammari and Das [33,34] designed a distributed hybrid forwarding protocol, which guaranteed k-coverage, and discussed the minimal deployment density of sensor nodes required in the three-dimensional k-coverage problem. In [35], a localized, pseudo-distributed scheme was further proposed to achieve k-coverage in 3D duty-cycled WSNs, but the diverse k-coverage requirements in different sub-regions were not mentioned.

The Existing Deployment Methods of UWSNs
The existing deployment methods of three-dimensional underwater sensor networks can be classified into two categories: the sea-bottom deployment methods and the sea-column deployment methods. Table 1 clearly displays the deployment classification of UWSN.
In the sea-bottom deployment, the seabed plane under the monitoring area is divided into many triangular grids, with nodes deployed in the intersections of the grids, so as to achieve full coverage of the monitoring area with a minimum number of nodes. Due to the 3D feature of the underwater environment, most applications require the network to collect subaqueous or aquatic data. Therefore, it is hard for the sea-bottom deployment methods to meet such demands.
In the sea-column deployment, nodes are deployed into 3D underwater space in order to achieve 3D coverage of the monitoring area. The existing sea-column deployment strategies fall into two categories: the uniform deployment and the non-uniform deployment. The uniform deployment requires the sensor nodes to be distributed uniformly in the monitoring area. Many uniform coverage optimization methods [36][37][38] have been proposed by continuously adjusting the diving depth of nodes and reducing sensing overlap areas of adjacent nodes. Alam et al. [31] proved the optimality of the truncated octahedron model in the problem of the three-dimensional coverage of sensor networks and provided various deployment strategies for practical applications. However, the coverage requirements of many practical applications are non-uniform. Some sub-regions of the whole monitoring area need to be covered by more than one sensor node; therefore, the uniform deployment-oriented optimization algorithms cannot achieve a satisfactory node deployment. The non-uniform deployment strategy, on the other hand, allows sensor nodes to be deployed non-uniformly according to the distribution states of underwater targets. Many research efforts working on the non-uniform deployment strategy are based on the concept of the "target interest event", which means the purpose of node deployment is to cover the targets and interesting events instead of the whole monitoring area [39] and to make the density distribution of sensor nodes be consistent with that of the target events. The non-uniform deployment strategy is more effective in practical applications and accords with the sparsity characteristic of sensor nodes deployed in UWSNs [6,21]. Aitsaadi et al. [40] proposed DDA (Differentiated Deployment Algorithm) for water quality monitoring in a closed lake. The algorithm took the distribution characteristics of pollutants into consideration and deployed nodes non-uniformly using mesh grid representation, thus achieving diverse coverage of the monitoring area. Golen et al. [41] divided the monitoring area into sub-regions based on environmental factors such as acoustic characteristics and then optimized node deployment using the Game Theory Field Design (GTFD) model. However, this method is difficult to implement due to its high complexity. Zakia et al. [42] proposed a heuristic deployment strategy based on sub-cube tessellation and mixed integer linear program optimization. Liu et al. [43] studied the topology control of UWSNs with diverse coverage requirements and proposed two algorithms in which the sensing radius of a sensor node is adjusted to achieve the diverse coverage requirements. Wang et al. [44] proposed a self-deployment algorithm for maintaining the maximum coverage and connectivity in underwater acoustic sensor networks based on an ant colony optimization in 2019. They carried out the greedy strategy and improved the path selection probability and pheromone update system. However, current non-uniform deployment strategies based on event-driven cannot fully meet the demands of underwater sensor network applications. The main reasons are listed as follows: • Most of the event-driven coverage algorithms obtained only one-coverage in the target event area, which might not meet the demands of practical applications. For instance, target localization tasks require that any point within the region must be covered by at least three sensor nodes.

•
Many methods mainly considered determined target events while ignoring the uncertain events in the open underwater environment. Not much attention was paid to the difference of the levels of coverage required in different sub-regions. Some research works considered diverse k-coverage requirements in different sub-regions, but only in two-dimensional situations [28], which needed to be extended to three-dimensional cases.
Current research works covered some aspects of the k-coverage problems in three-dimensional space, but few methods have been proposed to solve the diverse k-coverage problem in UWSNs. Some solutions for UWSNs were only extended from the general 3D k-coverage cases. Therefore, we propose an enhanced virtual force algorithm k-ERVFA (k-Equivalent Radius enhanced Virtual Force Algorithm) to solve effectively the diverse k-coverage problem in 3D UWSNs. The algorithm is based on the concept of "k-virtual force", and it studies the minimum deployment density of sensors required for 3D UWSNs' k-coverage according to the practical application. The proposed algorithm can realize different k-coverage requirements in 3D UWSNs and can be extended to general situations, such as the cases of node sparsity and k-coverage time-varying requirements.

Preliminaries and Problem Statement
In this paper, the proposed network model consists of anchors, underwater sensor nodes, and buoy nodes, which are connected to the sensor nodes via a rigid cable. Our goal is to achieve a desirable diverse k-coverage rate with minimal adjustment cost. We made the following assumptions: 1. The 3D underwater space A is a cube with side length L 0 (expressed as A = L 0 3 ), and its bottom surface is the plane of the seafloor with the X-Y axis, while the Z-axis extends forward to the sea surface. The i th sub-region with the k-coverage requirement is denoted by A k,i . The volume of A k,i is denoted as V A k,i , and A k,i should be k-covered. Both V A k,i and k are known in advance. 2. At the initial phase, a fixed number of sensor nodes and buoy nodes are sprinkled randomly by an airplane. The coordinates of sensor nodes are randomly distributed. The buoy nodes are fixed through anchors on the seafloor so that the (x,y)coordinates of the sensor nodes are also fixed. The diving depth of the sensor nodes can be adjusted by stretching the rigid cable between the sensor nodes and buoy nodes. Each underwater sensor node S i communicates with its buoy node B i via a wired cable and reports its current depth, and the perception data are collected. Each buoy is equipped with a GPS module in order to obtain its coordinate. In addition, the rigid wired cable guarantees that the buoy will not drift far. 3. The sink nodes collect the coordinates of buoy nodes and the corresponding underwater sensor nodes, then run the re-deployment algorithm, which calculates the re-deployment positions (only in the vertical direction) of sensor nodes to meet the k-coverage requirements of the sub-regions. When the calculation is complete, the sink nodes send a re-deployment message to the buoy nodes, and the buoy nodes adjust the length of the underwater wired cable; thus, the diving depths of sensor nodes are adjusted, and re-deployment is implemented. 4. Assume that the satisfactory k-coverage rate is denoted by η, whose value needs to be determined by practical applications, while a fixed value of the k-coverage rate is necessary for the more concise theoretical derivation and experimental simulations in this paper. In order to make the analysis more intuitive, the value of η will be set to 89% as a demonstration in this paper. Note that the input parameter η is still optional and adjustable. The proposed model and algorithm of this paper can adapt to different η values to meet various requirements.
Under the above assumptions, the proposed algorithm will calculate the number of sensors needed for a given diverse k-coverage problem and adjust the coordinates of sensor nodes automatically to achieve the required coverage rate. The explanations of the main notations are presented in Table 2.

Notation Description
r The coverage radius of sensor nodes k The multiplicity of the coverage requirement r k The k-equivalent radius Table 2. Cont.

Notation Description
The k-coverage rate of A k,i S x The x th sensor node L 0 The side length of the underwater monitoring area η The satisfactory coverage rate θ The node redundancy coefficient ρ min The minimum deployment density of sensor nodes n min The minimal number of sensor nodes K con f , K attr,k , K Ob,k The virtual force coefficients ∆L max The maximum moving distance of sensor nodes in each iteration N max The maximum iteration times in each round of k-ERVFA

Related Definitions of the Network Model
Definition 1 (Sub-regions with different k-coverage requirements). The whole underwater monitoring area consists of q sub-regions of various levels of coverage requirements. Namely, for any sub-region A k,i (i = 1, 2, ..., q, k ∈ Z + ), any point inside A k,i should be at least k-covered.
Definition 2 (Perceptual model). The Boolean perceptual coverage model is adopted in this paper. Let r be the coverage radius of the sensor node and (x i , y i , z i ) be the coordinates of sensor node S i . Then, the region covered by S i is a three-dimensional sphere with a radius of r and the center at (x i , y i , z i ). The region is called the perceptual sphere region, denoted as ε(i). Namely, it is the set of all points p that conform to Equation (1), in which d(S i , p) denotes the Euclidean distance between S i and p.
Definition 3 (k-equivalent radius). Let r be the coverage radius of the sensor node. Its k-equivalent radius is then defined as: Definition 4 (k-conflict radius). In the classic virtual force algorithm, two sensor nodes with Euclidean distance less than 2r are considered as conflicting nodes. Similarly, we consider two nodes to be k-conflicting nodes when the Euclidean distance between them is less than two-times the k-equivalent radius. The k-conflict radius r con f ,k is then defined in Equation (3): Definition 5 (k-coverage rate). With sub-region A k,i , we define its k-coverage rate (denoted as C k,i ) as the volume of the k-covered region within A k,i divided by the total volume of A k,i : where ∆(p k,i ) is the volume of point p k,i , which reflects the granularity of division.
Definition 6 (Sphere with the k-equivalent radius). The three-dimensional sphere with a radius of r k (r k = r k , k = 1, 2, 3...) is defined as a sphere with the k-equivalent radius. Note that the actual physics coverage radius is still r.

Analysis of the Deployment Density Satisfying Different k-Coverage Requirements
In this section, we consider the deployment density needed to meet the specific k-coverage requirement in a sub-region of UWSNs. Notice that [33,34] studied the minimal deployment density to meet the k-coverage requirement in 3D underwater region. In this paper, we obtain an improved minimum volume density of sensor nodes compared to that in [33,34]. The result is as follows: Claim 1. The minimum volume density of sensor nodes needed to obtain a desirable k-coverage rate (89%) is denoted as ρ kmin , and the value of ρ kmin is 3 √ 3/8r 3 , when k= 1; 3/πr 3 , when k = 2; 9/2πr 3 , when k = 3... (for other k values, refer to Equation (13)).
To reach the above conclusion, firstly, we divide the 3D underwater space into cubic grids as shown in Figure 1a (four grids are illustrated here). Let L be the side length of the small cubic grid and r be the coverage radius of the sensor node. Each sensor node is placed in the center of the corresponding cubic. We can see from Figure 1a that if the eight vertices are not covered by the centered sensor, it will not be covered by sensors from other grids either. In this case, similar coverage breaches will occur in every cubic grid in the 3D space. Namely, when r is not large enough (compared to L), the coverage rate will not be satisfactory. However, when r increases to the case shown in Figure 1b, all eight vertices of the cubic grid will be covered. Thus, in order to ensure a 100% one-coverage rate, we have: For the convenience of the solving process, L is represented as L = 2r , where k is the coverage multiplicity required in this given sub-region and m is an unknown parameter. According to the discussion above, we have: Sensor nodes are deployed uniformly in the whole 3D underwater space, so the density of sensor nodes is: Consider some point p inside the given sub-region and a sphere centered at p with a radius of r; the number of sensor nodes inside the sphere is: Any sensor S inside the sphere covers point p, while for any sensor S outside this sphere, the Euclidean distance between S and p is greater than r, which means that p will not be covered by S . In order to meet the k-coverage requirement, there shall be at least k sensors inside this sphere; thus, we have: Theoretically, k-coverage requirement will be met when both Equations (6) and (9) are satisfied simultaneously. However, the analysis above mainly considers the average number of sensors. In practice, the 100% k-coverage rate requirement cannot always be met due to randomness in sensor deployment.  For example, as shown in Figure 2, which is the 2D projection of the 3D case mentioned above, blue pentagrams stand for sensor nodes. Sphere 1 and Sphere 2 have the same volume, but Sphere 1 contains two sensor nodes, while Sphere 2 contains only one sensor node. Namely, all particles inside Sphere 1 are two-coverage, while the particles inside Sphere 2 and outside Sphere 1 are one-coverage.
Therefore, an adjustable parameter θ is introduced to ensure higher k-coverage in the actual deployment, which makes our solution more suitable for random deployment of initialization. Therefore, the above condition m 3 πk 6 ≥ k should to be adjusted to: where θ is introduced to provide some redundancy of nodes so that a high k-coverage rate can be guaranteed. θ is related to both k and η, expressed as θ = θ(k, η). The value of θ should be carefully chosen, for a higher θ value will lead to a higher sensor deployment cost, while the satisfactory coverage rate will not be obtained if θ is too small. We propose the Optimal θ Searching Algorithm (OθSA) to help decide the appropriate value of θ. In this paper, the coverage radius of sensor node r = 10 m, and the side length of the 3D underwater region L 0 = 100 m. Sensor nodes were deployed as shown in Figure 1, with a deployment interval of L = 2r/(m · 3 √ k). In order to make the analysis more intuitive, the value of η will be set to 89% as a demonstration in this paper.
In order to meet different requirements of the k-coverage rate in practical applications, we can obtain the optimal θ value through the optimal θ searching algorithm (Algorithm 1). When η is varied, the new θ value can be calculated by this algorithm to match it. In addition, we have displayed different values regarding θ when η = 88% and 90%, respectively, as shown in Figure 3 and Table 3 below.      Input: k, r, L and ∆x. Output: The coverage rate.
From Figure 3, we find that a higher value of θ (greater than 1.8) makes little contributions to the improvement of the coverage rate, while causing an unnecessary deployment cost. Thus, it is not recommended to adopt a higher value of θ. In addition, for higher values of k(k ≥ 5), the optimal value of θ also should be determined by OθSA.
It can be seen in Figure 3 and Table 3 that the simulation results are consistent with the theoretical derivation. The basic idea of the algorithm can be described as follows: for a given k, through traversing the value of θ(k, η) from 1.0 with an increment of 0.1, the value of m can be calculated according to Equation (10); thus, the quantitative relation between r and L is known, and the k-coverage rate in the sub-region can be obtained. Whenever the calculated coverage rate reaches η, the algorithm will be terminated at this moment, and the value of θ(k, η) is optimal. This algorithm belongs to the Monte Carlo method (statistical simulation algorithm). The shorter the step size of space particle coordinates, the higher the precision and the higher the algorithm complexity. The time complexity of OθSA is where ∆x is the step size of spatial variation and MAXN is ceil(L 0 /L). In summary, in order to achieve a satisfactory k-coverage rate of one given sub-region, both Equations (6) and (10)  and a greater volume density of sensor nodes, which means the deployment cost is more expensive. We should adopt the minimum m that satisfies both Equations (6) and (10). That is: The minimum m and corresponding volume density of nodes to ensure an 89% k-coverage rate can then be calculated and is shown in Table 4. Table 4. The minimum m and corresponding volume density to ensure an 89% k-coverage.
The minimal number of nodes in any given k-coverage sub-region can then be calculated as: where V k is the volume of the k-coverage sub-region.

The Enhanced Virtual-Force Algorithm k-ERVFA
In this section, we propose the enhanced virtual force algorithm, i.e., k-ERVFA to solve the diverse k-coverage problem of 3D UWSNs. The model of the sensor sphere with k-equivalent radius (seen in Section 3.2) is used to cover the entire 3D underwater area.

Repulsion between k-Conflicting Nodes Based on Coulomb's Law
According to Coulomb's Law, we define the k-repulsion between two sensor nodes i and j as Equation (15): where n is the total number of sensor nodes in the sensor network, d ij is the Euclidean distance between sensor node i and sensor node j, r con f ,k is the k-conflict radius (seen Definition 4 in Section 3.2), − → e ji is the unit vector in the direction from S j to S i , and K con f is the coefficient of the Coulomb repulsion between nodes, which is a constant independent of the value of k.
Note that for sub-regions with different k-values, the only difference in Equation (15) is r con f ,k , which is 2 · r/k (seen in Definitions 3 and 4 in Section 3.2). For greater k values, both r k and r con f ,k are smaller, i.e., the repulsion exists only when the Euclidean distance between neighboring nodes is rather small (compared to cases with small k values). A smaller distance between nodes means a higher volume density, which leads to a higher k-coverage rate. This explains why the k-equivalent radius was introduced and how it facilitates obtaining diverse k-coverage.
In this algorithm, the k-repulsion between nodes was adopted to minimize the overlapping coverage area of neighboring nodes.

Attraction from the k-Coverage Requirement Sub-Regions
If there were only repulsion between nodes, sensor nodes would just disperse as far as possible so there will be coverage breaches between sensors; thus, k-coverage will not be achieved. In an ideal deployment scheme, more sensor nodes should move towards sub-regions with k-coverage (k ≥ 2) requirements to guarantee better k-coverage. Therefore, we introduce the attraction from the k-coverage requirement sub-regions on nodes to "pull in" more sensor nodes. In order to achieve diverse k-coverage, the attraction coefficients of sub-regions with different k-coverage requirements should differ. The attraction from a sub-region on nodes is considered as the attraction from the region's centroid. For a given k (k ≥ 2), the k-attraction from the j th k-coverage requirement sub-region on sensor S i is: where K attr,k is the attraction coefficient of the k-coverage requirement sub-regions on sensor nodes and varies with different k values and d ikj is the equivalent distance between sensor S i and the k-coverage requirement region A k,j (the jth k-coverage requirement sub-region for a given k), which is considered as the Euclidean distance between S i and the centroid of A k,j : The direction of − → F attr,ikj is determined by the unit vector − → e ikj , which starts from S i and ends at the centroid of A k,j . Now, considering the value of K attr,k , naturally, for sub-regions with higher k-coverage requirements, the attraction should be greater so that more sensors can be pulled into these sub-regions; thus, we will set greater K attr,k for higher k values. The numerical relation of K attr,k for different k values is determined by (the ratio of the coverage requirement multiplicity in different sub-regions), where k 1 and k 2 represent the coverage requirement multiplicity in different sub-regions, respectively. The numerical relation between K attr,k and K con f will be stated later in the simulation section.

Obstacle Repulsion from the "Fixed" k-Coverage Requirement Sub-Regions
The k-ERVFA algorithm sorts all sub-regions in a descending order of k values, then considers the sensor deployment in each sub-region successively. Namely, the sub-regions with high k values are considered as "fixed" after their deployment rounds are over. The sensor nodes inside of the "fixed" sub-regions will not participate in the future deployment process (i.e., the deployment process in sub-regions with smaller k values). Besides, in the following deployment process, the aforementioned "fixed" sub-regions should be considered as obstacles to prevent the nodes outside the "fixed" sub-regions from re-entering. This kind of re-entering will result in extra deployment cost and can be avoided when the obstacle k-repulsion is adopted.
We set the obstacle region Ob k,i of "fixed" sub-region A k,i to be the combination of A k,i and a ring region formed by the spheres with a radius of r k and the centroid at the boundary of A k,i , as shown in Figure 4. Ob k,i can be defined as: In the current deployment round, if sensor node S j is outside the fixed region A k,i and inside the obstacle region Ob k,i , it will receive repulsion from Ob k,i : where − → F ob,ikj represents the repulsion of obstacle region Ob k,i on sensor S j with the direction from the centroid of region A k,i to S j , d ikj is the distance between Ob k,i and S j , and K ob,k is the repulsion coefficient of Ob k,i .

k-Resultant Force of the Sensor Node
Based on the virtual force model and the discussions above, the k-resultant force on node S i is the vector sum of the aforementioned k-repulsion and k-attraction, namely: In a single deployment round, the sink node receives the coordinate of S i (x i , y i , z i ) from the buoy, calculates the k-resultant force that acts on S i , and then decides the moving direction and distance of S i in this round. Since the underwater sensor nodes can only move vertically, we need to obtain the component force along the Z-axis by projecting − → F i onto the Z-axis: where − → e z is the unit vector in the positive direction of the Z-axis. The moving distance and direction of S i are then decided by the magnitude (|F iz |) and direction of − → F iz . Let the maximum moving distance (along the Z-axis) of all the sensor nodes in a single round be ∆L max (the suitable value of ∆L max is discussed in Section 6). After − → F iz of each S i is calculated in each round, the maximum of |F iz | can be obtained and be denoted by max{|F iz |} (i = 1, 2..., n). The moving distance of S i in this round is then: The normalized calculation in Equation (22) represents the correlation between the resultant force that acts on the sensor node and the sensor's moving distance in each round, namely greater force corresponds to longer moving distance under the premise that the maximal moving distance in a single round is ∆L max .

Motion Pattern of Boundary Nodes Based on the Ideal Elastic Collision
Any given closed region has its boundary; here, the word "region" means either one specific k-coverage requirement sub-region A k i ,j or the entire underwater monitoring area A. For one sub-region A k i ,j , when the "fixed and even" algorithm (seen in Sections 4.6 and 6.1) is performed, we need to make sure that all sensor nodes inside A k i ,j always stay in A k i ,j (that is why these nodes are considered to be "fixed"). For the entire monitoring area A, naturally, any sensor node should not move outside of the boundary of A. To prevent nodes from escaping a closed region, we propose the motion pattern of boundary nodes based on ideal elastic collision.
In this paper, sensor nodes can only move along the Z-axis, so we only consider the z-coordinate. If the calculated redeployment position z i + ∆z i is outside the boundary, then the actual redeployment position should be modified as follows according to ideal elastic collision: When When z i = z i + ∆z i < A: where A is the z-coordinate of the upper horizontal boundary of the 3D region A, A is that of the lower horizontal boundary, and z i is the modified redeployment position of node S i after collision.

The Fix and Even Redeployment Algorithm
In this Section, we only describe the procedure of the fix and even algorithm; a further discussion and explanation can be seen in Section 6.1.
1. The input parameter of this algorithm is the coverage multiplicity k i . Set the state of any sensor S x inside A k i ,j to be "fixed". Namely, for any sensor inside A k i ,j , its range of motion should be within A k i ,j . 2. Calculate the k-resultant force on S x and its moving vector ∆z x , based on the discussions in Sections 4.1-4.4. If the calculated redeployment position z x + ∆z x is outside A k i ,j , modify the actual redeployment position according to Equation (23) or (24).

Repeat
Step 2 until the iteration number reaches N max or the k-coverage rate of A k i ,j is met, that is C k i ,j ≥ η, where η is the satisfactory k-coverage rate for practical application. The fix and even algorithm is described in Algorithm 2.

Algorithm 2:
The "fix and even" algorithm.
Input: k i , the specific multiplicity of the coverage requirement. Output: The coverage rate inside A k i ,j is improved. 1 N already ← 0 repeat /* N already is used to count iteration times */ 2 for all S x in A k i ,j do 3 D x ← " f ixed /* Set the state of any sensor S x inside A k i ,j to be "fixed" */ 4 calculate the k-coverage parameters: − → F xz , ∆z x , z x = z x + ∆z x /* according to Equations (15) to (22) */

Steps of the k-ERVFA Algorithm
The procedure of k-ERVFA algorithm (described in Algorithm 3) is as follows: 1. Sort {k 1 , k 2 ..., k m } (k i is the coverage multiplicity in different k-coverage requirement sub-regions) in descending order, and get {k 1 , k 2 ..., k m }, where k 1 > k 2 > ... > k m . For each k i , calculate the corresponding k-equivalent radius r k i = r/k i . 2. Check the current iteration number N already ; if N already > N max , break out the loop in Step 2; otherwise, substitute the k-equivalent radius calculated in Step 1 into Equations (15)- (22). The sink nodes calculate the k-resultant force and the moving vector for all sensor nodes, except for those whose states are set as fixed in Step 3. The redeployment positions of sensors in this iteration round are then determined (for boundary nodes, the ideal elastic collision model will be adopted if necessary). Now, calculate the k-coverage rate C k i ,j of A k i ,j , (j = 1,2, ... , m k i , where m k i is the total number of sub-regions that need to be k i -covered). Break out of Step 2 if N already > N max or the k-coverage rate for all A k i ,j is met. Otherwise, repeat Step 2, and increase N already by one. 3. For all sub-regions that need to be k i -covered, the redeployment process is partially done. Now, either the k-coverage rate for all A k i ,j has been met or N already has reached N max . Then, call the fix and even algorithm for all A k i ,j . By doing this, all sensors within A k i ,j are set to be fixed and will no longer participate in the following redeployment rounds (smaller k i ). Besides, the sensors within A k i ,j are homogenized, and a better coverage rate is achieved.

Simulation Results and Performance Analysis
We performed simulations of the k-ERVFA algorithm using MATLAB R2015a. The entire three-dimensional underwater monitoring area was set to be a cube with side length of 100 m. The diverse k-coverage requirement was set as follows: • a cubic sub-region that is centered on (25 m, 25 m, 25 m) needs to be three-covered, and the size of the three-coverage requirement sub-region is 30 m * 30 m * 30 m; • a cubic sub-region that is centered on (70 m, 70 m, 70 m) needs to be two-covered, and the size of the two-coverage requirement sub-region is 40 m * 40 m * 40 m; • the rest of the monitoring area needs to be one-covered.
Sensor nodes were randomly distributed initially. We set the numerical relation of K con f : K attr,2 : K attr,3 : K ob,2 : K ob,3 to be 1:2:3:2:3. With a given maximum iteration number, this proportion helps to achieve a satisfactory k-coverage rate according to the simulation results. The detailed simulation parameters are shown in Table 5.
Here, r is set to be 10 m, and according to Equation (14), the minimal number of sensor nodes needed to meet the diverse k-coverage requirements (89% was considered as the satisfactory coverage rate) in different sub-regions is: 591 for the one-coverage sub-region; 62 for the two-coverage sub-region; 39 for the three-coverage sub-region; and 692 in total.  Max iteration times N max 100 K con f , K attr,k , K Ob,k K con f :K attr,2 :K attr,3 :K Ob,2 :K Ob,3 = 1:2:3:2:3 However, in the simulation experiments, we deployed no more than 650 sensor nodes into the monitoring area for the reasons below. Firstly, according to Equation (10) and Table 4, the calculated number of nodes will achieve a 100% one-coverage rate in the one-coverage requirement sub-regions, which is unnecessary. Secondly, θ was introduced, which leads to redundancy in node deployment. In practice, k-ERVFA achieved a satisfactory k-coverage rate with approximately 650 nodes, as shown in the simulation results.
In the simulation experiments, the number of nodes varied from 400 to 650 with an increment of 50 each time. We compared k-ERVFA with the following sensor deployment algorithms: RD (initial Random Deployment algorithm), VFA (classic Virtual Force Algorithm) [7], CLA-EDS (a Cellular Learning Automata-based Enhanced Deployment Strategy) [28], IDCA (Intelligent Deployment and Clustering Algorithm) [21]. Note that the IDCA algorithm does not take the k-coverage requirement into consideration, and it should be modified to be compared with other algorithms. We modified IDCA by replacing d avg (the expected distance between nodes, which was an input parameter of IDCA) with d avg /k (similar to the k-equivalent radius in k-ERVFA); thus, the modified IDCA can deal with diverse k-coverage problems. In addition, CLA-EDS and IDCA only study two-dimensional sensor networks, so we extended them to the 3D case while constraining sensor nodes to move along the Z-axis only. We compared five algorithms under the same conditions. Figure 5a-c show the k-coverage rate (k = 3, 2, 1) achieved by different algorithms with the number of sensor nodes varying from 400 to 650.
As shown in Figure 5a, with 450 sensor nodes, for the three-coverage requirement sub-region A 3,1 , the three-coverage rate was 23.43% by RD, 32.92% by VFA, 61.89% by CLA-EDS, 62.97% by IDCA, and 82.45% by k-ERVFA, which was the highest among the five algorithms. Similarly, in Figure 5b, with 450 sensor nodes, for the two-coverage requirement region A 2,1 , the two-coverage rate was 45.32% by RD, 56.65% by VFA, 83.17% by CLA-EDS, 84.92% by IDCA, and 86.44% by k-ERVFA, which was also the highest. Note that although 450 was far less than 692 calculated by Equation (14), to achieve over an 89% k-coverage rate, the two-coverage and three-coverage rate of k-ERVFA reached 82% and 86%, respectively. It is obvious that k-ERVFA improved the k-coverage rate significantly.
Similar results were obtained when n varied from 500 to 650. The two-coverage and three-coverage rates of k-ERVFA were higher compared to other algorithms. It can be found that there was a significant advantage in the three-coverage rate, while the advantage was rather slight for the two-coverage rate.  Figure 5c demonstrates the one-coverage rate achieved by the five algorithms. IDCA had the highest one-coverage rate, followed by the CLA-EDS algorithm, then the classic VFA and k-ERVFA algorithms. Note that the Y-axis started at 0.75, and the difference of the coverage rate between k-ERVFA and IDCA was from 3.33% to 9.11% as n varied from 400 to 650, which shows that k-ERVFA increased the two-and three-coverage rate significantly with a slight sacrifice in the one-coverage rate. In addition, the one-coverage rate of k-ERVFA reached 91.87% when n = 450, which was 12.64% higher than RD, showing that its one-coverage performance was acceptable.
The main reason was that the k-ERVFA algorithm sorted all diverse k-coverage requirement sub-regions in terms of k values and considered redeployment in high k-coverage sub-regions preferentially. Namely, the algorithm deployed nodes into high k-coverage requirement sub-regions first, followed by the redeployment in low k-coverage sub-regions. In particular, for some point in a sub-region that needs to be k-covered, if it can only be (k − 1) covered because of the lack of one sensor node, then the other existing (k − 1) sensor nodes are wasted for this point. For higher k values, this kind of squander of sensors was severe and should be eliminated with first priority, which explains why k-ERVFA set priorities according to the k values.
In this respect, k-ERVFA was superior to CLA-EDS, which did not take priority coverage for the higher k-coverage requirement into consideration (according to [28], the CLA-EDS algorithm showed an obvious decrease in the high k-coverage rate). IDCA introduced the k-expected distance d avg /k, which was similar to the k-equivalent radius in k-ERVFA, and both algorithms were based on VFA, while the k-ERVFA algorithm achieved a better two-and three-coverage rate by considering the k-resultant force comprehensively and applying the fix and even algorithm to redeploy nodes within sub-regions. Compared to the classic VFA algorithm in which nodes are deployed uniformly, k-ERVFA pulled more nodes into the two-and three-coverage requirement sub-regions to achieve more desirable two-and three-coverage rates with a slight decrease (<3%) in the one-coverage rate.
In addition, in order to prove that it is reasonable to deploy 650 nodes in practical applications, we added a set of experiments with 700 nodes, and the experimental results are shown in Figure 6. Deploying 700 sensor nodes could achieve a 100% one-coverage rate and at least an 89% two-and three-coverage rate. However, compared with the result of 650 nodes, the effect of the improvement was not obvious. It also caused the redundancy of nodes and the increase of the deployment cost. In fact, with dynamic adjustment of node position based on the k-resultant force, k-ERVFA achieved a satisfactory two-and three-coverage rate (though the one-coverage rate decreased slightly) using only 650 nodes. Therefore, it is reasonable and optimal to deploy 650 nodes in the real scenario.   Table 5. Figure 7a shows the initial Random Deployment (RD) where the green dots are the sensor nodes located in three-coverage requirement sub-region A 3,1 , the red triangles are the nodes in two-coverage requirement sub-region A 2,1 , and all other blue asterisks are sensors in one-coverage requirement sub-regions. It can be seen that sensor nodes were non-uniformly distributed with RD, and the number of nodes in A 3,1 and A 2,1 were few. Figure 7b shows the distribution of nodes after the first round of k-ERVFA, i.e., the round in which spheres with three-equivalent radius (r 3 = r/3) were applied to achieve high a three-coverage rate. Due to the k-attraction of sub-region A 2,1 , nodes would aggregate towards A 2,1 , and the two-coverage rate increased. However, the node distribution in A 1,1 was still non-uniform, which lead to many coverage breaches. Figure 7c is the ultimate distribution after k-ERVFA has been performed. Compared with Figure 7b, nodes were more uniformly distributed in A 3,1 and A 2,1 due to the effect of the fix and even algorithm. There were fewer coverage breaches in A 1,1 so the one-coverage rate increased accordingly. Eventually, the k-coverage rates achieved by k-ERVFA with 600 nodes were: 92.67% for one-coverage, 97.54% for two-coverage, and 95.22% for three-coverage, which was consistent with the theoretical analysis.

The Principle of the Fix and Even Algorithm
After completing Step 1 of the k-ERVFA algorithm, the redeployment process in the k i -coverage requirement sub-regions was only partially done. For now, we know that either the coverage rate in all k i -coverage requirement sub-regions reached η or the iteration number reached N max . The latter is the more common case in practice. It is worth emphasizing that the value of N max should not be too large due to restrictions of energy and running time. Namely, in most cases, the iteration number reached N max , but the k i -coverage rate inside A k i ,j was far from satisfaction.
However, if we continue the node redeployment process in Step 2 of k-ERVFA, some nodes already inside A k i ,j may move away from A k i ,j due to the attraction from other sub-regions, causing unnecessary moving and computing costs. In order to avoid this situation, here we consider site preservation and bound these nodes inside A k i ,j , i.e., these nodes cannot move out of A k i ,j and will not participate in future redeployment rounds. This is the operation of "fix". Now that the nodes are fixed inside A k i ,j , their distribution is still uneven, and the k i -coverage rate is not satisfactory. Naturally, we want to homogenize the node distribution inside A k i ,j . This is the operation of "even".
With the fix and even algorithm, we fixed the sensor nodes inside the sub-region and then homogenized the node distribution, for the purpose of achieving a higher k i -coverage rate inside A k i ,j .

Discussion of the Value of θ(k, η)
In this section, we discuss whether there is a rough acceptable range for the value of θ(k, η) to guide node deployment in practice.
As shown in Figure 3, when θ < 1.4, some individual points in the curve of the f-coverage rate were above the corresponding points in the curve of the three-coverage rate, and these two curves tended to merge when θ > 2.2. For a given θ, θ · k increased with k, which provided more extra nodes. Namely, if we used the same θ(k, η) for different k values, then the redundancy of nodes in cases with high k values was excessive and would result in extra deployment cost. Therefore, we conjectured that with the increment of k, the optimal θ(k, η) was no more than O(k) or even had an upper bound of a constant.
Proof. It is difficult to give the strict mathematical proof of the conjecture, so we discuss the problem by performing statistical simulations.
Let η be 89%. We ran the OθSA (Optimal θ Search Algorithm) simulation experiments with k values varying from 1 to 300. We had the maximum of θ(k, η) as max{θ(k, 89%)} k=1∼300 = 2.3 when k = 5. Namely, for k = 1 ∼ 300, the value of optimal θ(k, η) had a numerical upper bound of 2.3. The simulation result shown in Figure 8 also suggests that θ(k, η) decreased when k increased. In general, the value of θ(k, η) decreased when k increased from 10 to 250.
For higher k values, the node deployment interval (l = 2r/(m · 3 √ k)) was rather small, which led to a high time complexity of OθSA. Therefore, we only carried out the OθSA for several individual k values to testify to our hypothesis. We obtained θ(500, 89%) = 1.6 when k = 500, θ(700, 89%) = 1.6 when k = 700, and θ(900, 89%) = 1.6 when k = 900. None of the above θ values was larger than 2.3, which was consistent with our hypothesis. Note that θ(k, η) seemed to approximate to some lower bound when k was large, and this lower bound should be greater than one. When θ was set to be one, we performed simulations with different k values and calculated the k-coverage rate. We obtained a 28.79% k-coverage rate when k = 500, a 29.96% k-coverage rate when k = 700, and a 26.11% when k = 900. The average k-coverage rate obtained was 33.32%. This suggests that θ(k, η) should still be greater than one even when k is very large. In fact, when we increased θ to 1.1, the average k-coverage rate was then 61.71%, which almost doubled the k-coverage rate when θ = 1. This shows that the redundancy caused by θ plays a crucial role in the improvement of k-coverage. Even a slight increment of θ(k, η) over one could significantly increase the k-coverage rate.
The value of θ(k, η) was within a small constant range and had an upper bound, independent of k. It can be seen that the actual deployment volume density ensured the coverage was linearly related to k, and the deployment cost was acceptable even if the value of k was large. To study the influence of the value of r on the performance of k-ERVFA, we conducted simulation experiments with different r values. Most of the simulation parameters were the same as those in Table 5, except for the coverage radius r and the number of sensor nodes.
According to Equation (14), the total number of sensor nodes needed to obtain a satisfactory k-coverage rate (89% in this paper) is proportional to 1/r 3 . Figure 9 shows the minimum number of nodes needed in each sub-region to meet the corresponding k-coverage (k = 1, 2, 3) requirement with r varying from 5 m to 20 m.
As shown in Figure 9, with the increase of the node coverage radius r, the number of required nodes decreased rapidly at first and then became stable. This was because when the coverage radius was small, the coverage area was very limited, and a large number of nodes were needed to achieve the required coverage rate. As r increased, the coverage rate increased, and the number of nodes required decreased. When r reached a certain threshold, the number of required nodes tended to stabilize due to comprehensive coverage and network connectivity.
In addition, the coverage radius of the node was determined by its power. The increase of power would increase the energy consumption of the system and reduce the network life, so it was necessary to consider comprehensively the coverage and energy consumption of the system to choose the appropriate node radius in the practical applications.
In the case of r = 1 m, the calculated number of sensor nodes was: 5.9041 * 10 5 in the one-coverage requirement sub-region, 61,115 in the two-coverage requirement sub-region, and 38,675 in the three-coverage requirement sub-region. Note that the total number of sensor nodes needed when r = 10 m was 692 (seen in Section 5), and it is obvious that with different coverage radius, the number of sensor nodes needed differed greatly. Fewer sensors were needed when the coverage radius was large, which is in consistent with practical experience. Figure 10a-c show the curves of the k-coverage rate (k = 1, 2, 3) for different sensor numbers (400 and 500) when r varied from 5 m to 15 m.
It should be noted that when r was small (e.g., r = 5 m), the calculated total number of sensor nodes needed to obtain a satisfactory k-coverage rate was quite large (more than 5000), while we set n as 400 or 500 in the simulations. This is the problem of node deployment in sparse sensor networks. Because k-ERVFA adopts a greedy strategy that considers node deployment in high k-coverage requirement sub-regions as the first priority, the coverage rate in the low k-coverage requirement sub-regions achieved in sparse sensor networks was rather poor. Therefore, the Ak-ERVFA (seen in Section 6.4.1) was applied here to improve the performance in sparse sensor networks, especially for sub-regions with a low k-coverage requirement.
The coverage radius of the sensors (m)  For higher r values (r ≥ 9 m), we tended to use the original k-ERVFA algorithm to highlight its advantage of improving the k-coverage rate with high k values.
As shown from Figure 10a-c, the k-coverage rate in the corresponding area increased with the value of r, which is consistent with practical experience. When r = 5 m and n = 500, which is the aforementioned case of a sparse sensor network, the initial 1-, 2-, and 3-coverage rates in the corresponding sub-regions were 19.06%, 0.96%, and 0.37%. These coverage rates increased to 23.17%, 3.34%, and 1.68% in the simulation experiments of Ak-ERVFA.
The two-coverage rate and three-coverage rate increased significantly when 5 m ≤ r ≤ 10 m, as shown in Figure 10b,c. This increase tendency slowed down a bit when r was greater than 10 m. In particular, when r = 12 m and n = 500, the initial 1-, 2-, and 3-coverage rates were 94.36%, 89.7%, and 57.13% and were improved to 96.02%, 98.99%, and 97.89% by k-ERVFA. These k-coverage rates were rather impressive compared with the k-coverage rate obtained in Section 5 with 650 sensor nodes, considering only 500 sensors were deployed here. This indicates that a slight increment of r can lead to a significant improvement of the k-coverage rate. When r = 15 m and n = 500, the initial 1-, 2-, and 3-coverage rates were 99.13%, 99.94%, and 97.51%, and they all increased to 100% after redeployment with k-ERVFA. In conclusion, the value of r had a significant influence on the variation of the k-coverage rate.
Obviously, the k-coverage rate achieved by k-ERVFA was 100% for higher r values (r ≥ 15 m). Therefore, for the diverse 3D underwater k-coverage requirement scenarios described in this paper, 400∼500 sensor nodes with a physical coverage radius of 10 m∼15 m can guarantee a satisfactory k-coverage rate.

Discussion of the Value of N max
In this section, we study the influence of the maximum iteration number on the performance of the k-ERVFA algorithm. Similarly, the setting of the underwater monitoring area and all simulation parameters except for N max were the same as those in Section 5. We carried out simulations with n = 400 and N max increasing from 10-150 by steps of 10. Figure 11 shows the different k-coverage rates obtained by k-ERVFA with different values of N max . It should be mentioned that in order to eliminate the difference from the initial coverage rate caused by the initial random deployment, the initial sensor distributions in different rounds were set to be identical. The initial 1-, 2-, and 3-coverage rates were 77.51%, 41.96%, and 24.81%, respectively. Besides, we modified Step 2 of the original k-ERVFA so that the loop ends only if the iteration number has reached N max , i.e., the judging condition of whether the coverage rate has reached η was removed. The "fix" operation was also abolished; thus, sensor nodes moved more freely in order to demonstrate the effect of N max and the greedy strategy of k-ERVFA.
As expected, the two-coverage rate and three-coverage rate increased significantly with the growth of N max , as shown in Figure 11. When N max = 150, the two-coverage rate obtained was 91.33%, and the three-coverage rate was 83.17%. They increased by 49.37% and 58.36%, respectively, compared to the initial coverage rate.
The one-coverage rate increased smoothly when N max increased from 10 to 50, reaching a maximum coverage rate of 86.35% when N max = 50. However, it decreased gradually when N max was greater than 50 and reached a minimum of 77.96% when N max = 150. It is reasonable to predict that the two-and three-coverage rate will keep growing, while the one-coverage will decrease when N max keeps increasing. In fact, when N max was small, the homogenization effect from repulsion was the dominant factor to facilitate the increase of the one-coverage rate. When N max was large, the k-attraction of the sub-regions with the high k-coverage requirement was dominant, and more sensors moved into high k-coverage sub-regions, which led to the increment of the high k-coverage rate and the decline in the one-coverage rate. Larger N max also highlighted the characteristics of the greedy strategy adopted by the k-ERVFA algorithm, as mentioned above.
According to Figure 11, the two-coverage and three-coverage rate increased significantly when N max grew from 40 to 100, and this increasing tendency slowed down when N max was over 100. When N max was around 100, the one-coverage rate was also acceptable. Considering that the excessively large value of N max would result in unnecessary running time and energy consumption, it was appropriate to set N max as 80 to 100 in the underwater monitoring scenario described in this paper. In the practical scenario, an appropriate value of N max is important to achieve the better trade-off between coverage requirements and network cost.

Discussion of the Value of ∆L max
In this section, we study the effect of ∆L max , which is the maximum possible moving distance (along the Z-axis) of the sensor node in each iteration. Similarly, the parameters except for ∆L max were identical to those in Section 5. We carried out simulations with n = 400 and fixed the initial distribution of sensor nodes throughout the simulation process for the same reason mentioned in Section 6.3.2. The initial 1-, 2-, and 3-coverage rates were 77.51%, 41.96%, and 24.81%, respectively. Figure 12a shows the different k-coverage rates achieved by the k-ERVFA with ∆L max growing from 0.1 m to 1 m with a step size of 0.1 m. When ∆L max = 0.1 m, the final 1-, 2-, and 3-coverage rates were 83.72%, 58.78%, and 35.76%, respectively. Compared with the initial k-coverage rates, the improvement by k-ERVFA was inappreciable. The reason was that compared to the size of the underwater monitoring area (100 m * 100 m * 100 m), ∆L max was so small (0.1m) that sensor nodes could not reach the desirable destinations based on k-ERVFA within the maximum iteration times (N max = 100); thus, the k-coverage rate obtained was unsatisfactory. However, the effect of k-ERVFA improved with the increasing of ∆L max . The decline in the one-coverage rate shown in Figure 12a was due to the fact that more sensor nodes moved into the two-and three-coverage requirement sub-regions when ∆L max increased.
In Figure 12b, when ∆L max grew from 1 m to 20 m with a step size of 1 m, all the k-coverage rates oscillated with the growth of ∆L max . When the value of ∆L max was small (1m ≤ ∆L max ≤ 7 m), the moving range of sensor nodes in each iteration was longer compared with those in Figure 12a (∆L max < 1 m). With the same iteration number, the effect of the redeployment process was more significant; thus, all the k-coverage rates increased. However, when ∆L max kept increasing (∆L max > 7 m), the moving range of sensor nodes in each iteration was so long that their motion pattern turned into a certain kind of disordered oscillation. Hence, the k-coverage rates obtained by k-ERVFA could not satisfy the actual requirements.
In conclusion, the value of ∆L max should be moderate. If ∆L max is too small (≤1 m), then sensors cannot reach their final destinations before the redeployment process is terminated.
However, the excessively large value of ∆L max (≥10 m) will cause nodes to drift and oscillate; thus, the performance of k-ERVFA will be unpredictable. We can see from Figure 12b that when ∆L max = 7 m, the 1-, 2-, and 3-coverage rates all reached their maximum value. Therefore, ∆L max = 7 m was the optimal value under such simulation conditions (n = 400, N max = 100, r = 10 m).  Table 6 shows the optimal ∆L max for different total numbers of sensors (when r = 10 m). We can see that with the node number increasing, the optimal value of ∆L max decreased gradually. The reason was that with more sensor nodes deployed in the area, the average distance between the coverage breached, and the neighboring sensor nodes decreased; thus, the appropriate moving distance of the sensor node in each iteration (which was decided by ∆L max in Equation (22)) decreased as well. Therefore, the key to application implementation lies in how to obtain the optimal value of ∆L max according to the total number of nodes in the actual scenario.

The Improvement of k-ERVFA
The underwater sensor network is closely related to many marine environment factors, such as fish stock interference and tidal disturbances. For example, when the fish school destroy the nodes or the nodes and equipment age and fail, the previous deployment requirements will be changed to adapt to the new situation. The number of nodes may be insufficient in some scenarios, and even in some extreme conditions, the total number of sensor nodes is so small that the sensor network is quite sparse. In addition, in practical applications, k-coverage requirements may change over time according to seasonal changes or actual situations, which is also a real problem that needs to be considered. Aiming at these actual problems mentioned above, we discuss the Ak-ERVFA (Averaged k-ERVFA) algorithm for node sparsity and the Ck-ERVFA (Changed k-ERVFA) algorithm for time-variant demand, respectively.

Ak-ERVFA
Due to the strategy taken by k-ERVFA, which satisfied the k-coverage requirement with high k values preferentially and the fact that the total number of nodes was severely limited, the performance of k-ERVFA in sparse sensor networks was unsatisfactory. Especially, sub-regions with requirements of low coverage multiplicity tended to be ignored.

Modify the termination condition in
Step 2 of k-ERVFA; now, the iteration will also be terminated if the number of nodes inside A k,j reaches λ(r, k) · N total , where N total is the total number of sensor nodes.
We can see that with Ak-ERVFA, the number of nodes deployed in any sub-region was proportional to its expected node number calculated by Equation (14). In this way, sub-regions with different k-coverage requirements were treated equally to some extent. Thus, the k-coverage rate in sub-regions with the requirements of low coverage multiplicity became reasonable.

Ck-ERVFA
Ck-ERVFA was defined as the enhanced k-ERVFA for time-variant situations. In Section 4, the coverage requirement information of the underwater monitoring area was known in advance and was considered to be fixed throughout the life-cycle of the sensor network. However, the coverage requirement information is often time-variant in practical applications, and we should extend our k-ERVFA algorithm so that it may work in such time-variant situations. In fact, this extension is possible when the running time (response time) of k-ERVFA is far less than the time interval between changes in coverage requirement information, which is the case in most practical applications.
We propose the enhanced Ck-ERVFA as follows: 1. Sink nodes collect the latest coverage requirement information, i.e., the geometry attributes and the coverage multiplicities of all sub-regions. Then, the coordinates of the centroids of the sub-regions and the corresponding equivalent radius are calculated. All sub-regions are sorted in descending order of the k value. Set n as zero. 2. Sink nodes calculate the k-resultant force similar to Step 2 in the k-ERVFA algorithm. 3. Call the fix and even algorithm similar to Step 3 in the k-ERVFA algorithm. 4. Sink nodes send redeployment commands to sensors, and let n = n + 1. Check if n reaches k n , otherwise go to Step 2. 5. Sink nodes monitor the coverage requirement information constantly; go to step 1 if the coverage requirement information changes, otherwise go to Step 5.
In summary, Ck-ERVFA monitors changes of the coverage requirement information, and node redeployment will be performed according to the latest coverage requirement. In this way, k-ERVFA can adapt to underwater monitoring tasks with the time-variant coverage requirement. It should be noted that in Ck-ERVFA, sink nodes are capable of collecting information about the coverage requirements, as discussed in [28].

Conclusions and Future Work
In this paper, we studied the diverse k-coverage problem in three-dimensional underwater sensor networks. First, we analyzed the minimum density of sensors required to satisfy the different k-coverage requirements by dividing the monitoring area into grids. We introduced θ to provide the necessary redundancy of sensor nodes required in practical applications and proposed the OθSA algorithm to decide the appropriate value of θ. Then, we put forward the enhanced virtual force algorithm k-ERVFA to solve the diverse k-coverage problem. Both theoretical analysis and simulation experiments were carried out to demonstrate the effectiveness and feasibility of our proposed algorithm. Finally, we extended the k-ERVFA to Ak-ERVFA and Ck-ERVFA so that it could work in sparse sensor networks and time-variant situations.
Moving forward, we plan to combine our approach with methods based on the Voronoi diagram to better discover coverage breaches and further improve the diverse k-coverage rate in 3D UWSNs. We also plan to implement the proposed algorithm with a testbed to evaluate the performance of k-ERVFA in real-world application scenarios while taking the influence of marine environment factors and measurement error into account. In this subsection, we perform a theoretical analysis on the correctness and effectiveness of the k-ERVFA algorithm when k = 2.
According to the k-ERVFA algorithm, spheres with a two-equivalent radius (r 2 = r/2) are deployed in the two-coverage requirement sub-regions. Note that the actual coverage radius is still r, as mentioned in Section 3.2.
As shown in Figure A1, the larger spheres are the actual coverage region of the sensor nodes, while the smaller spheres are the two-equivalent radius spheres defined in Definition 6 (seen in Section 3.2). In Figure A1, the two-equivalent radius spheres are tangent to each other, which is considered as the critical case. When adopting the k-ERVFA algorithm, neighboring two-equivalent radius spheres tend to overlap with each other, which will result in a higher two-coverage rate than the critical case. For simplicity, assume that the sub-region is a cube with side length L (L ≥ 3r). If L is less than 3r, the number of nodes needed can be calculated by Equation (14) in the paper, then these sensor nodes are deployed into the cubic region. In this case, there will be overlapping areas between adjacent two-equivalent spheres, and the actual coverage rate will be greater than the critical case. Based on the discussions above, here we only need to study the critical case in Figure A1 when L ≥ 3r.
There are four kinds of three-dimensional regions that need to be discussed separately. These regions are marked as α, β, ε, and θ in Figure A1, where α ∩ β ∩ ε ∩ θ = ∅.
1. For any point within region α, it is covered by four sensor nodes because it is within the coverage range from Sensor 1 to Sensor 4. Namely, region α is at least four-covered. 2. For any point within region β (β includes the internal part and upper external part of the sphere with the two-equivalent radius of Sensor 4; here, the upper external part is inside the actual coverage sphere of Sensor 4), it is covered by Sensors 3 and 4 at least. Hence, region β is at least two-covered. 3. For any point within region ε (ε is located near the edge in the right front side of the cube), it can only be covered by Sensor 2. Thus, region ε is only one-covered. We need to know the total number of ε regions (denoted by n ε ) in the entire cubic in order to calculate the 2-coverage rate. Assume that L = K · r, K ≥ 3, K ∈ N + ; we derive the total number of region ε as follows (the detailed derivation process can be found in Appendix A.2 for k = 3): 4. For any point within region θ (θ is located in the vertices of the cube as shown in Figure A1; θ includes the external part of sphere with the two-equivalent radius of Sensor 1; here, the external part is inside the actual coverage region of Sensor 1), it can only be covered by Sensor 1; thus, region θ is only one-covered. There are eight θ regions in the cube, and each one corresponds to a vertex of the cube.
In order to calculate the two-coverage rate of the entire cubic, we need to estimate the volume of region ε and region θ, which are not two-covered.
The front, bottom, and right surfaces of region θ are flat, while the top, left, and back surfaces are spherical, which are formed by some parts of the spheres with a two-equivalent radius of Sensor 2, 4, and 6.
Based on the symmetry of region θ, we can estimate its volume as follows: where l θ 3 is the volume of the cubic with side length l θ and S is the area of the square with side length l θ (these squares share common vertices with the three spherical surfaces of region θ). 1 3 S (l θ − r 2 ) is the volume of the three pyramids that will be removed to estimate the volume of region θ. These pyramids are based on the aforementioned squares, and their apexes are at the center of the spherical surfaces of region θ. V θ is less than the value on the right side of Equation (A2) because the spherical surfaces are concave. With knowledge of solid geometry, we have l θ = 3− √ 2 2 r. There is a total of eight θ regions in the entire cube. Therefore, the percentage of V θ in the entire cube is: Similarly, we can estimate the volume of the region ε (denoted by V ε ) as follows: where l 2 θ l ε is the volume of the cuboid with height l ε . The terms after the minus sign are the volume of pyramids to be removed. Similarly, V ε is less than the value on the right side of Equation (A4) due to the concavity of region ε. From Equation (A1), we have obtained the total number of regions ε; therefore, the percentage of V ε in the entire cube is: The two-coverage rate of the entire cube is then: P reality,2 is the two-coverage rate obtained by k-ERVFA with enough sensor nodes and sufficient running time. P critical,2 is the two-coverage rate in the critical case in Figure A1 when spheres with two-equivalent radius are tangent to each other and L = 3r.
Thus far, we have proven that the k-ERVFA algorithm achieved over an 85.74% two-coverage rate in and underwater cube with side length L = 3r when there were enough sensor nodes and sufficient running time. The two-coverage rate increased with the increase of L. Especially, we achieved over a 98.31% two-coverage rate when L = 10r.
There is a similar analysis when k = 3 is presented in the next subsection (a 96.73% three-coverage rate was achieved). For k values greater than four, we propose a different approach to prove that the 100% k-coverage rate can be achieved.
When k = 1, the k-ERVFA algorithm will degenerate to the general virtual force algorithm for 3D UWSNs, for which the effectiveness and correctness have been extensively studied [7,12]. Now we analyze the three-coverage rate of the entire cube in the aforementioned critical case (L > 2r). As the physical coverage radius was still r and the distance between sensor nodes was 2 · r/3, any point in the interior part of the cube was at least four-covered. Therefore, we only need to study the coverage situation in the boundary part, which was near the surfaces and vertices of the cube. There are three kinds of three-dimensional regions that need to be discussed separately. In Figure A2, these regions are marked as α, β and ε, where α ∩ β ∩ ε = ∅.
1. For any point within region α, similar to the points in the interior part, it can be covered by Sensors 1, 2, 3, and 4, so region α is at least four-covered. 2. For any point within region β (around Sensor Node 4) that is near the boundary surface and not close to the edges of the cube, it is covered by Nodes 1, 3, and 4, so region β is at least three-covered. 3. For any point within region ε that is near the edges of the cube, it can only be covered by Node 1 or both Nodes 1 and 2, so region ε is not three-covered. We need to estimate the volume of region ε and the total number of such regions to calculate the three-coverage rate of the entire cube.
First, we estimated the volume of region ε. Figure A2b is the projection view of Figure A2a on the front surface of the cube. The radius of the circular section of Sphere 4 on the surface was then For the chord length l in Figure A2a, we have: The side length of the bottom surface of region ε is: We study the boundary situations with L = 2r in Figure A2b; thus, we obtain: Let h ε be the height of region ε on the edge; we see that the top vertex of region ε is the intersection point of the edge and Section Round 5. With the coordinate system shown in Figure A2b and knowledge of analytical geometry, we can derive the expression regarding Section Round 5 as (x − 5 3 r) 2 + (y − r 3 ) 2 = ( 2 √ 2 3 r) 2 . Let y be zero, and h ε is calculated as: Region ε is a concave polyhedron (similar to V θ in Figure A1) with concave surfaces, so its volume V ε is less than the volume of the corresponding cuboid, which has the bottom area of l ε 2 and height h ε .
Now, we consider the total number of regions like region ε. Note that region ε appears on the edges of the cube repeatedly with an interval length of r. The cube has 12 edges, so there are 12L/r of such regions in total.
Therefore, the total volume of ε regions is estimated as: Now, we consider the entire cubic region and set its side length as L = M · 2r/3, M ≥ 3, M ∈ N + , so the percentage of the total volume of all ε regions in the entire cubic region is: So far, our discussion has been based on the critical situation where spheres with k-equivalent radius are tangent to each other. The three-coverage rate of the critical case can be denoted by P critical,3 . In practice, with enough sensors and sufficient running time, the spheres with k-equivalent radius will overlap with each other so the three-coverage rate achieved by k-ERVFA (denoted as P reality,3 ) will be better than P critical,3 , i.e., P reality, 3 ≥ P critical, 3 = 1 − P ε >> 1 − 0.0327 = 0.9673 (A14) It can be found that in Equation (A13), P ε decreases when L increases. Therefore, the three-coverage rate in larger underwater areas will be greater. For example, when M = 10 (L = 20r/3), we have P ε ≤ 0.0029 and P reality,3 ≥ 0.9971. In the case of k ≥ 4 (k ∈ N + ), suppose that sensor nodes are uniformly deployed in a 3D underwater cube with side length L.
Without loss of generality, let L be the integral multiple of 2r k , i.e., L mod (2r k )= 0, where r k = r/k is the k-equivalent radius. We study the critical case where spheres with the k-equivalent radius are tangent to each other, as shown in Figures A1 and A2. Therefore, the total number of sensor nodes in the cube is (L/2r k ) 3 = L 3 /8r 3 k = L 3 /8(r/k) 3 = k 3 L 3 /8r 3 , and the average volume density of node deployment is ρ = k 3 L 3 8r 3 /L 3 = k 3 8r 3 . For some point p in the cubic region, it is covered by sensor S i only if the Euclidean distance between p and S i is less than the coverage radius of S i , i.e., d pS i ≤ r. We define the under-covered region of point p (denoted by R p ) as follows: 1. For any point in R p , its Euclidean distance to point p should be less than r; 2. The density of node deployment ρ in R p should be a non-zero constant.
Point p is k-covered if there are more than k sensor nodes in R p , i.e., ρ · V R p ≥ k. Based on this definition, we consider the k-coverage rate in the aforementioned cubic region. The problem can be divided into two different cases: 1. For any point p in the interior part of the cubic region, R p is a sphere center at p with the radius of r. The average number of sensor nodes in R p is then: 4 3 πr 3 · ρ = 4 3 πr 3 · k 3 8r 3 = πk 3 6 , so point p is πk 3 6 -covered. Let f (k) = πk 3 6 − k; we have its derivative as f (k) = π 2 k 2 − 1 > 0 when k > 1. Note that f (2) = 4π/3 − 2 > 0, and we have f (k) > 0( πk 3 6 > k) for ∀k ≥ 2; thus, point p is at least k-covered. Therefore, in the critical case where spheres with k-equivalent radius are tangent to each other, the k-coverage rate of the interior part of the cubic region is 100% when k ≥ 2. However, the discussion above is based on the average density, and the coverage rate in some sub-regions may not reach 100% in practice due to the non-uniform distribution of sensor nodes. 2. For any point p that is near the boundary surfaces and vertices of the cubic region, according to the definition, R p is no longer a complete sphere in this case. Here, the sphere is cut by the surface of the cube, and there are no sensors in the part outside the cube. Therefore, the volume of R p is now m · 4 3 πr 3 with 1/8 ≤ m ≤ 1. In particular, m = 1/8 is the case when p coincides with the eight vertices of the cube. The average number of sensor nodes in R p is then m · 4 3 πr 3 · ρ = 4mπr 3 3 · k 3 8r 3 = mπk 3 6 . Similarly, let f (k) = mπk 3 6 − k, and we get f (k) = πm 2 k 2 − 1, 1/8 ≤ m < 1. Let m be 1/8, which corresponds to the minimum of f (k) and f (k), and we have f (4) = 0.1888 > 0 and f (4) = π − 1 > 0. As f (k) and f (k) increase with k, it can be proven that f (k) > 0( mπk 3 6 > k), when k ≥ 4 and m = 1/8. Namely, point p can be k-covered even if it coincides with the eight vertices of the cubic region. For other cases with various m values (1/8 < m < 1), the volume of R p is greater than the case with m = 1/8, which means there are more sensor nodes in R p . Therefore, point p is also at least k-covered.
Based on the discussions above, points in the interior part of the cubic region are k-covered when k ≥ 2, and points that are near the boundary surfaces and vertices of the cubic region are k-covered when k ≥ 4. Theoretically, the k-coverage rate of the cubic region will reach 100% when k ≥ 4.
In fact, for the case of k ≥ 4, the same result can be obtained with the approach given in Appendix A.2.
In conclusion, with enough sensor nodes and sufficient running time, the k-ERVFA algorithm will achieve the following k-coverage rates in the critical case where spheres with k-equivalent radius are tangent to each other: • Over an 85.74% two-coverage rate for the two-coverage requirement sub-regions; • Over a 96.73% three-coverage rate for the three-coverage requirement sub-regions; • A 100% k-coverage rate for the k-coverage requirement sub-regions when k ≥ 4, k ∈ N + .

Appendix B. Selection of the Underwater Deployment Model
The deployment model we adopted in this paper is the 3D underwater sensor deployment model proposed in [45]. We also adopted the system proposed in [46], which can adjust the depth of sensors automatically. In general, a cloud-based architecture was proposed to provide a high quality of service for UWSNs. However, UWSNs also require low latency, and the central servers need to store and transmit tremendous data. Edge computing can overcome these drawbacks of traditional cloud computing. Therefore, we adopted the edge computing architecture, which performs the related tasks of UWSNs at the nearby edges of networks. The details of our model are shown in Figure A3. Initially, all sensor nodes were scattered into the target sea area from aircraft or ships, and they were then fixed by anchors on the seafloor. Rigid cables were used to prevent buoys from drifting, which would cause sensors to deviate from the target sea area. The length of the cable was adjustable; thus, the depth of the sensor node could be adjusted.
As for the transmission medium, the models proposed by [37,47,48] used an underwater acoustic link. However, due to the limitations of current technology, the underwater acoustic link suffers from issues like low bandwidth, high latency, path loss, noise [49,50], multi-path, etc. Here, we selected the wired cable as the transmission medium between underwater sensor nodes and buoy nodes. Buoy nodes were equipped with wireless communication modules, which received data collected by sensor nodes and transmitted the data to the edge servers. Edge servers became responsible for a subset of tasks in their close vicinity. The buoy nodes also acted as relays for communications between sensor nodes and were equipped with a positioning module (such as a GPS receiver) and an energy supplement module (such as solar panels). Figure A3. Underwater sensor deployment model.