A Co-Localization Algorithm for Underwater Moving Targets with an Unknown Constant Signal Propagation Speed and Platform Errors

Abstract

Underwater mobile acoustic source target localization encounters several challenges, including the unknown propagation speed of the source signal, uncertainty in the observation platform’s position and velocity (i.e., platform systematic errors), and economic costs. This paper proposes a new two-step closed-form localization algorithm that jointly estimates the angle of arrival (AOA), time difference of arrival (TDOA), and frequency difference of arrival (FDOA) to address these challenges. The algorithm initially introduces auxiliary variables to construct pseudo-linear equations to obtain the initial solution. It then exploits the relationship between the unknown and auxiliary variables to derive the exact solution comprising solely the unknown variables. Both theoretical analyses and simulation experiments demonstrate that the proposed method accurately estimates the position, velocity, and speed of the sound source even with an unknown sound speed and platform systematic errors. It achieves asymptotic optimality within a reasonable error range to approach the Cramér–Rao lower bound (CRLB). Furthermore, the algorithm exhibits low complexity, reduces the number of required localization platforms, and decreases the economic costs. Additionally, the simulation experiments validate the effectiveness of the proposed localization method across various scenarios, outperforming other comparative algorithms.

1. Introduction

Underwater mobile acoustic source target localization technology has been a significant research challenge in marine scientific research, resource exploration, and military applications [1,2,3]. It is well known that a GNSS provides 4D positioning (XYZT) on land and water and in the air and space. GNSSs are less suitable for underwater measurements XYZT due to the known limitations of the propagation of electromagnetic waves underwater. As a rule, underwater vehicles can determine their position from XYT acoustic signals from multiple surface sources. The receiver “listens” to radio acoustic repeaters, receives their messages, and determines the geographic coordinates of the acoustic speakers [4,5]. In [6], GIBs (GPS intelligent buoys) are a group of buoys that calculate their own GPS position and then send acoustic signals to a submersible as the basis for calculating its position. Furthermore, [7] uses buoys that act as mobile pseudolites and convert the GPS radio service into a service based on underwater acoustics. Using this service, an unlimited number of users can independently determine locations and navigate underwater. Briefly, a range of acoustic detection platforms (hereafter referred to as platforms), like underwater vehicles (UUVs) [8,9,10], underwater sensor networks (USNs) [11,12,13], and surface buoys [6,7], can use their position, velocity, and other information to provide the necessary assistance in the localization of underwater acoustic sources. These multiple underwater platforms enable the acquisition and processing of underwater acoustic source signals through acoustic sensors. The parameter information extracted from these signals, including the time, frequency, spatial, and energy domains, facilitates target localization. This information includes time of arrival (TOA) [14], time difference of arrival (TDOA) [10,12,15], frequency difference of arrival (FDOA) [16], angle of arrival (AOA) [11,17], signal strength (RSS) [18], and various combinations thereof [2,3,19,20,21,22]. However, the TOA method necessitates precise clock synchronization, which is both demanding and expensive in terms of the hardware equipment. Signal strength is influenced by water temperature, salinity, and flow velocity, resulting in low positioning accuracy. Combining the AOA and TDOA can only be used to estimate the position of the target [19,21], so currently, joint TDOA/FDOA positioning is often used for mobile target localization [8,22,23], which can achieve the estimation of the position and velocity to the target. However, this method demands the participation of at least five or more platforms in the localization, significantly elevating its economic costs and deployment complexities. Integrating bearing information can enhance the localization performance and reduce the number of required platforms. Thus, in this study, we will develop algorithms jointly utilizing AOA, TDOA, and FDOA observations.

The fusion of the above observations is typically nonlinear, rendering target localization a challenging problem. Consequently, numerous localization algorithms have been proposed to address this issue. The typical algorithms encompass closed-form and non-closed-form solution approaches. The closed-form solution algorithms primarily consist of the two-step weighted least squares (TS-WLS) algorithm [1,21,24,25], weighted spherical interpolation (WSI) [26], and similar techniques. The principal non-closed solution methods include the Taylor series iterative method and the semidefinite programming algorithm (SDP) [2,3,11,22]. Closed-form solutions exhibit low computational complexity. However, in the presence of significant measurement noise, the algorithm’s localization performance may suffer. Conversely, the Taylor series iterative method demands precise initial value selection and multiple iterations. On the other hand, the SDP algorithm, while demonstrating superior performance in high-error scenarios, imposes strict relaxation requirements and entails greater computational complexity. In this paper, we opt for the TS-WLS algorithm due to its closed-form solution and low computational complexity, and improvements have been made on that basis.

The aforementioned algorithms can attain the Cramér–Rao lower bound (CRLB) under reasonable observation error conditions [27]. However, due to the complexity of the marine environment, these localization algorithms are not directly applicable to underwater multi-platform localization scenarios. Figure 1 illustrates a real underwater multi-platform localization model. Firstly, sensors typically measure the platform’s underwater position, velocity, and other parameters, which suffer from inaccuracies and are inevitably influenced by measurement errors. We refer to these inaccuracies as the platform’s systematic errors. The neglect of such errors significantly degrades the localization performance [28]. Secondly, the propagation speed of the acoustic source signal is influenced by water temperature, water pressure, and seawater salinity [29,30], introducing uncertainty. In certain shallow water areas [31], it may even be entirely unpredictable. Therefore, in this study, we treat the sound source signal velocity as an unknown constant.

Some algorithms have been proposed to jointly estimate the propagation velocity of the sound source signal and the target source location in response to these issues. However, most of these existing algorithms are inadequate. Rui and Ho [32] proposed a four-step closed-form solution algorithm for localization, capable of estimating the sound velocity and target location. However, it is limited to two-dimensional scenarios and involves three pseudo-linear transformations to derive the closed-form solution, making the process cumbersome. Yang [31] employed the semidefinite relaxation technique (SDR) to convert a nonconvex localization problem into a semidefinite programming problem (SDP) to estimate the target position and sound velocity. Although it enhances the target localization accuracy in the presence of substantial sensor measurement errors, it demands significant computational resources. Fan [33] and Jia [34] enhanced the four-step method proposed by Rui and Ho [32] using the Lagrange multiplier method and the generalized trust region subproblem (GTRS), respectively, simplifying it into a two-step process to reduce the algorithmic complexity. However, these methods necessitate clock synchronization constraints. Sun [35] considers clock synchronization errors and employs a two-step weighting algorithm to estimate the target position and velocity, thereby enhancing the localization accuracy. However, these methods are tailored to multi-base static sonar scenarios and are solely capable of estimating the target position. Zhang [36] introduces a joint localization closure solution algorithm based on TDOA-FDOA, capable of estimating the position, velocity, and speed of sound of the moving target. However, it necessitates iterative solving, thereby increasing the algorithmic complexity. Additionally, the algorithm stipulates that the number of cooperative platforms involved in localization should not be fewer than five, potentially restricting its applicability in practical underwater localization scenarios. In conclusion, the current literature lacks a joint AOA/TDOA/FDOA method to address the unknown propagation speed of sound source signals and the uncertainty in the position and velocity information on moving platforms in underwater moving sound source target research.

This paper addresses the localization of underwater moving sound source targets using multiple mobile platforms with joint AOA/TDOA/FDOA and presents a two-step closed-form localization algorithm. Additionally, we derive and analyze the corresponding CRLB, conduct algorithm performance and complexity analyses, and verify through theoretical and simulation studies that the proposed method achieves the CRLB within an acceptable error margin. The primary contributions and innovations of this paper include:

• We propose a novel joint AOA/TDOA/FDOA target localization algorithm for underwater moving sound sources. This algorithm effectively estimates the target sound source position, velocity, and sound source signal velocity even under an unknown sound velocity and platform systematic error conditions. Additionally, it reduces the number of localization platforms, thus lowering the economic costs.
• We develop a new two-step closed-form solution that introduces fewer auxiliary variables, resulting in an exact solution containing only unknown variables. This approach avoids complex iterative operations and reduces the computational complexity of the algorithm.
• The proposed algorithm demonstrates asymptotic optimality within an acceptable error range. In five distinct localization scenarios, the root mean square error of the estimated parameters reaches the CRLB.
• The proposed method achieves effective localization even in scenarios involving fewer moving platforms, uncertain platform–target geometries, and far-field acoustic source targets. It also surpasses the other compared algorithms in terms of estimation performance.

The remainder of the paper is structured as follows: Section 2 provides a detailed description of the underwater localization model utilized in this paper. Section 3 derives the CRLB for this localization model and investigates the impact of an unknown sound speed on optimal estimation through simulation experiments. Section 4 outlines the localization method adopted in this paper. Section 5 thoroughly analyzes the performance and computational complexity of the algorithms proposed in this paper. The numerical simulation experiments in the five scenarios presented in Section 6 validate the correctness of the preceding theory and demonstrate the superiority of the algorithm proposed in this paper over other algorithms. Section 7 provides a summary of the entire paper.

The following notation will be used throughout this paper. Bold lowercase letters such as $\mathit{s}$ denote vectors, bold uppercase letters such as $\mathit{S}$ denote matrices, ${\mathit{S}}^T$ denotes the transpose matrix of $\mathit{S}$, ${(•)}^{\mathit{o}}$ denotes the actual value of $(•)$, $\parallel \mathit{a}\parallel$ denotes the 2-norm of the vector $\mathit{a}$, ${\mathit{O}}_{m\times n}$ denotes the zero matrix of $m\times n$, ${\mathit{I}}_{m\times m}$ denotes the unit matrix of $m\times m$, $\otimes$ is the Kronecker product, $\mathrm{diag}(\mathrm{a},\mathrm{b})$ is a diagonal matrix with the diagonal elements $\mathrm{a},\mathrm{b}$, $\mathrm{blakdiag}\left[R_1,R_2\right]$ is a matrix with $R_1$ and $R_2$ on the diagonal and 0 as the other elements, $\mathit{A}(n,:)$ denotes the $n_{th}$ row vector of the matrix $\mathit{A}$, $\mathit{G}\left(n,j:k\right)$ is a vector consisting of the $n_{th}$ row and the $j_{th}$ to $k_{th}$ column elements of the matrix $\mathit{G}$, and $\mathit{\Psi }(j:k)$ is a vector consisting of the $j_{th}$ to $k_{th}$ elements of the vector $\mathit{\Psi }$.

2. Localization Model

