A New Method for Joint Sparse DOA Estimation

Abstract

To tackle the issue of poor accuracy in single-snapshot data processing for Direction of Arrival (DOA) estimation in passive radar systems, this paper introduces a method for judiciously leveraging multi-snapshot data. This approach effectively enhances the accuracy of DOA estimation and spatial angle resolution in passive radar systems. Additionally, in response to the non-convex nature of the mixed norm, we propose a hyperbolic tangent model as a replacement, transforming the problem into a directly solvable convex optimization problem. The rationality of this substitution is thoroughly demonstrated. Lastly, through a comparative analysis with existing discrete grid DOA estimation methods, we illustrate the superiority of the proposed approach, particularly under conditions of medium signal-to-noise ratio, varying numbers of snapshots, and close target angles. This method is less affected by the number of array elements, and is more usable in practices verified in real-world scenarios.

1. Introduction

The Passive Bistatic Radar (PBR) system represents a variant of traditional radar systems with historical roots dating back decades [1,2,3]. PBR possesses advantages that active radar systems lack. Firstly, it can employ “non-cooperative radiation sources” as its energy transmission source, eliminating the need for dedicated transmitters [4,5,6]. This feature yields advantages such as cost-effectiveness, stealth, resistance to interference, and extensive application prospects. Secondly, PBR can operate across a wide frequency band without interfering with existing wireless systems, utilizing external radiation sources like television, radio stations, mobile phone base stations, and digital video broadcasting [7].

Recent research in passive radar has spawned various subfields, including direct wave purification and clutter suppression. Passive radar target localization and tracking are pivotal technologies, with localization methods categorized into direction-finding [8], frequency-azimuth joint [9], and direct localization methods [10]. Many localization methods require Direction of Arrival (DOA) information, and there are numerous DOA estimation methods. This paper primarily focuses on the DOA estimation method based on compressed sensing (CS) for multiple-shot data.

Numerous researchers have proposed algorithms based on sparse models, classifying sparse representations from different perspectives [11,12,13,14,15,16,17]. For instance, methods employed in sparse constraint with different norms can be categorized into five classes: $l_0$ norm minimization, $l_p$ norm minimization $\left(0<p<1\right)$, $l_1$ norm minimization, $l_2$ norm minimization, and $l_{2,1}$ norm minimization. Furthermore, various optimization algorithms for solving the above-mentioned sparse models can also be categorized, encompassing greedy strategies, constraint optimization strategies, approximation algorithms, and homotopy algorithms.

Addressing the challenge of developing new signal processing theories for spatial spectrum estimation and spatial filtering, some scholars have analyzed the interference effects of non-ideal factors in spatial spectrum estimation from different angles. They have developed a diagonal loading algorithm, directly seeking the optimal solution for the diagonal loading value. Moreover, recent Bayesian learning-based algorithms have been proposed [18]. This method excels in accurately reconstructing signals after sparse representation in array direction pattern synthesis and angle tracking.

In [19], Bin Hu et al. proposed an offline DOA estimation method based on compressed sensing in a multipath environment. The crucial step involves transforming the signal reception model into an Errors-in-Variables (EIV) model and then estimating the sparse coefficients using total least squares. The estimated sparse coefficients are utilized to update the grid, and subsequently, the updated grid and estimated sparse coefficients are employed to estimate the multipath coefficients. In [20], Linxi Liu et al. studied the DOA estimation problem of a one-dimensional uniform linear array system with unknown inter-antenna phase errors. Leveraging CS theory, they transformed the DOA estimation problem into a sparse reconstruction problem, capitalizing on the sparsity of signals in the spatial domain. Subsequently, they proposed an iterative DOA estimation method for updating sparse signals and phase errors. In [21], Berkan Kilic et al. introduced a hardware-efficient DOA estimation method based on compressed sensing. The hardware complexity of the DOA estimation system can be tailored by design. They also considered array imperfections in numerical analysis, anticipating that the proposed design would perform effectively in practical applications. Furthermore, the performance of their proposed method can be explored in various fields of signal processing where the number of samplers is less than the number of receivers.

2. Main Work and Symbol Description

The paper aims to address two key problems:

Utilization of Multi-Snapshot Data in Sparse Model-Based DOA Estimation: Many existing sparse model-based Direction of Arrival (DOA) estimation methods primarily focus on single-snapshot data. However, practical applications often involve multi-snapshot data, which tends to exhibit better sparsity recovery. Therefore, the paper emphasizes the need to find a reasonable approach for the effective utilization of multi-snapshot data.

Rationalization of Sparse Recovery Models for $l_{2,1}$ and $l_{2,0}$ norms: Joint sparse recovery problems, generalizations of single measurement vector problems in compressed sensing, aim to recover sets of jointly sparse vectors. The paper explores five norm-constrained sparse reconstruction algorithms, each with distinct characteristics. Notably, due to the NP-hard nature of the $l_{2,0}$ norm model, the paper emphasizes the importance of finding a replacement, potentially by deriving a theoretical analytical expression for the $l_{2,th}$ norm minimization problem addressed herein.

The main thrust of this paper centers on the rational utilization of multi-snapshot data in passive radar engineering. In addressing the non-convexity or non-differentiability issues inherent in sparse recovery models, particularly for $l_{2,1}$ and $l_{2,0}$ norms, the paper introduces a differentiable and optimizable surrogate model, the rationality of which is rigorously verified.

Symbol definitions used in this paper include:

$\parallel X{\parallel }_{g,h}$: The $g,h$ mixed norm of a matrix.

$X_{rowi}$: The $i$th row vector of a matrix.

$⟨x,y⟩$: The inner product of two vectors.

${‖X‖}_{2,th}$: The norm model replacing the $l_0$ norm with the hyperbolic tangent function.

${\displaystyle \frac{\partial {‖X‖}_{2,th}}{\partial x_{11}(u)}}$: The partial derivative of ${‖X‖}_{2,th}$ with respect to one of its variables.

