A Novel Monopulse Technique for Adaptive Phased Array Radar

Abstract

The monopulse angle measuring technique is widely adopted in radar systems due to its simplicity and speed in accurately acquiring a target’s angle. However, in a spatial adaptive array, beam distortion, due to adaptive beamforming, can result in serious deterioration of monopulse performance. In this paper, a novel constrained monopulse angle measuring algorithm is proposed for spatial adaptive arrays. This algorithm maintains the ability to suppress the unwanted signals without suffering from beam distortion. Compared with conventional adaptive monopulse methods, the proposed algorithm adopts a new form of constraint in forming the difference beam with the merit that it is more robust in most practical situations. At the same time, it also exhibits the simplicity of one-dimension monopulse, helping to make this algorithm even more appealing to use in adaptive planar arrays. The theoretical mean and variance of the proposed monopulse estimator is derived for theoretical analysis. Mathematical simulations are formulated to demonstrate the effectiveness and advantages of the proposed algorithm. Both theoretical analysis and simulation results show that the proposed algorithm can outperform the conventional adaptive monopulse methods in the presence of severe interference near the mainlobe.

1. Introduction

Monopulse is an established technique in radars for fast and precise estimation of direction of target [1,2]. Compared with other precise angle estimation methods, such as MUSIC, maximum likelihood (ML) estimation etc., it consumes much less time and demands much lower computation cost. It can take only a single snapshot to make a precise estimation while its competitors usually have to solve a non-linear optimization problem [3,4]. In the case of a single target, it is an approximate ML approach, as shown in [5]. However, monopulse technique is not applicable in the case of multiple targets. Large aperture radar nowadays can easily avoid the case of multiple targets with narrower beam width. Modern wide band radars can also employ this technique [6]. Thus, the monopulse technique has still been a promising candidate for angle estimation in modern radar systems. Historically, the monopulse technique is implemented by either constructing error voltages or phase angles from a two-antenna system exploiting a fixed functional relation with the angle of target [7]. In a typical monopulse system, two sets of identical antennas are either separated by some distance (phase monopulse) or located at the the same phase center but with a squint angle difference (amplitude monopulse). Their outputs are summed up to produce a sum beam and are subtracted yielding the difference beam. When it comes to two dimensions, the system needs four receiving channels placed in four quadrants symmetrically to measure azimuth and elevation angles separately. With the advent of phased array radar, operation in the element level is impractical due to the large number of elements. Thus, the monopulse technique is usually implemented with digitalized sub-array outputs. When interference is present, adaptive beamforming [8] is required to suppress the unwanted signals when forming sum and difference beams. More recently, the method of adaptive beamforming based on artificial neural network has been presented in [9].

However, adaptive beamforming changes the beam pattern dependent on the received signals and the desired direction for nulling out the unwanted signals while maximizing the gain in the desired direction. This change results in a deviation in the pointing direction [10] and deformation of the main beam which in turn deteriorates angle measurement accuracy. Additionally, the distortion of the beam pattern is a function of the external environment, such as the location of the interference, target signal energy, clutter energy, etc. [7].

To alleviate distortion in the sum beam, a technique called “diagonal loading” has been introduced [11,12,13]. This technique, in general, makes a compromise between beam distortion and nulling depth by adding a diagonal matrix to the data covariance matrix. It has been proved that diagonal loading is mathematically equivalent to adding constraints about the pointing direction when forming sum beam weights [12]. By choosing a proper loading coefficient, the shape of the main beam can be maintained to be approximately the same while imposing sufficient suppression of the interference.

However, beam distortion may be more severe in forming the difference beam where the “diagonal loading” technique is no longer applicable. The solution to this problem requires keeping sufficient null depth in the direction of interference and maintaining a constant monopulse ratio. Methods based on the ML estimator have been developed in [14,15,16,17] and generalized in [7]. These methods are also referred to as adaptive corrected monopulse as they correct the distortion and bias in the monopulse ratio formula. A different minimum variance adaptive monopulse (MVAM) solution was proposed in [18] and it was found that this solution is equivalent to the ML estimator when the sum beam weight is derived for optimal signal-to-noise ratio (SNR). All these methods have been extended to the space-time case. This has been presented in [19,20,21]. However, all these methods require calculation of first or second derivatives of array manifold with respect to angles of arrival, which makes these methods more sensitive to the error of array elements distribution. Furthermore, these derivatives are calculated at the radar look direction to approximate its true value at the unknown target direction, which would result in performance deterioration. Other ML-based estimation schemes for multi-dimensional parameter estimation were introduced in [22,23], but the solutions are non-linear and ignore the correlation between azimuth and elevation estimates [7]. Another approach to the problem is by adding additional constraints to clamp the values of monopulse ratio at several selected direction [24,25,26]. These constraints work by trading degrees of freedom (DOF) for the desired difference beam pattern around the look direction. The loss of DOF reduces the output signal to noise and interference ratio (SNIR). However, how to choose these constraints to optimize the trade-off remains an open question. Furthermore, the constraints used in [24] ignore the fact that output of the beamforming process in most of the radar systems is actually complex for the sake of measuring the Doppler effect. As a result, the constraints imposed may consume more DOF and lead to further losses in SNIR.

In this paper, a new constrained monopulse algorithm is proposed to calculate adaptive difference beam weights. By using the proposed algorithm, unnecessary loss of DOF is saved, thus improving the performance of monopulse angle measurement. The proposed algorithm is then extended to planar array application, in which the estimation of azimuth and elevation is correlated with each other. Therefore, decoupling constraint is added to the constrained problem, making azimuth and elevation estimation independent from each other and the two-dimensional estimation can be solved by two separate one-dimensional monopulse measurements. The theoretical mean and variance of the proposed monopulse estimator is then derived for analysis using the distribution of the monopulse ratio and the optimized parameters in the constraints are also derived. Theoretical analysis shows that the estimates using the proposed algorithm are unbiased and comparison with the MVAM solution shows that the proposed algorithm can achieve better performance with optimized parameters. In addition, for the first time, we demonstrated that the MVAM could be considered as a special case of the linear constraint monopulse technique. Numerical simulation is employed to demonstrate the effectiveness and advantage of the proposed algorithm and the theoretical variance derived fits the simulation result perfectly well.

The remainder of this paper is organized as follows: in Section 2, a data model for the adaptive monopulse technique is outlined. The proposed algorithm is introduced and a closed form solution generating the adaptive difference beam weights is derived in Section 3.1. The algorithm is then extended to the case of planar arrays in Section 3.2. In Section 4, the theoretical mean and variance of the estimator is derived and analyzed and comparison with the MVAM estimation is made. Mathematical simulations for both the linear array case and planar array case are formulated in Section 5. Finally, concluding remarks are given in Section 6.

2. Data Model

Consider a linear array consisting of m isotropic elements with spacing d. Assume a target signal impinges on the array with amplitude b and phase ξ and incidental angle ${\theta }_0$ which is the physical angle with respect to the line of the array. Besides the target, assume there are p strong sources of interference emitting jamming noise from different incidental angles $\{{\theta }_i:i=1,2\cdots p\}$ and with amplitudes $\{x_i:i=1,2\cdots p\}$ and phases $\{{\xi }_i:i=1,2\cdots p\}$ respectively. Let

$$\mathbf{a}(\theta )=[1e^{j\frac{2\pi d}{\lambda }cos(\theta )}\cdots e^{j\frac{2\pi (m-1)d}{\lambda }cos(\theta )}]$$

(1)

denote the array response to a plane wave of unit amplitude arriving at an incidental angle θ, j is the unit of imaginary numbers and λ is the wavelength of the transmitted signal in radar. Then, the array’s output at the kth moment can be written as

$$\mathbf{y}(k)=b(k)e^{j\xi (k)}\mathbf{a}({\theta }_0)+\sum _{i=1}^p{x}_i(k)e^{j{\xi }_i(k)}\mathbf{a}({\theta }_i)+\mathbf{n}(k),$$

(2)

