# Underwater Localization Using Differential Doppler Scale and TDOA Measurements with Clock Imperfection

Authors: Bingbing Zhang; Jiahua Zhu; Yanqun Wu; Weihua Zhang; Min Zhu
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## 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.


---

## Body

## 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;ai and Ai,j are the ith element of a and the i,jth element of A, respectively; ai:j denotes a subvector with the ith to the jth elements of a; Ai,: represents the ith row of A; · is the Euclidean distance norm; δ· denotes the Dirac delta function; ∗ denotes the convolution operator; superscript T denotes the transpose of a matrix (vector); e denotes the element-wise multiplication; 1 and 0 are vectors (or matrixes) of 1 and 0, O is zero matrix, I denotes the identity matrix (size indicated in the subscript if necessary); diag· and blkdiag· represent the diagonal and the block-diagonal matrices; E· denotes the expectation operator.


## 1.1. Notations

Column vectors and matrices are denoted by bold lower- and uppercase letters, respectively;ai and Ai,j are the ith element of a and the i,jth element of A, respectively; ai:j denotes a subvector with the ith to the jth elements of a; Ai,: represents the ith row of A; · is the Euclidean distance norm; δ· denotes the Dirac delta function; ∗ denotes the convolution operator; superscript T denotes the transpose of a matrix (vector); e denotes the element-wise multiplication; 1 and 0 are vectors (or matrixes) of 1 and 0, O is zero matrix, I denotes the identity matrix (size indicated in the subscript if necessary); diag· and blkdiag· represent the diagonal and the block-diagonal matrices; E· 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 Figure1, consider a three-dimensional underwater localization scenario where N moving anchors are used to determine the position u=x,y,zT and velocity u̇=ẋ,ẏ,żT 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 si=xi,yi,ziT and ṡi=ẋi,ẏi,żiT, i=1,⋯,N, respectively. Assume the anchors are synchronized and behave a common clock skew w and clock offset θ. Thus, the local time of the anchors with respect to (w.r.t.) a universal standard time t is give by [31].
(1)ct=wt+θ.Figure 1
The underwater localization scenario.We assume the local time of the target is the universal standard time, i.e., clock skew is1 and clock offset is 0.To implement a DDS/TDOA information-based localization algorithm, it is necessary to transmit signals from the target to anchors and/orvice 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 theith anchor that has the impulse response [6].
(2)hiτ,t=∑p=1MAi,ptδτ−τi,pt,where M is the number of paths, Ai,pt is the corresponding time-varying path attenuation, and τi,pt 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 ai,p as
(3)τi,pt=τi,p−ai,pt,and the path attenuation are assumed constant Ai,pt=Ai,p. Moreover, we assume all paths have a similar Doppler scale ai, which has already been justified in [6]. Finally, the underwater acoustic channel model is approximated as7
(4)hiτ,t=∑p=1MAi,pδτ−τi,p+ait.Letst be the transmitted passband signal at the target with duration T0, then the received passband signal at the ith anchor is given by
(5)yit=st∗hiτ,t=∑p=1MAi,ps1+ait−τi,p+εit,where εit is additive noise. Sampling yit 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 Ts′=Ts/w, where Ts is the reference sampling interval. Thus, the received signal at the ith anchor is discretized as
(6)yin=yitt=nTs/w=∑p=1MAi,ps1+aiwn−τi,p+εin,where n=wt/TsTs 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 T0 to Ti,r=wT0/1+ai,. The scaled signal duration can be estimated by cross-correlating the received signal with the known LFM preamble and postamble, denoted as T^i,r. Then, by knowing the original signal duration T0, the Doppler scale measurement of the ith anchor is estimated as a^i=T0/T^i,r−1. Thus, the Doppler scale measurement model can be defined as
(7)a^i=1+aiw−1+vi,i=1,⋯,N,where ai is modeled as the actual Doppler scale caused by the target/anchor motion and is given by
(8)ai=u−siTu̇−ṡicrio.c is the acoustic propagation speed, and rio=u−si; vi is the measurement noises, which are independent and identically distributed (i.i.d.) Gaussian random variables with zero-mean and variance σa,i2. 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 ith anchor
(9)τ^i=wt0+rio/c+θ+ni,i=1,⋯,N,where t0 is the unknown start transmission time of the target and ni is TOA estimation error and modeled by i.i.d. Gaussian random variables with zero-mean and variance στ,i2. 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 i=1 as the reference Doppler scale and TOA measurement and form the DDS and TDOA measurements, respectively, as
(10a)Δa^i1=a^i−a^1=1wai−a1+vi1,(10b)Δτ^i1=τ^i−τ^1=wcrio−r1o+ni1,where vi1=vi−v1, ni1=ni−n1, and i=2,⋯,N. It can be seen from (10b) that the parameters t0 and θ are cancelled out by the TDOA calculation while the clock skew w 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 3 paths is used.Figure 2
An example of extracting the Doppler scale and TOA measurements:T0 is the signal duration, w is the clock skew, ai is the Doppler scale, and the underwater acoustic channel has 3 paths with delay τi,1, τi,2, and τi,3, respectively.Stacking the DDS and TDOA measurements and fusing them into a2N−1×1 vector, we obtain
(11)m≜fT,τTT,where
(12)f=1wa2−a1,⋯,aN−a1T+v21,⋯,vN1T=fo+Γv,andΓ=−1N−1,IN−1,v=v1,⋯,vNT,τ=wcr2o−r1o,⋯,rNo−r1oT+n21,⋯,nN1T=τo+Γn,andn=n1,⋯,nNT.The measurement error vector isΔm≜ΓvT,ΓnTT, and its covariance matrix is Qm≜EΔmΔmT=blkdiagΓQvΓT,ΓQnΓT, where Qv=diagσa,12,⋯,σa,N2T and Qn=diagστ,12,⋯,στ,N2T. Note that the Doppler scale measurement noises v are assumed to be independent of the TOA measurement noises n. From (11), we can write the observed DDS and TDOA equation in matrix form as
(13)m=mo+Δm,where mo=foT,τoTT. The goal of underwater localization is estimating u, u̇, and w, from the noisy measurement vector m. Under the mutually independent Gaussian noise condition, the ML estimator of u, u̇, and w can be formulated as [36]
(14)u∧T,u̇∧T,w∧T=minu,u̇,wm−moTQm−1m−mo.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 m is Gaussian distributed as m~Nmo,Qm. The CRLB is the lowest possible variance that an unbiased estimator can achieve. For the unknown vector φ=u,u̇,wT, its CRLB is given by [36]
(15)CRLBφ=Jφ−1,where Jφ is the Fisher information matrix (FIM). Using the notation ∇ba=∂a/∂b, Jφ can be calculated as
(16)Jφ=∇φmoTQm−1∇φmo.The partial derivative in (16) can be expressed as
(17)∇φmo=∇ufo∇u̇fo∇wfo∇uτo∇u̇τo∇wτo,where
(18)∇wfoi−1=−w−2ai−a1,∇ufoi−1:=w−1u̇−ṡiTcrio−aiu−siTrio2−u̇−ṡ1Tcr1o+a1u−s1Tr1o2,∇u̇foi−1:=w−1cu−siTrio−u−s1Tr1o,∇wτoi−1=rio−r1oc,∇uτoi−1:=w2∇u̇foi−1:,∇u̇τo=ON−1×3,i=2,⋯,N.Based on the equations above, the FIM can be easily calculated, and hence, the CRLB is obtained.