3. Passive Radar DOA Estimation Model Based on Uniform Linear Array Antenna

In the Passive Radar System Model illustrated in Figure 1, the research scenario addressed in this paper involves a passive radar system designed to operate without emitting electromagnetic waves [22]. Instead, it harnesses the prevalent digital TV signals in space as the radiation source and employs a Uniform Linear Array (ULA) as the receiving station. The electromagnetic waves in space encompass not only target-reflected echoes, but also direct signals and multipath clutter, further compounded by noise within the receiving array. Given that most passive radar systems initially undertake clutter suppression and direct signal cancellation, the signal model proposed in this paper assumes a scenario where, post-elimination of multipath clutter and direct signals from the receiving array, only the target echoes and residual noise persist.

In Figure 2, a one-dimensional Uniform Linear Array (ULA) is depicted, comprising $M$ amplitude and phase-consistent elements selected as the reference element, with an element spacing of $d$. Considering $K$ far-field narrowband signals incident on the array from various directions, the signals received by the $m$th elements and the $k$th signals at time $t$ are expressed as [23]:

(1)$y_m={\displaystyle \sum _{k=0}^{K-1}{S}_k}\left(t-{\tau }_{m,k}\right)+n_m\left(t\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}m=1,2,\dots ,M$

Here, ${\tau }_{m,k}$ denotes the time delay of the signal $S_k\left(t\right)$ when it reaches the $m$th element relative to the main element. For the ULA, the delay is given by ${\tau }_{m,k}=\left(m-1\right)d\sin {\theta }_k/c$, where $c$ represents the speed of electromagnetic waves in a vacuum, generally considered as $3\times {10}^8,\mathrm{m}/\mathrm{s}$. The term $n_m\left(t\right)$ signifies the noise received by the $m$th element or generated by the array itself, with the assumption that all noises are independent Gaussian white noises.

The expressions for receiving signals at each element are:

(2)$y_m={\displaystyle \sum }_{k=0}^{K-1}{e}^{-j2\pi f_0{\tau }_{m,k}}{S}_k+n_m\left(t\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}m=1,2,\dots ,M$

where $f_0$ is the carrier frequency of the signal, $j=\sqrt{-1}$.

Writing the above equation in matrix form, we obtain the signals received by the uniform array at time [24]:

(3)$\mathit{y}\left(t\right)={\left[y_1\left(t\right),y_2\left(t\right),\dots ,y_M\left(t\right)\right]}^T$

(4)$\begin{array}{c}\mathit{y}\left(\mathit{t}\right)=\left[\begin{array}{cccc}1& 1& \dots & 1\\ {\mathit{e}}^{-j(2\pi d\sin {\theta }_1/\lambda )}& {\mathit{e}}^{-j(2\pi d\sin {\theta }_2/\lambda )}& \dots & {\mathit{e}}^{-j(2\pi d\sin {\theta }_K/\lambda )}\\ \dots & \dots & \dots & \dots \\ {\mathit{e}}^{-j(2\pi \left(M-1\right)d\sin {\theta }_1/\lambda )}& {\mathit{e}}^{-j(2\pi \left(M-1\right)d\sin {\theta }_2/\lambda )}& \dots & {\mathit{e}}^{-j(2\pi \left(M-1\right)d\sin {\theta }_K/\lambda )}\end{array}\right]\hfill \\ \ast \left[\begin{array}{c}{\mathit{S}}_{\mathbf{1}}\left(\mathit{t}\right)\\ {\mathit{S}}_{\mathbf{2}}\left(\mathit{t}\right)\\ \dots \\ {\mathit{S}}_{\mathit{K}}\left(\mathit{t}\right)\end{array}\right]+\left[\begin{array}{c}{\mathit{n}}_{\mathbf{1}}\left(\mathit{t}\right)\\ {\mathit{n}}_{\mathbf{2}}\left(\mathit{t}\right)\\ \dots \\ {\mathit{n}}_{\mathit{K}}\left(\mathit{t}\right)\end{array}\right]\hfill \end{array}$

The above equation can be expressed as:

(5)$\mathit{y}\left(t\right)=\mathit{A}\left(\theta \right)\mathit{s}\left(t\right)+\mathit{n}\left(t\right)$

where

(6)$\mathit{s}\left(t\right)={\left[S_1\left(t\right),S_2\left(t\right),\dots ,S_K\left(t\right)\right]}^T$

is the incident signal,

(7)$\mathit{n}\left(t\right)={\left[n_1\left(t\right),n_2\left(t\right),\dots ,n_K\left(t\right)\right]}^T$

is the noise signal,

(8)$\mathit{A}\left(\theta \right)={\left[\mathit{a}\left({\theta }_1\right),\mathit{a}\left({\theta }_2\right),\dots ,\mathit{a}\left({\theta }_K\right)\right]}^T$

is the flow matrix of the array, which is determined by the geometric structure of the array and the incident angle together, and

(9)$\mathit{a}\left(\theta \right)={\left[1,e^{-j(2\pi d\sin \theta /\lambda )},\dots ,e^{-j(2\pi d\sin \theta /\lambda )}\right]}^T$

is the steering vector of the array.

Practical array signal processing is basically a multi-snapshot situation, and the expression of the array received signal can be written as follows:

(10)$\mathit{Y}=\mathit{A}\left(\theta \right)\mathit{S}+\mathit{N}$

where

(11)$\mathit{Y}={\left[\mathit{y}\left(t_1\right),\mathit{y}\left(t_2\right),\dots ,\mathit{y}\left(t_T\right)\right]}^T$

is the array output signal,

(12)$\mathit{S}={\left[\mathit{s}\left(t_1\right),\mathit{s}\left(t_2\right),\dots ,\mathit{s}\left(t_T\right)\right]}^T$

is the incident signal matrix, and

(13)$\mathit{N}={\left[\mathit{n}\left(t_1\right),\mathit{n}\left(t_2\right),\dots ,\mathit{n}\left(t_T\right)\right]}^T$

is the noise signal matrix. $T$ is snapshots.