where $\mathbf{n}(k)$ is the time sample at the k-th moment from white Gaussian noise. With the array output, the monopulse technique requires forming a sum beam and a difference beam by cutting the array into its left half sub-array and right half sub-array. The direction of the mainlobe of the sum beam is determined by the radar look direction θ. However, in the adaptive monopulse radar, a slight deviation may exist due to the presence of interference [12]. If the target is located right at the radar look direction, the sum beam output will be maximum while the difference beam will be zero. If the target has a small offset angle from the look direction ${\theta }_0-\theta$, the ratio of the sum to difference beam, namely the monopulse ratio is employed to estimate this angle. The monopulse ratio is approximately independent of the target amplitude and has a good linear shape over the offset angle. For amplitude monopulse, one can easily show the following basic relation by Taylor expansion [27] with a fixed squint angle δ and antenna amplitude pattern h

$$\frac{h({\theta }_0-\theta -\delta )-h({\theta }_0-\theta +\delta )}{h({\theta }_0-\theta -\delta )+h({\theta }_0-\theta +\delta )}\approx -\frac{h^{\prime }(\delta )}{h(\delta )}·({\theta }_0-\theta ),$$

(3)

where ${()}^{\prime }$ is the first order derivative. However, the sum beam is incoherently formed in this way as we can see from (3), which leads to loss of gain. It is for this reason that we design our algorithm based on the phase monopulse technique. By ignoring noise and interference, the outputs of two sub-arrays have the form $b{e}^{j\xi }h({\theta }_0-\theta )e^{j\frac{2\pi \overline{d}}{\lambda }(cos({\theta }_0)-cos(\theta ))}$ and $b{e}^{j\xi }h({\theta }_0-\theta )e^{-j\frac{2\pi \overline{d}}{\lambda }(cos({\theta }_0)-cos(\theta ))}$ for phase monopulse where h denotes the common sub-array pattern and $\overline{d}$ denotes the distance from the phase center of each sub-array to the center of the whole array. Therefore, one can easily get the monopulse ratio $r=j·tan2\pi \frac{\overline{d}}{\lambda }(cos({\theta }_0)-cos(\theta ))$, which is again a linear function in the neighborhood of the look direction.

In the presence of the interference and noise, the aim is to suppress both the interference and noise as well as extracting the angle to the target. To fulfil this task, appropriate sum and difference beam weights needs to be found that are as orthogonal to the interference subspace as possible while maintaining linearity within the neighborhood of the look direction. Denoting such sum beam weights as ${\mathbf{w}}_{sum}$ and difference beam weights as ${\mathbf{w}}_{diff}$, we have

$$Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{y}(k)}{{\mathbf{w}}_{sum}^H\mathbf{y}(k)}\right)\approx \mathsf{\Gamma }({\theta }_0-\theta ),$$

(4)

where $Im\left(·\right)$ is the imaginary part of a complex variable, ${\left(·\right)}^H$ is the conjugate transpose and Γ is the slope of monopulse ratio. Consequently, the offset angle of target can be easily estimated using (4).

For a planar array, the array manifold becomes

$$\mathbf{a}(\theta ,\phi )=\left[\begin{array}{ccc}1& \cdots & e^{j\frac{2\pi }{\lambda }(m-1)dcos\theta +j\frac{2\pi }{\lambda }(n-1)dcos\phi }\end{array}\right],$$

(5)

where θ and φ are the incident angles w.r.t. the X and Y axis and named azimuth and elevation respectively, m and n are the number of elements along the corresponding dimension. In this case, the estimates of azimuth and elevation are correlated with each other. In fact, by using Taylor expansion, it can be seen that the monopulse ratio is correlated with two partial derivatives w.r.t. azimuth or elevation [18]. Thus, estimation of either azimuth or elevation would require knowledge of the other, which makes the problem harder to approach. However, with difference beam weights appropriately chosen, we can decouple this two-dimensional problem into two independent one-dimensional monopulse estimations and use (4) to solve them separately.

3. The Proposed Constrained Algorithm

3.1. Derivation of the Algorithm Used in a Linear Array

In this section, the proposed algorithm is introduced and the derivation of the constrained difference beam weights given. The procedure commences by first forming the adaptive sum beam weights using a minimum variance distortionless response beamformer (MVDR) [8], e.g.,

$${\mathbf{w}}_{sum}=\frac{{\mathbf{R}}^{-1}\mathbf{a}(\theta )}{{\mathbf{a}}^H(\theta ){\mathbf{R}}^{-1}\mathbf{a}(\theta )},$$

(6)

where $\mathbf{R}$ is the covariance of the interference plus noise and can be estimated by the ML estimator

$$\mathbf{R}=\frac{1}{N}\sum _{k=1}^N\tilde{\mathbf{y}}(k){\tilde{\mathbf{y}}}^H(k),$$

(7)

where N is number of snapshots used for estimation and $\tilde{\mathbf{y}}(k)$ is the sample of interference plus noise, which can be obtained using the interference covariance matrix reconstruction methods [28,29]. The diagonal loading technique can also be applied when forming the sum beam in case mismatch exists between the radar look direction and the target’s true direction.

The calculation of difference in beam weights is dependent on the sum beam weights obtained in (6). The difference beam weights are not only required to obtain sufficient suppression of interference and noise, but also to keep the approximation in (4) valid. Therefore, we formulate the following problem to derive the difference beam weight:

$$\underset{{\mathbf{w}}_{diff}}{min}\left({\mathbf{w}}_{diff}^H\mathbf{R}{\mathbf{w}}_{diff}\right),s.t.Im({\mathbf{w}}_{diff}^H\mathbf{C})=\mathbf{g}$$

(8)

where matrix $\mathbf{C}$ is

$$\mathbf{C}=\left[\frac{\mathbf{a}(\theta +\Delta \theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta +\Delta \theta )}\frac{\mathbf{a}(\theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta )}\right],$$

(9)

and row vector $\mathbf{g}$ is

$$\mathbf{g}=\left[\mathsf{\Gamma }\Delta \theta 0\right],$$

(10)

where $\Delta \theta$ is a small user-determined parameter whose value influences the performance of the proposed algorithm. The choice of $\Delta \theta$ will be discussed shortly.

Using the Lagrange Multiplier method to approach the constrained problem in (8) gives

$$\underset{{\mathbf{w}}_{diff},\mathbf{p}}{min}\left(f={\mathbf{w}}_{diff}^H\mathbf{R}{\mathbf{w}}_{diff}+Im({\mathbf{w}}_{diff}^H\mathbf{C})\mathbf{p})\right),$$

(11)

where $\mathbf{p}$ is a real column vector to be determined. To solve (11), a complex gradient is introduced as follows [30]:

$${\nabla }_{{\mathbf{w}}_{diff}^{*}}(f)=\left(\frac{\partial f}{\partial Re({\mathbf{w}}_{diff})}+j\frac{\partial f}{\partial Im({\mathbf{w}}_{diff})}\right)/2$$

(12)

where $Re(·)$ is the real part of a complex variable. With (11) and (12), we have

$$\begin{array}{c}\hfill {\nabla }_{{\mathbf{w}}_{diff}^{*}}(f)=\mathbf{R}{\mathbf{w}}_{diff}-\frac{j}{2}\mathbf{Cp},\end{array}$$

(13)

where we have used the fact that $\mathbf{R}$ is Hermitian. To find a minimum point for f, set the gradient to zero, then we have

$${\mathbf{w}}_{diff}=\frac{j}{2}{\mathbf{R}}^{-1}\mathbf{Cp}$$

(14)

Substituting (14) into the constraint of (8) gives

$$\mathbf{p}=-2{\left(Re\left({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}\right)\right)}^{-T}{\mathbf{g}}^T,$$

(15)

where ${(·)}^T$ is the transpose of a matrix and ${(·)}^{-T}$ is the transpose of the inverse matrix. Combining (14) and (15) together, we give the solution to (8) as follows:

$${\mathbf{w}}_{diff}=-j{\mathbf{R}}^{-1}\mathbf{C}{\left(Re\left({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}\right)\right)}^{-T}{\mathbf{g}}^T$$

(16)

From (4) and (16), it can be seen that the value of Γ has no effect on the estimate of ${\theta }_0$. Therefore, it can be any real number other than 0. For simplicity, we make $\mathsf{\Gamma }=1$ for the rest of the paper.

