This study considers a three-dimensional UWSN containing N anchor nodes with known locations and a target whose position needs determining. Suppose that the position of the i th anchor node is a i = [ a i 1 , a i 2 , a i 3 ] T , where i = 1 , ⋯ , N  and the target position is

x = [ x 1 , x 2 , x 3 ] T .

T represents the transpose operation. This study assumes that the target can transmit an acoustic signal with RSS information to the anchors (Mei et al., 2020a).

where P r i denotes the received signal power of the i th anchor node from the target, P 0 is the transmit power, α represents the path loss exponent, d i is the distance between the i th anchor and target, d 0 is the reference distance (1 m), γ i is the shadowing noise modeled as a Gaussian distribution with zero mean and variance σ i 2 , and α f is the absorption coefficient obtained from Thorp’s formula (Stojanovic and Preisig, 2009) according to

where f is the frequency.

A standard RSS-based terrestrial localization scheme is developed using α f = 0 . In practice, the frequency at which the acoustic system operates in UWSNs varies from 10 kHz to 1000 kHz, where the corresponding α f can be 0.001 dB/m to 0.32 dB/m (Stojanovic and Preisig, 2009).

In contrast to terrestrial environments, acoustic signals do not propagate straightly owing to stratification. The travel length is prone to an arc, as shown in Figure 2 (Yan et al., 2021). Moreover, the localization accuracy may degrade if directly using the straight propagation line (Mei et al., 2022b). Accordingly, a compensated model using Snell’s law is used (Zhang et al., 2020) such that

where θ x and θ i are the corresponding angles to the target and i th anchor, respectively, k is a constant, z represents the depth information, and

C ( · ) is the function of sound speed following the isogradient model (Su et al., 2022) such that

where ϑ is the steepness of the sound speed profile (SSP), and b is the speed of sound on the ocean surface.

From Figure 2 , we also have

where r denotes the horizontal range following r = r x − r i = ( x 1 − a i 1 ) 2 + ( x 2 − a i 2 ) 2 , l denotes the ray-tracing length, and t represents the travel time.

The following expression is obtained by exploiting derivates with respect to

z and

θ , i.e

Taking the integral and replacing (6) with (5), the underwater ray-tracing length, considering an isogradient SSP, can be expressed as

where θ x and θ i are further expressed as θ x = β 0 + α 0 and θ i = β 0 − α 0 with β 0 = arctan [ ( a i 3 − x 3 ) ( r x − r i ) ] and β 0 = arctan [ 0.5 ϑ ( r x − r i ) [ b + 0.5 ϑ ( a i 3 + x 3 ) ] ] , respectively.

The compensated model is then derived, and the distance between the

i th anchor and target can be acquired using the ray-tracing length. To this end, the RSS-based measurement model in (1) is then further transformed into

Let c i = α f l i , and provide the observation vector P = [ P r i ] T . Then the probability density function (PDF) can be expressed as

By maximizing the PDF, the maximum likelihood (ML) estimator can be derived as

where x ^ is the estimate.

However, solving the problem in (10) accurately is challenging owing to its high nonconvexity. Thus, the problem is then transformed into an ANCLS framework.

Let σ i = σ . By applying the transformation and linearization manipulation, the problem in (10) can be converted into

where κ i = P 0 − P r i − c i 10 α , and ∥ · ∥ is the ℓ 2 norm.

After squaring the range, (11) can be rewritten as

where χ = ∥ x ∥ 2 .

Given

ξ = [ x 1 , x 2 , x 3 , χ ] T as the estimated parameter, the localization problem can be transformed into

where

With the constraint of ξ , i.e., ξ ≥ 0 , the original problem in (10) can be reformed into the ANCLS framework, as shown in (13), with a single right-hand side (RHS) vector (Kannan et al., 2018).

Assuming that the matrix W has a full column rank, the problem in (13) is strictly convex because

W T W is positively definite. An optimal solution is acquired when it satisfies the Karush–Kuhn–Tucker (KKT) optimality condition, expressed as

where M is the number of estimated variables (blocks) in ξ .

To obtain a feasible solution in (13), the index set { 1 , ⋯ , M } is split into two groups, ℋ and ℂ , such that ℋ ∪ ℂ = { 1 , ⋯ , M } and ℋ ∩ ℂ = ∅ . In addition, ξ ℋ , ξ ℂ , y ℋ , and y ℂ are defined as subsets of variables with corresponding indices and W ℋ and W ℂ are the submatrices of W with corresponding column indices.

