Bispectrum- and Bicoherence-Based Discriminative Features Used for Classification of Radar Targets and Atmospheric Formations

Abstract

This chapter is dedicated to bispectrum-based signal processing in the surveillance radar applications. Detection, recognition, and classification of the targets by surveillance radars have various applications including security, military intelligence, battlefield purposes, boundary protection, as well as weather forecast. One of the particular and effective discriminative features commonly exploited in modern radar automatic target recognition (ATR) systems is the micro-Doppler (m-D) contributions extracted from joint time-frequency (TF) distribution. However, a common drawback of the energy-based strategy lies in the impossibility to retrieve additional particular information related to frequency-coupling and phase-coupling contributions containing in the radar backscattering. Phase coupling contains additional discriminative features related to individual target properties. Bispectrum-based strategy allows retrieving a phase-coupled data containing unique discriminative features related to individual target properties. Bispectrum tends to zero for a stationary zero-mean additive white Gaussian noise (AWGN), providing smoothing of AWGN in TF distributions. Hence, bispectrum-based approach allows improving extraction of robust discriminative features for ATR radar systems.

1. Introduction

Bispectral density or third-order cumulant spectrum, in opposite to the common energy spectrum, allows not only describing the statistical properties of an observed process more correctly and profoundly but also extracting novel information features such as the contributions caused by spectral component correlation relationships. Moreover, bispectrum estimate allows to extract the phase relationships existing between the spectral components contained in the observed process. Therefore, the main benefit of bispectrum in comparison with energy spectrum is in the preservation of phase information contained in a process under study and possibility to recover this important information.

Triple autocorrelation function (TAF) R_x(k, l) of a discrete process {x^(m)(i)}, m = 1, 2,…, M is related to the class of third-order statistics [1]. TAF is a function of two-shift variables, and it can be written in the discrete form as

$$R_x\left(k,l\right)={⟨\sum _{i=0}^{I-1}\left[x^{\left(m\right)}\left(i\right)-E\right]\left[x^{\left(m\right)}\left(i+k\right)-E\right]\left[x^{\left(m\right)}\left(i+l\right)-E\right]⟩}_{M\to \infty },$$ E1

where E is the mean value; k = −I + 1,…, I − 1 and l = −I + 1,…, I − 1 are the shift indices.

According to the definitions given in [1], bispectrum is the 2-D direct Fourier transform (DFT) of TAF. Unlike common real-valued spectral density, bispectrum is the complex-valued function ${\dot{B}}_x\left(p,q\right)$ of two independent frequency variables p and q defined as

$${\dot{B}}_x\left(p,q\right)=\sum _{k=-I+1}^{I-1}\sum _{l=-I+1}^{I-1}{R}_x\left(k,l\right)\text{exp}\left[-j2\pi \left(\mathit{kp}+\mathit{lq}\right)\right],$$ E2

or given in the form of

$${\dot{B}}_x\left(p,q\right)={⟨{\dot{X}}^{\left(m\right)}\left(p\right){\dot{X}}^{\left(m\right)}\left(q\right){\dot{X}}^{\ast \left(m\right)}\left(p+q\right)⟩}_{\infty },$$ E3

where ${\dot{X}}^{\left(m\right)}\left(p\right)$ is the DFT of mth realization {x^(m)(i)}, m = 1, 2,…, M; * denotes complex conjugation; ${\dot{B}}_x\left(p,q\right)=\left|{\dot{B}}_x\left(p,q\right)\right|\text{exp}\left[j{\gamma }_x\left(p,q\right)\right]$ ; $\left|{\dot{B}}_x\left(p,q\right)\right|$ and γ_x(p, q) are the magnitude bispectrum (bimagnitude) and phase bispectrum (biphase), respectively; p = −I + 1, …, I − 1 and q = −I + 1, …, I − 1 are the frequency indices.

According to (3), bispectral density ${\dot{B}}_x\left(p,q\right)$ is the ensemble averaging of triple product of three complex-valued functions related to three different frequencies: p, q, and p + q. It means that bispectrum can be defined as the correlation measure between the frequencies p, q, and p + q contained in the process. The main properties of bispectrum [1] are listed below.

1. Bispectrum tends to be zero for stationary zero-mean Gaussian process. Note, however, that for a process with non-symmetrical probability density function (PDF), its bispectrum differs from zero. This property allows to extract non-Gaussian contributions against the background of Gaussian noise.

2. TAF and bispectrum are equal to zero for deterministic signals having zero asymmetry. For example, TAF and bispectrum are of zero values for a simple signal x(i) = A0 cos(2πfi). However, when some little non-linear distortions appear or when a constant component is added to a harmonic signal, TAF and bispectrum become non-zero valued functions. The latter property can serve as a very sensitive tool to detect and evaluate non-linear distortions contained in the observed signal.

3. Bispectrum is the periodical function with period equal to 2π

   $${\dot{B}}_x\left(p,q\right)=\dot{B}\left(p+2\pi ,q+2\pi \right).$$ E4

4. Bispectral density has the following symmetry property

   $$\begin{array}{ll}\dot{B}\left(p,q\right)& =\dot{B}\left(q,p\right)={\dot{B}}^{\ast }\left(-p,-q\right)={\dot{B}}^{\ast }\left(-q,-p\right)\\ & =\dot{B}\left(-p-q,q\right)=\dot{B}\left(p,-p-q\right)=\dot{B}\left(-p-q,p\right)=\dot{B}\left(q,-p-q\right).\end{array}$$ E5

   By using (5), bispectrum can be defined by a hexagon given in the bifrequency plane. It follows from (5) that bispectrum of a real-valued process can be defined completely just only within the area of the main triangular domain given by the inequalities

   $$q\ge 0,p\ge q,p+q\le I–1.$$ E6

   It is sufficient to use the symmetry relationships (5) taking into account (6) to compute bispectrum function within all other parts of hexagonal domain. Note that the conditions (6) limiting the total number of bispectrum samples allow to decrease essentially the memory requirements and, hence, to reduce the computational time of bispectrum.

5. In many practical cases, it is necessary to detect and estimate the phase coupling contributions. Since the phase relationships are lost in power spectral density estimator, it is impossible to extract information about phase coupling from common power spectral density. Meanwhile, it is possible to do that by using bispectrum

   This property can be explained by considering the following two processes

   $$x_1\left(i\right)=\text{cos}\left(2\pi f_{1}i+{\phi }_1\right)+\text{cos}\left(2\pi f_{2}i+{\phi }_2\right)+\text{cos}\left(2\pi f_{3}i+{\phi }_3\right),$$ E7a
   $$x_2\left(i\right)=\text{cos}\left(2\pi f_{1}i+{\phi }_1\right)+\text{cos}\left(2\pi f_{2}i+{\phi }_2\right)+\text{cos}\left[2\pi f_{3}i+\left({\phi }_1+{\phi }_2\right)\right],$$ E7b

   where f₃ = f₁ + f₂ are the frequencies: φ₁, φ₂, φ₃ are the initial phases that are supposed to be independent random values varying within the limits of [0, 2π] and having uniform distribution law.

   Assume that the frequency f₃ in (7а) is an independent component, and an initial phase ϕ₃ corresponding to this frequency is an independent random value. Assume also that the frequency f₃ in (7b) is the result of phase coupling. It is evident that the power spectrum densities computed for the processes (7a) and (7b) are of the same shape because each considered process contains the contributions of the same frequencies f₁, f_2, and f₃, i.e., P_x1(p) = P_x2(p). At the same time, magnitude bispectrum of the process (7а) tends to be zero, i.e., ${\dot{B}}_{x1}\left(p,q\right)=0$ , and the magnitude bispectrum of the process (7b) is of non-zero value, i.e., ${\dot{B}}_{x2}\left(p,q\right)\ne 0$ .

