Tracking a Low-Angle Isolated Target via an Elevation-Angle Estimation Algorithm Based on Extended Kalman Filter with an Array Antenna

In a low-angle tracking situation, estimating the elevation angle is challenging due to the entrance of the multipath signals in the antenna’s main lobe. In this article, we propose two methods based on the extended Kalman filter (EKF) and frequency diversity (FD) process to estimate the elevation angle of a low-angle isolated target. In the first case, a simple weighting of the per-frequency estimates is obtained (termed WFD). Differently, in the second case, a matrix-based elaboration of the per-frequency estimates is proposed (termed MFD). The proposed methods are completely independent of prior knowledge of geometrical information and the physical parameters. The simulation results show that both methods have excellent performance and guarantee accurate elevation angle estimation in different multipath environments and even in very-low SNR conditions. Hence, they are both suitable for low-peak-power radars.


Introduction
For nearly half a century, researchers have been studying the problem of localizing low-elevation targets. Unfortunately, in the above scenario, the multipath effect causes a significant error in tracking the target elevation angle, even in the case of an isolated target. Accordingly, multipath represents the main obstacle in solving the aforementioned task accurately. Indeed, the antenna array receives more signals than one direct signal while tracking a low-angle target flying over the surface. Specifically, specular and diffuse components reflected from the surface are also received by the array and create a (ghost) image target below the surface.
These multipath components lead to an interfering effect in the elevation angle estimation process. Such an effect may be destructive or constructive, depending on the traveled distance from the target to the antenna and the phase difference between direct and multipath signals. However, the specular signal is more destructive due a relativelystronger effect and higher correlation with the direct signal. Accordingly, the conventional monopulse method fails in tracking a low-angle target in the presence of multipath. To mitigate such detrimental effects, White [1] enhanced the monopulse method by introducing the double null (DN) and fixed beam (FB) methods.
These two methods were improved in [2] by squinting the antenna's sum and difference patterns. Differently, in [3], the difference monopulse pattern was designed symmetrically via beamforming. The use of the frequency diversity (FD) in an optimal way to reduce the multipath effect is investigated in [4]. by predicting the value, the velocity and the acceleration of the target elevation angle based on previous estimations.
To further increase the accuracy of the proposed methods, we rely on a FD process. In the first method, we use a (rank-based) weighting of the per-frequency angle elevation estimates. Conversely, in the second method, we use a matrix-based processing of the aforementioned per-frequency estimates via a compact representation of the multiplemeasurement model.
The proposed methods can properly work in both rough and smooth surfaces and do not require prior knowledge of the target range, the surface reflection coefficient or sea state, the radar height, the polarization and the effective earth radius. The computational complexity of the proposed methods is relatively low, and yet they present significant performance in tracking low-angle targets even for a low number of antenna array elements. Still, the proposed methods show excellent performance even at very-low SNR (as opposed to other existing methods [15,16]), thus, being suitable for adoption in low-peakpowers radars.
The remainder of the manuscript is organized as follows. The propagation and multipath models are introduced in Section 2, and the proposed methods for elevation angle estimation are presented in Section 3. Section 4 presents (and discusses) the simulation results to assess the effectiveness of the proposed methods. Finally, our conclusions and future directions for research are given in Section 5.
Notation: in this paper, we use lower-case (resp. upper-case) bold letters denote vectors (resp. matrices), with a k (resp. a n,m ) representing the kth element (resp. (n, m)th element) of a (resp. A); i is used to denote the imaginary unit; the symbols (.) T and (.) H are used to indicate transpose and conjugate transpose, respectively; I m denotes the m-dimensional identity matrix; N (µ, Σ) denotes a normal distribution with mean vector µ and covariance matrix Σ; Rayleigh(ρ) denotes a Rayleigh distribution with parameter ρ; U (a, b) denotes a uniform distribution with support (a, b); and finally, the symbols ∼ and ≈ mean "distributed as" and "approximately equal to", respectively.

Propagation and Multipath Models
The present section describes the propagation and multipath models considered in this work. Specifically, in Section 2.1, we first describe the general antenna array model. Then, in Section 2.2, we detail the geometric parameters involved. Finally, in Sections 2.3 and 2.4, we report the considered specular and diffuse component modeling.

