Instantaneous altitude estimation of maneuvering target in over-the-horizon radar exploiting multipath Doppler signatures

Over-the-horizon radar systems are capable of localizing targets in range and azimuth but are unable to achieve reliable altitude estimation. Past work has shown that the time-varying Doppler signatures of micro-multipath signals provide rich information for reliable estimation of altitude changes. In this paper, we develop a new technique for the estimation of instantaneous altitude of maneuvering targets by exploiting the estimated multicomponent Doppler signatures. A key contribution of this paper is to apply effective non-stationary signal analysis techniques for estimating the time-varying Doppler signature of each individual multipath. This enables reliable target altitude estimation in an extended Kalman filter setup. The maximum a-posteriori criterion is used for the estimation of initial target states.

High-frequency (HF) over-the-horizon radar (OTHR) systems that exploit skywave propagation, i.e., reflection and refraction of radar signals from the ionosphere, provide wide area surveillance capabilities to detect and track targets at farther range (thousands of kilometers) [1, 2]. The capability of OTHR systems to cover a surveillance area beyond the range of conventional line-of-sight radars makes them uniquely important in a number of applications.

In OTHR operations, narrowband signals are used because the available signal bandwidth is constrained by the prevalent ionospheric conditions and the range extents. As such, the range resolution of an OTHR system is typically measured in tens of kilometers [2]. Traditional OTHR systems use a two-band linear array for transmit and a very long array on receive, which provide large apertures for high azimuth resolution. Recent developments in OTHR systems suggest the use of two-dimensional (2-D) transmit and receive array configurations, combined with advanced signal processing techniques, such as multiple-input multiple-output (MIMO) radar approaches, to effectively mitigate noise and spread-Doppler clutter sensitivity problems [3]. In these configurations, the aperture in the range direction provides observation and processing in the elevation dimension to support target altitude tracking.

The limitations in signal bandwidth and array aperture, as well as the low signal-to-noise ratio (SNR) make it difficult for an OTHR system to provide good range and cross-range resolutions. Because radar signals are reflected by multiple ionospheric layers, a number of techniques have been developed over the years to exploit the rich multipath propagation associated with moving targets to improve the target localization and tracking performance [4–6]. Multipath exploitation radar is also found in other applications, such as through-the-wall and urban terrain sensing [7, 8].

Practically, the range and azimuth resolutions may not necessarily require a very high resolution. This is not the case for target altitude. Target altitude provides important information for the determination of the type and flying course of a target. For this information to be useful, an accuracy of a few kilometers or less, rather than in the order of tens of kilometers, would be desirable. Toward this end, effective techniques based on micro-multipath model can be devised to provide accurate altitude information. These models make use of multipath returns due to ocean or ground reflections that are local to the target [4]. The micro-multipath propagation is illustrated in Figure 1. In essence, instead of using multipath propagation through different ionospheric layers, the micro-multipath propagation model focuses on the relatively stable E layer and earth surface reflection that is local to the target. For convenience of analysis and description, we use, in this paper, the simplified flat-earth ionosphere model, as depicted in Figure 2[9]. The approach presented for the flat-earth model can nevertheless be implemented in an actual propagation environment after the spheric earth model is properly considered.

Micro-multipath propagation in OTHR.

Full size image

Flat-earth model.

Full size image

In the flat-earth model, the reflected propagation paths can be equivalently considered as straight lines to their mirrored positions, where the top line corresponds to the micro-multipath (reflected at the ionosphere and the earth surface), and the second line from the top corresponds to the direct path (only reflected at the ionosphere). Assume that in addition to the target motion in the range direction, which generates a nominal Doppler signature, the target ascends in altitude. This occurs during the target’s departure and landing, as well as in its flying course. In this case, the direct path becomes longer, while the multipath becomes shorter. These paths behave conversely when the target descends. As such, changes in the target altitude will alter the distance of both direct path and multipath, generating Doppler variations with opposite signs around the nominal Doppler signature [4]. These propagation paths of the emitted/received signals associate themselves with different non-linear time-frequency trajectories, each corresponding to a Doppler signature of the target along a propagation path. For maneuvering targets, high-resolution time-frequency analyses have been shown to be effective in resolving the multicomponent Doppler signatures, and thus revealing rich and important information about the relative target altitude [9, 10]. Recently, new approaches for accurate parametric estimations of time-varying multicomponent signals with closely separated Doppler signatures encountered in OTHR systems have been developed and have yielded high accuracy estimation of the relative target altitude [11, 12]. We maintain, however, that the estimation of the actual instantaneous target altitude has not been considered within the non-stationary signal analysis framework.