6. Invariance property of bispectrum to a signal time delay or signal spatial shift: this property can be demonstrated by the following simple expression

   $$\begin{array}{ll}{\dot{B}}_{x\tau }\left(p,q\right)& ={\dot{X}}_{\tau }\left(p\right){\dot{X}}_{\tau }\left(q\right){\dot{X}}_{\tau }\left(-p-q\right)=\\ & =\dot{X}\left(p\right)\dot{X}\left(q\right)\dot{X}\left(-p-q\right)\text{exp}j2\pi \left\{-\mathit{\tau p}-\mathit{\tau q}-\tau \left[\left(-p-q\right)\right]\right\}={\dot{B}}_x\left(p,q\right),\end{array}$$ E8

   where ${\dot{X}}_{\tau }\left(p\right)=\dot{X}\left(p\right)\text{exp}\left(-j2\mathit{\pi \tau }\right)$ is the Fourier transform of a process shifted by a τ value in the temporal domain. It follows from (8) that both for the original process x^(m)(i) and for its replica x^(m)(i-τ) shifted by τ, the bispectra are congruent.

2. Extracting bispectral discriminative features from radar backscattering provoked by walking pedestrian

One of the most important problems arising in radar automatic target recognition (ATR) systems is the extraction of robust discriminative features contained in raw radar backscattering contaminated by clutter contribution. The Doppler velocity spectrum commonly serves for ground moving radar target detection, recognition, and classification. Often, the processed radar echo signals are of non-stationary behavior, whence the common notion of frequency becomes meaningless and a new parameter instantaneous frequency (IF) must be exploited in order to analyze how the Doppler spectral content varies with time. IF is defined as the frequency value existing within a short-time interval, and a given parameter allows one to estimate the distribution of spectral peaks varying with dwell time.

Note that the most frequently used tool for spectral analysis of non-stationary signals is an estimation of energy per unit time per unit frequency. Spectrogram is a typical representative of such energy-based estimators. Since the possible existing phase relationships between harmonics are lost in energy spectra, it is impossible to extract phase-coupled frequencies by means of time-varying energy spectrum. In other words, energy spectrum is unable to extract phase-coupled instantaneous frequencies (PCIFs) contained in time sequences processed. This peculiarity is the most serious shortcoming of energetic versions of evolutionary spectral estimators.

A walking person can be considered in the form of a complex physical phenomenon with simultaneous motion of different body parts: the torso, the legs, and the arms. Hence, for an extended human body, the backscattered field is a sum of multiple contributions.

Time-varying phase contributions are of paramount importance. Doppler frequencies provoked by swinging body parts may be phase coupled. In addition, radar returns obtained by applying vertical or horizontal polarizations differ one from another just by different spatial phase distribution of the reflection coefficients for corresponding vertical or horizontal polarizations.

The observed discrete-time received multicomponent and non-linearly FM radar returns corresponding to the vertical y^V(i) and horizontal y^H(i) polarizations can be written, respectively, as

$$y^V\left(i\right)=\sum _{m=1}^M{A}_m^V\left(i\right)\text{cos}\left[{\mathrm{\Phi }}_m^V\left(i\right)\right]=\sum _{m=1}^M{a}_m^V\left(i\right)\text{cos}\left[{\phi }_m^V\left(i\right)\right]F^2\left({\theta }_m\right)\text{cos}\left\{\frac{4\pi }{{\lambda }_0}\left[r_m\left(i\right)-r_0\left(i\right)\right]\right\},$$ E9a

$$y^H\left(i\right)=\sum _{m=1}^M{A}_m^H\left(i\right)\text{cos}\left[{\mathrm{\Phi }}_m^H\left(i\right)\right]=\sum _{m=1}^M{a}_m^H\left(i\right)\text{cos}\left[{\phi }_m^H\left(i\right)\right]F^2\left({\theta }_m\right)\text{cos}\left\{\frac{4\pi }{{\lambda }_0}\left[r_m\left(i\right)-r_0\left(i\right)\right]\right\},$$ E9b

where $a_m^V\left(i\right)$ and $a_m^H\left(i\right)$ are the time-varying magnitudes of the local reflection coefficients corresponding to an arbitrary mth object scattering center for vertical and horizontal polarizations, respectively; ${\phi }_m^V\left(i\right)$ and ${\phi }_m^H\left(i\right)$ are the local time-varying phases for vertical and horizontal polarizations, respectively; i = 1, 2, …, I is the sample index; F(θ) is the radar antenna pattern; θ_m is the angle location of the m-th object scattering center; r_m(i) and r₀(i) are the time-varying distances from antenna phase center to the arbitrary m-th moving object scattering center and object phase center, respectively; λ₀ is the radiated wavelength; A_m(i) and Φ_m(i) are the total radar echo signal magnitude and phase corresponding to the m-th object scattering center, respectively.

Two motion contributions are involved in the model (9a) and (9b): first, the translation motion of the human torso provoking approximately constant low Doppler frequency shift due to the relatively constant and low speed of the human body and, second, the swinging motion of the legs and arms provoking time-varying micro-Doppler frequency shifts caused by different non-uniform time-varying velocities of swinging body parts. These IFs may be related due to phase coupling.

Suggested approach takes into account the fact that the swinging torso, legs, and arms are not independent sources of Doppler frequency shifts but are related via the same “carrier” (the human torso) that can be considered as the “common basis” for the swinging legs and arms. Therefore, one can expect the presence of phase coupling contributions in radar returns (9a) and (9b).

Let us define the short-time bispectrum estimate of a transient signal y(i, n), whose time duration is significantly shorter than the total time observation interval of I samples contained in (9a) and (9b). Transient signal y(i, n) can be separated from the total observation (9a) and (9b) by a sliding window that moves step by step along the recorded process and takes n = 1, 2, …, N non-overlapping locations.

Time-varying bispectrum estimate can be computed according to (3) as follows

$$B\left(p,q;n\right)=Y\left(p,n\right)Y\left(q,n\right)Y^{\ast }\left(p+q,n\right)=\left|B\left(p,q;n\right)\right|\text{exp}\left[\mathit{j\beta }\left(p,q;n\right)\right],$$ E10

where Y(p, n) = |Y(p, n)| exp[jφ(p, n)] is the short-time Fourier transform of the transient signal y(i, n); |B(p, q; n)| and β(p, q; n) are the time-varying bimagnitude and biphase, respectively; p = 1, 2, …, P and q = 1, 2, …, P (P << I) are the frequency indices.

Coherent, homodyne, and surveillance radar operating in continuous-wave mode at the wavelengths of λ₀ = 2 cm, 3 cm, and 8 mm with the vertical polarization and polarimetric radar operating at the wavelength of λ₀ = 8.8 mm with the vertical and horizontal polarizations have been exploited for experimental study. The output of the 10 bits’ ADC allows collecting frequency radar returns containing Doppler frequency shifts. Experimental investigations have been carried out during the spring and summer periods. It is evidently seen from the histogram in Figure 1 the presence of a large number of the PCIFs contained in the radar return.

Figure 1.

The histogram of the PCIF distribution computed for a walking human.

First, we consider the case of vertical polarization. The radar returns were recorded during total data collection time of 63 s that corresponds to 200,000 digitized signal samples. The window width of 256 samples sliding along the signal provides frequency resolution of approximately 12.4 Hz. The sequence of time-varying bispectrum estimates contained 256 × 256 samples has been computed according to (10). Phase coupling contributions corresponding to the bimagnitude peaks were mapped from bifrequency plane onto the TF plane as two separate values on the frequency axis in the coordinates of p and q.