General Propagation Model
As shown in Figure 1, we consider a curved-earth model with both specular and diffuse propagation. The considered radar is equipped with a uniform linear array (ULA) made of N elements with inter-element spacing d. The received array signal (observation) vector x(n) x 1 (n) x 2 (n) · · · x N (n) T , including the direct and multipath signals, is given as where a(θ) 1 e −ikd sin(θ) · · · e −ik(N−1)d sin(θ) T denotes the steering vector, with k = (2π/λ) and λ representing the wave-number and the wavelength, respectively. In Equation (1), the terms θ t , θ r and θ d represent the arrival angles associated to the direct, specular and diffuse component, respectively. Analogously, s(n), s r (n) and s d (n) denote the amplitudes of the direct, specular and diffuse signals, respectively. Finally, w(n) is a zero-mean complex additive Gaussian observation noise vector, namely w(n) ∼ N (0, C).

Evaluation of Geometric Parameters
According to Figure 1, the reflecting specular and diffuse signals received by the antenna array are modeled based on a curved-earth surface assumption. Consequently, based on geometric considerations, we are able to obtain the explicit values of θ t and θ r , as well as that of the grazing angle ψ g (which are used in the following): where h r and h t are the radar and target heights, respectively. Additionally, r e denotes the effective earth's radius, which equals 4/3 times the earth's radius (≈8504 km). Furthermore, R, R 1 and R 2 represent (a) the radar to target range, (b) the radar to the surface reflection point range and (c) the reflection point to target range, respectively. Finally, r 1 and r 2 are the corresponding surface distances of R 1 and R 2 . Such distance pairs can be obtained through the following formulas: where f

Specular Component Modeling
In this study, we model the specular component as follows s r (n) = Γ D ρ s e −i∆φ s(n) (9) where Γ is the Fresnel reflection coefficient, D is the divergence coefficient, ρ s is the specular scattering coefficient, and ∆φ denotes the phase lag. We now provide details for each of these constituting parameters.
The Fresnel reflection coefficient, which determines the ratio of the signal amplitude after and before the incidence on a surface, is given for either vertical and horizontal polarization as follows [30] horizontal polarization (10) We observe that the coefficient Γ is a function of (a) the grazing angle ψ g and (b) the surface dielectric coefficient ε c . The latter parameter is expressed as ε c = ε − i60λσ, where ε is the dielectric constant and σ is the conductivity. Such dependence makes Γ a frequencydependent parameter.
Conversely, the divergence factor, accounting for the incidence of the signal on a curved earth surface, can be expressed approximately as where r 1 and r 2 are surface distances of R 1 and R 2 (cf. Equations (5)-(8)), respectively. We can extract the angles with respect to the curved earth and propagation model in Figure 1.
Furthermore, we use Northam's model [31] for accurately calculating the specular scattering coefficient, namely: where g = (σ h sin(ψ g )/λ) denotes the surface roughness coefficient, which depends on both the RMS of the surface roughness σ h and the grazing angle ψ g . Finally, the phase lag functional expression is ∆φ = −k(R − (R 1 + R 2 )). The latter quantity depends on the radar to target range (R), the radar to the surface reflection point range (R 1 ) and the reflection point to target range (R 2 ).

Diffuse Component Modeling
According to Beard's experimental data [32], the random diffuse component can be adequately modeled by considering its reflection point at the specular reflection point (i.e., θ r = θ d ).
Hence, the vector component associated to the diffuse signal is given as: where ρ d is used to indicate the (complex-valued) diffuse reflection coefficient. The latter parameter is modeled as a circularly-complex valued Gaussian random variable [31,33].
Thus, the diffuse coefficient ρ d can be equivalently expressed in the form ρ d = a d e iϕ d , where a d ∼ Rayleigh(σ d ) and ϕ d ∼ U (0, 2π), respectively. In the considered diffuse amplitude model, the Rayleigh parameter σ d assumes the following expression which depends on both the surface roughness coefficient (g) and the Fresnel reflection coefficient (Γ). To gain insight on the relative importance of specular and diffuse components in different scenarios, Figure 2 illustrates the specular scattering coefficient ρ s and the Rayleigh parameter σ d for different surface roughness conditions at an antenna height of 15 m. The considered target is assumed to be fixed at range 20 km and altitude 80 m. Specifically, in smooth surfaces (g < 0.1), the specular component has strong power, and the specular scattering coefficient reaches a unit value whereas the Rayleigh parameter of the diffuse component peak value equals 0.5, thus, showing that the specular signal has a more destructive effect. Therefore, this demonstrates that the angle error on the smooth surface is greater than the rough surface (0.1 < g < 0.5).