This section describes the utilization of a multi-system comprising $\mathit{N}$ underwater mobile platforms to localize a single moving sound source target. As depicted in Figure 1, first, we define the true position of the moving target source: ${\mathit{u}}^o={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$, ${\dot{\mathit{u}}}^o={\left[\begin{array}{ccc}\dot{x}& \dot{y}& \dot{z}\end{array}\right]}^T$. The true position and true velocity of the moving platform are, respectively, ${\mathit{s}}^{\mathit{o}}={\left[{({\mathit{s}}_1^{\mathit{o}})}^T,{({\mathit{s}}_2^{\mathit{o}})}^T,\dots ,{({\mathit{s}}_{\mathit{N}}^{\mathit{o}})}^T\right]}^T$, ${\dot{\mathit{s}}}^{\mathit{o}}={\left[{({\dot{\mathit{s}}}_1^{\mathit{o}})}^T,{({\dot{\mathit{s}}}_2^{\mathit{o}})}^T,\dots ,{({\dot{\mathit{s}}}_{\mathit{N}}^{\mathit{o}})}^T\right]}^T$, where ${\mathit{s}}_{\mathit{n}}^{\mathit{o}}={\left[x_n^o y_n^o{ \mathit{z}}_n^o\right]}^T$, ${\dot{\mathit{s}}}_{\mathit{n}}^{\mathit{o}}={\left[{\dot{x}}_n^o {\dot{y}}_n^o{\dot{\mathit{z}}}_n^o\right]}^T$, $n=1,2,\dots ,N$. In practice, the position and speed of underwater mobile platforms are also measured using sensors, so they will be affected by measurement errors, which are hereinafter collectively referred to as mobile platform systematic errors. Denoting ${\mathit{s}}_{\mathit{n}}={\mathit{s}}_{\mathit{n}}^{\mathit{o}}+\Delta {\mathit{s}}_{\mathit{n}}$ and ${\dot{\mathit{s}}}_{\mathit{n}}={\dot{\mathit{s}}}_{\mathit{n}}^{\mathit{o}}+\Delta {\dot{\mathit{s}}}_{\mathit{n}}$, $\Delta {\mathit{s}}_{\mathit{n}}$ and $\Delta {\dot{\mathit{s}}}_{\mathit{n}}$ are random measurement errors. The position and velocity parameters of each moving platform can be represented in vector form.

$$\mathit{S}={\mathit{S}}^o+\Delta \mathit{S}.$$

(1)

In (1), $\mathit{S}={[{\mathit{s}}^T,{\dot{\mathit{s}}}^T]}^T={\left[{({\mathit{s}}_1)}^T,{({\mathit{s}}_2)}^T,\dots ,{({\mathit{s}}_N)}^T,{({\dot{\mathit{s}}}_1)}^T,{({\dot{\mathit{s}}}_2)}^T,\dots ,{({\dot{\mathit{s}}}_N)}^T\right]}^T$ and ${\mathit{S}}^o={[{\mathit{s}}^{oT},{\dot{\mathit{s}}}^{oT}]}^T$ are used to denote the vector of the measured parameters of the moving platforms and the vector of the actual parameters, $\Delta \mathit{S}={[\Delta {\mathit{s}}^T,\Delta {\dot{\mathit{s}}}^T]}^T={\left[{(\Delta {\mathit{s}}_1)}^T,{(\Delta {\mathit{s}}_2)}^T,\dots ,{(\Delta {\mathit{s}}_N)}^T,{(\Delta {\dot{\mathit{s}}}_1)}^T,{(\Delta {\dot{\mathit{s}}}_2)}^T,\dots ,{(\Delta {\dot{\mathit{s}}}_N)}^T\right]}^T$, $\Delta \mathit{S}$ follows a Gaussian distribution with a mean of zero. Its covariance matrix is $E(\Delta {\mathit{S}}^T\Delta \mathit{S})={\mathit{Q}}_{\mathit{S}}$.

Figure 1 illustrates that the platform can measure the azimuth angle of the target as ${\mathit{\theta }}_n={\mathit{\theta }}_n^o+\Delta {\mathit{\theta }}_n$ and the pitch angle as ${\mathit{\beta }}_n{=\mathit{\beta }}_n^o+\Delta {\mathit{\beta }}_n$, where ${\theta }_n^o$ and ${\beta }_n^o$ denote the true angle; the precise observation equation for the AOA obtained by platform $n$ is

$$\begin{array}{l}{\mathit{\theta }}_n^o=\arctan \left(\left(\mathit{y}-y_n^o\right)/\left(\mathit{x}-x_n^o\right)\right),\\ {\mathit{\beta }}_n^o=\arctan \left((z-z_n^o)/\sqrt{{\left(x-x_n^o\right)}^2+{\left(y-y_n^o\right)}^2}\right).\end{array}$$

(2)

We can represent the angle information measured by each moving platform in vector form

$$\mathit{\theta }={\mathit{\theta }}^o+\Delta \mathit{\theta },\mathit{\beta }={\mathit{\beta }}^{\mathit{o}}+\Delta \mathit{\beta }.$$

(3)

In (3), $\mathit{\theta }={\left[{\mathit{\theta }}_1,{\mathit{\theta }}_2,\dots ,{\mathit{\theta }}_N\right]}^T$ and $\mathit{\beta }={\left[{\mathit{\beta }}_1,{\mathit{\beta }}_2,\cdots ,{\mathit{\beta }}_N\right]}^T$ denote the measured angle vectors, ${\mathit{\theta }}^o=[{\mathit{\theta }}_1^o,{\mathit{\theta }}_2^o,\dots {\mathit{\theta }}_N^o]$ and ${\mathit{\beta }}^o={\left[{\mathit{\beta }}_1^o,{\mathit{\beta }}_2^o,\cdots ,{\mathit{\beta }}_N^o\right]}^T$ denote the true angle vectors, and $\Delta \mathit{\theta }={\left[\Delta {\mathit{\theta }}_1,\Delta {\mathit{\theta }}_2,\cdots ,\Delta {\mathit{\theta }}_N\right]}^T$ and $\Delta \mathit{\beta }={[\Delta {\mathit{\beta }}_1,\Delta {\mathit{\beta }}_2,\cdots ,\Delta {\mathit{\beta }}_N]}^T$ conform to a Gaussian distribution with a mean of zero. Their covariance matrices are, respectively, $E(\Delta {\mathit{\theta }}^T\Delta \mathit{\theta })={\mathit{Q}}_{\mathit{\theta }}$ and $E(\Delta {\mathit{\beta }}^T\Delta \mathit{\beta })={\mathit{Q}}_{\mathit{\beta }}$.

Using the first moving platform as the reference point, the exact observation equation for the TDOA obtained for the nth platform is

$${\tau }_n^o=\frac{1}{c^o}(R_n^o-R_1^o)=\frac{1}{c^o}\left(\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_n^{\mathit{o}}\parallel -\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1^{\mathit{o}}\parallel \right),n=2,3,\dots ,N,$$

(4)

where $R_n^o=\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_n^{\mathit{o}}\parallel$ denotes the true distance from the $n_{\mathit{th}}$ platform to the target.

Delay measurements typically incur errors. We denote ${\tau }_n={\tau }_n^o+\Delta {\tau }_n$. Measurements from multiple platforms can be represented in vector form, as

$$\mathit{\tau }={\mathit{\tau }}^{\mathit{o}}+\Delta \mathit{\tau }.$$

(5)

In (5), $\mathit{\tau }={\left[{\tau }_2,{\tau }_3,\cdots ,{\tau }_N\right]}^T$,${\mathit{\tau }}^o={\left[{\tau }_2^o,{\tau }_3^o,\cdots ,{\tau }_N^o\right]}^T$, respectively, represent the measured delay parameter vector and the actual parameter vector. $\Delta \mathit{\tau }={\left[\Delta {\tau }_2,\Delta {\tau }_3,\cdots ,\Delta {\tau }_N\right]}^T$, $\Delta \mathit{\tau }$ is a Gaussian error vector with a mean of zero, and its covariance matrix is $E(\Delta {\mathit{\tau }}^T\Delta \mathit{\tau })={\mathit{Q}}_{\mathit{\tau }}$.

Deriving the observation equation for the TDOA yields the exact observation equation for the FDOA.

$${\dot{\tau }}_n^o=\frac{1}{c^o}({\dot{R}}_n^o-{\dot{R}}_1^o),n=2,3,\dots ,N,$$

(6)

where ${\dot{R}}_n^o=\frac{{\left({\mathbf{u}}^o-{\mathbf{s}}_n^o\right)}^T\left({\dot{\mathbf{u}}}^o-{\dot{\mathbf{s}}}_n^o\right)}{R_n^o}$ represents the rate of the change in the radial distance difference between the $n_{\mathit{th}}$ platform and the target.

Frequency difference measurements typically contain errors. We denote ${\dot{\tau }}_n={\dot{\tau }}_n^o+\Delta {\dot{\tau }}_n$; the frequency difference measurements from multiple platforms can be represented in vector form

$$\dot{\mathit{\tau }}={\dot{\mathit{\tau }}}^{\mathit{o}}+\Delta \dot{\mathit{\tau }}.$$

(7)

In (7), $\dot{\mathit{\tau }}={\left[{\dot{\tau }}_2,{\dot{\tau }}_3,\cdots ,{\dot{\tau }}_N\right]}^T$ and ${\mathit{\tau }}^o={\left[{\dot{\tau }}_2^o,{\dot{\tau }}_3^o,\cdots ,{\dot{\tau }}_N^o\right]}^T$, respectively, represent the measured frequency difference parameter vector and the actual frequency difference parameter vector, $\Delta \dot{\mathit{\tau }}={\left[\Delta {\dot{\tau }}_2,\Delta {\dot{\tau }}_3,\cdots ,\Delta {\dot{\tau }}_N\right]}^T$, $\Delta \dot{\mathit{\tau }}$ represents the Gaussian distribution delay error vector with a mean of zero, and its covariance matrix is $E(\Delta {\dot{\mathit{\tau }}}^T\Delta \dot{\mathit{\tau }})=Q_{\dot{\mathit{\tau }}}$.

Combining the measurements of the AOA, TDOA, and FDOA and integrating them into an $(4N-2)$ matrix, we obtain

$$\mathit{\gamma }={\left[{\mathit{\theta }}^T,{\mathit{\beta }}^T,{\mathit{\tau }}^T,{\dot{\mathit{\tau }}}^T\right]}^T={\mathit{\gamma }}^o+\Delta \mathit{\gamma },$$

(8)

$${\mathit{Q}}_{\mathit{\gamma }}=E(\Delta {\mathit{\gamma }}^T\Delta \mathit{\gamma })=\mathrm{bladiag}\left({\mathit{Q}}_{\mathit{\theta }},{\mathit{Q}}_{\mathit{\beta }},{\mathit{Q}}_{\mathit{\tau }},{\mathit{Q}}_{\dot{\mathit{\tau }}}\right),$$

(9)

where ${\mathit{\gamma }}^o={\left[{\mathit{\theta }}^{oT},{\mathit{\beta }}^{oT},{\mathit{\tau }}^{oT},{\dot{\mathit{\tau }}}^{oT}\right]}^T$,$\Delta \mathit{\gamma }={\left[\Delta {\mathit{\theta }}^T,\Delta {\mathit{\beta }}^T,\Delta {\mathit{\tau }}^T,\Delta {\dot{\mathit{\tau }}}^T\right]}^T$, and ${\mathit{Q}}_{\mathit{\gamma }}$ is the measurement error covariance matrix.

3. Cramér–Rao Lower Bound Analysis

The CRLB reflects the optimal statistical performance obtained using all the unbiased estimates. Therefore, before presenting the proposed algorithm, this section first provides the CRLB of the localization model described in Section 2 as a standard for evaluating the estimation performance. And through simulation experiments, we verify the impact of an unknown sound speed and the uncertainty of the positions and velocities of the mobile platforms on the optimal estimation accuracy.

3.1. Derivation of the CRLB

The unknown vector is represented as $\mathit{\xi }={\left[{\mathit{S}}^{oT},c^o,{\mathit{\phi }}^{oT}\right]}^T$, where ${\mathit{S}}^o$ is a random parameter, ${\mathit{\phi }}^o={[{\mathit{u}}^T,{\dot{\mathit{u}}}^T]}^T$, and $c^o$ represents the unknown parameters. By combining the observation vector $\mathit{\gamma }={\left[{\mathit{\theta }}^T,{\mathit{\beta }}^T,{\mathit{\tau }}^T,{\dot{\mathit{\tau }}}^T\right]}^T$ with the platform parameter vector $\mathit{S}$ containing systematic errors, we can obtain the log-likelihood function of the joint probability density function of the target state estimates under the Gaussian data model as

$$\ln \mathit{P}\left(\mathit{\gamma },\mathit{S};\mathit{\xi }\right)=\mathit{F}-\frac{1}{2}{(\mathit{\gamma }-{\mathit{\gamma }}^o)}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\mathit{\gamma }-{\mathit{\gamma }}^o\right)-\frac{1}{2}{({\mathit{S}}^o-\mathit{S})}^T{\mathit{Q}}_{\mathit{S}}^{-1}\left({\mathit{S}}^o-\mathit{S}\right),$$

(10)

where $\mathit{F}=-0.5\ln ({(2\pi )}^{4\times N}\left|{\mathit{Q}}_{\mathit{\gamma }}\right|)-0.5\ln ({(2\pi )}^{6\times N}\left|{\mathit{Q}}_{\mathit{s}}\right|)$, and $\mathit{F}$ is a constant term independent of $\mathit{\xi }$.

According to the estimation theory [27], the CRLB estimate with respect to the unknown vector $\mathit{\xi }$ is obtained by taking the second-order derivative of the likelihood function with respect to $\mathit{\xi }$ and taking the inverse of the expectation matrix, which partitions the CRLB into a 3 × 3 block matrix, namely