There, $\mathrm{A}\in {ℂ}^{M\times N}$ is the observation matrix $K<M\ll N$, that is, the sparsity of the signal to be recovered is much less than the size of the observation matrix, so the final received signal can be linearly represented by the following equation:

(14)$Y=\left[y\left(1\right),\dots ,y\left(k\right)\right]=A(\theta )X$

This is a norm-based solution method. When some measurement vectors share the same support set, we will recover a joint sparse signal from multiple measurement vectors.

(15)$\min \parallel X{\parallel }_{g,h}\mathrm{subject} \mathrm{to}Y=AX$

We can obtain (15), where $\parallel X{\parallel }_{g,h}$ is defined as follows:

(16)$\parallel X{\parallel }_{g,h}={\left[{{\displaystyle \sum }}_{\tilde{j}=1}^N{\left({{\displaystyle \sum }}_{\tilde{i}=1}^k{\left|X\left[\tilde{j},\tilde{i}\right]\right|}^g\right)}^{h/g}\right]}^{1/h}$

When $g=2$ and $h=0$ are equal, it is the $l_{2,0}$ norm, and when $g=2$ and $h=1$ are equal, it is the $l_{2,1}$ norm. Here arises a question, whether multi-snapshot data can use the idea of column blocking by restoring each shot’s original signal to calculate the joint sparse signal. The following example can fully negate this idea.

Assuming the observation matrix is

(17)$A=\left[\begin{array}{ccccc}2& 0& 0& 1& 0\\ 0& 0.5& 0& 1& 0\\ 0& 0& 1& 2& -0.5\\ 0& 0& 0& -1& 0.5\end{array}\right]$

and the signal received by the signal processor in the DOA estimation problem is

(18)$B={\left[\begin{array}{c}1 1 0 0 \\ 1 1 1 0 \end{array}\right]}^T$

The observation matrix partitions the original signal into spatial angle divisions, ranging from four to five divisions. This division process facilitates the recovery of the initial input signal $X$ sparsely through the equation $AX=B$. The matrix $B$ is subsequently column-blocked into two column vectors, denoted as $B=[b_1,b_2]$.

To elaborate, the individual components of $B$ are organized into two distinct vectors. Specifically, the two single-shot input signals, represented by $x_1={[0.5 2 0 0 0]}^T$ and $x_2={[0 0 1 2 0]}^T$, are determined by utilizing a sparse recovery method applied to the respective equations $A{x}_1=b_1$ and $A{x}_2=b_2$.

Subsequently, the $l_{2,0}$ norm of the concatenated vector $[x_1,x_2]$ is computed, resulting in a value of 4 (indicating the number of non-zero row vectors after the matrix is block diagonalized). However, it is noteworthy that there exists a more parsimonious solution with fewer sparse degrees. By denoting the five angles in space as $[1°,2°,3°,4°,5°]$, the estimated Direction of Arrival (DOA) of the target obtained through two single-shot data processing methods is denoted as single-shot $[1°,2°,3°,4°]$ and $[1°,2°,3°]$, while the estimated DOA obtained through joint sparse recovery is represented as joint recovery.

Upon analysis, it is evident that joint recovery provides a more accurate estimation compared to single-shot [25]. Consequently, this example highlights that relying solely on multiple single-shot data processing methods may not guarantee the desired solution. The specific issue demonstrated in this example could potentially arise in real Direction of Arrival (DOA) estimation problems. Hence, determining effective strategies for leveraging multiple fast-shot data becomes a crucial and nuanced aspect of addressing such challenges.

4. Proof of Model Equivalence and Proposed Algorithm

It can be seen that an underdetermined equation $AX=B$ has infinite solutions, but for DOA estimation in this paper, the sparsest solution is always wanted, and there are many methods to find this sparsest solution. This paper mainly studies the sparse solution of the equation under mixed norm constraints.

Since the five classifications based on the minimum norm model analyzed before each have advantages and disadvantages, the $l_{2,th}$ norm model is analyzed here. Inspired by $l_0$ minimization and $l_p$ minimization model transformation, this is a method for finding optimal solutions to MMV problems, which solves the following minimum optimization problems

(19)$\underset{X\in {ℝ}^{n\times r}}{\min }\parallel X{\parallel }_{2,th} \mathrm{subject} \mathrm{to} AX=B$

The $l_0$ minimization problem is a problem that has been shown to be NP-hard due to discreteness, so the $l_{2,0}$ minimization problem is also NP-hard. However, $\parallel X{\parallel }_{2,0}={\lim }_{\alpha \to \infty }\parallel X{\parallel }_{2,t{h}_{\alpha }}$, the $l_{2,th}$ minimization problem is theoretically solvable, and the proof of this problem has been reflected in the literature [26]. The question of whether the solutions obtained by $l_{2,th}$ minimizing and $l_{2,0}$ minimizing the two models are equivalent will be discussed next.

The $l_0$ norm of a given vector is defined as

