Inverse Synthetic Aperture Radar Imaging: A Historical Perspective and State-of-the-Art Survey

Radar can obtain high spatial resolution at large operating distances, and its operation is largely independent from weather and lighting conditions. This makes it an important instrument in an ever-increasing number of applications requiring robust and reliable sensing of the environment. Due to the advancements in semiconductor technology and digital computation power, radar can obtain increasingly better position estimation accuracy, detection sensitivity, and target resolution. High-resolution imaging is an example of a rapidly advancing field of radar technology. Inverse synthetic aperture radar (ISAR) is an imaging technique that uses the motion of a target to obtain a high-resolution radar image of it. This paper presents a comprehensive review of past work and the state of the art in ISAR imaging systems, algorithms, and applications. Starting from the early historical developments of simple monostatic range-Doppler imaging, we present the developments up to the current advanced spatially distributed imaging modalities and reconstruction algorithms. The list of references includes the most influential ISAR publications in the open literature up to early 2021.


I. INTRODUCTION A. BACKGROUND
Radar has several benefits compared to more traditional imaging sensors such as optical, infrared, or hyperspectral cameras. First, radar can provide a very high spatial resolution independent of the distance to the imaged target scene since the resolution only depends on the bandwidth, wavelength, and possible relative movement between the radar and the target. Second, radar can operate during day or night since it is an active sensor that is not dependent on daylight. Moreover, the electromagnetic waves at typical radar wavelengths can penetrate many visual obstacles, such as clouds, fog, or smoke, making the image quality largely independent of weather conditions.
High-resolution radar imaging using the synthetic aperture principle is a powerful tool for a large number of surveillance and environmental monitoring applications [1]. Synthetic aperture radar (SAR) is an imaging technique that uses the motion of the radar platform to create a large synthetic The associate editor coordinating the review of this manuscript and approving it for publication was Cheng Hu . antenna corresponding to a narrow beamwidth, i.e., a high angular (cross-range) resolution [1]- [3]. In SAR imaging, the radar typically collects the data while moving on board an air-or spaceborne platform and then focuses the data collected along the flight path to form a high-resolution range-cross-range image of the ground surface.
A high resolution in range is easily obtained by using a suitable high-bandwidth waveform. To get a focused target image with high cross-range resolution, a very narrow antenna beam is required. At large operating ranges, the diameter of a real aperture antenna would have to be impractically large to achieve this. Thus, an alternative solution-the so-called synthetic aperture principle-has to be adopted. By allowing the relative position between the target scene and the radar to change as a function of time, the received signal can be processed as if it had been received by a large synthetic array (the diameter of which is determined by the platform trajectory). This enables the creation of a very narrow antenna beam in the cross-range direction by digital signal processing after the signals have been recorded.
To create the synthetic aperture, only a suitable relative motion between the radar and target scene is required. Thus, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ either the radar platform or the target (or both) needs to move during the data collection. When the motion of a moving target (instead of just the radar with regard to the ground surface) is used to create the synthetic aperture, the imaging technique is called inverse synthetic aperture radar (ISAR). Depending on the sensor characteristics and the relative motion, the result of an ISAR system is a 1-, 2-, or 3-D description of the target's spatial electromagnetic reflectivity distribution. The radar images obtained by ISAR systems are commonly used to recognize unknown detected targets when cooperative identification is not possible. Alternatively, ISAR can be used as an additional independent information source even when cooperative identification is used. Non-cooperative target recognition (NCTR) techniques are used in maritime, airspace, near-space, and overland surveillance applications to identify targets such as vessels, aircraft, ballistic missiles, satellites, and ground vehicles.
Compared to SAR, the applications for ISAR are much more limited. The required signal processing algorithms in ISAR are often more complicated due to unknown target motion. Additionally, SAR can be used to image and monitor large areas with multiple targets, whereas ISAR is usually limited to imaging one target at a time.
From a historical perspective, military applications have been the driving force behind NCTR and the related ISAR techniques. However, an ever-increasing number of civilian applications is currently emerging. For example, ISAR has recently been used in automotive [4], [5] and space surveillance [6]- [8] use cases. The trend of ISAR imaging technology is moving toward higher spatial resolution and more advanced imaging modalities utilizing spatially distributed systems.

B. SCOPE AND STRUCTURE OF THE PAPER
The purpose of this paper is to present a comprehensive review of past developments and the state of the art in ISAR imaging. During recent years, specialized textbooks concerning general ISAR theory and applications have been published [9]- [12]. These references provide a detailed theoretical description and demonstrate several applications of ISAR imaging. This paper supplements these references by providing a more complete description of the historical developments of ISAR and a more extensive survey concerning the various research topics in the open literature. We focus on ISAR systems and algorithms, which aim to produce a radar image of the target. We do not cover techniques that use these images in a further application, such as NCTR (or automatic target recognition, ATR) and polarimetric ISAR.
We begin the article by introducing the basic ISAR terminology and a simple signal model in Section II. The presented material aids in understanding the significance of the various algorithms and concepts described in the later sections. Readers familiar with basic ISAR theory and terminology may skip this introductory section.
The early historical developments in ISAR imaging are reviewed in Section III. Starting from the first applications in range-Doppler (RD) imaging of the moon and the planets, we survey the early applications up to the introduction of turntable ISAR and imaging of air targets in the early 1990s.
Section IV reviews the developments in monostatic ISAR processing algorithms from the mid-1990s until today. The conventional assumption of dividing the target motion into translation and rotation as well as various motion compensation methods are reviewed.
Since the early 2000s, spatially distributed (bi-and multistatic) ISAR systems and the related processing algorithms have been of great research interest. Section V reviews the seminal work in these areas. Important concepts and applications such as passive bistatic ISAR and 3-D interferometric ISAR (InISAR) are surveyed.
Section VI surveys important recent applications of ISAR imaging. Ground moving target imaging from space-and airborne platforms is one use case receiving a great deal of interest [13]. Another emerging research topic is imaging non-rigid targets with rotating parts (such as drones and helicopters). Additionally, a number of space surveillance radars with ISAR imaging capabilities have been developed recently.
The advances in computer science and signal processing can also be seen in the current trends in ISAR research. Section VII reviews the most recent algorithmic developments in ISAR, focusing on compressed sensing (CS) and machine learning (ML) techniques.
Finally, Section VIII concludes the paper and outlines the future trends in ISAR research and applications.

II. BASIC PRINCIPLES OF ISAR A. BASIC DEFINITIONS
The term ISAR is defined in the IEEE standard radar definitions as [14] ''An imaging radar in which cross-range resolution (angular resolution) of a target (such as a ship, aircraft, or other reflecting object) is obtained by a synthetic aperture formed by the rotation or translation of the target, as opposed to motion of the radar. Cross-range resolution of target whose exact angular motion is unknown (e.g. non-cooperative targets such as ships at sea) is often achieved by resolving in the Doppler domain the different Doppler frequencies produced by echoes from the individual parts of the object.'' The basic ISAR operating principle for a monostatic system can be summarized as follows. First, the radar transmits a sequence of electromagnetic signals-in the form of a suitable high-bandwidth waveform-toward the imaged target. Commonly used waveforms include a pulsed linear frequency modulation (LFM), the stepped frequency waveform, or frequency modulated continuous wave (FMCW). Then, the radar receiver records the echoes from the target that are scattered toward the receiver. By performing range (fast-time) matched filtering (pulse compression) for the received signal, a series of high-resolution range profiles is obtained as a function of slow-time t (or pulse number m). Slow time is a time variable that remains constant during the formation of a single range profile and increases from one range profile to the next. In ISAR processing, the target is commonly modeled as a superposition of isolated and isotropic point scatterers (Born approximation). Figure 1 illustrates the basic ISAR geometry for a monostatic setup where the imaged target is constrained to move on a 2-D plane. The primed coordinate system-the origin of which is assumed to be located in the target's center of mass-is rigidly attached to the target. The spatial degrees of freedom for each instant of slow time include the line-ofsight vector r 0 , i.e., the location of the target's center, and the orientation of the primed coordinate system. The orientation is determined by the angle between the target velocity vector v and the line-of-sight vector r 0 . For simplicity, here we assume the target always to be oriented along v. Importantly, the target aspect angle θ changes as the object moves in a suitable manner. Consequently, the range R p (t) = ||R p (t)|| of every spatial position on the target changes in a unique way as a function of slow time. This creates the different Doppler frequencies for each scatterer required to resolve them in cross range.

