Fast Implementation of Approximated Maximum Likelihood Parameter Estimation for Frequency Agile Radar under Jamming Environment

Abstract

A computationally efficient target parameter estimation algorithm for frequency agile radar (FAR) under jamming environment is developed. First, the barrage noise jamming and the deceptive jamming are suppressed by using adaptive beamforming and frequency agility. Second, the analytical solution of the parameter estimation is obtained by a low-order approximation to the multi-dimensional maximum likelihood (ML) function. Due to that, fine grid-search (FGS) is avoided and the computational complexity is greatly reduced.

1. Introduction

With the development of new countermeasure technology, radar is faced with an increasingly complex electromagnetic environment [1]. Therefore, how to eliminate the negative influences using the complex environment is the precondition of successful target identification [2]. It is well known that the frequency agile radar (FAR) that adopts the pulse-to-pulse random frequency agility over shorter intervals can improve the performance in target detection [3,4], and its thumbtack ambiguity function [5] implies that one is able to obtain high range–velocity resolution and avoid range–velocity coupling problems which are present in traditional frequency agility waveforms [6]. Thus, the FAR is an important anti-electronic countermeasure [7,8].

As for the target detection, the early studies were mainly focused on incoherent integration such as the Hough transform [9] and the Radon transform based methods [10]. However, it has integration loss compared with the coherent method. Consequently, some researchers proposed the target detection method based on coherent integration [11]. In [12,13], Zhou et al. proposed an adaptive monopulse (AM) algorithm based on combining sum–difference and auxiliary beam at subarray level for anti-jamming. In [14], based on a wideband model, Qian et al. proposed a wideband scale-based Radon–Fourier transform (WSRFT) method for coherent integration of high-speed targets, but only for fixed carrier frequency radars. In [7], Wang et al. proposed an improved frequency agile coherent radon transform (FA-CRT) scheme to cope with the long-time coherent integration in frequency agile radar. In [15], P. Huang et al. proposed a method for removing the residual coupling between the agile carrier frequency and slow-time by using Keystone Transform (KT) between different bursts to achieve coherent accumulation. In [16], Li et al. proposed a combined range–velocity estimation method based on compressive sensing (CS) in the range–velocity plane, where the CS is employed to reduce the number of channels, therefore, the system complexity and the computational burden are effectively reduced [17].

However, the computational complexities of CRT and CS are very high which greatly limits their applications. In this communication, a computationally efficient algorithm for estimating target parameters of FAR under jamming environment is developed. The barrage noise jamming and the deceptive jamming are firstly suppressed by using adaptive beamforming and frequency agility [18], and then a low-order approximation of the multi-dimensional maximum likelihood (ML) function is employed to estimate target parameters [19]. Due to the fact that the fine grid-search is not necessary anymore, the computational complexity is greatly reduced relative to the conventional grid-search method. The experimental results show that this method has a lower computational complexity at the cost of trivial performance loss.

2. System Model Characterization

2.1. Signal Model

Assume that there are M elements and FAR transmits N pulses with random carrier frequencies. The frequency of the $n\mathrm{th}$ pulse can be denoted by $f_n=f_0+k_n\mathsf{\Delta }f$, where $f_0$ is the initial carrier frequency, $\mathsf{\Delta }f$ is the minimal frequency agility interval, and $k_n$ controls the rules of how carrier frequency varies. Based on the signal model given in [15], the received signal of the $n\mathrm{th}$ pulse at the $m\mathrm{th}$ array element after matched filtering can be expressed as

$$s_{mn}\left(t,\tau \right)=a\sin c\left(\mathsf{\Delta }f\left(t-{\tau }_{mn}\right)\right)\exp \left(-j2\pi f_n{\tau }_{mn}\right)$$

(1)

where $a$ is the complex amplitude of a moving point-like target, and ${\tau }_{mn}$ is the target time delay of the $n\mathrm{th}$ pulse at the $m\mathrm{th}$ array element, which can be written as