## 2.1. Problem Formulation

As depicted in Figure1, consider a three-dimensional underwater localization scenario where N moving anchors are used to determine the position u=x,y,zT and velocity u̇=ẋ,ẏ,żT 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 si=xi,yi,ziT and ṡi=ẋi,ẏi,żiT, i=1,⋯,N, respectively. Assume the anchors are synchronized and behave a common clock skew w and clock offset θ. Thus, the local time of the anchors with respect to (w.r.t.) a universal standard time t is give by [31].
(1)ct=wt+θ.Figure 1
The underwater localization scenario.We assume the local time of the target is the universal standard time, i.e., clock skew is1 and clock offset is 0.To implement a DDS/TDOA information-based localization algorithm, it is necessary to transmit signals from the target to anchors and/orvice 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 theith anchor that has the impulse response [6].
(2)hiτ,t=∑p=1MAi,ptδτ−τi,pt,where M is the number of paths, Ai,pt is the corresponding time-varying path attenuation, and τi,pt 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 ai,p as
(3)τi,pt=τi,p−ai,pt,and the path attenuation are assumed constant Ai,pt=Ai,p. Moreover, we assume all paths have a similar Doppler scale ai, which has already been justified in [6]. Finally, the underwater acoustic channel model is approximated as7
(4)hiτ,t=∑p=1MAi,pδτ−τi,p+ait.Letst be the transmitted passband signal at the target with duration T0, then the received passband signal at the ith anchor is given by
(5)yit=st∗hiτ,t=∑p=1MAi,ps1+ait−τi,p+εit,where εit is additive noise. Sampling yit 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 Ts′=Ts/w, where Ts is the reference sampling interval. Thus, the received signal at the ith anchor is discretized as
(6)yin=yitt=nTs/w=∑p=1MAi,ps1+aiwn−τi,p+εin,where n=wt/TsTs 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 T0 to Ti,r=wT0/1+ai,. The scaled signal duration can be estimated by cross-correlating the received signal with the known LFM preamble and postamble, denoted as T^i,r. Then, by knowing the original signal duration T0, the Doppler scale measurement of the ith anchor is estimated as a^i=T0/T^i,r−1. Thus, the Doppler scale measurement model can be defined as
(7)a^i=1+aiw−1+vi,i=1,⋯,N,where ai is modeled as the actual Doppler scale caused by the target/anchor motion and is given by
(8)ai=u−siTu̇−ṡicrio.c is the acoustic propagation speed, and rio=u−si; vi is the measurement noises, which are independent and identically distributed (i.i.d.) Gaussian random variables with zero-mean and variance σa,i2. 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 ith anchor
(9)τ^i=wt0+rio/c+θ+ni,i=1,⋯,N,where t0 is the unknown start transmission time of the target and ni is TOA estimation error and modeled by i.i.d. Gaussian random variables with zero-mean and variance στ,i2. 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 i=1 as the reference Doppler scale and TOA measurement and form the DDS and TDOA measurements, respectively, as
(10a)Δa^i1=a^i−a^1=1wai−a1+vi1,(10b)Δτ^i1=τ^i−τ^1=wcrio−r1o+ni1,where vi1=vi−v1, ni1=ni−n1, and i=2,⋯,N. It can be seen from (10b) that the parameters t0 and θ are cancelled out by the TDOA calculation while the clock skew w 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 3 paths is used.Figure 2
An example of extracting the Doppler scale and TOA measurements:T0 is the signal duration, w is the clock skew, ai is the Doppler scale, and the underwater acoustic channel has 3 paths with delay τi,1, τi,2, and τi,3, respectively.Stacking the DDS and TDOA measurements and fusing them into a2N−1×1 vector, we obtain
(11)m≜fT,τTT,where
(12)f=1wa2−a1,⋯,aN−a1T+v21,⋯,vN1T=fo+Γv,andΓ=−1N−1,IN−1,v=v1,⋯,vNT,τ=wcr2o−r1o,⋯,rNo−r1oT+n21,⋯,nN1T=τo+Γn,andn=n1,⋯,nNT.The measurement error vector isΔm≜ΓvT,ΓnTT, and its covariance matrix is Qm≜EΔmΔmT=blkdiagΓQvΓT,ΓQnΓT, where Qv=diagσa,12,⋯,σa,N2T and Qn=diagστ,12,⋯,στ,N2T. Note that the Doppler scale measurement noises v are assumed to be independent of the TOA measurement noises n. From (11), we can write the observed DDS and TDOA equation in matrix form as
(13)m=mo+Δm,where mo=foT,τoTT. The goal of underwater localization is estimating u, u̇, and w, from the noisy measurement vector m. Under the mutually independent Gaussian noise condition, the ML estimator of u, u̇, and w can be formulated as [36]
(14)u∧T,u̇∧T,w∧T=minu,u̇,wm−moTQm−1m−mo.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 m is Gaussian distributed as m~Nmo,Qm. The CRLB is the lowest possible variance that an unbiased estimator can achieve. For the unknown vector φ=u,u̇,wT, its CRLB is given by [36]
(15)CRLBφ=Jφ−1,where Jφ is the Fisher information matrix (FIM). Using the notation ∇ba=∂a/∂b, Jφ can be calculated as
(16)Jφ=∇φmoTQm−1∇φmo.The partial derivative in (16) can be expressed as
(17)∇φmo=∇ufo∇u̇fo∇wfo∇uτo∇u̇τo∇wτo,where
(18)∇wfoi−1=−w−2ai−a1,∇ufoi−1:=w−1u̇−ṡiTcrio−aiu−siTrio2−u̇−ṡ1Tcr1o+a1u−s1Tr1o2,∇u̇foi−1:=w−1cu−siTrio−u−s1Tr1o,∇wτoi−1=rio−r1oc,∇uτoi−1:=w2∇u̇foi−1:,∇u̇τo=ON−1×3,i=2,⋯,N.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
(19a)ṙi1=ṙi1o+ε̇i1=1wṙio−ṙ1o+ε̇i1,(19b)ri1=ri1o+εi1=wrio−r1o+εi1,where ṙi1=cΔa^i1, ri1=cΔτ^i1, ε̇i1=cvi1, εi1=cni1, and
(20)ṙio=u−siTu̇−ṡirio.Considering the noise free TDOA measurement in (19b), which can be written as
(21)ri1o=wrio−r1o.Upon rewriting (21) as ri1o+wr1o=wrio, squaring both sides, and substituting rio=u−si for rio2 and r1o2, we have a set of TDOA equations for i=2,⋯,N.
(22)ri1o2+2wri1or1o=w2siTsi−s1Ts1−2si−s1Tu.Taking the time derivative of (22) results in a set of DDS equations for i=2,⋯,N.
(23)wri1oṙi1o+w2ṙi1or1o+ri1oṙ1o=wsiTṡi−s1Tṡ1−ṡi−ṡ1Tu−si−s1Tu̇.In terms of the noisy quantities by puttingri1o=ri1−εi1 and ṙi1o=ṙi1−ε̇i1 into (22) and (23) and ignoring the second-order error terms, we arrive at
(24)e1=h1−G1φ1,where φ1=u,η1,η2,u̇,η3,w1T is defined as the parameter vector, in which η1=r1ow−1, η2=w−2, η3=ṙ1ow−1, and w1=wr1o. The vector h1 and matrix G1 are defined as
(25)h1=s1Ts1−s2Ts2⋮s1Ts1−sNTsN2r21ṙ21+s1Tṡ1−s2Tṡ2⋮2rN1ṙN1+s1Tṡ1−sNTṡN,G1i−1,:=−2si−s1T,ri1,0.5ri12,01×5,G1i+N−2,:=−2ṡi−ṡ1T,01×2,si−s1T,ri1,ṙi1,i=2,⋯,N,and the noise vector e1 is defined as
(26)e1=B1ε1ε2,B1=OBDC,ε1=ε̇21,⋯,ε̇N1T,ε2=ε21,⋯,εN1T,B=2w−1diagr2o,⋯,rNo,C=2w−1diagṙ2o,⋯,ṙNo,D=2wdiagr2o,⋯,rNo.Thanks to the nuisance variablesη1, η2, η3, and w1 introduced in the parameter vector φ1, it makes (24) become a set of linear equations with respect to φ1. As a result, φ1 can be estimated by the WLS method, whose solution is given as [36]
(27)φ^1=G1TW1G1−1G1TW1h1,where W1 is the weighting matrix chosen as
(28)W1=Ee1e1T−1=B1QmB1T−1,where Qm=c2Qm. Note that the weighting matrix is dependent on the true values of clock skew, source position, and velocity through B1. To cope with this, W1=Q~m−1 is first employed to calculate an initial estimate from (27), which will be used back into (28) for an improved version of W1, then leading to a better estimate of φ1. The estimation error in φ^1 can be calculated as
(29)Δφ1=φ^1−φ1=G1TW1G1−1G1TW1e1.The covariance matrix of φ^1 is, therefore, assuming small measurement noise so that the noise in G1 can be ignored [36].
(30)covφ^1≈G1TW1G1−1.