B. SIGNAL MODEL
In practical ISAR scenarios, the targets are objects such as cars, ships, satellites, or aircraft. The distance between the radar and the target is usually several kilometers, whereas the dimensions of the targets are much smaller. Thus, we can make a far-field approximation r 0 = r 0 r p = r p . This can be used to derive an expression for the distance R p between an arbitrary point on the object and the radar as [15] R p (t) ≈ r 0 (t) + x p cos θ(t) + y p sin θ(t), where x p = r p sin θ p , and y p = r p cos θ p . This is a key result because the motion can be divided into two parts: translation and rotation. The translation r 0 (t) is the same for every target point p (spatially invariant). Thus, it is not useful for ISAR imaging since it gives each scatterer the same term in the Doppler frequency where λ is the carrier wavelength. The rotational motion θ(t) produces the unique (spatially variant) Doppler frequencies for different target points p required to obtain the cross-range resolution.
In both SAR and ISAR, the so-called start-stop approximation is commonly employed. This means that we assume the target to be stationary during the formation of a single range profile. Assuming that the radar is illuminating a single point target at (x p , y p ), the range profile can be expressed mathematically as where B is the frequency bandwidth of the signal, c is the speed of light, and σ is the radar cross section (RCS) distribution of the target. Furthermore, a is the waveform's autocorrelation function. For example, in the case of LFM, we have a(r) = sinc(r). Using the Born approximation, we can write the reflectivity function of the target as a sum of P point scatterers as where δ is the Dirac delta function. As the target moves as a function of t, the distance R p of each scatterer changes. Furthermore, the range profile of the target reflection can be expressed as a convolution between the point target response (3) and the target reflectivity function g. Mathematically, the ISAR signal model can thus VOLUME 9, 2021 be expressed as where is the rectangle function, and T is the length of the coherent processing interval (CPI). It is helpful to consider the signal model in the range spatial frequency k r domain by taking the Fourier transform of (5) with regard to the range variable r. Assuming that a(r) = sinc(r), using the convolution theorem, and plugging in R p from (1), we get where k c = (2π)/λ, and G(k x , k y ) is the two-dimensional Fourier transform of the reflectivity function g. The expression in (6) is a key result for both ISAR and spotlight mode SAR [3], [16]. It can be interpreted as follows: the range-compressed radar signal is a series of phase-modulated slices of the 2-D Fourier transform of the reflectivity function g. The rectangle function in (6) represents the fact that the signal is band limited in k r (which is an approximation, since the range profile is assumed to be of finite support). Consequently, the slices span an annulus sector in the 2-D spatial frequency domain. To successfully reconstruct the target image g, we need to compensate the phase modulation caused by r 0 and know the angles θ(t) to calculate the inverse Fourier transform.
In non-cooperative ISAR, difficulty arises from the fact that both r 0 and θ are unknown a priori. To reconstruct a properly focused image, their values as a function of t have to somehow be estimated. It is also evident from (6) that r 0 is a nuisance that does not provide any useful information for the image reconstruction. Its estimation is referred to as translational motion compensation. Correspondingly, the estimation of θ is referred to as rotational motion compensation.

C. IMAGE RECONSTRUCTION AND MOTION ESTIMATION
The results in (5) and (6) aid in understanding the ISAR image reconstruction process. The problem is to estimate the spatial reflectivity distribution g(x, y) given the received signal (5). Essentially, we need to apply a slow-time-dependent matched filter h for each image position to recover g. Mathematically, this can be expressed as where I denotes the ISAR image. Ideally, the filter h is the impulse response of a point target located at position (x, y). Thus, the phase behavior of h is determined by the relative motion (r 0 and θ). For non-cooperative targets, a prerequisite for calculating (7) is target motion estimation. Since the matched filtering is based on coherently combining the range profiles, the motion has to be estimated very precisely (within a fraction of the wavelength λ) to obtain well-focused images.
In SAR imaging, very accurate information about the relative motion is usually provided by a global positioning system and inertial sensors. However, since in non-cooperative ISAR the target motion is unknown and unpredictable, the estimation problem is much more difficult.
In the 2-D model described above, we have two degrees of freedom (r 0 and θ) for each instant of slow time. In case the target moves in 3-D, additional degrees of freedom need to be introduced. In general, these are two additional angles that describe the orientation of the primed coordinate frame (now including a third coordinate z ) in the unprimed frame. Because the number of range profiles M is typically in the order of hundreds, the total number of unknown degrees of freedom 2M is often prohibitively large to allow a computationally feasible solution.
Thus, many assumptions and approximations have to be made in order to successfully estimate the target motion. These are application dependent, and no successful general method to solve the motion estimation problem has thus far been presented. This has led to a large number of different motion compensation and image reconstruction algorithms in the ISAR literature. The different methods contain various trade-offs between estimation accuracy (robustness), image resolution (information content), computational requirements, and system complexity.
The simplest approach to ISAR is the RD algorithm, where the Fourier transform is used as the matched filter for each range bin of the signal (5). However, this method is not often applicable for high-resolution imaging since it requires many assumptions to be satisfied. First, the range R p of each scatterer should not change more than the range resolution r = c/(2B) during the CPI. Otherwise, scatterers migrate from one range bin to another, causing smearing in range and a reduced Doppler resolution. Moreover, for the Doppler frequency f p (see (2)) of each scatterer to be constant, we need to make a small angle approximation, and the aspect angle needs to change linearly.
Whereas the RD technique represents the most simple and crude reconstruction method, back projection (BP) is the most accurate and general way to reconstruct the image. In BP, the point target response (4) is used as the matched filter h for each image position (x, y). The difficulty in using BP lies in the fact that the motion parameters r 0 and θ need to be known a priori to construct the matched filter h. Moreover, calculating the 2-D correlation integral (7) for each image pixel position (x, y) is computationally very intensive. This makes it very difficult to incorporate a motion estimation approach with BP.
The flow chart in Figure 2 summarizes the typical ISAR processing steps described in this section. The model in (5) represents the signal after range compression. For clarity, the flow chart depicts motion compensation and image reconstruction as separate steps. However, in some cases they are performed jointly, and it is not possible to distinguish them into separate processing steps.   The plot on the left shows signal intensity after range compression (range profiles). The movement of the car during the CPI manifests itself as the changing location of the intensity distribution inside the range gate. The plot on the right shows the reconstructed ISAR image with a spatial resolution of approximately 15 cm. In this example, the radar and the car are located on the same ground plane. For this reason, the car corner closest to the radar produces the strongest reflection, whereas the opposite corner is shadowed.

D. SPATIALLY DISTRIBUTED SYSTEMS
Bi-and multistatic ISAR systems have spatially separated transmitter(s) and receiver(s). This setup can be used to obtain 3-D images of targets, more useful information about the target's reflectivity, and more degrees of freedom for the motion estimation. Moreover, they can be used as passive ISAR sensors. Importantly, spatially distributed ISAR is more robust against unfortunate imaging geometries-for a monostatic sensor, a radially moving target is impossible to resolve in cross range (since the Doppler shift of each scatterer is the same).
Even though the above-described signal model only holds for a monostatic geometry, the basic concepts are essentially the same for the spatially distributed case. A detailed analysis of the spatially distributed ISAR signal model is outside the scope of this paper. References [17], [18] provide a detailed description of the bi-and multistatic ISAR models, respectively. Spatially distributed ISAR methods are surveyed and analyzed in more detail in Section V of this paper.

III. EARLY HISTORICAL DEVELOPMENTS
We begin our survey by providing a historical perspective of ISAR imaging. This section surveys the early ISAR developments from the late 1950s up to the early 1990s. The important early milestones illustrated in the timeline of Figure 4 are described next.