$${\tau }_{mn}=2\left(r-vnT-\frac{1}{2}md\sin \theta \right)/c$$

(2)

where $r$ is the target range, $v$ is the target velocity, $d$ is the array element spacing and $\theta$ is the target angle. For simplicity, we define $\mathsf{\eta }={\left[\begin{array}{ccc}r& v& \theta \end{array}\right]}^{\mathrm{T}}$ as the target parameter vector. The target data can be assembled in the form of a space–time snapshot

$$\overline{s}\left(\mathsf{\eta }\right)=as\left(\mathsf{\eta }\right)=a{\left[\begin{array}{cccc}{s}_0^{\mathrm{T}}& s_1^{\mathrm{T}}& \dots & s_{N-1}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$$

(3)

where $s_n={\left[\begin{array}{cccc}\exp \left(-j2\pi f_n{\tau }_{0n}\right)& \exp \left(-j2\pi f_n{\tau }_{1n}\right)& \dots & \exp \left(-j2\pi f_n{\tau }_{\left(M-1\right)n}\right)\end{array}\right]}^{\mathrm{T}}$ is an $M\times 1$ spatial steering vector of the $n\mathrm{th}$ pulse, the superscript T means the matrix transpose. Meanwhile, the deceptive jamming, which is usually obtained by delaying the transmitting signal for multiple pulse repetition intervals, can be easily suppressed by matched filtering. Thus, the component of deceptive jamming is omitted after matched filtering. Furthermore, the barrage noise jamming often emerges in various radar scenes and should be taken into account here. According to the model given in [20], the jammer space–time snapshot may be written as $x_J={\alpha }_J\otimes a_J$ and the jammer space–time covariance matrix is $R_J=E\left(x_J{x}_J^H\right)=I_N\otimes {\Phi }_J$, where ${\alpha }_J$ is a random vector containing the jammer amplitudes, $a_J$ is the jamming steering vector, and ${\Phi }_J={\sigma }_J^2{a}_J{a}_J^H$ means the jamming spatial covariance matrix, and ${\sigma }_J^2$ is the jamming power. In addition, the noise space–time snapshot is denoted by $x_{\mathrm{n}}$ and the noise space–time covariance matrix is $R_{\mathrm{n}}=E\left[x_{\mathrm{n}}{x}_{\mathrm{n}}^H\right]={\sigma }^2{I}_{MN}$, where ${\sigma }^2$ is the noise power. Finally, the echo data can be expressed as

$$x_1=as\left(\mathsf{\eta }\right)+x_0$$

(4)

where $x_0=x_n+x_J$ encompasses the barrage noise jamming and noise, which has the covariance matrix

$$R_0=E\left[x_0{x}_0^H\right]=R_{\mathrm{J}}+R_{\mathrm{n}}=I_N\otimes {\Phi }_{\mathrm{Jn}}$$

(5)

where ${\Phi }_{\mathrm{Jn}}={\Phi }_{\mathrm{J}}+{\sigma }^2{I}_M$.

2.2. ML Parameter Estimation

In this section, a computationally efficient algorithm for target parameters estimation based on the ML method is developed. The joint probability density function of the observed data is given by

$$f\left(x_1|\mathsf{\eta }\right)=\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left|R_0^{}\right|}^{\frac{1}{2}}}\exp \left(-{\left(x_1-as\left(\mathsf{\eta }\right)\right)}^H{R}_0^{-1}\left(x_1-as\left(\mathsf{\eta }\right)\right)\right)$$

(6)

According to the definition of the likelihood function and ignoring the constant term, the corresponding log-likelihood function is

$$L\left(a,\mathsf{\eta }\right)={\left(x_1-as\left(\mathsf{\eta }\right)\right)}^H{R}_0^{-1}\left(x_1-as\left(\mathsf{\eta }\right)\right)$$

(7)