Elevation Angle Estimation
This section describes the proposed methodology for low-angle tracking in multipath conditions. First, in Section 3.1, the EKF filter is derived, and the problem of unknown parameters is tackled. Then, in Section 3.2, the capitalization of FD is explained, and the fusion of per-frequency estimates is integrated within the considered EKF workflow.

EKF-Based Tracking
In this section, we propose an EKF method for low-elevation-angle estimation. The (linear) KF is not valid for the low-angle target estimation because the observation relation with the elevation angle in Equation (1) is nonlinear. As mentioned earlier, we also use the predictions of the velocityθ t (n) and the accelerationθ t (n) of target elevation angle to improve the estimation accuracy. Accordingly, the state vector considered herein is The considered state-evolution and measurement models are summarized as follows where u(n) ∼ N (0, Q) denotes the system (viz. process) noise (assumed to be independent from sample to sample) and In the above definitions, the term T represents the sampling period. Conversely, h(·) denotes the N-dimensional measurement mapping, expressed in Equation (1), which relates θ t (i.e., the elevation angle of the target) to the noise-free received signal vector. Hence, the observation vector does not (explicitly) contain the velocity (θ t ) and the acceleration (θ t ) of the elevation angle. Finally, w(n) is the measurement noise vector as defined in Equation (1).
Thus, the prediction equation for obtainingθ t (n|n − 1) and the corresponding minimum prediction MSE matrix M(n|n − 1) can be expressed aŝ whereθ t (n − 1|n − 1) and M(n − 1|n − 1) are the final (elevation angle) target vector estimate at previous sample and the corresponding MSE matrix, respectively. Both these quantities are obtained based on the correction step driven by the considered measurement model. Regarding such a model, we first analyze, in what follows, the single-frequency (narrowband) case. Then, in Section 3.2, we generalize the correction step to the multifrequency case. Starting from Equation (1), we ignore at the design phase the unknown terms a(θ r ) s r (n) and a(θ d ) s d (n), and we approximate the (non-linear) function h(·) by linearizing it around the current estimate. In order to do so, we resort to the calculation of the Jacobian matrix of the mapping h(·) w.r.t. the state vector θ t (n), denoted with H(n): We recall that H(n) is an N × 3-dimensional matrix. Additionally, from inspection of h(·), the observation does not contain the velocity and the acceleration of the elevation angle. Hence, the second and third columns of H(n) are always null. From inspection of Equation (19), it is apparent that the presence of the unknown scattering parameter s(n) precludes the exact calculation of both h(n) and H(n) at each time step n. Accordingly, we propose to useĥ(n) andĤ(n), where s(n) is replaced by a corresponding suitable estimatê s(n).
Specifically, given the received array signal x(n), we can estimate the complex amplitude of the direct signal from the following least squares (LS) problem s(n) arg min s(n) x(n) − a(θ t (n)) s(n) 2 (21) Thus, the estimated complex amplitude of the direct signalŝ(n) is obtained in closed form asŝ In the above equation, the unknown AoA θ t (n) is in turn replaced with the corresponding prediction one-step predictionθ t (n|n − 1). As a result, the equations of Kalman gain matrix K(n), correction (estimation)θ t (n|n) and filtering MSE matrix M(n|n) for the complex-valued EKF are given as [34] where I 3 is the 3 × 3 identity matrix, and C is the observation noise covariance matrix.

Proposed FD-Based Smoothing
When capitalizing the diversity arising from wideband processing, multiple measurements (at different frequencies) are available for the direct AoA θ t (n) at the same time. Accordingly, the state-evolution and measurement models generalize as follows: where f = 1, . . . , F is used to denote the frequency bins, with F being the total number of frequencies considered. The above measurement model at different frequencies can be also collected in compact form as follows: where we have used the definitions x fd (n) Herein, we first observe that, with respect to the single-frequency case, the prediction step is unvaried, i.e.,θ Conversely, in the following we present two FD-based methods for improving the performance of the correction step. Specifically, one relies on simple weighting of per-frequency elevation estimates (WFD), whereas the other provides a matrix-based elaboration of the aforementioned estimates (MFD). The final elevation-angle estimate provided by the proposed FD-based methods has very high accuracy due to smoothing and the multipath mitigation effect. This is achieved due to the combined use of EKF-based tracking (which leverages the time-dependence of the AoA) and multiple-frequency capitalization.