In this paper, we develop a robust altitude estimation technique for maneuvering targets in an MIMO radar environment. Specifically, we focus on the array apertures in the range direction that allow for spatial processing in the elevation dimension. The proposed technique is based on recent advances in non-stationary signal analyses for instantaneous multicomponent Doppler signature estimations. Because Doppler information alone does not provide sufficient information of the absolute target altitude as well as the elevation maneuvering direction (ascending or descending), the target positions are estimated using the extended Kalman filer that exploits different hypotheses of the initial conditions. We then utilize the maximum a-posteriori (MAP) criterion to estimate the initial target altitude and the maneuvering direction. Reliable target altitude estimation is enabled using the initial estimate and the knowledge guided by the time-varying Doppler signature of each individual multipath.

The following notations are used in this paper. A lower (upper) case bold letter denotes a vector (matrix). $\mathbb{E}[·]$ represents the statistical mean operation. (·)^∗, (·)^T, and (·)^H respectively denote complex conjugation, transpose, and conjugate transpose (Hermitian) operations. $\dot{a}$ denotes the derivative of a with respect to time. I _N expresses the N×N identity matrix. In addition, $\mathbb{R}$ denotes the complete set of real scalars, whereas ${ℂ}^{N\times M}$ denotes the complete set of N×M complex entries.

2.1 OTHR micro-multipath model

Consider an OTHR system as illustrated in Figure 2. Our main interest lies in the target altitude estimation. Therefore, for simplicity and without loss of generality, we only consider 2-D position and velocity (range direction x and altitude direction z) in Cartesian earth-centered earth-fixed coordinates. The target states at slow-time index t and is then described by the 2-D position and velocity as:

Correspondingly, only the array apertures that lie in the range direction are considered, and the cross-range array apertures are ignored. As such, both transmit and receive arrays are considered to be linear, and their apertures extend along the x-axis.

Due to the presence of micro-multipath propagation, the combination of the direct path and the multipath in both forward and return links yields the following four combinations of the two-way propagation: path I (l ₁:l ₁), path II (l ₂:l ₂), path III (l ₁:l ₂), and path IV(l ₂:l ₁). Among them, paths III and IV yield virtually identical two-way slant range.

To further describe the advantages of exploiting micro-multipath information, we can represent the multipath propagation in an alternative view, as shown in Figure 3. The micro-multipaths result in two sets of virtual sensors that are respectively located on the top and bottom in this figure. The very large separation of these virtual sensors achieves a large virtual array aperture in the elevation dimension and thus permits accurate target altitude estimations.

Extended array aperture due to micro-multipath propagation.

Full size image

2.2 MIMO signal model

Consider a monostatic MIMO radar system consisting of N _t closely spaced transmit antennas and N _r closely spaced receive antennas. Denote $\mathbf{\text{S}}\in {ℂ}^{N_t}\times T$ as the narrowband waveform matrix which contains orthogonal waveforms to be transmitted from the N _t antennas over a pulse repetition period of T fast-time samples. We assume that the waveform orthogonality is achieved in the fast-time domain, i.e., by denoting s _i as the i th row of matrix S, s _i and s _m are orthogonal for any i≠m with different delays, and s _i is orthogonal to its delayed versions. We also assume that s _i has a unit norm, i.e., $\mathbf{\text{S}}{\mathbf{\text{S}}}^H={\mathbf{I}}_{N_t}$.

Thus, for the k th slow-time pulse, the signal received at the m th receive array element is expressed as:

where I=4 is the number of multipath components, ρ _i,k is the complex response of the i th multipath that accounts for target reflection and propagation delays and will be discussed later in Section 3. Superscripts [D] and [A] respectively denote the departure and the arrival modes. The phase term ${\varphi }_{i,m,k}^{\left[A\right]}$ in (2) is determined by the m th receive antenna location, $d_m^{\left[A\right]}$, and is expressed as:

where ${\theta }_{i,k}^{\left[A\right]}$ is the elevation angle of arrival of the i th path. Similarly, ${\varphi }_{i,n,k}^{\left[D\right]}$ is determined by the n th transmit array sensor location, $d_n^{\left[D\right]}$, and is expressed as:

where ${\theta }_{i,k}^{\left[D\right]}$ is the elevation angle of departure of the i th path. To simplify processing, we assume that the clutter is sufficiently removed through, e.g., notch filtering of the low-frequency components around the direct current (DC) region. As such, w _m,k in (2) at the m th receive element can be considered as additive noise, which is independent and identically distributed (i.i.d.) white complex Gaussian $\mathcal{C}\mathcal{N}(0,{\sigma }_n^2)$ and is independent of each other and of the target returns.

For path I, the departure and arrival elevation angles of the target are identical, i.e.,

Similarly, for path II, the departure and arrival elevation angles of the target share the following expression:

For path III, the elevation angle of the departure path is ${\theta }_{3,k}^{\left[D\right]}={\theta }_{1,k}$ and that of the return path is ${\theta }_{3,k}^{\left[A\right]}={\theta }_{2,k}$. For path IV, the elevation angles corresponding to the departure and return paths are ${\theta }_{4,k}^{\left[D\right]}={\theta }_{2,k}$ and ${\theta }_{4,k}^{\left[A\right]}={\theta }_{1,k}$, respectively.

Multiplying u _m,k by ${\mathbf{\text{s}}}_n^H$, we obtain:

where $w_{m,n,k}={\mathbf{w}}_{m,k}{\mathbf{\text{s}}}_n^H$.

Denote ${\mathbf{\text{a}}}_{i,k}={[e^{j({\varphi }_{i,1,k}^{\left[A\right]}+{\varphi }_{i,1,k}^{\left[D\right]})},\mathrm{...},e^{j({\varphi }_{i,N_r,k}^{\left[A\right]}+{\varphi }_{i,N_t,k}^{\left[D\right]})}]}^T\in {ℂ}^{N\times 1}$ as the steering vector of the virtual array corresponding to the i th path, where N=N _t N _r , and let A _k =[a _1,k ,...,a _I,k ]. Further, denote ρ _k = [ρ _1,k ,...,ρ _I,k ]^T, and ${\mathbf{w}}_k={[w_{1,1,k},\mathrm{...},w_{N_t,N_r,k}]}^T$. Then, the measurement data vector is expressed as:

3.1 Doppler frequency model

Now, we consider the Doppler frequency contained in ρ _i,k in (7). Toward this end, we express ρ _i,k as

for i=1,...,I, where λ=c/f _c is the wavelength corresponding to carrier frequency f _c , and σ _i,k represents the propagation channel coefficients representing the combined effect of transmit power, target reflection, which is a function of radar cross section, and the ionosphere as well as surface reflections of the i th path at time index k. In addition, η _i,k is the two-way slant range. For the ease of presentation, we approximate the actual propagation environment with a flat-earth model as illustrated in Figure 2. In practice, the transmit and receive arrays are separated by a distance, which is much smaller than the range, to mitigate power leakage. We assume that the difference of the range due to the array positions is negligible. Thus, η _i,k is related to the multipath lengths as:

where $l_{u_{T,i},k}$ and $l_{u_{R,i},k}$, respectively denote the slant range of the forward and return links associated with path i at time instant k. In particular,

The slant ranges l _1,k and l _2,k can be expressed in terms of the ground range x _k , the ionosphere altitude H, and the target altitude z _k , as

For convenience of computation, we use the following two variables as our observations:

Equations 12 and 13 are used for numerical evaluations of the slant ranges. Nevertheless, for illustrative purposes, we also use the following approximations, which hold true for x _k ≫H≫z _k , to reveal a clear relationship between the slant ranges and the other parameters:

The average and difference Doppler frequencies are then obtained as:

where $K_k=1-2H^2/x_k^2$.

From the above discussion, it is evident that, while the dominant Doppler component $f_{\text{ave},k}=2K_k{\dot{x}}_k/\lambda$ is shared by all the four paths and reveals the target velocity in the range direction, the small Doppler difference between the paths, $f_{\text{diff},k}=4H{Ż}_k/\left(x_k\lambda \right)$, is a function of ${Ż}_k$. Effective time-frequency analysis allows separation of the multicomponent Doppler signatures based on the Doppler frequency difference [11]. In this paper, the resolved Doppler signatures are used for improved target altitude tracking.

In practice, there is ambiguity in the sign of the estimated Doppler difference, i.e., the Doppler frequency difference f _diff,k which by itself does not reveal whether a target is ascending and descending. This ambiguity will be considered and resolved in the target tracking process.

3.2 Instantaneous frequency estimation

Similar to many other narrowband radar systems, an OTHR system heavily relies on Doppler analysis to separate targets from clutter and to reveal many important information. Because of the low SNR, OTHR systems often utilize a long coherent integration time (CIT) to achieve a high processing gain. The use of chirplet transform allows a longer CIT for weak target detection when the target returns can be modeled as linear frequency modulated signals [13]. In this paper, we use the signal stationarization technique to achieve further enhancement of non-stationary signals that demonstrate high-order time-varying characteristics [10, 11]. In the following, the instantaneous frequency (IF) estimation and signal stationarization techniques are briefly summarized. The proposed technique is based on the local analysis of time, frequency, and phase coherence, and uses this information to merge local components in order to estimate the global time-frequency structures characterizing the signal.

The primary challenges of this problem lie in the difficulty of separating and resolving the multicomponent returns that are very close in their time-varying Doppler signatures as well as their spatial signatures. We use a two-step procedure for an improved IF estimation, i.e., coarse signal stationarization and fine IF estimations.

In the first step of our approach, we assume that the multicomponent signals are difficult to separate and only a coarse estimation of their direction, ${\widehat{\theta }}_{\text{coarse}}$, can be obtained by treating the received signal as a single-component signal. For this purpose, ${\widehat{\theta }}_{\text{coarse}}$ is considered invariant for the time period being considered because of its small variation. The corresponding steering vector is denoted by ${\mathbf{\text{a}}}_0=\mathbf{\text{a}}\left({\widehat{\theta }}_{\text{coarse}}\right)$. Capon beamforming is applied to the array output to form a single signal stream with array gain achieved:

where R is the expected covariance matrix of r _k , i.e., $\mathbf{R}={\mathbb{E}}_k\left[{\mathbf{r}}_k{\mathbf{r}}_k^H\right]$. In practice, R is obtained through sample averaging.

The next step is to obtain the IF estimates of the time-varying Doppler signatures from the combined signal ${\tilde{r}}_k$. The known information that the three distinct Doppler signatures are closely separated allows us to stationarize all the Doppler components with the IF estimate of a single signal component. The process is illustrated in Figure 4. Denote the resulting Doppler signature estimate of one component (which can be any of the three distinct components) as f _D,coarse and the corresponding phase signature as ψ _coarse. The procedure of estimating the IF estimation of a single-component signal is summarized in the Appendix. By multiplying ${\tilde{r}}_k$ by exp(−j ψ _coarse), all the three components ${\tilde{r}}_{k,\text{coarse}}={\tilde{r}}_{k}exp(-j{\psi }_{\text{coarse}})$ become approximately stationarized and are located in the low frequency region (Figure 4b). With the waveform being more stationary, the IF estimation procedure can be performed again to obtain finer IF and phase estimates (denoted as f _D,fine and ψ _fine, respectively). By multiplying ${\tilde{r}}_{k,\text{coarse}}$ by exp(−j ψ _fine), improved signal stationarization is achieved (Figure 4c). As a result, one of the signal components is stationarized with a high accuracy, i.e., it approximately becomes a DC component, whereas the Doppler signatures of the other two components have non-zero frequencies. The original Doppler signature of the stationarized component can be determined from the IF estimates used in the stationarization, i.e., f _D,coarse+f _D,fine.