TF distribution of the PCIFs and bimagnitude estimates are demonstrated in Figure 2.

Figure 2.

TF distribution of the PCIFs for a walking human, vertical polarization, λ₀ = 2 cm. (a) Bimagnitude: vertical (b) and horizontal (c) polarization.

The analysis of the PCIF evolutionary behavior in Figure 2a indicates the contributions provoked by both the swinging pedestrian torso in the low-frequency range (30 Hz) and the legs and arms in the high-frequency range (from 120 through 240 Hz).

The following question arises: is there a difference between the radar Doppler TF signatures obtained with the vertical and horizontal polarizations and how much they differ one from another? In fact, according to the signal model (9a) and (9b), the last terms (electrical distances) in the quadruple products in the sums are the same for the vertical and horizontal polarizations. The distinction is due to the different behavior of the local time-varying phases ${\phi }_m^V\left(i\right)$ and ${\phi }_m^H\left(i\right)$ . It can be explained by the distinction in-phase spatial distributions in a human body for the vertical and horizontal polarizations. In order to confirm this hypothesis experimentally, we pay attention to the difference between the bimagnitude estimates computed for the vertical and horizontal polarizations for a pedestrian walking toward the radar.

The plots of the bimagnitudes are shown in Figure 2b and c for vertical and horizontal polarization, respectively. The difference between these bimagnitude estimates is clearly seen.

TF distributions computed by common energy-based strategy do not provide discriminating contributions provoked by the vertical and horizontal polarizations. Using proposed bispectrum-based approach permits recognizing the difference between PCIF TF distributions computed for the vertical and horizontal polarizations. This important peculiarity of the suggested approach can serve as a new perspective classification feature for solving classification problems in radar ATR systems.

Let us define the short-time cross-bispectrum estimate B_VH(p, q; n) in the form of

$$B_{\mathit{VH}}\left(p,q;n\right)=S_V\left(p,n\right)S_V\left(q,n\right)S_H^{\ast }\left(p+q,n\right)=\left|B_{\mathit{VH}}\left(p,q;n\right)\right|\text{exp}\left[j{\beta }_{\mathit{VH}}\left(p,q;n\right)\right],$$ E11

where S_V (p, n) and S_H (q, n) are discrete time-varying Fourier transforms of the transient signals y^V(i, n) and y^H(i, n), respectively; |B_VH (p, q; n)| and β_VH (p, q; n) are the time-varying cross-bimagnitude and cross-biphase, respectively.

Note that if there exists a phase coupling between two certain frequencies contained in the signals (9a) and (9b), a peak in the cross-bimagnitude distribution |B_VH (p,q;n)| will appear at the intersection of two frequency samples p and q in the bifrequency plane. Otherwise, for the Gaussian process provoked by a backscattering from vegetation, the cross-bimagnitude tends to be zero. Therefore, cross-bispectrum estimates (11) will contain the peaks caused by phase coupling. This important characteristic can serve for retrieving novel detection and recognition features for radar objects moving in a growing vegetation clutter.

The radar returns corresponding to the V and H polarizations both for a walking pedestrian and vegetation stirred by a light breeze were measured and recorded in the PC memory with the time interval of approximately 1 min. The sequence of N cross-bispectrum estimates (11) contained 256 × 256 samples each has been computed. The peaks of cross-bimagnitude |B_VH(p, q; n)| have been analyzed, and their values were projected into the TF domain. The plots of TF distributions computed for a human walking in the vegetation clutter and the vegetation clutter itself are shown in Figures 3 and 4, respectively.

Figure 3.

TF distribution for a human walking in vegetation clutter.

Figure 4.

TF distribution for vegetation clutter.

The number of pronounced cross-bimagnitude peaks is observed in Figure 3: the torso provokes the peaks in the low-frequency range, and the responses of swinging arms and legs are observed in a high-frequency range. The TF distribution computed for the vegetation clutter differs from that of a walking human substantially: the cross-bimagnitude peaks are observed only in the low-frequency range, and their values are considerably smaller in comparison to those of a walking human. It permits one to single out and identify a human walking in a vegetation clutter by the TF signature proposed.

Note that detection and classification of moving objects in through-foliage environment is one of the most important problems arising in surveillance radar systems serving for state border protection in the forest environment.

3. Classification of aerial radar targets by using bicoherence-based features extracted from micro-Doppler content

The problem of classification of various types of aerial targets by extraction of the information features contained in radar returns is of particular interest for ATR systems [2]. High-resolution radar range profile (HRRP) can be considered as a projection of a target spatial intensity distribution of the backscattered electromagnetic field onto the radar line of sight. Therefore, HRRP contains certain information about geometry of aircraft. Contributions caused by rotating turbine or propeller blades and contained in radar backscattering in the form of micro-Doppler data [3] also can contain supplementary classification information about aircraft’s particular characteristics. Both HRRPs and jet engine modulation (JEM) allow extracting classification features conforming to the aircraft geometry and particular engine characteristics.

As a result of aerial target travel over a large azimuth angle, its HRRPs measured and recorded during this travel considerably suffer by three main reasons: speckle, rotational range migration, and translational range migration. Due to the speckle effect, HRRP can fluctuate essentially even for slight aspect variations (about tenths of one angle degree), i.e., only a slight rotation of aerial target in elevation or aspect azimuth is enough to provoke considerable variations in the HRRPs. A typical example demonstrating considerable variability of HRRPs for the aircrafts AH-64, An-26, and B-52 under consideration is shown in Figure 5 for three fixed aspect angles of 185°, 190°, and 195°. These HRRPs have been computed by using radar backscattering models described in [2].

Figure 5.

Examples of dependence of HRRP on aspect angle.

It is well-known that performance achieved in ATR systems largely depends both on robustness of information features and dimensionality of classification feature vector. Consequently, the most serious limitation for using HRRPs for ATR is in their extreme variability. In addition, when the signal-to-noise ratio (SNR) is of rather low value, classification probability rate considerably degrades due to the contamination of radar signature by the noise. Therefore, search for both noise-robust and aspect-angle robust classification features is of paramount interest in ATR radar systems.

Phase information contained in radar returns can provide additional knowledge about the target. Unfortunately, phase information contained in radar returns in the form of phase coupling of micro-Doppler harmonics is irretrievably lost in common HRRP. In addition, micro-Doppler content is less sensitive to aspect angle variations compared to the HRRPs. This peculiarity can give certain advantages over HRRP-based approaches for ATR.

The approach which is most frequently cited for comparison of the performances of aircraft classification has been proposed by Zyweck and Bogner in [4]. We will refer to it as “Zyweck and Bogner technique” (ZBT). ZBT contains the following data processing steps.

• Computation of the following sequence of N HRRPs as

$$Y_n\left(m\right)={\left|K_n\left(m\right)\right|}^2,n=1,2,\dots ,N;m=1,2,\dots ,M,$$ E12

where K_n = S(:, n) ∈ Cⁿ is the matrix representation related to the sequence of N radar returns accumulated in the form of M in-phase (I) and M quadrature-phase (Q) digital samples, respectively.

• Each HRRP must be normalized to provide level invariance property. Aligning consecutive HRRPs (12) by using correlation procedure and following averaging of normalized HRRPs ${\widehat{\mathrm{Y}}}_{\mathrm{n}}\left(\mathrm{m}\right)$ as

$$R\left(m\right)=\sum _{n=1}^P{\widehat{Y}}_n\left(m\right),$$ E13