### 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φ1. In fact, they are related to each other through the following equations:
(31a)η1w1=r1o2=u−s12,(31b)η12=η2r1o2=η2u−s12,(31c)η1η3=η2ṙ1or1o=η2u−s1Tu̇−ṡ1,(31d)η3w1=ṙ1or1o=u−s1Tu̇−ṡ1.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φ^1 is Δφ1. Expressing φ^11:3=u+Δφ11:3 and subtracting both sides by s1, we have
(32)φ^11:3−s1⊙φ^11:3−s1≈u−s1⊙u−s1+2u−s1⊙Δφ11:3.where the second-order error terms are ignored. Similarly, expressing φ^16:8=u̇+Δφ16:8 and combing with the position estimates, we have
(33)φ^11:3−s1⊙φ^16:8−ṡ1≈u−s1⊙u̇−ṡ1+u̇−ṡ1⊙Δφ11:3+u−s1⊙Δφ16:8.Substitutingη1=φ^14−Δφ14, η2=φ^15−Δφ15, η3=φ^19−Δφ19, and w1=φ^110−Δφ110 into the equations (31a), (31b), and (31d), we obtain
(34a)φ^110Δφ14+φ^14Δφ110≈φ^14φ^110−u−s12,(34b)2φ^14Δφ14−r1o2Δφ15≈φ^142−φ^15u−s12,(34c)φ^110Δφ19+φ^19Δφ110≈φ^19φ^110−u−s1Tu̇−ṡ1.In order to provide an estimation ofw in the second stage, we incorporates an additional equation
(35)Δφ15=φ^15−w−2.Staking (32)–(35), we arrive at
(36)e2=h2−G2φ2,where
(37)h2=φ^11:3−s1⊙φ^11:3−s1φ^11:3−s1⊙φ^16:8−ṡ1φ^14φ^110φ^142φ^19φ^110φ^15,G2=I3O3×303×1O3×3I303×111×301×30φ^1511×301×3001×311×3001×301×31,φ2=u−s1⨀u−s1u−s1⨀u̇−ṡ1w−2.On the left side of (45), the noise vector e2 is defined as
(38)e2=B2Δφ1,B2=2C1O3×2O3×3D1O3×2C1O4×3E1TO4×3O3×2O3×2E2T,C1=diagu−s1,D1=diagu̇−ṡ1,E1=wr1o2w−1r1o0−r1o20001,E2=00w−1r1o0wr1o0w−1ṙ1o0.Equation (36) is a set of linear equations with respect to φ2; its WLS solution is given by [36]
(39)φ^2=G2TW2G2−1G2TW2h2,where W2 is the weighting matrix, which is given as
(40)W2=Ee2e2T−1=B2covφ∧1B2T−1.To examine the covariance of the second-stage solutionφ^2, subtracting both sides of (39) by φ2 and using (36) gives
(41)Δφ2=φ2−φ^2=G2TW2G2−1G2TW2e2.Hence, under small noise conditions, the covariance of the second stage solutionφ^2 is given by [36]
(42)covφ^2≈G2TW2G2−1.Finally, the clock skew, source position, and velocity estimates can be deduced from the definition ofφ2.
(43a)u^=Πφ^21:3+s1,(43b)u̇^=φ^24:6.u−s1+ṡ1,(43c)w^=1φ^27,where Π=diagsignφ^11:3−s1 is used to avoid the sign ambiguity caused by the square root operation.Similar toW1, the weighting matrix W2 is also dependent on the true values of clock skew, source position, and velocity through B2. In practice, these true values can be substituted by the solution in φ^1 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 Algorithm1.Algorithm 1
The proposed estimatorInput: Anchor’s parameter, DDS/FDOA measurements, the measurement noise covariance.First stage processing:1:Initialization:W1=Q~m−1.2:Forl=1toNiter (Niter is the number of iterations)3: computingφ^1 from (27);4: substitutingφ^1 into (28) to update W1;5:end ForSecond stage processing:6: Computingcovφ^1 using (30).7: Usingφ^1 to calculate B2 and obtaining W2 using (40).8:Forl=1toNiter9: computingφ^2 from (39);10: applying (43a), (43b), and (43c) to generate the estimates;11: substituting the estimates from (43a), (43b), and (43c) in B2 and updating W2 using (40) accordingly;12:end ForOutput: the target position and velocity, and the clock skew estimates.


