Underwater Target Localization and Synchronization for a Distributed SIMO Sonar with an Isogradient SSP and Uncertainties in Receiver Locations

A distributed single-input multiple-output (SIMO) sonar system is composed of a sound source and multiple underwater receivers. It provides an important framework for underwater target localization. However, underwater hostile environments bring more challenges for underwater target localization than terrestrial target localization, such as the difficulties of synchronizing all the underwater receiver clocks, the varying underwater sound speed and the uncertainties of the locations of the underwater receivers. In this paper, we take the sound speed variation, the time synchronization and the uncertainties of the receiver locations into account, and propose the underwater target localization and synchronization (UTLS) algorithm for the distributed SIMO sonar system. In the distributed SIMO sonar system, the receivers are organized in a star topology, where the information fusion is carried out in the central receiver (CR). All the receivers are not synchronized and their positions are known with uncertainties. Moreover, the underwater sound speed is approximately modeled by a depth-dependent sound speed profile (SSP). We evaluate our proposed UTLS algorithm by comparing it with several benchmark algorithms via numerical simulations. The simulation results reveal the superiority of our proposed UTLS algorithm.


Introduction
As a fundamental function for sonar systems, target localization based on multistatic sonar has gained significant attention. Multistatic sonar can provide superior localization accuracy and detection performance as well as improved robustness and flexibility over its monostatic counterpart [1]. There are several typical multistatic systems enabling target localization, such as multi-input multiple-output (MIMO) sonar/radar [2,3], single-input multiple-output (SIMO) sonar/radar [4,5], wireless sensor networks [6] and multistatic radar [7,8]. Multiple trials conducted by the Centre for Maritime Research and Experimentation [9] have shown that multistatic sonar systems can provide better detection, localization and tracking performances than monostatic sonars. In this paper, we focus on the underwater target localization based on the SIMO sonar system which consists of a sound source and multiple receivers organized in a star topology.
Various localization algorithms have been proposed for wireless networks [10][11][12], and terrestrial SIMO or MIMO radar systems [13,14]. However, they cannot be applied to the distributed sonar system directly due to the special challenges posed by underwater environments. Firstly, the underwater sound ray (USR) bends according to the varying underwater sound speed. The curve of the USR gives rise to ranging errors thereby localization errors. Secondly, The mobility of the underwater receivers results in uncertainties of the receiver locations. Therefore, the real locations of the receivers may not be the same as what we have known. Thirdly, time synchronization has to be carried out due to the absence of GPS and clock differences of underwater receivers. Equipping the underwater receivers with accurate atomic clocks [15] is a way to keep all the underwater receiver clocks in synchronization. However, the atomic clock is quite expensive [16]. Various time synchronization algorithms have been proposed for underwater devices [17,18]. However, most of them require a number of signal transmissions which are energy intensive. In addition, it is difficult to replace batteries for underwater receivers due to the hostile underwater environments. Therefore, it is necessary to perform target localization considering the variation of underwater sound speed, the uncertainties of the receiver locations and the time synchronization.
The sonar researchers have proposed a lot of localization algorithms for various sonar systems [19,20]. The authors of [21] perform a fundamental investigation on the performance of elliptic localization. In [22], the authors utilize the time difference of arrival-angle of arrival (TDOA-AOA) measurements obtained by two stations to localize a point source. A Bayesian localization method for multistatic active sonar is proposed in [23]. The authors of [24] take the mobility of the target into consideration and propose a moving target localization algorithm for multistatic sonar. The multi-target localization problem is considered in [25], the authors propose a multi-target positioning algorithm via bistatic range space projection. Efficient closed-form estimators are proposed for multistatic sonar localization in [26].
Some works take the location uncertainties of the anchors into consideration [31][32][33][34][35][36][37][38] in source localization. The authors of [31][32][33] analyse the effect of anchor location errors on source localization. In [34,35], the calibration emitters are utilized to alleviate the effect of anchor location errors on source localization. The authors of [36] consider the problem of locating multiple disjoint sources using TDOAs and frequency differences of arrival (FDMA) in the presence of anchor position and velocity errors, and propose an asymptotically efficient estimator. The synchronization clock bias and sensor position errors are considered in the TDOA source localization in [37]. In [38] a robust localization algorithm using geographic information in bistatic sonar is proposed for source localization. Some researchers attempt to find a close-form solution for underwater source localization [39][40][41][42][43]. Some closed-form solutions for source localization are obtained by utilizing the TDOA measurements [42] or the combination of TDOA and FDMA measurements [40]. The problem of source localization is extended to multiple disjoint sources localization in [41], and an approximately efficient closed-form solution is proposed. In [43], the authors proposed a closed-form solution for joint source and sensor localization. It locates multiple disjoint sources and refines erroneous sensor positions simultaneously.
All the above underwater target localization algorithms assume that the underwater sound speed is constant during the process of localization. The constant sound speed results in the straight-line trajectory of the USR. However, the assumption of the constant underwater sound speed is validated only in some shallow sea environments. Generally, the underwater sound speed varies with depth, temperature, salinity, etc. Especially, it can be modeled by an isogradient depth-dependent sound speed profile (SSP) approximately in deep sea environments [44] as the salinity and temperature almost remain unchanged. Therefore the depth is the only factor that affects the SSP in deep sea. In deep sea environments, the variation of the sound speed makes the trajectory of the USR a curve. Therefore, the constant sound speed assumption makes the performance degradation of underwater localization in deep sea environments.
The authors of [45,46] proposed algorithms to perform localization and tracking of a mobile target with an isogradient SSP. however, the time synchronization is ignored by them. A target localization algorithm is proposed for a distributed SIMO sonar with an isogradient SSP in our previous work [4]. A delay model is presented in an explicit form with respect to (w.r.t.) the locations of receivers and target. The clock parameters of the receivers and the location of the target are estimated with the assumption that the locations of the receivers are accurate. Actually, in order to monitor efficiently, the receivers have to be anchored at some special depths with anchor chains. In this case, the receivers drift with ocean currents. Therefore, the locations of the receivers vary with ocean currents, which means that there are uncertainties in the receiver locations.
In this work, with the consideration of the USR curve, the time synchronization and the uncertainties of the receivers, we propose the underwater target localization and synchronization (UTLS) algorithm for the distributed SIMO sonar system. The receivers consist of one central receiver (CR) and several normal receivers (NRs) organized in a star topology. Assuming that the receivers are anchored at some special depths with anchor chains, all the receivers are not synchronized, and the underwater sound speed is modeled by an isogradient SSP.
We assume that all the receivers are equipped with acoustic modems [47] to communicate with each other. The CR starts the synchronization and target localization process by broadcasting an initial signal to the NRs after the reception of the direct and reflected signals. In the UTLS algorithm, we first employ the expectation maximization (EM) algorithm [48] to estimate the clock skews of the receivers, and apply the weighted least squares (WLS) to estimate the clock offsets sequentially with the clock skews estimated via the EM algorithm. During the estimation of the clock skews, the receiver locations are estimated as latent variables with the consideration of the uncertainties of the receiver locations. We apply the nonlinear weighted least squares (NWLS) to perform target localization with the estimates of the clock parameters and the receiver locations. Simulation results show the superiority of the UTLS algorithm compared to the tailored benchmark algorithms, such as the tailored NWLS algorithms, the algorithms which combine the tailored WLS and NWLS algorithms.
The rest of the paper is organized as follows. In Section 2, we describe the system models which include the sound speed and time delay models, the clock model, the uncertainty model of the receiver locations and the measurement models. In Section 3, the UTLS algorithm is introduced. It includes the EM algorithm, the WLS algorithm and the NWLS algorithm to achieve time synchronization and target localization. We evaluate the performance of the UTLS algorithm by comparing it with the tailored benchmark algorithms in Section 4 through several simulations. Finally, we conclude this paper in Section 5.