A. RADAR ASTRONOMY
The seminal paper by Ausherman et al. [19] summarized the earliest ISAR developments concisely. The concept was formulated for the first time in the 1950s in the context of radar astronomy [20] with the aim of improving the resolution of the radars measuring the surfaces of the moon and the planets of the solar system. In these cases, the rotation of the planets and the moon was used to obtain the differential Doppler resolution capability [21]- [23], though the term ISAR was introduced only later by the radar community.
The first radar images of the moon were obtained by the Millstone Hill radar operating at 440 MHz center frequency [24]. The images had a range resolution of 75 km and a Doppler resolution of 0.1 Hz, obtained with a CPI of 10 seconds. During the following years, the same concept was used to obtain images of different planets. The Millstone radar was used to acquire images of Venus and Mercury [25]- [27]. The Arecibo telescope was used to measure and analyze the Doppler spectrum of Mercury, Mars, and Jupiter [28]- [31]. However, the Arecibo instrument was not able to properly resolve the planets in range due to its lowbandwidth waveform. Thus, the Doppler spectra were used to estimate the planets' rotation parameters, not to map their surfaces. Nevertheless, the Doppler analysis technique used in this case was essentially based on the ISAR principle.

B. SPACE TARGET IMAGING
The next significant developments occurred in the late 1960s, when it was recognized that the RD approach could be applied to imaging of satellites orbiting the Earth in space. For the purpose of space target identification, a radar operating at 94 GHz with a bandwidth of 1 GHz using LFM was developed by the Aerospace Corporation [32]. The first radar images of space targets were obtained by the Advanced Research Projects Agency (ARPA) Lincoln C-band Observables Radar (ALCOR) in the early 1970s [33]. The range resolution of ALCOR was initially 50 cm, and the system was later developed-due to the success of early ALCOR results-to have a higher pulse repetition frequency (PRF) and more robust data acquisition for imaging purposes [19], [33].
In the 1970s, the next significant advancement in this field was the development of the long-range imaging radar (LRIR) by the Lincoln Laboratory [34]. LRIR had a range resolution of 25 cm and a PRF of 1000 Hz at X-band, allowing for high-quality imaging of rapidly rotating space targets. In this context, significant progress was also made in the signal processing algorithms used to reconstruct the target images, leading to increased spatial resolution. Concepts for using longer CPIs were developed to achieve wide-angle and 3-D imaging capabilities [19], [34].

C. TURNTABLE IMAGING
The first developments for imaging rotating ground-based targets began in the 1960s at the Willow Run Laboratories. Scientists recognized that the concept was essentially equivalent to SAR imaging, which had been introduced during the previous decade [35]. The first experiments were carried out using a rotating platform (turntable) and a coherent ground-based radar system. The turntable setup provided a well-controlled environment, allowing for the development and demonstration of increasingly advanced imaging algorithms. Initially, the image reconstruction was based on the Fourier transform and a simple procedure to compensate for the scatterer migration through range resolution cells (MTRC) [36].
A more general formulation of the rotating target imaging problem was given by Walker [37], [38]. He introduced the polar format data storage and processing concept, which solved the scatterer MTRC problem. MTRC occurs when the scatterer ranges R p change more than the range resolution r during the CPI, which is common in highrange resolution and wide-angle imaging. The algorithms and the results of the turntable experiments were reported in the seminal paper [38]. Early theoretical developments of ISAR were also reported in [39], [40]. At about the same time, Mensa et al. [41], [42], Wehner et al. [43], and Chen et al. [44], [45] also reported their first concepts and results in RD imaging of rotating targets.

D. OTHER EARLY DEVELOPMENTS
In the early 1970s, Raney demonstrated imaging ground moving targets using an airborne SAR system [46]. The concept was further developed and demonstrated in the 1980s [47], [48] and early 1990s [49]- [51] and also applied to imaging maritime targets [52]. Motion estimation and compensation for non-cooperative targets was achieved using the dominant scatterer technique (i.e., estimating the phase of prominent point scatterers on the target as a function of slow time) [2], [49] or by using time-frequency analysis techniques [50], [51].
In the late 1980s, the first results from imaging airborne targets with ground-based ISAR systems were reported. Steinberg [53] presented X-band images of an aircraft with a 1-m resolution. The first concepts related to estimating the target motion without a priori knowledge were also presented for this purpose [54], [55]. They were also based on estimating the phase of multiple dominant scatterers on the target.
Since the early 1990s, the number of reported applications for ISAR has steadily increased, with a trend toward higher spatial resolution and more sophisticated imaging modalities. A difficult challenge in these early applications was to image maneuvering targets with non-linear motion. A significant development to address this problem was the introduction of time-frequency techniques for motion compensation and imaging. A detailed description of these methods is presented in Section IV-B.
The earliest applications for imaging ground-based moving targets from air-or spaceborne platforms as well as imaging space targets at the low Earth orbit region remain very active application and research areas. This is due to the need for active space surveillance to prevent space debris collisions with active satellites [8]. Modern applications and developments in these areas are reviewed in more detail in Section VI. The important later developments described in the next sections are summarized in the timeline seen in Figure 5.

IV. MONOSTATIC ISAR PROCESSING
Following the developments described in the previous section, the motion compensation and image reconstruction algorithms required for successful ISAR imaging have received a great deal of interest in the open literature. Until the early 2000s, most of the systems and methods were developed in the context of monostatic ISAR (all of the early systems described in the previous section were VOLUME 9, 2021 FIGURE 6. Methods used for translational motion compensation in conventional monostatic ISAR imaging. The two-step approach is computationally more tractable but starts to break down for long CPIs with significant target rotation. monostatic). This section reviews the various ISAR motion compensation and image reconstruction algorithms that have been developed and demonstrated for monostatic ISAR systems. It should be noted that many of these methods can be used for bi-and multistatic processing as well-albeit usually with some modifications.
Instead of going through the developments in chronological order, we review specific problems and the methods used to solve them in a sequential manner. Each of the next two subsections deals with a distinct aspect of ISAR processing. In Section IV-A, we first consider the methods to estimate and compensate the translation r 0 of the target. Then, we proceed to review the approaches used for estimating the target rotation θ and performing the image reconstruction in Section IV-B.

A. TRANSLATIONAL MOTION COMPENSATION
Translational motion compensation is commonly divided into two parts. The first part performs an initial coarse estimation and is called range alignment. It aims to compensate the target translation to an accuracy comparable to the range resolution r, i.e., to ''eqnarray'' the range profiles. It first estimates r 0 (t) and then shifts the range profiles to compensate for the changing range. The second part, called autofocus, performs a fine estimation to accurately remove the residual spatially invariant phase errors remaining after range alignment.
The division into two steps is usually done for computational reasons. However, it is also possible to estimate and compensate the translation using a single-step approach. Figure 6 summarizes these two classes of different methods for translational motion compensation. Next, we review the specific techniques in more detail.