(20)${‖x‖}_0=\left\{\begin{array}{l}1⟨x,x⟩\ne 0\\ 0⟨x,x⟩=0\end{array}\right.$

Obviously, this piecewise function is not differentiable, so this paper approximates the $l_0$ norm with the following functions, where $\alpha >0$.

(21)${‖x‖}_0\approx {\displaystyle \sum _{i=1}^n\tanh \alpha x_i}={\displaystyle \sum _{i=1}^n\frac{e^{\alpha x_i}-e^{-\alpha x_i}}{e^{\alpha x_i}+e^{-\alpha x_i}}}$

The graph of the hyperbolic tangent function is shown in the Figure 3, and it can be seen that the larger the value $\alpha$, the closer the function of the $l_0$ norm.

The $l_2$ norm of a vector is defined as

(22)${‖x‖}_2=\sqrt{{\displaystyle \sum _{i=1}^n{\left(x_i\right)}^2}}$

The mixed $l_{2,0}$ norm of a matrix is defined as

(23)${‖X‖}_{2,0}={‖\left[{‖X_{row1}‖}_2;{‖X_{row2}‖}_2\dots ;{‖X_{row\mathrm{m}}‖}_2\right]‖}_0$

Bring in the $l_0$ norm model approximated by the hyperbolic tangent function and define it as

(24)${‖X‖}_{2,0}\approx {\displaystyle \sum _{\mathrm{i}=1}^m\tanh \alpha {‖x_{rowi}‖}_2}={‖X‖}_{2,th}$

Some ${‖X‖}_{2,th}$ properties discussed below illustrate their equivalence with the $l_{2,0}$ norm model.

Theorem 1.

For a given matrix $X\in {ℝ}^{m\times n}$, ${‖X‖}_{2,th}$ there is an upper bound $M$ and satisfied $M\le m$, there is a lower bound $0$.

A lemma needs to be proven before we can prove this theorem.

Lemma 1.

$f(\alpha x)={\displaystyle \frac{e^{\alpha x_i}-e^{-\alpha x_i}}{e^{\alpha x_i}+e^{-\alpha x_i}}} ,\alpha >0$ is a monotonically increasing function under $x\in [0,+\infty ]$ and its value is given by $[0,1)$.

Without loss of generality, consider a function $f(x)={\displaystyle \frac{e^x-e^{-x}}{e^x+e^{-x}}}$ of one variable, where $[-\infty ,+\infty ]$ is a derivable function, whose derivative is,

(25)$f^{\prime }(x)=1-f^2(x)$

Since $\underset{x\to +\infty }{{\displaystyle \lim }}f(x)=1$ and $\max \left|f(x)\right|=1$, it follows that $f^{\prime }(x)\ge 0$; therefore, $f(x)$ is a monotone nondecreasing function with upper bound 1, and according to $\underset{x\to 0}{{\displaystyle \lim }}f(x)=0$, $0\le f(x)\le 1x\in [0,+\infty ]$.

And then according to ${‖X‖}_{2,th}={\displaystyle \sum _{\mathrm{i}=1}^m\tanh {‖x_{rowi}‖}_2}$, considering ${‖x_{rowi}‖}_2$ is the $l_2$ norm of the row vector of a matrix $X$, and by the definition of $l_2$ norm, therefore ${‖x_{rowi}‖}_2$ is a bounded nonnegative number. And then according to the property of $f(x)$, for any $x_{rowi}$, the following equations are

(26)$0\le \tanh {‖x_{rowi}‖}_2\le 1$

Since $X\in {ℝ}^{m\times n}$, when all ${‖x_{rowi}‖}_2\to +\infty$, ${‖X‖}_{2,th}\to m$, when all ${‖x_{rowi}‖}_2=0$, ${‖X‖}_{2,th}=0$. There is an ${‖X‖}_{2,th}$ upper bound $M$ and satisfies $M\le m$, and there is a lower bound 0. Proof done.

Theorem 2.

For a given matrix $X\in {ℝ}^{m\times n}$, ${‖X‖}_{2,th}$ is a continuously derivable function for every element in the matrix.

To prove Theorem 2, we need to assume a time-varying matrix and, without loss of generality, let the matrix have only one variable, so we can assume that the matrix is as follows,

(27)$X=\left[\begin{array}{cccc}{x}_{11}(u)& x_{12}& \dots & x_{1n}\\ x_{21}& \dots & \dots & \dots \\ \dots & \dots & \dots & \dots \\ x_{m1}& \dots & \dots & x_{mn}\end{array}\right]$

And then compute it ${‖X‖}_{2,th}$,

(28)${‖X‖}_{2,th}=\tanh (\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}})+{\displaystyle \sum _{i=2}^m\tanh (\sqrt{{\displaystyle {\sum }_{k=1}^n{x}_{ik}^2}})}$

Since ${(•)}^n$ and $\tanh (•)$ are both continuous functions, according to Lemma 2 below,

Lemma 2.

A composition of continuous functions is a continuous function in its domain.

Therefore, ${‖X‖}_{2,th}$ is a continuous function, and according to the derivative formula of the composite function, its derivative is

(29)${\displaystyle \frac{\partial {‖X‖}_{2,th}}{\partial x_{11}(u)}}=\left[1-{\tanh }^2(\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}})\right]\left[{\displaystyle \frac{x_{11}(u)}{\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}}}}\right]$

Proof is completed.

Theorem 3.

For a given matrix $X\in {ℝ}^{m\times n}$, when ${‖X‖}_{2,th} x_{ij}\in [0,+\infty )$ is a monotone nondecreasing function for every element of the matrix.

To prove that Theorem 3 holds, we need to show that the ${‖X‖}_{2,th}$ derivative of Equation (29) exists ${\displaystyle \frac{\partial {‖X‖}_{2,th}}{\partial x_{11}(u)}}\ge 0$ for any $x_{11}(u)\in [0,+\infty )$. Analyze the second term in Equation (29), when $x_{11}(u)\in (0,+\infty )$ $\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}}>0$, and when $x_{11}(u)=0$ ${\displaystyle \frac{x_{11}(u)}{\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}}}}=0$, analyze the first term in Equation (29) by the properties of Equation (25) in Lemma 1; we know that when $x_{11}(u)\in [0,+\infty )$, $1-{\tanh }^2\left(\sqrt{{(x_{11}(u))}^2+{\displaystyle {\sum }_{i=2}^n{x}_{1i}^2}}\right)\ge 0$.

To sum up, when ${\displaystyle \frac{\partial {‖X‖}_{2,th}}{\partial x_{11}(u)}}\ge 0$ holds if $x_{11}(u)\in [0,+\infty )$ and only if $x_{11}(u)=0$, the equality sign holds. So, for the given matrix $X\in {ℝ}^{m\times n}$, when ${‖X‖}_{2,th} x_{ij}\in [0,+\infty )$, it is a monotone nondecreasing function for every element of the matrix. Proof done.

According to the proofs of Theorems 1–3, it is obvious that the $l_{2,th}$ model proposed in this paper is equivalent to the $l_{2,0}$ model.