At the initiation, let

ℋ = ∅ ,

ℂ = { 1 , ⋯ , M } , and

ξ ℂ = y ℋ = 0 . Accordingly, (17) holds for any values of

ξ ℋ and

y ℂ . A complementary basic solution

( ξ ℋ , y ℂ ) that can be acquired by (18) and (19) is feasible only when

ξ ℋ ≥ 0 and

y ℂ ≥ 0 . Otherwise,

( ξ ℋ , y ℂ ) is infeasible.

If

( ξ ℋ , y ℂ ) is infeasible, an updated procedure must be manipulated by exchanging variables for which (16) does not hold. This study defines an index set as

Accordingly, a variable

ξ j with

j ∈ ℱ is an infeasible variable. The updated rule of

ℋ and

ℂ can be expressed as

where ℛ is the nonempty subset, such that ℛ ⊂ ℱ .

Notably, a cycle or failure to find an optimal solution may occur during the updated procedure. Therefore, a backup rule should be adopted to ensure a finite termination, which can be expressed as

In addition, this study defines variable Δ as a buffer in the updated procedure. The value of Δ is reduced by 1 if the updated round increases the number of infeasible variables. When Δ is zero, the backup rule is used, such that the number of infeasible variables is smaller than the lowest value and stored in β . The procedure for coarse estimation is as follows: |F| denotes the number of variables at each iteration, and Δ = 3 is suggested by (Cimini and Bemporad, 2017) to accelerate the convergence.

Although CEA can obtain a relatively feasible solution, it can easily drop to a local optimum owing to its nature. This section proposes a Taylor series expansion-based algorithm to refine the solution further.

Recalling the problem in (11), the equation is expanded, and the constant item is ignored.

(24) can be further converted into a surrogate function, such that

where f ( x ) is convex, and g ( x ) is concave.

Theoretically, the concave function should have an upper bound based on the first-order Taylor expansion (Zhang and Zheng, 2010). Assuming that

x ℘ is the estimate in the

℘ th iteration,

where ∇ g ( · ) denotes the gradient.

Substituting (26) for

g ( x ) in (25) yields

Notably, the surrogate function is bounded below and can converge to a limit point (Sun et al., 2017). Therefore, given the ℘ th location estimate, the refined solution can be obtained by an iterative operation, such that

where the initial guess is obtained through CEA; that is x 0 = x ^ CEA . The flowchart of the proposed method is depicted in Figure 3 , while the pseudocode of the CFLM is summarized as follows:

The CRLB using the RSS-based technique under the stratification effect was derived as a benchmark for the unbiased estimator (Sengijpta, 1993). Given the Gaussian noise assumption, such that σ i = σ , the CRLB is expressed as

where Tr is the trace function, and F I M is the Fisher information matrix.

From (1), the partial derivatives can be expressed as

Inserting (6) into (5) and taking the integral operation yields

The following equations are derived after obtaining the partial derivatives of (33) and (34) with respect to

x 3 and

r x , respectively:

According to (7), the derivatives of the ray-tracing length with respect to

x are

Where

Following (30)–(43), ∂ P ∂ x can be further expressed as

where ϖ = 10 α ln 10 .

Let

The CRLB using RSS-based technique under the stratification effect is expressed as

Figure 4 shows the results of cumulative distribution function (CDF) when σ i = 3  dB , s i d e l e n g t h = 50  m , and ϑ = 0.1 . Figures 4A–D ) show the schemes with N = 4 , N = 6 , N = 8 , and N = 10 , respectively. Notably, the number of available measurements increases as the anchor grows. In this case, the localization accuracy improves, as shown in Figure 4 . For example, the probability of ∥ x ^ − x ∥ ≤ 5.43  m reaches 95% for the CFLM in Figure 4A ), where the CFLM error is 3.61 m at the same probability as in Figure 4D ). The same situation is observed for WLS, PPSL, CEA, PSO, DE, and TLLA-RSS. Table 2 summarizes the corresponding errors and probabilities for Figures 4A, D ). Among these methods, the CFLM performs better than the WLS, PPSL, CEA, PSO, DE, and TLLA-RSS. Its performance is close to that of CRLB, although the margin remains remarkable when N > 4 .