1) RANGE ALIGNMENT
The first and simplest range alignment method is based on choosing a reference range profile (usually corresponding to the first pulse of the CPI) and calculating the cross-correlation between it and the other range profiles to determine the range shifts [53], [56]- [58]. However, the accuracy of this method suffers if the target reflectivity distribution (RCS) changes during the CPI due to rotation. A way to avoid this problem is to correlate each adjacent range profile (and thus estimate the differential range shifts). However, this method is vulnerable to error accumulation since the final estimation is obtained by integrating the estimated differential range shifts. More recently, extensions to the cross-correlation method have been developed to partially overcome the above-mentioned problems [59], [60].
Methods based on the Hough transform [61] and keystone formatting [62]- [65] have also been successfully used for range alignment. Of these two methods, keystone formatting also plays an important role in rotational motion compensation, as described later. Both of these methods assume a simple model (first-or second-order polynomial) for r 0 as a function of slow time. Keystone formatting is based on eliminating the RD coupling by re-scaling the slow-time variable in the (k r , t) domain. The Hough transform is a wellknown method that can be used to find specific shapes in images. In this case, it is used to find the first-or secondorder curves traced by the target in the range-compressed (r, t) domain.
Several other range alignment methods are based on optimizing a quality measure calculated from the amplitude of the range profiles. The minimum entropy method [66]- [68] is based on minimizing the entropy of the sum of two adjacent profiles as a function of the range shift applied to one of the profiles. The global range alignment method [69], [70] estimates the translation by minimizing a loss function, the value of which quantitatively measures the quality of the range alignment for the entire CPI. Several different quality measures have been proposed [71], [72], such as the contrast and entropy of the average range profile and the sum of the squared differences of adjacent range profiles. The global method uses the entire CPI to estimate the motion, whereas the minimum entropy method operates on a pulse-to-pulse basis. A disadvantage of the global method is an increased computational complexity, which means that usually a parametric model needs to be used for r 0 (t). The minimum entropy method is essentially nonparametric. Methods to speed up the computation of the global method were proposed in [72], [73]. An efficient hybrid between the global and local pulse-to-pulse methods uses the average range profile (i.e., data from the entire CPI) in the loss function but does the estimation on a pulse-topulse basis [74]- [76]. This method is more robust when the target exhibits rapid motions not easily captured by a loworder parametric model.
Currently, the most important limitations of range alignment methods are related to poor performance under low signal-to-noise ratio (SNR) and sensitivity to significant changes in the target RCS during the CPI. The optimizationbased approaches suffer the most in the presence of low SNR because the quality metrics only use the absolute value of the range-compressed signal-i.e., a possible coherent integration gain is not utilized. An underlying basic assumption behind all of the methods is that the target RCS is constant during the CPI. In practice, this means that the target should not rotate much during the CPI. The global method provides the most robustness against target rotation, whereas the cross-correlation method is the most sensitive to it. Table 1 summarizes the most important strengths and weaknesses of the range alignment methods described above.

2) AUTOFOCUS
After range alignment is performed, a slow-time-dependent phase error motion may remain in the data. This is because range alignment is essentially a non-coherent procedure: it only takes into account the amplitude envelopes of the signal. Moreover, some sources of phase errors (e.g., target vibrations or atmospheric distortions) may be so small that they do not cause a noticeable range shift while still significantly distorting the slow-time phase-necessitating an additional phase compensation step.
In the SAR (and ISAR) literature, autofocus refers to the residual translational motion compensation that is applied to each range bin signal after range alignment. Autofocus has been extensively studied in the context of airborne spotlight mode SAR [2], [3] and synthetic aperture sonar imaging [77] since the problem is essentially the same as in ISAR. Two well-known classes of autofocus algorithms are based on image sharpness maximization (a general term that includes contrast optimization and entropy minimization) [67], [78]- [82] and phase gradient autofocus (PGA) [83], [84].
In sharpness maximization autofocus, the objective function is calculated from the image intensities. For example, it can be the coefficient of variation of the image intensities or amplitudes [79], the sum of the squared intensities [82], or the entropy of the intensity image [67]. The sharpness maximization autofocus can be considered a local optimization problem [80]- [82]. Thus, it can be solved efficiently by first-or second-order local optimization algorithms [15], [82], [85]. Since the required numerical optimization is generally a computationally intensive task, several approaches have been proposed to improve the efficiency [15], [75], [85]- [93]. One option is to use a parametric model for the phase errors as a function of slow time. Additionally, different image sharpness metrics have been analyzed and compared [94]- [96].
PGA uses an iterative procedure to estimate the unknown phase error [3], [84]. After isolating dominant scatterers in the image (range-cross-range) domain, the signal is transformed back into the range-slow-time domain to perform the phase error estimation. The process is repeated iteratively since it becomes easier to isolate the point scatterers when the image focus increases. PGA is computationally more efficient than sharpness maximization. The conventional PGA [84] was limited to Fourier transform-based image reconstruction using a single strong scatterer per range bin for estimation. For this reason, extensions to PGA have been considered since its conception in the 1990s. More accurate and robust estimators were proposed in [97]- [105]. The most recent advances [106], [107] proposed an extension of PGA to allow its application for diverse imaging geometries and different image reconstruction methods.
Another well-established autofocus method called multichannel autofocus is based on a linear algebraic formulation of the phase error estimation problem [108]- [110]. It uses an assumption that the image support is limited (when the image is correctly focused, there is a region of approximately zero intensity pixels) to estimate the phase errors. The above-mentioned references have demonstrated its superior focusing performance compared to PGA and sharpness maximization for spotlight mode SAR images.
Sharpness maximization works well when the ideally focused image has a high contrast, i.e., one or multiple strong scattering centers on a weak background. However, when the ideally focused image has a relatively uniform intensity distribution, sharpness maximization may produce erroneous results. In this case, PGA is usually more robust. Even though PGA relies on isolating strong scatterers, it averages the estimation result from each range bin to obtain the final result. In most ISAR applications, the target can be considered a sum of point-like reflections with weak background (clutter) reflections. Thus, both sharpness maximization and PGA are widely used solutions in ISAR. The relative strengths and weaknesses of the above-mentioned autofocus algorithms are summarized in Table 2.

3) SINGLE-STEP COMPENSATION
The above-described two-step approach to translational motion compensation is commonly used due to its computational efficiency. However, it has some drawbacks. Most importantly, the autofocus approach requires successful range alignment (the residual error has to be smaller than the range resolution). If the range alignment fails to adequately compensate the translation, autofocus cannot remedy the situation. Thus, the two-step approach is ultimately limited by the accuracy of the range alignment algorithm. As mentioned previously, long CPIs present a problem even for the most robust range alignment methods. Increased cross-range resolution requires a longer CPI (larger target rotation), limiting the applicability of the two-step approach in very high-resolution ISAR imaging.
To overcome these limitations, several methods requiring only a single step of translational motion estimation and compensation have been considered in the literature. An early example of such a method is the centroid tracking approach [111]. In this approach, the translation is compensated by tracking and making the target centroid (defined as the weighted average of the target range and Doppler) constant as a function of range and slow time. Another early method was based on optimizing a burst derivative quality measure, which describes the image sharpness [112], [113]. Compared to sharpness maximization autofocus, the advantage of this method is that image reconstruction does not have to be performed to evaluate the burst derivative metric.
The sharpness maximization concept can also be used as a single-step approach. In this case, the estimation is extended from the residual slow-time phase to the full translation r 0 . The main drawback of this approach is the resulting computational intensity-the autofocus problem can be considered a 1-D (cross-range) focusing problem, whereas now the focusing needs to be done in 2-D (range and cross-range). Berizzi et al. [79] optimized the coefficients of a parametric translation model using the image contrastdefined as the ratio between the standard deviation and the mean of image intensities-as the objective function. The same approach has since been successfully used with different optimization algorithms, objective functions, and target motion models [114]- [119].
An approach based on the CLEAN algorithm, which uses multiple independent scattering centers on the target to estimate the translation, was proposed in [120]. More recently, the Doppler centroid estimation method was extended (to include estimating the Doppler rate) in [121]. In [122], a method based on using an improved version of the cross-correlation range alignment method with phasederived velocity estimates was used to perform single-step compensation.

B. ROTATIONAL MOTION COMPENSATION
After the translation r 0 has been compensated, the target rotation θ needs to be estimated to properly perform image reconstruction. The rotation determines the phase progression of the slow-time matched filter h. The process of estimating the unknown rotational motion and using it in image reconstruction is generally more complicated than translational motion compensation. Many different methods for performing slow-time filtering and estimating the target rotation have been proposed in the ISAR literature.
The most common approach is to divide the problem into two distinct parts: image reconstruction filtering and crossrange scaling. The purpose of the first part is to design such a filter that the target reflections are properly compressed into point-like responses in the ISAR image domain. Since the slow-time signal can be considered a superposition of signals with slow-time varying Doppler frequencies f p , methods from time-frequency signal processing are commonly used for this purpose.
However, the problem with time-frequency filtering methods is that they only estimate the Doppler frequency of the scatterers. The Doppler frequency is directly proportional R. Vehmas, N. Neuberger: ISAR Imaging: Historical Perspective and State-of-the-Art Survey to the cross-range position, but there is an unknown scale factor depending on the aspect angle change (target rotation) during the CPI. Moreover, the direction of the target rotation determines the image projection plane. If the target is not constrained to move rigidly on a 2-D plane (see Figure 1), knowing the projection plane is crucial for interpreting the image and estimating the size of the target.
Methods for obtaining an estimate for the rotational angle are called cross-range scaling since the Doppler axis can be properly scaled into spatial cross-range units when the rotation is estimated. This step is crucial for image interpretation and exploitation in NCTR applications. Figure 7 summarizes the methods used for the abovementioned two steps (image reconstruction filtering and cross-range scaling) into different classes. Next, we review these algorithms in more detail.