After getting the sum beam weight from (6) and the difference beam weight from (16), the monopulse ratio can be easily calculated and therefore the target angle.

The optimization problem in (8) seeks the solution that minimizes the power of interference plus noise under two hard constraints. Since the interference and noise are the source of variance of monopulse estimation, minimizing their power is equal to minimizing the variance of monopulse estimation. The hard constraints clamp the values of monopulse ratio at two close positions, forcing the imaginary part of the monopulse ratio to satisfy (4) within the neighborhood of the radar look direction θ. One of the constraints guarantees that the difference beam has a null at θ, which eliminates estimation bias. The other ensures the approximate linearity within $[\theta -\Delta \theta ,\theta +\Delta \theta ]$, if $\Delta \theta$ is small. The reason for this can be explained as follows. Suppose the target angle being ${\theta }_0=\theta +\tilde{\lambda }\Delta \theta ,\left|\tilde{\lambda }\right|\le 1$ where || means the absolute value or the modulus. By using Taylor expansion for the array manifold and omitting infinitesimal of higher order, we have

$$\mathbf{a}(\theta +\tilde{\lambda }\Delta \theta )\approx \mathbf{a}(\theta )+\tilde{\lambda }(\mathbf{a}(\theta +\Delta \theta )-\mathbf{a}(\theta )),$$

(17)

After beamforming, the output SNIR of both sum and difference beam should be high enough, thus the monopulse ratio can be written as

$$Im(r)\approx Im\left(\frac{{\mathbf{w}}_{diff}^H(\mathbf{a}(\theta )+\tilde{\lambda }(\mathbf{a}(\theta +\Delta \theta )-\mathbf{a}(\theta )))}{{\mathbf{w}}_{sum}^H(\mathbf{a}(\theta )+\tilde{\lambda }(\mathbf{a}(\theta +\Delta \theta )-\mathbf{a}(\theta )))}\right),$$

(18)

Considering $\Delta \theta$ is small and the numerator is much smaller than the denominator in (18), we have the following approximation:

$${\mathbf{w}}_{sum}^H\mathbf{a}(\theta )\approx {\mathbf{w}}_{sum}^H\mathbf{a}(\theta +\Delta \theta )$$

(19)

when calculating the imaginary part of monopulse ratio. From the constraints, we have

$$\begin{array}{c}Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{a}(\theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta )}\right)=0,\hfill \\ Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{a}(\theta +\Delta \theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta +\Delta \theta )}\right)=\Delta \theta .\hfill \end{array}$$

(20)

Combining (18) and (20), one can easily get

$$Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{y}}{{\mathbf{w}}_{sum}^H\mathbf{y}}\right)\approx \tilde{\lambda }\Delta \theta ,\left|\tilde{\lambda }\right|\le 1.$$

(21)

which is equal to (4) within $[\theta -\Delta \theta \theta +\Delta \theta ]$.

From the above analysis, it is clear that the constraints in (20) are able to clamp the monopulse ratio. For convenience, we have assumed that the angle of target lies in $[\theta -\Delta \theta ,\theta +\Delta \theta ]$. However, the angle of target does not have to be confined in this range. As long as $\tilde{\lambda }\Delta \theta$ is small enough to keep (17) and (19) valid, a reasonable estimation of the target’s offset angle can still be obtained by (4). However, with the target moving further away from the radar look direction, the approximation error (non-linearity error) will grow, thus the performance of estimation can be expected to deteriorate. Therefore, as long as $\Delta \theta$ is small, the range of angle estimation does not rely on $\Delta \theta$.

Indeed, the choice of $\Delta \theta$ should be able to keep (17) and (19) valid which demands that $\theta +\Delta \theta$ lies within the neighborhood of the look direction. Denote the range of this neighborhood as ε. In practice, ε is usually considered to be the 3 dB beamwidth of the sum beam. With the target also lying within the neighborhood of the look direction, the approximation error is usually negligible compared to the error incurred by noise and interference. Therefore, the choice of $\Delta \theta$ should be made to minimize the power of noise and interference

$$\Delta \theta =\underset{\Delta \theta }{min}\left({\mathbf{w}}_{diff}^H\mathbf{R}{\mathbf{w}}_{diff}\right),s.t.\Delta \theta \in [\theta -\epsilon /2,\theta +\epsilon /2]\backslash \left\{0\right\}.$$

(22)

With (16), we have

$$\begin{array}{c}{\mathbf{w}}_{diff}^H\mathbf{R}{\mathbf{w}}_{diff}\hfill \\ =\mathbf{g}Re{\left({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}\right)}^{-1}{\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}Re{\left({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}\right)}^{-1}{\mathbf{g}}^H.\hfill \end{array}$$

(23)

Combining (22) and (23), we can get

$$\begin{array}{c}\Delta \theta =\underset{\Delta \theta }{min}\left(F=\frac{{\Delta \theta }^2}{det(Re({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}))}\right),\hfill \\ s.t.\Delta \theta \in [\theta -\epsilon /2\theta +\epsilon /2]\backslash \left\{0\right\}.\hfill \end{array}$$

(24)

which is a non-linear optimization problem. However, using Taylor expansion for the array manifold and omitting infinitesimal of higher order, we can rewrite (24) as

$$F\approx \frac{{\mu }^2}{{\left|\mu +\eta \Delta \theta \right|}^2(\mu \nu -{\left|\eta \right|}^2)},$$

(25)

where $\mu ={\mathbf{a}}^H(\theta ){\mathbf{R}}^{-1}\mathbf{a}(\theta )$, $\nu ={\mathbf{a}}^{\prime H}(\theta ){\mathbf{R}}^{-1}{\mathbf{a}}^{\prime }(\theta )$, $\eta ={\mathbf{a}}^{\prime H}(\theta ){\mathbf{R}}^{-1}\mathbf{a}(\theta )$. The denominator of (25) is a quadratic function of $\Delta \theta$. Its axis of symmetry is at $\Delta \theta =-\frac{Re(\eta )\mu }{{\left|\eta \right|}^2}$. From Cauchy–Schwarz inequality, $\mu \nu \ge {\left|\eta \right|}^2$. Therefore, the choice of $\Delta \theta$ depends on the sign of $Re(\eta )$:

$$\Delta \theta =\left\{\begin{array}{cc}\frac{\epsilon }{2}& \mathrm{if}Re(\eta )>0\\ -\frac{\epsilon }{2}& \mathrm{if}Re(\eta )<0\end{array}\right.$$

(26)

It should be noted that the performance improvement brought by optimizing $\Delta \theta$ depends on the value of $-\frac{Re(\eta )\mu }{{\left|\eta \right|}^2}$. If $\left|\frac{Re(\eta )\mu }{{\left|\eta \right|}^2}\right|\gg \epsilon /2$, the improvement will be more prominent. However, if $\left|\frac{Re(\eta )\mu }{{\left|\eta \right|}^2}\right|\to 0$, this improvement would be less obvious.

Function F is not continuous at $\Delta \theta =0$, but we can make it continuous by assigning value to it. Denote

$$\begin{array}{cc}\hfill F(0)& =\underset{\Delta \theta \to 0}{lim}\frac{{\Delta \theta }^2}{det\left(Re\left({\mathbf{C}}^H{\mathbf{R}}^{-1}\mathbf{C}\right)\right)}\hfill \\ & \approx \frac{1}{\mu \nu -{\left|\eta \right|}^2}.\hfill \end{array}$$

(27)

In fact, as $\Delta \theta$ approaches zero from positive axis, we have

$$\underset{\Delta \theta \to 0}{lim}\frac{Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{a}(\theta +\Delta \theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta +\Delta \theta )}\right)-Im\left(\frac{{\mathbf{w}}_{diff}^H\mathbf{a}(\theta )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta )}\right)}{\Delta \theta }=1.$$

(28)

The constraints in (20) are equal to the demanding derivative of the imaginary part of the monopulse ratio at the radar look direction, to be one with a bias of zero. This case, as we will show in next section, is approximately equal to the MVAM.