System Models
The scenario considered in this research consists of one target whose coordinates are x = [x, y, z] T , one sound source and M + 1 distributed receivers {r i , i = 1, 2, · · · , M, h} as shown in Figure 1. The receivers organized in the fashion of a star topology. The receiver r h is referred as CR, while the receivers r i , i = 1, 2, · · · , M are referred as NRs. The coordinates of the receivers are indicated as {x i , i = 1, 2, · · · , M, h}. We do not know the accurate locations of the receivers as they move with ocean currents. The signals depart from the sound source directly to the receivers r i , i = 1, 2, · · · , M, h are referred as direct signals s d = {s d 1 , · · · , s d M , s d h }; while the signals which depart from the sound source, reflected by the target and received by the receiver r i , i = 1, 2, · · · , M, h are referred as the reflected signals s = {s 1 , · · · , s M , s h }. There are a large amount of research focus on the detection of the direct signals and the reflected signals [49][50][51]. We assume that the detection problem can be solved by the existing approaches, and focus on the target localization and time synchronization algorithm.
In order to estimate the target location x, there are three aspects have to be considered. The first one is the varying underwater sound speed which results in the curve of USR. The second one is time synchronization. The third one is the uncertainties of the locations of the receivers. With the consideration of these three aspects, we employ four types of models: the sound speed and time delay models, the clock model, the uncertainty model of the receiver locations, and the time-of-arrival (TOA) measurement models.

