ISAR Processing

Authored by: M. Martorella , F. Berizzi , E. Giusti , A. Bacci

Radar Imaging for Maritime Observation

Print publication date:  June  2016
Online publication date:  June  2016

Print ISBN: 9781466580817
eBook ISBN: 9781315374253
Adobe ISBN:

10.1201/9781315374253-5

 

Abstract

In this section the ISAR signal model for a monostatic configuration will be presented.

 Add to shortlist  Cite

ISAR Processing

3.1  ISAR Signal Model

In this section the ISAR signal model for a monostatic configuration will be presented.

ISAR geometry

Figure 3.1   ISAR geometry

Consider the typical ISAR geometry shown in Figure 3.1 in which the transmitter and the receiver are colocated in the origin of the Cartesian reference system Tξ. As stated in Chapter 1, the signal in the frequency/slow time domain at the output of the matched filter can be expressed as

3.1 S ( f , n ) = C rect ( n N ) | S T ( f , n ) | 2 V f ( x ) e j 2 π f τ ( x , n ) d x
where in conventional active radar | S T ( f ) | 2 rect ( f f 0 B ) and rect ( n N ) = u ( n + N 2 ) u ( n N 2 1 ) defines the signal support in the frequency/slow time domain. The function u(·) yields 1 when (·) > 0. τ (x, n) is the delay time of a scatterer placed in x = (x1, x2, x3)T with respect to the reference system Tx and is expressed as follows:
3.2 τ ( x, n ) = 2 R ( x, n ) c
where R(x, n) is the distance between the radar and the scatterer placed in x in the nth sweep, c is the speed of light in a vacuum and ΩT(n) and Ω(n) are, respectively, the total target rotation vector and the effective target rotation vector, which is obtained as the projection of ΩT(n) onto the plane orthogonal to the radar line of sight (LoS).

By assuming that the radar target distance is much greater than the size of the target (straight iso-range approximation), the term R(x, n) can be approximated as follows:

Rotation along the LoS

Figure 3.2   Rotation along the LoS

3.3 R ( x, n ) R 0 ( n ) + x T i L o S ( n )
where R0(n) is the distance between the radar and a reference point on the target at the time instant nTR, x is the column vector that identifies a scatterer on the target and iLoS is the column unit vector that identifies the radar LoS.

It is worth pointing out that Equation (3.3) holds true only when considering the target as a rigid body. In this case, its motion with respect to the radar can be considered as the superimposition of the translational motion of a reference point and a rotation vector applied to this reference point. The rotation vector takes into account both the rotation due to the translation and the proper target rotation. It is worth pointing out that, since the radar can measure only the distance variation along the LoS, any rotation around the y2 axis does not produce any effect. Then, the only component to be considered is Ω(n), as it contributes to the synthetic aperture formation. This result is evident by observing Figure 3.2 in which the case in which the target undergoes a rotation around the LoS is represented. In such a case, the distance between the radar and each scatterer on the target does not change, R(x, t1) = R(x, t2). Then, such a rotation does not produce any effect from the radar point of view.

Let us assume that the effective rotation vector is constant during the overall observation time, i.e., Ω(n) ≈ Ω for |nTR| < Tobs, and consider the geometry depicted in Figure 3.1 in which the axis x3 is oriented along Ω. Under this assumption, the inner product y2(x, n) = x · iLoS(n) can be expressed in a closed form by solving the differential equation system expressed as follows:

3.4 y ˙ ( n ) = Ω T × y ( n )
with the initial condition y(0) = x.

The resulting closed form solution is

3.5 y ( n ) = a + b cos ( Ω T n T R ) + c Ω T sin ( Ω T n T R )
where
a = Ω T x Ω T 2 Ω T b = x Ω T x Ω T 2 Ω T c = Ω T × x and  Ω T = | Ω T | Ω T = | Ω T |
with the above-mentioned choice of the reference system, ΩT = (0, Ω T 2 , Ω) (with respect to Ty) where the component Ω T 2 does not produce any effect as stated above.

Equation (3.5) can be reasonably approximated by its first order Taylor series around t = 0 as follows:

3.6 y ( n ) a + b + c ( n T R ) = x + c ( n T R )
resulting in:
3.7 y 2 ( x , n ) = x 2 + Ω n T R x 1
where Ω = |Ω|.

By substituting Equation (3.6) in Equation (3.1) and by defining W ( f , n ) = C rect ( f f 0 B )

rect (n/N) (as defined in Equation (1.17)), the received signal can be rewritten as
3.8 S ( f , n ) = W ( f , n ) V f ( x ) e j 4 π f c ( R 0 ( n ) + x 2 + Ω x 1 n T R ) d x

The same result can be obtained in a different way considering that the LoS unit vector in the nth sweep in the Tx reference system is given by