1) RANGE-INSTANTANEOUS-DOPPLER IMAGE RECONSTRUCTION
The simplest way to perform image reconstruction is to use a 1-D Fourier transform, i.e., to assume a constant f p (linear slow-time phase progression) for each scatterer during the CPI. However, when the range resolution is extremely high and a similarly high resolution is desired in the crossrange direction, this simple RD algorithm rarely produces an acceptable result. As the CPI length and the rotational angle increase, scatterer MTRC and the time dependence of f p cause defocusing and thus degraded image resolution (in both range and cross-range dimensions).
Moreover, if the rotational motion is not uniform, f p will not be constant even during a short slow-time interval. As a solution to the MTRC problem, a technique called keystone formatting can be used [15], [62]- [64], [123]- [130]. This approach was initially introduced by Perry et al. [62] in the context of airborne SAR imaging of moving objects. As previously mentioned, it can also be used to achieve range alignment. However, keystone formatting can only compensate for MTRC of a specific order.
Commonly, either the linear or quadratic component is compensated.
To mitigate the effect of slow-time-dependent scatterer Doppler frequencies f p (t) (non-linear phase histories), time-frequency representations (TFRs) can be used as the slow-time matched filter. This approach was first introduced by Chen [9], [131], [132] in the 1990s. Since then, it has matured into a well-established technique [133], [134], and many different TFRs have been successfully used in ISAR imaging. Theoretical analyses of various trade-offs in the TFR-based imaging approach-concerning metrics such as resolution and SNR-were provided in [133]- [137]. To distinguish this approach from the conventional Fourier transform-based RD algorithm, it is usually referred to as the range-instantaneous-Doppler (RID) method. The conventional RD approach uses the entire CPI to reconstruct a single ISAR image, whereas the RID approach produces a set of ISAR images (one for each slow-time instant) for the same CPI.
In an attempt to eliminate the harmful cross terms inherent in the Wigner-Ville TFR, several variations have been successfully used to reconstruct ISAR images [138]- [142]. The so-called adaptive joint time-frequency transform [143] has also been demonstrated in numerous studies [144]- [148]. Furthermore, the S-method [149] has been shown to greatly surpass the conventional RD approach [15], [150]- [153]. The fractional and local polynomial Fourier transforms represent a compromise between the pure RD and RID methods: instead of transforming the signal into the full slow-time-Doppler domain, they only calculate the TFR along a certain line (linear combination of t and f p ) in the TFR plane [154], [155].
More recently, time-frequency signal processing methods such as the synchrosqueezing transformation [156] and timefrequency reassignment [157] have been demonstrated in ISAR applications [130], [158], [159]. These methods have the potential to further improve the image focus and mitigate the effect of sidelobes in the image domain. References [160]- [163] used the chirplet decomposition to perform slow-time filtering. In these methods, the signal in each range bin is considered to be a superposition of ''chirplets,'' which are basis functions exp iφ(t), the slowtime phase history φ(t) of which follows a chosen model (most commonly quadratic or cubic). A number of different methods based on this assumption (superposition of quadratic or cubic phase signals) were developed and demonstrated in [164]- [185]. These methods are based on estimating the Doppler modulation parameters (second-and third-order phase coefficients), which can then be used to estimate the rotation and perform slow-time filtering.
The large number of different RID methods demonstrates the difficulty of formulating a general method for ISAR imaging. Choosing a suitable method largely depends on a priori knowledge about the possible target types and their motion dynamics, RCS distributions, etc. The RID method is more robust and produces a better cross-range resolution than RD in cases where the target exhibits complex motions. However, this robustness comes at the cost of having to tune the parameters of a chosen TFR to a specific use case, limiting general applicability. Moreover, the well-known time-frequency uncertainty relation of TFRs also presents a performance limitation for RID imaging. In other words, a compromise between the CPI length, TFR cross-term suppression, and cross-range resolution is unavoidable in one form or another.

2) OTHER METHODS
Even though the TFR-based RID approach mitigates the effects of non-linear slow-time phase histories (changing f p ), it is not a complete solution to the problem. To mitigate the above-mentioned limitations, a number of other rotation estimation and image reconstruction techniques have been proposed.
A general method for reconstructing the image without any approximations is the BP algorithm. However, it is computationally very intensive when combined with a data-driven motion estimation method (usually including numerical optimization). Thus, BP has so far been mostly limited to imaging of cooperative targets, where the target motion is known a priori [186]- [190].
More recently, BP algorithms with integrated motion compensation have been considered [191], [192]. However, these methods require approximations concerning the target motion to be made to realize a computationally efficient implementation. The BP approach has had more success for a multistatic setup, where the motion estimation is easier due to the increased degrees of freedom provided by the spatial diversity [18]. This method is later described in Section V.
Several non-linear filtering methods have also been used in ISAR. The most popular approach is to use sparse reconstruction algorithms from CS. These methods, which have seen significantly higher use in recent years, are reviewed in Section VII-A. Other examples of applying non-linear slow-time filtering in ISAR include the Capon method [193], the iterative minimum means squares error method [194], the weighted least squaresbased iterative adaptive approach [195], and super-resolution techniques [196]. These methods have been shown to surpass the RD and RID methods in terms of image resolution in specific examples. However, their increased performance comes at a cost. The methods require a high SNR to work well, the target reflectivity often needs to follow a specific model, and the computational cost is high. Table 3 summarizes the strengths and weaknesses of the different image reconstruction filters described above. Whereas the RID approach remains the most widely used concept due to its computational efficiency, BP-and CSbased methods are becoming more attractive as the available computation power continues to increase.

3) TIME WINDOW OPTIMIZATION
In a typical operational scenario, radar measures the noncooperative target as long as it remains in the field of view. This often results in a large amount of recorded data, and thus it becomes necessary to choose a suitable subset (CPI) to perform the ISAR imaging. If the CPI is too short (and the target barely rotates during the interval), motion compensation is easy, but the cross-range resolution will be poor. If the CPI is long and the target rotates significantly during it, high cross-range resolution can be obtained. However, motion compensation becomes more difficult due to scatterer MTRC and slow-time-dependent Doppler frequencies f p .
To choose a suitable CPI for image reconstruction, a socalled time window optimization procedure can be used. Its purpose is to locate a slow-time window during which the target rotational motion is as uniform as possible. In its simplest form, it chooses a single time window (length and location) from the entire recorded slow-time signal, which is then used in the RD imaging method. More generally, the optimization can be seen as a method for optimizing the parameters of the TFR that are used in the RID image reconstruction [15], [130].
Martorella et al. proposed choosing the slow-time window by maximizing the image contrast when using the RD method [197]. Another possibility is to estimate the instantaneous Doppler of different scatterers on the target, as demonstrated in [198]. These procedures are useful in cases where the target exhibits complex motions and the direction of the rotation changes, such as in imaging maritime targets [199]- [205]. More recently, the time window optimization procedure has been used to image ground and air targets as well [15], [130], [206].
Time window optimization is technically not a motion compensation procedure, since it does not produce estimates for the target motion. However, it has been shown to be a useful technique to avoid more complicated estimation strategies and thus to produce a computationally feasible approach for reconstructing images with sufficient crossrange resolution and focusing quality.