The Sound Speed and Time Delay Models
The first aspect relates to the USR bend. As we know, the underwater sound speed is not constant but varies with temperature, salinity, pressure (depth), etc. The SSP can be approximated as a piece-wise linear function of the depth [52]. In our scenario, the assumption of isogradient SSP is made. Therefore, the SSP c(z) can be expressed as where a is a constant, and b is the sound speed at the surface. Both of them are assumed to be known as a priori. For a single ray such as the ray from the target to the receiver r i as shown in Figure 1, the time-of-flight (ToF) [4] τ i can be expressed in an explicit form w.r.t.
x and x i where 2 denoting the distance between the receiver r i and the target. Similarly, the delays τ hi between the CR and the NRs can be expressed as where d hi is the distance between the CR and the ith NR. The delay τ s between the sound source and the target can be written as where d s is the distance between the sound source and the target. The delays τ d i between the sound source and the receivers can be expressed as (5)

The Clock Model
The second aspect is the time synchronization. In this work, we assumed that all the NRs were fully asynchronous. The clock model can be stated as where µ i and o i are the unknown clock skew and offset of the receiver r i , respectively. The symbol T i stands for the timestamp from the clock of the receiver r i w.r.t. standard time t. For simplicity, we stack all the clock skews and clock offsets into the vectors

The Uncertainty Model of the Receiver Locations
The uncertainties of the receivers are caused by ocean currents and measurement noise. The ocean current pushes the receivers to move within some ranges. It makes the real locations of the receivers unknown. We obtain the receiver locations by measuring the ocean current. However, the measurements are corrupted by noise. Therefore we have to carry out time synchronization and target localization with the corrupted locations of the receivers.

The Ocean Current Model
The double-gyre phenomena in large-scale ocean circulation is typical for the northern mid-latitude ocean basins [53]. It is quite dominant and persistent in oceans and consists of a sub-polar and a sub-tropical gyre. However, the double-gyre model is only suitable for large-scale ocean circulation. It cannot resolve local, small-scale flow disturbances properly. In addition to the double-gyre model, a small-scale gyre model is adapted in [53] to depict the local, small-scale background flow velocity (BFV) disturbances. In this paper, we adopt the ocean current model used in [53] v where v 0 (t) and v 1 (t) are the large-scale and small-scale BFV model, respectively. We refer to [53] for more details of the ocean current model.

The Positions of Receivers
Note that the BFV relates to the position. We denote the BFV at the position of the ith receiver as v i (t). In order to keep expressions concise, the time parameters of the BFV and the positions of the receivers are omitted when it dose not cause ambiguity. There are some assumptions given by

•
In order to monitor the ocean more efficiently, the receivers need deploying at some special depths. Therefore, we assume that all the receivers are anchored at the bottom of the sea with anchor chains. The length of the i-th receiver chain is ζ i .

•
All the receivers are equipped with buoyancy balls to keep them suspending underwater and straightening the anchor chains.

•
All the receivers are equipped with pressure sensors to determine their depths.