The ML estimate of $a$ can be obtained by setting the derivative of the log-likelihood function with respect to $a$ equal to zero. After some manipulations, we have

$$\widehat{a}=\frac{s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}{s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}$$

(8)

Substituting (8) into (7) yields the ML estimate of $\mathsf{\eta }$, that is

$$\widehat{\mathsf{\eta }}=\arg \underset{\mathsf{\eta }}{\max }L\left(\mathsf{\eta }\right)=\frac{{\left|s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{x}_1\right|}^2}{s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}$$

(9)

In order to obtain the ML estimate of $\mathsf{\eta }$, it is necessary to perform a search over the whole parameter grids. However, the three-dimensional fine grid-search brings large computational cost. To reduce the computational complexity, we can resort to a solving equation. As we know, the ML solution should satisfy the following condition: the partial derivatives of the log-likelihood function with respect to these parameters are equal to zero, that is

$$\left[\begin{array}{l}{\dot{L}}_{\theta }\left(\mathsf{\eta }\right)\\ {\dot{L}}_v\left(\mathsf{\eta }\right)\\ {\dot{L}}_r\left(\mathsf{\eta }\right)\end{array}\right]=\left[\begin{array}{l}0\\ 0\\ 0\end{array}\right]$$

(10)

where ${\dot{L}}_{\theta }\left(\mathsf{\eta }\right)$, ${\dot{L}}_v\left(\mathsf{\eta }\right)$, ${\dot{L}}_r\left(\mathsf{\eta }\right)$ represents the first-order partial derivative of $L\left(\mathsf{\eta }\right)$ with respect to $\theta$, $v$ and $r$, respectively. After some manipulations, ${\dot{L}}_{\theta }\left(\mathsf{\eta }\right)$ can be written as

$${\dot{L}}_{\theta }\left(\mathsf{\eta }\right)=\frac{\left(\mathsf{\Delta }{\Sigma }^{\mathrm{H}}+{\mathsf{\Delta }}^{\mathrm{H}}\Sigma \right)c_{11}-\Sigma {\Sigma }^{\mathrm{H}}\left(c_{12}+c_{21}\right)}{c_{11}^2}=0$$

(11)