$$\mathit{CLRB}\left(\mathit{\xi }\right)=-E{\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial \mathit{\xi }\partial {\mathit{\xi }}^T}\right]}^{-1}={\left[\begin{array}{ccc}{\mathit{X}}_{11}& {\mathit{X}}_{12}& {\mathit{X}}_{13}\\ {{\mathit{X}}_{12}}^T& {\mathit{X}}_{22}& {\mathit{X}}_{23}\\ {{\mathit{X}}_{13}}^T& {{\mathit{X}}_{23}}^T& {\mathit{X}}_{33}\end{array}\right]}^{-1}.$$

(11)

From (11), each block matrix corresponds to the estimated performance of ${\mathit{\phi }}^o,{\mathit{S}}^o \mathrm{and}{ c}^o$, and their expressions are

$$\begin{array}{l}{\mathit{X}}_{11}=-E\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial {\mathit{\phi }}^o\partial {\mathit{\phi }}^{oT}}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{\phi }}^o})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{\phi }}^o}\right),{\mathit{X}}_{12}=-E\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial {\mathit{\phi }}^o\partial c^o}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{\phi }}^o})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o}\right),\\ {\mathit{X}}_{22}=-\mathit{E}\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial c^o\partial c^o}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o}\right),{\mathit{X}}_{13}=-E\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial {\mathit{\phi }}^o\partial {({\mathit{S}}^{\mathit{o}})}^T}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{\phi }}^o})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}\right),\\ {\mathit{X}}_{23}=-E\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial c^o\partial {({\mathit{S}}^{\mathit{o}})}^T}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}\right),{\mathit{X}}_{33}=-E\left[\frac{{\partial }^2\mathit{ln}P\left(\mathit{\gamma };\mathit{\xi }\right)}{\partial {\mathit{S}}^{\mathit{o}}\partial {({\mathit{S}}^{\mathit{o}})}^T}\right]={(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}})}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left(\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}\right)+{\mathit{Q}}_{\mathit{S}}^{-1}.\end{array}$$

(12)

A more detailed derivation of the expression of each partial derivative in (12) is provided in Appendix A.

Next, we examine how the speed of sound propagation and perturbations in the position and velocity of the moving platform affect the localization performance. To emphasize the impact of two conditions—the unknown speed of sound and the presence of systematic errors—on the localization performance, we explore the CRLB with and without systematic errors and with knowledge of the unknown speed of sound. We examine this under two scenarios.

1. $c^o$ is known: When the speed of sound is known, the partial derivative of (11) with respect to $c^o$ is zero, as follows:

$$\mathrm{CRLB}\left({\mathit{\phi }}^o\right)={\left[\begin{array}{cc}{\mathit{X}}_{11}& {\mathit{X}}_{13}\\ {{\mathit{X}}_{13}}^T& {\mathit{X}}_{33}\end{array}\right]}^{-1}.$$

(13)

Using our knowledge of matrix chunking theory [37], the CRLB of the unknown vector ${\mathit{\phi }}^o$ can be obtained as:

$$\begin{array}{ll}\mathrm{CLRB}\left({\mathit{\phi }}^o\right)& ={({\mathit{X}}_{11}-{\mathit{X}}_{13}{\mathit{X}}_{33}^{-1}{\mathit{X}}_{13}^T)}^{-1}\\ & ={{\mathit{X}}_{11}}^{-1}+{{\mathit{X}}_{11}}^{-1}{\mathit{X}}_{13}{({\mathit{X}}_{33}-{\mathit{X}}_{13}^T{{\mathit{X}}_{11}}^{-1}{\mathit{X}}_{13})}^{-1}{\mathit{X}}_{13}^T{{\mathit{X}}_{11}}^{-1}.\end{array}$$

(14)

In (14), ${{\mathit{X}}_{11}}^{-1}$ represents the CRLB when the speed of sound is known and there are no systematic errors, and the second term on the right-hand side indicates the performance loss due to the presence of systematic errors.

2. $c^o$ is unknown: When the speed of sound is completely unknown, it is possible to jointly estimate the unknown information on the speed of sound and the observed target, and we denote ${\mathit{\nu }}^o={\left[{\mathit{\phi }}^{oT},c^o\right]}^T$, obtaining the CRLB of ${\mathit{\nu }}^o$ as

$$\begin{array}{ll}\mathrm{CRLB}({\mathit{\nu }}^o)& ={\left(\left[\begin{array}{ll}{\mathit{X}}_{11}\hfill & {\mathit{X}}_{12}\hfill \\ {\mathit{X}}_{12}^{\mathit{T}}\hfill & {\mathit{X}}_{22}\hfill \end{array}\right]-\left[{\mathit{X}}_{13}{\mathit{X}}_{23}\right]{\mathit{X}}_{33}^{-1}\left[\begin{array}{l}{\mathit{X}}_{13}^{\mathit{T}}\hfill \\ {\mathit{X}}_{23}^{\mathit{T}}\hfill \end{array}\right]\right)}^{-1}\\ & ={({\mathit{X}}_c-{\mathit{X}}_d{\mathit{X}}_{33}^{-1}{\mathit{X}}_d^T)}^{-1}={{\mathit{X}}_c}^{-1}+{{\mathit{X}}_c}^{-1}{\mathit{X}}_d{({\mathit{X}}_{33}-{\mathit{X}}_d^T{{\mathit{X}}_c}^{-1}{\mathit{X}}_d)}^{-1}{\mathit{X}}_d^T{{\mathit{X}}_c}^{-1},\end{array}$$

(15)

where ${\mathit{X}}_c={\left[\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial \mathit{\phi }} \frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o}\right]}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\left[\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial \mathit{\phi }} \frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o}\right]$, ${\mathit{X}}_d={\left[\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial \mathit{\phi }} \frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial c^o}\right]}^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}$.

The first term on the right-hand side of (15) incorporates the effect of the unknown acoustic propagation velocity on the CRLB of $\mathit{\phi }$. The second term then contains the effect of the unknown acoustic signal propagation velocity and the platform system errors on the CRLB of ${\mathit{\nu }}^o$, which is clearly not generated in an additive manner. If the effect of systematic errors is not considered, i.e., $\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}=0$, then the first term on the right-hand side of (15) is the CRLB estimate of ${\mathit{\nu }}^o$. At this point, the $\mathrm{CRLB}$ values of the unknown vectors ${\mathit{\phi }}^o$ and $c^o$ are

$$\left\{\begin{array}{l}\mathrm{CLRB}\left({\mathit{\phi }}^o\right)={{\mathit{X}}_{11}}^{-1}+{{\mathit{X}}_{11}}^{-1}{\mathit{X}}_{12}\mathrm{CLRB}\left(c^o\right){\mathit{X}}_{12}^T{{\mathit{X}}_{11}}^{-1}\\ \mathrm{CLRB}\left(c^o\right)={({\mathit{X}}_{22}-{\mathit{X}}_{12}^T{{\mathit{X}}_{11}}^{-1}{\mathit{X}}_{12})}^{-1}\end{array}\right..$$

(16)

The ultimate CRLB estimates for the unknown vectors ${\mathit{u}}^o,{\dot{\mathit{u}}}^o \mathrm{and}{ c}^o$ are derived as follows

$${\mathrm{CRLB}}_{\mathit{u}}={\displaystyle \sum }_{k=1}^3\mathrm{CLRB}{\left({\mathit{\nu }}^o\right)}_{k,k},{\mathrm{CRLB}}_{\dot{\mathit{u}}}={\displaystyle \sum }_{k=4}^6\mathrm{CLRB}{\left({\mathit{\nu }}^o\right)}_{k,k},{\mathrm{CRLB}}_{c^o}=\mathrm{CLRB}{\left({\mathit{\nu }}^o\right)}_{7,7}.$$

(17)

3.2. Comparison of CRLB Simulation Experiments

To analyze the impact of the known/unknown $c^o$ and the systematic errors of the mobile platform on the CRLB, we conduct an analysis here by adjusting various sensor measurement errors and mobile platform systematic errors. The simulation utilizes the positioning system described in reference [36], comprising underwater mobile platforms, to estimate the position, velocity, and unknown speed of sound for the underwater targets, as illustrated in Figure 2. Here, with $N=10$, the actual position and velocity parameters of the platform are presented in

Table 1. Real position and velocity parameters of moving platforms.

Platform $\mathit{no}.\mathit{n}$	1	2	3	4	5	6	7	8	9	10
$x_n^o$	0	0	0	0	1000	1000	1000	1000	500	500
${\mathit{y}}_n^o$	100	0	0	1000	1000	1000	1000	1000	500	500
$z_n^o$	0	0	1000	1000	0	0	1000	1000	1000	0
${\dot{x}}_n^o$	3	−3	1	1	2	2	1.2	−1.5	1.3	2.5
${\dot{\mathit{y}}}_n^o$	−2	1	−2	2	−1	−1	−1.5	1.2	1.3	2.5
${\dot{z}}_n^o$	2	−2	1	3	1	1	1.5	−1.2	1.3	2.5

.

The simulation results are presented in Figure 3 and Figure 4. Figure 3a,b depict the variation in the CRLB for the target’s position and velocity, respectively, with sensor measurement errors under $c^o$ unknown/known conditions. Setting the magnitude of the systematic errors to ${\sigma }_s=1 \mathrm{m}$, Figure 3a also illustrates the change in the CRLB of the target’s position with sensor measurement errors when systematic errors are absent. Regardless of the presence of systematic errors, the CRLB at an unknown sound speed is higher than that with the speed of sound known. For instance, when the perturbation parameter of the observables $\sigma =1$, i.e., $10\lg {\sigma }^2=0$, ${\sigma }_{\mathrm{aoa}}=0.1 \mathrm{deg}$, ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$, a difference of 1.12 dB is observed between the two in terms of the position estimation performance and 1.11 dB in terms of the speed estimation performance in the presence of systematic errors. In the absence of systematic errors, there is also a difference of 1.25 dB in the position estimation performance. Additionally, it can be observed from Figure 3a that the presence of systematic errors also diminishes the target localization accuracy when the sound speed is known/unknown.

Figure 4 shows the variation in the CRLBs incorporating the systematic error of the moving platform. The measurement error variance ${\sigma }_{\mathrm{aoa}}=1 \mathrm{deg}$, ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$; again, it can be seen that the CRLB at an unknown sound speed is higher than the CRLB with the speed of sound known. It can also be seen from Figure 4b that when the systematic error is less than −17 dB, the CRLB with systematic error and with the speed of sound known is higher than the position estimation error at an unknown sound speed and with no systematic error. This indicates that the smaller the error, the greater the effect of the speed of sound being unknown on localization.

In order to better verify the localization performance of moving targets when the speed of sound is unknown/known, it is assumed that the parameters of the underwater mobile observation platform are the same as above, ${\sigma }_{\mathrm{aoa}}=1 \mathrm{deg}$, ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$, ${\sigma }_s=1 \mathrm{m}$. The unknown moving target moves from a water area of ${[0,0,200]}^T\mathrm{m}$ from left to right at a speed of $1 \mathrm{m}/\mathrm{s}$ until it reaches a water area of ${[1000,1000,200]}^T\mathrm{m}$ with a localization period of 1 s. Figure 5 gives the CRLBs in the two cases with the moving targets at different positions. It can be seen that at the same position, the CRLB values in Figure 5b are generally larger than the CRLB values corresponding to Figure 5a.

In summary, the simulation also verifies the previous statement that the value of the estimated CRLB for the position and velocity of the target increases when the speed of sound is unknown and there is uncertainty in the information on the position and velocity of the moving platform. Therefore, it is necessary for us to design a joint AOA/TDOA/FDOA-based localization algorithm to achieve an asymptotically optimal statistical performance under these two conditions.

4. Positioning Algorithms

In this section, we will discuss the underwater mobile target localization problem with inaccurate mobile platform position and velocity information and completely unknown sound source signal velocity information. We propose a two-step closed-form localization algorithm. In the first stage, we minimize the introduction of additional unknowns to formulate pseudo-linear equations. After linearization, we obtain an initial solution. In the second stage, we derive an exact solution containing only unknown variables by considering the relationship between the source parameters and redundant variables. The specific process will be described below.

4.1. The First-Step Calculation Principle

First, we need to linearize the three joint localization equations so that we can construct a linear data model and thus roughly estimate the estimated target position ${\mathit{u}}^{\mathit{o}}$, velocity ${\dot{\mathit{u}}}^{\mathit{o}}$, and speed of sound $c^{\mathit{o}}$.

1.

For AOA: First, (2) can be transformed into a pseudo-linear equation

$$\left\{\begin{array}{l}\left(x-x_n^o\right)\sin \left({\theta }_n^o\right)-\left(y-y_n^o\right)\cos \left({\theta }_n^o\right)=0\Rightarrow {\mathit{\omega }}_1^{\mathit{T}}\left({\theta }_n^o\right)\mathit{u}={\mathit{\omega }}_1^{\mathit{T}}\left({\theta }_n^o\right){\mathit{s}}_{\mathit{n}}^{\mathit{o}},\hfill \\ =\left(\left(x-x_n^o\right)\cos \left({\theta }_n^o\right)+\left(y-y_n^o\right)\sin \left({\theta }_n^o\right)\right)\sin \left({\beta }_n^o\right)-\left(z-z_n^o\right)\cos \left({\beta }_n^o\right)=0\hfill \\ \Rightarrow {\mathit{\omega }}_2^{\mathit{{\rm T}}}\left({\theta }_n^o,{\beta }_n^o\right)\mathit{u}={\mathit{\omega }}_2^{\mathit{T}}\left({\theta }_n,{\beta }_n^o\right){\mathit{s}}_{\mathit{n}}^{\mathit{o}},\hfill \end{array}\right.$$

(18)

where ${\mathit{\omega }}_1\left({\theta }_n^o\right)={\left[-\sin \left({\theta }_n^o\right),\cos \left({\theta }_n^o\right),0\right]}^{\mathit{T}}$, ${\mathit{\omega }}_2\left({\theta }_n^o,{\beta }_n^o\right)=\left[-\sin \left({\mathit{\beta }}_n^o\right)\cos \left({\theta }_n^o\right),\right.$ ${\left.-\sin \left({\mathit{\beta }}_n^o\right)\sin \left({\theta }_n^o\right),\cos \left({\mathit{\beta }}_n^o\right)\right]}^T$.

The first-order Taylor expansion of $\sin \left({\theta }_n^o\right)$ and $\cos \left({\theta }_n^o\right)$ at ${\theta }_n$ and $\sin \left({\beta }_n^o\right)$ and $\cos \left({\beta }_n^o\right)$ at ${\mathit{\beta }}_n$, respectively, leads to the following equation

$$\left\{\begin{array}{l}\sin \left({\theta }_n^o\right)\approx \sin \left({\theta }_n\right)-\cos \left({\theta }_n\right)\Delta {\theta }_n,\cos \left({\theta }_n^o\right)\approx \cos \left({\theta }_n\right)+\sin \left({\theta }_n\right)\Delta {\theta }_n\\ \sin \left({\beta }_n^o\right)\approx \sin \left({\beta }_n\right)-\cos \left({\beta }_n\right)\Delta {\beta }_n,\cos \left({\beta }_n^o\right)\approx \cos \left({\beta }_n\right)+\sin \left({\beta }_n\right)\Delta {\theta }_n\end{array}\right..$$

(19)

Substituting (19) into (18) and performing a simple multiplicative transformation yields

$$\begin{array}{l}{\mathit{\omega }}_1{\left({\theta }_n\right)}^T{\mathit{s}}_{\mathit{n}}-{\mathit{\omega }}_1{\left({\theta }_n\right)}^T\mathit{u}=R_n\cos ({\beta }_n)\Delta {\theta }_n+{\mathit{\omega }}_1{\left({\theta }_n\right)}^T\Delta {\mathit{s}}_{\mathit{n}},\\ {{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right){\mathit{s}}_{\mathit{n}}-{{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right)\mathit{u}=R_n\Delta {\beta }_n+{{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right)\Delta {\mathit{s}}_{\mathit{n}}.\end{array}$$

(20)

2.

For TDOA: Instead of (4), we use the available moving platform position ${\mathit{s}}_{\mathit{n}}$, which can be obtained by performing a first-order Taylor expansion of $\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}}^{\mathit{o}}\parallel$ and $\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1^{\mathit{o}}\parallel$ at ${\mathit{s}}_{\mathit{n}}$ and ${\mathit{s}}_1$, respectively:

$$\begin{array}{l}\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}}^{\mathit{o}}\parallel \approx \parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}}\parallel +{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_{\mathit{n}}}^T\Delta {\mathit{s}}_{\mathit{n}},\\ \parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1^{\mathit{o}}\parallel \approx \parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1\parallel +{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}^T\Delta {\mathit{s}}_1.\end{array}$$