•
All the receivers can measure the BFVs around their positions by using an acoustic Doppler current profiler.
• The initial position of the i-th receiver is the receiver's position when the anchor chain is perpendicular. The initial positions are known as prior information.
With these assumptions, the feasible positions of the ith receiver compose a hemispherical surface as shown in Figure 2. Actually, the receivers drift with the BFV, their positions determined by the BFV. As the buoyancy balls always keep the receivers suspending underwater, the depths of the receivers are determined by the strength provided by the BFV. Their horizontal coordinates are determined by the directions of the BFV. The receiver r i first measures its depth, then utilizes the depth measurement to calculate the horizontal distance between its initial position and its current position. The horizontal distance given by where [x i , y i , z i ] T and [x i (0), y i (0), z i (0)] T are the positions of the receiver r i at time t and 0, respectively. The parameter r i is the horizontal distance between the positions of the receiver r i at time t and 0. With the horizontal distance, the feasible positions ith receiver are limited to a circle whose center is [x i (0), y i (0), z i ] and radius is r i . As shown in Figure 3, the receiver drifts with the ambient BFV, the direction of the BFV is the direction of the horizontal coordinates of the receiver. Therefore the horizontal coordinates can be given by [ where v i is the ambient BFV of the receiver r i . In conclusion, by combining the horizontal coordinates with the depth coordinate, we have Actually, it is hard to obtain the real BFV and the real depth of the receiver. We only have the corresponding measurements which are corrupted by measurement noise. Taking the measurement noises into account we havex where n x i = [n ix , n iy , n iz ] T is assumed as a zero-mean Gaussian noise with variance σ 2 Figure 4 reveals the signal time sequence of the target localization process. Firstly, a signal is transmitted by the sound source at time t 0 , and arrives at the target at time t 0 + τ s . We assume that the signal is reflected by the target as soon as it arrives at the target. Secondly, the direct signal s d i received by receiver r i at time t d i = t 0 + τ d i . The receiver r i records the arrival time of the direct signal s d i as T d is . The measurement model of the direct signals can be stated as

The TOA Measurement Models
where n d i ∼ N (0, σ 2 i ), i = 1, 2, · · · , M, h are Gaussian measurement noises of the direct signals s d . Thirdly, the reflected signal s i received by the receiver r i at time t i = t 0 + τ s + τ i . The receiver records the arrival time of the reflected signal as T is . The measurement model of the reflected signal can be represented as where n i ∼ N (0, σ 2 i ) is the Gaussian measurement noise of the reflected signal s i with σ 2 i known as a priori. The measurement noises of each receivers are assumed to be independent with each other. Therefore, the elements of the set {n i , n d i , i = 1, 2, · · · , M, h} are independent with each other.
The source The target The central receiver propagation delay between the CR and the NR r i . The NR r i records the signal arrival time and denotes it as T ih . The measurement models of the signals from the CR to the NRs are expressed as where n ih ∼ N (0, σ 2 i ), i = 1, 2, · · · , M are independent Gaussian measurement noises. Fifthly, the NR r i replies a signal to the CR at time t ib after receiving the broadcast signal from the CR. The signal contains all the recorded timestamps T d is , T is , T ih and T ib of itself. The CR records the arrival time of the signal as T hi . The measurement models of the signals from the NRs to CR are represented as where n hi ∼ N (0, σ 2 h ), i = 1, 2, · · · , M are independent and identically distributed (IID) Gaussian measurement noises.

The UTLS Algorithm
We considered the problem of underwater target localization and synchronization for the distributed SIMO sonar system. Firstly, we combined the direct signals and the broadcast signals by the receivers to carry out time synchronization. Sequentially, we localized the target with the estimated clock parameters and the TOAs of reflected signals.

Time Synchronization
In this section, we perform time synchronization by estimating the clock skews of the receivers firstly, and estimating the clock offsets of the receivers sequentially.

Estimation of the Clock Skews
Different from the situation of our previous work [4], there are uncertainties in the locations of the receivers. In order to reduce the influence of the uncertainties on the target localization, we estimate the clock skews of the receivers in conjunction with the locations of the receivers with EM algorithm. By combining (13)-(16), we have which give that whereε i = µ h µ i n ih + n hi denoting the noise. It follows from (11) and (17b) that which can be combined as whereε i = µ h n d i µ i + n hi − n d h denoting the noise. We combine (18) and (20) and represent it in the vector form as where Θ = [ϑ T , µ h ] T denoting the parameter vector, the vector ϑ = [µ h /µ 1 , · · · , µ h /µ M ] T , and the vector ε = [ε 1 , · · · ,ε M ,ε 1 , · · · ,ε M ] T denoting the noises. The measurements related to the target position are not included in (21). Therefore, (21) is independent on the target position. The covariance matrix of ε is of the form where The vector ϕ ∈ R 2M is the observation which can be stated as The entries of h are related to the locations of the receivers which are not exactly known. In order to obtain a fixed observation matrix, we linearize the entries of h by the Taylor expansion as The vector h can be approximated as 2τ hM ] T denoting the Taylor expansion point of h. The vector γ is named as latent vector whose entries are latent variables. It can be expressed as where collecting the partial derivatives of h w.r.t. the locations of the receivers. Its submatrices can be expressed as By substituting (23) into (21), we can modify (21) to whereH = [Ah], the set {ϕ, γ} is termed complete data set. The parameter Θ can be estimated if the complete data set is known. The EM algorithm can be utilized to estimate the latent vector γ firstly and the parameter vector Θ sequentially. We employ the EM algorithm to estimate the parameter vector Θ by maximizing the expectation of the complete data logarithm likelihood function under the posterior distribution of the latent variables. The process consists of two steps named as the expectation step (E-step) and the maximization step (M-step). The purpose of the E-step is to estimate the latent variables by calculating the conditional expectation. The M-step aims to estimate the parameters with the estimated latent variables. They are explained as follows.
Since the latent vector γ and the parameter vector Θ are independent, the second term in the right hand side of (27) can be omitted as it is independent on Θ. Therefore the function Q(Θ, Θ (l−1) ) can be represented as Through the derivation (see Appendix A), we rewrite the objective function as where the variablesγ andV s can be calculated according to (A2) and (A4), respectively.