where $\{\begin{array}{l}\mathsf{\Delta }={\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1\\ \Sigma =s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{x}_1\\ c_{11}=s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{s}^{\mathrm{H}}\left(\mathsf{\eta }\right)\\ c_{12}=s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)\\ c_{21}={\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\end{array}$

where ${\dot{s}}_{\theta }^{}\left(\mathsf{\eta }\right)$ represents the first-order partial derivative of $s\left(\mathsf{\eta }\right)$ with respect to $\theta$.

Simplifying (11), ${\dot{L}}_{\theta }\left(\mathsf{\eta }\right)$ can be expressed as

$${\dot{L}}_{\theta }\left(\mathsf{\eta }\right)=\mathrm{Re}\left(\frac{{\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}\right)-\mathrm{Re}\left(\frac{{\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}\right)=0$$

(12)

Rearranging (12), the equation can be written as

$$\mathrm{Re}\left(\frac{w_{\theta }^{\mathrm{H}}\left(\mathsf{\eta }\right)x_1}{w_{}^{\mathrm{H}}\left(\mathsf{\eta }\right)x_1}\right)=\mathrm{Re}\left(\frac{w_{\theta }^{\mathrm{H}}\left(\mathsf{\eta }\right)s\left(\mathsf{\eta }\right)}{w_{}^{\mathrm{H}}\left(\mathsf{\eta }\right)s\left(\mathsf{\eta }\right)}\right)$$

(13)

where $w_{\theta }^{}\left(\mathsf{\eta }\right)=R_0^{-1}{\dot{s}}_{\theta }\left(\mathsf{\eta }\right)$ and $w\left(\mathsf{\eta }\right)=R_0^{-1}s\left(\mathsf{\eta }\right)$ can be considered as the parameter-dependent adaptive sum and difference weight vectors.

Since the equations in (13) are complex nonlinear functions, the solution cannot be solved in analytical form. Observing (13), we can see that the complexity of (13) mainly results from the fact that the sum and difference weight vectors are functions of $\mathsf{\eta }$. By searching the values of $\mathsf{\eta }$, the weight vectors can obtain more matching results with the echo data, and thus the optimal solution is guaranteed. Although (13) can obtain the optimal solution of parameter estimation, the computational complexity resulted from fine grid-search is unacceptable. In practice, the performance improvement due to the variation of $\mathsf{\eta }$ only in one resolution is trivial. Thus, we can approximate the ML solution by setting the parameter-independent adaptive sum and difference weight vectors in (13), that is, $w_{\theta }^{}\left({\mathsf{\eta }}_0\right)=R_0^{-1}{\dot{s}}_{\theta }\left({\mathsf{\eta }}_0\right)$ and $w\left({\mathsf{\eta }}_0\right)=R_0^{-1}s\left({\mathsf{\eta }}_0\right)$, then (13) can be simplified as

$$\mathrm{Re}\left(\frac{w_{\theta }^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}\right)=\mathrm{Re}\left(\frac{w_{\theta }^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}\right)$$

(14)

Similarly, ${\dot{L}}_v\left(\mathsf{\eta }\right)$ can be written as

$${\dot{L}}_v\left(\mathsf{\eta }\right)=\mathrm{Re}\left(\frac{{\dot{s}}_v^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}\right)-\mathrm{Re}\left(\frac{{\dot{s}}_v^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}\right)=0$$

(15)

Approximating (15) by using the parameter-independent adaptive sum and difference weight vectors, (15) can be approximated by

$$\mathrm{Re}\left(\frac{w_v^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}\right)=\mathrm{Re}\left(\frac{w_v^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}\right)$$

(16)

where $w_v^{}\left({\mathsf{\eta }}_0\right)=R_0^{-1}{\dot{s}}_v\left({\mathsf{\eta }}_0\right)$.

Meanwhile, ${\dot{L}}_r\left(\mathsf{\eta }\right)$ can be written as

$${\dot{L}}_r\left(\mathsf{\eta }\right)=\mathrm{Re}\left(\frac{{\dot{s}}_r^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}{x}_1}\right)-\mathrm{Re}\left(\frac{{\dot{s}}_r^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}{s_{}^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)}\right)=0$$

(17)

and (17) can be approximated by

$$\mathrm{Re}\left(\frac{w_r^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}\right)=\mathrm{Re}\left(\frac{w_r^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)}\right)$$

(18)

where $w_r^{}\left({\mathsf{\eta }}_0\right)=R_0^{-1}{\dot{s}}_r\left({\mathsf{\eta }}_0\right)$.

Finally, the ML solution can be obtained by solving Equations (14), (16) and (18). Define $F_i\left(\mathsf{\eta }\right)=\mathrm{Re}\left(w_i^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)/\left(w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)s\left(\mathsf{\eta }\right)\right)\right),i=\theta ,r,v$ as the response curve of target parameters and ${\mathsf{\eta }}_i=\mathrm{Re}\left(\frac{w_i^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}{w_{}^{\mathrm{H}}\left({\mathsf{\eta }}_0\right)x_1}\right),i=\theta ,r,v$ as the difference–sum–ratio output. Then (14), (16) and (18) can be written as

$$\left[\begin{array}{l}{F}_{\theta }\left(\mathsf{\eta }\right)\\ F_v\left(\mathsf{\eta }\right)\\ F_r\left(\mathsf{\eta }\right)\end{array}\right]=\left[\begin{array}{l}{\mathsf{\eta }}_{\theta }\\ {\mathsf{\eta }}_v\\ {\mathsf{\eta }}_r\end{array}\right]$$

(19)