Now, in order to design an algorithm to estimate the joint sparse matrix $X$, the following two theorems are introduced and proven to show the rationality of the sparse reconstruction of the matrix by solving the $l_{2,th}$ norm minimization.

Theorem 4.

For the matrix $Y\in {ℝ}^{m\times p}$ under the constraint $Y=AX$, it satisfies:

(30)$X=\Lambda \left(X\right)A^T{\left(A\Lambda \left(X\right)A^T\right)}^{†}Y$

where $\Lambda \left(X\right)\in {ℝ}^{m\times m}$ is a diagonal matrix whose entries on the diagonal $\Lambda {\left(X\right)}_{i,i}$ are:

(31)$\Lambda {\left(X\right)}_{i,i}={\displaystyle \frac{{‖X_{rowi}‖}_2}{{\tanh }^{\prime }\left({‖X_{rowi}‖}_2\right)}}$

The proof of Theorem 4 reconsiders the following questions:

(32)$\underset{X}{\min }\parallel X{\parallel }_{2,th}\mathrm{subject} \mathrm{to} Y=AX$

$\parallel X{\parallel }_{2,th}$ derivability with respect to each variable has already been proven in Theorem 2, so when using the KKT condition on Equation (32), a Lagrange function is defined here:

(33)$L\left(X,\lambda \right)=\parallel X{\parallel }_{2,th}-\lambda \left(AX-Y\right)$

It can be seen that the estimation algorithm is a form of fixed point iteration, so it is necessary to discuss the convergence of the algorithm. Namely, Theorem 5 below.

Theorem 5.

For a sequence of matrices $\left\{X_t\right\}$ satisfying the following recurrence relation:

(34) $X_{t+1}=\Lambda \left(X\right)A^T{\left(A\Lambda \left(X_t\right)A^T\right)}^{†}Y$

Satisfies the following inequality:

(35)$\parallel X_{t+1}{\parallel }_{2,th}\le \parallel X_t{\parallel }_{2,th}$

And $\underset{t\to \infty }{\lim }{X}_t=X_{end}$; this limit value is the final estimated value of Equation (34).

Now, Theorem 5 is proven. According to the Formula (30), the following conditions can be obtained:

(36)$\left\{\begin{array}{l}{\Lambda }^{†}\left(X\right)X_{t+1}=A^T{\left(A\Lambda \left(X_t\right)A^T\right)}^{†}Y\\ Y=A{X}_{t+1}\end{array}\right.$

The Formula (36) can be transformed into the following minimum quadratic optimization problem:

(37)$\min X^T{\Lambda }^{†}\left(X_t\right)X \mathrm{subject} \mathrm{to} Y=AX$

Namely,

(38)${\displaystyle \sum _{i=1}^m\frac{\tanh \left({‖X_{t,rowi}‖}_2\right)}{{‖X_{t,rowi}‖}_2}}{\left({‖X_{t+1,rowi}‖}_2\right)}^2\le {\displaystyle \sum _{i=1}^m\tanh \left({‖X_{t,rowi}‖}_2\right)}{‖X_{t,rowi}‖}_2$

From Theorem 3 obtained in the previous proof, we know that ${‖X‖}_{2,th}$ in $x_{ij}\in [0,+\infty )$ is a monotone nondecreasing function, so we can obtain:

(39)$\begin{array}{l}\tanh \left({‖X_{t+1,rowi}‖}_2\right)-{\displaystyle \frac{{‖X_{t+1,rowi}‖}_2^2\tan h^{\prime }\left(X_t\right)}{2{‖X_{t,rowi}‖}_2}}\\ \le \tanh \left({‖X_{t,rowi}‖}_2\right)-{\displaystyle \frac{{‖X_{t,rowi}‖}_2\tan h^{\prime }\left(X_t\right)}{2}}\end{array}$

According to Theorem 1, the iteration result has a lower bound 0, so the Formula (38) is proven and Theorem 5 holds. In this paper, a DOA estimation method based on $l_{2,th}$ minimization is proposed as Algorithm 1.

Algorithm 1. Second mixed norm constraint

Input: observed signal $Y$, measurement matrix $A$, signal sparsity $K$.Output: An estimate of the signal $X^{*}$.Initialization: Number of iterations $t=1$, let $M_1=I_n$.In the $t$th iteration: Step 1: Update the estimate of the original signal $X_{t+1}^{*}=M_t^{-1}{A}^T{(A{M}_t^{-1}{A}^T)}^{-1}Y$.Step 2: Use the ${\displaystyle \frac{{‖X_{rowi}‖}_2}{{\tanh }^{\prime }\left({‖X_{rowi}‖}_2\right)}} i=1,2,\dots ,n$ updated $M_{t+1}^{-1}$ main diagonal elements Step 3: Stop iterating on number of $t>K$ iterations or convergence

5. Numerical Experiment

The radiation source employed in this study is based on the digital TV signal DVB-T standard, characterized by a carrier frequency of 522 MHz and a signal bandwidth of 100 kHz. Figure 4 depicts a simulation of the corresponding baseband signal.

Furthermore, in Figure 5, the antenna pattern function of a single array element is illustrated, showcasing a consistent gain across the angular space of $[1°,90°]$. Where the blue line is the normalized gain of the antenna on all directions. This observation emphasizes the uniformity in gain within this specific angle range.

The experimental scenario in the simulation experiment unfolds in a two-dimensional plane space spanning an area of $2500\times 1000 {\mathrm{km}}^2$. The simulation parameters are delineated as Table 1:

In the first simulation experiment as Figure 6, the spatial angle is discretized into 891 grid points within the range of $[1°,90°]$. The anticipated location of the target is set at $30°$. The experimental findings reveal that, under a Signal-to-Noise Ratio (SNR) of 10 dB, all three Direction of Arrival (DOA) estimation algorithms consistently and accurately estimate the DOA of the target.

In the second simulation experiment as Figure 7, the number of targets in the spatial domain was increased, and the assumed Directions of Arrival (DOA) for the target were set at $[10°,20°,30°]$. They are displayed as the three colored circles in the figure. All other conditions remaining are unchanged. The obtained result figures indicate that all three methods accurately estimate the DOA of the targets.

In the actual passive radar DOA estimation, due to the limitation of hardware conditions, it is often impossible to have a lot of array elements, so the performance of the proposed DOA estimation method is particularly important when the array elements are insufficient.

In the simulation experiment in Figure 8, the number of elements is set to 8, and it can be seen that the $l_{2,0}$ M-OMP method cannot correctly estimate the target DOA at all, and the $l_{2,1}$ base-matching tracking method can estimate the target DOA, but the main lobe of the spectrum is too wide, while the method proposed $l_{2,th}$ in this paper can accurately estimate the target DOA.

In the simulation experiment in Figure 9, the number of elements is further reduced, and the number of elements is set to 4. It can be seen that the M-OMP method still fails to correctly estimate the target DOA. Although the base matching tracking method can estimate the target DOA, the error of DOA estimation will be large as the main lobe of the spectrum becomes wider. However, the proposed method can still accurately estimate the target DOA.

DOA estimation methods based on sparse recovery generally require prior information about the number of targets. However, in the actual passive radar system, this prior information cannot be obtained most of the time, and the sparse recovery is generally carried out using a fixed number of iterations. Therefore, whether the algorithm can stably estimate the DOA of the target with the increase in the number of iterations is particularly important.

Figure 10 shows the result of setting the number of iterations to 60. It can be seen that the M-OMP method cannot correctly estimate the target DOA at all, while the base-matching tracking method can estimate the target DOA but has some secondary peaks, while the method proposed in this paper can accurately estimate the target DOA.

We want to apply the algorithms proposed in the theory to practical problems, so it is necessary to conduct simulation experiments that approximate the real environment. We simulate a real passive radar scenario by setting the number of array elements to 8, the number of sources to 1, and the signal-to-noise ratio to 0 dB. Considering that phase error often exists in real equipped receiving systems, this phase error may be due to the spacing of antenna array elements presenting an unknown non-uniformity in the production or use, or it may also be caused by the amplitude of different channels of the hardware for the conversion of analog to digital signals in the receivers; phase inconsistencies result. Therefore, the next simulation experiments consider the effect of random phase errors with a 60° Gaussian distribution for each array element on the DOA estimation.

Figure 11 shows the simulation results of the proposed method and the two methods compared in the relatively realistic scenario in the above. We set the iteration termination condition to 15 times and all three algorithms have outputs at the real DOA position of the target. The M-OMP results have many false peaks with an amplitude of more than 0.5, which are difficult to identify due to the high accuracy requirements of the M-OMP method on the antenna and receiver. The basis pursuit results have relatively few false peaks with an amplitude of more than 0.5, but have too many burrs, which is not conducive to the actual DOA estimation. The proposed method can accurately estimate the DOA, but there are fluctuations around the true DOA of the target, which is generally better than the two compared methods and is more likely to be applied in practice.

This study conducts a comparative analysis of the root mean square error (RMSE) in Angle estimation for three sparse recovery algorithms. The evaluation is performed under various signal-to-noise ratios and differing sparsity and snapshot numbers. The RMSE is defined as follows:

(40)$E_{RMS}=\left({\displaystyle \frac{1}{M}}\sqrt{{\displaystyle \sum _{i=1}^M{\left(\stackrel{\wedge }{\theta }-{\theta }_{real}\right)}^2}}\right)$

In the context of the following experiments, $\stackrel{\wedge }{\theta }$ represents the outcome of Direction of Arrival (DOA) estimation, ${\theta }_{real}$ signifies the true angle of arrival, and $M$ stands for the number of Monte Carlo experiments, consistently set to 100 in all subsequent trials.

In the first comparative experiment as Figure 12, the RMSE curves for the Direction of Arrival (DOA) estimation of three targets are presented over a signal-to-noise ratio range of $[-30 \mathrm{dB},0 \mathrm{dB}]$. The overall RMSE of the three methods shows a decreasing trend with the signal-to-noise ratio. The accuracy of the M-OMP method improves rapidly as the signal-to-noise ratio increases. The basis pursuit method improves accuracy slowly. The figure in this paper illustrates that, as the Signal-to-Noise Ratio (SNR) increases, the mean square error of the proposed second-norm constraint algorithm surpasses that of the $l_{2,0}$ and $l_{2,1}$ minimization methods, particularly in low SNR scenarios. However, in high SNR situations, the performance of the proposed algorithm does not outperform that of the M-OMP algorithm. This can be attributed to the fact that the norm model is a discontinuous function and experiences rapid changes in ideal conditions.

In the second comparative experiment as Figure 13, we investigated the impact of the angular separation between two targets on the measurement error. Specifically, the reference angle of one target is fixed at $10°$, while the angle of the other target is set as $(10+\Delta \theta )°$. Simulation results reveal that the proposed algorithm exhibits superior resolution and yields the smallest Root Mean Square Error (RMSE) for the target positioned closer in angular space. It is worth mentioning that the basis pursuit method performs poorly when the target angle difference is small, but performs well when the target angle difference is large, and is suitable for DOA estimation when the targets are far apart in azimuthal angle.

In the third comparison experiment as Figure 14, we set up a change in the number of snapshots. Although the three DOA estimation methods are based on different algorithmic implementations, they may all be subject to the same constraints when these methods are compared at different numbers of snapshots. Under these constraints, the performance differences between the methods may be attenuated, leading to convergence of results. For example, if all the methods are constrained by the array element or the signal-to-noise ratio, their errors may fluctuate within a similar range. However, the proposed method fluctuates less and is relatively stable compared to the basis pursuit method. The good performance of the proposed method compared to the two methods compared at low snapshots is due to the nonlinear convex function.

In the fourth comparison experiment as Figure 15, we set up different numbers of array elements to analyze the RMSE of the three methods. Other conditions remain unchanged; the goniometric error usually tends to decrease with the increase in the number of array elements because more array elements can provide better spatial resolution and thus improve the accuracy of the measurement. Significant differences in goniometric errors may exist between different methods. M-OMP is more sensitive to the number of array elements and has a large error when the number of array elements is small, which decreases significantly after the number of array elements increases. This is related to its non-convexity. The basis pursuit method and the proposed method are not sensitive to the number of array elements, and overall, the proposed method performs better than the basis pursuit method, and the error of the proposed method decreases significantly when the number of array elements is higher than 12, which indicates that it is more usable in the actual low number of array element devices. This indicates that it is more usable in real devices with a low number of array elements. Different methods or algorithms may have significant differences in the goniometric errors under the same number of array elements, and the optimal method or algorithm should be selected according to the actual application scenarios and requirements.

As shown in Figure 16a, the antenna used in this field experiment is an eight-array antenna composed of two four-array antennas, with the main flap width recorded as [1°,90°], placed on the roof of the Hangzhou Research Institute of Xidian University, with the normal direction recorded as 45°, and the three locations of the Beigaofeng TV tower in Hangzhou, Hangzhou Research Institute of Xidian University, and Hangzhou Xiaoshan Airport are calibrated, as shown in Figure 16b. We set the normal direction of the antenna to the east direction (the direction of the black line in the figure) and use the 530 MHZ DVB-T signal transmitted from Hangzhou Beifeng TV tower to estimate the DOA of the echoes of the civil aviation targets taking off and landing at Hangzhou Xiaoshan Airport, under which the target can be regarded as a point target. As shown in Figure 17, the azimuth angle of the target that should appear according to the geographic location is 48°, while the azimuth angle measured by using the method proposed in this paper is 47.3°, which is 0.7° different from the real position.

6. Conclusions

To facilitate the judicious application of multi-snapshot data in the context of Direction of Arrival (DOA) estimation in passive radar, this paper undertakes a comprehensive exploration. Initially, it examines the correctness and efficacy of multi-snapshot data in contrast to single-snapshot data. Subsequently, acknowledging the non-convex nature of the norm model, an alternative model is introduced. The correctness and rationality of this alternative model are demonstrated, effectively addressing the NP-hard problem associated with norm optimization.

In conclusion, the proposed method’s rationality is substantiated through simulation experiments. Comparative analysis with classical sparse methods reveals that, under specific signal-to-noise ratio conditions, the proposed method exhibits smaller errors, is less sensitive to the number of snapshots, and demonstrates enhanced resolution in the DOA estimation of targets with similar spatial angles. The robustness of the proposed method in this paper is greatly improved compared to M-OMP and is more practical in real-life scenarios. We have completed experiments in real-life scenarios and achieved certain results.

In this research article, we designed a new algorithm for DOA estimation based on the proposed new model. Currently, within the more ideal case analysis, we believe that there is still a lot of room for improvement in the method. The current work is only a phase, and we will conduct more experiments in real environments so that the method can really be applied in practice.

Although this study tries to be rigorous and meticulous, due to the limitations of time, resources, and personal knowledge, there are inevitably omissions and shortcomings. It is hoped that this study will have the effect of drawing attention to these shortcomings.

Author Contributions

Software, J.H.; Conceptualization, C.W.; methodology, C.W.; Formal analysis, H.Z.; project administration, F.Z.; Simulation and verification, Z.Z.; Writing—original draft preparation, Z.Z.; Writing—review and editing, Z.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

The research in this paper does not involve ethical issues.

Informed Consent Statement

The research in this paper does not involve ethical issues.

Data Availability Statement

Due to some data privacy issues, the data used in this article cannot be provided.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Passive radar system model.

Figure 2. Uniform array antenna model.

Figure 3. Hyperbolic tangent functions with different parameters.

Figure 4. Digital TV signal simulation.

Figure 5. Single array antenna pattern.

Figure 6. Single target effect experiment.

Figure 7. Three-objective effect experiment.

Figure 8. Eight matrix DOA estimation results.

Figure 9. DOA estimation results of four array elements.

Figure 10. Results of 60 iterations.

Figure 11. Simulation analysis similar to the real environment.

Figure 12. Relationship between direction finding error and SNR.

Figure 13. Relationship between direction finding error and angular interval.

Figure 14. Algorithm performance versus number of snapshots.

Figure 15. Algorithm performance versus number of array elements.

Figure 16. Antennas and map used in the experiment.

Figure 17. Experimental results of the measured data.

Table 1

Experiment settings.

Radiation Source Location	(0, 0)
Receiving Station Location	(10, 0)
Target Position	(13.5, 50)
Number of Array Elements	64
Snapshot Count	1024
Spatial Meshing	$0.1°$

References

1. Lombardo, P.; Colone, F. Advanced processing methods for passive bistatic radar systems in Principles of Modern Radar. Adv. Radar Tech.; 2012; 43, pp. 739-821.

2. Palmer, J.; Cristallini, D.; Kuschel, H. Opportunities and current drivers for passive radar research In Proc. IEEE Radar Conf.; 2015; 14, pp. 145-150.

3. Klemm, R. Novel radar techniques and applications in Part III. Passiv. Multistatic Radar; 2017; 36, pp. 457-462.

4. Ulander, L.M.H.; Frolind, P.; Gustavsson, A.; Ragnarsson, R.; Stenstrom, G. VHF/UHF bistatic and passive SAR ground imaging. IEEE Radar Conf.; 2015; 6, pp. 669-673.

5. Dawidowicz, B.; Samczynski, P.; Malanowski, M.; Misiurewicz, J.; Kulpa, K.S. Detection of moving targets with multichannel airborne passive radar. IEEE Aerosp. Electron. Syst. Mag.; 2012; 27, pp. 42-49. [DOI: https://dx.doi.org/10.1109/MAES.2012.6380825]

6. Wojaczek, F.; Cristallini, D.J.; Schell, J.; Hagan, D.O. Polarimetric antenna diversity for improved reference signal estimation for airborne passive radar. IEEE Radar Conf.; 2020; 5, pp. 1-6.

7. Wan, L.T.; Kong, X.J.; Xia, F. Joint Range-Doppler-Angle Estimation for Intelligent Tracking of Moving Aerial Targets. IEEE Internet Things J.; 2018; 5, pp. 1625-1636. [DOI: https://dx.doi.org/10.1109/JIOT.2017.2787785]

8. Zeng, W.; He, H.Y.; Li, H.L. A dual-satellite passive direction finding and location method for aerial targets. J. Nanjing Univ. Aeronaut. Astronaut.; 2022; 54, pp. 234-248.

9. Han, K.; Huang, K.; Ko, S.W.; Lee, S.; Ko, W.S. Joint Frequency-and-Phase Modulation for Backscatter-Tag Assisted Vehicular Positioning. J. Stage; 2019; 1905, pp. 114-119.

10. Zhu, K.H.; Jiang, H.; Huo, Y.C.; Yu, Q.; Li, J.F. A Direct Position Determination Method Based on Subspace Orthogonality in Cross-Spectra under Multipath Environments. Sensors; 2022; 22, 7245. [DOI: https://dx.doi.org/10.3390/s22197245]

11. Donoho, D.L. Compressed sensing. IEEE Trans. Inf. Theory; 2006; 52, pp. 1289-1306. [DOI: https://dx.doi.org/10.1109/TIT.2006.871582]

12. Qaisar, S.; Bilal, R.M.; Iqbal, W.; Naureen, M.; Lee, S. Compressive sensing: From theory to applications, a survey. J. Commun. Netw.; 2013; 15, pp. 443-456. [DOI: https://dx.doi.org/10.1109/JCN.2013.000083]

13. Hu, L.; Shi, Z.G.; Zhou, J.X.; Fu, Q. Compressed sensing for complex sinusoids: And approach based on dictionary refinement. IEEE Trans. Signal Process.; 2012; 60, pp. 3809-3822.

14. Lin, J.; Ma, X.; Yan, S.; Hao, C.; Lin, G. Off-grid fast relevance vector machine algorithm for direction of arrival estimation. IET Radar Sonar Navig.; 2016; 10, pp. 718-725. [DOI: https://dx.doi.org/10.1049/iet-rsn.2015.0304]

15. Fang, Y.; Wu, J.; Huang, B. 2D sparse signal recovery via 2D orthogonal matching pursuit. Sci. China Inf. Sci.; 2012; 55, pp. 889-897. [DOI: https://dx.doi.org/10.1007/s11432-012-4551-5]

16. Zhang, R.; Wang, S.; Lu, X.; Duan, W.; Cai, L. Two-Dimensional DoA Estimation for Multipath Propagation Characterization Using the Array Response of PN-Sequences. IEEE Trans. Wirel. Commun.; 2016; 15, pp. 341-356. [DOI: https://dx.doi.org/10.1109/TWC.2015.2473156]

17. Wang, L.; Zhao, L.; Bi, G.; Wan, C.; Zhang, L.; Zhang, H. Novel Wideband DOA Estimation Based on Sparse Bayesian Learning with Dirichlet Process Priors. IEEE Trans. Signal Process.; 2016; 64, pp. 275-289. [DOI: https://dx.doi.org/10.1109/TSP.2015.2481790]

18. Wang, Q.S.; Yu, H.; Li, J.; Dong, C.; Ji, F.; Chen, Y.K. Spatial proximate Signal DOA estimation algorithm based on sparse Bayesian learning. J. Electron. Inf. Technol.; 2022; 43, pp. 1145-1153.

19. Hu, B.; Shen, X.Y.; Jiang, M.; Zhao, K.H. Off-Grid DOA Estimation Based on Compressed Sensing on Multipath Environment. Int. J. Antennas Propag.; 2023; 10, pp. 1246-1249. [DOI: https://dx.doi.org/10.1155/2023/9315869]

20. Liu, L.X.; Zhang, X.; Chen, P. Compressed Sensing-Based DOA Estimation with Antenna Phase Errors. Electronics; 2019; 8, 294. [DOI: https://dx.doi.org/10.3390/electronics8030294]

21. Gungor, A.; Kilic, B. Hardware-Efficient Direction of Arrival Estimation using Compressive Sensing. IEEE Int. Symp. Phased Array Syst. Technol.; 2022; 6, pp. 1123-1134.

22. Blasone, G.P.; Colone, F.; Lombardo, P.; Wojaczek, P.; Cristallini, D. Passive radar DPCA schemes with adaptive channel calibration. IEEE Trans. Aerosp. Electron. Syst.; 2020; 56, pp. 4014-4034. [DOI: https://dx.doi.org/10.1109/TAES.2020.2987478]

23. Meng, X.T.; Cao, B.X.; Yan, F.G.; Greco, M.; Gini, F.; Zhang, Y. Real-valued MUSIC for efficient direction of arrival estimation with arbitrary arrays: Mirror suppression and resolution improvement. Signal Process.; 2023; 202, 108766. [DOI: https://dx.doi.org/10.1016/j.sigpro.2022.108766]

24. Liu, Y.; Wan, X.; Tang, H.; Yi, J.; Cheng, Y.; Zhang, X. Digital television based passive bistatic radar system for drone detection. IEEE Radar Conf.; 2017; 51, pp. 1493-1497.

25. Tao, R.; Gao, Z.W.; Wang, Y. Side Peaks Interference Suppression in DVB-T Based Passive radar. IEEE Trans. Aerosp. Electron. Syst.; 2012; 48, pp. 3610-3619. [DOI: https://dx.doi.org/10.1109/TAES.2012.6324746]

26. Wang, C.L.; Peng, J.G. Exact Recovery of Sparse Multiple Measurement Vectors by l_2,p-minimization. J. Inequalities Appl.; 2018; 17, pp. 45-60. [DOI: https://dx.doi.org/10.1186/s13660-017-1601-y]