Illustration of the IF estimation process. (a) Original Doppler signature, (b) after coarse stationarization, (c) Stationarization with accurate Doppler estimate, (d) after removal of the stationarized component, (e) stationarized with one of the remaining components, and (f) after removal of the stationarized component.

Full size image

Now, we can filter out the stationarized signal component by removing the DC component, yielding the two remaining components (Figure 4d). We then repeat the IF estimation and stationarization process to obtain the IF estimate of the second signal (Figure 4e) and remove the stationarized signal component to estimate the last signal component (Figure 4f). After all the three components are estimated, they are sorted such that f _D,1,k and f _D,2,k have the highest and the lowest values, respectively, whereas f _D,1,k takes value between them.

3.3 Signal filtering

With the separation of the IF signatures, each multipath signal component can be separated as well. In particular, we are interested in signal components corresponding to paths I and II, which allow us to analyze the target maneuvering, in terms of the Doppler frequency and elevation angle, for each individual path.

By multiplying the vector signal r _k by the conjugate of the phase estimate of the i th path, exp(−j ψ _i,k ) for i=1 and 2, the signal associated with the i th path is concentrated around the DC component. We can design a narrowband filter to keep the i th component and to filter out the other multipath components. Note that unlike in the direction-of-arrival estimation approaches [14] where only a single zero-frequency bin is used, we need to capture multiple frequency bins of the stationarized i th component so as to keep the time-varying Doppler information, which is particularly important to determine the direction of the elevation velocity of the target. The captured signal vector is multiplied by exp(j ψ _i,k ) to restore the original Doppler information in each signal component., i.e.,

where $\mathcal{P}$ denotes the filtering processing. In this paper, the filtering is implemented in the frequency domain by masking the Fourier transform coefficients, i.e., $\mathcal{P}(·)={\mathcal{F}}^{-1}\mathcal{M}\mathcal{F}(·)$, where $\mathcal{F}$ denotes the Fourier transform, and $\mathcal{M}$ is a proper binary mask with ones in the passband around the DC component.

By eliminating the effect of interactions between different paths, ${\mathbf{r}}_k^{\left[i\right]}$ enables better association of the measurement data with the target maneuvering.

Naturally, the array data vector, r _k , is considered as the observation vector. Target geolocation based on extended Kalman filtering exploiting the array data vector, however, does not yield satisfactory target altitude information. To take advantage of the resolved estimates of the Doppler signatures, the following steps are performed to yield high-accuracy target altitude estimates.

• The average and difference Doppler signatures, denoted as a vector f _D,k =[f _ave,k ,f _diff,k ]^T, are used as additional observations. The incorporation of the instantaneous Doppler estimates generally provides good estimation of relative target altitude, but the instantaneous target altitude is still very sensitive to the initial target position vector ${\tilde{\mathbf{x}}}_0$, particularly the initial altitude z ₀, assumed in time k=0. In addition, as we discussed before, the Doppler difference has an ambiguity in the target direction of its elevation maneuvering.

• To overcome these problems, we use multiple hypotheses of the initial target position and vertical orientation and find the best solution that maximizes the MAP criterion. Note that due to the low SNR involved in this problem, the a-posteriori probability offered by the measured data at each time instant is not reliable enough to provide meaningful information. Rather, we use the a-posteriori probability of all the observed time instants so that a reliable MAP metric is achieved.

In the following, we introduce the target state model, and the estimation for instantaneous target estimation, assuming an initial target position, is described. The MAP-based estimation of the initial target altitude and the elevation orientation are then addressed.

4.1 Target state model

The target state vector, x _k , consists of the target position vector, ${\tilde{\mathbf{x}}}_k$, and the complex response vector σ _k =[σ _1,k ,σ _2,k ,σ _3,k ]^T. Note that σ _4,k is excluded from σ _k because it equals to σ _3,k due to propagation reciprocity. Therefore, we have:

Then, the target state evolves according to the following linear stochastic model:

where

and Δ is the pulse repetition interval.

The process noise v _k is assumed to be zero-mean white Gaussian noise with covariance matrix Q. In this paper, we define Q to have the form:

where ${\sigma }_p^2$, ${\sigma }_v^2$, and ${\sigma }_g^2$ are, respectively, the process noise variance for the position, velocity, and mode power level components of the state vector. Under these assumptions, the state transition probability density function (pdf) has the following form:

4.2 Instantaneous target altitude estimation

In this step, the target altitude will be estimated with hypotheses of the initial target position, x ₀, and the direction of the target elevation velocity, i.e., ascending or descending. The ambiguities will be solved in the following subsection.

In addition to the array data, the resolved Doppler estimates, which are related to the target states through (13), can be used as additional constraints. Because of the ambiguity of target vertical orientation as represented by the absolute value operation in (13), we divide the problem into two cases, i.e., when the target ascends or when it descends. For the former,

whereas when the target descends, we have:

Such relationships can be treated as additional constraints or penalty functions to improve tracking performance [6, 15]. In this paper, however, we simply augment the Doppler estimates, f _D,k =[f _ave,k ,f _diff,k ]^T, to be part of the observation vector of the maneuvering target at time instant k. As a result, the observation vector is defined as ${\mathbf{z}}_k={\left[{\left({\mathbf{r}}_k^{\left[1\right]}\right)}^T{\left({\mathbf{r}}_k^{\left[2\right]}\right)}^T{\mathbf{f}}_{D,k}^T\right]}^T$, which contains the array data individually separated for the two multipath modes and the Doppler signatures. The corresponding observation equation is expressed as:

where h _k (x _k ) collectively define the non-linear relationship between z _k and x _k as described in different parts of Section 3, and ${\check{\mathbf{w}}}_k$ represents the observation noise which is assumed to be zero mean and with covariance matrix R _n . More specifically, h _k (x _k ) is divided into four components, i.e., ${\mathbf{h}}_k\left({\mathbf{x}}_k\right)={\left[{\mathbf{h}}_{1,k}^T\right({\mathbf{x}}_k\left){\mathbf{h}}_{2,k}^T\right({\mathbf{x}}_k\left)h_{3,k}\right({\mathbf{x}}_k\left)h_{4,k}\right({\mathbf{x}}_k\left)\right]}^T$. These four components are respectively expressed as:

4.3 MAP-based initial state estimation

We attempt to solve the aforementioned two problems, i.e., the initial target altitude estimation and target elevation movement direction (ascending or descending), by making multiple hypotheses of the initial target position and vertical orientation and find the best solution that maximizes the MAP criterion.

Define

to be the collection of the state and observation vectors, and

to be the collection of the estimated state under the assumption of initial target position x ₀ and vertical orientation ν, where ν∈[1,2] with 1 denoting target descending and 2 target ascending. In MAP estimation, the goal is to find $\widehat{\mathbf{X}}$ that maximizes the posterior density, f(X|Z). Based on the Bayesian theorem, the $\widehat{\mathbf{X}}$ that maximizes f(X|Z) also maximizes lnf(Z , X)= lnf(Z|X)+ lnf(X). With the initial target altitude and moving direction in mind, we have

where q ₀(x ₀) is the pdf of the initial target altitude. On the other hand, $lnf\left(\mathbf{Z}\right|{\mathbf{X}}^{[{\mathbf{x}}_0,\nu ]})$ is given by

Our objective is to find x ₀ and ν such that $lnf\left({\mathbf{X}}^{[{\mathbf{x}}_0,\nu ]}\right|\mathbf{Z})$ is maximized.

The MAP-based initial state estimation can be expressed as:

As an example, we consider a maneuvering target which makes a 360^o circular turn of radius 5 km in approximately T=179.5 s to descend by approximately 2,250 meters. At a waveform repetition frequency of f _s =40 Hz, the total number of data sampled is K=7,180. The centers of the transmit and receive arrays are set as the coordinate origins. The other key simulation parameters are listed in Table 1. All the multipath signals are considered to fall within the same range cell. The target maintains a horizontal velocity of 175 m/s (630 km/hr). Its elevation velocity varies sinusoidally, and the corresponding target altitude is expressed as:

To account for the Doppler shift due to ionosphere, we use a general Gaussian scatter model proposed in the ITU [16] for the simulation of HF ionospheric channels. This model is based on Watterson’s channel model [17]. The Matlab code and parameters provided in [18] are used to generate the propagation channel coefficients, σ _i , where i=1,2, and 3. The correlation between σ ₃ and the other two coefficients is assumed to be 0.64 and that between σ ₁ and σ ₂ is assumed to be 0.49. The yielding Doppler signatures are depicted in Figure 5.

Three distinct Doppler signatures due to micro-multipath propagation.

Full size image

We assume an OTHR system that uses 2-D arrays in the cross-range and range directions. For target altitude estimation purpose, we are only concerned with the array aperture in the range direction. As such, by assuming that the signals across the cross-range dimension are coherently combined, we consider a small-size MIMO array that consists of six transmit antennas and ten receive antennas. Both arrays are linear, and the antennas are extended in the range direction. The minimum redundant array configurations [19, 20] are used at the transmit and receive arrays, with the unit separation being one wavelength. The transmit and receive arrays are separated by a 100-km cross-range distance. The typical ground range is in the order of thousands of kilometers. Therefore, the radar can be considered monostatic. In particular, the elevation angles are virtually identical for both arrays. Notice that because of the multipath propagation, the co-array equivalence between the transmit and receive arrays is only achieved at each resolved path (i.e., path I and path II).

Figure 6 shows the spectrogram of the Capon beamformer output (${\tilde{r}}_k$), where the input SNR at each antenna is -10 dB, and a Hamming window of length 255 is applied when performing the short-time Fourier transform. Signals falling in between -1 and 1 Hz are filtered out for clutter removal. Figure 7a shows the estimated IF signatures of the three Doppler components. Figure 7b shows the estimation frequency error for the three components. The estimation error is in the order of a few hundredth of a hertz. It is relatively larger in the two edges due to the lack of paring segment and around the two instants where the data is suppressed due to clutter removal.

Spectrogram of the received signal at the Capon beamformer output.

Full size image

Doppler frequency estimation results. (a) Estimated Doppler frequency. (b) Estimated Doppler error.

Full size image

The estimated target positions are depicted in Figure 8. In determining the target initial states, we assume that the initial target altitude is uniformly distributed between 5,000 and 15,000 m, and a step size of 50 m is used. Because the initial target altitude and the elevation motion direction are coupled, these two variables are jointly determined. As a result, the radar system can correctly identify that the target flies in a descending mode. In 100 Monte Carlo trials with different noise realizations, the motion directions are all correctly identified. The average value of the estimated initial altitude is 9,825 m, and the standard deviation is 219 m.

Tracking performance. (a) Target ground range. (b) Target altitude. (c) Range-direction velocity. (d) Elevation velocity.

Full size image

Figure 8a shows the ground range estimate (x _k ), whereas the corresponding range velocity (${\dot{x}}_k$) is plotted in Figure 8c. The initial target altitude estimate in this example is 9,850 m. Very high accuracy estimation is achieved for both parameters largely because of the small Doppler estimation error relative to the overall Doppler frequencies. Figure 8b shows the target altitude (z _k ), and the corresponding elevation velocity (${Ż}_k$) is shown in Figure 8d. Overall, good instantaneous target altitude estimation is achieved.

In this paper, we have developed a novel technique that provides an accurate estimation of the instantaneous altitude of a maneuvering target in an over-the-horizon radar system. The proposed method utilizes the micro-multipath Doppler signatures due to earth surface reflection, which are analyzed and separated by advanced non-stationary signal analysis techniques. Accurate estimation of the instantaneous frequency of such Doppler signatures allows high-quality target trajectory estimation through an extended Kalman filter with hypotheses of target elevation motion direction and initial target altitude. These hypotheses are then estimated based on the maximum a-posteriori criterion. Simulation results verified the high-accuracy target position estimation results.