## 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
(19a)ṙi1=ṙi1o+ε̇i1=1wṙio−ṙ1o+ε̇i1,(19b)ri1=ri1o+εi1=wrio−r1o+εi1,where ṙi1=cΔa^i1, ri1=cΔτ^i1, ε̇i1=cvi1, εi1=cni1, and
(20)ṙio=u−siTu̇−ṡirio.Considering the noise free TDOA measurement in (19b), which can be written as
(21)ri1o=wrio−r1o.Upon rewriting (21) as ri1o+wr1o=wrio, squaring both sides, and substituting rio=u−si for rio2 and r1o2, we have a set of TDOA equations for i=2,⋯,N.
(22)ri1o2+2wri1or1o=w2siTsi−s1Ts1−2si−s1Tu.Taking the time derivative of (22) results in a set of DDS equations for i=2,⋯,N.
(23)wri1oṙi1o+w2ṙi1or1o+ri1oṙ1o=wsiTṡi−s1Tṡ1−ṡi−ṡ1Tu−si−s1Tu̇.In terms of the noisy quantities by puttingri1o=ri1−εi1 and ṙi1o=ṙi1−ε̇i1 into (22) and (23) and ignoring the second-order error terms, we arrive at
(24)e1=h1−G1φ1,where φ1=u,η1,η2,u̇,η3,w1T is defined as the parameter vector, in which η1=r1ow−1, η2=w−2, η3=ṙ1ow−1, and w1=wr1o. The vector h1 and matrix G1 are defined as
(25)h1=s1Ts1−s2Ts2⋮s1Ts1−sNTsN2r21ṙ21+s1Tṡ1−s2Tṡ2⋮2rN1ṙN1+s1Tṡ1−sNTṡN,G1i−1,:=−2si−s1T,ri1,0.5ri12,01×5,G1i+N−2,:=−2ṡi−ṡ1T,01×2,si−s1T,ri1,ṙi1,i=2,⋯,N,and the noise vector e1 is defined as
(26)e1=B1ε1ε2,B1=OBDC,ε1=ε̇21,⋯,ε̇N1T,ε2=ε21,⋯,εN1T,B=2w−1diagr2o,⋯,rNo,C=2w−1diagṙ2o,⋯,ṙNo,D=2wdiagr2o,⋯,rNo.Thanks to the nuisance variablesη1, η2, η3, and w1 introduced in the parameter vector φ1, it makes (24) become a set of linear equations with respect to φ1. As a result, φ1 can be estimated by the WLS method, whose solution is given as [36]
(27)φ^1=G1TW1G1−1G1TW1h1,where W1 is the weighting matrix chosen as
(28)W1=Ee1e1T−1=B1QmB1T−1,where Qm=c2Qm. Note that the weighting matrix is dependent on the true values of clock skew, source position, and velocity through B1. To cope with this, W1=Q~m−1 is first employed to calculate an initial estimate from (27), which will be used back into (28) for an improved version of W1, then leading to a better estimate of φ1. The estimation error in φ^1 can be calculated as
(29)Δφ1=φ^1−φ1=G1TW1G1−1G1TW1e1.The covariance matrix of φ^1 is, therefore, assuming small measurement noise so that the noise in G1 can be ignored [36].
(30)covφ^1≈G1TW1G1−1.