To acquire the solution of (19), one-order Taylor approximation to $F_i\left(\mathsf{\eta }\right)$ can be employed

$$F_i\left(\mathsf{\eta }\right)=p_i+{\delta }_{i\theta }\left(\theta -{\theta }_0\right)+{\delta }_{iv}\left(v-v_0\right)+{\delta }_{ir}\left(r-r_0\right),i=\theta ,r,v$$

(20)

where $p_i=F_i\left(\mathsf{\eta }\right){|}_{\mathsf{\eta }={\mathsf{\eta }}_0}$, ${\delta }_{i\theta }=\partial F_i\left(\mathsf{\eta }\right)/\partial \theta {|}_{\mathsf{\eta }={\mathsf{\eta }}_0}$, ${\delta }_{iv}=\partial F_i\left(\mathsf{\eta }\right)/\partial v{|}_{\mathsf{\eta }={\mathsf{\eta }}_0}$, ${\delta }_{ir}=\partial F_i\left(\mathsf{\eta }\right)/\partial r{|}_{\mathsf{\eta }={\mathsf{\eta }}_0}$, and $\partial$ means the partial derivative operator. After some manipulations, we have

$\{\begin{array}{l}{\delta }_{i\theta }=\mathrm{Re}\left(\frac{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\ddot{s}}_{i\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+{\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_{\theta }^{}\left(\mathsf{\eta }\right)\right)-\left({\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\dot{s}}_{\theta }^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_{\theta }^{}\left(\mathsf{\eta }\right)\right)}{{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)}^2}\right)\\ {\delta }_{iv}=\mathrm{Re}\left(\frac{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\ddot{s}}_{iv}^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+{\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_v^{}\left(\mathsf{\eta }\right)\right)-\left({\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\dot{s}}_v^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_v^{}\left(\mathsf{\eta }\right)\right)}{{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)}^2}\right)\\ {\delta }_{ir}=\mathrm{Re}\left(\frac{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\ddot{s}}_{ir}^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+{\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_r^{}\left(\mathsf{\eta }\right)\right)-\left({\dot{s}}_i^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)\left({\dot{s}}_r^H\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)+s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}{\dot{s}}_r^{}\left(\mathsf{\eta }\right)\right)}{{\left(s^{\mathrm{H}}\left(\mathsf{\eta }\right)R_0^{-1}s\left(\mathsf{\eta }\right)\right)}^2}\right)\end{array}$ and ${\ddot{s}}_{ij}^H\left(\mathsf{\eta }\right)$ represents the second-order partial derivative of $s\left(\mathsf{\eta }\right)$ with respect to $i$ and $j$.

Reformulating (20) into a matrix form yields

$$\left[\begin{array}{l}{F}_{\theta }\left(\mathsf{\eta }\right)\\ F_v\left(\mathsf{\eta }\right)\\ F_r\left(\mathsf{\eta }\right)\end{array}\right]=\left[\begin{array}{l}{p}_{\theta }\\ p_v\\ p_r\end{array}\right]+\left[\begin{array}{lll}{\delta }_{\theta \theta }& {\delta }_{\theta v}& {\delta }_{\theta r}\\ {\delta }_{v\theta }& {\delta }_{vv}& {\delta }_{vr}\\ {\delta }_{r\theta }& {\delta }_{rv}& {\delta }_{rr}\end{array}\right]\cdot \left[\begin{array}{l}\theta -{\theta }_0\\ v-v_0\\ r-r_0\end{array}\right]$$

(21)

Substituting (21) into (19) yields

