Real-time synthetic aperture radar (SAR) processing ... - IEEE Xplore
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 30, NO. 4, JULY 1992
714
Real-Time Synthetic Aperture Radar (SAR) Processing with a New Subaperture Approach Albert0 Moreira
Abstruct- A time domain subaperture approach has been developed at DLR which is suitable for real-time synthetic aperture radar ( S A R ) processing and produces high quality, full resolution images. The real-time subaperture algorithm is based on an approximation of the phase history correction in each subaperture with a simple linear correction, which can be carried out by an up/down-conversion of the received signal followed by a moving average operation. The characteristics of the impulse response function are improved considerably by means of a frequency overlap of the subapertures and become comparable to a conventional matched filter response. Special algorithms for high quality SAR processing, such as auto-focus, Doppler centroid estimation, and range migration correction, can be efficiently implemented. High flexibility is achieved for multilook processing and no secondary range compression is needed, even for high squint angles. Dedicated hardware is under development for the DLR airborne SAR (E-SAR) System. Keywords-Synthetic Aperture Radar Processing (SAR), Subaperture Processing, Real-Time SAR Processing.
I. INTRODUCTION
S
AR processing is based on a two-dimensional correlation of the backscattered signal with a space-variant reference function. The image coordinates are azimuth (along track) and range (across track). The commonly used processing algorithms [3], [6], [9], [18], [19], [20], [23] use a two-dimensional Fast Fourier Transform (FFT) (w - IC approach) or 1-D FFT’s (range-Doppler approach) to perform the correlation process efficiently and to obtain a well focused image. In these cases, a long FFT is needed for efficient block processing and large memory capacity is required, especially in 2-D FFT processing. Since the azimuth modulation of the backscattered signal is dependent on range and follows the range migration line, an interpolation and a reference function update or phase correction must be carried out in order to achieve high quality imagery. In addition to this, the platform motion errors (e.g., velocity and/or squint angle variation) may require an update of the reference function in azimuth direction, which will increase the computation time considerably. Basic subaperture processing consists of dividing the received SAR signal (in range or in azimuth) into subapertures, each corresponding to a small part of the total signal bandwidth. The data in each subaperture is correlated with the matched reference function in order to obtain the low resoluManuscript received August 23, 1991; revised February 24, 1991. The author is with the German Aerospace Research Establishment (DLR), Institute for Radio Frequency Technology, D-8031 Oberpfaffenhofen, Germany. IEEE Log Number 9200972.
tion images. By adding all the subapertures coherently, a high resolution image is obtained. Several subaperture techniques have been used in the past [2], [17], [21], [22], [25], and [26]. The most common technique is the step transform. In this algorithm, the subapertures are formed by multiplying the received SAR signal with a set of overlapped complex conjugated reference signals and performing a frequency analysis with short FFT’s on the results of the multiplications. After applying a time delay for the subaperture signals, in order to correct their relative positions, and after a phase multiplication to provide a phase continuity along the subapertures, a second FFT is used to obtain the high resolution image. One disadvantage of the step transform algorithm occurs due to the undersampling of the subapertures that deteriorates the image quality. If the subapertures are not undersampled, the computational efficiency of this algorithm is decreased considerably. The new Real-Time Subaperture algorithm (RTS algorithm), in its original form, does not make use of the FFT and requires a reduced number of subapertures to achieve the desired geometric resolution. The memory and computational requirements can be reduced by undersampling the subapertures without deteriorating the image quality. Since the bandwidth of the subapertures is much smaller than the bandwidth of the original signal, the ur,dersampling is easily performed by picking one sample from each group of Nunder samples, where NTLnder is the undersampling factor. Additional flexibility to multilook processing is obtained due to the subaperture structure and special algorithms for high quality SAR processing (such as auto-focus and Doppler centroid estimation) can be efficiently implemented. A detailed comparison between the computational effort of the RTS algorithm and other conventional SAR processing algorithms can be found in [14]. Briefly, the computational requirement of the RTS-algorithm is comparable to the basic range-Doppler algorithm using FFT for a time-bandwidth product ( T B P ) greater than 100. For a T B P less than 100 (some airborne cases), the RTS algorithm requires less computational effort than the range-Doppler algorithm. The basic subaperture principle for SAR processing will be introduced in Section 11. The subaperture processing will then be simplified in order to reduce the computational requirements (Section 111). A linear approximation for correction of the phase variation in the subapertures is used, corresponding to unfocused subaperture processing. A periodic phase error in the correlation process arises from the linear phase approximation. This leads to paired echoes in the final, high resolution impulse response function. Through an optimization of the
0196-2892/92$03.00 0 1992 IEEE
715
MOREIRA: REAL-TIME SAR PROCESSING
processing parameters (frequency overlap , integration time and amplitude weighting of the subapertures) the paired echoes are suppressed more than 35 dB. The basic equations for one-dimensional SAR processing using the RTS algorithm are given in Section IV. The use of improved multilook techniques [14] with the RTS algorithm, as well as the implementation of auto-focus and Doppler centroid estimation into the RTS algorithm, are briefly described in this section. Processing with the RTS algorithm will be extended in Section V to the twodimensional case, in order to accommodate range migration. It will be shown that the linear part of range migration can be corrected prior to azimuth processing, so that only the quadratic term of range migration must be considered during azimuth processing. The simulation results in Section VI, using the parameters for the DLR airborne real-time SAR processor, show that no visible deterioration of the impulse response function is produced, even for high squint angles. Two raw data sets of the E-SAR System were used for processing with the RTS algorithm. The results of the ensuing image analysis are presented in Section VII. The image quality requirement for the E-SAR system is achieved using autofocus. The principal advantage of the hardware implementation with the RTS algorithm is the subaperture structure with reduced computation, which allows a parallel hardware configuration using commercial digital signal processors to achieve real-time SAR processing.
for the correlation with the reference functions. From (1)-(3) we obtain: I
S I
sinc[b. TI . ( t - t o ) ]
(4)
where sinc(n:) is the sin(z)/z-function and to is a time constant corresponding to the peak position of the IRF. The sinc-function has a low geometric resolution corresponding to the integration time of each subaperture. Performing the coherent addition of the subapertures with the phase correction 4 = 2.71.b.Tl.t(first term in the summation of (4)) results in:
where B is the processed signal bandwidth and T the total integration time (T = T l . ( 2 . N l + l ) )The . same IRF is achieved as in conventional matched filter processing. Equations (1)-(3) and 5 provide the basis for typical subaperture processing. Variations of this basic processing are used in the polyphase [2] and in the step transform [21], [22], [25], [26] processing algorithms. Other subaperture processing algorithms for medium geometric resolution and reduced image quality (e.g., [5],[7]) have been developed but will not be considered here, since the main objective is to achieve the same image quality as with the traditional matched filter approach. Three principal advantages are obtained by using subaperture processing: 11. THE BASK SUBAPERTURE PRINCIPLE 1. Flexibility to multilook processing. A group of subaperConsider the one-dimensional frequency modulated signal tures can be added incoherently in order to form the s ( t ) and the corresponding matched filter h , ( f ) . The imlooks without changing the basic processing parameters. pulse response function (IRF) using conventional time domain 2. Parallel hardware structure. This enables the use of subaperture processing (coherent look processing) can be a dedicated processing unit for each subaperture, thus formulated as follows: decreasing processing time. 3. Auto-focus capability. In azimuth processing, the time delay between the signals of adjacent subapertures can Ifo(t>l = hOn( t )@ s n ( t ) (1) be measured by means of a correlation of the amplitude ~ n ~ ~ . ~ l of the subaperture signals. The deviation of the measured where ( 2 . N1 1) is the number of subapertures, each with time delay from the nominal value T1 is then used to an integration time TI and @ stands for the one-dimensional determine the variation of the platform velocity. correlation. Assuming the signal s ( t ) with a linear modulation rate and h, as the complex conjugate, time reversed of s ( t ) , 111. ANALYSIS OF SUBAPERTURE PROCESSING the corresponding signals s n ( t ) and ho,l( t )in each subaperture The basic objective of this analysis is to simplify the n are given by processing in each subaperture, so that the computational requirement can be reduced. An optimization of the processing S , ( t ) = s ( t ) ( n - +) 5 t 5 ( n +) T l , (2) parameters should be carried out in order to achieve the same s n ( t )= 0 elsewhere image quality as with a matched filter. and similarly Performing a linear approximation to the phase of the subapertures results in: h,,(t) = h ( t ) . ( n - +) 5 t 5 (71 $) TI, (3) elsewhere. hOn( t )= 0. 4,,(t) = Tl . b . [ 2 . n . ( t - t o )- n2 . TI] A bank of FIR (finite impulse response) filters can be used 71= - s, to obtain the signals s, ( t ) , while the reference functions h,,(t) for each subaperture can be obtained directly from the full aperture reference function h,(t). The advantage of using the FIR filters is that the subaperture signals can be where (ij,, is the phase history of the reference function (see undersampled after the filtering, thus reducing the data amount Fig. 1) and rect[t/Tl] is the rectangular function with duration
1
+
i
+
i
+
'
'
1
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 30, NO. 4, JULY 1992
716
2.5L'
3
*ol
a
O L 0
'
"
I
'
"
'
"
"
"
'4
i
I
1
0.5 0.0
200
400 600 TIME (arbitrary scale)
800
1000
Fig. 1. Linear approximation of the phase correction for subaperture processing. The integration time of each subaperture is T I .
7'1 . This approximation leads to a phase error that has its maximum (q5e7naz) at t = n + to. For phase errors & ,,, smaller than 180", the IRF can be approximated to
0
{
7r
. B . [(t - to) - k . ( 2 . N,+l)]
1,
. (7)
The IRF consists of a main response ( k = 0) with high geometric resolution superposed by an infinite summation of paired echoes, whose amplitudes are given by the Bessel function J k of the first kind and kth order. From this equation, the loss, peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) can be calculated as a function of the phase error. For qbemaz = go', a PSLR of -9.6 dB and a ISLR of -4.6 dB are achieved, both values not sufficient for high image quality. Decreasing the phase error to 3" will almost eliminate the paired echoes but the number of subapertures will be increased approximately by a factor of 6, according to the following relation:
where fix(.) is the integer value operator. Due to the great number of subapertures for a requirement of qh,, = 3", the processing using this scheme becomes inefficient. In the following, a more effective way will be derived for improving the PSLR and ISLR. Starting from (1) and using the linear approximation for the phase history (6), a modified expression of (4) is obtained [14]: If(t)l =
'
I4t)l.
IR(t)l is a periodic sin(:c)/z-function and is given by:
The following properties can be associated with lR(t)l :
(9)
200
400 600 TIME ( a r b i t r a r y scale)
800
1000
Fig. 2. IRF and product terms for subaperture processing with periodical phase error. I.i(t)l is the low resolution response of the subapertures, IR(t)l is the periodical & i ( . i , ) / . r -function and I f ( t ) l is the high resolution impulse response function.
+
1. Period of ( 2 . NI l ) / B . 2. High resolution in time (zl/B). 3 . Independence of the phase error d,,,,. The second term of (9) is given by
IA(t)l = sinc
I
i
! ul I
"
linear approximation
~
d m .{[C(tl)+ C(t2)I2+ [S(tl)+ S ( ~ P2 })11]2
(11) where S(E)and C(t) are the complex Fresnel integrals. Some important properties of IA(t)l are: 1. Low resolution in time (= TI). 2. Independence of the number of subapertures. 3. PSLR and ISLR are related with 4,,,. Fig. 2 shows the IRF for YO" phase error as well as the product terms IR(t)l and IA(t)l.As can be seen from this figure, the amplitude of the paired echoes is given by the amplitude of the minima of the low resolution function IA(t)l.Two solutions were found to be efficient for decreasing the paired echoes. Improvement of the shape of the low resolution function IA(t)l. This can be achieved by using an amplitude weighting in the correlation of each subaperture, so that the effective amount of phase error is reduced and the amplitude of the minima of IA(t)l is attenuated. By using a strong amplitude weighting, the geometric resolution of IA(t)l is deteriorated by M 30% and the distance of the minima of IA(t)l apart from the mainlobe is increased proportionally. Use of a frequency overlap between the subapertures, so that the period of IR(t)l is increased and the position of the first paired echo around the mainlobe coincides with the first minimum of IA(t)l (obtained by amplitude weighting). Using a triangular weighting for the correlation in the subapertures and allowing a phase error 4,,, in each subaperture of 28.8"', the characteristics of the IRF are optimized if 'The phase error in each subaperture should be limited to 28.8' in order to maintain processing loss lower than 0.6 dB. The triangular weighting for the correlation in each subaperture will be adopted for the RTS approach, since it can he easily implemented by two moving average operations (correlation of two retangular functions). The adoption of the triangular weighting and the condition of 0, ,. 5 28.8' imposes an overlap of the subapertures that must be greater than 21% for sufficient attenuation of the sidelobes.
717
MOREIRA: REAL-TIME SAR PROCESSING "
0
"
"
'
I
-20
0
Fig. 3.
"
PiLR
=
ISLR
=
Loss
34 dB 14.4 dB
= 0.6 d B
1000 TIME (arbitrary scale)
500
1500
IRF of the RTS-algorithm for 21% frequency overlap (the response of one subaperture is represented by the dotted line).
the subapertures are overlapped by 21%, as shown in Fig. 3. A n additional weighting w(n) (Hamming, CY = 0.6) was introduced in the coherent summation of the subapertures in order to decrease the sidelobes of the IRF as much as the amplitude of the paired echoes. If better values of PSLR and ISLR are required, an overlap of 57% will place the first paired echo in the second amplitude minimum of IA(t)l and the PSLR will be attenuated to -45 dB for a Hamming weighting with a = 0.54.
IV. THERTS ALGORITHM Assuming a linear approximation of the phase history, a triangular function for correlation of the subapertures and a frequency overlap ovl, we obtain:
If ( t ) /=
(12) where iVlobl = N1/(1 - ovl) and n' = n . (1 - owl). The double integral in (12) corresponds to a double correlation of the signal s ( t ) with a rectangular function, where a frequency shift (up/down conversion of the subapertures) is performed before and after the correlation. We define the first frequency shift ( g l l ) and the second frequency shift (h,) before coherent summation of the subapertures as follows:
g n ( t ) = exp [ - j . 2 . 12'. TI . b . (t - t o ) ] .
(13)
Fig. 4. Basic block diagram of the RTS algorithm. In this figure, five subapertures are added coherently, whereby the high resolution IRF is obtained (no multilook processing is used).
Fig. 4 shows the basic block diagram of the impulse compression using the RTS algorithm. First, the received signal s ( t ) is multiplied with the frequency shifts g, , which shift the signal s ( t ) from the desired band for each subaperture up or down to baseband. The results of the frequency shift are correlated twice with a rectangular function by means of two moving average operations, where only two real additions and two real subtractions are needed per correlated point. This operation can also be interpreted as an integration with duration TI, as represented in Fig. 4. Before the second frequency shift is carried out, a time shift ri' . Tl is performed with each subaperture. This shift compensates for the time displacement of the subapertures, which are derived from different parts of the full aperture. Then the subapertures are up and down converted to the initial frequency value by means of hlL.An additional phase term b . a'* . T: is provided by h, (see 14) in order to ensure a phase continuity over the subapertures. Finally, the subapertures are added coherently, whereby the high resolution IRF is obtained. When designing an RTS processor, optimized values should be used for the maximum phase error (4e,naz5 28.8') and for the minimum overlap (ovl 2 0.21). Obviously, these requirements can be met for range processing and also for azimuth processing over all distances (including a possible squint angle). Then the integration time TI of each subaperture and the total number of subapertures ( 2 . Nloz,l 1) are calculated as follows:
+
h,(t) = ~ ( n. e)z p [ j . b . TL'. ( 2 . TI . ( t - t o )- 72' . T:)]. (14) Using (12)-(14) given by
the final equation of the RTS-algorithm is where B is the desired signal bandwidth to be processed. The processing itself is performed according to (15). For instance, the time integration of each subaperture is done over , a discrete number of points, i.e. fix[Tl . f s / ( N u n d e r ) ]where
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 30, NO. 4, JULY 1992
718
fs is the sampling frequency. The discrete integration does not influence the processing performance since the integer number of points is greater than 15 for all available SAR systems. Although the RTS algorithm works in the time domain, the computational requirement is low, since the time correlation in each subaperture is substituted by a moving average operation. The determination of the reference function for processing is simplified to the calculation of the functions gT1, h,, and the parameters TI and n’ . T I . The computational requirement can be reduced further, if the frequency shifts gn and h, are set to a submultiple of the sampling frequency. Then, short FFT’s can be used for the formation and the summation of the subapertures [13]. For a time-bandwidth product greater than 100 (more than % 10 subapertures are needed), the use of FFT’s in the RTSalgorithm is preferred. In the following, three specific algorithms for high resolution SAR processing will be discussed: multilook, autofocus and Doppler centroid estimation. We show that these algorithms can be incorporated into the RTS algorithm without additional operations, such as FFT.
enough accuracy to focus the SAR image [lo]. In the case of subaperture processing, a cross-correlation between adjacent subapertures is carried out to determine the relative time shift of the subapertures. If no velocity variations occur, the nominal time shift 71’ . TI will cause the proper positioning of the subapertures. The measured displacement Ax between the subapertures is then used to estimate the correct Doppler rate IC, = for the azimuth processing as follows
4
IC, =
+ v,,,
. Ax ’
(1 - owl) . TI
where V,,, and ICatir are the erroneous velocity and Doppler rate, respectively, assumed for processing. In comparison with the traditional auto-focus with multilook [lo], a better sensitivity to high frequency velocity variations is achieved by the subaperture processing due to the low integration time of the subapertures. Obviously, the auto-focus output degrades proportionally to the integration time, since the crosscorrelation result broadens as the integration time decreases. In an airborne SAR system, the auto-focus technique can be used to determine the aircraft velocity variations in realA. Multilook Processing For RTS processing with Lo looks, a number of (2.N1,,1+ time and to update the PRF (pulse repetition frequency) in a l ) / L o subapertures should be added coherently in order to continuous manner. This ensures a well focused image with obtain each look response. Then, the results of these additions geometric fidelity, since the azimuth pixel spacing is held are detected and added incoherently to improve the radiometric constant by changing the PRF. A prediction algorithm for the resolution. A n optimization of the look overlap to achieve an determination of the PRF (possibly based on a Kalman filter) improved radiometric resolution can be found in 1131. Only is required due to the delay associated with the autofocus the last stage of the RTS processor (coherent and incoherent computation. addition of the subapertures) must be changed to vary the number of looks and the look overlap. C. Doppler Centroid Estimation An IML (improved multilook) technique [13] can also Doppler centroid estimation can be carried out directly in be implemented efficiently in the RTS processor. The IML azimuth processing because the subaperture domain (index technique is based on the formation of looks with different bandwidths. The looks with larger bandwidth contribute to an 7~’)corresponds to a frequency domain with a low frequency improvement of the overall geometric resolution, while the resolution = P R F / ( 2 . NIo7./ 1). Then, by averaging looks with smaller bandwidth improve the overall radiomet- the signal intensity in each subaperture over the azimuth ric resolution. Multilook processing with the IML technique time or over a group of range bins, an estimation of the improves the radiometric resolution of SAR images, since it azimuth antenna pattern can be made [14]. It is assumed makes a more efficient use of the available signal bandwidth. that the number of subapertures is large enough (> 10) For processing with the IML technique using the correlation that the azimuth antenna pattern can be considered constant in frequency domain (e.g., 1111, [24]), approximately 50% within a subaperture. After interpolating the results of the more computation time is required than with the traditional averaging in order to obtain the desired frequency accuracy for multilook technique. This occurs, since a number of L L looks determination of the Doppler centroid, conventional methods with larger bandwidth and a number of Ls looks with smaller (e.g., energy balance method [lo]) can be used for efficient bandwidths must be formed in the frequency domain, followed estimation of the Doppler centroid. The time domain structure of the RTS algorithm (no block by ( L L L s ) inverse FFT’s. Then L.5 extra inverse FFT’s processing is used) imposes no restrictions on the update rate are needed, leading to the increased computation time. Using of the Doppler centroid. This feature is an important point for the RTS algorithm, the extra computation time is reduced to a minimum. Only extra additions are needed to create the looks the accurate processing of future spaceborne SAR systems. with larger and smaller bandwidth (coherent addition) and to For the X-SAR system, which will be launched on the Space Shuttle in 1993, an azimuth Doppler centroid variation of up form the final image (incoherent addition). to 280 Hz/s will occur due to the instability of the platform [4]. Since the required accuracy of the Doppler centroid is f 5 0 Hz, B. Auto-Focus an update rate of 0.2 s is necessary to obtain a high quality Auto-focus is commonly implemented in azimuth process- image. The RTS algorithm can accommodate these updates ing in order to determine the Doppler frequency rate with without decreasing its computational efficiency.
+
+
MOREIRA: REAL-TIME SAR PROCESSING
V. EXTENSION TO 2-D
719
TABLE I PROCESSING PARAMETERS FOR THE SIMULATION OF THE IRF IN C-BAND(SQUINT ANGLE= 8", AIRCRAFT VELOCITY = 70 mis AND SHORTEST RANGETO POINTTARGET= 3000 m)
PROCESSING
The extension of the 1-D RTS algorithm to a 2-D RTS algorithm means that the functions g, and h, are twodimensional and that the triangular function for the correlation should be substituted by a pyramidal function. However, this 2-D implementation will have limitations on the correction of the linear and quadratic range migration and will require a large memory for the correlation process. Therefore, the 1-D RTS algorithm (as described in the previous section), with a time domain correction of the range migration, will be used to accomplish the 2-D nature of the azimuth reference function. The adopted processing steps are the same as for a 1-D time domain SAR processor (range compression, corner turn operation, azimuth processing, and detection). The correction of range migration is divided into two steps: A) range walk correction, which is performed during range processing and corner turn operation and B) range curvature correction, which is carried out by an interpolation of samples prior to the coherent summation of the subapertures in the azimuth processing.
~~
Simulation Parameters
range orocessine
azimuth processing
processed signal bandwidth
100 MHz
100 Hz
_____
reference function duration
5
1.7 s 0.104 s
PS
integration time of each subaperture
0.178 p s
subaperture overlap
57%
number of subapertures
63
21% 20
"I
A. Range Walk Correction
The fractional part of the range walk can be easily corrected during range processing by adding a linear phase term before the coherent summation of subapertures is carried out. Considering the periodic sin(z)/z-function R ( t )in (9) and including the linear phase term, 41 = 2 . n' . (1 - out) . TI . b . At, where At is the desired time shift, we obtain: sinj7r . B . ( t - At - to)] IR(t)l = sin ?r-B.(t-At-t,)]
[
Fig. 5. Simulated point target response for the E-SAR system of DLR (C-band, 8' squint angle).
VI.
SIMULATION OF THE
IRF
(19)
By using the parameters for the C-band mode of the airborne
Since the IRF is given by IR(t)l . IA(t)l, where IA(t)l is a slowly varying function (refer to (9) and (ll)), the obtained time shift of the range compressed data is At. This first step of the range migration correction ensures an exact interpolation of the data in range direction which can be updated from sample to sample, since no block processing is used by the RTS algorithm. The update of the range walk becomes necessary when the Doppler centroid varies with range (e.g., large swath width and/or small look angle). The integer part of the range walk is corrected by a shift of range samples during the corner turn operation. An additional deskew of the data is performed with the final image (possible during slant to ground range conversion) in order to compensate for the integer sample shift after the range processing.
integrated, a point target simulation and processing was carried out in order to estimate the performance of the RTS algorithm. The critical points in processing the E-SAR data are the correction of the motion errors due to the instability of the small DO-228 aircraft and the high range walk that occurs when flying with strong crosswinds. The motion errors of the E-SAR system can be corrected efficiently by the reflectivity displacement method [15] and will not be covered here. We analyze below the performance of the RTS algorithm for high drift angles. A detailed hardware simulation of the real-time processor including quantization and linearity errors due to the finite bit length in the processing can be found in [12]. A squint angle of 8' was used for the simulation. In this case, the total range walk within the azimuth integration time is greater than 30 range bins, causing an azimuth time-bandwidth product
714
Real-Time Synthetic Aperture Radar (SAR) Processing with a New Subaperture Approach Albert0 Moreira
Abstruct- A time domain subaperture approach has been developed at DLR which is suitable for real-time synthetic aperture radar ( S A R ) processing and produces high quality, full resolution images. The real-time subaperture algorithm is based on an approximation of the phase history correction in each subaperture with a simple linear correction, which can be carried out by an up/down-conversion of the received signal followed by a moving average operation. The characteristics of the impulse response function are improved considerably by means of a frequency overlap of the subapertures and become comparable to a conventional matched filter response. Special algorithms for high quality SAR processing, such as auto-focus, Doppler centroid estimation, and range migration correction, can be efficiently implemented. High flexibility is achieved for multilook processing and no secondary range compression is needed, even for high squint angles. Dedicated hardware is under development for the DLR airborne SAR (E-SAR) System. Keywords-Synthetic Aperture Radar Processing (SAR), Subaperture Processing, Real-Time SAR Processing.
I. INTRODUCTION
S
AR processing is based on a two-dimensional correlation of the backscattered signal with a space-variant reference function. The image coordinates are azimuth (along track) and range (across track). The commonly used processing algorithms [3], [6], [9], [18], [19], [20], [23] use a two-dimensional Fast Fourier Transform (FFT) (w - IC approach) or 1-D FFT’s (range-Doppler approach) to perform the correlation process efficiently and to obtain a well focused image. In these cases, a long FFT is needed for efficient block processing and large memory capacity is required, especially in 2-D FFT processing. Since the azimuth modulation of the backscattered signal is dependent on range and follows the range migration line, an interpolation and a reference function update or phase correction must be carried out in order to achieve high quality imagery. In addition to this, the platform motion errors (e.g., velocity and/or squint angle variation) may require an update of the reference function in azimuth direction, which will increase the computation time considerably. Basic subaperture processing consists of dividing the received SAR signal (in range or in azimuth) into subapertures, each corresponding to a small part of the total signal bandwidth. The data in each subaperture is correlated with the matched reference function in order to obtain the low resoluManuscript received August 23, 1991; revised February 24, 1991. The author is with the German Aerospace Research Establishment (DLR), Institute for Radio Frequency Technology, D-8031 Oberpfaffenhofen, Germany. IEEE Log Number 9200972.
tion images. By adding all the subapertures coherently, a high resolution image is obtained. Several subaperture techniques have been used in the past [2], [17], [21], [22], [25], and [26]. The most common technique is the step transform. In this algorithm, the subapertures are formed by multiplying the received SAR signal with a set of overlapped complex conjugated reference signals and performing a frequency analysis with short FFT’s on the results of the multiplications. After applying a time delay for the subaperture signals, in order to correct their relative positions, and after a phase multiplication to provide a phase continuity along the subapertures, a second FFT is used to obtain the high resolution image. One disadvantage of the step transform algorithm occurs due to the undersampling of the subapertures that deteriorates the image quality. If the subapertures are not undersampled, the computational efficiency of this algorithm is decreased considerably. The new Real-Time Subaperture algorithm (RTS algorithm), in its original form, does not make use of the FFT and requires a reduced number of subapertures to achieve the desired geometric resolution. The memory and computational requirements can be reduced by undersampling the subapertures without deteriorating the image quality. Since the bandwidth of the subapertures is much smaller than the bandwidth of the original signal, the ur,dersampling is easily performed by picking one sample from each group of Nunder samples, where NTLnder is the undersampling factor. Additional flexibility to multilook processing is obtained due to the subaperture structure and special algorithms for high quality SAR processing (such as auto-focus and Doppler centroid estimation) can be efficiently implemented. A detailed comparison between the computational effort of the RTS algorithm and other conventional SAR processing algorithms can be found in [14]. Briefly, the computational requirement of the RTS-algorithm is comparable to the basic range-Doppler algorithm using FFT for a time-bandwidth product ( T B P ) greater than 100. For a T B P less than 100 (some airborne cases), the RTS algorithm requires less computational effort than the range-Doppler algorithm. The basic subaperture principle for SAR processing will be introduced in Section 11. The subaperture processing will then be simplified in order to reduce the computational requirements (Section 111). A linear approximation for correction of the phase variation in the subapertures is used, corresponding to unfocused subaperture processing. A periodic phase error in the correlation process arises from the linear phase approximation. This leads to paired echoes in the final, high resolution impulse response function. Through an optimization of the
0196-2892/92$03.00 0 1992 IEEE
715
MOREIRA: REAL-TIME SAR PROCESSING
processing parameters (frequency overlap , integration time and amplitude weighting of the subapertures) the paired echoes are suppressed more than 35 dB. The basic equations for one-dimensional SAR processing using the RTS algorithm are given in Section IV. The use of improved multilook techniques [14] with the RTS algorithm, as well as the implementation of auto-focus and Doppler centroid estimation into the RTS algorithm, are briefly described in this section. Processing with the RTS algorithm will be extended in Section V to the twodimensional case, in order to accommodate range migration. It will be shown that the linear part of range migration can be corrected prior to azimuth processing, so that only the quadratic term of range migration must be considered during azimuth processing. The simulation results in Section VI, using the parameters for the DLR airborne real-time SAR processor, show that no visible deterioration of the impulse response function is produced, even for high squint angles. Two raw data sets of the E-SAR System were used for processing with the RTS algorithm. The results of the ensuing image analysis are presented in Section VII. The image quality requirement for the E-SAR system is achieved using autofocus. The principal advantage of the hardware implementation with the RTS algorithm is the subaperture structure with reduced computation, which allows a parallel hardware configuration using commercial digital signal processors to achieve real-time SAR processing.
for the correlation with the reference functions. From (1)-(3) we obtain: I
S I
sinc[b. TI . ( t - t o ) ]
(4)
where sinc(n:) is the sin(z)/z-function and to is a time constant corresponding to the peak position of the IRF. The sinc-function has a low geometric resolution corresponding to the integration time of each subaperture. Performing the coherent addition of the subapertures with the phase correction 4 = 2.71.b.Tl.t(first term in the summation of (4)) results in:
where B is the processed signal bandwidth and T the total integration time (T = T l . ( 2 . N l + l ) )The . same IRF is achieved as in conventional matched filter processing. Equations (1)-(3) and 5 provide the basis for typical subaperture processing. Variations of this basic processing are used in the polyphase [2] and in the step transform [21], [22], [25], [26] processing algorithms. Other subaperture processing algorithms for medium geometric resolution and reduced image quality (e.g., [5],[7]) have been developed but will not be considered here, since the main objective is to achieve the same image quality as with the traditional matched filter approach. Three principal advantages are obtained by using subaperture processing: 11. THE BASK SUBAPERTURE PRINCIPLE 1. Flexibility to multilook processing. A group of subaperConsider the one-dimensional frequency modulated signal tures can be added incoherently in order to form the s ( t ) and the corresponding matched filter h , ( f ) . The imlooks without changing the basic processing parameters. pulse response function (IRF) using conventional time domain 2. Parallel hardware structure. This enables the use of subaperture processing (coherent look processing) can be a dedicated processing unit for each subaperture, thus formulated as follows: decreasing processing time. 3. Auto-focus capability. In azimuth processing, the time delay between the signals of adjacent subapertures can Ifo(t>l = hOn( t )@ s n ( t ) (1) be measured by means of a correlation of the amplitude ~ n ~ ~ . ~ l of the subaperture signals. The deviation of the measured where ( 2 . N1 1) is the number of subapertures, each with time delay from the nominal value T1 is then used to an integration time TI and @ stands for the one-dimensional determine the variation of the platform velocity. correlation. Assuming the signal s ( t ) with a linear modulation rate and h, as the complex conjugate, time reversed of s ( t ) , 111. ANALYSIS OF SUBAPERTURE PROCESSING the corresponding signals s n ( t ) and ho,l( t )in each subaperture The basic objective of this analysis is to simplify the n are given by processing in each subaperture, so that the computational requirement can be reduced. An optimization of the processing S , ( t ) = s ( t ) ( n - +) 5 t 5 ( n +) T l , (2) parameters should be carried out in order to achieve the same s n ( t )= 0 elsewhere image quality as with a matched filter. and similarly Performing a linear approximation to the phase of the subapertures results in: h,,(t) = h ( t ) . ( n - +) 5 t 5 (71 $) TI, (3) elsewhere. hOn( t )= 0. 4,,(t) = Tl . b . [ 2 . n . ( t - t o )- n2 . TI] A bank of FIR (finite impulse response) filters can be used 71= - s, to obtain the signals s, ( t ) , while the reference functions h,,(t) for each subaperture can be obtained directly from the full aperture reference function h,(t). The advantage of using the FIR filters is that the subaperture signals can be where (ij,, is the phase history of the reference function (see undersampled after the filtering, thus reducing the data amount Fig. 1) and rect[t/Tl] is the rectangular function with duration
1
+
i
+
i
+
'
'
1
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 30, NO. 4, JULY 1992
716
2.5L'
3
*ol
a
O L 0
'
"
I
'
"
'
"
"
"
'4
i
I
1
0.5 0.0
200
400 600 TIME (arbitrary scale)
800
1000
Fig. 1. Linear approximation of the phase correction for subaperture processing. The integration time of each subaperture is T I .
7'1 . This approximation leads to a phase error that has its maximum (q5e7naz) at t = n + to. For phase errors & ,,, smaller than 180", the IRF can be approximated to
0
{
7r
. B . [(t - to) - k . ( 2 . N,+l)]
1,
. (7)
The IRF consists of a main response ( k = 0) with high geometric resolution superposed by an infinite summation of paired echoes, whose amplitudes are given by the Bessel function J k of the first kind and kth order. From this equation, the loss, peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) can be calculated as a function of the phase error. For qbemaz = go', a PSLR of -9.6 dB and a ISLR of -4.6 dB are achieved, both values not sufficient for high image quality. Decreasing the phase error to 3" will almost eliminate the paired echoes but the number of subapertures will be increased approximately by a factor of 6, according to the following relation:
where fix(.) is the integer value operator. Due to the great number of subapertures for a requirement of qh,, = 3", the processing using this scheme becomes inefficient. In the following, a more effective way will be derived for improving the PSLR and ISLR. Starting from (1) and using the linear approximation for the phase history (6), a modified expression of (4) is obtained [14]: If(t)l =
'
I4t)l.
IR(t)l is a periodic sin(:c)/z-function and is given by:
The following properties can be associated with lR(t)l :
(9)
200
400 600 TIME ( a r b i t r a r y scale)
800
1000
Fig. 2. IRF and product terms for subaperture processing with periodical phase error. I.i(t)l is the low resolution response of the subapertures, IR(t)l is the periodical & i ( . i , ) / . r -function and I f ( t ) l is the high resolution impulse response function.
+
1. Period of ( 2 . NI l ) / B . 2. High resolution in time (zl/B). 3 . Independence of the phase error d,,,,. The second term of (9) is given by
IA(t)l = sinc
I
i
! ul I
"
linear approximation
~
d m .{[C(tl)+ C(t2)I2+ [S(tl)+ S ( ~ P2 })11]2
(11) where S(E)and C(t) are the complex Fresnel integrals. Some important properties of IA(t)l are: 1. Low resolution in time (= TI). 2. Independence of the number of subapertures. 3. PSLR and ISLR are related with 4,,,. Fig. 2 shows the IRF for YO" phase error as well as the product terms IR(t)l and IA(t)l.As can be seen from this figure, the amplitude of the paired echoes is given by the amplitude of the minima of the low resolution function IA(t)l.Two solutions were found to be efficient for decreasing the paired echoes. Improvement of the shape of the low resolution function IA(t)l. This can be achieved by using an amplitude weighting in the correlation of each subaperture, so that the effective amount of phase error is reduced and the amplitude of the minima of IA(t)l is attenuated. By using a strong amplitude weighting, the geometric resolution of IA(t)l is deteriorated by M 30% and the distance of the minima of IA(t)l apart from the mainlobe is increased proportionally. Use of a frequency overlap between the subapertures, so that the period of IR(t)l is increased and the position of the first paired echo around the mainlobe coincides with the first minimum of IA(t)l (obtained by amplitude weighting). Using a triangular weighting for the correlation in the subapertures and allowing a phase error 4,,, in each subaperture of 28.8"', the characteristics of the IRF are optimized if 'The phase error in each subaperture should be limited to 28.8' in order to maintain processing loss lower than 0.6 dB. The triangular weighting for the correlation in each subaperture will be adopted for the RTS approach, since it can he easily implemented by two moving average operations (correlation of two retangular functions). The adoption of the triangular weighting and the condition of 0, ,. 5 28.8' imposes an overlap of the subapertures that must be greater than 21% for sufficient attenuation of the sidelobes.
717
MOREIRA: REAL-TIME SAR PROCESSING "
0
"
"
'
I
-20
0
Fig. 3.
"
PiLR
=
ISLR
=
Loss
34 dB 14.4 dB
= 0.6 d B
1000 TIME (arbitrary scale)
500
1500
IRF of the RTS-algorithm for 21% frequency overlap (the response of one subaperture is represented by the dotted line).
the subapertures are overlapped by 21%, as shown in Fig. 3. A n additional weighting w(n) (Hamming, CY = 0.6) was introduced in the coherent summation of the subapertures in order to decrease the sidelobes of the IRF as much as the amplitude of the paired echoes. If better values of PSLR and ISLR are required, an overlap of 57% will place the first paired echo in the second amplitude minimum of IA(t)l and the PSLR will be attenuated to -45 dB for a Hamming weighting with a = 0.54.
IV. THERTS ALGORITHM Assuming a linear approximation of the phase history, a triangular function for correlation of the subapertures and a frequency overlap ovl, we obtain:
If ( t ) /=
(12) where iVlobl = N1/(1 - ovl) and n' = n . (1 - owl). The double integral in (12) corresponds to a double correlation of the signal s ( t ) with a rectangular function, where a frequency shift (up/down conversion of the subapertures) is performed before and after the correlation. We define the first frequency shift ( g l l ) and the second frequency shift (h,) before coherent summation of the subapertures as follows:
g n ( t ) = exp [ - j . 2 . 12'. TI . b . (t - t o ) ] .
(13)
Fig. 4. Basic block diagram of the RTS algorithm. In this figure, five subapertures are added coherently, whereby the high resolution IRF is obtained (no multilook processing is used).
Fig. 4 shows the basic block diagram of the impulse compression using the RTS algorithm. First, the received signal s ( t ) is multiplied with the frequency shifts g, , which shift the signal s ( t ) from the desired band for each subaperture up or down to baseband. The results of the frequency shift are correlated twice with a rectangular function by means of two moving average operations, where only two real additions and two real subtractions are needed per correlated point. This operation can also be interpreted as an integration with duration TI, as represented in Fig. 4. Before the second frequency shift is carried out, a time shift ri' . Tl is performed with each subaperture. This shift compensates for the time displacement of the subapertures, which are derived from different parts of the full aperture. Then the subapertures are up and down converted to the initial frequency value by means of hlL.An additional phase term b . a'* . T: is provided by h, (see 14) in order to ensure a phase continuity over the subapertures. Finally, the subapertures are added coherently, whereby the high resolution IRF is obtained. When designing an RTS processor, optimized values should be used for the maximum phase error (4e,naz5 28.8') and for the minimum overlap (ovl 2 0.21). Obviously, these requirements can be met for range processing and also for azimuth processing over all distances (including a possible squint angle). Then the integration time TI of each subaperture and the total number of subapertures ( 2 . Nloz,l 1) are calculated as follows:
+
h,(t) = ~ ( n. e)z p [ j . b . TL'. ( 2 . TI . ( t - t o )- 72' . T:)]. (14) Using (12)-(14) given by
the final equation of the RTS-algorithm is where B is the desired signal bandwidth to be processed. The processing itself is performed according to (15). For instance, the time integration of each subaperture is done over , a discrete number of points, i.e. fix[Tl . f s / ( N u n d e r ) ]where
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 30, NO. 4, JULY 1992
718
fs is the sampling frequency. The discrete integration does not influence the processing performance since the integer number of points is greater than 15 for all available SAR systems. Although the RTS algorithm works in the time domain, the computational requirement is low, since the time correlation in each subaperture is substituted by a moving average operation. The determination of the reference function for processing is simplified to the calculation of the functions gT1, h,, and the parameters TI and n’ . T I . The computational requirement can be reduced further, if the frequency shifts gn and h, are set to a submultiple of the sampling frequency. Then, short FFT’s can be used for the formation and the summation of the subapertures [13]. For a time-bandwidth product greater than 100 (more than % 10 subapertures are needed), the use of FFT’s in the RTSalgorithm is preferred. In the following, three specific algorithms for high resolution SAR processing will be discussed: multilook, autofocus and Doppler centroid estimation. We show that these algorithms can be incorporated into the RTS algorithm without additional operations, such as FFT.
enough accuracy to focus the SAR image [lo]. In the case of subaperture processing, a cross-correlation between adjacent subapertures is carried out to determine the relative time shift of the subapertures. If no velocity variations occur, the nominal time shift 71’ . TI will cause the proper positioning of the subapertures. The measured displacement Ax between the subapertures is then used to estimate the correct Doppler rate IC, = for the azimuth processing as follows
4
IC, =
+ v,,,
. Ax ’
(1 - owl) . TI
where V,,, and ICatir are the erroneous velocity and Doppler rate, respectively, assumed for processing. In comparison with the traditional auto-focus with multilook [lo], a better sensitivity to high frequency velocity variations is achieved by the subaperture processing due to the low integration time of the subapertures. Obviously, the auto-focus output degrades proportionally to the integration time, since the crosscorrelation result broadens as the integration time decreases. In an airborne SAR system, the auto-focus technique can be used to determine the aircraft velocity variations in realA. Multilook Processing For RTS processing with Lo looks, a number of (2.N1,,1+ time and to update the PRF (pulse repetition frequency) in a l ) / L o subapertures should be added coherently in order to continuous manner. This ensures a well focused image with obtain each look response. Then, the results of these additions geometric fidelity, since the azimuth pixel spacing is held are detected and added incoherently to improve the radiometric constant by changing the PRF. A prediction algorithm for the resolution. A n optimization of the look overlap to achieve an determination of the PRF (possibly based on a Kalman filter) improved radiometric resolution can be found in 1131. Only is required due to the delay associated with the autofocus the last stage of the RTS processor (coherent and incoherent computation. addition of the subapertures) must be changed to vary the number of looks and the look overlap. C. Doppler Centroid Estimation An IML (improved multilook) technique [13] can also Doppler centroid estimation can be carried out directly in be implemented efficiently in the RTS processor. The IML azimuth processing because the subaperture domain (index technique is based on the formation of looks with different bandwidths. The looks with larger bandwidth contribute to an 7~’)corresponds to a frequency domain with a low frequency improvement of the overall geometric resolution, while the resolution = P R F / ( 2 . NIo7./ 1). Then, by averaging looks with smaller bandwidth improve the overall radiomet- the signal intensity in each subaperture over the azimuth ric resolution. Multilook processing with the IML technique time or over a group of range bins, an estimation of the improves the radiometric resolution of SAR images, since it azimuth antenna pattern can be made [14]. It is assumed makes a more efficient use of the available signal bandwidth. that the number of subapertures is large enough (> 10) For processing with the IML technique using the correlation that the azimuth antenna pattern can be considered constant in frequency domain (e.g., 1111, [24]), approximately 50% within a subaperture. After interpolating the results of the more computation time is required than with the traditional averaging in order to obtain the desired frequency accuracy for multilook technique. This occurs, since a number of L L looks determination of the Doppler centroid, conventional methods with larger bandwidth and a number of Ls looks with smaller (e.g., energy balance method [lo]) can be used for efficient bandwidths must be formed in the frequency domain, followed estimation of the Doppler centroid. The time domain structure of the RTS algorithm (no block by ( L L L s ) inverse FFT’s. Then L.5 extra inverse FFT’s processing is used) imposes no restrictions on the update rate are needed, leading to the increased computation time. Using of the Doppler centroid. This feature is an important point for the RTS algorithm, the extra computation time is reduced to a minimum. Only extra additions are needed to create the looks the accurate processing of future spaceborne SAR systems. with larger and smaller bandwidth (coherent addition) and to For the X-SAR system, which will be launched on the Space Shuttle in 1993, an azimuth Doppler centroid variation of up form the final image (incoherent addition). to 280 Hz/s will occur due to the instability of the platform [4]. Since the required accuracy of the Doppler centroid is f 5 0 Hz, B. Auto-Focus an update rate of 0.2 s is necessary to obtain a high quality Auto-focus is commonly implemented in azimuth process- image. The RTS algorithm can accommodate these updates ing in order to determine the Doppler frequency rate with without decreasing its computational efficiency.
+
+
MOREIRA: REAL-TIME SAR PROCESSING
V. EXTENSION TO 2-D
719
TABLE I PROCESSING PARAMETERS FOR THE SIMULATION OF THE IRF IN C-BAND(SQUINT ANGLE= 8", AIRCRAFT VELOCITY = 70 mis AND SHORTEST RANGETO POINTTARGET= 3000 m)
PROCESSING
The extension of the 1-D RTS algorithm to a 2-D RTS algorithm means that the functions g, and h, are twodimensional and that the triangular function for the correlation should be substituted by a pyramidal function. However, this 2-D implementation will have limitations on the correction of the linear and quadratic range migration and will require a large memory for the correlation process. Therefore, the 1-D RTS algorithm (as described in the previous section), with a time domain correction of the range migration, will be used to accomplish the 2-D nature of the azimuth reference function. The adopted processing steps are the same as for a 1-D time domain SAR processor (range compression, corner turn operation, azimuth processing, and detection). The correction of range migration is divided into two steps: A) range walk correction, which is performed during range processing and corner turn operation and B) range curvature correction, which is carried out by an interpolation of samples prior to the coherent summation of the subapertures in the azimuth processing.
~~
Simulation Parameters
range orocessine
azimuth processing
processed signal bandwidth
100 MHz
100 Hz
_____
reference function duration
5
1.7 s 0.104 s
PS
integration time of each subaperture
0.178 p s
subaperture overlap
57%
number of subapertures
63
21% 20
"I
A. Range Walk Correction
The fractional part of the range walk can be easily corrected during range processing by adding a linear phase term before the coherent summation of subapertures is carried out. Considering the periodic sin(z)/z-function R ( t )in (9) and including the linear phase term, 41 = 2 . n' . (1 - out) . TI . b . At, where At is the desired time shift, we obtain: sinj7r . B . ( t - At - to)] IR(t)l = sin ?r-B.(t-At-t,)]
[
Fig. 5. Simulated point target response for the E-SAR system of DLR (C-band, 8' squint angle).
VI.
SIMULATION OF THE
IRF
(19)
By using the parameters for the C-band mode of the airborne
Since the IRF is given by IR(t)l . IA(t)l, where IA(t)l is a slowly varying function (refer to (9) and (ll)), the obtained time shift of the range compressed data is At. This first step of the range migration correction ensures an exact interpolation of the data in range direction which can be updated from sample to sample, since no block processing is used by the RTS algorithm. The update of the range walk becomes necessary when the Doppler centroid varies with range (e.g., large swath width and/or small look angle). The integer part of the range walk is corrected by a shift of range samples during the corner turn operation. An additional deskew of the data is performed with the final image (possible during slant to ground range conversion) in order to compensate for the integer sample shift after the range processing.
integrated, a point target simulation and processing was carried out in order to estimate the performance of the RTS algorithm. The critical points in processing the E-SAR data are the correction of the motion errors due to the instability of the small DO-228 aircraft and the high range walk that occurs when flying with strong crosswinds. The motion errors of the E-SAR system can be corrected efficiently by the reflectivity displacement method [15] and will not be covered here. We analyze below the performance of the RTS algorithm for high drift angles. A detailed hardware simulation of the real-time processor including quantization and linearity errors due to the finite bit length in the processing can be found in [12]. A squint angle of 8' was used for the simulation. In this case, the total range walk within the azimuth integration time is greater than 30 range bins, causing an azimuth time-bandwidth product