4) CROSS-RANGE SCALING
When using the above-mentioned RD, RID, or any other Doppler frequency estimation-based approach for image reconstruction, the target rotational motion parameters (θ as a function of slow time) are usually not explicitly estimated. To properly interpret the ISAR image, it is necessary to know the image projection plane and the spatial scale of the crossrange axis. The rotational motion parameters determine the relationship between the Doppler frequency and the crossrange position of the scatterers. The procedure for estimating the target rotation during the CPI-referred to as cross-range scaling-is a very difficult task with no general solution.
Several cross-range scaling methods based on estimating the slow-time phase history of multiple dominant point scatterers were proposed in [17], [207]- [218]. These methods utilize the relationship between the higher order slow time phase coefficients, the cross-range location, and the rotational velocity of the scatterers. The second-or third-order phase coefficients are obtained by isolating point scatterers and estimating their phase as a function of slow time. Because the phase estimation approaches are not usually directly applicable to RID image reconstruction, a suitable approach for incorporating them was recently considered in [219].
Another class of algorithms is based on optimizing a suitable objective function to estimate the target rotation. These optimization-based methods divide the CPI into two or more parts to produce sub-aperture images of the target.
These images can be considered scaled and rotated versions of each other. Thus, maximizing their correlation can be used to estimate the rotation between them. This provides an estimate for the target rotation between the sub-aperture images. For example, the rotation correlation and polar mapping methods [220], [221] are based on maximizing the correlation between two sub-aperture images as a function of the rotation applied to one of them.
More recently, approaches that attempt to extract features (e.g., dominant point scatterers or lines) from sub-aperture images and correlate them instead of the entire image have been proposed [222]- [228]. The sub-aperture technique has also been used to obtain 3-D locations for the target scatterers [229], and cross-range scaling approaches have recently been proposed for this application [230], [231]. Another possible method is to maximize the image contrast based on the rotational motion model used in the slow-timematched filtering [232]- [234].
Of the above-mentioned two classes, the phase estimationbased approaches require usually at least three range bins that contain a single isolated dominant scatterer with high SNR. When this assumption is not fulfilled, the phase estimation accuracy (and consequently also the cross-range scaling accuracy) significantly degrades. The optimization-based approaches are more generally applicable since they do not require the isolation of dominant scatterers. On the contrary, they work best when the target reflectivity distribution is more even, and highly dominant scatterers may even degrade their performance. As a drawback, the required optimization algorithms are often computationally intensive and do not guarantee that the global optimum is found. Moreover, the use of sub-aperture images (or features) works ideally only when the rotation is uniform, since non-uniform rotation causes a different scaling factor and defocusing between the subaperture images.
Thus, the choice of a suitable cross-range scaling method depends on the expected target type and motion dynamics. The strengths and weaknesses of the above-described two algorithm classes are summarized in Table 4.

C. DISCUSSION
Even though the single-step approach to translational motion compensation is generally more robust than range alignment followed by autofocus, its use has been limited mainly due to its higher computational cost. However, as the desired image resolution improves, it becomes necessary to abandon the more efficient two-step approach. Moreover, the single-step approach is often more straightforward to implement because VOLUME 9, 2021 each of the range alignment and autofocusing algorithms have their own tunable parameters and limitations.
The two-step approach (image reconstruction filter and cross-range scaling) in rotational motion compensation is fairly robust against different scenarios. It should be noted that the RID imaging approach offers some robustness even against residual translational motion errors. Cross-range scaling is performed independently and usually does not rely on the previous motion compensation steps to work properly. BP performs the image reconstruction and crossrange scaling in one step (since the rotation is explicitly used in the image reconstruction filter). However, currently there are no generally applicable computationally feasible approaches incorporating data-driven motion compensation in BP.
Reference [15] demonstrated the full ISAR processing chain consisting of range alignment, autofocus, RID imaging, and cross-range scaling using the state-of-the-art algorithms in each of the aforementioned steps. The results indicated that an image resolution of about 10-15 cm was possible to achieve in X-band, even in complicated target motion scenarios including highly non-linear translation and rotation. To increase the resolution beyond 10 cm, the approximations of the step-wise approach start to break down, and a more general method such as BP is required.

V. SPATIALLY DISTRIBUTED ISAR
Spatially distributed systems can be used to solve certain problems that are inherent in monostatic ISAR. For example, a distributed configuration is robust against unfortunate imaging geometries: a target that approaches a monostatic radar radially cannot be imaged in cross range because all the scatterers have the same Doppler frequency f p . Moreover, obtaining images from multiple aspect angles gives important additional information about the target reflectivity. A 3-D imaging capability is also obtained by using receivers at different heights. This information is of crucial importance in NCTR applications. However, multistatic systems are more complex and require additional considerations such as precise synchronization between the radar stations. These challenges, as well as additional comparisons between monostatic and multistatic ISAR, were analyzed in [235].
This section reviews the most important spatially distributed ISAR system concepts and the associated processing algorithms. Section V-A surveys the early references for biand multistatic systems. An interesting and important special case of bistatic systems is passive bistatic ISAR, which is reviewed in Section V-B. Section V-C investigates spatially distributed systems capable of locating target scatterers in 3-D using interferometric techniques.

A. BI-AND MULTISTATIC ISAR
Early theoretical developments, including a bistatic ISAR signal model and analysis concerning the resolution and geometrical distortions, were presented in [236]- [238]. One important finding was that under certain conditions, the bistatic configuration can be approximated as an equivalent monostatic setup, which greatly simplifies the ISAR processing [237]. However, as a drawback, this approximation causes geometric distortions and thus makes the interpretation of the bistatic ISAR images more difficult.
More recent developments in bistatic ISAR include dedicated motion compensation methods [239], a system with a ship-borne receiver and a shore-based transmitter with associated imaging and motion estimation algorithms [240], and techniques to mitigate the distortions caused by the monostatic approximation [241]- [244].
Early experimental results obtained with multistatic ISAR were reported in [245]. Theoretical developments of multistatic ISAR image reconstruction and resolution analysis along with experimental demonstrations followed in [246]. Further theoretical considerations related to multistatic motion compensation and image reconstruction were later published in [247], [248]. Brisken et al. developed the Fraunhofer experimental multistatic motion estimation system (FEMMES) and demonstrated the related imaging techniques in a series of publications [18], [249]- [253]. The FEMMES system is capable of forming multistatic images in X-band with approximately 10-cm resolution. This research demonstrated the entire multistatic processing chain, starting from synchronization using the direct path signal all the way up to multilateration-based motion estimation and BP-based image reconstruction.
The multistatic ISAR concept makes the motion compensation problem more tractable since the additional receivers provide new information for the data-driven estimation algorithms. However, multistatic systems are currently not very common due to the required system complexity. With developing radar technology and the increased interest in fully coherent multiple-input multiple-output techniques, multistatic ISAR has the potential to become ubiquitous in new ISAR applications.

B. PASSIVE BISTATIC ISAR
In recent years, passive radars have received a great deal of interest in the scientific community. Passive systems have several advantages over active systems, such as reduced cost, low vulnerability to electronic countermeasures, low probability of intercept, wide coverage areas, and continuous signal acquisition. These systems use different types of electromagnetic sources called illuminators of opportunity. Examples include active radar, communications, navigation, or digital broadcasting transmitters. As a drawback, the waveform is not under the operator's control, and thus the achievable range resolution may be inadequate for highresolution imaging.
Problems specific to passive ISAR are the often relatively poor range resolution (due to the limited bandwidth of conventional illuminators of opportunity) and the estimation of complex target motions. Due to the poor range resolution, small targets typically cannot be resolved in range-their response falls into a single range bin. Thus, the ISAR images obtained by passive ISAR systems are usually in the form of 1-D cross-range profiles. The low carrier frequency of typical illuminators requires a very long CPI to be processed in order to obtain good cross-range resolution. This poses strict requirements for motion estimation and compensation, which is a prerequisite for successful cross-range imaging.
The first demonstrations of passive radar imaging were in the context of airborne SAR. References [254]- [256] presented methods for obtaining images of moving targets in passive SAR images. In the early 2010s, the first passive ISAR applications using the Global System for Mobile and digital broadcasting signals for imaging air and maritime targets were published in [257]- [261]. As an early theoretical development, [262] proposed a passive ISAR imaging method for point targets relying on a test of a binary hypothesis, which was used to determine the target position and velocity.
Martorella and Giusti [263] presented the first rigorous theoretical development of passive bistatic ISAR imaging. A more recent pair of papers by Garry et al. described several additional theoretical and practical considerations, such as motion compensation and using different illuminators of opportunity [264], [265]. Colone et al. [266], [267] published the first comprehensive demonstration of the passive ISAR processing chain with WiFi-and very high frequency based experimental data. The processing chain for passive radar is in many ways more complicated than for active monostatic radars. This is because range compression and clutter cancellation require a more sophisticated approach. More recently, a number of passive ISAR experiments using digital television and digital video broadcasting and signals were reported in [190], [268]- [273]. Additionally, global navigation satellite systems signals have also been used for passive ISAR [274]- [276].