Since the constraints are enough to clamp the monopulse ratio, adding extra constraint has very little effect on improving linearity, but it might incur extra loss of DOF. For this reason, the constraint at $\theta -\Delta \theta$ as in [24] has little effect on improving linearity but could result in more loss of DOF. Furthermore, the constraint used in [24] could not be applied to a complex signal. Therefore, the proposed algorithm is superior to the algorithm in [24].

3.2. Extension to Planar Array Application

In a typical planar phased array, the sum and difference beamforming is usually done at the sub-array level using the digital output of the sub-arrays to reduce the dimension of computation. When adaptive beamforming is applied, a serious problem can arise in that the monopulse ratio curve of one dimension may be dependent on the value of the other. This problem can be severe enough to make the monopulse technique fail. Consider a planar array placed in the X–Y plane of a Cartesian system with the center of array overlapping with the origin. The array manifold is shown in (5).

Unlike the linear array case, to estimate either azimuth or elevation angle, additional constraint is required to make the estimation of one angle independent on the value of the other. Therefore, the constrained problem that generates difference beam weights takes the form:

$$\begin{array}{c}\underset{{\mathbf{w}}_{az}}{min}\left({\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{az}\right),\underset{{\mathbf{w}}_{el}}{min}\left({\mathbf{w}}_{el}^H\mathbf{R}{\mathbf{w}}_{el}\right)\hfill \\ s.t.\hfill \\ \mathrm{Im}\left({\mathsf{\Omega }}^H{\mathbf{C}}_2\right)={\mathbf{G}}_2,\mathsf{\Omega }=\left(\begin{array}{cc}{\mathbf{w}}_{az}& {\mathbf{w}}_{el}\end{array}\right)\hfill \end{array}$$

(29)

where ${\mathbf{w}}_{az}$ and ${\mathbf{w}}_{el}$ are the azimuth and elevation difference beam adaptive weights and

$$\begin{array}{c}{\mathbf{C}}_2=\left(\begin{array}{ccc}\frac{\mathbf{a}(\theta ,\phi )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta ,\phi )}& \frac{\mathbf{a}(\theta +\Delta \theta ,\phi )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta +\Delta \theta ,\phi )}& \frac{\mathbf{a}(\theta ,\phi +\Delta \phi )}{{\mathbf{w}}_{sum}^H\mathbf{a}(\theta ,\phi +\Delta \phi )}\end{array}\right),\hfill \\ {\mathbf{G}}_2=\left(\begin{array}{ccc}0& \Delta \theta & 0\\ 0& 0& \Delta \phi \end{array}\right).\hfill \end{array}$$

(30)

Using the Lagrange Multiplier method gives

$$\mathsf{\Omega }=-j{\mathbf{R}}^{-1}{\mathbf{C}}_2{\left(Re\left({{\mathbf{C}}_2}^H{\mathbf{R}}^{-1}{\mathbf{C}}_2\right)\right)}^{-T}{\mathbf{G}}_2^T.$$

(31)

The reason for the constraints keeping linearity within the neighborhood of the radar look direction can be demonstrated in a way similar to the case of linear arrays. Let $(\theta ,\phi )$ denote the radar look direction. The area of the neighborhood is composed by azimuth ranging from $\left(\theta -\epsilon /2\right)$ to $\left(\theta +\epsilon /2\right)$ and elevation ranging from $\left(\phi -\delta /2\right)$ to $\left(\phi +\delta /2\right)$. Practically, this area is also determined by the 3 dB width of the mainlobe. Assume the following approximation holds within the neighborhood of the radar look direction:

$$\begin{array}{cc}\hfill \mathbf{a}\left(\tilde{\theta },\tilde{\phi }\right)& =\left(1-{\lambda }_{\theta }\right)\left(1-{\lambda }_{\phi }\right)\mathbf{a}\left(\theta ,\phi \right)\hfill \\ & +\left(1-{\lambda }_{\theta }\right){\lambda }_{\phi }\mathbf{a}\left(\theta ,\phi +\Delta \phi \right)\hfill \\ & +{\lambda }_{\theta }\left(1-{\lambda }_{\phi }\right)\mathbf{a}\left(\theta +\Delta \theta ,\phi \right)\hfill \end{array}$$

(32)

$${\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\approx {\mathbf{w}}_{sum}^H\mathbf{a}\left(\tilde{\theta },\tilde{\phi }\right)$$

(33)

where $\tilde{\theta }=\theta +{\lambda }_{\theta }\Delta \theta \in \left[\begin{array}{cc}\theta -\epsilon /2,& \theta +\epsilon /2\end{array}\right]$ and $\tilde{\phi }=\phi +{\lambda }_{\phi }\Delta \phi \in \left[\begin{array}{cc}\phi -\delta /2,& \phi +\delta /2\end{array}\right]$. Then, with the constraints in (29), one can easily have

$$Im\left(\frac{{\mathsf{\Omega }}^H\mathbf{a}\left(\tilde{\theta },\tilde{\phi }\right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\tilde{\theta },\tilde{\phi }\right)}\right)\approx \left(\begin{array}{c}{\lambda }_{\theta }\Delta \theta \\ {\lambda }_{\phi }\Delta \phi \end{array}\right)$$

(34)

Optimization of parameters $\Delta \theta$ and $\Delta \phi$ requires solving the following non-linear problem:

$$\begin{array}{c}\underset{\Delta \theta ,\Delta \phi }{min}\left(tr\left(\mathbf{K}{\mathsf{\Omega }}^H\mathbf{R}\mathsf{\Omega }\right)\right),\mathbf{K}=\left(\begin{array}{cc}{k}_{az}& 0\\ 0& k_{el}\end{array}\right)\hfill \\ s.t.\Delta \theta \in \left[\begin{array}{cc}-\epsilon /2,& \epsilon /2\end{array}\right]\backslash \left\{0\right\},\Delta \phi \in \left[\begin{array}{cc}-\delta /2,& \delta /2\end{array}\right]\backslash \left\{0\right\}\hfill \end{array}$$

(35)

where $tr(·)$ is the trace of a matrix, $\mathbf{K}$ is a weighting matrix that weights the variance of azimuth and elevation difference beam outputs. Unfortunately, there is no easy way to solve this problem as in the linear array case. Thus numerical methods such as the Nelder–Mead method are suggested for this problem. However, as long as the constrained directions $\theta +\Delta \theta$ and $\phi +\Delta \phi$ lie within the neighborhood of the radar look direction, the proposed algorithm with optimized $\Delta \theta$ and $\Delta \phi$ can always achieve a better performance than the previous method, e.g., the MVAM.

3.3. Summary and Computational Complexity of the Proposed Algorithm

We conclude the procedure of the proposed monopulse technique and discuss the computational complexity as follows:

• calculate the adaptive sum beam weights using (6). In this step, calculation of the sample matrix inversion (SMI) is the most expensive. Fortunately, we can use the recursive matrix inversion formula for the one rank updated sample covariance matrix [31] which can reduce the computational complexity to the level of $\mathrm{O}\left(l^2\right)$ where l is the dimension of the array manifold. The other matrix multiplication is also in the order of $\mathrm{O}\left(l^2\right)$ [32,33]. Therefore, the total computational complexity in this step is $\mathrm{O}\left(l^2\right)$.
• Determine the parameter $\Delta \theta$ or the parameter pair $(\Delta \theta ,\Delta \phi )$ as follows: in the linear array case, use (26) to determine $\Delta \theta$, while in the planar array case, solve (35) for $(\Delta \theta ,\Delta \phi )$ or directly make $\Delta \theta$ and $\Delta \phi$ close to zero to approximate the performance of the MVAM if seeking to reduce computational cost. In this step, determining the parameter $\Delta \theta$ according to (26) or directly choosing a small value close to zero does not require any computation. However, in the planar array case, optimization of the parameters requires recursive iteration to solve (35), thus it is not suggested for real-time applications.
• Use (16) for linear array applications or (31) for planar array applications to calculate the constrained difference beam weights. The sample matrix inversion is already calculated in the first step and the other matrix inversion in (16) and (31) is usually negligible because $l\gg 3$ in practice. Therefore, the total computational complexity in this step is in the order of $\mathrm{O}\left(l^2\right)$.
• Perform beamforming with the beam weights calculated in the previous steps and then calculate the monopulse ratio along with the angle estimates. This last step has the computational complexity of $\mathrm{O}\left(2l\right)$.

In general, the total computational complexity of the above procedure is $\mathrm{O}\left(l^2\right)$ which is about the same as the MVAM. In practice, such computational cost is usually affordable for modern radar systems.

It is worth noting that before forming the monopulse ratio, any linear processing, such as coherent integration, can be applied to both beam outputs to improve SNIR as the processing effect on the target signal is identical in both beam outputs.

4. Performance Analysis of the Proposed Monopulse Estimator

4.1. Mean and Variance of the Proposed Estimator

In this section, the theoretical performance of the proposed estimator is analyzed. The derivation follows the ideas developed in [17,34,35] using conditional distributions. All the beam weights are considered to be given quantities. Let $d_x$ and $d_y$ denote the output of the azimuth difference beam and elevation difference beam respectively. s is the output of the sum beam. The beam outputs ${\left(\begin{array}{ccc}{d}_x& d_y& s\end{array}\right)}^T$ are assumed to be circularly-symmetric complex Gaussian distributed with mean $\mathbf{u}$ and covariance matrix $\mathbf{Q}$. The mean vector and covariance matrix can be partitioned as

$$\mathbf{u}=\left(\begin{array}{c}{\mathbf{u}}_d\\ u_s\end{array}\right),\mathbf{Q}=\left(\begin{array}{cc}{\mathbf{Q}}_d& {\mathbf{q}}_{ds}\\ {\mathbf{q}}_{ds}^H& q_s\end{array}\right),$$

(36)

where ${\mathbf{u}}_d$ is the mean of $\mathbf{d}={\left(\begin{array}{cc}{d}_x& d_y\end{array}\right)}^T$, ${\mathbf{Q}}_d$ is its covariance, ${\mathbf{q}}_{ds}=E\left((\mathbf{d}-{\mathbf{u}}_d){(s-u_s)}^H\right)$ is the covariance with s and $E(·)$ means the expectation value. Here, we only consider the case of the deterministic target. Therefore, we have

$${\mathbf{u}}_d=b{e}^{j\xi }{\mathsf{\Omega }}^H\mathbf{a}({\theta }_0,{\phi }_0),u_s=b{e}^{j\xi }{\mathbf{w}}_{sum}^H\mathbf{a}({\theta }_0,{\phi }_0);$$

(37)

where $({\theta }_0,{\phi }_0)$ is the direction of the target and

$${\mathbf{Q}}_d={\mathsf{\Omega }}^H\mathbf{R}\mathsf{\Omega },q_s={\mathbf{w}}_{sum}^H\mathbf{R}{\mathbf{w}}_{sum},{\mathbf{q}}_{ds}={\mathsf{\Omega }}^H\mathbf{R}{\mathbf{w}}_{sum}.$$

(38)

The conditional distribution of $\mathbf{d}$ conditioned on s is also complex Gaussian with mean ${\mathbf{u}}_{d|s}={\mathbf{u}}_d+{\mathbf{q}}_{ds}{q}_s^{-1}(s-u_s)$ and covariance ${\mathbf{Q}}_{d|s}={\mathbf{Q}}_d-{\mathbf{q}}_{ds}{q}_s^{-1}{\mathbf{q}}_{ds}^H$. Rewrite the imaginary part of the monopulse ratio as

$$\mathbf{r}=Im\left(\frac{\mathbf{d}{s}^{*}}{\left|s^2\right|}\right).$$

(39)

Then, the conditional distribution of $\mathbf{r}$ conditioned on s is Gaussian with mean $Im\left({\mathbf{q}}_{ds}{q}_s^{-1}+\left({\mathbf{u}}_d-{\mathbf{q}}_{ds}{q}_s^{-1}{u}_s\right)s^{-1}\right)$ and covariance $Re\left(\frac{{\mathbf{Q}}_{d|s}}{2{\left|s\right|}^2}\right)$. Let $u_{az|s}$ and $q_{az|s}$ denote the mean and variance of the azimuth estimator conditioned on s respectively. With (37) and (38), they can be written as

$$\begin{array}{c}{u}_{az|s}=Im\left(\beta +\left({\mathbf{w}}_{az}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)-\beta {\mathbf{w}}_{sum}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)\right)b{e}^{j\xi }{s}^{-1}\right),\hfill \\ q_{az|s}=\left({\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{az}-\frac{{\left|{\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{sum}\right|}^2}{{\mathbf{w}}_{sum}^H\mathbf{R}{\mathbf{w}}_{sum}}\right){\left(2{\left|s\right|}^2\right)}^{-1},\hfill \\ \beta =\frac{{\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{sum}}{{\mathbf{w}}_{sum}^H\mathbf{R}{\mathbf{w}}_{sum}}.\hfill \end{array}$$

(40)

The performance of the estimator depends on the sum beam weights. In the presence of interference, the MVDR beamformer, as in (6), is usually a practical method for calculating sum beam weights and the proposed monopulse technique is based on it.

Using (6) and the constraint $Im\left(\beta \right)=0$, we can get

$$\begin{array}{c}{u}_{az|s}\hfill \\ =Im\left(\frac{\left({\mathbf{w}}_{az}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)-{\mathbf{w}}_{az}^H\mathbf{a}\left(\theta ,\phi \right){\mathbf{w}}_{sum}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)\right)b{e}^{j\xi }}{s^{-1}}\right),\hfill \\ q_{az|s}\hfill \\ =\left({\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{az}-\frac{{\left|{\mathbf{w}}_{az}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)}\right){\left(2{\left|s\right|}^2\right)}^{-1}.\hfill \end{array}$$

(41)

By integrating over the distribution of s, the mean and covariance of $\mathbf{r}$ can then be obtained. The sum beam output power ${\left|s\right|}^2$ cannot be too low, otherwise the monopulse estimation would fail. Therefore, a threshold α is set for the lower bound of the integration and the mean and variance derived in the following is still conditioned on ${\left|s\right|}^2>\alpha$ [17]. This threshold is relevant with the scenario where the radar system is working, including SNIR, false alarm rate, etc. Using the results from [17], the mean of the azimuth estimator $u_{az}$ is given by

$$u_{az}=Im\left(\frac{{\mathbf{w}}_{az}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)}\right)\left(1-\frac{a_1}{p_d}\right)$$

(42)

where

$$\begin{array}{c}{a}_1=e^{-(\alpha +{\left|u_s\right|}^2)/q_s}·I_0\left(2\sqrt{\alpha }\left|u_s\right|/q_s\right),\hfill \\ p_d=\frac{1}{q_s}{\int }_{\alpha }^{\infty }{e}^{-(t+{\left|u_s\right|}^2)/q_s}·I_0\left(2\sqrt{t}\left|u_s\right|/q_s\right)dt\hfill \end{array}$$

(43)

where $I_0$ denotes the modified Bessel function of order zero. In the derivation of (42), we have used the constraint $Im\left(\beta \right)=0$. The mean of the elevation estimator can be derived analogously:

$$u_{el}=Im\left(\frac{{\mathbf{w}}_{el}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left({\theta }_0,{\phi }_0\right)}\right)\left(1-\frac{a_1}{p_d}\right)$$

(44)

As α approaches zero, $a_1$ is approaching a small value determined by the SNR of the sum beam and $p_d$ is approaching unity. Thus, the proposed estimator is approximately unbiased.

With the results from [17], the variance of the azimuth and elevation estimator is

$$\begin{array}{cc}\hfill q_{az}& =\frac{1}{2}\left({\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{az}-\frac{{\left|{\mathbf{w}}_{az}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)}\right)\frac{a_2}{p_d}\hfill \\ & =\frac{1}{2}\left({\mathbf{w}}_{az}^H{\mathbf{R}}^{1/2}{\mathbf{P}}_{null}{\mathbf{R}}^{1/2}{\mathbf{w}}_{az}\right)\frac{a_2}{p_d},\hfill \\ \hfill q_{el}& =\frac{1}{2}\left({\mathbf{w}}_{el}^H\mathbf{R}{\mathbf{w}}_{el}-\frac{{\left|{\mathbf{w}}_{el}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)}\right)\frac{a_2}{p_d}\hfill \\ & =\frac{1}{2}\left({\mathbf{w}}_{el}^H{\mathbf{R}}^{1/2}{\mathbf{P}}_{null}{\mathbf{R}}^{1/2}{\mathbf{w}}_{el}\right)\frac{a_2}{p_d},\hfill \end{array}$$

(45)

where ${\mathbf{P}}_{null}$ is the projection matrix to the null space of ${\mathbf{R}}^{-1/2}\mathbf{a}\left(\theta ,\phi \right)$ and

$$a_2={\int }_{\alpha }^{\infty }\frac{1}{q_s}{e}^{-(t+{\left|u_s\right|}^2)/q_s}·I_0\left(2\sqrt{t}\left|u_s\right|/q_s\right)t^{-1}dt.$$

(46)

With (31), we can rewrite $Re\left({\mathsf{\Omega }}^H\mathbf{a}\left(\theta ,\phi \right)\right)$ as

$$\begin{array}{c}Re\left({\mathsf{\Omega }}^H\mathbf{a}\left(\theta ,\phi \right)\right)\hfill \\ ={\mathbf{G}}_2\left(Re{\left({{\mathbf{C}}_2}^H{\mathbf{R}}^{-1}{\mathbf{C}}_2\right)}^{-1}\right)Im\left({{\mathbf{C}}_2}^H{\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)\right).\hfill \end{array}$$

(47)

With the constraints in (29), one can easily get $Im\left({{\mathbf{C}}_2}^H{\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)\right)=\mathbf{0}$. Therefore, we have

$$\begin{array}{c}{q}_{az}=\frac{1}{2}\left({\mathbf{w}}_{az}^H\mathbf{R}{\mathbf{w}}_{az}\right)\frac{a_2}{p_d},\hfill \\ q_{el}=\frac{1}{2}\left({\mathbf{w}}_{el}^H\mathbf{R}{\mathbf{w}}_{el}\right)\frac{a_2}{p_d}.\hfill \end{array}$$

(48)

As we can see from (48), $\frac{a_2}{p_d}$ is determined by the SNR of the sum beam output and the threshold α, the other term is the variance of the difference beam output. Therefore, given the sum beam weights and α, the variance of the estimator only depends on the variance of the difference beam output which hence explains (22) and (35) for the optimization of parameter $\Delta \theta$.

4.2. Comparison with MVAM

In the MVAM [18], the adaptive difference beam weights were derived under different objective function and different constraints from (29). However, as will be proved in the following, the adaptive difference beam weights calculated by the MVAM are equivalent with that of the proposed algorithm in a special case where $\Delta \theta \to 0$. Therefore, the proposed algorithm can outperform the MVAM with the parameter $\Delta \theta$ optimized using (26) or in the planar array case, using numerical methods. Let $\tilde{\mathsf{\Omega }}=\left(\begin{array}{cc}{\tilde{\mathbf{w}}}_{az}& {\tilde{\mathbf{w}}}_{el}\end{array}\right)$ denote the adaptive difference beam weights given by the MVAM. These weights are calculated from [18]:

$$\begin{array}{c}\underset{{\tilde{\mathbf{w}}}_{az}}{min}\left({\tilde{\mathbf{w}}}_{az}^H\tilde{\mathbf{R}}{\tilde{\mathbf{w}}}_{az}\right),\underset{{\tilde{\mathbf{w}}}_{el}}{min}\left({\tilde{\mathbf{w}}}_{el}^H\tilde{\mathbf{R}}{\tilde{\mathbf{w}}}_{el}\right)\hfill \\ s.t.Re\left({\tilde{\mathbf{C}}}^H\tilde{\mathsf{\Omega }}\right)={\mathbf{I}}_2\hfill \end{array}$$

(49)

where ${\mathbf{I}}_2$ is the unity matrix of rank 2 and

$$\begin{array}{cc}\tilde{\mathbf{R}}\hfill & =\mathbf{R}-\frac{\mathbf{a}\left(\theta ,\phi \right){\mathbf{w}}_{sum}^H}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)}\mathbf{R}-\mathbf{R}\frac{{\mathbf{w}}_{sum}\mathbf{a}{\left(\theta ,\phi \right)}^H}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{w}}_{sum}}\hfill \\ & +\frac{\mathbf{a}\left(\theta ,\phi \right){\mathbf{a}}^H\left(\theta ,\phi \right)}{{\left|{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2}{\mathbf{w}}_{sum}^H\mathbf{R}{\mathbf{w}}_{sum}\hfill \\ \tilde{\mathbf{C}}\hfill & =\left(\begin{array}{cc}{\mathbf{c}}_{az}& {\mathbf{c}}_{el}\end{array}\right),\hfill \\ {\mathbf{c}}_{az}\hfill & =\frac{\left({\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\right){\mathbf{d}}_{\theta }-\left({\mathbf{w}}_{sum}^H{\mathbf{d}}_{\theta }\right)\mathbf{a}\left(\theta ,\phi \right)}{{\left|{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2},\hfill \\ {\mathbf{c}}_{el}\hfill & =\frac{\left({\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\right){\mathbf{d}}_{\phi }-\left({\mathbf{w}}_{sum}^H{\mathbf{d}}_{\phi }\right)\mathbf{a}\left(\theta ,\phi \right)}{{\left|{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)\right|}^2},\hfill \end{array}$$

(50)

where ${\mathbf{d}}_{\theta }$ and ${\mathbf{d}}_{\phi }$ are respectively the partial derivatives of $\mathbf{a}\left(\theta ,\phi \right)$ w.r.t. θ and φ at the look direction. The singularity of $\tilde{\mathbf{R}}$ prevents solving (49) directly with the Lagrange method. In [18], this problem was solved by adding an additional constraint ${\tilde{\mathsf{\Omega }}}^H{\mathbf{w}}_{sum}=\mathbf{0}$.

Because ${\mathbf{R}}^{-1}=\mathbf{M}\mathbf{P}{\mathbf{M}}^H+\frac{1}{{\sigma }^2}\mathbf{I}$ where $\mathbf{M}$ is the interference eigenvector matrix and ${\sigma }^2$ is the power of the noise, it follows that

$${\tilde{\mathsf{\Omega }}}^H{\mathbf{w}}_{sum}=\frac{\left({\tilde{\mathsf{\Omega }}}^H\mathbf{MP}{\mathbf{M}}^H\mathbf{a}\left(\theta ,\phi \right)+\frac{1}{{\sigma }^2}{\tilde{\mathsf{\Omega }}}^H\mathbf{a}\left(\theta ,\phi \right)\right)}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)}.$$

(51)

Because the power of interference is assumed to be much stronger than the noise, the maximum eigenvalue of $\mathbf{P}$ is much smaller than $\frac{1}{{\sigma }^2}$ and the difference beam weights should be orthogonal to the sub-space of interference. Therefore, we have

$${\tilde{\mathsf{\Omega }}}^H{\mathbf{w}}_{sum}\approx \frac{1}{{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)}\left(\frac{1}{{\sigma }^2}{\tilde{\mathsf{\Omega }}}^H\mathbf{a}\left(\theta ,\phi \right)\right)$$

(52)

and further the constraint ${\tilde{\mathsf{\Omega }}}^H{\mathbf{w}}_{sum}=\mathbf{0}$ approximately equals to ${\tilde{\mathsf{\Omega }}}^H\mathbf{a}\left(\theta ,\phi \right)=\mathbf{0}$. Thus, with the additional constraint ${\tilde{\mathsf{\Omega }}}^H\mathbf{a}\left(\theta ,\phi \right)=\mathbf{0}$ added to (49), the solution of the MVAM is given by:

$$\begin{array}{c}\tilde{\mathsf{\Omega }}={\overline{\mathbf{R}}}^{-1}\overline{\mathbf{C}}{\left(Re\left({\overline{\mathbf{C}}}^H{\overline{\mathbf{R}}}^{-1}\overline{\mathbf{C}}\right)\right)}^{-1}\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right),\hfill \\ \overline{\mathbf{R}}=\tilde{\mathbf{R}}+a_3\mathbf{a}\left(\theta ,\phi \right){\mathbf{a}}^H\left(\theta ,\phi \right),\overline{\mathbf{C}}=\left(\begin{array}{cc}\tilde{\mathbf{C}}& \mathbf{a}\left(\theta ,\phi \right)\end{array}\right)\hfill \end{array}$$

(53)

where $a_3$ is a non-zero loading factor that makes $\overline{\mathbf{R}}$ invertible. The effect of the loading matrix is cancelled by the additional constraint.

Let $a_3=1/{\mathbf{a}}^H\left(\theta ,\phi \right){\mathbf{R}}^{-1}\mathbf{a}\left(\theta ,\phi \right)$, then the weight solution becomes

$$\tilde{\mathsf{\Omega }}={\mathbf{R}}^{-1}\overline{\mathbf{C}}{\left(Re\left({\overline{\mathbf{C}}}^H{\mathbf{R}}^{-1}\overline{\mathbf{C}}\right)\right)}^{-1}\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right)$$

(54)

Let us rewrite $\overline{\mathbf{C}}$ as follows:

$$\begin{array}{c}\overline{\mathbf{C}}=\underset{\begin{array}{c}\Delta \theta \to 0\\ \Delta \phi \to 0\end{array}}{lim}\left(\Pi \right),\Pi =\left(\begin{array}{ccc}{\mathbf{z}}_1& {\mathbf{z}}_2& {\mathbf{z}}_3\end{array}\right),\hfill \\ {\mathbf{z}}_1=\frac{1}{\Delta \theta }\left(\frac{\mathbf{a}\left(\theta +\Delta \theta ,\phi \right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta +\Delta \theta ,\phi \right)}-\frac{\mathbf{a}\left(\theta ,\phi \right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)}\right),\hfill \\ {\mathbf{z}}_2=\frac{1}{\Delta \phi }\left(\frac{\mathbf{a}\left(\theta ,\phi +\Delta \phi \right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi +\Delta \phi \right)}-\frac{\mathbf{a}\left(\theta ,\phi \right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)}\right),\hfill \\ {\mathbf{z}}_3=\frac{\mathbf{a}\left(\theta ,\phi \right)}{{\mathbf{w}}_{sum}^H\mathbf{a}\left(\theta ,\phi \right)}.\hfill \end{array}$$

(55)

Combining (54) and (55), we have

$$\begin{array}{c}\tilde{\mathsf{\Omega }}=\underset{\begin{array}{c}\Delta \theta \to 0\\ \Delta \phi \to 0\end{array}}{lim}\left({\mathbf{R}}^{-1}\Pi {\left(Re\left({\Pi }^H{\mathbf{R}}^{-1}\Pi \right)\right)}^{-1}\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right)\right)\hfill \\ =\underset{\begin{array}{c}\Delta \theta \to 0\\ \Delta \phi \to 0\end{array}}{lim}\left(j\mathsf{\Omega }\right).\hfill \end{array}$$

(56)

where Ω is calculated by (31). It can be seen from (56) that the MVAM corresponds to a special case of the proposed algorithm where $\Delta \theta \to 0$ and $\Delta \phi \to 0$. In the case of linear array, one can prove that the variance of estimation using the MVAM is ${\left(\mu \nu -{\left|\eta \right|}^2\right)}^{-1}$ which is (27).

Denote $G(\Delta \theta ,\Delta \phi )=tr\left(\mathbf{K}{\mathsf{\Omega }}^H\mathbf{R}\mathsf{\Omega }\right)$ and define $G(0,0)$ by its limit value. Then, function G is continuous almost everywhere within the domain $\Delta \theta \in \left[\begin{array}{cc}-\epsilon /2,& \epsilon /2\end{array}\right],\Delta \phi \in \left[\begin{array}{cc}-\delta /2,& \delta /2\end{array}\right]$. Therefore, by optimizing the choice of $\Delta \theta$ and $\Delta \phi$, the proposed algorithm could achieve better performance than the MVAM or Nickel’s corrected adaptive monopulse [18] which gives the same performance as the MVAM. In the case of linear array, the strategy for choosing the parameter is given by (26).

Furthermore, the proposed algorithm does not require knowledge of the derivatives of array response. Thus, it is more robust than the MVAM in the case where errors about the array response exist.

5. Numerical Examples and Applications

In order to examine the performance of the proposed monopulse algorithm, mathematical simulations are conducted using an echo simulation software programmed in MATLAB. Comparison is made between the proposed algorithm, traditional constrained method [24] and the MVAM. The simulation is composed of two parts. In the first part, different methods are compared using a uniform linear array, while in the second they are applied to a rectangular planar array for comparison.

5.1. Simulation in Linear Array

Consider a linear array with 16 elements equally spaced by half of the wavelength. The desired radar look direction is 20${}^{\circ }$ away from the boresight. A deterministic target is located at the direction 20.8${}^{\circ }$ off the boresight with the power of 10 dB relative to internal noise. Two noise jammers with the power of 30 dB relative to internal noise are incident from 10${}^{\circ }$ and 14${}^{\circ }$ off the boresight respectively. The sum beam is formed adaptively using (6). The beam patterns using the proposed algorithm and the MVAM are plotted in Figure 1. As a reference, the quiescent beam patterns are also plotted. For the proposed algorithm, $\Delta \theta$ is set to be 4${}^{\circ }$ (half of the 3 dB beamwidth), because $Re\left(\eta \right)>0$. It can be seen from Figure 1 that the mainlobe of the sum beam pattern is slightly nudged away from the desired look direction due to the strong jammers. Two clear nulls can be observed at the directions of jammers for both of the methods, providing jammer rejection. However, the proposed algorithm places another null at the radar look direction because of the constraint while the MVAM only tends to have a null at the look direction [18].

In Figure 2, the standard deviations of the proposed algorithm, the MVAM and the traditional constrained method are plotted with signal to noise ratio (SNR) changing from −3 dB to 13 dB and Root-Mean-Square Errors (RMSE) of these three methods are plotted in Figure 3. The covariance matrix $\mathbf{R}$ is estimated from 100 snapshots. For the traditional constrained method, the constrained points are respectively −4${}^{\circ }$, 0${}^{\circ }$ and 4${}^{\circ }$ away from the radar look direction. It can be seen that the proposed algorithm outperforms both the MVAM and the traditional constrained method, especially in the case of low SNR. Comparing the standard deviation and RMSE at each SNR scenario, we can tell that all these three methods give approximately unbiased estimates of the target direction. Additionally, the theoretical standard deviation calculated by (48) predicts the simulation results well. However, at low SNR, due to the disturbance of noise and a finite number of samples, the simulation results slightly deviate from the theoretical curve.

The effect of choosing different $\Delta \theta$ is evaluated in Figure 4 where the MVAM is used as a comparison. SNR is set to be 10 dB. It can be seen that the RMSE of the proposed algorithm slowly decreases with $\Delta \theta$ growing until $\Delta \theta$ reaches half of the beamwidth and the RMSE near 0 is approximately equal to that of the MVAM. This result shows that by optimizing $\Delta \theta$, the proposed algorithm can achieve better performance than the MVAM. It can also be observed that the standard deviation continues to decrease with $\Delta \theta$ growing beyond the half of beamwidth which is analyzed in section III. However, the RMSE deteriorates after that point which suggests that, in this case, the bias of the estimation cannot be ignored.

Finally, the RMSE versus target angle is plotted in Figure 5. When the target is very close to the radar look direction, the proposed algorithm with $\Delta \theta =4$ achieves much better performance than the MVAM and the traditional constrained method. However, its performance deteriorates fast with the target offset angle growing and would be worse than the MVAM when the offset angle continues growing. This deterioration is due to the bias increase of the estimation. As can be seen in Figure 5, when $\Delta \theta$ is adjusted to 1, the proposed algorithm has better performance than the MVAM at every offset angle of target. However, compared with the case where $\Delta \theta =4$, its RMSE is worse at positions close to the radar look direction.

From the above simulation results, we can see that the proposed algorithm always outperforms the traditional constrained method and when $\Delta \theta$ goes to 0, it gives the same performance as the MVAM. With proper choice of $\Delta \theta$, the proposed algorithm can outperform the MVAM.

5.2. Simulation in Planar Array

In this section, a simulation is conducted based on a $12\times 12$ rectangular planar array with elements horizontally and vertically spaced by half of the wavelength. This whole array is divided into 16 $3\times 3$ sub-arrays all steered to look at the direction $({70}^{\circ },{90}^{\circ })$. For simplicity of comparison and discussion, no amplitude weighting is considered at the array elements. The target is located at $({69}^{\circ },{89}^{\circ })$ with power equal to the internal noise. Two noise jammers with power of 30 dB relative to the internal noise are incident from the direction $({77}^{\circ },{90}^{\circ })$ near the mainlobe and the direction $({84}^{\circ },{90}^{\circ })$ in the sidelobe region. The time series data is passed through a 100-point fast Fourier transform (FFT) filter before applying the adaptive weighting. The covariance matrix $\mathbf{R}$ is estimated from 100 snapshots.

In the first scenario, instead of solving (35), we first evaluate two choices of $\Delta \theta$: $\Delta \theta =3^{\circ }$ and $\Delta \theta =-3^{\circ }$, whose absolute value is almost half of the beamwidth. The other parameter $\Delta \phi$ is chosen to be 0.1. The azimuth and elevation difference beam pattern cuts for the case of $\Delta \theta =-3^{\circ }$ are plotted in Figure 6. It can be observed that nulls are placed at the directions of jammers for both the proposed algorithm and the MVAM and like the case of linear array, the MVAM tends to place a null at the radar look direction.

The RMSE of azimuth estimates versus the SNR are plotted in Figure 7. It can be seen that the proposed algorithm outperforms the MVAM with $\Delta \theta =-3^{\circ }$ but performs worse than the MVAM in the case of $\Delta \theta =3^{\circ }$. Therefore, to achieve the best performance, the parameters need to be optimized. Combining with the simulation results in the case of linear array, it can be observed that in all these scenarios, the proposed algorithm with optimized parameters outperforms the MVAM.

The RMSE of elevation estimates versus the SNR are plotted in Figure 8. From Figure 8, we can see that the proposed algorithm and the MVAM have the same performance at each SNR. Therefore, we demonstrate that by choosing $\Delta \phi$ close to 0, the proposed algorithm gives the same performance as the MVAM.

In the second scenario, we evaluate the performance w.r.t. different target locations. Assume the target is incident from the direction $({71}^{\circ },{91}^{\circ })$ while the rest of the conditions remain the same as the first scenario. Again, two sets of parameters $\Delta \theta =3^{\circ },\Delta \phi ={0.1}^{\circ }$ and $\Delta \theta =-3^{\circ },\Delta \phi ={0.1}^{\circ }$ are evaluated. The RMSE of estimates versus the SNR are plotted in Figure 9 and Figure 10. It can be seen from Figure 9 that in the estimation of azimuth, the proposed algorithm with $\Delta \theta =-3^{\circ }$ still outperforms the MVAM. In Figure 10, it can be seen that the estimation of elevation is irrelevant with the choice of $\Delta \theta$ and with $\Delta \phi ={0.1}^{\circ }$, the proposed algorithm is equivalent with the MVAM. More results with different target locations are shown in

Table 1. RMSE of azimuth estimates versus target locations.

Method	Case 1	Case 2	Case 3	Case 4
Proposed algorithm	${0.79}^{\circ }$	${0.60}^{\circ }$	${1.10}^{\circ }$	${0.64}^{\circ }$
MVAM	${0.93}^{\circ }$	${0.72}^{\circ }$	${1.20}^{\circ }$	${0.79}^{\circ }$

and

Table 2. RMSE of elevation estimates versus target locations.

Method	Case 1	Case 2	Case 3	Case 4
Proposed algorithm	${0.27}^{\circ }$	${0.22}^{\circ }$	${0.72}^{\circ }$	${0.37}^{\circ }$
MVAM	${0.27}^{\circ }$	${0.22}^{\circ }$	${0.72}^{\circ }$	${0.37}^{\circ }$

. In these cases, the SNR is fixed at −6 dB. From case 1 to 4, the target is incident from $({71}^{\circ },{89}^{\circ })$, $({69}^{\circ },{91}^{\circ })$, $({72}^{\circ },{92}^{\circ })$ and $({68}^{\circ },{88}^{\circ })$ respectively.

In the third scenario, we evaluate the performance of the proposed algorithm with different jammer locations. The target is located at $({69}^{\circ },{89}^{\circ })$ with the SNR being −5 dB. We consider three different cases. In the first case, assume two jammers are incident from $({77}^{\circ },{90}^{\circ })$ and $({70}^{\circ },{97}^{\circ })$ respectively. In the second case, they are $({63}^{\circ },{90}^{\circ })$, $({56}^{\circ },{90}^{\circ })$ and in the third case, $({63}^{\circ },{90}^{\circ })$, $({70}^{\circ },{83}^{\circ })$. The rest of the conditions are the same as scenario one. The RMSE of the proposed algorithm and the MVAM are shown in

Table 3. RMSE of azimuth estimates versus jammer locations.

Method	Case 1	Case 2	Case 3
Proposed algorithm	${0.33}^{\circ }$	${0.81}^{\circ }$	${0.42}^{\circ }$
MVAM	${0.33}^{\circ }$	${0.96}^{\circ }$	${0.47}^{\circ }$

and

Table 4. RMSE of elevation estimates versus jammer locations.

Method	Case 1	Case 2	Case 3
Proposed algorithm	${0.32}^{\circ }$	${0.21}^{\circ }$	${0.37}^{\circ }$
MVAM	${0.33}^{\circ }$	${0.28}^{\circ }$	${0.40}^{\circ }$

. For the proposed algorithm, parameters are determined by solving (35) with grid search within 3 dB beamwidth around the radar look direction. The weighting matrix $\mathbf{K}$ is chosen to be the unity matrix. Since the range of searching area is small, the computation cost is reasonable. For case 1, the optimum parameters are $\Delta \theta =-3^{\circ }$ and $\Delta \phi =-3^{\circ }$. For case 2, the optimum parameters are $\Delta \theta =3^{\circ }$ and $\Delta \phi =-3^{\circ }$ while for case 3, they are $\Delta \theta =3^{\circ }$ and $\Delta \phi =3^{\circ }$. It can be seen that the proposed algorithm outperforms the MVAM in all three cases. However, the performance of the proposed algorithm depends on the locations of jammers. In the case where both jammers have the same elevation angle (case 2), the advantage of the proposed algorithm is more obvious than the other cases.

6. Conclusions

In this paper, a novel constrained monopulse algorithm has been proposed for phased array radar. In the real battlefield, airborne target signals are usually masked by strong noise jammers which are close to the target in space. By using novel adaptive difference beam weights, the difference beam adaptively provides suppression against unwanted signals and maintains the beam shape around the look direction. Therefore, the monopulse ratio formed by the proposed algorithm is able to maintain the linearity around the look direction and minimize the effect from noise jamming. The theoretical mean and variance of the proposed monopulse estimator is derived for performance evaluation. From the theoretical analysis and mathematical simulation, it can be seen that the estimation of the target’s direction is unbiased and accurate. In comparison with the previous methods, it is found that the MVAM is equivalent with a special case of the proposed algorithm where the constrained point is chosen to be close to the look direction. By optimizing the constrained point, the proposed algorithm can outperform the existing adaptive monopulse techniques. The optimized constrained point is derived in the linear array case. In the planar array case, the optimization requires recursive iteration, thus is not suggested for real-time application. However, in such case, the proposed algorithm can still achieve the same performance as the MVAM by making the constrained points close to the look direction. Mathematical simulation results in scenarios where severe jamming sources occur near the mainlobe, demonstrate that the proposed monopulse technique can accurately acquire target direction and outperform the previous adaptive monopulse techniques. In addition, the computational complexity of the proposed algorithm is found to be $\mathrm{O}\left(l^2\right)$ which is about the same as the traditional adaptive monopulse techniques and usually affordable in modern phased array radar. All in all, the proposed algorithm proves to be effective in the real scenarios where strong jammers occur near the target.