M-Step:
The maximization step attempts to find an optimal Θ which maximizes the objective function Q(Θ, Θ (l−1) ). We take the first order derivative of the object function Q(Θ, Θ l−1 ) w.r.t. Θ and set the result to zero, which gives Define and where The Equations (30) and (31) can be combined into the following vector form which gives the estimate of Θ as According to (23), errors have been introduced by the linearization of h. In order to reduce the errors, we need to updateh by updating the Taylor expansion pointx r . Assuming thatx r and the real locations x r of the receivers. Actually, during the l-th iteration, the latent vectorγ (l) collects the location offsets of the receivers, whilex (l) r is the estimate of the receivers' locations. Therefore, we can estimate the locations of the receivers after the last iteration byx r =x (l) r +γ (l) . DenoteΘ as the estimate of Θ, the clock skews of the receivers can be calculated as followsμ The details of the EM based clock skews estimation algorithm are stated in Algorithm 1.

Algorithm 1 : Expectation maximization (EM) based clock skews estimation algorithm.
Require: T ib , T hi , T ih , i = 1, 2, · · · , M, T hb , a, b, x h , L, Nx i , T d is , i = 1, 2, · · · , M, h, Ensure:μ i , i = 1, 2, · · · , M, h,x r 1: l = 1, index = 1,x linearize h atx r according to (23) 5: E-Step: 6: calculateγ (l) according to (A2) 7: calculateV (l) s according to (A4) 8: M-Step: 9: calculate Θ (l) according to (35) r +γ (l) , l = l + 1 14: end if 15: end while 16 We remark that the uncertainties of the receiver locations are considered in two steps. In the first step, the statistical properties of the receiver locations are utilized to estimate the location offsets and obtain the estimateγ. In the second step, the estimateγ is used to update the receiver locations, which can reduces the influence of the uncertainties of the receiver locations.

Estimation of the Clock Offsets
We employ the WLS algorithm to estimate the clock offsets of the receivers. By substituting the estimates of the clock skews into (11), (17a) and (17b), we have where the vector ϕ o ∈ R 3M+1 can be stated as = [n d 1 , · · · , n d M , n d h , n 1h , · · · , n Mh , n h1 , · · · , n hM ] T ∈ R 3M+1 denoting the noise. Its covariance matrix can be expressed as . By applying the WLS algorithm, we can express the estimates of the clock offsets aŝ where the weighted matrix is defined as W o = Σ −1 .