C. 3-D InISAR
Spatially distributed ISAR systems can be used as interferometric systems to estimate the 3-D positions of the target scatterers. In these systems, the setup is either bi-, multi-, or quasi-monostatic. The important aspect needed for the 3-D localization capability is having at least two receive (RX) antennas, which are separated in the elevation direction. This way, the phase differences between the RX channels can be used to estimate the height of the scatterers-using the same principle that is widely used in cross-track interferometric SAR systems, e.g., to produce digital elevation models of the ground surface. An important limitation of interferometry is that only a single scatterer can occupy the same range-crossrange resolution cell.
The aforementioned references mostly dealt with the quasi-monostatic setup, where closely spaced RX antennas are used to obtain 3-D estimates of scatterer positions. The images from the individual RX channels are combined to achieve better estimates of the target's dimensions. As previously mentioned, the quasi-monostatic setup has certain drawbacks (e.g., related to unfortunate imaging geometries). For this reason, 3-D InISAR has more recently been considered for spatially distributed bi-and multistatic setups [298]- [300]. These techniques are based on fusing the information contained in the 2-D images obtained from different RX stations of the multistatic setup to estimate the 3-D scatterer locations.
Besides their 3-D capability, these InISAR systems have an additional attractive ability: they can be used to estimate the direction of the target rotation vector. Monostatic systems can only resolve the magnitude of the rotation, not the direction. Because the direction of the rotation vector determines the image projection plane in 2-D imaging, its estimation is highly important for image interpretation and NCTR. References [301]- [306] described techniques for estimating the rotation vector from multisensor measurements. Table 5 summarizes the various aspects of using spatially distributed ISAR systems. The listed advantages make these techniques an important trend in ISAR that are receiving increasing research interest. The disadvantages related to system complexity are becoming less of an issue due to the developing technology.

VI. RECENT ISAR APPLICATIONS
Since the early days of ISAR, the most important applications have been in the military context, e.g., imaging ships, ground vehicles, and aircraft. Nowadays, the developing radar technology makes ISAR a more attractive concept in other applications as well. This section surveys a number of recent important application areas of ISAR: ground moving target VOLUME 9, 2021 imaging from an airborne platform, imaging targets with internal motion, and space surveillance.
The timeline in Figure 5 summarizes the most significant ISAR applications and algorithmic developments since the 1990s, including the methods described in Sections V-VII.
A. GROUND MOVING TARGET IMAGING ISAR can be used to image targets moving on Earth's surface using a space-or airborne SAR sensor. This concept was first demonstrated in Raney's paper [46], and it has attracted a great deal of interest more recently [13]. For this application, there is an additional challenge besides motion compensation: the strong ground clutter returns may mask the moving target reflections. For this purpose, multichannel systems that are able to cancel the ground clutter returns and detect the moving targets are commonly used.
An early method for performing motion compensation in this application was the prominent point processing method [2], [49], where the slow-time phase of multiple dominant scatterers was estimated to compensate the target motion. Another method proposed by Fienup [307] used a sharpness maximization procedure to focus patches of the SAR image that contained smeared moving target returns. Keystone formatting was initially proposed to compensate the moving target translation in SAR images [62]. Combined with an autofocus algorithm, this approach was demonstrated to produce focused moving target images in [126], [308]. Another more general motion compensation approach was proposed and demonstrated in [309]. All of these methods were designed for single-channel systems. Moreover, they assumed the moving target reflections to be much stronger than the ground clutter returns.
More recently, this problem has been approached by using the conventional ISAR processing scheme after isolating and projecting the moving target signal from the SAR image into the range-slow-time domain. Experimental demonstrations with Cosmo-SkyMed data using this approach were reported in [310]- [312]. Furthermore, a frequency-domain focusing algorithm and experimental demonstration with airborne data were reported in [313], [314].
When imaging weakly reflecting targets, the cancellation of the ground clutter signal is a prerequisite for successful ISAR image reconstruction. A common method for this purpose is space-time adaptive processing (STAP) [315], which can be used for systems with multiple RX channels. An approach combining STAP and ISAR processing was developed and demonstrated in [316]- [321]. Additionally, a virtual multichannel approach, i, e., using a single-channel system to emulate a multichannel system, was used and compared to actual multichannel systems in [322].
Ground moving target detection with airborne radar using STAP is a well-established technique. Adding the capability of ISAR imaging to this technology offers the possibility to obtain a significantly improved situational awareness. The targets can not only be detected and located but possibly also classified or recognized based on their ISAR signatures.
This application presents a novel algorithmic challenge because the SNR (or more precisely, the signal-to-clutter ratio) is often much lower than in conventional ISAR processing. Even though the STAP-ISAR approach shows promise, much work remains to be done to make this concept more robust and generally applicable.

B. IMAGING TARGETS WITH ROTATING PARTS
In conventional ISAR processing, the target is assumed to be a rigid body. Under this assumption, the position of each target scatterer depends only on the rotation θ(t) of the rigid target. However, certain common ISAR targets, such as helicopters, aircraft, and drones, contain rapidly rotating parts, violating the rigid body assumption. For example, the extremely fast rotation of rotor blades causes highly non-stationary Doppler frequency modulation, which makes it very difficult to obtain a focused ISAR image of the target. Even if the translation and rotation of the rigid part of the target are properly estimated and compensated, the micro-Doppler caused by the non-rigid rotating part may blur and defocus the entire image. Moreover, the echo of the blurred rotating part may mask important rigid body target features.
To overcome the aforementioned problem, several methods for separating the rigid body and rotating part reflections have been considered in the ISAR literature. On the timefrequency plane, the rigid body echoes follow straight lines whereas the rotating parts follow a sinusoidal curve. Fundamentally, all of the separation methods are in some way based on using this qualitative difference in the timefrequency behavior of the echoes.
An early development based on using a TFR and order statistics to obtain the rotating part signals was proposed in [323]. This approach was further developed, analyzed, and demonstrated in later publications [324]- [326]. An approach based on this concept but using a procedure to improve the time-frequency resolution was proposed in [327].
Other separation methods have since been proposed, e.g., based on the Hough transform [328], [329] and the Radon transform [330]. Various approaches based on different decomposition methods from linear algebra-principal component analysis [331], empirical mode decomposition [332], complex local mean decomposition [333], and bivariate variational mode decomposition [334]-have also been used for this purpose.
The methods for imaging non-rigid targets with rotating parts are an important extension of ISAR imaging. Besides enabling successful imaging of new targets in old applications, these techniques may open up new potential application areas, such as imaging moving persons.