Smoothing via Weighting-Based Frequency Diversity (WFD)
In the first approach, we arrange the estimated elevation angles for all frequencies by increasing magnitude, namelŷ After that, we multiply each of these angles (as well as its corresponding MSE matrix M f ) by an appropriate weighting coefficient and add them [4]. By doing so, we can eliminate large peak errors by always giving the lowest proportion to the largest elevation angle estimation. On the other hand, intermediate estimations have lower error probabilities. Then, weighting by rank is expected to yield better estimate performance than a simple average. The rank-based weighting formula considered in this work is the following where r 1 (resp. r F ) represents the biggest (resp. smallest) coefficient. Accordingly, we assign these coefficients and derive the final angle estimation (θ wfd t ) and averaged MSE matrix (M wfd ) as followŝ where M ( f ) , f = 1, . . . , F is the filtering MSE matrix at f th frequency (after re-ordering). For brevity, we have used the short-hand notationsθ  (28) and (29)) for the next time step.

Smoothing via Matrix-Based Frequency Diversity (MFD)
In such a case, the update step in the presence of FD is obtained by the following set of correction equations exploiting the compact representation reported in Equation (27): Analogously, the final angle estimationθ mfd t and filtering MSE matrix M mfd are used to update the prediction step (cf. Equations (28) and (29)) for the next time step. In the above formulation, we capitalized the stacked observation covariance matrix, e.g., C fd = diag(C 1 , . . . , C F ) and the stacked non-linear mapping and Jacobian matrix. The former is clearly equal to: . . .
whereas the latter is given as: Since the signal over different frequencies is unknown, a similar approximation can be performed by plugging the estimatesŝ(n, f ) in H fd (n) and h fd (n), obtained as: Capitalizing the above estimates provides the quantitiesĤ fd (n) andĥ fd (n), respectively.