Target Localization
By utilizing the estimates of the clock parameters, (12) can be rewritten as T and ω = [n 1 /μ i , · · · , n M /μ M , n h /μ h ] T as the observation vector, the propagation delay vector and the noise vector, respectively. We rewrite (39) as The covariance matrix of the noise ω is . It is hard to give an analytical solution of the estimation of the target's location as the propagation delays f s (x, x s ) and f i (x, x i ) are nonlinear w.r.t. the variable x. We employ the NWLS algorithm for the target localization. The iterative process can be expressed as where the weight matrix W = Σ −1 ω , and the matrix For the l-th iteration, the Jacobi matrix D (l−1) can be stated as

Numerical Simulations
We evaluate the performance of our proposed algorithms via the software "Matlab" with an isogradient SSP and different variances of the locations of the receivers. At first, we introduce the simulation setup. Then we show the estimation performance including the parameters of the receivers and the location of the target.

Simulation Setup
As the clock skews of the receivers are designed to be 1, we assume that the real value of the clock skews follow the Gaussian distributions  Table 1. As shown in Table 1, the competitors include the NWLS-UTL algorithm, the WLS-NWLS algorithm, the AWLS-NWLS algorithm, the ANWLS algorithm, and the NWLS-RP algorithm. The NWLS-UTL algorithm is proposed in [4] for underwater target localization. In our proposed UTLS algorithm, the NWLS algorithm is employed to estimate the location of the target with the values of the clock skews and clock offsets estimated by the EM and WLS algorithms, respectively. For comparison, the WLS-NWLS algorithm employs the NWLS algorithm to localize the target with the estimates of the clock skews and clock offsets obtained by the WLS algorithms. The difference between the AWLS-NWLS and WLS-NWLS is that the values of the clock skews, which are utilized in the estimations of the clock offsets and target location, are estimated by an AWLS algorithm in the AWLS-NWLS algorithm. In the ANWLS algorithm, the clock skews and clock offsets of the receivers are assumed to be "1" and "0", respectively. The NWLS-RP algorithm acts as a benchmark since the real values of the clock parameters and the locations of the receivers are used. While the measurements of the receiver locations are used in the NWLS-UTL, the WLS-NWLS, the AWLS-NWLS and ANWLS algorithms.
Firstly, we illustrate the necessity of the consideration of the uncertainties of the receiver locations. We show the performance of the estimation of the clock skews by the comparison of the proposed EM algorithm with the WLS algorithm of the WLS-NWLS algorithm, the AWLS algorithm of the AWLS-NWLS algorithm (whose details are shown in Appendix B) and the tailored WLS of the NWLS-UTL algorithm. We show the performance of the estimations of the clock offsets by comparing the WLS algorithms of the UTLS, the WLS-NWLS and the AWLS-NWLS algorithms. Sequently, we show the necessarity of taking time synchronization into account. We reveal the target localization performance of our proposed UTLS algorithm by comparing the NWLS algorithm of the UTLS with the NWLS algorithms of the WLS-NWLS, the AWLS-NWLS, the ANWLS and the NWLS-RP algorithms. In addition, we compare the UTLS algorithm with the NWLS-UTL algorithm which is proposed in our previous work [4] to show the improvement.
For ease of presentation, we replace the root mean square errors (RMSEs) of the estimations of the clock parameters and the locations of the NRs with the corresponding average value which can be stated as where the variable i can be chosen as µ i , o i or x i , the variable N denotes the number of the variable i . For example, N = M if the variable i stands for µ i , i = 1, 2, · · · , M. For simplicity, we utilize "the RMSE of µ i \o i \x i " to stand for "the average RMSE of µ i \o i \x i , i = 1, 2, · · · , N", from now on.