## 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φ1. In fact, they are related to each other through the following equations:
(31a)η1w1=r1o2=u−s12,(31b)η12=η2r1o2=η2u−s12,(31c)η1η3=η2ṙ1or1o=η2u−s1Tu̇−ṡ1,(31d)η3w1=ṙ1or1o=u−s1Tu̇−ṡ1.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φ^1 is Δφ1. Expressing φ^11:3=u+Δφ11:3 and subtracting both sides by s1, we have
(32)φ^11:3−s1⊙φ^11:3−s1≈u−s1⊙u−s1+2u−s1⊙Δφ11:3.where the second-order error terms are ignored. Similarly, expressing φ^16:8=u̇+Δφ16:8 and combing with the position estimates, we have
(33)φ^11:3−s1⊙φ^16:8−ṡ1≈u−s1⊙u̇−ṡ1+u̇−ṡ1⊙Δφ11:3+u−s1⊙Δφ16:8.Substitutingη1=φ^14−Δφ14, η2=φ^15−Δφ15, η3=φ^19−Δφ19, and w1=φ^110−Δφ110 into the equations (31a), (31b), and (31d), we obtain
(34a)φ^110Δφ14+φ^14Δφ110≈φ^14φ^110−u−s12,(34b)2φ^14Δφ14−r1o2Δφ15≈φ^142−φ^15u−s12,(34c)φ^110Δφ19+φ^19Δφ110≈φ^19φ^110−u−s1Tu̇−ṡ1.In order to provide an estimation ofw in the second stage, we incorporates an additional equation
(35)Δφ15=φ^15−w−2.Staking (32)–(35), we arrive at
(36)e2=h2−G2φ2,where
(37)h2=φ^11:3−s1⊙φ^11:3−s1φ^11:3−s1⊙φ^16:8−ṡ1φ^14φ^110φ^142φ^19φ^110φ^15,G2=I3O3×303×1O3×3I303×111×301×30φ^1511×301×3001×311×3001×301×31,φ2=u−s1⨀u−s1u−s1⨀u̇−ṡ1w−2.On the left side of (45), the noise vector e2 is defined as
(38)e2=B2Δφ1,B2=2C1O3×2O3×3D1O3×2C1O4×3E1TO4×3O3×2O3×2E2T,C1=diagu−s1,D1=diagu̇−ṡ1,E1=wr1o2w−1r1o0−r1o20001,E2=00w−1r1o0wr1o0w−1ṙ1o0.Equation (36) is a set of linear equations with respect to φ2; its WLS solution is given by [36]
(39)φ^2=G2TW2G2−1G2TW2h2,where W2 is the weighting matrix, which is given as
(40)W2=Ee2e2T−1=B2covφ∧1B2T−1.To examine the covariance of the second-stage solutionφ^2, subtracting both sides of (39) by φ2 and using (36) gives
(41)Δφ2=φ2−φ^2=G2TW2G2−1G2TW2e2.Hence, under small noise conditions, the covariance of the second stage solutionφ^2 is given by [36]
(42)covφ^2≈G2TW2G2−1.Finally, the clock skew, source position, and velocity estimates can be deduced from the definition ofφ2.
(43a)u^=Πφ^21:3+s1,(43b)u̇^=φ^24:6.u−s1+ṡ1,(43c)w^=1φ^27,where Π=diagsignφ^11:3−s1 is used to avoid the sign ambiguity caused by the square root operation.Similar toW1, the weighting matrix W2 is also dependent on the true values of clock skew, source position, and velocity through B2. In practice, these true values can be substituted by the solution in φ^1 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 Algorithm1.Algorithm 1
The proposed estimatorInput: Anchor’s parameter, DDS/FDOA measurements, the measurement noise covariance.First stage processing:1:Initialization:W1=Q~m−1.2:Forl=1toNiter (Niter is the number of iterations)3: computingφ^1 from (27);4: substitutingφ^1 into (28) to update W1;5:end ForSecond stage processing:6: Computingcovφ^1 using (30).7: Usingφ^1 to calculate B2 and obtaining W2 using (40).8:Forl=1toNiter9: computingφ^2 from (39);10: applying (43a), (43b), and (43c) to generate the estimates;11: substituting the estimates from (43a), (43b), and (43c) in B2 and updating W2 using (40) accordingly;12:end ForOutput: the target position and velocity, and the clock skew estimates.