where P is the number of HRRPs selected for extraction of classification features.

• Direct Fourier transform is computed for averaged HRRP (13). In order to achieve translation invariance property, magnitude Fourier spectrum is used for extraction of classification features.

In practice, HRRPs consecutively accumulated according to the chain of received pulses have variable spatial shifts observed within fixed range strobe. These shifts are caused by translation of aircraft, and hence, migration of intensity contributions from one range cell to another is observed. ZBT exploits aligning procedure based on correlation analysis to provide translation invariance.

One more approach using classification feature extraction from HRRP is proposed in [5] by Kim et al. Below it will be referred as “Kim technique” (KT). The main idea behind KT is to use the first 20 central moments as the features. KT contains the following data processing steps:

1. - computation of range profiles according to (12) and their normalization;

2. - aligning consecutive HRRPs by using correlation procedure and averaging of aligned range profiles according to (13); and

3. - using the first 20 central moments related to averaged HRRP as classification features.

Micro-Doppler content can be extracted from the following accumulation of radar returns by the procedure represented in [3] as:

$$D\left(n\right)=\sum _{m=1}^M{K}_n\left(m\right).$$ E14

Joint TF distribution of radar echo-signal can be defined by the following short-time Fourier transform:

$$\mathit{TF}\left(f,t\right)=\sum _{l=1}^{L}D\left(l+\left(t-1\right)\left(L-Q\right)e^{\frac{j2\mathit{\pi ft}}{L}}w\left(l\right)\right),$$ E15

where L is the length of suitable smoothing window w(l); Q is the number of overlapping samples; f and t are the IF and time indices, respectively. The dimensionality of TF distribution (15) can be defined as C^LxG, where G = (P−Q)/ (L−Q).

Examples of spectrograms computed for three types of aerial targets are shown in Figure 6.

Figure 6.

Spectrograms of aerial targets computed for aspect angle equal to 185°: (a) helicopter AH-64 Apache; (b) military transport aircraft An-26; and (c) B-52 Stratofortress.

Each of these radar signatures contains micro-Doppler contributions inherent in particular aerial target. These radar signatures are of periodical behavior, and a period depends on the radial velocity of corresponding rotating parts. The contributions caused by front and back helicopter blades can be recognized in the spectrogram shown in Figure 6a. The front blade provokes larger radar cross-section (RCS) value and, in the position perpendicular to the radar line of sight, excites harmonics in the whole spectrum. The back helicopter blade provokes smaller RCS value and excites harmonics with smaller frequency bandwidth. The patterns referred to the aircrafts An-26 and B-52 and shown in Figure 6b and c, respectively, demonstrate higher periodicity frequency, and they look like more intricate radar signatures.

An approach proposed in [4] provides a translation invariance property for moving radar target classification dealing with cepstrum-based features. Rather low-dimensional and shift-invariant classification features can be extracted by using the cepstrum coefficients C(.). They are computed by using TF distribution (15) by the following transform

$$C\left(.\right)={\left|\mathit{\text{IDTF}}\left({\mathit{log}}_{10}\left(\frac{1}{T}\sum _{t=1}^T{\left|\mathit{TF}\left(:,t\right)\right|}^2\right)\right)\right|}^2,$$ E16

where IDFT denotes the indirect discrete Fourier transform.

Note that a dimensionality of the cepstrum-based classification features extracted by using transform (16) depends on required frequency resolution. It is relevant to note here that frequency resolution depends on the window length equal to L samples (see (15)).

Suggested feature extraction strategy is based on bicoherence estimation of radar data [6]. In order to compute bicoherence estimates from finite data records, first, the following evolutionary sequence of short-time bispectrum estimates B_t(f₁, f₂) corresponding to TF distributions (15) are computed as:

$$B_{\text{t}}(f_1,f_2)=TF(f_1,t)TF(f_2,t)T{F}^{*}(f_1+f_{2}t).$$ E17

Estimation of the bimagnitude of radar returns allows extracting information about presence of phase coupling in micro-Doppler content. It makes it possible to extract the features related to novel particular properties of target under classification.

Squared bicoherence version defined in [1] can serve as a quantitative measure of phase coupling. Squared bicoherence ${\widehat{b}}^2\left(f_1,f_2\right)$ can be interpreted as the proportion of signal energy corresponding to the frequency (f₁ + f₂) which is coupled with the frequencies f₁ and f₂ as follows:

$${\widehat{b}}^2\left(f_1,f_2\right)=\left({\left|\frac{1}{T}\sum _{t=1}^T{B}_t\left(f_1,f_2\right)\right|}^2\right)/\left[\left(\widehat{X}\left(f_1,f_2\right)\widehat{P}\left(f_1,f_2\right)\right)\right],$$ E18

where $\widehat{P}\left(f_1+f_2\right)=\frac{1}{T}{\sum }_{t=1}^T{\left|\mathit{TF}\left(f_1+f_2\right)\right|}^2$ is the power spectrum estimate averaged over T finite data samples; $\widehat{X}\left(f_1,f_2\right)=\frac{1}{T}{\sum }_{t=1}^T{\left|\mathit{TF}\left(f_1+f_2\right)\right|}^2$ .

The bicoherence (18) values vary within the limits of [0, 1]. If bicoherence tends to ${\widehat{b}}^2\left(f_1,f_2\right)\ne 0$ , it means that there exists a phase coupling at some pairs of frequencies, and if ${\widehat{b}}^2\left(f_1,f_2\right)=0,$ it means that there is no phase coupling in the signal under study.

Note that, unfortunately, the value of bicoherence estimate (18) depends on variations of target velocity. In order to provide target velocity invariance property, the following bicoherence-based feature (referred below as Bic) is introduced as

$$\mathit{Bic}=\left|F_{2D}\left({\widehat{b}}^2\left(f_1,f_2\right)\right)\right|,$$ E19

where F_2D denotes the 2-D Fourier transform.

Examples of bicoherence estimates (19) are illustrated in Figure 7. The bicoherence signature computed for helicopter AH-64 corresponds to the contribution of numerous phase-coupled harmonics, and ${\widehat{b}}^2\left(f_1,f_2\right)$ tends to the unity for a multiple number of points distributed in bifrequency plane. In contrast, bicoherence estimate computed for aircraft An-26 has the small number of phase-coupled harmonics. Moreover, only two points can be found with contribution larger than value of 0.1. The bicoherence estimate computed for aerial target B-52 has a number of frequency coupled responses with phase coupling coefficient larger than 0.6.

Figure 7.

Bicoherence estimates.

Two different types of classifiers are studied in simulations [7, 8, 9]. The first classifier is the support vector machine (SVM) belonging to the non-probabilistic linear classifiers, and the second type is the Naive Bayes Classifier (NBC) belonging to the probabilistic linear classifiers. Linear discriminant analysis (LDA) [5] is applied for reduction of feature vector dimensionality.

Modeled radar echo-signals backscattered by different flying aircrafts are computed in the sequences of I and Q samples by using the software described in [2]. Bandwidth of chirp radar signal equal to Δf = 150 МHz provides range resolution of ΔR = c/2Δf = 1 m. Pulse repetition frequency and central wavelength are equal to 2 KHz and 3 cm, respectively. The coherent chain of N = 2000 chirp pulses is received within observation time interval equal to 1 s.

Aerial targets of three different types are studied: helicopter AH-64 (the speed is of 160 km/h and height is of 200 m) and airplanes An-26 (the speed is of 400 km/h and height is of 2000 m) and B-52 (the speed is of 800 km/h, height is of 2000 m). For each target, data are available in the form of N = 2000 realizations related to the radar returns recorded within 1 s. Each radar return corresponding to reception of one pulse contains M = 1200 I and Q digital samples. Three fixed values of aspect angles equal to A = {185°, 190°, or 195°} are considered in computer simulations. The aspect angle equal to α = 180° corresponds to the situation, when the aircraft is flying toward the radar.

Classification probability rates are examined depending on the following three different scenarios. Under Scenario #1, a half of received signal samples related to each of three abovementioned aspect angles are used for training, and the remainder half is used as testing data. For extraction the TF-based classification features, the parameters used in TF distribution (19) are selected to be equal to L = 56 and Q = 46, and Hamming window is exploited in computer simulations.

According to Scenario #2, radar data conforming to one of three aspect angles equal to A = {185°, 190°, or 195°} are used for training database, and the records related to other two aspect angles are used for testing procedure. Step by step, each radar record from three available aspect angles is used one time for training and other two times for testing procedure. Under Scenario #3, data accumulated for two aspect angles are used for training procedure, and the remainder is used for testing. In the latter case, three sets of radar data corresponding to different aspect angles are examined, and three combinations of testing and training sets are possible for study. For each of them, classification probability rate is estimated, and final result is evaluated by averaging.

Table 1 shows the results obtained for three types of aerial targets under abovementioned three scenarios. It can be seen from Scenario #1 that both SVM and NBC classifiers provide approximately similar classification probability rates. The best performance is achieved by ZBT and Bic technique. HRRP-based schemes, i.e., ZBT and KT outperform the performances of TF-based techniques (Bic and Cepstrum) approximately by 1–4% for SVM classifier. This can be explained by a fact that the optimal values of parameters L, G, and Q given in (15) are not estimated. The best performance between the features related to the TF-based techniques is achieved by bicoherence-based (Bic) feature with NBC.

	Classifier	Feature extraction technique
		ZBT	KT	Cepstrum	Bic
Scenario #1	SVM	99.99	99.98	95.51	98.98
	NBC	99.97	93.06	98.67	100
Scenario #2	SVM	63.65	76.5	76.89	83.9
	NBC	71.8	79.32	77.38	–
Scenario #3	SVM	78.52	85.63	90.97	87.14
	NBC	75.92	79.46	88.15	82.78

Table 1.

Classification probability rate computed in percent.

Within Scenario #2, the best classification probability rate is achieved by using Bic technique. It allows to note that bicoherence-based classification feature demonstrates better robustness regarding the aspect of angle variations than other features. The Bic technique outperforms KT approximately by 4–8% and ZBT approximately by 12–20%.

Feature spaces of highest rank related to Scenario #2 are illustrated in Figure 8.

Figure 8.

Feature space after LDA, Scenario #2. Aspect angle of 185° is given for training data set, and aspect angles of 190° and 195° are selected for testing data set: (a) KT, (b) Cepstrum, and (c) Bic.

At the same time, the highest probability of correct classification is achieved by cepstrum-based features and SVM classifier for Scenario #3. Within Scenario #3, Bic technique takes the second place having the lags equal approximately to 4 or 6% depending on SVM or NBC classifier, respectively.

Summarizing the results presented in Table 1, the best robustness referred to aspect angle variations is demonstrated by TF-based schemes and extraction of cepstrum-based and bicoherence-based discriminative features. Taking into account the results obtained for Scenario #2, the information about phase coupling provides the best immunity regarding the aspect angle variations.

It can be seen from Figure 8c that variance of training samples is rather small compared to the distances observed between classes. Consequently, a problem occurs for computing the classification conditional probabilities for testing data. It follows from Figure 8a and b that testing sequences of samples are good when separated between each other. However, a variation of aspect angle leads to shifting of the features and misclassification in the final analysis. In contrast, the Bic-based features (see Figure 8c) belonging to testing sequence are overlapped for classes related to An-26 and B-52. As a result, classification errors can appear.

In order to study the influence of atmosphere turbulence, the same strategy is used as approach described for Scenario #1 with small difference. The training set is constructed in the same way by using clean, i.e., computed with less turbulence range profiles. At the same time, the testing set is constructed from the range profiles corrupted by a turbulence effect.

Classification probability rates computed under atmospheric turbulence influence are represented in Table 2. The following peculiarities should be emphasized by the comparison of the results demonstrated in Tables 1 and 2. Classification probability rates computed for ZBT, Bic-, and cepstrum-based approaches tend to decreasing approximately by 1–2% in a turbulence environment. Classifier using KT- and SVM-based decision making demonstrates the best robustness regarding influence of atmosphere turbulence.

Classifier	Feature extraction technique
	ZBT	KT	Cepstrum	Bic
SVM	98.22	99.99	97.09	97.82
Bayes	99.47	91.79	97.55	98.08

Table 2.

Classification probability rate computed for turbulence environment.

Classification performance evaluated for various SNRs is demonstrated in Figure 9. As can be seen from Figure 9a, the best performance is achieved by KT technique. Bicoherence-based classification feature becomes comparable to KT starting from SNR = 1 dB. The cepstrum-based approach outperforms Bic method only for SNR smaller than 2 dB. The worst performance is achieved by ZBT, and the loss is rather high comparing to other techniques. Concluding analysis of the results represented in Figure 9a, the best performance of 72% is achieved by KT even for SNR = −10 dB, whereas other techniques provide the probability of random guess of 33%.

Figure 9.

Classification probability rates depending on SNR value.

Bayes classifier demonstrates quite different results shown in Figure 9b. The Kim and Cepstrum techniques provide the same result as the results obtained with SVM classifier. The main difference is achieved by Bic method. For SNR ≥ −1 dB, the Bic method outperforms all other considered techniques. Errorless classification can be achieved starting from SNR = 2 dB for bicoherence-based features.

Thus, the best robustness to the aspect angel variability is provided by bicoherence- and cepstrum-based classifiers. These classifiers can be recommended for ATR radar systems.

4. Bispectrum-based radar range profile evaluation for the naval targets embedded in sea clutter environment

First, let us focus on the peculiarities characteristic for spatial-temporal model of surveillance radar signals received in sea clutter environment. Such signal model includes data contribution related to backscattering of naval surface target and interference contribution caused by a sea clutter. Assume that radar antenna irradiates the chain of M RF pulses of rectangular shapes and without any carrier modulation. Radar echo signal ${\dot{s}}_{\mathit{ik}}\left(t\right)$ can be written in the form of

$${\dot{s}}_{\mathit{ik}}\left(t\right)={\dot{S}}_{\mathit{ik}}\left(t\right)\text{exp}\left(j2\pi f_{0}t\right),$$ E20

where t ∈ [−T/2, T/2]; i, k = {H, V} are the indices corresponding to horizontal H and vertical V polarization, respectively; ${\dot{S}}_{\mathit{ik}}\left(t\right)$ is the complex-valued envelope; τ_p is the pulse length; T = T_rM is the total radar signal processing time, T_r is the pulse repetition period; m = 1, 2, 3,…, M defines the pulse number index; f₀ is the central frequency given in radar signal, f ∈ [f₀−ΔF/2, f₀−ΔF/2], and the frequency bandwidth ΔF is assumed to be of ΔF << f₀; $j=\sqrt{-1}$ .

During the time interval T, the naval target location on the sea surface can vary, and both its orientation and visibility can also be changed due to the wind’s influence. In order to obtain a smooth estimation of the radar range profile (RP) contaminated by contribution of the unpredictable variability of maritime environment, the common strategy presumes dividing time interval T into N separate short-time segments and performing the averaging for N observed realizations [2].

Assuming that contribution of additive noise is negligible, the nth (n = 1, 2, …, N) arbitrary echo signal realization (further, the nth scan) observed at the amplitude detector output can be written as

$$S_{\mathit{ik}}^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t\right)=\text{Re}\{\underset{\Omega }{\int }\dot{F}\left(\theta \right){\dot{\epsilon }}_{\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)S\left[l\left(\Delta \mathit{t}+T_r\right)-\tau \left(\theta \right)\right]\mathit{d\theta }$$ E21

where $\dot{F}\left(\theta \right)$ is the complex-valued antenna pattern; ${\dot{\epsilon }}_{\mathit{ik}}^{\left(n\right)}\left(\theta \right)$ is the nth relative (normalized for surface unit equal to 1 m²) complex backscattering coefficient; τ(θ)=2R(θ)/c denotes the time signal delay; R(θ) and c are the slant range and speed of light, respectively; l = 1, 2, …, L is the index of current temporal digital sample; L = M/N is the total sample number related to one scan; Δt is the pulse repetition period related to the range bin ΔR = cΔt. Eq. (21) is valid under spatial-temporal and band-limitedness condition, i.e., for 2ΔF/f₀ < λ₀/D, where D is the size of radar antenna aperture.

The complex-valued coefficients ${\dot{\epsilon }}_{\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)$ in (21) contain the sum of two contributions. First of them is random backscattering signal components ${\dot{\epsilon }}_{\mathit{Tik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)$ corresponding to the contributions of the echoes caused by small surfaces Δθ of the randomly moving target and second is the random component ${\dot{\epsilon }}_{\mathit{Sik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)$ caused by a sea clutter. These two contributions are statistically independent. Hence, ${\dot{\epsilon }}_{\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)$ can be expressed as

$${\dot{\epsilon }}_{\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)={\dot{\epsilon }}_{\mathit{Tik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)+{\dot{\epsilon }}_{\mathit{Sik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right),$$ E22

where the first term is given by signal polarization matrix contained the coefficients ${\dot{\epsilon }}_{T\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)=\left|{\dot{\epsilon }}_{T\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)\right|\text{exp}\left[j{\xi }_{T\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)\right]$ , and the second term is provoked by interference that can be represented by a polarization matrix containing the random coefficients ${\dot{\epsilon }}_{S\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)=\left|{\dot{\epsilon }}_{S\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)\right|\text{exp}\left[j{\xi }_{S\mathit{ik}}^{\left(n\right)}\left(\theta ,\mathrm{\Delta }\mathit{t}\right)\right]$ .

Commonly, the RP estimates are obtained by using the averaging strategy performed over N observed realizations. The mixture of signal and interference can be written as

$${\stackrel{⌢}{S}}_{\mathit{ik}}\left(\mathit{l}\mathrm{\Delta }t\right)=<{\tilde{S}}_{\mathit{ik}}^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t\right)>=<{\tilde{S}}_T^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t-{\tau }_T^{\left(n\right)}\right)>+<{\tilde{S}}_S^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t-{\tau }_S^{\left(n\right)}\right)>,$$ E23

where ${\tilde{S}}_T^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t-{\tau }_T^{\left(n\right)}\right)$ and ${\tilde{S}}_S^{\left(n\right)}\left(\mathit{l}\mathrm{\Delta }t-{\tau }_T^{\left(n\right)}\right)$ are the nth echo envelopes smoothed by antenna pattern and corresponding to both target and sea clutter backscattering, respectively; ${\tau }_T^{\left(n\right)}$ is the time lag integrated for all object backscattering centers during the nth scan; ${\tau }_S^{\left(n\right)}$ is the time lag integrated for all sea backscattering elements during the nth scan; <… > denotes the averaging for N observed scans. Both ${\tau }_T^{\left(n\right)}$ and ${\tau }_S^{\left(n\right)}$ are the random values that vary from one realization to another.

Due to random motions of target on the sea waves, the backscattering signal is a fluctuation process. Different wave propagation paths lead to different time lags ${\tau }_T^{\left(n\right)}$ and ${\tau }_S^{\left(n\right)}$ contained in (23), and these lags vary randomly. The target and sea radar responses are independent processes for which the corresponding PDF and correlation intervals differ from each other. The correlation intervals for processes that correspond to backscattering from sea surface are comparable to one scan processing time. Therefore, the contribution caused by sea clutter is an unpredictable fluctuating process that rapidly varies from one scan to another. Therefore, the time-varying sea echo radar responses that usually overlap with the signal target response can destroy the estimate (23).

Let us consider bispectrum-based data processing approach aimed for smoothing maritime interference contribution. Bispectrum estimate $\widehat{{\displaystyle \dot{B}}}\left(p,q\right)$ of the process (20) can be written as

$$\widehat{{\displaystyle \dot{B}}}\left(p,q\right)=\left|\widehat{B}\left(p,q\right)\right|\text{exp}\left[j{\displaystyle \widehat{\gamma }}\left(p,q\right)\right]=<X^{\left(n\right)}\left(p\right)X^{\left(n\right)}\left(q\right)X^{\left(n\right)\ast }\left(p+q\right)>,$$ E24

where $\left|\widehat{B}\left(p,q\right)\right|$ and $\widehat{\gamma }\left(p,q\right)$ are the bimagnitude and biphase estimates, respectively; p = 1, 2, …, L and q = 1, 2, …, L are the frequency indices; ${\dot{X}}^{\left(n\right)}\left(\dots \right)$ is the direct Fourier transform of (20) equal to ${\dot{X}}^{\left(n\right)}\left(p\right)={\dot{S}}_T^{\left(n\right)}\left(p\right)\text{exp}\left(j2{\mathit{\pi \tau }}_T^{\left(n\right)}p\right)+{\dot{S}}_S^{\left(n\right)}\left(p\right)\text{exp}\left(j2{\mathit{\pi \tau }}_S^{\left(n\right)}p\right)$ ; ${\dot{S}}_T^{\left(n\right)}\left(p\right)$ and ${\dot{S}}_S^{\left(n\right)}\left(p\right)$ are the Fourier transforms of the object ${\tilde{S}}_T^{\left(n\right)}\left(\dots \right)$ and sea ${\tilde{S}}_S^{\left(n\right)}\left(\dots \right)$ contributions given in (23), respectively.

Bispectrum estimate $\widehat{{\displaystyle \dot{B}}}\left(p,q\right)$ (24) contains eight terms. The first term in the form of ${\widehat{{\displaystyle \dot{B}}}}_T\left(p,q\right)=<{\dot{S}}_T^{\left(n\right)}\left(p\right){\dot{S}}_T^{\left(n\right)}\left(q\right){\dot{S}}_T^{\left(n\right)\ast }\left(p+q\right)>$ corresponds to the contribution of original target bispectral estimate. The other seven terms represent the interference contribution. Theoretically, considerably good accuracy of signal contribution ${\widehat{{\displaystyle \dot{B}}}}_T\left(p,q\right)$ can be obtained under traditional assumptions that interference is of zero mean, and its PDF is close to a symmetric law. However, in number of real-life situations, interference does not obey the latter symmetric law. Therefore, the real-life case of sea clutter with non-zero mean value and with long-tail PDFs needs to be studied, and this problem is one of the subjects of our consideration and discussion in this chapter.

Final step of the bispectrum-based forming of RP estimate can be represented as the following inverse Fourier transform (IFT)

$${\widehat{S}}_{\mathit{\text{range}}}\left(l\right)=\left|\mathit{IFT}\left\{\left|{\widehat{\dot{S}}}_{\mathit{\text{bisp}}}\left(r\right)\right|e^{j{{\displaystyle \widehat{\varphi }}}_{\mathit{\text{bisp}}}\left(r\right)}\right\}\right|r=1,2,\dots ,L,$$ E25

where $\left|{\widehat{\dot{S}}}_{\mathit{\text{bisp}}}\left(r\right)\right|$ and ${\widehat{\varphi }}_{\mathit{\text{bisp}}}\left(r\right)$ are the magnitude and phase target RP Fourier spectrum estimates, respectively. They can be recovered from (24) by recursive algorithms described in details in [10].

The main difference between traditional signal processing (23) and suggested bispectrum-based approaches can be explained by the difference in averaging procedures necessary to smooth the sea clutter contribution. Traditional signal processing (23) using direct ensemble averaging is performed in the signal space, i.e., in time domain. Due to the random signal lags ${\tau }_T^{\left(n\right)}$ contained in (23), averaged received response estimate (23) becomes considerable spread shape. It provokes decreasing the performance in ATR system by using common RP evaluation. At the same time, suggested bispectrum-based technique exploits coherent averaging of non-translated bispectrum estimates in bifrequency domain. As a result, due to the invariance property of bispectrum to a signal lags ${\tau }_T^{\left(n\right)}$ , suggested strategy provides coherent ensemble averaging of bispectrum estimates. Therefore, the shape of RP transformed considerably less as compared to the common signal averaging procedure (23). Hence, better performance of bispectrum-based ATR system can be achieved.

Experimental study was performed during summer period using X-band polarimetric surveillance radar operating at the central frequency equal to f₀ = 9.370 MHz. The fixed transmitting/receiving radar antenna was located on the sea shore at the height of y₀ = 8 m over sea level.

After passing through analog IF amplifier and amplitude detector, the received radar echo signals were digitized and accumulated in the processor memory. The sampled experimental data were accumulated and recorded in the form of the sets of scans (realizations) contained L = 32 samples in each scan for each HH, HV, VH, and VV polarization, respectively. The scan duration for each polarization was set to 320 ms.

The anchored metallic buoy was served as the naval surface target. Its size was considerably smaller than the radar range bin. The radar echoes were recorded under grazing angles of about 0.4°.

We compare the bispectrum-based RP recovery with RP estimations evaluated by the common direct averaging procedure. The number of N scans accumulated for each separate polarization mode was equal to N = 256.

Consecutive two arbitrary scans experimentally measured for HH polarization mode are demonstrated in Figure 10.

Figure 10.

Consecutive scans.

The plots of normalized range profile (NRP) are represented in the form of function of range sample index l. A priory known slant range from antenna to buoy was equal to R = 1500 m. Therefore, the buoy location corresponds approximately to the range sample index l = 26. Note that random buoy translation on sea waves can be expected to become apparent during considerably long observation time interval. Hence, the target ${\dot{\epsilon }}_{\mathit{Tik}}^{\left(n\right)}\left(\theta ,\mathit{l}\mathrm{\Delta }t\right)$ and sea clutter ${\dot{\epsilon }}_{\mathit{Sik}}^{\left(n\right)}\left(\theta ,\mathit{l\Delta t}\right)$ contributions given in (22) are the values that are randomly varying during time interval related to each separate observed scan.

As can be seen from Figure 10, both radar backscattering response corresponding to buoy RP and to sea clutter have random nature, and their appearance and shapes considerably vary from scan to scan. Particularly, the buoy is not visible for some scans at all, e.g., for scan in Figure 10b. Though the target RP should appear itself as single peak because its size is much smaller than the radar range bin, the target response shape is considerably distorted. Due to sea clutter contribution, several peaks are observed in the neighborhood of the range index number l = 26 related to the buoy location. Such distortions in the NRP prevent target classification in ATR system.

In order to improve estimator accuracy, the average procedure performed over N scans is commonly used according to (23). Example of averaged NRP is demonstrated in Figure 11a. The bispectrum-based NRP signature is shown in Figure 11b. Note that due to bispectrum shift invariance property, the latter NRP is centered with respect to the center of gravity.

Figure 11.

Averaged RP, VV polarization, wind speed 7–10 m/s; sea state 2–2.5 (a). Bispectral NRP, VV polarization, wind speed 7–10 m/s; sea state 2–2.5 (b).

As clearly seen from the Figure 11a, the NRP ensemble averaged according to (23) is considerably distorted. The shape of averaged NRP is spread considerably. This value is essentially larger than the radar range bin equal to 75 m that can be expected as the target response width. This spreading results in decreasing of radar range resolution caused by sea clutter influence. It is clearly seen that range resolution obtained with bispectral NRP (see Figure 11b) is considerably better as compared with averaged RP (see Figure 11a).

In Tables 3 and 4, the summary of the key parameters of the NRP estimation, namely the range resolution data under sea clutter level, is represented for two considered techniques.

Absence of wind			Wind speed of 7–10 m/s
HH	HV	VV	HH	HV	VV
NRP width at the half-amplitude level, m
450	300	300	900	300	300
Sea clutter level, dB
−13	−6	−6	−13	−8	−13

Table 3.

Experimental results obtained for common averaging procedure according (23).

Absence of wind			Wind speed of 7–10 m/s
HH	HV	VV	HH	HV	VV
NRP width at the half-amplitude level, m
150	150	150	450	150	150
Sea clutter level, dB
−24	−20	−26	−17	−28	−26

Table 4.

Experimental results obtained for suggested bispectrum-based data processing.

As can be seen from Table 3, the range resolution is worsen from 2 to 12 times as compared to the theoretical radar range bin equal to 75 m.

The results represented in Table 4 demonstrate the stability of the range resolution that is worsen only by 2 times in comparison to the above mentioned theoretical radar range bin except the case of wind speed of 7–10 m/s, as well as for HH polarization.

The values of the sea clutter levels given in Tables 3 and 4 depend both on polarization mode and sea state. As clearly seen from comparison, the data contained in the Tables 3 and 4, bispectrum-based technique provides suppression of sea clutter level from 11 dB (HH polarization mode) to 20 dB (VV polarization mode) in the case of absence of wind and from 4 dB (HH polarization mode) to 20 dB (HV polarization mode) in the case of wind speed equal to 7–10 m/s. Hence, the bispectrum-based technique provides considerable improvement of range resolution in sea clutter environment.

Consequently, conventional RP estimation technique based on the averaging of the received signal envelopes has low-range resolution and low robustness to sea clutter.

5. Classification of the atmospheric formations by using bicoherence-based features

Recently, weather radars play an important role in the world for atmosphere surveillance and weather forecasts. Radar techniques dedicated to study the meteorological formations are based on extraction of the data from backscattering of electromagnetic waves by the particles concentrated in the clouds and precipitations. Spatial distribution of the numerous backscattering centers contained in a cloud creates the contribution into total radar echo signal. Parameters of this backscattering signal contain meteorological information about atmosphere formations [11, 12].

Reliable meteorological data serve as one of the main and important component necessary for securing flight safety in modern aviation. According to statistical data, dangerous phenomena like thunderstorm frontal passages, squalls, and wind fluctuations cause approximately 70% of the aviation crashes. Turbulent movement of atmospheric flows as opposed to laminar air movement is specified with random variability of wind-rate field, presence of fluctuation of air heterogeneity, or so-called turbulent atmospheric vortexes that cause mixing of the airflows. Atmospheric turbulence causes abrupt movement of aircraft called as bumpy flight. As the bumpy flight piloting becomes a difficult and dangerous procedure, additional mechanical loads related to the elements of aircraft construction appear and even destruction of aircraft construction may happen. Consequently, reliable recognition of dangerous turbulent formations and areas is an important problem in modern aviation radar meteorology and aviation safety of the flights.

Power spectral density of radar backscattering signals is one of commonly used and widespread information features extracted from “angel echo.” Usually, the severity of turbulence is accessed from the measurements of the radial velocity variance observed in the air flows. The power spectral density describes the dissipation of turbulence energy caused by decay of large atmospheric vortexes into smaller fractions. In order to recognize the type of cloud formation, it is commonly accepted that it is enough to estimate the power spectral width of radar backscattering signal. In order to make a decision about atmospheric turbulence state and, hence, about dangerous areas located at the course of aircraft, the value referred to the power spectrum spreading must be estimated. After that, the pilot must make a decision in order to avoid a weather hazard area. However, the spreading of power spectrum can be related to sufficiently rough classification feature serving for estimation of the scale of turbulence contained in the atmosphere formations. Solving the problems of discrimination of turbulent and laminar atmosphere flows by using spectral density width can lead to ambiguous interpretation of measured data and, consequently, to provoke the errors related to the correction of the aircraft course in order to round the dangerous atmospheric turbulent formations and areas.

In order to improve reliability of recognition and discrimination of laminar and turbulent meteorological formations, it is necessary to seek novel information features contained in radar backscattering signals. Spatial-temporal, frequency-phase relationships, as well as phase-coupling phenomena may serve as novel discrimination and classification features. These frequency-phase relationships can be caused by physical nature of radar signal scattering with laminar aerial flow contained within a pulse volume.

In order to extract phase coupling contained and prevail in laminar atmospheric formation, we suggest exploiting bispectrum-based strategy to the weather radar signal processing.

The approach suggested in this section is based on the hypothesis that radar signal backscattered by atmospheric formation can contain phase-coupled spectral contributions carrying useful data for extraction of the classification features. This phase coupling contains novel classification features used for discrimination and classification of turbulent and laminar atmospheric flows and, hence, information about dangerous areas located in the atmosphere. As novel information feature and measure, we suggest exploiting bicoherence in the form defined in [13].

Bicoherence b²(p, q) [13] can be computed in the form of statistical averaging of the bispectral estimate and further normalization procedure as

$$b^2\left(p,q\right)={\left|<B\left(p,q\right)>\right|}^2/\left\{<{\left|X\left(p+q\right)\right|}^2<{\left|X\left(p\right)X\left(q\right)\right|}^2>\right\}$$ E26

Bicoherence b²(p, q) given in the form of (26) can be interpreted as squared version related to the normalized bispectrum. The dimensionless value (26) varies within the limits from zero to unity, i.e. b²(p,q) ∈ [0, 1]. So, bicoherence is the value serving as a quantitative measure of phase relationships between spectral components contained in a signal under study.

If the electromagnetic waves are backscattered by volume-distributed turbulent meteorological formation, all the phases in backscattering signals are of random and independent values. In such case, bicoherence estimate (26) tends to be of zero value. Phase-coupled spectral components are contained in the radar backscattering caused by volume-distributed laminar formation. In latter case, bicoherence estimate (26) tends to be a non-zero value.

Experimental study was performed by upgraded weather radar of the type of MRL-1 equipped by additional facilities necessary for control of orientation of maximum of antenna main lobe pattern, calibration of radar energy potential, expansion of dynamic range in receiving device, and digital signal processing. Weather radar operates at the wavelengths of 8 mm and 3.2 cm.

The power spectral density and bicoherence of the echo signal measured for the case of antenna elevation angle equal to 10° and the downrange of 3750 m are plotted in the Figure 12. These demonstrative examples correspond to the stratocumulus “non-translucent” clouds #122.

Figure 12.

Power spectrum (a) and bicoherence estimate (b) for “non-translucent” clouds #122.

Due to stochastic behavior of the echo-signals, their power spectral density function is of pronounced irregular shape. Because of this, in order to smooth power spectral density estimate, statistical averaging has been performed over each four realizations. After digitalization procedure and averaging performed over four realizations, 512 digital samples corresponding to each range bin were accumulated in PC memory. These data are used to compute both power spectral density and bicoherence estimate.

Computation of smoothed bicoherence estimate has been accomplished by averaging over eight segments contained 64 received signal samples with 50% overlap in each segment. Additional smoothing of the bicoherence estimate has been carried out under fixed antenna elevation angle equal to 10°. Ensemble averaging has been performed over 100 observed realizations and with time interval equal to 0.85 s. Therefore, total time interval exploited for computation of bicoherence estimate was equal to 100 × 0.85 s = 85 s.

Decision about belonging of meteorological formation to laminar or turbulent kind is made by comparison the bicoherence value with threshold given a priory. In the case, when the bicoherence samples do not take the values belonging to the given threshold, decision is made in favor of turbulent formation. In opposite case, one made decision in favor of laminar formation.

Two different threshold values equal to 0.7 and 0.5 were examined (see Table 5).

Type of cloudiness	Slant range (m)	Index of excess
		Threshold 0.7 (%)	Threshold 0.5 (%)
Clouds #122	3750	0.02	0.76
White cumulus clouds	3150	0.77	11.2
Flocculent clouds T-36	5100	16.19	95.0
	3300	18.4	96.7
	1500	0.6	3.7
“White flaky”	2400	0.1	4.07

Table 5.

The indexes of exceeding of the threshold computed in percent.

Data represented in Table 5 indicate the following results. For the stratocumulus “non-translucent” clouds #122, the percentage of the bicoherence values exceeding both threshold values of 0.7 and 0.5 are equal to 0.02 and 0.76%, respectively. For the flocculent altocumulus clouds T-36, the percentage of the bicoherence values acceding the thresholds of 0.7 and 0.5 are equal to 16.19 and 95%, respectively. It is also seen from the data contained in Table 5 that for the clouds #122, the turbulent contributions are predominated at the slant range of 3750 m. At the same time, the laminar content is predominated for the T-36 clouds at the slant range equal to 5100 m.

Thus, using common classification features evaluated in the form of width of power spectral density leads to certain ambiguity. At the same time, suggested bicoherence-based approach proposed in this chapter eliminates this ambiguity.

6. Conclusions

Bispectrum density estimator in opposite to the common energy spectrum estimator allows not only describing statistical characteristics of a process more correctly and more profoundly but also to detect and extract a new class of dependences contained in the data under study. These dependences can exist in the form of spectral component correlation relationships and phase coupling between the spectral component pairs. Therefore, the main difference between bispectral and energy spectral strategies is in the preservation of phase information and possibility of extracting this phase contributions. Bispectrum-based signal processing allows extracting novel information features providing signal detection and discrimination, as well as object recognition and classification. The benefits of the suggested bispectrum-based data processing techniques were demonstrated by experimental study of radar target detection, classification and identification for naval, aerial, and ground moving objects. Experimental results represented in this chapter demonstrate sea clutter suppression and improving of naval object range resolution provided by polarimetric X-band radar. Essential reduction of speckle distortions in aerial high resolution range profiles obtained by using suggested bicoherence-based classification features was shown. Experimental results of time-frequency analysis of backscattered signals recorded by ground surveillance Doppler radar were represented and discussed. Suggested bispectrum-based approach can serve for improving the detection and recognition performances in radar ATR systems operating in vegetation clutter. Bicoherence estimates were proposed for detection and discrimination of the atmospheric turbulence formations performed by weather surveillance radar. Results of experimental study indicate discrimination of the turbulent and laminar air flows by exploiting bicoherence classification features.