$$\left[\begin{array}{l}\theta \\ v\\ r\end{array}\right]=\left[\begin{array}{l}{\theta }_0\\ v_0\\ r_0\end{array}\right]+{\left[\begin{array}{lll}{\delta }_{\theta \theta }& {\delta }_{\theta v}& {\delta }_{\theta r}\\ {\delta }_{v\theta }& {\delta }_{vv}& {\delta }_{vr}\\ {\delta }_{r\theta }& {\delta }_{rv}& {\delta }_{rr}\end{array}\right]}^{-1}\cdot \left(\left[\begin{array}{l}{\mathsf{\eta }}_{\theta }\\ {\mathsf{\eta }}_v\\ {\mathsf{\eta }}_r\end{array}\right]-\left[\begin{array}{l}{p}_{\theta }\\ p_v\\ p_r\end{array}\right]\right)$$

(22)

To improve the accuracy of Taylor approximation, the iterative technique can be employed, which uses the estimated target parameters instead of the assumed target parameters, i.e., $[{\theta }_0,v_0,r_0]=[\widehat{\theta },\widehat{v},\widehat{r}]$. It is shown from (22) that target parameter estimation can be realized by a simple expression, where the computational complexity is $O\left(12{\left(NM\right)}^2+12NM+3^3\right)$ complex multiplications. Due to the fact that the computational complexity of complex addition is less than that of complex multiplication, complex additions are not considered here. In contrast, the conventional grid-search technique has the computational complexity of $O\left(\left({\left(NM\right)}^2+2NM\right)J^3\right)$ complex multiplications, where $J$ is the number of grids in one dimension. To guarantee the accuracy, the number of grids is usually larger than 100. To demonstrate the advantage of the computational complexity, the computational complexities of the presented method and the grid-search method as a function of the number of array elements are plotted in Figure 1, where the number of grids is 100, and the number of pulses is 192. It is observed that the presented method owns lower computational complexity than the fine grid-search method.

3. Numerical Experiments

In this section, the effectiveness of the proposed method is verified by simulation experiments. The radar parameters are designed as follows: a uniform linear array with 8 elements is employed, the transmitter transmits a multi-cycle inter-pulse stepped frequency signal by 192 pulses, each group has 16 pulses, and the carrier frequency hops from 2 to 2.016 GHz with the frequency increment of 1 MHz, the sampling frequency is 16 MHz, and the pulse repetition frequency is 4000 Hz.

3.1. Deceptive Interference Suppression and Target Parameter Estimation for Frequency Agile Radar

A frequency agile radar is designed based on the above parameters, and the synthesized bandwidth is increased by 16 times. Therefore, under the same parameters, the range resolution of the frequency agile radar is 16 times that of the fixed frequency radar. By using the bandwidth synthesis function of FAR [13], 16 fine range gates in each coarse range unit can be obtained. Thus, distance information can be estimated accurately.

To analyze the performance of suppressing deceptive jamming, the experiment randomly simulates multiple targets and multiple deceptive jamming located in various coarse range bins, velocities and distances. Figure 2 gives the echo plane after matched filtering for the frequency non-agile case and the frequency agile case, where multiple target signals and deceptive jamming signals with random parameters are simulated, where the targets are located in the range bins of 76, 131, 176 and 201, respectively, and the jamming signals are located in the range bins of 99, 126, 156 and 191, respectively. It is observed that the target signals and the deceptive jamming signals are all preserved for the frequency non-agile case in Figure 2a, whereas only the target signals are well preserved for the frequency agile case in Figure 2b.

In order to see the performance of suppressing deceptive jamming and range resolution enhancement resulted from frequency agility, three target signals are added in the echo data, where the two targets are set at the 71st coarse range bin and another one is located in the coarse range bin of 114. Meanwhile, four deceptive jamming signals are added, where two deceptive jamming signals are set at the 110th coarse range bin and others are located in the coarse range bin of 149. Additionally, the signals at the same range bin have different speeds and angles.