C. SPACE SURVEILLANCE
As the number of operating and defunct satellites orbiting the Earth is continuously increasing, the risk of collisions between active satellites and space debris becomes more pronounced. Thus, it is increasingly important to have an accurate situational awareness of near-Earth space. A number of radar systems have recently been developed to improve space surveillance capabilities. Some of these systems are low resolution and only provide the ability to detect and track targets. Others are capable of high-resolution ISAR imaging, providing an important information source for monitoring and characterizing the observed satellites.
The earliest space surveillance systems LRIR and ALCOR were developed in the 1970s, as previously described in Section III. Since then, the Lincoln Laboratory has developed the Ku-band Haystack Auxiliary system with 13-cm range resolution and the Ka-band Millimeter-Wave radar with 7-cm resolution [335]. Most recently, the W-band Haystack Ultrawideband Satellite Imaging Radar was put into operational use [335]. It can produce ISAR images with an extremely high resolution of approximately 2-cm, representing the absolute state of the art in long-range noncooperative ISAR imaging performance.
Another high-performance space surveillance radar system is the Tracking and Imaging Radar (TIRA) developed by Fraunhofer FHR [6], [336]. This Ku-band system is capable of producing decimeter resolution ISAR images of satellites. The reconstructed images from TIRA have recently been used to analyze the rotation of ENVISAT [7] and to monitor the uncontrolled re-entry of the Tiangong-1 satellite [337]. The most recent new space surveillance radar with high-resolution ISAR imaging capability is the X-band Imaging of Satellites in Space system [8] developed by the German Aerospace Center. While still in its commissioning phase, the system has already produced high-resolution ISAR images of the International Space Station [8], [338].
As space surveillance is becoming an increasingly important task, ISAR techniques will be pushed forward by the requirements to detect and identify very small objects (various space debris particles) with extremely high resolution. The fast motion and uncontrolled rotation of space debris particles pose a novel challenge for ISAR motion compensation algorithms.

VII. RECENT ALGORITHMIC DEVELOPMENTS
During recent years, there has been great research interest in CS and ML techniques. In addition to theoretical developments in these fields, they have been used in an increasing number of engineering applications. In this section, we survey their use in ISAR imaging.

A. CS IN ISAR
The ISAR imaging problem can be formulated as a linear inverse problem where the unknown target reflectivity distribution (RCS as a function of spatial position) is reconstructed from noisy measurements of scattered electromagnetic signals. CS techniques exploit the sparsity of the target signal to find sparse solutions to the under-determined linear equations that constitute the imaging problem. The simplest form of sparsity is to consider the target as a small number of dominant point scatterers-often an accurate assumption in ISAR.
CS algorithms can reconstruct the target image with suppressed sidelobes and increased spatial resolution of point scatterers (super resolution). Furthermore, CS provides increased robustness to limitations in the data quality and quantity. Typical data limitations include missing pulses and sparse frequency spectrum waveforms. Thus, CS can significantly increase the ISAR image quality, providing enhanced target information for a subsequent NCTR or ATR stage.
The use of CS techniques in SAR and other radar imaging applications was reviewed in [339], [340]. These papers provide a concise theoretical description of the image reconstruction procedure. Moreover, they discuss various practical considerations and CS-based radar imaging applications. However, both of these papers are mainly focused on SAR-only a few ISAR references are given. To fill in this gap, we proceed to survey the CS-based methods used in the context of ISAR imaging.
The CS-based approach is useful in situations where the available radar data are incomplete or non-uniformly sampled-e.g., pulses are missing, the target rotates rapidly and non-uniformly, and frequencies are missing due to jamming. References [341]- [349] demonstrate ISAR image reconstruction using CS when a number of pulses (range profiles) are missing, i.e., in the case of sparse slowtime sampling. References [350]- [352] proposed CS-based methods to deal with the sparse angular sampling caused by a rapid target rotation. Another possibility is to use sparse sampling in the waveform by using a sparse steppedfrequency waveform [353]- [355]. This approach can be beneficial in applications where covert radar operation with a low probability of intercept is required.
CS-based image reconstruction has been found to be useful in 3-D InISAR processing [356]. Moreover, it is useful in the case of multistatic 3-D ISAR imaging because it offers a way to avoid the coherent combination of the signals from the different RX stations [253], [357]. Examples of CSbased 3-D InISAR imaging methods can be found in [358], [359]. Furthermore, specific methods and results have been presented for space [352], [360], maritime [356], and general maneuvering target [361], [362] imaging.
An important feature of CS is its ability to enhance the image resolution. Approaches in which the achievable image resolution no longer only depends on the signal bandwidth and target rotation are commonly referred to as super-resolution techniques. Reference [363] presented comparisons between CS and other super-resolution techniques in ISAR. A number of different CS methods have been used for this purpose: sparse Bayesian learning [364]- [371], the RELAX method [372], matching pursuit [373], and dictionary learning [374], [375].
In ISAR processing, the unknown target motion usually hinders the direct application of CS-based reconstruction techniques. Moreover, the motion compensation approaches described in the previous sections are usually not directly applicable: they have been designed for linear reconstruction methods and Nyquist sampled data. Thus, several approaches for performing data-driven motion compensation for CSbased ISAR have been developed. For example, the image sharpness maximization can be included in the CS optimization problem [376]. It is also possible to include the unknown phase errors in addition to the sparse scattering coefficients in the CS optimization problem formulation. This approach, demonstrated in references [377]- [383], is equivalent to joint image reconstruction and autofocusing.
Despite the potential and voluminous literature on the subject, CS-based ISAR imaging in operational systems is still very limited. Most of the reported results dealt with idealized simulation scenarios or simple, closely controlled target scenarios. The reviewed methods require stringent approximations to hold in order to produce meaningful ISAR images. Thus, their performance heavily depends on a favorable target scenario. Moreover, CS-based approaches generally require a high SNR to provide useful results, and their computational complexity is often much higher than for the linear image reconstruction and motion compensation methods described in the previous sections. Table 6 presents a summary of the aforementioned advantages and disadvantages of the CS-based ISAR imaging approaches.

B. ML IN ISAR
The recent developments in computing power and ML techniques are making their way into ISAR. Previously, the use of ML in the context of ISAR was limited to NCTR and ATR, i.e., classifying targets based on the reconstructed images. This subject is outside the scope of this survey (for a relatively recent review of ATR techniques, see [384]).
During the last few years, a number of methods using ML techniques in ISAR image reconstruction have been proposed. In [385]- [387], the authors used a convolutional neural network (CNN) to reconstruct ISAR images from sparsely sampled data and demonstrated their improved performance compared to CS-based ISAR. Furthermore, an autofocus approach was integrated into the CNN framework in [388]. Neural networks have also recently been used to enhance image resolution [389], [390]. Additionally, the CNN-based approach has recently been demonstrated in the ground moving target imaging application [391].
The inherent challenge in ML-based imaging approaches is the need for a large amount of representative training data. In the radar application, it is notoriously difficult to obtain comprehensive training data: the RCS modeling and simulation for different aspect angles and frequencies is cumbersome and requires accurate target models, which are not readily available for all possible targets. Moreover, sharing and obtaining data is quite restricted for military targets.
As an advantage, proper ML architecture avoids all the conventional challenges with ISAR processing: by incorporating the unknown motion dynamics in the training data, the neural network performs all the required processing without the need to explicitly design separate motion estimation algorithms. In conclusion, despite the promising concept, ML-based ISAR imaging is still in its infancy.

VIII. CONCLUSION
This paper presented a comprehensive literature survey of ISAR imaging techniques and applications from the late 1950s to the present day. The timelines in Figures 4 and 5 summarize the important technical milestones described in this paper. After introducing the basic ISAR principles and terminology, we described the early historical developments, starting from the imaging of the moon and planets and continuing to early algorithmic developments and imaging airborne targets. Then, we proceeded to analyze the various methods used for motion compensation and image reconstruction in monostatic ISAR processing. We continued by analyzing spatially distributed ISAR concepts, which offer various advantages and hold great potential for future ISAR developments. We proceeded by describing important current applications, including ground moving target imaging, nonrigid target imaging, and space surveillance. Finally, we surveyed the ISAR literature using the latest developments from CS and ML.
In general, non-cooperative ISAR imaging is a very challenging task with no general purpose solution-hence the voluminous literature on the subject. However, the increasing computational resources has allowed for the development of more generally applicable algorithms and the use of CS and ML techniques. To provide images with increased information content, the trend of ISAR is moving toward spatially diverse systems and higher spatial resolution. The rapid development of affordable high-resolution radar systems continually opens up new possible application areas for ISAR imaging, transforming this from a niche field to the mainstream.