Performance of the Estimations of the Clock Skews and Locations of the Receivers
In this subsection, we show the performance of the clock skew estimations carried out by the algorithms listed in the second column of Table 1. The estimation performance of the receiver locations via the EM algorithm of the UTLS algorithm is shown in this subsection too. The simulation results of the estimation of the clock skews and the receiver locations are shown in Figures 5-7. For the estimation performance of the clock skews and the receiver locations, we have the following remarks 1. The UTLS algorithm is the best one, the WLS-NWLS algorithm is the second one and the AWLS-NWLS is the worst one. 2. The performance of the estimation of the clock skews can be improved by updating the locations of the receivers. The main difference between the EM algorithm of the UTLS algorithm and the WLS algorithm of the WLS-NWLS algorithm is that the EM algorithm estimates the locations of the receivers by calculating the conditional expectation and utilizes the estimates of the locations for the estimation of the clock skews, while the WLS algorithm of the WLS-NWLS algorithm only utilizes the statistical properties for the estimation of the clock skews. As shown in Figure 5, the UTLS algorithm is better than the WLS-NWLS algorithm and the difference increases with the decreasing of the variance σ 2 i . This is because that with a smaller variance σ 2 i , the UTLS algorithm obtains more accurate estimates of the locations of the receivers as shown in the upper subplot of Figure 7, with which the estimation accuracy of the clock skews can be improved. 3. By utilizing the statistical properties of the locations of the receivers, the estimation performance of clock skews can be improved as shown in Figure 5. The main difference between the WLS algorithm of the WLS-NWLS algorithm and the AWLS algorithm of the AWLS-NWLS algorithm is that the WLS algorithm utilizes the statistical properties of the receivers' locations, while the AWLS dose not. The performance superiority of the WLS-NWLS algorithm compared with the AWLS-NWLS algorithm arises from the utilizing of the statistical properties. 4. In the NWLS-UTL algorithm, the clock offsets are eliminated by a designed subtraction of the time delay. The clock skews are estimated by a tailored WLS algorithm. Its. estimation performance of the clock skews is similar to the WLS-NWLS and the AWLS-NWLS algorithms and poorer than the UTLS algorithm. 5. As shown in Figure 6, the larger the variance σ 2 x of the receiver location noises, the greater the superiority of the UTLS algorithm. 6. The UTLS algorithm is more robust against the location uncertainties of the receivers than the WLS-NWLS, the AWLS-NWLS and the NWLS-UTL algorithms. As shown in Figure 6 the performance of the UTLS algorithm decreases slowly with the increasing of the variance σ 2 x . While the WLS-NWLS, the AWLS-NWLS and the NWLS-UTL algorithms decrease faster than the UTLS algorithm with the increasing of the variance σ 2 x . The main reason is that the accuracy of the estimates of the receivers' locations via the EM algorithm of the UTLS algorithm declines slowly with the increasing of the variance σ 2 x as shown in the lower subplot of Figure 7.

Performance of the Estimations of the Clock offsets of the Receivers
In this subsection, we show the performance of the clock offset estimations carried out by the algorithms listed in the third column of Table 1. We show the simulation results of the estimation of the clock offsets in Figures 8 and 9. For the estimate of the clock offsets, we have the following remarks 1. For the estimation of the offsets of the receivers, the UTLS algorithm is the best and the differences between the UTLS algorithm and the other two algorithms increase with the decreasing of the variance σ 2 i as shown in Figure 8. There are two reasons. The first one is that the differences of the clock skew estimates between the EM algorithm of the UTLS algorithm and the other two algorithms increase with the decreasing of the variance σ 2 i . The second one is that the EM algorithm of the UTLS algorithm provides more accurate locations of the receivers for the WLS algorithm of the UTLS algorithm, and the estimation accuracy of the locations of the receivers increases with the decreasing of the variance σ 2 i . 2. The performance superiority of the WLS algorithm of the WLS-NWLS algorithm compared with the WLS algorithm of the AWLS-NWLS algorithm arises from the utilizing of the statistical properties of the receivers' locations as shown in Figure 8. 3. As shown in Figure 9, the WLS algorithm of the UTLS algorithm is more robust against the uncertainty of the receivers' locations. This characteristic is inherited from the EM algorithm of the UTLS algorithm.