The matched filtering results for the frequency agile case and the frequency non-agile case are given in Figure 3 and Figure 4. Comparing the plane of coarse range bin and velocity given in Figure 3a and Figure 4a, it is clearly observed that the target signal is preserved, and the deceptive jamming is suppressed due to the agile carrier frequency. Figure 3b and Figure 4b are the range cuts of Figure 3a and Figure 4a at the velocity of −3.5 m/s, which show a perfect performance of deceptive jamming suppression for the frequency agile case. Furthermore, a plane of coarse range bin and fine range bin given in Figure 3c indicates that more accurate distance information can be obtained since there is a fine resolving within coarse range bin. In addition, Figure 3d shows the estimation of target speed and accurate distance.

3.2. Estimation of Target Parameters in Various Conditions

In this section, the simulation results in various conditions are carried out to validate the effectiveness of the proposed algorithm. Experiments of the fine grid-search method [21] and the adaptive monopulse method [12] are also performed to prove the superiority of the proposed method. All results are obtained by running a Monte Carlo simulation, where the number of Monte Carlo experiments is 500. For simplicity, the normalized root mean square error (NRMSE) is employed here, which is defined as ${\overline{\sigma }}_i={\sigma }_i/{\rho }_i,i=[\theta ,v,r]$, where ${\sigma }_i$ is the classical RMSE, and ${\rho }_i$ is the resolution. In this simulation, the resolution of ${\rho }_{\theta }$ is 14.32°, the resolution of ${\rho }_v$ is 1.56 m/s, and the resolution of ${\rho }_r$ is 9.37 m. The target and deceptive jamming are located at the range bins from 1500 to 3000 randomly, the velocity of them are from −6 to 6 m/s.

A. NRMSE in Signal Noise Ratio (SNR) Cases

In Figure 5, the NRMSEs as a function of SNR are given, where the SNR varies from 15 to 35 dB, the line of sight is 113° and the angle of the target is 113.5°. Parameter estimates of angle, range and velocity by three methods are also incorporated. It is shown that RMSEs of target parameters decrease with the SNR increasing, and the RMSEs of our presented method approach those of the fine grid-search method in most of SNR cases. Meanwhile, our presented method shows better estimated accuracy than the adaptive monopulse method.

B. NRMSE in the Number of Barrage Noise Jamming Cases

In Figure 6, the NRMSEs of angle, range and velocity as a function of the number of barrage noise jamming are given for the three methods, the line of sight is 113° and the angle of the target is 113.5°. It should be noted that as the number of barrage noise jamming varies from one to four, the RMSEs of target parameters will increase, which means the estimated accuracy is decreasing. Certainly, the target parameters estimated by the fine grid-search method is the closest to the actual values, followed by the ML method, and the adaptive monopulse method.

C. NRMSE in the Angle Deviation from the Line of Sight Cases

In Figure 7, the NRMSEs of the target angle versus angle deviation from the line of sight are depicted using three methods. The line of sight is 113°, where the deviation range between the target and the look direction is from 0° to 5°. It is observed that the RMSEs increase and the estimated performance decreases as the target angle deviates from the look direction. As shown in Figure 6, it can be clearly seen that the target parameters are finely estimated by the fine grid-search method and roughly estimated by the adaptive monopulse method, the ML method may achieve a good balance between them.

D. NRMSE in the Line of Sight Cases

In Figure 8, the NRMSEs as a function of the line of sight are given, where the look direction varies from 90° to 125° and the deviation between the target and the look direction is 0.5°, it is shown that as the look direction deviates from the normal direction, the RMSEs increase, and the parameter estimation performance decreases. Although the performance of the presented method is worse than that of the FGS method, it has a better performance compared to the adaptive monopulse method.

4. Conclusions

In this communication, fast approximated ML estimation for FAR under the sophisticated jamming environment including the deceptive jamming and the barrage noise jamming is considered. The presented method suppresses the barrage jamming by adaptive beamforming and the deceptive jamming by frequency agility. To reduce the computational complexity, we estimate the target parameters by a low-order approximation of the multi-dimensional ML function. Thus, the proposed method may achieve a good balance between the computational cost and the target parameter estimation performance. Numerical examples are given to demonstrate the effectiveness of the presented method.