IF estimation of single-component polynomial phase signal

We first consider a single-component polynomial phase signal (PPS) that is characterized by the IF:

where ϕ(t) is continuous on [0,T]. A number of techniques have been developed for the estimation of the polynomial phase coefficients, $\left\{a_m^{\left[w\right]}\right\},m=0,\mathrm{...},M$ (see e.g., [21, 22] and references therein). Below, we briefly summarize the concept of IF estimation based on multi-lag high-order ambiguity function (mlHAF) [22], which was developed based on the HAF, or polynomial phase transform, concept presented in [21]. Define the M th-order multi-lag high-order instantaneous moment (mlHIM) of signal s(t) as:

where τ _M is the set of lags, and ${\text{mlHIM}}_1\left[s\right(t\left)\right]=s\left(t\right)$. The corresponding mlHAF of s(t) is the Fourier transform of HIM[s(t);τ _M−1], i.e.,

It can be shown that |mlHAF[s(t);ω,τ _M ]| has a global maximum at:

Therefore, an estimate of a _M can be obtained from the peak position of |mlHAF[s(t);ω,τ _M ]|. Then, the M th-order polynomial phase component can be removed by multiplying the original signal s(t) with the conjugate of the estimated M th order polynomial phase signal as $s^{[M-1]}\left(t\right)=s\left(t\right)exp(-j{â}_M{t}^M)$, where ${â}_M$ denote the estimate of a _M . This procedure can be repeated by estimating a _m−1 from s ^[m−1](t) for m=M,...,2.

Now, we consider a single-component non-stationary signal which is characterized by its IF, but the IF law is rather complicated and is difficult to be represented by a PPS with a reasonable polynomial order. We divide the signal into W half-overlapped adjacent segments, each with data length T, such that the segmented signal can be considered as a PPS with a low polynomial order M (We use M=3 as the highest order of the PPS in each segment). In the w th segment, w=1,...,W, we consider that the signal is represented as a PPS of order M, i.e.,

where $t_0^{\left[w\right]}=(w-1)T/2$ is the starting time of segment w, and T is assumed to be even. The polynomial phase coefficients, $a_m^{\left[w\right]}$, can be determined using the procedure described above.

To improve the robustness of IF estimation against noise and other possible perturbation sources, multiple estimates of the polynomial phase coefficient sets are obtained for each segment. Denote $a_m^{[w,n_c^{\left[m\right]}]}$, $n_c^{\left[w\right]}=1,\mathrm{...},N_c$, as the N _c coefficient sets estimated in the w th segment. Because of the half overlapping between neighboring segments, the second half of the w th segment and the first half of the (w+1)th segment reflect the same portion of the signal and thus, ideally, should be identical. The coefficient set is chosen at each segment to minimize the following penalty function

Finally, the IF at an arbitrary time instant t ^′, w T−T/2≤t ^′≤w T, is expressed as:

where

is a triangular window function that uses a higher weight towards the center of each segment, and $t_0^{\left[w\right]}=(w-1)T$ is the start time of the w th segment. Note that h(τ) takes a value of one in edge segments when a pairing counterpart is not available.

The work of Y. D. Zhang and M. G. Amin was supported in part by a subcontract with Dynetics, Inc. for the research sponsored by the Air Force Research Laboratory (AFRL) under Contract FA8650-08-D-1303. Y. D. Zhang was also supported in part by the Air Force Office of Scientific Research (AFOSR) through the Air Force Summer Faculty Fellowship Program under contract number FA9550-09-C-0114.

Authors and Affiliations

1. Wireless Communications and Positioning Laboratory, Center for Advanced Communications, Villanova University, Villanova, PA, 19085, USA

   Yimin D. Zhang & Moeness G. Amin

2. Electrical and Computer Engineering Department, University of Denver, Denver, CO, 80210, USA

   Jun Jason Zhang

3. Air Force Research Laboratory, AFRL/RYMD, Dayton, OH, 45433, USA

   Braham Himed

Corresponding author

Correspondence to Yimin D. Zhang.

Competing interests

The authors declare that they have no competing interests.

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions