I. INTRODUCTION
In the field of non-cooperative target recognition with radar, several different techniques have been introduced for the classification of air targets. All of these techniques extract information from the scattered signal, like micro Doppler modulations from moving parts, especially jet engine modulation (JEM) produced by the fan or compressor of the engine [Reference Chadwick and Williams1]. Other techniques, like high-resolution range profiles [Reference Bieker2, Reference Haumtratz, Worms and Schiller3] or inverse synthetic aperture radar (ISAR) [Reference Chadwick and Williams1, Reference Rosenbach and Schiller4], exploit the imaging capabilities of radar. For these techniques a waveform with a high bandwidth, like a linearly frequency modulated (LFM) chirp, is necessary [Reference Chen and Martorella5]. Most of the high-resolution-based classifiers are trained on a limited amount of similar targets and are used to distinguish between them. This procedure is based on the assumption that the class such as airliner, business jet, fighter or transport aircraft can be discriminated by a previously applied classifier.
This paper is an extension of [Reference Wagner, Dommermuth and Ender6] and a technique, which can be used as part of a pre-classification stage to eliminate a huge number of potential targets, is presented. With the proposed decomposition of ISAR images, the number and position of jet engines of targets visible in ISAR images can be determined. With these features a separation of targets with wing- and fuselage-mounted engines is possible, further is the number of engines significant for a list of possible targets used for the final classification [Reference Rihaczek and Hershkowitz7]. The efforts presented in the past to locate jet engines in ISAR images are based on the Doppler spread and the resulting JEM lines. In [Reference Tait, Griffiths and Blacknell8], range-Doppler signatures with short dwell times were used to detect jet engines along the target. In [Reference Ricardi, Dell'Acqua and Aprile9], an algorithm to detect the range position of jet engines in ISAR images was introduced that utilizes the occurring JEM lines. A drawback of both given examples is that only the range position of the jet engines can be determined with these methods. With the method proposed here also the cross range position and the number of engines can be extracted.
For the feature extraction, a sparse decomposition of the ISAR image is performed. To achieve a separation between engine scatterers and normal scattering centers, two different signal models are necessary to form the dictionary for the sparse decomposition. As signal model for the normal scattering centers, the common isotropic point scattering model is used. For the engine scattering, a waveguide scattering model, which is very common to model jet engines [Reference Anastassiu10] is utilized. This model describes a dispersive behavior of the cavity of the engine, such that the LFM chirp is extended in time and produces artifacts in the image. For these artifacts, a mathematical description that is used as model for the decomposition algorithm was introduced by Borden [Reference Borden11]. Both models are reviewed in Section II and the sparse decomposition framework is presented in Section III. An example of the decomposition applied to a real ISAR image, measured with the tracking and imaging radar (TIRA) of Fraunhofer FHR, is shown in Section IV. Finally, Section V concludes the paper and points out some open problems for future research.
II. ISOTROPIC AND WAVEGUIDE SCATTERING SIGNAL MODELS
For the proposed framework two distinguishable signal models of jet engines and point targets are necessary. The general idea of the engine model is to model the entrance structure as a circular terminated waveguide. Within this model the scattered field
$\vec E^S$
of the jet engines can be considered as a sum of three independent parts, which are depicted in Fig. 1. These scattering mechanisms are a strong reflection
$\vec E_{rim}^S $
from the edge at the opening of the engine, weak reflections
$\vec E_{ext}^S $
from the exterior of the engine and the reflection
$\vec E_{wg}^S $
from the fan or the first stage of the compressor inside the engine [Reference Moghaddar and Walton12]
$$\vec E^S = \vec E_{rim}^S + \vec E_{ext}^S + \vec E_{wg}^S. $$
Fig. 1. Different scattering mechanisms of a jet engine.
The first two elements of the sum are non-dispersive and can be modeled as common isotropic frequency-independent scatterers. With y as the range and x as the cross-range direction this model is written as
$$s_p(k_x,\,k_y,\,x,\,y) = e^{ - j2(k_xx + k_yy)},$$
which is used to describe a scattering center in the two-dimensional (2D) k-space domain [Reference Chen and Martorella5]. k x and k y are the wavenumbers in there specific directions, calculated with the small aperture approximation by k x = ksin(Θ) ≈ kΘ and k y = kcos(Θ) ≈ k, where Θ is the azimuth angle of the target and k is the wavenumber of the transmitted instantaneous frequency. With a band-limited 2D inverse Fourier transform we get the point spread function (PSF) in the image domain, a 2D sinc-function
$$\eqalign{PSF_p({x}^{\prime},{y}^{\prime}) = & \int\!\!\!\!\!\int_{\opf K} e^{ - j2(k_xx + k_yy)}e^{j2(k_x{x}^{\prime} + k_y{y}^{\prime})}dk_xdk_y \cr \approx & A \cdot {\rm sinc}(\bar k\Delta \Theta ({x}^{\prime} - x)) \, {\rm sinc} (\Delta k({y}^{\prime} - y)),} $$
which will be used in the dictionary for the standard scatterer at position (x′,y′). The inverse Fourier transform is band-limited because the covered area
$\opf{K}$
of k
x
and k
y
is limited. ΔΘ is the totally covered aperture during the imaging process, Δk is the used bandwidth in k-space k
max − k
min and
$\bar k = (k_{{\rm min}} + k_{{\rm max}})/2$
. The complex scaling factor A can be ignored in the proposed framework, since this is working on the magnitude image and uses normalized PSFs. An example of the PSF, weighted with a Hamming window, is shown in Fig. 3(a).
Inside a waveguide, the propagation speed of the electromagnetic wave depends on the frequency, which is called waveguide dispersion. This dispersion will spread the transmitted LFM waveform in time and thus the pulse compression, which is fitted to a certain pulse length, will lead to delayed returns in the received signal. These duct delayed returns produce artifacts in the range direction of the image behind the physical position of the waveguide structure. An example of these delayed returns is depicted in Fig. 2(a), which is a nose on ISAR image of a small business jet. The duct delayed returns can be seen behind the engines, of which the position is visible in the target shape in Fig. 2(b).
Fig. 2. Image and shape of the target. (a) Testimage, (b) shape of the target.
To model this waveguide effect, the standard model s p in (2) must be corrected with a dispersion term including a wavenumber k y in range, which is nonlinearly dependent on the frequency ω
$$s_{wg}(k_y,L,\omega ) = e^{ - j2Lk_y(\omega )}. $$
This term represents the propagation of the incident wave to the terminated end of the waveguide over its length L and back. The general form of the wavenumber in dispersive media is given by [Reference Orfanidis13]
$$k(\omega ) = \displaystyle{{ - j} \over {c_0}}\sqrt {j\omega + a + b} \sqrt {j\omega + a - b}. $$
The appropriate choice of a and b in a hollow metallic waveguide are a = 0 and b = jω c , which is the cut-off frequency of the propagating mode. With this substitution, (5) results in
$$k_{wg}(\omega ) = \displaystyle{1 \over {c_0}}\sqrt {\omega ^2 - \omega _c^2}. $$
Using (6) as k y (ω) in (4), the corrective term becomes
$$s_{wg}(k_y,L,d) = e^{ - j2L\sqrt {k_y^2 - \kappa _m^2}}, $$
where κ m is the cut-off wavenumber of the mode m. It is known from electromagnetic theory that inside waveguides electromagnetic waves propagate in discrete modes κ m , which depend on the diameter d of the opening [Reference Orfanidis13]. Thus, (7) must be replaced by a sum over the propagating modes if an inlet is located at the corresponding position [Reference Borden14],
$$s_{wg}(k_y,\Theta, L,d) = \displaystyle{1 \over k_y}\sum\limits_{m = 1}^M S_m(\Theta, k_y)e^{ - j2L\sqrt {k_y^2 - {\kappa _m^2}}}. $$
The term S m (Θ, k) in (8) represents the energy spread of the waveform. Since the opening of the engine is large compared with the wavelength of the radar, many modes will propagate inside the waveguide and the energy of the waveform is spread over them. The number of modes M is limited by the largest transmitted frequency k max,
$$M = \max \left\{ {m \vert \kappa _m \le k_{{\rm max}}} \right\}.$$
To calculate the PSF of a waveguide scatterer, the inverse Fourier transformation of the corrected range term with the standard model in cross-range must be determined, which leads to [Reference Borden11]
$$\matrix{ \matrix{PSF(x,y) \hfill \cr \quad \approx \underbrace{{{\rm sinc}(\bar k\Delta \Theta x)}}_{{PSF_x}} \;\underbrace{{\sum\limits_{m = 1}^M S_m\int_{2k_{{\rm min}}}^{2k_{{\rm max}}} \displaystyle{1 \over {k_y}}e^{ - j2L\sqrt {k_y^2 - \kappa _m^2}} e^{jk_yy}dk_y}}_{{PSF_y}}. \hfill} \cr} $$
Equation (10) is a multiplication of two separated imaging kernels in the x- and in the y-direction. The kernel in the x-direction only depends on the carrier frequency
$\bar k$
and the complete rotation angle ΔΘ of the target and is the same as for the standard signal model. The kernel in range direction is a weighted sum of integrals, that still has to be examined and for a more detailed version of the following deviation of the PSF, the reader is referred to the original paper of Borden [Reference Borden11] or his textbook [Reference Borden14]. The first step is the substitution
$$\matrix{ \matrix{\displaystyle e^{ - j2L\sqrt {k_y^2 - {\kappa _m^2}}} = 2j\sqrt {k_y^2 - \kappa _m^2} \int_L^\infty J_0\left( {2\kappa _m\sqrt {{{y}^{\prime}}^2 - L^2}} \right)e^{ - j2k_y{y}^{\prime}}d{y}^{\prime}, \hfill} \hfill \cr} $$
where J 0( · ) is the Bessel function of the first kind. With this substitution and a rearrangement of k y , the integral in (10) can be written as
$$\eqalign{& 2j\int_L^\infty {J_0} \left( {2\kappa _m\sqrt {y^{\prime 2} - L^2} } \right) \cr & \times \int_{2\max (k\min ,\kappa _m)}^{2k_{\max }} {\sqrt {1 - {\left( {\displaystyle{{\kappa _m} \over {k_y}}} \right)}^2} } e^{j2k_y(y - y^{\prime})}dk_ydy^{\prime},} $$
where the lower bound of the initial integral in (10) is changed since components below the cut-off frequency of the specific mode will be exponentially damped. The expression in (12) can be interpreted as a convolution of the Bessel function with the second integral, which can be approximated as
$$\eqalign{& \int_{2\max (k_{\min}, \kappa _m)}^{2k_{\max}} \sqrt {1 - {\left( {\displaystyle{{\kappa _m} \over {k_y}}} \right)}^2} e^{j2k_y(y - {y}^{\prime})}dk_y \cr & \quad \approx \alpha \Delta k_{(1,c)}{\rm sinc}\left( {\Delta k_{(1,c)}(y - {y}^{\prime})} \right)e^{j2{\bar k}_{(1,c)}(y - {y}^{\prime})}.} $$
In (13), the width of the sinc-function depends on the bandwidth of this specific mode, which is defined as Δk
(1,c) = k
max − max (k
min, κ
m
).
$\bar k_{(1,c)}$
is the center frequency of the mode defined as (max (k
min, κ
m
) + k
max)/2. With these definitions the final version of the imaging kernel in range direction is
$$\eqalign{PSF_y = & \sum\limits_{m = 1}^M 2jS_m\alpha \Delta k_{(1,c)}\int_L^\infty J_0\left( {2\kappa _m\sqrt {{{y}^{\prime}}^2 - L^2}} \right) \cdots \cr & \times {\rm sinc} (\Delta k_{(1,c)}(y - {y}^{\prime}))e^{j2{\bar k}_{(1,c)}(y - {y}^{\prime})}d{y}^{\prime}.} $$
From (14) it is visible that the down-range effect in the image domain has the form of the function J
0(y), shifted according to
$y \to \sqrt {y^2 - L^2} $
dilated by y → 2κ
m
y and blurred by the sinc(Δk
(1,c)(y − y
′)) function. In this first approach the amplitudes S
m
are given by 1/m and the energy of the PSF is normalized to one.
A point not considered by this model is the length of the artifacts in the ISAR image. Since the Bessel function is damped with a factor of y −1/2 for large y, the tail would be visible over the complete image, if the image is depicted in an appropriate dynamic range in dB. In the example of Fig. 2(a), it can be seen that the damping is much stronger than the term of the Bessel function. Because of this, the Bessel function is truncated after a certain reasonable length. An example of a 2D waveguide PSF is depicted in Fig. 3(b).
Fig. 3. Point spread functions for different scattering mechanisms. (a) PSF point scatterer, (b) PSF waveguide scatterer.
III. SPARSE DECOMPOSITION FRAMEWORK
The sparse decomposition framework is based on the assumption that the ISAR image I(x, y) under test can be represented as a sum of an Image I p (x, y) with point-like scatterers and an image I wg (x, y) with waveguide scatterers,
$$I(x,y) = I_p(x,y) + I_{wg}(x,y).$$
This framework is closely related to the one in [Reference Wagner15], where SAR images were decomposed in clutter and target content by a technique called morphological component analysis [Reference Starck, Elad and Donoho16]. To separate the image in the two components, a so called dictionary Φ is needed for each of the components. These dictionaries must represent the corresponding content sparsely and should be inefficient to represent contents of the other component. In the proposed framework, the dictionaries contain the different PSFs of point-like and waveguide scatterers at all possible positions in the image. In the literature on sparse representation, the single functions of the dictionary are called atoms. With these dictionaries both components can be described by a matrix vector multiplication of the dictionary with a corresponding coefficient vector, α:
$$I_p(x,y) = \Phi _p\alpha _p,$$
$$I_{wg}(x,y) = \Phi _{wg}\alpha _{wg}.$$
To form the dictionaries in (16) and (17) the 2D PSFs are vectorized into 1D vectors. Thus the images I p and I wg , marked in the equations as 2D, are also 1D vectors. For the following description of the algorithm, the so-called l p -norm
$$ \vert \vert x \vert \vert_p = \left( {\sum\limits_{i = 1}^N \vert x_i \vert ^p} \right)^{1/p}$$
is used regularly to measure different properties of length N vectors x. For p = 2 this expression is equivalent to the Euclidean norm, the limiting case for p → 0 gives the number of elements not equal to zero and is called the l 0-norm
$$ \vert \vert x \vert \vert_0 = \# \{ x_i \vert x_i \ne 0\}, $$
although it is not a norm in the mathematical meaning. Nevertheless, using this definition, the l 0-norm becomes a direct measure for the sparsity of a signal.
As already mentioned in Section II, the PSFs of both the point-like and waveguide scatterers are normalized to an energy (l 2-norm) of one, leading to a dictionary with normalized atoms. With these dictionaries, an optimization problem can be formulated that minimizes the number of coefficients in α, while the error of the representation is below a certain value σ
$$\eqalign{& \mathop {\min }\limits_{\alpha _p,\alpha _{wg}} \vert \vert \alpha_p\vert \vert_0 + \vert \vert \alpha_{wg}\vert \vert_0\quad s.t. \cr & \quad \vert \vert I - \Phi _p\alpha _p - \Phi _{wg}\alpha _{wg}\vert \vert_2 \lt \sigma .} $$
The problem stated in (20) is NP-hard and thus not efficiently solvable, but several algorithms have been proposed over the last years using either a greedy approach or a convex relaxation of the problem [Reference Tropp and Wright17]. In this paper, a very basic greedy approach called matching pursuit [Reference Mallat and Zhang18] is used to perform the sparse decomposition and is implemented by matrix vector multiplications. In the first iteration of this algorithm, the input image is correlated with all atoms in the dictionary and the maximum correlation index is used as coefficient in α. Then the selected atom is weighted with the correlation index and subtracted from the input image, leading to the residual of the current iteration. After that, this procedure is repeated with the residual instead of the input image until a stopping criteria is met. As stopping criteria, a maximum number of iterations or the approximation error σ can be used.
IV. EXPERIMENTAL RESULTS
In this section, the decomposition result of a real ISAR image, measured by the TIRA of Fraunhofer FHR, is presented. The ISAR image of a small business jet is depicted in Fig. 2(a), the shape of the target can be seen in Fig. 2(b).
For the model, the diameter of the engine and the assumed length of the waveguide structure are needed. In the first experiment, an opening diameter of 90 cm and a length of 30 cm are used to generate the model. The length of the waveguide is approximated by the distance from the entrance to the first compressor stage. As mentioned at the end of Section II, the length of the Bessel function is important to model the damping in range. To consider this damping, the length of the Bessel function is limited to 4.5 m. This cut of the Bessel function is necessary because the model used in Section II assumes a perfectly conducting material for the waveguide. In reality the material is not a perfect conductor and the corrective term in (4) must be extended with a real exponential term that represents the damping and the decay of the artifacts in the image.
The result of the decomposition can be seen in Fig. 4 in the form of detected coordinates for the two different scattering mechanisms. The coordinates of the detected waveguide scatteres are labeled by triangles in the plots on the right side. The matching pursuit algorithm was stopped after 150 iterations and in Fig. 4 the detected coordinates are shown after 50, 100 and 150 iterations. The last two images thus represent the final result.
Fig. 4. Detected point-like (left) and waveguide scatterer (right). (a) 50 iterations, (b) 50 iterations, (c) 100 iterations, (d) 100 iterations, (e) 150 iterations, (f) 150 iterations.
In the images on the right side of Fig. 4, the detected waveguide scattering centers clearly form two pointclouds at the position of the jet engines. Two false alarms are visible in the center of the target, but with a clustering algorithm an automatic detection and location of the engines should be possible. The two false alarms are marked by a red circle in Fig. 4(f).
V. CONCLUSION AND FUTURE WORK
The sparse decomposition of ISAR images is a promising tool for feature extraction and a successful example of a model based signal processing approach. In this paper, we validated the model presented by Borden and showed a successful detection of jet engines. Nevertheless, this work can only be seen as a first step toward feature extraction via a model-based sparse decomposition; there remain many open questions and challenges for the future. One challenge is the shift from the image domain to the raw data domain, and thus incorporate the feature extraction into the imaging process. A framework that has the capability to fulfill this task was presented in [Reference Samadi, Cetin and Masnadi-Shirazi19]. Inside this framework, the mentioned damping term might be included in the model and thus give a better representation of the reality. Another open point is the use of more sophisticated algorithms to perform the sparse decomposition, like the morphological component analysis in [Reference Wagner15].
Simon Wagner received the B.Eng. and the M.Sc. degree in Electrical Engineering from the University of Applied Sciences in Trier, Germany, in 2010 and 2013, respectively. He joined the Fraunhofer FHR in 2012 for his Master Thesis and since 2013 he is a research scientist in the NCI group of the cognitive radar department. His research interests include deep learning and sparse representation for target classification purposes. He is currently pursuing the Ph.D. degree at the University of Siegen in Germany.
Joachim Ender held the positions as the Director of the Fraunhofer-Institute for High Frequency Physics and Radar Technology FHR at Wachtberg, Germany, and as Chair for High Frequency Sensors and Radar Techniques at the University of Siegen until July 2016. After his retirement he is still active at both institutions as senior scientist. After his diploma in mathematics/physics he performed research on various topics of radar science since more than forty years. After the start of his carrier as young scientist in 1976 he became the Head of a department and in 2003 the Director of the FHR. Joachim Ender is author and co-author of numerous papers. Among other prizes, he received the “Group Technical Achievement Award – For contributions to Array Signal Processing and Multichannel Synthetic Aperture Radar” from EURASIP. In 2014, he was named as an IEEE-fellow “for contributions to multi-channel synthetic aperture radar and radar array signal processing”. He was one of the founder members of the “European Conference on Synthetic Aperture Radar” (EUSAR), which takes place every two years since 1996.