An error analysis is performed to further evaluate the performance of each Monte Carlo trial, as shown in Figure 5 . Figures 5A, B ) show the scenarios with N = 6 and N = 10 , respectively. As shown in Figure 5 , outliers can occur in some situations. However, the CFLM performs better than the other methods. Its third quartile, median, and first quartile errors in Figure 5A ) are 3.12 m, 2.51 m, and 1.96 m, respectively. This outperformance is remarkable with N = 10 . The median errors of the CFLM, TLLA-RSS, CEA, PPSL, PSO, DE, and WLS in Figure 5B ) are 2.14 m, 2.54 m, 2.71 m, 2.84 m, 2.81 m, 2.55 m, and 3.33 m, respectively. Moreover, the corresponding results show that the proposed CEA obtains a better solution than the WLS and PPSL in some situations. The refined step proposed in Section 3.3 can further improves the localization accuracy.

Figure 6 shows the CDF of variable noises under N = 8 , s i d e l e n g t h = 50  m , and ϑ = 0.1 . Figures 6A–D show the schemes with σ i = 1  dB , σ i = 3  dB , σ i = 5  dB , and σ i = 7  dB , respectively. The measurement error increases with σ i . Therefore, the performance degrades with increased noise. For example, the probability of ∥ x ^ − x ∥ ≤ 2.15  m reaches 95% for the CFLM in Figure 6A ), exhibiting a 3.86 m error at the same probability in Figure 6B ). It is worth noting that the performance is close between CFLM and PSO when the noise enlarges. However, the PSO localization achieves the accuracy at cost of numerous particles, which increases the computational complexity. On the contrary, the proposed method can get a satisfactory performance at a relatively low computational complexity. Table 3 lists the results. The CFLM outperforms the other methods, which is within 3.86 m when σ i = 3  dB , reaching 95% compared with that of the WLS (5.47 m), PPSL (5.30 m), CEA (4.73 m), PSO (4.35 m), DE (4.16 m), and TLLA-RSS (4.35 m). However, the margin between the CFLM and CRLB remains remarkable.

Figures 7A, B ) show the scenarios with σ i = 3  dB and σ i = 7  dB , respectively. The PPSL is sensitive to noise. Its performance fluctuates significantly compared with that of the other methods. Moreover, the outliers of the PPSL increases with noise. The median errors of the CFLM when σ i = 3  dB and σ i = 7  dB are 2.24 m and 3.49 m, respectively, better than that of the other methods.

Figure 8 shows the CDF of the variable side lengths with N = 8 , σ i = 3  dB , and ϑ = 0.1 . Figures 8A–D ) show the schemes with the sidelength at 25 m, 50 m, 75 m, and 100 m, respectively. Notably, the signal attenuation becomes severe with enlarged area of interest. Accordingly, the localization error increases with the sidelength. Among the evaluated methods, the CFLM reaches ∥ x ^ − x ∥ ≤ 2.73 m and ∥ x ^ − x ∥ ≤ 5.48 m at 95% with sidelength at 25 m and 100 m, respectively. By contrast, the other methods could not reach the same localization accuracy. Table 4 lists the performances of the methods, whereas Figure 9 shows the corresponding box plot. The results in Figure 9 ) further proves the superior performance of the CFLM. Its median errors are 1.58 m and 3.17 m, as shown in Figures 9A, B ), respectively.

Steepness depends on the ocean environment characteristics, such as the temperature (Ramezani et al., 2013). Accordingly, simulations are performed to evaluate whether the proposed method still works in variable steepness. Figure 10 shows the results with N = 8 , σ i = 3  dB , and s i d e l e n g t h = 50  m . The WLS, PPSL, CEA, TLLA-RSS, PSO, DE, and CFLM are robust to the steepness, as shown in Figure 10 and Table 5 . Conversely, the CRLB is sensitive to the change in steepness. Its error increases with the steepness. Moreover, the CFLM performs better than the other methods, regardless of the changes in steepness. The CFLM achieves ∥ x ^ − x ∥ ≤ 3.87 m and ∥ x ^ − x ∥ ≤ 3.89 m at 95% with ϑ = 0.1 and ϑ = 0.7 , respectively. The box plot of the deviation in Figure 11 also proves the superior performance of the CFLM. Its median errors are 2.24 m and 2.26 m, as shown in Figures 11A, B , respectively.