## 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φ2 defined below (45), we can relate the estimation error of (43a), (43b), and (43c) with that of φ^2 as
(44)Δφ=Δu,Δu̇,ΔwT=B3−1Δφ2,B3=2diagu−s1O3×303×1diagu̇−ṡ1diagu−s103×101×301×3−2w−3.The bias of the final solution is given by taking expectation ofΔφ. Obviously, it can be seen that Δφ is linearly related to ε1,ε2T through the definitions of e1, Δφ1, e2, and Δφ2. Since ε1,ε2T 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
(45)covφ^=B3−1covφ^2B3−T.After substituting the corresponding covariance matrices in (30) and (42), (45) becomes
(46)covφ^=G3TQm−1G3−1,where
(47)G3=c−1B1−1G1B2−1G2B3.Note that (46) has the same form as the CRLB given in (15) and (16). We shall compare covφ^ with the CRLB under the case of far-field target. It has been shown that comparing covφ^ with the CRLB for near-field target is not easy due to the tedious form of covφ^; 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
(48)ir1o≈r2o≈⋯≈rNo,iiṙio/rio≈0,i=1,⋯,N.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 matrixG3. Substituting the relevant matrices into G3, after some straightforward algebraic manipulation and appropriately using the conditions in (48), we can show that the elements of G3 can be approximated as
(49)G3i−1,1:3≈w−1ṡ1−ṡicrio−ṙios1−sicrio2,G3i−1,4:6≈w−1cs1−sirio,G3i−1,7≈−w−1cr1oṙi1rio,G3N+i−2,1:3≈wcs1−sirio,G3N+i−2,4:6≈01×3,G3N+i−2,7≈w−1cr1ori1rio,where i=2,⋯,N. Then, using the variable relationships that ṙio=cai, wṙi1o=ṙio−ṙ1o, and ri1o=wrio−r1o and the approximations that ṙi1≈ṙi1o, ri1≈ri1o and rio≈r1o,i≠1, we arrive at
(50)G3≈∇φmo.This completes the proof that the proposed solution in Section3 can attain the CRLB accuracy for small Gaussian noise and far-field target. Due to the tedious form of covφ^, 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 1/Nexp∑j=1Nexpξ∧j−ξ2, where ξ∧j is the estimate of ξ∈u,u̇,w obtained in the jth trial. Each simulation result is averaged over Nexp=2000 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 1000×1000×1000 cube centered at 500,500,500T, the velocities of the nodes are randomly drawn from −3,3m/s, and the clock skew is randomly drawn from 0.995,1.005. The sound propagation speed is set as 1500 m/s. The errors εi1 and ε̇i1 are zero-mean white Gaussian processes with identical variances of σ2 and 0.01σ2, respectively. In scenario 2, the target is placed in −300,−300,−300T 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.Figure3 plots the RMSE of the considered localization algorithms versus σ∈1,5ms 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 σ≤3ms. 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 8.9×10−4 s, 10×10−4 s, and 7.8×10−4 s, respectively. This indicates the computational attractiveness of the proposed method.Figure 3
RMSE performance versusσ under scenario 1: (a) clock skew; (b) target position; (c) target velocity.
(a)(b)(c)Figure4 plots the RMSE of the considered localization algorithms versus σ∈1,5ms 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 σ≤2ms. 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.Figure 4
RMSE performance versusσ under scenario 2: (a) clock skew; (b) target position; (c) target velocity.
(a)(b)(c)Figure5 plots the RMSE of the considered localization algorithms versus clock skew ∈0.96,1.04 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.Figure 5
RMSE performance versusσ under scenario 2: (a) clock skew; (b) target position; (c) target velocity.
(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.



---

*Source: 6597132-2022-03-21.xml*