3.9 i L o S ( n ) = [ sin ( Ω n T R ) cos ( Ω n T R ) 0 ]

As a consequence, the inner product in Equation (3.3) results in

3.10 y 2 ( x T , n ) = x i L o S ( n ) = x 2 cos ( Ω n T R ) + x 1 sin ( Ω n T R )

ISAR processing chain

Figure 3.3   ISAR processing chain

Equation (3.1) can therefore be written as

3.11 S ( f , n ) = W ( f , n ) x 1 , x 2 f ( x 1 , x 2 ) e j 4 π f c [ R 0 ( n ) + x 2 cos ( Ω n T R ) + x 1 sin ( Ω n T R ) ] d x 1 d x 2
where f ( x 1 , x 2 ) = x 3 f ( x ) d x 3 is the projection of the target's reflectivity function onto the image projection plane (IPP).

By assuming the total aspect angle variation is small during the observation time, i.e., Δϑ = ΩTobs ≪ 1, Equation (3.11) becomes Equation (3.8).

3.2  ISAR Image Formation Chain

The steps to be performed to obtain an ISAR image are summarized in the functional block diagram shown in Figure 3.3. A time interval of the overall observation within which the effective rotation vector can be assumed constant time is selected first by means of the time windowing selection. This operation is important because when the radar observes the target for long time intervals, this assumption does not hold true and the model in Equation (3.8) cannot be applied. Through the time windowing process, a subset of data for which the effective rotation vector is constant is properly selected to form the ISAR image. After that, the motion compensation operation is performed. This operation aims at removing the phase term R0(n) that accounts for the target radial motion. After motion compensation, the image formation step is performed via a range-Doppler (RD) approach. This operation provides the ISAR image in the delay time/Doppler coordinates (τ, ν). Finally, in order to better asses the target size, a scaling operation to spatial coordinates is needed.

All these processing steps will be extensively described in the following sections.

3.2.1  Image Formation

The ISAR image is formed via the RD approach. Consider the received signal expressed in Equation (3.11). Assuming that the phase term due to the translational motion R0(n) is perfectly compensated during the ISAR image formation processing, the motion compensated signal is expressed as follows:

Spatial frequency domain

Figure 3.4   Spatial frequency domain