Simulation Results
In this section, we evaluate the tracking performance of the proposed approach via simulations. In detail, a moving target from 5 to 20 km at 80 m altitude is considered. The target emits a constant signal during the whole observation window. The noise vector is assumed as w(n) ∼ N (0, σ 2 c I N ). Two multipath environments are considered to examine the proposed method at a satisfactory depth: (i) a smooth surface scenario (with dominance of the specular reflection) and (ii) a rough surface (with dominance of the diffuse reflection). The SNR is set to 10 dB, and five operational frequencies for the WFD process are considered, taken around the center frequency of 15 GHz. The sampling period equals T = 0.01 s. The process noise vector within EKF is modeled as u(n) ∼ N (0, σ 2 q I N ), where σ q = 0.005.
The ULA with the 15 m height is made of N = 10 elements positioned vertically and spaced by the center wavelength. The number of the snapshots used in the simulations is 10 for each frequency and the number of Monte Carlo trials is 100. More details about the antenna specifications and environment parameters are listed in Table 1. We compare the proposed FD-based methods with the double unitary transform based GMUSIC (DU-GMUSIC) and the polarization smoothing GMUSIC (PS-GMUSIC) methods [15,16]. Since both methods require a high number of snapshots to estimate the covariance matrix of received signals accurately, S = 256 snapshots are considered. In what follows, we first focus on one instance of the considered trajectory, i.e., by analyzing only a single realization of the measurement noise.
As shown in Figure 3a,b, the proposed methods are tracking very accurately the (true) target elevation angle in both smooth and rough surfaces. Indeed, both the proposed approaches suppress the diffuse fluctuations very well based on EKF smoothing ability and with the capitalization of FD. As depicted in Figure 4a,b, the proposed methods are also capable of providing explicit velocity and acceleration estimates of the elevation angle. Conversely, MUSIC-based baselines (PS-GMUSIC and DU-GMUSIC) are not designed to estimate the velocity/acceleration explicitly. For the sake of a complete comparison, we obtained their velocity/acceleration estimates as the first derivative of the angle estimates at consecutive time indices. Accordingly, the proposed methods are tracking the velocity/acceleration in a relatively-accurate fashion in both surfaces and over all the trajectory, while the other two competitor methods are providing out-of-scale velocity/acceleration estimates.  To provide a complementary (finer-grained) error analysis, Figure 5a,b reports the scatter plots for the elevation angle estimate, which obtained over 100 Monte Carlo runs of the same trajectory. As apparent from the results, the estimates for each range bin are mostly distributed around the true angle for both the proposed methods, with insignificant outlier points. Accordingly, the proposed estimators exhibit both a negligible bias and a very limited variance. As expected, proposed methods performed better over rough surfaces than over smooth ones. Indeed, the latter causes a stronger reflection power (according to Figure 2). On the contrary, both PS-GMUSIC and DU-GMUSIC present highly-biased angle estimates with a significant (higher) standard deviation.
To assess the estimation accuracy of the elevation angle estimation, we now assess it via the well-known root mean square error (RMSE). Specifically, Figure 6 shows the RMSE of the proposed method, which is much lower than (about one tenth of) the RMSE of the DU-GMUSIC and PS-GMUSIC methods in both cases (smooth and rough surfaces). Accordingly, the error of the proposed methods is negligible. However, the other two com-petitors use more snapshots than the proposed method, but they are unable to track the target with respect to their severe fluctuations in estimation caused by the multipath effect.
The performance of these methods are degraded in the very-low angle situation. Then, to stress the performance of the proposed approaches in a wide range of relevant scenarios, Figures 7 and 8 represent their performance against the SNR and the number of snapshots, respectively. First, by looking at Figure 7, the proposed method shows the excellent performance in tracking the low-angle targets. The performance of both proposed methods, even at very-low SNR, are appreciably higher (viz. lower in terms of RMSE) with respect to the other two (MUSIC-based) baselines considered.
As a result, the proposed FD-based methods can be used in very-low SNR situation based on their unique high performance peculiarity in the above condition. Moreover, the proposed methods are capable of ensuring satisfactory performance even in snapshot-limited (high-variability) scenarios, whereas MUSIC-originated methods require a considerably-higher number of snapshots, as evident from Figure 8. Furthermore, to assess the stability of the proposed trackers to non-nominal conditions, the impact of the mismatch of observation noise variance σ 2 c between measurements and EKF is assessed in Figure 9. In the above case, a mismatch factor is considered in the observation noise covariance matrix (i.e., C = σ 2 c I N ) within the EKF. As expected, Figure 9 reveals that a high value of the mismatch factor has a slight negative impact on the proposed method overall performance and that the proposed method is strongly robust against observation noise mismatch over both FD algorithms.
Interestingly, in the rough surface scenario, the proposed method exhibits a higher performance over undervalued mismatches (i.e., below 10 0 ). This mainly occurs since reducing the mismatch factor implies that the EKF is weighting more the correction step than the prediction step (cf. Equations (24) and (35)), which, for the rough surface scenario, better represents the observation noise due to random behavior of the diffuse component.
Finally, Figure 10 shows the RMSE of the target elevation angle estimate against the RMS of the surface roughness (σ h ). The performance results are reported for both the proposed approaches and the considered MUSIC-based baselines. As expected, the RMSE of the proposed methods are significantly lower in different multipath situation, especially for high values of σ h .

Conclusions and Future Directions
In this paper, we tackled the problem of low-angle isolated-target tracking in multipath environments. Specifically, we proposed an elevation angle estimation method via the EKF and promoted this method by two FD processes. Our method does not require any additional prior knowledge of geometrical information and physical parameters and, equally important, enjoys a low computational load. The simulation results showed that the proposed methods have a precise (close) elevation angle estimation in any multipath situation whereas the considered baselines (i.e., DU-GMUSIC and PS-GMUSIC methods) are not able to track the target accurately.
Indeed, due to the EKF properties and the proposed FD processes, our approaches are able to considerably suppress the diffuse and specular components, thus, mitigating their corresponding fluctuations. Both proposed methods maintain their excellent performance even in very-low SNR situation, while the performance of many other relevant methods drops off significantly. These methods can be used in the low probability of intercept or low-peak-power radars due to their high performance in very low SNR. Future works will account for multiple-target tracking in the considered setup, privacy issues associated to tracking estimates [35] and assessing the performance of the proposed approaches with real data.

Conflicts of Interest:
The authors declare no conflict of interest.