(21)

where ${\mathit{\rho }}_{\mathit{a},\mathit{b}}=(\mathit{a}-\mathit{b})/\Vert \mathit{a}-\mathit{b}\Vert$.

Substituting (21) into (4), the TDOA observation equation can be expressed as

$$c^o{\tau }_n=\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}}\parallel -\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1\parallel +\epsilon ,$$

(22)

where $\epsilon ={\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_{\mathit{n}}}^T\Delta {\mathit{s}}_{\mathit{n}}-{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}^T\Delta {\mathit{s}}_1+c^o\Delta {\tau }_n$. Moving $\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1\parallel$ from the right-hand side of (22) to the left side and squaring both sides simultaneously while neglecting second-order and higher-order error terms yields the following TDOA equation.

$${\mathit{s}}_1^T{\mathit{s}}_1-{\mathit{s}}_n^T{\mathit{s}}_n+2{({\mathit{s}}_n-{\mathit{s}}_1)}^T{\mathit{u}}^{\mathit{o}}+{\tau }_n^2{c}^{o2}+2{\tau }_n{R}_1{c}^o=2\parallel {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}}\parallel \epsilon .$$

(23)

The simplification can be obtained

$${\mathit{s}}_1^T{\mathit{s}}_1-{\mathit{s}}_n^T{\mathit{s}}_n+2{({\mathit{s}}_n-{\mathit{s}}_1)}^T{\mathit{u}}^{\mathit{o}}+{\tau }_n^2{c}^{o2}+2{\tau }_n{R}_1{c}^o=2{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}})}^T\Delta {\mathit{s}}_{\mathit{n}}+2R_n{c}^o\Delta \tau -2\frac{R_n}{R_1}{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1)}^T\Delta {\mathit{s}}_1.$$

(24)

3.

For FDOA: Taylor expansion of ${\dot{R}}_1^o={\left({\mathbf{u}}^o-{\mathbf{s}}_1^o\right)}^T\left({\dot{\mathbf{u}}}^o-{\dot{\mathbf{s}}}_1^o\right)/R_1^o$ at ${\mathit{s}}_1$ and ${\dot{\mathit{s}}}_1$ can be obtained

$${\dot{R}}_1^o\approx {\dot{R}}_1+{\mathit{q}}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}^T\Delta {\mathit{s}}_1+{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}^T\Delta {\dot{\mathit{s}}}_1,$$

(25)

where ${\dot{R}}_1=\left(\frac{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1){({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_1)}^T}{\Vert {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1\Vert }\right)$, ${\mathit{q}}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}=\left({\mathit{I}}_3-\left(\frac{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1){({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1)}^T}{\Vert {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1{\Vert }^2}\right)\right)\frac{({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_1)}{\Vert {\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1\Vert }$.

Then, we move $R_1^o$ and ${\dot{R}}_1^o$ to the left-hand side of Equations (4) and (6), respectively, to obtain

$$c{\tau }_n^o+R_1^o=R_i^o,c{\dot{\tau }}_n^o+{\dot{R}}_1^o={\dot{R}}_i^o.$$

(26)

Multiplying both sides of Equation (26) separately leads to the transformed FDOA observation equation system as

$$c^{o2}{\tau }_n^o{\dot{\tau }}_n^o+c^o{\dot{\tau }}_n^o{R}_1^o+c^o{\tau }_n^o{\dot{R}}_1^o={\dot{R}}_i^o{R}_i^o-{\dot{R}}_1^o{R}_1^o.$$

(27)

We substitute ${\tau }_n={\tau }_n^o+\Delta {\tau }_n$, ${\dot{\tau }}_n={\dot{\tau }}_n^o+\Delta {\dot{\tau }}_n$, ${\dot{\mathit{s}}}_{\mathit{n}}={\dot{\mathit{s}}}_{\mathit{n}}^{\mathit{o}}+\Delta {\dot{\mathit{s}}}_{\mathit{n}}$, ${\mathit{s}}_{\mathit{n}}={\mathit{s}}_{\mathit{n}}^{\mathit{o}}+\Delta {\mathit{s}}_{\mathit{n}}$ along with (21) and (25) into Equation (27) and neglect second-order errors to obtain

$$\begin{array}{l}{c}^o{R}_n\Delta {\dot{\tau }}_n+c^o{\dot{R}}_n\Delta {\tau }_n+{({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_n)}^T\Delta {\mathit{s}}_n-{({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_1+c^o{\tau }_n{\mathit{q}}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}+c^o{\dot{\tau }}_n{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1})}^T\Delta {\mathit{s}}_1\\ +{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}})}^T\Delta {\dot{\mathit{s}}}_{\mathit{n}}-{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1+c^o{\tau }_n{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1})}^T\Delta {\dot{\mathit{s}}}_1\\ \approx {\mathit{s}}_1^T{\dot{\mathit{s}}}_1-{\mathit{s}}_n^T{\dot{\mathit{s}}}_n+{({\dot{\mathit{s}}}_n-{\dot{\mathit{s}}}_1)}^T{\mathit{u}}^{\mathit{o}}+{({\mathit{s}}_n-{\mathit{s}}_1)}^T{\dot{\mathit{u}}}^{\mathit{o}}+{\tau }_n{\dot{\tau }}_n{c}^{o2}+{\dot{\tau }}_n{R}_1{c}^o+{\tau }_n{\dot{R}}_1{c}^o.\end{array}$$

(28)

Please note that it is widely assumed in the literature that neglecting second-order compound error terms is effective under small error conditions and has been confirmed in references [31,32,33,35,36].

4.

Fusion and pseudo-linear transformation of AOA/FDOA/TDOA measurements: First, collect the unknown target source parameters and two introduced auxiliary variable parameters into a vector ${\mathit{\psi }}_1^o$, where ${\mathit{\psi }}_1^o=\left[{\mathit{u}}^{\mathit{o}},{\dot{\mathit{u}}}^{\mathit{o}},c^{o2},R_1{c}^o,{\dot{R}}_1{c}^o\right]$; then, establish a pseudo-linear equation based on this parameter vector.

First, representing Equations (20), (24), and (28) in matrix form, we have

$${\mathit{h}}_1-{\mathit{A}}_1{\mathit{\psi }}_1^o={\mathit{B}}_1\Delta \mathit{\gamma }+{\mathit{D}}_1\mathit{\Delta }\mathit{S}.$$

(29)

In (29), ${\mathit{h}}_1=[{\mathit{h}}_{11}^T,{\mathit{h}}_{12}^T,{\mathit{h}}_{13}^T,{\mathit{h}}_{14}^T]$ represents a $(4N-2)\times 1$ observation vector, ${\mathit{A}}_1={[{\mathit{A}}_{11}^T,{\mathit{A}}_{12}^T,{\mathit{A}}_{13}^T,{\mathit{A}}_{14}^T]}^T$ represents a $(4N-2)\times 9$ observation vector coefficient matrix, and ${\mathit{D}}_1={[{\mathit{D}}_{11}^T,{\mathit{D}}_{12}^T,{\mathit{D}}_{13}^T,{\mathit{D}}_{14}^T]}^T$ represents a $(4N-2)\times 6N$ mobile platform system error coefficient matrix.

• where ${\mathit{h}}_{11}$ and ${\mathit{h}}_{12}$ are an $N\times 1$ matrix, ${\mathit{A}}_{11}$ and ${\mathit{A}}_{12}$ are an $N\times 9$ matrix, and ${\mathit{D}}_{11}$ and ${\mathit{D}}_{12}$ are an $N\times 6N$ matrix, $n=1,2,\dots N$. Their $n_{th}$ row matrices are