3.12 S C ( X 1 , X 2 ) = W ( X 1 , X 2 ) x 1 , x 2 f ( x 1 , x 2 ) e j 2 π ( x 1 X 1 + x 2 X 2 ) d x 1 d x 2
where (X1, X2) are the spatial frequencies defined as
3.13 { X 1 ( f , n ) = 2 f c sin ( Ω n T R ) X 2 ( f , n ) = 2 f c sin ( Ω n T R )

From Equation (3.13) it is obvious that the motion compensated signal corresponds to the windowed Fourier transform (FT) of the projected target's reflectivity function. The term W(X1, X2) defines the region (a sector of annulus) in the spatial frequencies domain in which the signal is defined, see Figure 3.4.

As stated above, by assuming that the total aspect angle variation is small, i.e., Δϑ = ΩTobs ≪ 1, the spatial frequencies can be approximated as follows:

3.14 { X 1 ( n ) = 2 f 0 c ( Ω n T R ) X 2 ( f ) = 2 f c

It is worth noting that for X1 the frequency f has been substituted by the central frequency f0, being the result of the approximation of the polar domain in Figure 3.4 with a rectangular window that intersects the angular sector at the coordinates X 2 = 2 f 0 c

. It should be noted that this approximation is the one that leads to the minimum error as it can be inferred by examining Figure 3.5 where both the polar and the rectangular domains are represented.

Spatial frequency domain with small aspect angle variation

Figure 3.5   Spatial frequency domain with small aspect angle variation

By assuming that the target is composed by K point-like scatterers and neglecting all the interactions among the scatterers, the target reflectivity function can be expressed as follows:

3.15 f ( x 1 , x 2 ) = k = 1 K σ k δ ( x 1 x 1 k , x 2 x 2 k )
where x 1 k , x 2 k and σk are the cross-range coordinate, the range coordinate and the reflectivity value of the k(th) scatterer, respectively.

By substituting Equation (3.15) in Equation (3.12) the received signal after motion compensation results

3.16 S C ( X 1 , X 2 ) = W ( X 1 , X 2 ) k = 1 K σ k e j 2 π ( X 1 x 1 k + X 2 x 2 k )

The ISAR image can be obtained by applying a 2D-FT to Equation (3.16) and can be expressed as

3.17 I ( x 1 , x 2 ) = w ( x 1 , x 2 ) f ( x 1 , x 2 ) = k = 1 K σ k w ( x 1 x 1 k , x 2 x 2 k )

By substituting Equation (3.14) in Equation (3.16), the received signal after motion compensation becomes:

3.18 S C ( f , n ) = W ( f , n ) k = 1 K σ k e j 2 π ( f τ k + n T R v k )
where the delay-time and the Doppler coordinates are defined as:
3.19 τ = 2 x 2 c v = 2 f 0 Ω x 1 c

The ISAR image in the delay time/Doppler domain, I(τ, ν) can then be obtained via a two-dimensional Fourier transform of Equation (3.18) as follows:

3.20 I ( τ , v ) = w ( τ , v ) k = 1 K σ k δ ( τ τ k , v v k ) = k = 1 K σ k w ( τ τ k , v v k )
where
3.21 w ( τ , v ) = 2 D-FT [ W ( f , t ) ] = B T o b s sin c ( τ B ) sin c ( v T o b s ) e j 2 π f 0 τ
and 2D-FT[·] is the two-dimensional Fourier transform operator, sin c ( x ) = sin ( π x ) π x and ν [ P R F 2 , P R F 2 ] . w ( τ , ν ) defines the system point spread function (PSF) and can be used to find the resolution in the delay time/Doppler dimensions, i.e., δτ and δν, defined as the first null of the two sinc functions in both directions
3.22 δ τ = 1 B δ v = 1 T o b s

It is worth pointing out that Equation (3.20) represents the ISAR image in the delay time/Doppler domain. To obtain the target ISAR image in the spatial coordinates (range/cross-range) a scaling operation is needed. This scaling operation can be performed by exploiting Equation (3.19) as follows:

3.23 { x 1 = c 2 f 0 Ω v x 2 = c 2 τ
and, as a consequence, the spatial resolutions can be expressed as:
3.24 { δ x 1 = c 2 f 0 Ω T o b s δ x 2 = c 2 B

Equation (3.23) shows that the scaling operation from the delay time to the range coordinate is straightforward while this is not true for the scaling operation from the Doppler frequency to the cross-range coordinates. In fact, this coordinates conversion involves the knowledge of the modulus of the effective rotation vector, Ω, which depends on the motion of the non-cooperative moving target and, as a consequence, cannot be assumed known in the ISAR scenario.

As can be noted from Equation (3.24), the resolution becomes finer as the aspect angle variation increases. This is in contrast with the small aspect angle variation assumption which is mandatory for the applicability of the range Doppler image formation algorithm. This hypothesis determines an upper bound in the cross-range resolution that can be achieved with this algorithm.

Another important consideration concerns the definition of the IPP where the target is represented. As stated above, the reference system Tx has been chosen so that the x2 axis is oriented along the LoS for t = 0 and the x3 axis is parallel to the vector Ω that produces the aspect angle variation. Since any rotation around the LoS does not produce any variation in the aspect angle, it is quite obvious that Ω is the projection onto the plane orthogonal to the LoS of the total rotation vector ΩT and can be expressed as follows:

3.25 Ω ( n ) = i L o S × [ Ω T ( n ) × i L o S ]

Therefore, the cross-range axis x1 is defined as:

3.26 x 1 = Ω × i L o S
where the dependency on n is dropped out because the IPP can be defined only if Ω is constant during the observation time. It is evident that, as it happens for the cross-range scaling operation, the IPP depends on the unknown quantity so that the IPP cannot be a priori predicted leading to some difficulties in the ISAR image interpretation.

It is worth pointing out that the above-mentioned RD technique is based on the assumption that the Doppler frequency of each scatterer, relative to the focusing point, is constant during the CPI (coherent processing interval). This hypothesis usually holds true for low spatial resolution (of the order of meter) and in the case in which the target does not undergo fast maneuvers and/or is affected by significant oscillation such as roll, pitch and yaw. In the case in which very high resolutions are required, longer CPIs are needed and the Doppler frequency becomes time-varying. This can happen also when the target undergoes angular motions, as in the case of ships. In this case the first order Taylor approximation in Equation (3.6) and (3.7) is not valid. To solve this problem, the time frequency transforms (TFTs) is used instead of the RD approach to form ISAR images. TFTs are, in fact, suitable for the analysis of non-stationary signals. An example of a TFT-based technique is the range instantaneous Doppler (RID) approach described in [2].

3.2.2  Time Window Selection

As stated before, the applicability of the RD technique relies on the assumption that the effective rotation vector is constant during the observation time, i.e., Ω(n) Ω. However, the target's own motion may induce a non-uniform target rotation vector. In particular, maritime targets can undergo very complex motions (roll, pitch and yaw) induced by the sea surface. In order to minimize the target's rotation variations, the CPI can be controlled via a time-windowing approach.

Typically, a fixed window length can be set and a sliding window processing along the whole observation time is applied to obtain a sequence of ISAR images of the targets. Among the obtained frames, the more suitable for target classification and recognition is selected.

It is worth pointing out that, in order to obtain good classification and recognition performances, the most important requirement is a good level of focus, as the target details would then appear sharper than in defocused images. As stated above, the image cross-range resolution is inversely proportional to the total aspect angle variation and, as a consequence, to the considered CPI so that longer CPI involves a finer resolution in the cross-range dimension. On the other hand, longer CPIs increase the chance that the target rotation vector may not be considered constant and therefore producing a defocused image. It is quite obvious that a trade-off must be identified to obtain a well focused and a high resolved image at the same time.

An automatic time windowing technique is proposed in [9] and will be here recalled for the sake of clarity. Specifically, the position of the window along the whole observation time and its length are automatically selected in order to obtain the image with the highest focus. The concept can be explained referring to Figure 3.6, in which the data is refereed to be distributed along the time axis n (n denotes the sweep index at time t = nTR). The considered temporal window is defined by two parameters, namely nw that identifies the window's position and Nw that denotes the window's length in number of sweeps, so that the temporal length of the windows is Tw = NwTR.

The measure of the image focus is made through the image contrast (IC) which is assumed to be maximum when the window position (nw) and the window length (Nw) identify a subset of data where the conditions of constant rotation vector and resolution are optimal in terms of image focus.

In particular, the optimal time window is obtained by maximizing the IC with respect to the pair (nw, Nw). This problem can be mathematically formulated as follows:

Time windowing concept

Figure 3.6   Time windowing concept

Optimal time window estimator

Figure 3.7   Optimal time window estimator

3.27 ( n w , o p t , N w , o p t ) = arg max ( n w , N w ) ( I C ( n w , N w ) )

It should be noted that the variables (nw, Nw) are both discrete so that the problem in Equation (3.27) is a discrete optimization problem. Specifically, such a problem can be classified as a non-linear Knapsack Problem [11]. The solution proposed in [9] is based on a double linear search, which can be briefly described by the following steps:

  1. Contrast maximization with respect to nw for a given guess Nw,in. Let nw,opt be the solution of such maximization;
  2. Optimization with respect to Nw with nw = nw,opt.

The procedure is depicted in Figure 3.7 for the sake of clarity. The justification of this procedure is heuristic. It can be observed both in simulated and real ISAR data of several types of targets that the position of the optimal time-window is almost independent on the window length and, as a consequence, the IC peak position along the slow-time axis, n, does not change by changing Nw. This can be physically explained by considering that the target own motion is characterized by regular motions at given times and less regular motion at other times. For example, a ship usually undergoes regular pitch and roll motions due to the sea surface state. This regular motion may be disturbed by an incoming wave that generates complex motions which cause the rotation vector to change rapidly. For this reason, the estimation of nw,opt and Nw,opt can be performed separately.

3.2.3  Motion Compensation

As stated in Section 3.2.1, the radial motion compensation is a fundamental step in forming well-focused ISAR images. This step is typically termed image autofocus and consists of removing the phase term 4 π f c R 0 ( n )

in Equation (3.8). This term depends on the radial motion of the reference point. When no external data are available, this operation must be performed only by using the radar received signal and, for this specific reason, the image focusing process is called ISAR autofocus.

Several techniques have been developed by the ISAR research community, each of them showing pros and cons. Such techniques can be classified as parametric and non-parametric methods [1]. Parametric methods need a parametric model of the radar received signal while non-parametric techniques do not make any assumption. Some autofocus techniques are detailed in the next subsections.

3.2.3.1  Image Contrast Based Autofocus

The image contrast based autofocus (ICBA) is a parametric technique that is based on the obvious concept that a higher value of IC corresponds to a more focused image [10]. This autofocus technique aims at removing the term R0(n) via an iterative estimation of the motion parameters. For a relatively small observation time, Tobs = NTR, and relatively smooth target motions, the radar-target residual distance can be expressed by exploiting a Lth polynomial model as follows:

3.28 R 0 ( n ) = l = 0 L 1 l ! α l ( n T R ) l
in which αl are the above-mentioned focusing parameters which can be grouped in a vectorial form as α = [α0, α1, …, αL−1]T.

The estimation of R0(n) results in the estimation of the target motion parameters via the maximization of the IC with respect to α as:

3.29 α ^ = arg max α [ I C ( α ) ]
where the IC is defined as follows:
3.30 I C ( α ) = A v [ [ I ( τ , v , α ) A v [ I ( τ , v , α ) ] ] 2 ] A v [ I ( τ , v , α ) ]
and where Aυ[·] denotes the average operation over the variables τ and ν, I(τ, ν, α) is the ISAR image magnitude or intensity (power) after motion compensating the target's translational motion by using α as focusing parameters and is expressed as follows:
3.31 I ( τ , v , α ) = | 2 D-FT [ S ( f , n ) e j 4 π f c R 0 ( α , n ) ] | p

Referring to Equation (3.31), p = 1 in the case of image amplitude and p = 2 in the case of image intensity. R0(α, n) is the radial motion model evaluated with the α motion parameters.

3.2.3.2  Image Entropy Based Autofocus

A similar approach to the one proposed in Section 3.2.3.1 can be devised by substituting the IC with the image entropy (IE) expressed as follows:

3.32 I E ( α ) = ln ( I ¯ ( τ , v , α ) ) I ¯ ( τ , v , α ) d τ d v
where I ¯ ( τ , ν , α ) = I ( τ , ν , α ) A υ [ I ( τ , ν , α ) ] .

As for the IC, the IE is a good indicator of the image focus. Unlike the IC, the smaller IE value the better the image focus.

3.2.3.3  Dominant Scatterer Autofocus

The dominant scatterer autofocus (DSA), also known as the hot spot (HS), is a non-parametric two stages technique. The first stage consists of a rough alignment before applying the phase compensation as a second step. This technique is derived from the studies in two different areas of research, namely the time delay estimation [3] and the adaptive beamforming [12]. A complete description of this algorithm can be found in [4] while a brief overview is reported here.

After measuring and storing the complex envelopes of the echo samples, high resolution range profiles in the nth sweep, S(τ, n), can be obtained via Fourier transforming Equation (3.1) along the frequency dimension. A cross-correlation shift operation can be performed between successive range profiles in order to obtain a rough range profile alignment. Let Sa(τ, n) be the roughly aligned profiles expressed as follows:

3.33 S α ( τ , n ) = A ( τ , n ) e j ϕ ( τ , n )
where τ denotes a range cell (delay time) and n is the pulse number and A(τ, n) ∝ |Sa(τ, n)|ejϕA, where ejϕA(τ) is a phase term independent on the slow-time variable, n.

Then, a search along the range coordinate is performed in order to find a dominant and stable scatterer. The range cell where such scatterer is located is the reference range denoted by τ0 and the corresponding range profile is given by

3.34 S α ( τ 0 , n ) = S 0 ( n ) = A ( τ 0 , n ) e j ϕ 0 ( n )

The τ0 value is found by measuring the normalized echo variance in each range cell and is determined to be the range for which the variance value is minimum. This approach relies on the assumption that a dominant scatterer with a large radar cross-section exists and, therefore, the measured phase can be attributed to the phase generated by one point scatterer. The next step performs a phase compensation by using the phase history of S0(n). For the range cell data corresponding to τ0 the result is expressed as follows:

3.35 S 0 , c ( n ) = A ( τ 0 , n ) A

Applying the same operation to the other range cells, motion compensation is achieved and the result is expressed as

3.36 S C ( τ , n ) = A ( τ , n ) e j [ ϕ ( τ , n ) ϕ 0 ( n ) ]

This algorithm is also known as minimum variance algorithm (MVA) due to the criterion used to choose the reference range cell.

A more robust version of this algorithm, which combines the echoes from several range cells has been proposed in [4] and it is named multiple scatterer algorithm (MSA). It is based on an average operation performed on the phase differences of Mc reference scatterer (after phase unwrapping) in order to provide the phase correction. Typically, three reference range bins (Mc = 3) are sufficient to produce a focused image. Mathematical details are given below. Let the mth reference cell be expressed as:

3.37 S m ( n ) = A ( τ m , n ) e j ϕ m ( n )
and let Mc be the number of selected range cells. The estimation of the phase history to be compensated, ϕ0(n), is
3.38 ϕ 0 ( n ) = 1 M c m = 1 M c ϕ m ( n )

Phase compensation is then performed as expressed in Equation (3.35) and (3.36). The DSA algorithm is summarized in the functional block depicted in Figure 3.8.

3.2.3.4  Phase Gradient Autofocus

The Phase Gradient Autofocus (PGA) can be seen as an extension of the DSA algorithm proposed in [4]. It also exploits the information that remains in those range cells discarded by the minimum variance (MV) range selection. The PGA algorithm replaces the phase compensation approach used in DSA with a solution based on the maximum likelihood (ML) estimator. The ML estimator theoretically uses the information contained in all the range cells. However, only the range cells in which the SNR is high enough are considered to improve the PGA performance.

DSA algorithm functional block diagram

Figure 3.8   DSA algorithm functional block diagram

PGA functional block diagram

Figure 3.9   PGA functional block diagram

The functional block diagram that summarizes the PGA algorithm is depicted in Figure 3.9.

A range Doppler image is formed and a rough range alignment is performed. Then, the peak value in each range cell, which is supposed to be the return from the dominant scatterer, is found along the Doppler dimension and centershifted and filtered in the Doppler domain (low-pass filter). Each range cell is then transformed back in the pulse domain via an inverse Fourier transform (IFT) to obtain the phase shifted and filtered time histories.

This operation corresponds to isolating the strongest scatterer contribution for each range cell and forcing it to zero-Doppler, which can be interpreted as a way to remove the scatterer radial motion.

Then the phase gradient is estimated by using the ML approach as described as follows.

Let the arbitrary mth range cell time history be represented as follows, where two consecutive samples are considered:

3.39 S α ( τ m , n 1 ) = A ( τ m ) e j ϕ ( n 1 ) + S n ( τ m , n 1 ) S α ( τ m , n ) = A ( τ m ) e j ϕ ( n ) + S n ( τ m , n )
where τm denotes the mth range cell, A(τm) = |Sa (τm, n − 1)| ≃ |Sa(τm, n)| denotes the amplitude, which is assumed to be constant for two consecutive samples, ϕ(n) is the phase at time n and Sn(τ, n) is an additive white Gaussian noise sample.

The phase gradient between two consecutive samples can be defined as

3.40 Δ ϕ = ϕ ( n ) = ϕ ( n 1 )

The derivation of the ML estimator of the phase gradient is given in [5]. The solution is expressed as follows

3.41 Δ ϕ ^ ( n ) = ( m = 1 M c S α ( τ m , n ) S α * ( τ m , n 1 ) )
where Mc is the number of range cells used in the estimation process and where the symbol ∠(·) indicates the phase of a complex number.

The phase correction term is then obtained by integrating the estimated phase gradient, as follows:

3.42 ϕ ^ ( n ) = l = 1 n Δ ϕ ^ ( l ) ϕ ^ ( 0 ) = 0

Finally, the phase correction term is used to compensate the phase error and obtain the compensated signal.

3.2.4  Image Scaling

The ISAR processing generates a two-dimensional high-resolution image of a target in the delay-time/Doppler domain. In order to determine the size of the target it is, however, required to have a fully scaled image in the range/cross-range domain. The scaling along the range dimension can be performed by using the well-known relationship x 2 = c τ 2

, where x2 is the slant-range coordinate and τ is the delay-time. On the other hand, cross-range scaling requires the estimation of the modulus of the target's effective rotation vector, Ω, since the cross-range scaling is performed as x 1 = c 2 f 0 Ω ν . In this section an algorithm to obtain a fully scaled ISAR image is described.

Under the assumption that the target rotation vector is constant within the CPI, the chirp rate (which corresponds to the variation rate of the Doppler frequency of a scatterer) produced by the scattering centers is related to the modulus of the target rotation vector. Therefore, each scattering center carries information about the modulus of the target rotation vector through its chirp rate. To demonstrate this we will consider the radial motion compensated signal expressed in Equation (3.11).

The phase term in the integral can be approximated with a second order Taylor polynomial expansion, as follows:

3.43 ϕ ( f , n ; x 1 , x 2 ) = j 4 π f c ( x 2 + x 1 Ω n T R 1 2 x 2 Ω 2 ( n T R ) 2 )
so that the echo of an ideal scatterer located at ( x 1 k , x 2 k ) , with finite reflectivity function f ( x 1 , x 2 ) = σ k δ ( x 1 x 1 k , x 2 x 2 k ) , can be written as follows:
3.44 S C ( f , n ; x 1 k , x 2 k ) = W ( f , n ) σ k e ( j 4 π f c ( x 2 k + x 1 k Ω n T R 1 2 x 2 k Ω 2 ( n T R ) 2 ) )

After applying an FT to the signal in Equation (3.44) along the frequency domain, the range compressed signal is obtained, as shown in Equation (3.45):

3.45 S C ( τ , n ; x 1 k , x 2 k ) = B σ k sin c ( B ( τ 2 c x 2 k ) ) rect ( n N ) e j 4 π f 0 c ( x 2 k + x 1 k Ω n T R 1 2 x 2 k Ω 2 ( n T R ) 2 )

It should be noted that the phase term in Equation (3.45) represents the phase of a chirp signal, where μ k = 2 f 0 c x 2 k Ω 2

is usually referred to as chirp rate. If a method for perfectly estimating the chirp rate of a given scatterer was available, the following equation could be exploited to estimate the modulus of the effective target rotation vector:
3.46 μ = 2 f 0 c x 2 Ω 2

In fact, as the scatterer range coordinate x2 can be readily obtained by measuring the delay-time τ, the effective rotation vector can be obtained by inverting Equation (3.46), as follows:

3.47 Ω = c 2 f 0 x 2 μ

In practice, a scatterer chirp rate, as well as its range, must be estimated from the received data. Therefore, the estimation of the effective rotation vector magnitude would, in general, be affected by an error. Different techniques for estimating target scattering center chirp rates making use of atomic decomposition [15], CLEAN technique [13, 8] or based on the IC method have been proposed in [7].

Specifically, the method proposed in [7] is considered here, as it has been proven to be the most effective. According to such a method, the signal components of the scattering centers are processed by means of a polynomial Fourier transform (PFT) in order to estimate the chirp rate. The use of such a method, namely local polynomial Fourier transform (LPFT), requires the solution of an optimization problem [14]. It has been largely proven in literature that the image contrast (IC) and the image entropy (IE) are good indexes to assess the quality of an image. Therefore, the proposed algorithm estimates the scatterers chirp rate that maximize the IC.

To make the estimation more accurate and robust, the chirp rates of a number of target scatterers can be measured together with their ranges. The problem of estimating the effective rotation vector magnitude is then transformed into a problem of estimating the slope of a straight line that fits the scatterplot generated by the set of range and chirp rate estimates. One way of solving this problem is to apply a least square error (LSE) approach [6]. The mathematical problem can be set as follows:

3.48 μ ^ k = a x 2 k + b + ε k
where a = 2 f 0 c Ω 2 , b is the μ-intercept of the line and μ ^ k , x 2 k and k are the chirp rate estimate, the range estimate and the estimation error for the kth scatterer, respectively. Both a and b are unknown values. The LSE problem and its solution for the estimation of a are given in Equation (3.49) and Equation (3.50).
3.49 ( a ^ , b ^ ) = arg min a , b k = 1 K ε k 2
3.50 a ^ = K K = 1 K μ ^ k x 2 k k = 1 K μ ^ k k = 1 K x 2 k K k = 1 K ( x 2 k ) 2 ( k = 1 K x 2 k ) 2

For the sake of clarity, the algorithm steps can be summarized as follows:

  1. Image segmentation
  2. Chirp rate estimation via LPFT
  3. Effective rotation vector estimation

In detail, after motion compensation, an ISAR image obtained by means of the RD technique can be modeled as the composition of several sinc-like functions centered at each scatterer location in the IPP, as shown in Equation (3.51).

3.51 I ( τ , v ) = C k = 1 K σ k sin c [ B ( τ 2 x 2 k c ) ] sin c [ T o b s ( v 2 f 0 Ω x 1 k c ) ]
where K is the number of the scattering centers. After taking the absolute value of the ISAR image, a threshold is applied in order to extract the scattering centers. A number of connected groups of pixels are then found by means of a clustering procedure (segmentation task). In this way the backscattered signal coming from a specific scatterer can be selected. Its expression can be written as follows:
3.52 S C ( τ k , n ) = σ k e j 4 π f 0 c x 2 x e [ j 4 π f 0 c ( x 1 k Ω n T R x 2 k 2 Ω 2 ( n T R ) 2 ) ]

By applying a second order LPFT [14] a complex three variables ISAR image can be obtained in the delay time/Doppler/chirp rate domain as follows:

3.53 S C ( τ k , v , μ ) = σ k e j 4 π f 0 c x 2 x n e [ j π ( ( v + 2 f 0 c x 1 k Ω ) n T R + 1 2 ( μ 2 f 0 c x 2 k Ω 2 ) ( n T R ) 2 ) ]

From Equation (3.53) it is possible to estimate the Doppler frequency, νk, and the chirp rate, μk, by maximizing the IC of Sc(τk, ν, μ) with respect to ν and μ.

From the estimated value of μk and the scatterer location τk the modulus of the effective rotation vector can be estimated via the LSE approach described in Equation (3.50).

3.3  ISAR Parameters Setting

In a realistic scenario, the signal to work with is fully digital so all the variables introduced in the signal modeling are discrete and can be expressed as:

3.54 f = f 0 + m f Δ f m f = 0 , , M f 1 τ = m Δ τ m = 0 , , M 1 v = n v Δ v n v = 0 , , N v 1

Since the relationship between (f, n) and (τ, ν) is given by a Fourier transformation, the sampling intervals in the image domain depend on the ISAR signal spectral occupancy, as follows:

3.55 δ τ = 1 M Δ f = 1 B δ v = 1 N T R = 1 T o b s
where M = Mf and N = Nν. As a consequence, the sampling interval along the range and cross-range axes can be derived as in Equation (3.56)
3.56 δ r = c 2 B δ c r = c 2 f 0 Ω T o b s

Moreover, the above-mentioned sampling operation in the ISAR signal domain leads to an ambiguity in the ISAR image domain that must be taken into account in the radar design. Specifically, the non-ambiguity region in the cross-range and range dimension, i.e., Δ x 1

and Δ x 2 , can be expressed by considering the discrete spatial frequency domain shown in Figure 3.10, as follows:
3.57 Δ x 1 = c 2 f 0 Ω T R Δ x 2 = c 2 Δ f

In order to avoid ambiguity, the relations in Equation (3.57) define the maximum sizes of the target. Any target bigger than the non-ambiguity region will be seen as folded within the image making target recognition more difficult. It is worth pointing out that Δ x 1

depends also on the modulus of the effective rotation vector, Ω, which, as stated above, cannot be a priori predicted leading to a more difficult design of an ISAR system.

Non-ambiguity regions

Figure 3.10   Non-ambiguity regions

Moreover, considering Equation (3.24) some remarks about spatial resolutions can be made. Specifically, the range resolution depends on the transmitted signal bandwidth and can be directly controlled in the radar design process, since it can be considered as a known parameter. On the other hand, it is evident that the spatial resolution in the cross-range dimension depends on the total aspect angle variation, Δϑ = ΩTobs. This results in a parameter that is not controllable and cannot be selected during the radar design process, unless making some hypothesis about the target motion

3.4  Conclusions

The main concepts and theory behind ISAR imaging have been treated in this chapter. More specifically, the ISAR geometry and the received signal model have been defined. Then, mathematical details have been introduced in this chapter to support the main concepts behind ISAR technique and to provide the basis for implementing basic ISAR imaging algorithms. The concepts presented in this chapter mainly apply to monostatic active radar system. Nevertheless, the ISAR signal model formulation presented in this chapter is independent of the transmitted waveform; therefore, as it will be clarified in the following chapters, the same concepts apply to active and passive bistatic radar systems.

References

F. Berizzi, M. Martorella, B. Haywood, E. Dalle Mese, and S. Bruscoli. A survey on ISAR autofocusing techniques. Image Processing, 2004. ICIP ’04. 2004 International Conference on, 1:9–12, Oct. 2004.
F. Berizzi, E.D. Mese, M. Diani, and M. Martorella. High-resolution ISAR imaging of maneuvering targets by means of the range instantaneous doppler technique: modeling and performance analysis. Image Processing, IEEE Transactions on, 10(12):1880–1890, Dec 2001.
G. Clifford Carter. Time delay estimation for passive sonar signal processing. Acoustics, Speech and Signal Processing, IEEE Transactions on, 29(3):463–470, Jun 1981.
B. Haywood and R. J. Evans. Motion compensation for ISAR imaging. Proceedings of ASSPA 89, pages 113–117, 1989.
C.V.J. Jakowatz, D.E. Wahl, P.H. Eichel, D.C. Ghiglia, and P.A. Thompson. Spotlight-Mode Synthetic Aperture Radar: A Signal Processing Approach: A Signal Processing Approach. Springer, 2011.
Steven M. Kay. Fundamentals of Statistical Signal Processing: Estimation Theory. Signal Processing. 1993.
M. Martorella. Novel approach for ISAR image cross-range scaling. Aerospace and Electronic Systems, IEEE Transactions on, 44(1):281–294, 2008.
M. Martorella, N. Acito, and F. Berizzi. Statistical CLEAN technique for ISAR imaging. IEEE Tr. on Geoscience and Remote Sensing, 45(11):3552–3560, 2007.
M. Martorella and F. Berizzi. Time windowing for highly focused ISAR image reconstruction. Aerospace and Electronic Systems, IEEE Transactions on, 41(3):992–1007, 2005.
M. Martorella, F. Berizzi, and B. Haywood. Contrast maximisation based technique for 2-D ISAR autofocusing. Radar, Sonar and Navigation, IEE Proceedings, 152(4):253–262, 2005.
R. Gray Parker and Ronald L. Rardin. Discrete Optimization. Academic Press Professional, Inc., San Diego, CA, 1988.
B.D. Steinberg. Radar imaging from a distorted array: The radio camera algorithm and experiments. Antennas and Propagation, IEEE Transactions on, 29(5):740–748, Sept. 1981.
J. Tsao and B.D. Steinberg. Reduction of sidelobe and speckle artifacts in microwave imaging: The CLEAN technique. IEEE Tr. on Antennas and Propagation, 36:543–556, 1988.
Yongmei Wei and Guoan Bi. Efficient analysis of time-varying multicomponent signals with modified LPTFT. EURASIP J. Adv. Sig. Proc., 2005(8):1261–1268, 2005.
O.A. Yeste-Ojeda, J. Grajal, and G. Lopez-Risueno. Atomic decomposition for radar applications. Aerospace and Electronic Systems, IEEE Transactions on, 44(1):187–200, Jan. 2008.
Search for more...
Back to top

Use of cookies on this website

We are using cookies to provide statistics that help us give you the best experience of our site. You can find out more in our Privacy Policy. By continuing to use the site you are agreeing to our use of cookies.