Abstract
Underwater acoustic localization is an important, yet challenging problem: (1) node mobility issue, (2) Doppler effect, and (3) clock imperfection. To be specific, underwater nodes are not stationary in real-life due to unpredictable currents. Relative motion between a transmitter and a receiver causes the time scaling problem on the received signals, where the time scaling factor is termed as Doppler scale. Then, due to the slow acoustic signal propagation speed, the underwater Doppler scale becomes more severe compared with the one in terrestrial environments. Thus, the differential Doppler scale (DDS) measurements should also be collected, other than the time measurements like time-difference-of-arrival (TDOA), for enhancing the underwater localization. Since DDS/TDOA measurements and clock skew are tightly coupled, clock synchronization is essential for accurate localization. However, due to the stringent cost and power constrains of underwater nodes, low-cost clocks with relative low precision are normally employed, which makes it even more difficult to guarantee a perfect clock synchronization between transmitter/receiver pairs. In order to cope with those issues, we propose an algebraic underwater localization method using the hybrid DDS/TDOA measurements, which is particularly robust against the node clock imperfection. A new DDS/TDOA measurement model with clock imperfection is first presented by analyzing the received signals over underwater acoustic channels. Then, we devise a two-step weighted least square-based estimator, and the analytical study shows that our estimator can achieve the Cramer-Rao lower bound (CRLB) accuracy under small noise. Simulations corroborate the theoretical results and the good performance of the proposed method.
1. Introduction
Underwater localization has been an active area in recent years owing to its extensive applications such as data collection, environment monitoring, military surveillance, and assisted navigation [1–5]. Basically, the underwater localization process follows two steps, i.e., measurement collection and measurement fusion. To be specific, the measurements between the target and the predeployed anchors with prior known locations are first collected, from which we infer the target location. Since the Global Positioning System (GPS) signals are not available in underwater scenarios due to the severe power attenuation of electromagnetic waves, we can do nothing but collecting measurements from the acoustic signals [6]. For acoustic communication, there are many modulation methods, such as frequency-shift keying (FSK), phase-shift keying (PSK), and orthogonal frequency-division multiplexing (OFDM) [7]. For the low complexity of receivers, which is required to deal with highly dispersive channels, in our study, we consider the OFDM scheme for data exchange [8]. There are different types of measurements that can be extracted from data exchanges, such as time-of-arrival (TOA), time-difference-of-arrive (TDOA) [9, 10], received signal strength (RSS), differential RSS, and angle-of-arrival (AOA) [11]. Among those, the TDOA-based method becomes the primary concern of most engineers in designing the underwater localization owing to its relative high-accuracy ranging results and the relaxation of clock synchronization requirements [12]. We should also mention here the kind of differential Doppler scale (DDS) measurement, which results from the underwater nodes’ relative motion that always exist in practice [13]. This kind of measurement will be explained in details later. With all the collected measurements, the target’s position can be estimated by different families of measurement fusion methods, which are maximum likelihood (ML), semidefinite programming (SDP) based, least squares (LS) based, and alternating direction method of multipliers (ADMM) based [14, 15]. The ML method is asymptotically optimal, but the formed ML optimization problem is highly nonlinear and nonconvex, and thus, a closed-form solutions does not exist [16]. Although it can be solved approximately by iterative methods, they involve intensive computations and cannot guarantee the convergence to the correct solution unless the initial guess is close enough to it. The SDP-based method relaxes the nonconvex optimization problem to a convex one such that a global minimum can be effectively found [17, 18]. However, this method still requires a high complexity as well as a tight relaxation to guarantee an accurate estimate. To allow low complexity implementation as well as to ensure global convergence, (weighted) LS-based methods have been proposed in [19–21]. The LS-based method rearrange the nonlinear equations into a set of linear equations by introducing extra variables, which are functions of the target parameters (position and velocity). Furthermore, the relation between the extra variables and the unknown parameters of the target can also be utilized to improve the estimation accuracy, which reaches the Cramer-Rao lower bound (CRLB) under Gaussian noise at moderate to high signal-to-noise [19]. Therefore, due to the attractive advantages of the LS-based method over other methods, we choose it as our measurement fusion method.
Besides choosing the appropriate measurement fusion method, it is worth noting that when the target is moving, the DDS measurements can be explored together with the TDOA measurements to further improve the localization accuracy. Actually, underwater nodes can hardly maintain stationary due to the unpredictable ocean currents. A relative transmitter/receiver motion results into the time scaling problem on the received signals [22]. The time scaling factor is conventionally called as Doppler scale. Moreover, the low propagation speed of acoustic waves (about 1500 m/s) makes the acoustic signals very susceptible to the Doppler scale. Thus, the DDS measurements have to be effectively extracted from the received signals and exploited for underwater localization. It is worth noting that a similar kind of measurement called frequency-difference-of-arrival (FDOA) has been extensively studied for terrestrial localization problem [19–21]. One might ask what is the difference between the DDS and FDOA measurements. In fact, the DDS is the difference in received Doppler scales in time-domain, while the FDOA is the difference in received Doppler frequency offsets in frequency-domain, and both of them results from the Doppler effect. However, from the measurement point of view, narrowband signals will have more precise FDOA measurements compared with wideband signals [23]. Though underwater channels are wideband in nature because the signal bandwidth is not negligible compared to the carrier frequency [24]. To be specific, taking OFDM signals as an example, each subcarrier experiences a Doppler-induced frequency offset, which depends on the frequency of the subcarrier. This kind of Doppler shifts is called as nonuniform Doppler frequency offsets, and directly estimating it in frequency domain is intractable [6]. As a result, precisely estimating the Doppler frequency offsets in frequency-domain from the received signals over underwater channels might be difficult. Therefore, we alternatively choose to measure the time scaling factor of received signals in time-domain, which results in DDS measurements.
It is well known that time-based localization is very susceptible to the clock imperfection that always exists in practice. For example, due to the stringent cost and power constrains of underwater nodes, low-cost clocks are normally employed. This implies that the clock parameters of underwater nodes, i.e., clock skew and clock offset, might drift away over time. Although the anchor nodes can be synchronized with a reference clock by precalibration, the target node is usually very difficult to be guaranteed the same particularly for underwater scenarios. Even though the target is synchronized with the anchors, the clock imperfection might still exist due to the fact that the synchronization performance might significantly deteriorate in severe underwater communication environments. This will certainty incur the clock imperfection for the target node. Thus, taking into account the clock imperfection in underwater localization becomes an important task. In this work, we assume an independent clock for the target node and synchronized clocks among all the anchors. Traditionally, clock imperfection are usually considered while developing TOA- or TDOA-based localization algorithm [25, 26]. Compared with the TOA-based localization, TDOA technique resolves the clock offset ambiguity, though it can still suffer from the clock skew [27]. Other than the TDOA measurements, we also use the DDS measurements for enhancing the underwater localization. Obviously, the DDS measurement provides more information for localization. As a trade-off, when considering the clock imperfection, using the DDS is at price of a complicated measurement model since the measured DDS at the receiver is actually a combination of the Doppler effect and the clock skew between the transmitter/receiver pair, which will exacerbate the nonlinearity issue. However, the coupling nature between the DDS and clock skew is often overlooked in Doppler measurement based localization approaches. Recently, some pioneering research works noticed the potential coupling relationship between the Doppler scale measurement and clock skew [28–30]. However, they only explored the time synchronization problem in underwater sensor networks (UWSNs) using the clock skew-interfered Doppler scale measurements. In a nutshell, the clock imperfection will significantly degrade the DDS/TDOA-based localization performance, and hence, the localization methods that are robust against to this imperfection are ungently required.
To tackle the aforementioned clock imperfection problem, the first contribution of this work is proposing a new DDS/TDOA model with the clock imperfection by analyzing the received signals over underwater acoustic channels. Based on this model, we then contribute to localization algorithm development by exploiting the linearization approach. Based on the pseudolinear equations with nuisance variables, two weighted least square (WLS) estimators are devised. The first one ignores the built-in relationships between the unknown parameters and gives a coarse solution. The second one improves the solution by exploiting the known relationships between the estimates resulted from the first WLS estimator. Since only two WLS estimators are involved, our proposed method is computationally attractive. We compare analytically the location accuracy of the proposed estimator to the Cramér-Rao lower bound (CRLB) for Gaussian measurement noise. We also conduct simulations to compare the performance of the proposed localization method with the corresponding CRLB and that of the WLS method assuming perfect clock [19] and the proposed method using a mismatched DDS model under different noise conditions. The numerical results verify the effectiveness of the proposed localization method.
1.1. Notations
Column vectors and matrices are denoted by bold lower- and uppercase letters, respectively; and are the th element of and the th element of , respectively; denotes a subvector with the th to the th elements of ; represents the th row of ; is the Euclidean distance norm; denotes the Dirac delta function; denotes the convolution operator; superscript denotes the transpose of a matrix (vector); denotes the element-wise multiplication; and are vectors (or matrixes) of 1 and 0, is zero matrix, denotes the identity matrix (size indicated in the subscript if necessary); and represent the diagonal and the block-diagonal matrices; denotes the expectation operator.
2. Problem Formulation and CRLB
In this section, we first show how the underwater acoustic channel and clock imperfection affect the DDS and TDOA measurement model, then the localization problem is formulated. The Cramér-Rao lower bound (CRLB) is also derived.
2.1. Problem Formulation
As depicted in Figure 1, consider a three-dimensional underwater localization scenario where moving anchors are used to determine the position and velocity of a moving target using the DDS and TDOA measurements. The anchors are predeployed in the interested monitoring area, and their positions and velocities are known to the localization algorithm as and , , respectively. Assume the anchors are synchronized and behave a common clock skew and clock offset . Thus, the local time of the anchors with respect to (w.r.t.) a universal standard time is give by [31].
We assume the local time of the target is the universal standard time, i.e., clock skew is and clock offset is .
To implement a DDS/TDOA information-based localization algorithm, it is necessary to transmit signals from the target to anchors and/or vice versa. For simplicity, we here assume that the target radiates a signal at a single time instant and received by each anchor after a propagation delay. Generally, some preprocessing steps including detection, synchronization, and Doppler scale estimation are required for underwater acoustic communication systems [24]. Several structures of the transmitted signal can be employed for the preprocessing steps, such as linear-frequency-modulated (LFM) signal, cyclic-prefixed (CP) OFDM signal, and m-sequence [32]. In this paper, as originally suggested in [33], we adopt an LFM preamble and an LFM postamble around each data frame to estimate the DDS and TDOA. The main reason for choosing an LFM signal rather than other signals is its robustness against the Doppler effect as well as its good cross-correlation performance in environments corrupted by white Gaussian noise.
Consider a multipath underwater channel between the target and the th anchor that has the impulse response [6]. where is the number of paths, is the corresponding time-varying path attenuation, and is the time-varying path delay. Assuming that the duration of transmitted signal is short enough, the relative movement between the target and anchor is small. Thus, the time variation of the path delay can be reasonably approximated by a Doppler scale as and the path attenuation are assumed constant . Moreover, we assume all paths have a similar Doppler scale , which has already been justified in [6]. Finally, the underwater acoustic channel model is approximated as7
Let be the transmitted passband signal at the target with duration , then the received passband signal at the th anchor is given by where is additive noise. Sampling yields its discrete time sequence which is used for extracting the DDS and TDOA measurements. Due to the clock imperfection, the actual sampling interval of the anchor is , where is the reference sampling interval. Thus, the received signal at the th anchor is discretized as where is the time index and denotes the upward rounding operator. From (6), we observe that the joint effect of Doppler scaling and clock skew manifests itself in scaling the signal duration from to . The scaled signal duration can be estimated by cross-correlating the received signal with the known LFM preamble and postamble, denoted as . Then, by knowing the original signal duration , the Doppler scale measurement of the th anchor is estimated as . Thus, the Doppler scale measurement model can be defined as where is modeled as the actual Doppler scale caused by the target/anchor motion and is given by
is the acoustic propagation speed, and ; is the measurement noises, which are independent and identically distributed (i.i.d.) Gaussian random variables with zero-mean and variance . Furthermore, based on the LFM preamble that inserted in the transmitted signal, each anchor is able to measure the arrival time of the first path by using matched filtering [33]. Considering the clock model defined in (1), the TOA measurement for the signal transmitted from the target at the th anchor where is the unknown start transmission time of the target and is TOA estimation error and modeled by i.i.d. Gaussian random variables with zero-mean and variance . Note that the ocean ambient noise is very complicated; and it is hard to model it accurately [34, 35]; we hence use the Gaussian noise model for the convenience of derivation. Without loss of generality, we choose as the reference Doppler scale and TOA measurement and form the DDS and TDOA measurements, respectively, as where , , and . It can be seen from (10b) that the parameters and are cancelled out by the TDOA calculation while the clock skew still affects the TDOA measurements. So far, we have elaborated the method for extracting the DDS and TDOA measurements. For better understanding, this method is illustrated in Figure 2 in which an underwater acoustic channel example with paths is used.
Stacking the DDS and TDOA measurements and fusing them into a vector, we obtain where
The measurement error vector is , and its covariance matrix is where and . Note that the Doppler scale measurement noises are assumed to be independent of the TOA measurement noises . From (11), we can write the observed DDS and TDOA equation in matrix form as where . The goal of underwater localization is estimating , , and , from the noisy measurement vector . Under the mutually independent Gaussian noise condition, the ML estimator of , , and can be formulated as [36]
The ML problem is nonconvex, implying that there exist multiple local minima, and the global minimum can hardly be obtained. We shall develop an efficient estimator by converting the nonconvex ML problem to a linear one.
2.2. CRLB
As observed in (13), the DDS/TDOA measurement vector is Gaussian distributed as . The CRLB is the lowest possible variance that an unbiased estimator can achieve. For the unknown vector , its CRLB is given by [36] where is the Fisher information matrix (FIM). Using the notation , can be calculated as
The partial derivative in (16) can be expressed as where
Based on the equations above, the FIM can be easily calculated, and hence, the CRLB is obtained.
3. Proposed Methods
The proposed method has two stages. The first stage creates a set of pseudolinear equations for the nonlinear DDS and TDOA measurements by introducing nuisance variables. Then, an initial solution is obtained through WLS optimization. The second stage utilizes the relationship between the nuisance variables and the interested parameters to refine the first stage solution.
3.1. First Stage: Transforming the Nonlinear Measurement Equation to a Pseudolinear One
To fuse the DDS and TDOA measurements, we start from transforming (10a) and (10b) as where , , , , and
Considering the noise free TDOA measurement in (19b), which can be written as
Upon rewriting (21) as , squaring both sides, and substituting for and , we have a set of TDOA equations for .
Taking the time derivative of (22) results in a set of DDS equations for .
In terms of the noisy quantities by putting and into (22) and (23) and ignoring the second-order error terms, we arrive at where is defined as the parameter vector, in which , , , and . The vector and matrix are defined as and the noise vector is defined as
Thanks to the nuisance variables , , , and introduced in the parameter vector , it makes (24) become a set of linear equations with respect to . As a result, can be estimated by the WLS method, whose solution is given as [36] where is the weighting matrix chosen as where . Note that the weighting matrix is dependent on the true values of clock skew, source position, and velocity through . To cope with this, is first employed to calculate an initial estimate from (27), which will be used back into (28) for an improved version of , then leading to a better estimate of . The estimation error in can be calculated as The covariance matrix of is, therefore, assuming small measurement noise so that the noise in can be ignored [36].
3.2. Second Stage: Refining the Estimate Obtained in the First Stage
In this stage, we shall refine the estimate obtained in the first stage by utilizing the relationship between the parameters in . In fact, they are related to each other through the following equations:
Note that there are only three independent equations among the equations in (31a), (31b), (31c), and (31d), since any one of the equations in (31a), (31b), (31c), and (31d) can be interpreted from the others. Without loss of generality, we use the relationships in (31a), (31b), and (31d) to perform the second stage estimation.
To start with, recall that the estimation error in is . Expressing and subtracting both sides by , we have where the second-order error terms are ignored. Similarly, expressing and combing with the position estimates, we have
Substituting , , , and into the equations (31a), (31b), and (31d), we obtain
In order to provide an estimation of in the second stage, we incorporates an additional equation
Staking (32)–(35), we arrive at where
On the left side of (45), the noise vector is defined as
Equation (36) is a set of linear equations with respect to ; its WLS solution is given by [36] where is the weighting matrix, which is given as
To examine the covariance of the second-stage solution , subtracting both sides of (39) by and using (36) gives
Hence, under small noise conditions, the covariance of the second stage solution is given by [36]
Finally, the clock skew, source position, and velocity estimates can be deduced from the definition of . where is used to avoid the sign ambiguity caused by the square root operation.
Similar to , the weighting matrix is also dependent on the true values of clock skew, source position, and velocity through . In practice, these true values can be substituted by the solution in and then updated by the values in (43a), (43b),and (43c). We find that iterating one or two times leads to a good solution that meets the CRLB performance.
We summarize the prototype of our proposed estimator in Algorithm 1.
| ||||||||||||||||||||||||||||||||||
4. Performance Analysis
In this section, we shall analyze the theoretical covariance matrix of the proposed solution and compare it with the CRLB. By taking the differential of defined below (45), we can relate the estimation error of (43a), (43b), and (43c) with that of as
The bias of the final solution is given by taking expectation of . Obviously, it can be seen that is linearly related to through the definitions of , , , and . Since is zero mean (when the noise is small), is also zero mean, which implies that the solution estimate is unbiased over a small noise region. Multiplying (44) by its transpose and taking expectation yields
After substituting the corresponding covariance matrices in (30) and (42), (45) becomes where
Note that (46) has the same form as the CRLB given in (15) and (16). We shall compare with the CRLB under the case of far-field target. It has been shown that comparing with the CRLB for near-field target is not easy due to the tedious form of ; however, the theoretical result drawn from the case of far-field target is also valid for near-field target in most cases [20]. For far-field target, we have the following two conditions
The first condition indicates that the distances of the target to different receivers are approximately the same since the target is very far away from the anchors. The second condition implies that the velocities of underwater nodes are ignorable compared to the distances. This is valid due to two reasons: (1) the distances are large for far-field target; (2) underwater objects usually move slowly (several meters per second).
We now evaluate the matrix . Substituting the relevant matrices into , after some straightforward algebraic manipulation and appropriately using the conditions in (48), we can show that the elements of can be approximated as where . Then, using the variable relationships that , , and and the approximations that , and , we arrive at
This completes the proof that the proposed solution in Section 3 can attain the CRLB accuracy for small Gaussian noise and far-field target. Due to the tedious form of , so far, we are not able to study the performance of the proposed solution under the case of near-field target. However, the simulation results in the next section show that the proposed solution can reach the CRLB for near-field target as well.
5. Numerical Examples
In this section, Matlab simulations are carried out to verify the effectiveness of the proposed localization algorithm (denoted by “proposed method-case 1”) by comparing with the CRLB, the WLS method assuming perfect clock [19], and the proposed method considering the combination model of clock skew-free DDS model and clock skew-involved TDOA model (denoted by “proposed method-case 2”) (Note that the case of only considering the clock skew in DDS measurements is not included for comparison. This is because that the DDS measurement is much less than the TDOA measurement (about two orders of magnitude) so that the clock skew marginally affects the DDS measurement or even be overwhelmed by the noise. Thus, the estimates under this case is of poor accuracy as we observed in the simulations. As a result, the estimation performance of this case is not included.). The algorithm derivation of the proposed method, case 2, is similar to the procedure given in Section 3, and we would not repeat the derivation for simplicity. The performance criterion is the average root mean square error (RMSE), which can be expressed as , where is the estimate of obtained in the th trial. Each simulation result is averaged over Monte Carlo trials. Two scenarios are considered in our simulation with the target being a near-field one and a far-field one. In scenario 1, 10 anchors and 1 target are randomly placed in a cube centered at , the velocities of the nodes are randomly drawn from , and the clock skew is randomly drawn from . The sound propagation speed is set as m/s. The errors and are zero-mean white Gaussian processes with identical variances of and , respectively. In scenario 2, the target is placed in while the other settings are the same as scenario 1. Scenario 1 is a near-field target case while the scenario 2 is a far-field target case. Both scenarios are used to test the RMSE performance of the considered localization methods.
Figure 3 plots the RMSE of the considered localization algorithms versus ms under scenario 1. It should be noted that perfect clock is assumed to be available in [19]. Therefore, only the clock skew estimation of our proposed method is presented in Figure 3(a). It is seen that the accuracy of the proposed method-case 1 approaches the CRLB in the whole noise range while the one that is using a mismatched DDS model (i.e., “proposed method-case 2”) is optimal only when ms. This corroborates the theoretical analysis in Section 4. As expected, the WLS method assuming a perfect clock yields the worst performance in the whole noise range. The performance gap between the WLS method assuming perfect clock and the proposed method-case 1 increases as the measurement noise increases. This demonstrates the significance of taking the clock imperfection into account. Regarding implementation complexity, the average computation times per trial for the proposed method-case 1, proposed method-case 2, and the WLS method assuming perfect clock [19] are measured as s, s, and s, respectively. This indicates the computational attractiveness of the proposed method.
(a)
(b)
(c)
Figure 4 plots the RMSE of the considered localization algorithms versus ms under scenario 2. We again see the superiority of our proposed method over the WLS method without considering the clock imperfection. However, the location accuracy is generally worse for a far-field target than a near-field target as shown in Figures 4(a) and 4(b). To be specific, the proposed method reaches the CRLB only when ms. Despite the optimality, the RMSE of the target position and the velocity are larger than the one under scenario 1 by around 10 times as illustrated in Figures 4(b) and 4(c). This is because that the far-field target case implies a bad localization geometry, which leads to a significant decrease in target location accuracy. What is more, as shown in Figure 4(b), the position RMSEs of the proposed method-case 1 and proposed method-case 2 almost overlap. This indicates that the contribution of explicitly considering the clock skew in DDS model to the position estimation is marginal for a far-field target. To conclude the effect of the target position on the considered localization methods, a good localization geometry is very significant for all the methods. Although the RMSE performance of our proposed method is heavily deteriorated under a far-field target case, it still outperforms the localization method ignoring the clock imperfection.
(a)
(b)
(c)
Figure 5 plots the RMSE of the considered localization algorithms versus clock skew under scenario 1. The noise standard deviation is set to be 3 ms. In this simulation, we aim to study the impact of the value of clock skew on the considered algorithms. Note that the range of clock skew adopted here is based on the clock synchronization performance of widely used underwater acoustic modem (S2CR series) [37]. As it can be seen from Figure 5, along with the increase of the clock skew (both in positive or negative direction), the localization performance of the WLS method which have no regard for the clock skew gets worse. However, it has almost no effect on our proposed methods. This validates that the proposed localization algorithm is robust to the clock skew variation.
(a)
(b)
(c)
6. Conclusion
A new model for DDS/TDOA-based underwater localization with considering the clock imperfection has been proposed. Based on such a model, two WLS estimators are devised for underwater target parameters (position and velocity) and the clock skew estimation. The first one introduces nuisance variables to eliminate the coupling relationships between parameters; as a result, we obtain a pseudolinear estimation model for calculating a coarse estimate. The second WLS estimator refines the estimate by exploiting the coupling relationships between the target parameter, clock skew, and nuisance variables. The RMSE performances of the proposed estimator are compared with other competitive estimators by computer simulations. The performance of the proposed method is shown in theory and by simulations to reach the CRLB accuracy under sufficiently small noise conditions. As a future direction, we shall validate the performance of the proposed method using real underwater DDS/TDOA measurements.
With regard to the scalability issue, we notice that the proposed algorithm is designed to localize a single target node by using multiple anchor nodes. However, when there presents multiple target nodes, the localization scheme proposed in this work may turn to be inapplicable. In the future work, we aim to develop a multistage underwater node localization scheme to achieve large-scale underwater network positioning. The key idea of multistage scheme lies in that the localized target nodes can be used as anchor nodes to localize other unlocated nodes. As performing iteratively, the localization range is gradually increasing.
Data Availability
All the data used to support the findings of this study are available in this paper.
Conflicts of Interest
The authors declare that there is no conflict of interest.
Acknowledgments
This work was supported in part by the National Nature Science Foundation of China under the grants 52101391, 62101573, and 62101578 and the Project of National University of Defense Technology under grant nos. ZK20-35 and ZK19-36.