Performance of the Estimation of the Location of the Target
In this subsection, we show the performance of the target localization carried out by the algorithms listed in the last column of Table 1. The simulation results of the target localization are shown in Figures 10 and 11. We have the following remarks on the simulation results of the target localization.
1. The performance of the UTLS is the best compared with that of the WLS-NWLS, the AWLS-NWLS and the ANWLS algorithms as shown in Figure 10. The advantages of the UTLS arise from two reasons. First, the estimation accuracies of the clock skews and offsets carried out by the EM and WLS algorithms of the UTLS algorithm, respectively, are higher than that of the WLS-NWLS, the AWLS-NWLS and the ANWLS algorithms. Second, the estimation accuracy of the receivers' locations by the EM algorithm of the UTLS algorithm is higher than the corresponding measurements of the receivers' locations. 2. The performance of the WLS-NWLS algorithm is better than the AWLS-NWLS and the ANWLS algorithms. This superiority of the WLS-NWLS arises from the WLS algorithms of the WLS-NWLS algorithm. 3. As shown in Figure 10, the performance variation trend of the UTLS is almost the same as the NWLS-RP in which the real clock parameters and locations of the receivers are used. While the performance of the other algorithms almost do not vary with the variance of σ 2 i when the variance σ 2 i is small. 4. The performance of the UTLS algorithm is almost the same as the WLS-NWLS algorithm when the variance σ 2 x is small as shown in Figure 11. The reasons can be stated as follows. As shown in the lower subplot of Figure 7, the difference between the estimated locations of the receivers and the measurements of the locations is small when the variance σ 2 x is small. As shown in Figures 6 and 9, this phenomenon makes the performance superiority of the EM and WLS algorithms of our proposed UTLS algorithm small when compared to the WLS algorithms of the WLS-NWLS algorithm. 5. The UTLS algorithm is still better than the AWLS-NWLS and ANWLS algorithms when the variance σ 2 x is small. The larger the variance σ 2 x , the greater the superiority of the UTLS algorithm. 6. The localization performance of the NWLS-UTL algorithm is poor as shown in Figure 10.
This phenomenon may arise from the fact that the NWLS-UTL algorithm is susceptible to the location uncertainties of the receivers as shown Figure 11. The clock offsets have been eliminated by the NWLS-UTL algorithm. Therefore, the computational cost is smaller because of the absence of the clock offset estimates.
7. The localization performance of the NWLS-RP algorithm is better the our proposed UTLS algorithm. The reason is that the real clock parameters and receiver locations are employed in the NWLS-RP algorithm, which is impractical. In summary, the superiority of the EM algorithm of our proposed UTLS algorithm arises from the consideration of the uncertainties of the locations of the receivers and the estimation of the receivers' locations. With the estimates of the EM algorithm, the WLS algorithm of our proposed UTLS algorithm obtains more accurate clock offsets. The NWLS algorithm of our proposed UTLS algorithm achieves more accurate target localization with the estimates of the clock skews and offsets obtained via the EM and WLS algorithms of our proposed UTLS algorithm, respectively. The superiority of the UTLS algorithm increase with the increasing of the variance σ 2 x .

Conclusions
In this work, we propose the target localization and time synchronization (UTLS) algorithm for the distributed SIMO sonar system with an isogradient SSP and the uncertainties of the locations of the receivers. The UTLS algorithm contains two steps. In the first step, we carry out time synchronization via the EM and the WLS algorithms of our proposed UTLS algorithm. The locations of the receivers are also estimated by the EM algorithm in this step. In the second step, the NWLS algorithm of our proposed UTLS algorithm is employed to perform target localization with the estimates of the clock parameters and the locations of the receivers. We compared the UTLS algorithm with the NWLS-UTL, the WLS-NWLS, the AWLS-NWLS, the ANWLS and the NWLS-RP algorithms. The simulation results show that our proposed UTLS algorithm is better than the NWLS-UTL, the WLS-NWLS, the AWLS-NWLS and the ANWLS algorithms, and the advantage of the UTLS increases with the increasing of the variance σ 2 x . An isogradient SSP is an approximation of the real underwater SSP. In order to obtain a more accurate measurement of distance, we will take a more elaborate SSP into consideration for target localization in our future work.