$${\mathit{h}}_{11}(n,1)=\left[{\mathit{\omega }}_1{\left({\theta }_n\right)}^T{\mathit{s}}_n\right],{\mathit{h}}_{12}(n,:)=\left[{{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right){\mathit{s}}_n\right],$$

(30)

$${\mathit{A}}_{11}(n,:)=\left[{\mathit{\omega }}_1{\left({\theta }_n\right)}^T,{\mathit{O}}_{1\times 6}\right],{\mathit{A}}_{12}(n,:)=\left[{{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right),{\mathit{O}}_{1\times 6}\right],$$

(31)

$${\mathit{D}}_{11}(n,:)=\left[{\mathit{O}}_{1\times 3(n-1)},{\mathit{\omega }}_1{\left({\theta }_n\right)}^T,{\mathit{O}}_{1\times (6N-3n)}\right],{\mathit{D}}_{12}(n,:)=\left[{\mathit{O}}_{1\times 3(n-1)},{{\mathit{\omega }}_2}^T\left({\theta }_n,{\beta }_n\right),{\mathit{O}}_{1\times (6N-3n)}\right],$$

(32)

where ${\mathit{h}}_{13}$ and ${\mathit{h}}_{14}$ are an $(N-1)\times 1$ matrix, ${\mathit{A}}_{13}$ and ${\mathit{A}}_{14}$ are an $(N-1)\times 9$ matrix, and ${\mathit{D}}_{13}$ and ${\mathit{D}}_{14}$ are ab $(N-1)\times 6N$ matrix, $n=2,3,\dots N$. Their ${(n-1)}_{th}$ row matrices are

$${\mathit{h}}_{13}(n-1,:)=\left[{\mathit{s}}_1^T{\mathit{s}}_1-{\mathit{s}}_n^T{\mathit{s}}_n\right],{\mathit{h}}_{14}(n-1,:)=\left[{\mathit{s}}_1^T{\dot{\mathit{s}}}_1-{\mathit{s}}_n^T{\dot{\mathit{s}}}_n\right],$$

(33)

$$\begin{array}{l}{\mathit{A}}_{13}(n-1,:)=-\left[2{\left({\mathit{s}}_n-{\mathit{s}}_1\right)}^T,{\mathit{O}}_{1\times 3},{\tau }_n^2,2{\tau }_n,0\right],\\ {\mathit{A}}_{14}(n-1,:)=-\left[{\left({\dot{\mathit{s}}}_n-{\dot{\mathit{s}}}_1\right)}^T,{\left({\mathit{s}}_n-{\mathit{s}}_1\right)}^T,{\tau }_n{\dot{\tau }}_n,{\dot{\tau }}_n,{\tau }_n\right],\end{array}$$

(34)

$$\begin{array}{ll}{\mathit{D}}_{13}(n-1,:)=& [-2\frac{R_n}{R_1}{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1)}^T,{\mathit{O}}_{1\times 3(n-2)},2{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}})}^T,{\mathit{O}}_{1\times (6N-3n)}],\\ {\mathit{D}}_{14}(n-1,:)=& \left[-{({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_1+c^o{\tau }_n{\mathit{q}}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1}+c^o{\dot{\tau }}_n{\mathit{\rho }}_{{\mathit{u}}^{\mathit{o}},{\mathit{s}}_1})}^T,{\mathit{O}}_{1\times 3(n-2)},{({\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_n)}^T,{\mathit{O}}_{1\times 3(N-n)},\right.\\ & \left.-\frac{R_n}{R_1}{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_1)}^T,{\mathit{O}}_{1\times 3(N-n)},{({\mathit{u}}^{\mathit{o}}-{\mathit{s}}_{\mathit{n}})}^T,{\mathit{O}}_{1\times (3N-3n)}\right].\end{array}$$

(35)

Similarly, in Equation (29), we have a $(4N-2)\times (4N-2)$ matrix of observation vector error coefficients, which can be represented in block matrix form as

$${\mathit{B}}_1=\mathrm{blkdiag}\left({\mathit{B}}_{11},{\mathit{B}}_{12},{\mathit{B}}_{13}\right),$$

(36)

where

$$\begin{array}{l}{\mathit{B}}_{11}=\mathrm{diag}\left(\right[R_1\cos ({\beta }_1),R_2\cos ({\beta }_2),\dots ,R_N\cos ({\beta }_N)\left]\right),\\ {\mathit{B}}_{12}=\mathrm{diag}\left(\right[R_1,R_2,\dots ,R_N\left]\right),\\ {\mathit{B}}_{13}=\left[\begin{array}{cc}\mathrm{diag}\left(2\right[R_2{c}^o,R_3{c}^o,\dots ,R_N{c}^o\left]\right)& {\mathit{O}}_{(N-1)\times (N-1)}\\ \mathrm{diag}\left(\right[c^o{\dot{R}}_2,c^o{\dot{R}}_3,\dots ,c^o{\dot{R}}_N\left]\right)& \mathrm{diag}\left(\right[R_2{c}^o,R_3{c}^o,\dots ,R_N{c}^o\left]\right)\end{array}\right].\end{array}$$

(37)

5.

Initial solution: Solving Equation (29) using weighted least squares (WLS) yields the estimated value ${\mathit{\psi }}_1$.

