An Offgrid DOA Estimation Method for Passive Sonar Detection Based on Iterative Proximal Projection
https://doi.org/10.1007/s11804024004190

Abstract
Traditional direction of arrival (DOA) estimation methods based on sparse reconstruction commonly use convex or smooth functions to approximate nonconvex and nonsmooth sparse representation problems. This approach often introduces errors into the sparse representation model,necessitating the development of improved DOA estimation algorithms. Moreover,conventional DOA estimation methods typically assume that the signal coincides with a predetermined grid. However,in reality,this assumption often does not hold true. The likelihood of a signal not aligning precisely with the predefined grid is high,resulting in potential grid mismatch issues for the algorithm. To address the challenges associated with grid mismatch and errors in sparse representation models,this article proposes a novel highperformance offgrid DOA estimation approach based on iterative proximal projection (IPP). In the proposed method,we employ an alternating optimization strategy to jointly estimate sparse signals and grid offset parameters. A proximal function optimization model is utilized to address nonconvex and nonsmooth sparse representation problems in DOA estimation. Subsequently,we leverage the smoothly clipped absolute deviation penalty (SCAD) function to compute the proximal operator for solving the model. Simulation and sea trial experiments have validated the superiority of the proposed method in terms of higher resolution and more accurate DOA estimation performance when compared to both traditional sparse reconstruction methods and advanced offgrid techniques.Article Highlights● A highprecision DOA estimation method is proposed based on the offgrid model. A comparison with other offgrid DOA estimation methods reveals that the proposed method exhibits superior sparse signal recovery performance, resulting in heightened DOA estimation accuracy.● The incorporation of offgrid models in the proposed method facilitates the mitigation of grid mismatch issues, thereby enhancing the accuracy of DOA estimation.● Simulation results and sea trial analysis indicate that the proposed method outperforms its counterparts in low SNR, a limited number of snapshots, and a reduced number of array elements. 
1 Introduction
The effectiveness of sonar in remotely detecting, locating, extracting, and recognizing target signals in underwater environments greatly depends on the signal processing technology employed by underwater acoustic arrays. Among these technologies, DOA estimation plays a pivotal role in the realm of underwater acoustic array signal processing. In recent decades, there has been a proliferation of superresolution DOA estimation algorithms tailored for underwater detection. Notable examples include the multiple signal classification (MUSIC) method (Wagner et al., 2021) and the estimation of signal parameter via rotational invariance technique (ESPRIT) method (Li et al., 2020).
While subspacebased DOA estimation algorithms like MUSIC can achieve superresolution DOA estimation, they are not without their limitations. For instance, these algorithms rely on eigenvalue decomposition of the data covariance matrix and the exhaustive search of all spatial angles, resulting in prohibitively high computational complexity (Guo et al., 2022). Furthermore, the majority of subspacebased DOA estimation algorithms exhibit limited adaptability when confronted with scenarios involving small snapshots, low SNR, and correlated signals. In contrast, sparse reconstructionbased DOA estimation algorithms excel in adapting to such challenging conditions.
In recent years, the application of sparse signal representation (Chen et al., 2001; Donoho et al., 2006) and compressed sensing (Donoho, 2006) has emerged as a powerful mathematical tool in the realm of DOA estimation. Notably, Malioutov et al. (2005) introduced a groundbreaking sparse representation model based on the l_{1}norm penalty for DOA estimation, commonly known as the l_{1}SVD method. This method entails performing singular value decomposition (SVD) on the data matrix to reduce its dimensionality and subsequently solving it through a secondorder cone (SOC) program. In parallel, Stoica et al. (2011) proposed the sparse iterative covariance basis estimation (SPICE) method for DOA estimation, centered around the minimization of covariance matching criteria. These sparse reconstructionbased DOA estimation methods exhibit superior performance compared to traditional DOA estimation algorithms. However, a noteworthy limitation is their presumption that the DOAs of signals align precisely with predetermined angular grids, which is often overly idealistic. Consequently, these methods frequently encounter grid mismatch issues, leading to a reduction in DOA estimation accuracy.
In addressing the challenge of grid mismatch, researchers have introduced DOA estimation methods based on an offgrid model (Zhu et al., 2011; Gretsistas and Plumbley, 2012; Jagannath and Hari, 2013; Tan et al., 2014; Wu et al., 2018; Yang et al., 2013; Zhang et al., 2019). Zhu et al. (2011) employed the perturbation parameter method in DOA estimation, resulting in the development of a sparse total least squares method. Nevertheless, this method presents difficulties in parameter selection. Yang et al. (2013) leveraged firstorder Taylor expansion to construct an offgrid model and introduced an offgrid sparse Bayesian inference (SBI) DOA estimation method. Tan et al. (2014) incorporated the structured error of the array manifold dictionary set into DOA estimation, proposing a model capable of jointly estimating sparse signals and grid offset errors. Zhang et al. (2019) adopted the iterative sparse projection (ISP) model (Sadeghi and BabaieZadeh, 2016) as a sparse signal recovery model, coupling it with the firstorder Taylor expansion offgrid model to create a highperformance DOA estimation method. Recently, Dai et al. (2024) introduced a secondorder Taylor expansion offgrid model, demonstrating enhanced DOA estimation accuracy under low SNR conditions when compared to the firstorder Taylor expansion offgrid model.
With the evolution of sparse signal recovery theory, there is potential to devise more advanced DOA estimation methods that exhibit superior performance. The sparse signal representation model introduced by Sadeghi and BabaieZadeh (2016) lacks extrapolation steps and employs smoothing functions to address nonconvex and nonsmooth problems, leaving room for enhancing its sparse signal representation model. In this regard, Ghayem et al. (2018) adopted nonconvex and nonsmooth functions to portray the sparse signal model. They employed the IPP method to resolve the sparse signal representation model. The IPP method includes an extrapolation step, further augmenting its performance. Chen et al. (2020) applied the IPP method as a sparse signal representation model in the context of DOA estimation for monostatic MultipleInput MultipleOutput (MIMO) radar, yielding commendable DOA estimation performance. However, it is worth noting that this method, as described, does not integrate with the offgrid model and consequently does not address grid mismatch issues.
In this study, we introduce the SCAD function as a nonsmooth sparse promotion function into the sparse representation model for DOA estimation in underwater detection. The devised DOA estimation algorithm leverages the IPP method to solve the sparse representation model. Hence, the sparse signal recovery method employed in this study exhibits superior capability in recovering sparse signals, leading to a notable enhancement in DOA estimation accuracy compared to ISP method employed in Zhang et al. (2019). Subsequently, we propose a novel offgrid DOA estimation algorithm by combining the sparse representation model with an offgrid model constructed through firstorder Taylor expansion. Simultaneously, the offgrid model, constructed through firstorder Taylor expansion, effectively addresses the grid mismatch problem not accounted for in the framework proposed in Chen et al. (2020). This significantly contributes to the ultimate achievement of accurate DOA estimation.
2 Offgrid DOA estimation model
Assuming the presence of Q farfield narrowband signals, denoted as s_{q}(t) for q = 1, 2, …, Q from Q distinct directions (θ_{1}, θ_{2}, …, θ_{Q}), and incident upon a uniform linear array (ULA) consisting of M elements. In this ULA, the separation between adjacent elements is set to half the wavelength of the signal. The kth snapshot signal data vector received by the ULA can be expressed as:
$$ \boldsymbol{x}(k)=\sum\limits_{q=1}^Q \boldsymbol{a}\left(\theta_q\right) s_q(k)+\boldsymbol{n}(k)=\boldsymbol{A}(\boldsymbol{\theta}) \boldsymbol{s}(k)+\boldsymbol{n}(k) $$ (1) where $\boldsymbol{x}(k)=\left[x_1(k), \cdots, x_M(k)\right]^{\mathrm{T}}$ represents the kth snapshot signal data vector received by ULA, $\boldsymbol{\theta}=\left[\theta_1, \theta_2, \cdots\right., \left.\theta_Q\right]^{\mathrm{T}}$, and n(k) signifies an independent identically distributed Gaussian white noise signal with a zero mean. $\boldsymbol{A}(\boldsymbol{\theta})=\left[\boldsymbol{a}\left(\theta_1\right), \cdots, \boldsymbol{a}\left(\theta_Q\right)\right]$ represents the array manifold matrix, where the qth steering vector of A(θ) is denoted as $\boldsymbol{a}\left(\theta_1\right)=\left[1, e^{j {\rm{ \mathsf{ π} }} \sin \theta_q}, \cdots, e^{j(M1) {\rm{ \mathsf{ π} }} \sin \theta_q}\right]^{\mathrm{T}}$.
The spatial range of [−90°, 90°] is divided into a set of search angles $\overline{\boldsymbol{\theta}}=\left\{\bar{\theta}_1, \cdots, \bar{\theta}_G\right\}$ with equally spaced angular intervals denoted as ‘b’, where $G \gg Q$ signifies the number of grid points. In cases where the qth DOA of the signal does not precisely align with a grid point, that is, $\theta_q=\bar{\theta}_{g q}+\beta_q \notin \overline{\boldsymbol{\theta}}$ and $\bar{\theta}_{g q} \in \overline{\boldsymbol{\theta}}$, we employ firstorder Taylor expansion to approximate the true steering vector of the qth signal. This approach enables us to:
$$ \boldsymbol{a}\left(\theta_q\right) \approx \boldsymbol{a}\left(\bar{\theta}_{g q}\right)+\boldsymbol{a}^{\prime}\left(\bar{\theta}_{g q}\right) \beta_q $$ (2) where $\left\beta_q\right \leqslant \frac{b}{2}$ signifies the absolute value of the associated grid offset, with θ_{gq} representing the nearest grid point to the DOA.
Let’s define Φ as a diagonal matrix, denoted as Φ = diag(β), where $\boldsymbol{\beta}=\left[\beta_1, \cdots, \beta_G\right]^{\mathrm{T}}$. Furthermore, there is $\boldsymbol{B}(\overline{\boldsymbol{\theta}})=\left[\boldsymbol{a}^{\prime}\left(\bar{\theta}_1\right), \cdots, \boldsymbol{a}^{\prime}\left(\bar{\theta}_N\right)\right]$, and the resulting modified array manifold dictionary set matrix is then:
$$ \hat{\boldsymbol{A}}(\boldsymbol{\beta})=\boldsymbol{A}(\overline{\boldsymbol{\theta}})+\boldsymbol{B}(\overline{\boldsymbol{\theta}}) \boldsymbol{\varPhi} $$ (3) When the array receives data from K snapshots, the matrix representing the received data by the array can be expressed as:
$$ \boldsymbol{X}=\hat{\boldsymbol{A}}(\boldsymbol{\beta}) \tilde{\boldsymbol{S}}+\boldsymbol{N} $$ (4) where $\boldsymbol{X} \boldsymbol{\epsilon} \mathbb{C}^{M \times K}, \tilde{\boldsymbol{S}} \boldsymbol{\epsilon} \mathbb{C}^{N \times K}$ denotes the unknown sparse signal matrix, and $\boldsymbol{N} \boldsymbol{\epsilon} \mathbb{C}^{M \times K}$ is the noise data matrix.
Here, we define an indicator function:
$$ \mathrm{I}(s)=\left\{\begin{array}{l} 0, s=0 \\ 1, \text { others } \end{array}\right. $$ (5) and the mixed norm of $\\tilde{\boldsymbol{S}}\_{2, 0}$ is
$$ \\tilde{\boldsymbol{S}}\_{2, 0}=\sum\limits_{g=1}^G I\left(\left\\tilde{\boldsymbol{S}}_{g, :}\right\_2\right) $$ (6) It is evident that among all grid points, only Q grid points exhibit signal presence. As a result, $\tilde{\boldsymbol{S}}$ constitutes a rowsparse signal matrix, which can be reconstructed employing sparse signal recovery theory. In other words:
$$ \min\limits_{\tilde{\boldsymbol{S}}, \boldsymbol{\beta}}\\tilde{\boldsymbol{S}}\_{2, 0} \text { s.t. }\\boldsymbol{X}[\boldsymbol{A}(\overline{\boldsymbol{\theta}})+\boldsymbol{B}(\overline{\boldsymbol{\theta}}) \boldsymbol{\varPhi}] \tilde{\boldsymbol{S}}\_F^2 \leqslant \varepsilon $$ (7) where ε represents the admissible noise threshold and ‖∙‖_{F} stands for the Frobenius norm. It is wellrecognized that minimizing the l_{0}norm poses an NPhard problem. To address this challenge, various approaches have been employed, such as the use of smooth convex functions like the l_{1}norm or nonconvex smoothing functions like the l_{p}norm to approximate it (Elad 2010). However, employing convex or smooth functions to represent nonconvex and nonsmooth sparse signal recovery problems can introduce model inaccuracies. Fortunately, recent advances in sparse signal recovery have indicated that the selection of appropriate nonconvex and nonsmooth sparse promotion functions can yield superior sparse recovery results compared to other relaxation techniques (Ghayem et al., 2018).
3 Proposed offgrid DOA estimation method for passive sonar detection based on IPP
In this section, we will employ the concept of alternating optimization iterations in conjunction with the IPP method to solve equation (7). Prior to commencing the i + 1th iteration, we have acquired (or set) the values of $\tilde{\boldsymbol{S}}_i$ and β_{i}. Subsequently, we can determine the values of $\tilde{\boldsymbol{S}}_{i+1}$ and β_{i + 1} for the i + 1th iteration through the following process:
$$ \tilde{\boldsymbol{S}}_{i+1}=\underset{\overline{\boldsymbol{S}}_i}{\arg \min }\left\\overline{\boldsymbol{S}}_i\right\_{2, 0}+\frac{1}{2}\left\\boldsymbol{X}\left[\boldsymbol{A}(\overline{\boldsymbol{\theta}})+\boldsymbol{B}(\overline{\boldsymbol{\theta}}) \boldsymbol{\varPhi}_i\right] \tilde{\boldsymbol{S}}_i\right\_F^2 $$ (8) and
$$ \boldsymbol{\beta}_{i+1}=\underset{\boldsymbol{\beta}_i}{\arg \min }\left\\boldsymbol{X}\left[\boldsymbol{A}(\overline{\boldsymbol{\theta}})+\boldsymbol{B}(\overline{\boldsymbol{\theta}}) \boldsymbol{\varPhi}_i\right] \tilde{\boldsymbol{S}}_{i+1}\right\_F^2 $$ (9) 3.1 Restore the unknown sparse signal matrix
To attain the sparse solution as indicated in equation (8), the DOA estimation problem in passive sonar can be reformulated into an optimization problem involving the proximal function (Parikh and Boyd, 2014). This transformation yields:
$$ \min\limits_{\bar{\boldsymbol{S}}_i, \boldsymbol{Z}} \mathrm{H}(\boldsymbol{Z})+\Gamma_{\mathrm{A}_{\varepsilon}}\left(\overline{\boldsymbol{S}}_i\right) s.t. \boldsymbol{Z}=\overline{\boldsymbol{S}}_i $$ (10) where H(Z) represents a nonconvex and nonsmooth sparse promotion function, while Γ_{Αε} denotes the indicator function (Eftekhari et al., 2009) of set $\mathrm{A}_{\varepsilon} \triangleq\left\{\overline{\boldsymbol{S}}_i:\right.\left.\left\\boldsymbol{X}\left[\boldsymbol{A}(\overline{\boldsymbol{\theta}})+\boldsymbol{B}(\overline{\boldsymbol{\theta}}) \boldsymbol{\varPhi}_i\right] \overline{\boldsymbol{S}}_i\right\_{\boldsymbol{F}}^2 \leqslant \varepsilon\right\}$. This leads to:
$$ \Gamma_{\mathrm{A}_{\varepsilon}}= \begin{cases}0, & \overline{\boldsymbol{S}}_i \in \mathrm{A}_{\varepsilon} \\ +\infty, & \text { others }\end{cases} $$ (11) By incorporating the norm penalty method (Nocedal and Wright, 1999) into the sparse signal representation model presented in equation (10), we have
$$ \min\limits_{\overline{\boldsymbol{S}}_i, \boldsymbol{Z}} \mathrm{H}(\boldsymbol{Z})+\Gamma_{\mathrm{A}_{\varepsilon}}\left(\overline{\boldsymbol{S}}_i\right)+\frac{1}{2 v}\left\\overline{\boldsymbol{S}}_i\boldsymbol{Z}\right\_F^2 $$ (12) in which υ > 0 denotes a penalty parameter. The alternating minimization method can solve (12), there is
$$ \left\{\begin{array}{l} \boldsymbol{Z}_{i+1}=\underset{\boldsymbol{Z}}{\arg \min } v \mathrm{H}(\boldsymbol{Z})+\frac{1}{2}\left\\overline{\boldsymbol{S}}_i\boldsymbol{Z}\right\_F^2 \\ \overline{\boldsymbol{S}}_{i+1}=\underset{\overline{\boldsymbol{S}}_i}{\arg \min } \Gamma_{\mathrm{A}_{\varepsilon}}\left(\overline{\boldsymbol{S}}_i\right)+\frac{1}{2}\left\\overline{\boldsymbol{S}}_i\boldsymbol{Z}_{i+1}\right\_F^2 \end{array}\right. $$ (13) In line with the proximal projection definition (Parikh and Boyd, 2014), equation (13) can be streamlined to:
$$ \left\{\begin{array}{l} \boldsymbol{Z}_{i+1}=\operatorname{prox}_{v H}\left(\overline{\boldsymbol{S}}_i\right) \\ \overline{\boldsymbol{S}}_{i+1}=\operatorname{prox}_{\Gamma_{\mathrm{A}_{\varepsilon}}}\left(\boldsymbol{Z}_{i+1}\right) \end{array}\right. $$ (14) where $\operatorname{prox}_{v H}\left(\overline{\boldsymbol{S}}_i\right)$ represents the proximal operator of υH(Z). Furthermore, the inclusion of an extrapolation step (Nesterov, 2004) significantly enhances the performance of the IPP algorithm, resulting in:
$$ \overline{\boldsymbol{S}}_i=\tilde{\boldsymbol{S}}_i+a\left(\tilde{\boldsymbol{S}}_i\tilde{\boldsymbol{S}}_{i1}\right) $$ (15) in which a ≥ 0 denotes a weighting constant. To express equation (14) in a more concise and compact form, we have:
$$ \tilde{\boldsymbol{S}}_{i+1}=\operatorname{prox}_{\Gamma_{\mathrm{A}_{\varepsilon}}}\left(\operatorname{prox}_{\mathrm{vH}}\left(\overline{\boldsymbol{S}}_i\right)\right)=\operatorname{P}_{\mathrm{A}_{\varepsilon}}\left(\operatorname{prox}_{\mathrm{vH}}\left(\overline{\boldsymbol{S}}_i\right)\right) $$ (16) where P_{Αε} represents the projection of set Α_{ε}.
In a study by Ghayem et al. (2018), it was showcased that the fusion of the SCAD function (Fan and Li, 2001) as a sparse promotion function with IPP methods outperforms other nonsmooth functions in terms of performance. For computational convenience, we employ $\overline{\boldsymbol{S}}_i^{l 2}$ to denote the column vector formed by the l_{2}norm of each row vector in matrix $\overline{\boldsymbol{S}}_i$. Consequently, we have:
$$ \operatorname{prox}_{\mathrm{vH}}\left(\overline{\boldsymbol{S}}_i^{l 2}\right)=\mathrm{F}_{\tau, w}^{\mathrm{scad}}\left(\overline{\boldsymbol{S}}_i^{l 2}\right) $$ (17) in which the $\boldsymbol{F}_{\tau, w}^{\mathrm{scad}}\left(\overline{\boldsymbol{S}}_i^{l 2}\right)$ represents the proximal mapping of SCAD function. And there is
$$ \mathrm{F}_{\tau, w}^{\text {scad }}\left(\overline{\boldsymbol{S}}_i^{l 2}\right) \triangleq \begin{cases}\operatorname{sign}\left(\overline{\boldsymbol{S}}_i^{l 2}\right)\left(\left\overline{\boldsymbol{S}}_i^{l 2}\right\tau\right)_{+} & \left\overline{\boldsymbol{S}}_i^{l 2}\right \leqslant 2 \tau \\ \frac{(w1) \overline{\boldsymbol{S}}_i^{l 2}\operatorname{sign}\left(\overline{\boldsymbol{S}}_i^{l 2}\right) w \tau}{w2} & 2 \tau <\left\overline{\boldsymbol{S}}_i^{l 2}\right \leqslant w \tau \\ \overline{\boldsymbol{S}}_i^{l 2} & \left\overline{\boldsymbol{S}}_i^{l 2}\right>w \tau\end{cases} $$ (18) in which w > 2, τ > 0 represents a threshold constant and $(s)_{+} \triangleq \max (s, 0)$. Upon obtaining the sparse vector through the solution of equation (16), the DOAs of the target can be determined by locating the position of the spectral peak.
Considering that the primary emphasis of this article lies in the application of the IPP algorithm to devise DOA estimation methods with robust sparse signal recovery capabilities, and given that prior mathematical studies (Ghayem et al., 2018) have extensively detailed the utilization of the SCAD function and proximal operator for sparse signal recovery within the IPP algorithm, this paper places greater focus on the conclusive nature of employing the SCAD function and proximal operator.
3.2 Update grid offset β
In order to obtain the grid offset parameter β_{i + 1} after acquiring the sparse signal matrix $\overline{\boldsymbol{S}}_{i+1}$, we solve equation (9). This can be expressed as follows:
$$ \boldsymbol{\beta}_{i+1}=\underset{\boldsymbol{\beta}_i}{\arg \min }\left\\boldsymbol{X}\left[\boldsymbol{A}_1+\boldsymbol{B}_1 \boldsymbol{\varPhi}_i\right] \tilde{\boldsymbol{S}}_{i+1}\right\_F^2 $$ (19) in which A_{1} denotes A(θ), B_{1} represents B(θ). Let’s define
$$ \Pi\left(\boldsymbol{\beta}_i\right)=\left\\mathbf{X}\left[\boldsymbol{A}_1+\boldsymbol{B}_1 \boldsymbol{\Phi}_i\right] \tilde{\boldsymbol{S}}_{i+1}\right\_F^2 $$ (20) and (19) is equivalent to:
$$ \boldsymbol{\beta}_{i+1}=\underset{\boldsymbol{\beta}_i}{\arg \min } \Pi\left(\boldsymbol{\beta}_i\right) $$ (21) In (21), Π(β_{i}) represents a convex function. In accordance with the principles of convex optimization theory, for a continuously quadratic differentiable and convex function Π(β_{i}), its minimum value is attained when
$$ \Delta \Pi\left(\boldsymbol{\beta}_i\right)=0 $$ (22) where Δ(∙) denotes derivative symbol. And easy to obtain
$$ \begin{gathered} \boldsymbol{\beta}_{i+1}= \\ \operatorname{Re}\left\{\left[\boldsymbol{B}_1^H \boldsymbol{B}_1 \circ\left(\tilde{\boldsymbol{S}}_{i+1} \tilde{\boldsymbol{S}}^H{ }_{i+1}\right)^T\right]^{1} \operatorname{diag}\left(\tilde{\boldsymbol{S}}_{i+1}\left(\mathbf{X}\boldsymbol{A}_1 \tilde{\boldsymbol{S}}_{i+1}\right)^H \boldsymbol{B}_1\right)\right\} \end{gathered} $$ (23) in which Re(.) represents the operator for extracting the real part, with the symbol○indicating the Hadamard product.
Table 1 presents a summary of the proposed offgrid DOA estimation method based on IPP. The key parameters outlined in Table 1 are as follows: a = 0.9, w = 30, I = 3, K = 0.9, $\tau_0=5 \cdot \max \left\tilde{\boldsymbol{S}}_0(:, l)\right$, in which $\left\tilde{\boldsymbol{S}}_0(:, l)\right$ denotes the modulus vector comprising all elements in the lth column vector of the matrix $\tilde{\boldsymbol{S}}_0$. It is crucial to note, as indicated in Table 2, the presence of an inner loop between the ith iteration and the i + 1th iteration. The steps involved in this inner loop are outlined in Table 1.
1: Input: $\boldsymbol{X}, \boldsymbol{A}_1, \boldsymbol{B}_1, Q, \tau_0, a, w, I, K$ 2: Initialization: $\tilde{\boldsymbol{S}}_0=\boldsymbol{A}_1^{+} \boldsymbol{X}, \tau=\tau_0, i=0$ 3: while i ≤ I do 4: i = i + 1 5: while τ > τ_{f} do 6: for n = 1, …, N do 7: $\overline{\boldsymbol{S}}_n=\tilde{\boldsymbol{S}}_n+a\left(\tilde{\boldsymbol{S}}_n\tilde{\boldsymbol{S}}_{n1}\right)$ 8: $\tilde{\boldsymbol{S}}_{n+1}=\mathbf{P}_{\mathrm{A}_{\varepsilon}}\left(\boldsymbol{F}_{\tau, w}^{\mathrm{scad}}\left(\overline{\boldsymbol{S}}_n^{l 2}\right)\right)$ 9: end for 10: τ = K∙τ 11: end while 12: $\tilde{\boldsymbol{S}}_{i+1}=\tilde{\boldsymbol{S}}_{n+1}$ 13: Closed form solution (23) of β_{i + 1} 14: end while 15: Output: $\tilde{\boldsymbol{S}}, \boldsymbol{\beta}$ 4 Analysis of simulation results
In this section, a significant number of simulation experiments were conducted to assess the performance of the proposed DOA estimation algorithm and to compare it with several stateoftheart algorithms, including l_{1}SVD (Malioutov et al., 2005), OGSBI (Yang et al., 2013), OGISP (Zhang et al., 2019), and IPP method (Chen et al., 2020). For all simulation experiments, an ULA with M = 18 sensors was employed to receive target signals. Unless otherwise specified, the number of signal snapshots in each simulation experiment was fixed at T = 100. Additionally, we considered Q = 2 uncorrelated farfield narrowband signals incident onto the ULA from angles of −10.3°and 10.3°, respectively. For all the simulation experiments, we maintained a grid spacing of b = 1° and conducted a total of MC=200 Monte Carlo operations. To quantitatively assess the DOA estimation accuracy of various algorithms, the root mean square error (RMSE) was employed as a performance metric, which is
$$ \mathrm{RMSE}=\frac{1}{\mathrm{MC}} \sum\limits_{p=1}^{\mathrm{MC}} \sqrt{\frac{1}{Q}\left\\overline{\boldsymbol{\theta}}^{(p)}\boldsymbol{\theta}^{(p)}\right\_2^2} $$ (24) where MC represents the number of Monte Carlo operations, θ^{(p)} represents the estimated DOAs, and θ^{(p)} represents the true DOAs in the pth experiment, respectively.
In Simulation 1, Figure 1 displays the normalized spatial spectrum of l_{1}SVD, OGSBI, OGISP and the proposed OGIPP methods. For this simulation, we considered Q = 2 uncorrelated farfield narrowband signals incident onto an ULA from angles of − 6.5°and 6.5°, respectively. The grid spacing was set to b = 1°, and the orientation of the target signal deviated by 0.5° from the nearest grid point. As illustrated in Figure 1, the l_{1}SVD method exhibited the poorest performance due to its lack of integration with the offgrid model, resulting in a deviation of 0.5° from the true signal direction. In comparison, the OGSBI method showed improved performance, with an estimated target signal orientation deviating by 0.1° from the true orientation. Furthermore, both the proposed method and OGISP accurately estimated the target orientation, aligning perfectly with the true orientation of the target signal, thereby demonstrating excellent DOA estimation performance.
In Simulation 2, the RMSE performance of l_{1}SVD, OGSBI, OGISP and the proposed OGIPP methods under various input SNR is depicted in Figure 2 through Figure 5. In Figures 2 and 3, the number of snapshots for each algorithm was set to T=50, with corresponding array elements of M=18 and M=20, respectively. From the observations made in Figure 2 and Figure 3, it is evident that the l_{1}SVD algorithm inherently suffers from grid bias, leading to an inability to enhance DOA estimation accuracy as SNR increases. Conversely, offgrid DOA estimation methods can overcome grid mismatch issues, with RMSE approaching 0 as SNR continues to rise. Notably, among the three offgrid DOA estimation methods, the proposed method exhibits the lowest RMSE, indicating the highest DOA estimation accuracy. In Figures 4 and 5, each algorithm was provided with T=100 snapshots, and the corresponding number of array elements was M=18 and M=20, respectively. Similar to Figures 2 and 3, the proposed method consistently demonstrated the lowest RMSE among the three offgrid DOA estimation methods in Figure 4 and Figure 5. Across these four distinct conditions, involving varying input SNR, snapshot number, and array element number, it can be concluded that the proposed method consistently exhibits the lowest RMSE curve, signifying superior DOA estimation accuracy.
In Simulation 3, Figure 6 illustrates the RMSE performance of the proposed method as a function of SNR with different numbers of array elements. For Simulation 3, the number of snapshots was set to T=50. As evident from Figure 6, an increase in the number of array elements results in a gradual improvement in the DOA estimation accuracy of the proposed method. In essence, the greater the number of hydrophone arrays, the larger the array aperture that can be formed, leading to improved spatial resolution and, consequently, enhanced DOA estimation performance. This phenomenon can be attributed to the benefits of waveform diversity gain.
In Simulation 4, the RMSE performance of IPP DOA estimation method (Chen et al., 2020) and the proposed OGIPP method under various input SNR is depicted in Figure 7. From Figure 7, it is evident that as the SNR increases, the accuracy of DOA estimation improves for both algorithms. Nonetheless, because the IPP method did not incorporate an offgrid model, the RMSE of DOA estimation remained constant even as the input SNR exceeded −8dB. In contrast, the proposed OGIPP method, which integrates offgrid models, demonstrates a reduced RMSE and the capability to mitigate grid mismatch issues. With a continued increase in the input SNR, the RMSE of DOA estimation gradually converges towards 0.
5 Data processing
To validate the algorithm’s effectiveness, underwater target detection experiments were conducted in the South China Sea in September 2020. A fiber optic hydrophone array system, featuring an ULA comprising twenty sensors, was deployed on the seabed. The ULA has an array element spacing of 6.25 meters and was positioned at a depth of approximately 120 meters within the designated sea area. The observation range of the ULA extended from [−90°, 90°]. According to GPS navigation data, the ship receiving the sound source is positioned at the azimuth of 65° relative to the array. Additionally, another type of experimental cooperative target ship is located around the azimuth of 50° relative to the array.
In this underwater target detection experiment, the processed data spanned a duration of 400 s, which was divided into 20 segments, each lasting 20 s. The dataset used in the experiment was gathered between 13:23 and 13:29 on September 24, 2020, during the midday hours. Throughout the data processing in the sea trial, a farfield narrowband signal model was employed to both describe and process the data. In the course of the experiment, the cooperative target ship at 50° transmitted a singlefrequency signal with a signal frequency of 73 Hz. Consequently, the narrowband signal frequency chosen for processing was set to 73 Hz. The array receiving data was sampled at a rate of 32 000 Hz, necessitating FFT transformations every 32 000 data points. For DOA estimation purposes, the algorithm generated spatial spectra plots every 20 snapshots in the frequency domain and produced BearingTime Records (BTRs) every 2 s.
During the processing of the sea trial data, four DOA estimation methods—CBF, OGSBI, OGISP, and OGIPP—were employed to analyze and determine the DOAs of targets within the sea trial area. Figure 8 presents the results of CBF processing on the sea trial data. It is evident that the main lobe of the CBF spatial spectrum is excessively broad, allowing only a rough determination of the target’s azimuth range. This limitation makes it challenging to precisely measure the direction of the target. Hence, the CBF algorithm is unable to differentiate between cooperative targets at 50° and ship targets receiving sound sources at 65°. Furthermore, the CBF algorithm detected a noncooperative target in the vicinity of the −20 degrees orientation. Figures 9 to 11 display the BTRs generated by the OGSBI, OGISP, and OGIPP algorithms, respectively. From Figures 9 to 11, it becomes apparent that all three algorithms accurately detected the strong interference targets at − 20° and the sound source receiving ship targets at 65°. Moreover, they also successfully detected a weaker target at 50°, a task for which the CBF algorithm demonstrated subpar azimuth estimation accuracy, failing to detect any target at 50°. Furthermore, it is worth noting that the proposed algorithm’s target tracking trajectory appears the narrowest among the three algorithms, indicating the highest DOA estimation performance among them.
6 Conclusions
This paper introduces a novel highperformance offgrid DOA estimation method based on IPP. In the proposed method, an alternating optimization strategy is utilized to jointly estimate sparse signals and grid offset parameters. We employ a proximal function optimization model to address the nonconvex and nonsmooth sparse representation challenges in DOA estimation. Additionally, the SCAD function is employed to derive the proximal operator for solving the model. As a result, the proposed algorithm exhibits robust sparse signal recovery capabilities and delivers more accurate DOA estimation results. Through simulations and sea trial experiments, the proposed approach has demonstrated superior resolution and more precise DOA estimation performance in comparison to traditional sparse classification methods and advanced offgrid methods. While the proposed algorithm demonstrates high accuracy in DOA estimation, it comes with a high computational complexity. There remains a gap in achieving realtime processing for seamless integration into the detection system. Additionally, there is potential for improvement in the firstorder Taylor expansion offgrid model. A future direction could involve exploring the use of multiorder Taylor expansion terms to design offgrid models, aiming to further enhance the accuracy of DOA estimation.
Competing interest Jingwei Yin is an editorial board member for the Journal of Marine Science and Application and was not involved in the editorial review, or the decision to publish this article. All authors declare that there are no other competing interests. 
Table 1 Offgrid DOA Estimation Based on IPP
1: Input: $\boldsymbol{X}, \boldsymbol{A}_1, \boldsymbol{B}_1, Q, \tau_0, a, w, I, K$ 2: Initialization: $\tilde{\boldsymbol{S}}_0=\boldsymbol{A}_1^{+} \boldsymbol{X}, \tau=\tau_0, i=0$ 3: while i ≤ I do 4: i = i + 1 5: while τ > τ_{f} do 6: for n = 1, …, N do 7: $\overline{\boldsymbol{S}}_n=\tilde{\boldsymbol{S}}_n+a\left(\tilde{\boldsymbol{S}}_n\tilde{\boldsymbol{S}}_{n1}\right)$ 8: $\tilde{\boldsymbol{S}}_{n+1}=\mathbf{P}_{\mathrm{A}_{\varepsilon}}\left(\boldsymbol{F}_{\tau, w}^{\mathrm{scad}}\left(\overline{\boldsymbol{S}}_n^{l 2}\right)\right)$ 9: end for 10: τ = K∙τ 11: end while 12: $\tilde{\boldsymbol{S}}_{i+1}=\tilde{\boldsymbol{S}}_{n+1}$ 13: Closed form solution (23) of β_{i + 1} 14: end while 15: Output: $\tilde{\boldsymbol{S}}, \boldsymbol{\beta}$ 
Chen J, Zheng Y, Zhang T, Chen S, Li J (2020) Iterative reweighted proximal projection based DOA estimation algorithm for monostatic MIMO radar. Signal Processing, 172: 107537. https://doi.org/10.1016/j.sigpro.2020.107537 Chen SS, Donoho DL, Saunders MA (2001) Atomic decomposition by basis pursuit. SIAM Rev 43(1): 129158. https://doi.org/10.1137/S003614450037906X Dai ZH, Zhang L, Wang C, Han X, Yin JW (2024) Enhanced secondorder offgrid DOA estimation method via sparse reconstruction based on extended coprime array under impulsive noise. IEEE Transactions on Instrumentation and Measurement 73: 117, 2024, Art No. 8500417. https://doi.org/10.1109/TIM.2023.3328069 Donoho DL (2006) Compressed sensing. IEEE Transactions on Information Theory 52(4): 12891306. https://doi.org/10.1109/TIT.2006.871582 Donoho DL, Elad M, Temlyakov VN (2006) Stable recovery of sparse overcomplete representations in the presence of noise. IEEE Transactions on Information Theory 52(1): 618. https://doi.org/10.1109/TIT.2005.860430 Eftekhari A, BabaieZadeh M, Jutten C, Moghaddam HA (2009) RobustSL0 for stable sparse representation in noisy settings. IEEE International Conference on Acoustics, Speech and Signal Processing, Taipei, China, 34333436 Elad M (2010) Sparse and Redundant Representations. New York, NY, USA: Springer Fan J, Li R (2001) Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association 96(456): 13481360. https://doi.org/10.1198/016214501753382273 Ghayem F, Sadeghi M, BabaieZadeh M, Chatterjee S, Skoglund M, Jutten C (2018) Sparse signal recovery using iterative proximal projection. IEEE Transactions on Signal Processing 66(4): 879894. https://doi.org/10.1109/TSP.2017.2778695 Gretsistas A, Plumbley MD (2012) An alternating descent algorithm for the offgrid DOA estimation problem with sparsity constraints. 2012 Proceedings of the 20th European Signal Processing Conference (EUSIPCO), Bucharest, Romania, 874878 Guo K, Guo LX, Li YS, Zhang L, Dai ZH, Yin JW (2022) Efficient DOA estimation based on variable least Lncosh algorithm under impulsive noise interferences. Digital Signal Processing 122: 103383. https://doi.org/10.1016/j.dsp.2021.103383 Jagannath R, Hari KVS (2013) Block sparse estimator for grid matching in single snapshot DOA estimation. IEEE Signal Processing Letters 20(11): 10381041. https://doi.org/10.1109/LSP.2013.2279124 Li W, Liao W, Fannjiang A (2020) Superresolution limit of the ESPRIT algorithm. IEEE Transactions on Information Theory 66(7): 45934608. https://doi.org/10.1109/TIT.2020.2974174 Malioutov D, Cetin M, Willsky AS (2005) A sparse signal reconstruction perspective for source localization with sensor arrays. IEEE Transactions on Signal Processing 53(8): 30103022. https://doi.org/10.1109/TSP.2005.850882 Nesterov Y (2004) Introductory Lectures on Convex Optimization: A Basic Course. Boston, MA, USA: Kluwer Nocedal J, Wright SJ (1999) Numerical Optimization. New York, NY, USA: Springer Parikh N, Boyd S (2014) Proximal algorithms. Foundations and Trends in Optimization 1(3): 127239. https://doi.org/10.1561/2400000003 Sadeghi M, BabaieZadeh M (2016) Iterative sparsificationprojection: fast and robust sparse signal approximation. IEEE Transactions on Signal Processing 64(21): 55365548. https://doi.org/10.1109/TSP.2016.2585123 Stoica P, Babu P, Li J (2011) SPICE: A sparse covariancebased estimation method for array processing. IEEE Transactions on Signal Processing 59(2): 629638. https://doi.org/10.1109/TSP.2010.2090525 Tan Z, Yang P, Nehorai A (2014) Joint sparse recovery method for compressed sensing with structured dictionary mismatches. IEEE Transactions on Signal Processing 62(19): 49975008. https://doi.org/10.1109/TSP.2014.2343940 Wagner M, Park Y, Gerstoft P (2021) Gridless DOA estimation and rootMUSIC for nonuniform linear arrays. IEEE Transactions on Signal Processing 69: 21442157. https://doi.org/10.1109/TSP.2021.3068353 Wu X, Zhu WP, Yan J, Zhang Z (2018) Two sparsebased methods for offgrid directionofarrival estimation. Signal Processing 142: 8795. https://doi.org/10.1016/j.sigpro.2017.07.004 Yang Z, Xie LH, Zhang C (2013) Offgrid direction of arrival estimation using sparse Bayesian inference. IEEE Transactions on Signal Processing 61(1): 3843. https://doi.org/10.1109/TSP.2012.2222378 Zhang XW, Jiang T, Li YS, Liu X (2019) An offgrid DOA estimation method using proximal splitting and successive nonconvex sparsity approximation. IEEE Access 7: 6676466773. https://doi.org/10.1109/ACCESS.2019.2917309 Zhu H, Leus G, Giannakis GB (2011) Sparsitycognizant total leastsquares for perturbed compressive sampling. IEEE Transactions on Signal Processing 59(5): 20022016. https://doi.org/10.1109/TSP.2011.2109956