$${\mathit{\psi }}_1={({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1{\mathit{h}}_1,$$

(38)

where the weighting matrix is ${\mathit{W}}_1={({\mathit{B}}_1{\mathit{Q}}_{\mathit{\gamma }}{\mathit{B}}_1^T+{\mathit{D}}_1{\mathit{Q}}_S{\mathit{D}}_1^T)}^{-1}$.

According to (38), the calculation of the weighting matrix requires the use of an unknown parameter; we can replace ${\mathit{W}}_1$ with ${\mathit{Q}}_{\mathit{\gamma }}^{-1}$ to obtain the initial estimate and then update the weighted matrix using (32)–(37) to provide a more accurate estimate.

Then, we evaluate the covariance matrix of ${\mathit{\psi }}_1$, setting ${\mathit{\psi }}_1^o={({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1{\mathit{\psi }}_1$, with $\Delta {\mathit{\psi }}_1$ as the estimation error of the first-step algorithm, to obtain

$$\begin{array}{ll}\Delta {\mathit{\psi }}_1& ={\mathit{\psi }}_1-{\mathit{\psi }}_1^o={({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1({\mathit{h}}_1-{\mathit{A}}_1{\mathit{\psi }}_1^o)\\ & ={({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1({\mathit{B}}_1\Delta \mathit{\gamma }+{\mathit{D}}_1\mathit{\Delta }\mathit{S}).\end{array}$$

(39)

If the sensor measurement error and platform system error are sufficiently small to be disregarded, the estimation of ${\mathit{\psi }}_1^o$ in the initial step of the proposed algorithm can be deemed unbiased [27]. Subsequently, the error covariance matrix of the estimation value obtained ${\mathit{\psi }}_1$ in the first step of the algorithm can be approximated as

$$\mathrm{cov}({\mathit{\psi }}_1)=E(\Delta {\mathit{\psi }}_1\Delta {\mathit{\psi }}_1)={({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}.$$

(40)

4.2. The Second-Step Calculation Principle

In the first step of the algorithm, due to neglecting the relationship between the introduced auxiliary variables and the target variable, information is lost. Therefore, in this step, we improve the accuracy of the target estimation obtained in the first-step solution by utilizing the correlation among the elements in the vector of unknown parameters ${\mathit{\psi }}_1^o$. Initially, we express the true values of the unknown parameters in terms of the initial estimates in the first step as follows

$${\mathit{u}}^{\mathit{o}}={\mathit{\Psi }}_1(1:3)-\mathit{\Delta }{\mathit{\Psi }}_1(1:3),$$

(41)

$${\dot{\mathit{u}}}^o={\mathit{\Psi }}_1(4:6)-\mathit{\Delta }{\mathit{\Psi }}_1(4:6),$$

(42)

$$c^{o2}={\mathit{\Psi }}_1(7)-\mathit{\Delta }{\mathit{\Psi }}_1(7),$$

(43)

$$R_1{c}^o={\mathit{\Psi }}_1(8)-\mathit{\Delta }{\mathit{\Psi }}_1(8),$$

(44)

$${\dot{R}}_1{c}^o={\mathit{\Psi }}_1(9)-\mathit{\Delta }{\mathit{\Psi }}_1(9).$$

(45)

Moreover, there exists the following nonlinear relationship among the unknown vectors

$$c^o‖\mathit{u}-{\mathit{s}}_1‖=R_1{c}^o,$$

(46)

$$\frac{c^o{\left({\mathit{u}}^o-{\mathit{s}}_1\right)}^T\left({\dot{\mathit{u}}}^o-{\dot{\mathit{s}}}_1\right)}{R_1}={\dot{R}}_1{c}^o.$$

(47)

Substituting Equations (43) and (44) into Equation (46) and squaring both sides, while neglecting second- and higher-order terms, we obtain

$$\left({\mathit{u}}^T\mathit{u}+{\mathit{s}}_1^{\mathit{T}}{\mathit{s}}_1\right){\mathit{\Psi }}_1(7)-{\mathit{\Psi }}_1^2(8)-2{\mathit{s}}_1^{\mathit{T}}{\mathit{\Psi }}_1(7)\mathit{u}=\parallel \mathit{u}-{\mathit{s}}_1{\parallel }^2\mathit{\Delta }{\mathit{\Psi }}_1(7)-2{\hat{\mathit{\Psi }}}_1(8)\mathit{\Delta }{\mathit{\Psi }}_1(8).$$

(48)

Expanding ${\mathit{u}}^T\mathit{u}$ and $\parallel \mathit{u}-{\mathit{s}}_1{\parallel }^2$ around ${\mathit{\Psi }}_1(1:3)$ using Taylor series yields

$$\begin{array}{c}{\mathit{u}}^T\mathit{u}\approx {\mathit{\Psi }}_1^{\mathit{T}}(1:3){\mathit{\Psi }}_1(1:3)-2{\mathit{\Psi }}_1^{\mathit{T}}(1:3)\mathit{\Delta }{\mathit{\Psi }}_1(1:3),\\ \parallel \mathit{u}-{\mathit{s}}_1{\parallel }^2\approx \parallel {\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1{\parallel }^2-2{\left({\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1\right)}^T\mathit{\Delta }{\mathit{\Psi }}_1(1:3).\end{array}$$

(49)

Substituting (49) into (48) yields

$$\begin{array}{l}\left({\mathit{\Psi }}_1^{\mathit{T}}(1:3){\mathit{\Psi }}_1(1:3)+{\mathit{s}}_1^{\mathit{T}}{\mathit{s}}_1\right){\mathit{\Psi }}_1(7)-{\mathit{\Psi }}_1^2(8)-2{\mathit{s}}_1^{\mathit{T}}{\mathit{\Psi }}_1(7)\mathit{u}\\ =2{\mathit{\Psi }}_1(7){\mathit{\Psi }}_1^{\mathit{T}}(1:3)\mathit{\Delta }{\mathit{\Psi }}_1(1:3)+\parallel {\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1{\parallel }^2\mathit{\Delta }{\mathit{\Psi }}_1(7)-2{\mathit{\Psi }}_1(8)\mathit{\Delta }{\mathit{\Psi }}_1(8).\end{array}$$

(50)

Similarly, we multiply both sides of (47) by $c^o$ at the same time and substitute (43)–(45) into (47) to obtain:

$${\mathit{\Psi }}_1(7)\left({\mathit{u}}^{\mathit{o}T}{\dot{\mathit{u}}}^{\mathit{o}}-{\dot{\mathit{s}}}_1^{\mathit{T}}\mathit{u}-{\mathit{s}}_1^{\mathit{T}}{\dot{\mathit{u}}}^{\mathit{o}}+{\mathit{s}}_1^{\mathit{T}}{\dot{\mathit{s}}}_1\right)={\left({\mathit{u}}^o-{\mathit{s}}_1\right)}^T\left({\dot{\mathit{u}}}^o-{\dot{\mathit{s}}}_1\right)\mathit{\Delta }{\mathit{\Psi }}_1(7)-{\mathit{\Psi }}_1(9)\mathit{\Delta }{\mathit{\Psi }}_1(8)-{\mathit{\Psi }}_1(8)\mathit{\Delta }{\mathit{\Psi }}_1(9).$$

(51)

Expanding ${\mathit{u}}^{\mathit{o}T}{\dot{\mathit{u}}}^{\mathit{o}}$ and ${\left({\mathit{u}}^o-{\mathit{s}}_1\right)}^T\left({\dot{\mathit{u}}}^o-{\dot{\mathit{s}}}_1\right)$ around ${\mathit{\Psi }}_1(1:3)$ and ${\mathit{\Psi }}_1(4:6)$ using Taylor series yields

$$\begin{array}{c}{\mathit{u}}^{\mathit{o}T}{\dot{\mathit{u}}}^{\mathit{o}}\approx {\mathit{\Psi }}_1^T(1:3){\mathit{\Psi }}_1(4:6)-{\mathit{\Psi }}_1^T(4:6)\mathit{\Delta }{\mathit{\Psi }}_1(1:3)-{\mathit{\Psi }}_1^T(1:3)\mathit{\Delta }{\mathit{\Psi }}_1(4:6),\\ \begin{array}{l}{\left({\mathit{u}}^o-{\mathit{s}}_1\right)}^T\left({\dot{\mathit{u}}}^o-{\dot{\mathit{s}}}_1\right)\approx {\left({\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1\right)}^T\left({\mathit{\Psi }}_1(4:6)-{\dot{\mathit{s}}}_1\right)-\\ {\left({\dot{\mathit{u}}}^o-{\dot{\mathit{s}}}_1\right)}^T\mathit{\Delta }{\mathit{\Psi }}_1(1:3)-{\left({\mathit{u}}^o-{\mathit{s}}_1\right)}^T\mathit{\Delta }{\mathit{\Psi }}_1(4:6).\end{array}\end{array}$$

(52)

Substituting (52) into (51) yields

$$\begin{array}{l}{\mathit{\Psi }}_1(7)({\mathit{\Psi }}_1^T(1:3){\mathit{\Psi }}_1(4:6)+{\mathit{s}}_1^{\mathit{T}}{\dot{\mathit{s}}}_1)-{\mathit{\Psi }}_1(9){\mathit{\Psi }}_1(8)-{\mathit{\Psi }}_1(7)({\dot{\mathit{s}}}_1^{\mathit{T}}{\mathit{u}}^{\mathit{o}}+{\mathit{s}}_1^{\mathit{T}}{\dot{\mathit{u}}}^{\mathit{o}})={\mathit{\Psi }}_1(7){\mathit{\Psi }}_1^T(4:6)\mathit{\Delta }{\mathit{\Psi }}_1(1:3)+\\ {\mathit{\Psi }}_1(7){\mathit{\Psi }}_1^T(1:3)\mathit{\Delta }{\mathit{\Psi }}_1(4:6)+{\left({\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1\right)}^T\left({\mathit{\Psi }}_1(4:6)-{\dot{\mathit{s}}}_1\right)\mathit{\Delta }{\mathit{\Psi }}_1(7)-{\mathit{\Psi }}_1(9)\mathit{\Delta }{\mathit{\Psi }}_1(8)-{\mathit{\Psi }}_1(8)\mathit{\Delta }{\mathit{\Psi }}_1(9).\end{array}$$

(53)

Combining (41)–(43) and (50) and (53), the matrix expression for the second step of the algorithm can be obtained

$${\mathit{h}}_2-{\mathit{A}}_2{\mathit{\psi }}_2^o={\mathit{B}}_2\mathit{\Delta }{\mathit{\Psi }}_1.$$

(54)

In (54), the unknown vector parameter is ${\mathit{\psi }}_2^o={[{\mathit{u}}^{\mathit{o}T},{\dot{\mathit{u}}}^{\mathit{o}T},c^{o2}]}^T$, ${\mathit{h}}_2$ is the $9\times 1$ measurement vector, ${\mathit{A}}_2$ is the $9\times 7$ measurement matrix, and ${\mathit{B}}_2$ is the $9\times 9$ error coefficient matrix, which is defined as

$$\begin{array}{l}{\mathit{h}}_2=\left[{\hat{\mathit{\Psi }}}_1^T(1:7),\left({\hat{\mathit{\Psi }}}_1^{\mathit{T}}(1:3){\hat{\mathit{\Psi }}}_1(1:3)+{\mathit{s}}_1^{\mathit{T}}{\mathit{s}}_1\right){\hat{\mathit{\Psi }}}_1(7)-{\hat{\mathit{\Psi }}}_1^2(8),\right.\\ \left.{\mathit{\Psi }}_1(7)\left({\mathit{\Psi }}_1^T(1:3){\mathit{\Psi }}_1(4:6)+{\mathit{s}}_1^{\mathit{T}}{\dot{\mathit{s}}}_1\right)-{\mathit{\Psi }}_1(9){\mathit{\Psi }}_1(8)\right],\end{array}$$

(55)

$${\mathit{A}}_2=\left[\begin{array}{ccc}{\mathit{I}}_{3\times 3}& {\mathit{O}}_{3\times 3}& {\mathit{O}}_{3\times 1}\\ {\mathit{O}}_{3\times 3}& {\mathit{I}}_{3\times 3}& {\mathit{O}}_{3\times 1}\\ {\mathit{O}}_{1\times 3}& {\mathit{O}}_{1\times 3}& 1\\ 2{\mathit{s}}_1^{\mathit{T}}{\hat{\mathit{\Psi }}}_1(7)& {\mathit{O}}_{1\times 3}& 0\\ {\dot{\mathit{s}}}_1^{\mathit{T}}{\hat{\mathit{\Psi }}}_1(7)& {\mathit{s}}_1^{\mathit{T}}{\hat{\mathit{\Psi }}}_1(7) & 0\end{array}\right],$$

(56)

$${\mathit{B}}_2=\left[\begin{array}{lllll}& \hfill {\mathit{I}}_{7\times 7}& & \hfill {\mathit{O}}_{7\times 2}& \\ 2{\hat{\mathit{\Psi }}}_1(7){\hat{\mathit{\Psi }}}_1^{\mathit{T}}(1:3)& \hfill {\mathit{O}}_{1\times 3}\hfill & \hfill \parallel {\hat{\mathit{\Psi }}}_1(1:3)-{\mathit{s}}_1{\parallel }^2\hfill & -2{\hat{\mathit{\Psi }}}_1(8)& \hfill 0\hfill \\ {\mathit{\Psi }}_1(7){\mathit{\Psi }}_1^T(4:6)& {\hat{\mathit{\Psi }}}_1(7){\hat{\mathit{\Psi }}}_1^{\mathit{T}}(1:3)& {\left({\mathit{\Psi }}_1(1:3)-{\mathit{s}}_1\right)}^T\left({\mathit{\Psi }}_1(4:6)-{\dot{\mathit{s}}}_1\right)& -{\mathit{\Psi }}_1(9)& -{\hat{\mathit{\Psi }}}_1(8)\end{array}\right].$$

(57)

Solving (50) using the WLS algorithm yields the estimate ${\mathit{\psi }}_2$ of ${\mathit{\psi }}_2^o$ as follows

$${\mathit{\psi }}_2={({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2{\mathit{h}}_2,$$

(58)

where the weighting matrix is ${\mathit{W}}_2={({\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T)}^{-1}$.

Let ${\mathit{\psi }}_2^o={({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2{\mathit{\psi }}_2^o$; subtracting ${\mathit{\psi }}_2^o$ from ${\mathit{\psi }}_2$, we obtain the estimation error $\Delta {\mathit{\psi }}_2$ of the second-step algorithm

$$\Delta {\mathit{\psi }}_2={({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2{\mathit{B}}_2\mathit{\Delta }{\mathit{\Psi }}_1.$$

(59)

From (59), it is evident that $\Delta {\mathit{\psi }}_2$ and $\Delta {\mathit{\psi }}_1$ are linearly correlated, implying unbiased estimation of the unknown variable ${\mathit{\psi }}_2^o$. Thus, the covariance matrix of the estimate ${\mathit{\psi }}_2$ obtained using the second-step algorithm can be approximated as follows

$$\mathrm{cov}({\mathit{\psi }}_2)=E(\Delta {\mathit{\psi }}_2\Delta {\mathit{\psi }}_2^T)={({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}.$$

(60)

Finally, the valuation of the position and velocity of the target source as well as the propagation speed of the signal are

$$\left\{\begin{array}{c}\mathit{u}={\mathit{\psi }}_2(1:3)\\ \dot{\mathit{u}}={\mathit{\psi }}_2(4:6)\\ c=\sqrt{{\mathit{\psi }}_2(7)}\end{array}\right..$$

(61)

4.3. Algorithmic Steps

We summarize the proposed algorithm in Algorithm 1.

Algorithm 1: The underwater moving target localization algorithm proposed in this paper

Input:  Information on the position and speed parameters of the moving platforms $\mathit{S}$, the observation vectors $\mathit{\gamma }={\left[{\mathit{\theta }}^T,{\mathit{\beta }}^T,{\mathit{\tau }}^T,{\dot{\mathit{\tau }}}^T\right]}^T$, the measurement error covariance matrix ${\mathit{Q}}_{\mathit{\gamma }}$, the platform systematic error covariance matrix ${\mathit{Q}}_{\mathit{s}}$.

First-step algorithm:

Compute the observation vector matrix ${\mathit{h}}_1$ and the pseudo-linear observation vector coefficient matrix ${\mathit{A}}_1$ using (32)–(35). Initialize ${\mathit{W}}_1={\mathit{Q}}_{\gamma }^{-1}$, and compute an initial estimate of ${\mathit{\Psi }}_1$ using (38). Calculate ${\mathit{B}}_1$ and ${\mathit{D}}_1$ using an initial estimate of ${\mathit{\Psi }}_1$ and (36)–(37). Recalculate the weighting matrix ${\mathit{W}}_1$. Calculate the new estimate of ${\mathit{\Psi }}_1$ using (38).

Second-step algorithm:

6. Calculate the error covariance matrix $\mathrm{cov}({\mathit{\psi }}_1)$ using (40). 7. Calculate ${\mathit{B}}_2$ using ${\mathit{\psi }}_1$ and (57); then, obtain the weighting matrix ${\mathit{W}}_2$. 8. Calculate ${\mathit{h}}_2$ using (55). 9. Calculate the estimation vector ${\mathit{\psi }}_2$ using (58). 10. Calculate $u,\dot{u},c$ using (61). Output: Target position and velocity and acoustic source signal velocity.

5. Algorithm Analysis

5.1. Algorithm Performance Analysis

This summary will testify to the theoretical performance of the proposed algorithm in realizing the CRLB. We will derive the covariance matrix of the unknown vector estimates ${\mathit{\psi }}_2$ using the proposed algorithm and compare and analyze it with (15) to determine whether it is equal. Firstly, the algorithm in this paper is applied under the conditions that the position error, velocity error, and sensor measurement error of the mobile platform are small, i.e., small errors, and therefore satisfy

$$\begin{array}{l}\mathit{C}1:‖\Delta {\theta }_n‖\ll ‖{\theta }_n^o‖,‖\Delta {\beta }_n‖\ll ‖{\beta }_n^o‖,\\ \mathit{C}2:‖\Delta {\tau }_n‖\ll ‖{\tau }_n^o‖,\\ \mathit{C}3:‖\Delta {\dot{\tau }}_n‖\ll ‖{\dot{\tau }}_n^o‖,\\ \mathit{C}4:‖\Delta s_n‖\ll ‖s_n^o‖,‖\Delta {\dot{s}}_n‖\ll ‖{\dot{s}}_n^o‖.\end{array}$$

(62)

Conditions C1 and C4 are achievable with small bearing and platform systematic errors; C2 can be met with small time delay differences or distant targets; and C3 is attainable with small frequency difference errors or slow-moving targets. In underwater environments, object speeds typically correspond to a few meters per second, thus making Conditions C2 and C3 suitable for underwater localization within a few kilometers. Furthermore, these error conditions are commonly employed in the existing literature [31,32,33,35,36].

We define the estimate of the unknown parameter as $\mathit{\nu }$ and the estimate of the target unknown parameter as $\mathit{\phi }$. Then, the unknown vector estimation error is $\Delta \mathit{\nu }={\left[\Delta {\mathit{\phi }}^T,\Delta c\right]}^T$, where $\Delta \mathit{\phi }=\mathit{\phi }-{\mathit{\phi }}^o$, $\Delta c=c-c^o$, which is carried into (61), and $\Delta \mathit{\phi }$ can be expressed according to the estimation error $\Delta {\mathit{\psi }}_2$ in the second step of the proposed algorithm as

$$\Delta \mathit{\phi }={\mathit{\psi }}_2(1:6)-{\mathit{\psi }}_2^o(1:6)=\Delta {\mathit{\psi }}_2(1:6).$$

(63)

Expanding ${\mathit{\psi }}_2^o(7)$ around $c$ using Taylor series yields

$$\Delta {\mathit{\psi }}_2(7)={\mathit{\psi }}_2(7)-{\mathit{\psi }}_2^o(7)=2c\Delta c.$$

(64)

Therefore, the unknown variable error $\Delta \mathit{\nu }$ can be related to $\Delta {\mathit{\psi }}_2$ as follows

$$\Delta \mathit{\nu }={\mathit{B}}_3\Delta {\mathit{\psi }}_2,$$

(65)

where ${\mathit{B}}_3=\mathrm{diag}\left(1,1,1,1,1,1,\frac{1}{2c}\right)$.

Under small error conditions, $\Delta \mathit{\nu }$ is linearly correlated with $\Delta {\mathit{\psi }}_2$ from (65), and the estimate of ${\mathit{\psi }}_2^o$ is unbiased, so the final estimate of the unknown parameter $\mathit{\nu }$ is unbiased. The covariance matrix of $\mathit{\nu }$ can be approximated as

$$\mathrm{cov}(\mathit{\nu })=E(\Delta \mathit{\nu }\Delta {\mathit{\nu }}^T)={\mathit{B}}_3\mathrm{cov}({\mathit{\psi }}_2){\mathit{B}}_3^T={\mathit{B}}_3{({\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T)}^{-1}{\mathit{B}}_3^T.$$

(66)

By substituting ${\mathit{W}}_1={({\mathit{B}}_1{\mathit{Q}}_{\mathit{\gamma }}{\mathit{B}}_1^T+{\mathit{D}}_1{\mathit{Q}}_S{\mathit{D}}_1^T)}^{-1}$ into (62) and utilizing the matrix inversion lemma [37], we obtain

$$\mathrm{cov}(\mathit{\nu })={({\mathit{B}}_3^{-T}{\mathit{A}}_2^T{\mathit{B}}_2^{-T}{\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1{\mathit{B}}_2^{-1}{\mathit{A}}_2{\mathit{B}}_3^{-1})}^{-1}.$$

(67)

According to the Sherman–Morrison–Woodbury formula, expanding ${\mathit{W}}_1$ yields

$${\mathit{W}}_1={\mathit{B}}_1^{-T}{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{B}}_1^{-1}-{\mathit{B}}_1^{-T}{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{B}}_1^{-1}{\mathit{D}}_1{({\mathit{Q}}_{\mathit{S}}^{-1}+{\mathit{D}}_1^T{\mathit{B}}_1^{-T}{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{B}}_1^{-1}{\mathit{D}}_1)}^{-1}{\mathit{D}}_1^T{\mathit{B}}_1^{-T}{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{B}}_1^{-1}.$$

(68)

Substituting (68) into (67) to obtain

$$\mathrm{cov}(\mathit{\nu })={({\mathit{G}}_1^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{G}}_1-{\mathit{G}}_1^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{G}}_2{({\mathit{Q}}_{\mathit{S}}^{-1}+{\mathit{G}}_2^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{G}}_2)}^{-1}{\mathit{G}}_2^T{\mathit{Q}}_{\mathit{\gamma }}^{-1}{\mathit{G}}_1)}^{-1},$$

(69)

where ${\mathit{G}}_1={\mathit{B}}_1^{-1}{\mathit{A}}_1{\mathit{B}}_2^{-1}{\mathit{A}}_2{\mathit{B}}_3^{-1}$,${\mathit{G}}_2={\mathit{B}}_1^{-1}{\mathit{D}}_1$.

Comparing (69) and (15), it is seen that they have the same form, so it is necessary to compare the relationship between ${\mathit{G}}_1$ and $\left[\frac{\partial {\mathit{\gamma }}^o}{\partial {\mathit{\phi }}^o} \frac{\partial {\mathit{\gamma }}^o}{\partial c^o}\right]$ and ${\mathit{G}}_2$ and $\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}$, which can be obtained using (58). We can obtain

$${\mathit{G}}_1\approx \left[\frac{\partial {\mathit{\gamma }}^o}{\partial {\mathit{\phi }}^o} \frac{\partial {\mathit{\gamma }}^o}{\partial c^o}\right],{\mathit{G}}_2\approx -\frac{\partial {\mathit{\gamma }}^{\mathit{o}}}{\partial {\mathit{S}}^{\mathit{o}}}.$$

(70)

According to (66), it can be shown that $\mathrm{cov}(\mathit{\nu })\approx \mathrm{CLRB}\left({\mathit{\nu }}^o\right)$, and the target position, velocity, and speed of sound estimated using the proposed algorithm can reach the CRLB under this localization model. The detailed calculation procedure for (66) is given in Appendix B.

5.2. Algorithmic Complexity Analysis

This subsection provides a theoretical analysis of the computational complexity of the proposed algorithm, focusing on the required number of real multiplications. Assuming $m\times p$ matrices are multiplied by $p\times q$ matrices, the computational effort is denoted as $mpq$. The computation of $n\times n$ matrix inversion is denoted by $n^3$. The computational complexity of the first step of the algorithm proposed in this paper mainly consists of calculating the weighting matrix ${\mathit{W}}_1$ and the initial solution ${\mathit{\psi }}_1$. Calculating ${\mathit{W}}_1$ requires the computation of several arithmetic units, which are ${\mathit{B}}_1{\mathit{Q}}_{\mathit{\gamma }}$, ${\mathit{B}}_1{\mathit{Q}}_{\mathit{\gamma }}{\mathit{B}}_1^T$, ${\mathit{D}}_1{\mathit{Q}}_{\mathit{S}}$, ${\mathit{D}}_1{\mathit{Q}}_{\mathit{S}}{\mathit{D}}_1^T$, and ${({\mathit{B}}_1{\mathit{Q}}_{\mathit{\gamma }}{\mathit{B}}_1^T+{\mathit{D}}_1{\mathit{Q}}_{\mathit{S}}{\mathit{D}}_1^T)}^{-1}$; solving for ${\mathit{\psi }}_1$ consists of ${\mathit{A}}_1^T{\mathit{W}}_1$, ${\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1$,${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}$, ${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1$, and ${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1{\mathit{h}}_1$. The computational complexity of the second step mainly focuses on the computation of ${\mathit{W}}_2$, which includes ${\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1)$, ${\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T$, ${({\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T)}^{-1}$, and ${\mathit{\psi }}_2$, which include ${\mathit{A}}_2^T{\mathit{W}}_2$, ${\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2$, ${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}$, ${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2$, and ${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2{\mathit{h}}_2$.

The complexity of each unit of the proposed algorithm combining the AOA/TDOA/FDOA is given in

Table 2. Complexity of each unit of the proposed algorithm.

Calculation Unit		Unit Complexity	Total Unit Complexity
${\mathit{W}}_1$	${\mathit{B}}_1{\mathit{Q}}_{\gamma }$	${(4N-2)}^3$	$\begin{array}{l}3{(4N-2)}^3+\\ 36(4N-2)N^2+\\ 6{(4N-2)}^{2}N\end{array}$
	${\mathit{B}}_1{\mathit{Q}}_{\gamma }{\mathit{B}}_1^T$	${(4N-2)}^3$
	${\mathit{D}}_1{\mathit{Q}}_S$	$36(4N-2)N^2$
	${\mathit{D}}_1{\mathit{Q}}_S{\mathit{D}}_1^T$	$6{(4N-2)}^{2}N$
	${({\mathit{B}}_1{\mathit{Q}}_{\gamma }{\mathit{B}}_1^T+{\mathit{D}}_1{\mathit{Q}}_S{\mathit{D}}_1^T)}^{-1}$	${(4N-2)}^3$
${\mathit{\psi }}_1$	${\mathit{A}}_1^T{\mathit{W}}_1$	$9{(4N-2)}^2$	$\begin{array}{l}9{(4N-2)}^2+\\ 171(4N-2)+\\ 729\end{array}$
	${\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1$	$81(4N-2)$
	${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}$	729
	${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1$	$81(4N-2)$
	${({\mathit{A}}_1^T{\mathit{W}}_1{\mathit{A}}_1)}^{-1}{\mathit{A}}_1^T{\mathit{W}}_1{\mathit{h}}_1$	$9(4N-2)$
${\mathit{W}}_2$	${\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1)$	729	2187
	${\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T$	729
	${({\mathit{B}}_2\mathrm{cov}({\mathit{\psi }}_1){\mathit{B}}_2^T)}^{-1}$	729
${\mathit{\psi }}_2$	${\mathit{A}}_2^T{\mathit{W}}_2$	567	1855
	${\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2$	441
	${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}$	343
	${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2$	441
	${({\mathit{A}}_2^T{\mathit{W}}_2{\mathit{A}}_2)}^{-1}{\mathit{A}}_2^T{\mathit{W}}_2{\mathit{h}}_2$	63

, and all the $(4N-2)$ matrices in the table need to be replaced by $(2N-2)$ matrices if only the computational complexity of the TDOA/FDOA observations is utilized.

Compared to references [22,32,36], our algorithm is optimized for computational complexity. Compared to [36], both algorithms have the same computational complexity for ${\mathit{W}}_1$ and ${\mathit{\psi }}_1$, both being ${(2N-2)}^2(12N+3)$ and $(2N-2)(36N^2+171)$. However, [36] requires iterations in solving ${\mathit{W}}_2$ and introduces more unknown variables, resulting in increased computational complexity in the second step compared to our algorithm. Ref [32] proposes a four-step solution using only TDOA observations, allowing estimation of the target position, but the four-step calculation process increases the overall computational load. The SDP solution proposed in [22] introduces $N$ additional unknowns and requires multiple iterations, making it the most computationally demanding. In summary, the computational complexity of the second-step correction in the proposed algorithm is lower than that of existing closed-form algorithms. This is because the additional variable count is reduced to 2 in the first step of our algorithm, the second step only considers the relationship between additional variables and unknowns, and the matrix dimension is independent of the number of platforms, thereby reducing the computational complexity.

6. Simulation Experiments

In this section, the theoretical analysis of the proposed algorithm is further validated through Monte Carlo (MC) simulation experiments. We compare the proposed algorithm with three classical methods: the SDP method [22], four-step WLS methods (FS-WLS) [32], and two-step WLS methods (TS-WLS) [36]. Although the localization scenarios for some algorithms differ from those in this paper, these three localization algorithms cover a wide range of scenarios, including the speed of sound being known or unknown and taking systematic errors into account or not, and possess a high degree of novelty themselves. Hence, they are highly comparable to the algorithms in this paper according to several key dimensions, facilitating a comprehensive and in-depth evaluation of the algorithms’ performance. The propagation speed of sound in seawater generally ranges between $1420 \mathrm{m}/\mathrm{s}$ and $1560 \mathrm{m}/\mathrm{s}$ [29]. Hence, the unknown speed of sound in each of the following simulation scenarios is randomly generated within the range of $(1420,1560) \mathrm{m}/\mathrm{s}$ with a uniform distribution. Additionally, all the simulation experiments are evaluated using the root mean square error (RMSE) in each simulation, calculated as follows:

$$\mathrm{RMSE}(\mathit{\nu })=\sqrt{{\displaystyle \sum _{n=1}^M{‖{\mathit{\nu }}^n-{\mathit{\nu }}^o‖}^2/M}},$$

(71)

where ${\mathit{\nu }}^n$ is the estimate of ${\mathit{\nu }}^o$ from the $n_{\mathit{th}}$ simulation experiment, and the total number of MC simulation experiments is $M=5000$. All the RMSE results are shown in dB; in addition, the square root CRLBs as the performance bounds are also given in the simulation. In addition, the simulation experiments are carried out using MATLAB2022 software.

6.1. Effect of Measurement Errors on Positioning Performance

This subsection will investigate the performance of the proposed algorithm as the sensor measurement error varies, consistent with the target’s position in Section 2. To emphasize the variation in the velocity, the speeds of high-speed underwater targets such as submarines and torpedoes are adopted as ${\dot{\mathit{u}}}^{\mathit{o}}={\left[-20,-15,-20\right]}^T\mathrm{m}/\mathrm{s}$ [31], with the systematic error set to ${\sigma }_s=1 \mathrm{m}$, and the observation disturbance parameter varies from −5 dB to 20 dB. Since the FS-WLS algorithm only uses TDOA observations, it estimates only the target position and sound speed. In the SDP method, the speed of sound is known. In the simulation, we take $c=c^o+20$ as the known sound speed and use the MATLAB toolbox CVX to solve the SDP problems [38]. Figure 6a–c, respectively, depict the RMSE curves of the position estimation, velocity estimation, and speed of sound estimation as the sensor measurement error varies.

Based on the simulation results, the positioning algorithm proposed in this paper can achieve estimates of the target position, velocity, and speed of sound within the entire range of ${\sigma }^2$ with corresponding accuracy in the CRLB, outperforming the other three algorithms. The SDP algorithm performs lower in its target position and velocity estimations compared to the other algorithms, indicating the significance of estimating the speed of sound to enhance the algorithm performance. From Figure 6a,b, it is evident that the FS-WLS algorithm exhibits a lower estimation RMSE under small error conditions because it only utilizes the TDOA for positioning and does not consider system errors. In Figure 6b, the velocity estimation RMSE of this algorithm is slightly higher than that of the FS-WLS algorithm due to the target velocity estimation being primarily influenced by the FDOA equation, with the target velocity being significantly smaller than the position value, resulting in a minimal difference in the velocity estimation between the two algorithms. Since the FS-WLS algorithm can estimate three parameters, the data from Figure 6a–c are extracted into

Table 3. Estimated RMSE of the TS-WLS algorithm and the algorithm in this paper.

Estimated Parameters	Position		Velocity		${\mathit{c}}^{\mathit{o}}$
${\sigma }^2$ (dB)	−5	20	−5	20	−5	20
TS-WLS (dB)	4.12	10.08	2.08	5.05	6.27	12.87
Proposed (dB)	0.66	6.48	2.06	4.99	5.10	11.39

for comparison with this algorithm. At an SNR of −5 dB, the TS-WLS algorithm differs from this algorithm by 3.46 dB in its position estimation RMSE, 0.02 dB in its sound speed estimation RMSE, and 1.17 dB in its velocity estimation RMSE; meanwhile, at an SNR of 20 dB, the difference in the position estimation RMSE between the two algorithms increases to 3.60 dB, the velocity estimation RMSE difference rises to 0.06 dB, and the sound speed estimation RMSE difference reaches 1.48 dB. As the measurement errors increase, the differences in the estimations between the two algorithms gradually grow, indicating that this algorithm demonstrates stronger robustness in response to varying measurement errors.

6.2. Algorithm for the Effect of Systematic Errors on Positioning Performance

In this subsection, we will consider the variation in the proposed algorithm’s performance with the systematic error of the mobile platform by setting the sensor observation perturbation parameter ${\sigma }_{\mathrm{aoa}}=1 \mathrm{deg}$ and ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$. The systematic error is varied from −5 dB to 20 dB, and the rest of the simulation conditions remain unchanged; the simulation results are shown in Figure 7a–c.

Similarly, the estimated RMSE of the proposed algorithm in this paper achieves the CRLB with variations in the systematic error and outperforms the other three algorithms. From Figure 7a, it is evident that when ${\sigma }_s^2\le 0 \mathrm{dB}$, the position estimation RMSE of the algorithm proposed in this paper is at least 6 dB lower than that of the SDP algorithm, which does not consider errors in the speed of sound, and at least 2 dB lower than that of the TS-WLS algorithm. Comparing Figure 7b,c, it can be seen that as the systematic error increases, the trends in the growth of the speed estimation RMSE of the SDP algorithm and the sound speed estimation RMSE of the FS-WLS algorithm are more significant than those of the other two algorithms. This is mainly due to the fact that the SDP algorithm ignores the effect of sound velocity conditions in its operation, while the FS-WLS algorithm does not fully consider systematic errors. It is also noted that the performance gap between the TS-WLS algorithm and the algorithm in this paper decreases because both algorithms ignore second-order errors. When the error increases to a certain point, this approximation method will no longer be accurate. Nevertheless, under small error conditions, the algorithm in this paper still outperforms the TS-WLS algorithm.

6.3. Error Cumulative Distribution Function (CDF) Comparison

The interaction between the mobile platform and the shape of the target typically affects the positioning performance. We will analyze in detail the performance of the proposed algorithm compared to existing algorithms using different position distributions of the target and the mobile platform by evaluating the cumulative distribution function (CDF) of the RMSE. We assume 10 mobile platforms and 1 target source are randomly distributed in an $\left(\right[-1000,1000]\times [-1000,1000]\times [-1000,1000\left]\right) {\mathrm{m}}^2$ underwater cubic space, with the speed parameters of the mobile platforms and the target randomly generated within the range of $[-5,5] \mathrm{m}/\mathrm{s}$ to simulate uncertainty in real-world environments. The experiment randomly generates 1000 scenarios, with the RMSE of the target parameter estimates in each scenario used as the simulation data. The sensor measurement error is set to ${\sigma }_{\mathrm{aoa}}=1 \mathrm{deg}$, ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$, the system error is set to ${\sigma }_s=1 \mathrm{m}$, and the MC simulation is carried out 5000 times for each scenario in the simulation.

Figure 8 depicts the cumulative distribution function (CDF) plots of the RMSEs for the target source position, velocity, and sound speed estimation. Observing Figure 8a, it is evident that the abscissa corresponding to the curve of the proposed algorithm is significantly smaller than that of the other three algorithms, even at the same CDF value. This indicates that in the simulation environment in this section, the RMSE in the target position estimation of the proposed algorithm is significantly smaller than that of the other three algorithms, even at the same probability level. Similarly, the conclusions drawn from Figure 8b,c indicate that the proposed algorithm performs comparably to the TS-WLS algorithm in velocity estimation but outperforms the SDP algorithm. Additionally, in sound speed estimation, the RMSE of the proposed algorithm is superior to that of the TS-WLS and FS-WLS algorithms.

6.4. Effect of the Number of Mobile Platforms on Localization Performance

The number of mobile platforms also affects the target localization performance. Therefore, this section evaluates the performance of the proposed algorithm under different numbers of mobile platforms. The simulation conditions are the same as those set in Section 6.3, with a randomly selected localization scenario. The number of mobile platforms is set to $N=4,5,\dots 10$.

Figure 9 depicts how varying the number of mobile platforms affects the localization performance in random scenarios. The simulation results show that the proposed algorithm closely approximates the CRLB and outperforms the other three algorithms. As illustrated in Figure 9a, the performance of the proposed algorithm is comparable to that of other algorithms with eight platforms when there are only four mobile platforms. With fewer mobile platforms, the performance gap between the proposed algorithm and the other algorithms widens. This is because the other algorithms introduce more auxiliary variables and have fewer measurements, which makes it impossible to obtain a unique solution with fewer observation platforms, leading to a greater performance loss. The algorithm proposed in this study utilizes AOA observations and only needs to meet a certain condition to achieve target localization.

6.5. Effect of Target Position on Localization Performance

The localization performance is significantly affected by the target’s location, especially its distance from the moving platform. Since the experiments focused on near-field targets, this subsection will investigate the impact of far-field targets on the proposed algorithm’s performance. Assuming the localization scenario described below, where the moving platform’s position and velocity parameters match those in Table 1, we define a far-field cubic region with dimensions within a range of $[2000,3000] \mathrm{m}$ along the x, y, and z axes and randomly place 50 target sources within this region. The systematic error is set to ${\sigma }_s=1 \mathrm{m}$, while ${\sigma }^2$ is adjusted from −5 dB to 20 dB. Additionally, 5000 MC simulations are conducted at each location. Figure 10a,b present boxplots illustrating the estimated RMSEs of the target source at 50 random positions and the estimated RMSEs of the corresponding speed of sound considering measurement errors, respectively.

Figure 10 illustrates that the box maximum and minimum values for the estimated RMSEs and CRLBs from our algorithm are remarkably close. For instance, when ${\sigma }^2=20 \mathrm{dB}$, the maximum RMSE for position estimation and its corresponding CRLB are 16.4537 m and 16.0758 m, respectively, with a difference of merely 0.3779 m. Similarly, the minimum values are 8.5482 m and 7.5654 m, differing by 0.9828 m, while the average values are 12.0175 m and 12.0354 m, differing by only 0.0179 m. In contrast, for the TS-WLS algorithm, the maximum, minimum, and mean values are significantly higher, at 90.643 m, 26.5601 m, and 62.9327 m, respectively, compared to 74.18 m, 18.0119 m, and 50.9152 m for our algorithm. Furthermore, Figure 10b reveals that our algorithm’s sound speed estimation RMSE deviates by 0.0863 m/s and 1.5868 m/s from the maximum values of the corresponding CRLB and TS-WLS algorithms, respectively. It also deviates by 0.0230 m/s and 0.0206 m/s from the minimum values and by 0.0719 m/s and 0.3002 m/s from the mean values. The analysis above clearly demonstrates that within a specific range of measurement errors, for various far-field target positions, our algorithm achieves estimated RMSE values for the target position and speed of sound that closely match the corresponding CRLB and outperform the TS-WLS algorithm.

The sensor’s error is set to ${\sigma }_{\mathrm{aoa}}=1 \mathrm{deg}$, ${\sigma }_{\mathrm{d}}=1 \mathrm{m}$, the systematic error is varied from −20 dB to 5 dB, and the other simulation parameters are unchanged.

Figure 11 depicts a boxplot of the target source parameter estimation RMSE as a function of the systematic error variation. It can be observed from Figure 11 that within the entire ${\sigma }_{\mathit{s}}^2$ range of variation, the maximum, minimum, and mean values of the boxplot of our algorithm are close to the corresponding CRLB, and all are lower than those of the TS-WLS algorithm. Therefore, we can conclude that within a reasonable range of systematic errors, the proposed algorithm demonstrates asymptotic optimality in its position and sound speed estimation for far-field target sources at different positions, and it significantly outperforms the TS-WLS algorithm.

7. Conclusions

In this paper, a new two-step closed-form localization algorithm with a joint AOA/TDOA/FDOA method is proposed to estimate the position and velocity of a target sound source and its signal propagation velocity given its uncertainty in an underwater environment, as well as the position and velocity errors for moving platforms. The algorithm first introduces auxiliary variables in the first stage, obtaining initial solutions by constructing pseudo-linear equations; then, in the second stage, based on the relationship between the auxiliary variables and unknown variables, it further obtains exact solutions containing only unknown variables. Through in-depth theoretical analysis and simulation experiments, the proposed algorithm achieves a localization accuracy comparable to the precision of the CRLB when the sensor measurement errors and mobile platform system errors are within reasonable ranges. The experiments also demonstrate that even in scenarios with fewer mobile platforms, uncertain platform–target geometries, and far-field sound source targets, the algorithm can still achieve effective localization. Moreover, the root mean square errors of the estimated target sound source position, velocity, and sound propagation speed are superior to those of existing algorithms, demonstrating its estimation accuracy.

This study considers the target sound source as the point source. In reality, the physical properties of the sounding object, like its shape and size, often influence the sound wave propagation, impacting the localization accuracy. Future research will delve into factors such as the shape and size of the target object, examining their effects on the localization performance and refining the localization algorithm accordingly.