IRCI‐free multiple subband MIMO SAR

Mohammed AlShaya, Institute for Digital Communications, University of Edinburgh, Edinburgh, UK. Email: mohammed.alshaya@ed.ac.uk Abstract An inter‐range cell interference (IRCI)‐free multiple subband multiple‐input multiple‐ output (MIMO) synthetic aperture radar (SAR) algorithm is proposed to obtain a high‐resolution wide swath (HRWS) imaging. The received signals of the proposed configuration are formulated as multiple‐input single‐output (MISO) system identification problems using the principle of the displaced phase centre in a way that all the subband waveforms are processed simultaneously without the need to separate them at the receiver. This would maximise the bandwidth utilisation efficiency. The channel impulse response in the range dimension is identified using a frequency domain system identification‐‐based estimation algorithm, which is free from IRCI, instead of using conventional matched filters. A pulse repetition frequency lower than the Doppler bandwidth is used to obtain the HRWS imaging, and the resulted azimuth ambiguity is removed using a set of spatial filters applied on the received signals formulated as separate MISO system identification problems. Finally, both simulated and constructed raw data are used to validate the efficiency of the proposed algorithm.


| INTRODUCTION
The attention of researchers has been drawn recently to the application of a multiple-input multiple-output (MIMO) concept to synthetic aperture radar (SAR) imaging due to its powerful applications in both military and civilian sectors [1]. Two different configurations, namely, single-input multipleoutput (SIMO) SAR [2] and MIMO SAR [3][4][5] are introduced in the literature to overcome the contradicting requirements between obtaining the desired wide swath and high cross-range resolution in conventional single-input single-output (SISO) stripmap SAR. A SIMO SAR with a single-phase centre and multiple azimuth beam is introduced in [6,7] in which multiple contiguous receiving beams with different squint angles are used. The Doppler spectra of all narrow receiving beams are then combined to form a wide Doppler spectrum.
It is stated in [8], where the concept of MIMO SAR is introduced, extra phase centres can be obtained, which allows to form a high-resolution wide swath (HRWS) image. MIMO SAR can be categorised into orthogonal waveform encoding MIMO SAR [9], beam-space MIMO SAR [10] and multiple subband MIMO SAR [11], which is considered here.
Frequency division multiple access waveforms are used in multiple subband MIMO SAR in which each subband waveform is emitted using a different carrier frequency. According to the literature [12,13], returns that correspond to each subband waveform are separated at the receiver using a bank of bandpass filters. The leakage among the subband spectra (i.e. overlapping) degrades the performance of the estimated range profile. Therefore, guard bands are added between the spectra of the adjacent subbands to reduce the overlapping, but this affects the bandwidth utilisation efficiency. Multiple subband MIMO SAR is used to improve the range resolution because it is possible to transmit more bandwidth while the narrowband assumption corresponds to the bandwidth of the individual transmitted subband waveform. It is worth pointing out that the number of effective phase centres (EPCs) in multiple subband MIMO SAR is the same as the one obtained in SIMO SAR [11].
In [14], we first proposed a multiple subband MIMO SAR configuration shown in Figure 1. The beams of the transmitters are wide and multiple adjacents received subbeams with different phase centres and squint angles are generated in the azimuth direction. The system pulse repetition frequency (PRF) only needs to satisfy the Nyquist rate of the partial Doppler bandwidth corresponding to each receiving subbeam. Wide Doppler bandwidth is then synthesised to obtain wide swath high-resolution images. This configuration has a drawback that the receive beams are contiguous, as shown in Figure 2, which makes the return from the covered area of one receive 3 dB beam badly contaminated by the returns from the covered area by the other receive 3 dB beams. This results in a bad azimuth ambiguity which will degrade the quality of the formed SAR image.
It should be noted that the range profile in SAR consists of a large number of scattering points rather than a small number as in the case of the classic radar mode. Accordingly, it is desirable to have a high range resolution with low sidelobe levels in the estimated SAR range profile so that the scattering points with low reflectivity are not overwhelmed with bright scattering points. The effects of the sidelobes of a scattering point on the adjacent range cells are known as inter-range cell interference (IRCI), which increases proportionally with the number of range cells and reduces the resultant range resolution accordingly [15]. Cyclic prefix-based orthogonal frequency division multiplexing waveforms are proposed in the literature [15,16], which has IRCI-free property . This class of waveforms is impractical in some applications, such as stripmap SAR, as the waveform length should be at least equal to the channel to be estimated.
Here, we extend the study presented in [14] in a way that the beamwidths of all transmitters and receivers are the same and they simultaneously illuminate the same imaging area as shown in Figure 3 (i.e. the beams of all transmitters and receivers are the same as the Tx beam shown in Figure 2), which allows to overcome the receiving interbeams overlapping problem. In addition, the received signal in each receiver is formulated as a multiple-input single-output system identification problem using the principle of displaced phase centre (DPC), which is then used to estimate impulse response in the range dimension using our proposed frequency domain system identification (FDSI) algorithm which has the IRCI-free property; although the transmitted waveforms are partially correlated. In addition, the proposed configuration does not require the separation of the subbands at the receiver as they are processed jointly without a need to add guard bands between the adjacent subbands. Also, utilising the whole bandwidth, linear frequency modulated (LFM) waveforms are used for transmission, which have many desirable properties such as the constant envelop and unity peak to average power. Finally, a wide swath image is formed by utilising a PRF lower than the Doppler bandwidth, and the resulted azimuth ambiguity is removed using a set of spatial filters.
Notations. Matrices are denoted by upper case bold letters (e.g. A) while scalars and vectors are denoted by upper∖lower case letters (e.g. a and A) and bold lower case letters (e.g. a), respectively. In addition, (.) T and (.) H denote the transpose and hermitian operators, respectively. diag(a) denotes a matrix with the vector a as main diagonal. C represents the field of complex numbers.

| MIMO SAR SYSTEM MODEL
In this section, the multiple subband MIMO SAR system is assumed to be narrowband and the numbers of antenna elements at the transmitter and receiver are denoted by M and N, F I G U R E 1 MIMO SAR configuration proposed in [14]. MIMO, multiple-input multiple-output F I G U R E 2 Transmit and receive antenna beam patterns which are generated using a raised cosine function F I G U R E 3 Multiple subband MIMO SAR configurations in which the beam patterns of the transmitters and receivers are the same. MIMO, multiple-input multiple-output; SAR, synthetic aperture radar ALSHAYA ET AL. -687 respectively. In addition, the transmitted waveforms are assumed to be subband waveforms (i.e. LFM waveforms) which can be expressed as the following:

| Received signal model
The baseband received signal at the nth receiving channel can be written as the following: F c is the carrier frequency, t and η denote the fast time and slow time, respectively. x m (t) denotes the mth subband waveform, σ l is the radar cross-section of the lth scattering point, w n (t) is the additive white noise generated by the nth receiver and R l,i denotes the slant range from the ith transmitter/receiver to the lth scattering point and can be approximated under the assumption of a straight sensor trajectory as the following [17] (∀n = 1, 2, …, N ): where v p is the platform velocity and R l denotes the slant range of the closest approach to the lth scatterer. The azimuth impulse response is a function of the slant ranges and can be expressed as the following for the case when the signal is transmitted by the mth transmitter and received by the nth receiver [2]: where Δy mn denotes the separation between the mth transmitter and the nth receiver. The following can be expressed using Equations (5) and (6) where h az,t1 and h az,t2 denote the azimuth impulse response whose phase centres' locations are as indicated by y t1 and y t2 in Figure 4, respectively. The following can be expressed as well using the principle of DPC for the general case when N ≥ 3 (i.e. without loss of generality it is assumed M = N = 3 as shown in Figure 3): The received signal in Equation (2) after sampling in both range and azimuth dimensions by the sampling frequencies of f s = 1/T s and f p = PRF, respectively, can be expressed as the following (∀n = 1, 2, …, N ): where h mn ½l; n a � ¼ σ l e −j2πF c τ l;mn ðn a Þ ¼ σ l e jϕ l;mn ðn a Þ ð17Þ n t and n a denote the range and azimuth time indices, respectively.

F I G U R E 4
The locations of the phase centres for the case when M = N = 2. y t1 and y t2 indicate the locations of the phase centres whose azimuth impulse responses are denoted by h az,t1 and h az,t2 , respectively 688 -ALSHAYA ET AL.

| Azimuth ambiguity removal
The received signals are aliased in the azimuth dimension because the PRF used is assumed to be lower than the Doppler bandwidth. The resulted azimuth ambiguity will be removed using the algorithm proposed in [2], which states that it is possible to reconstruct the original signal from N independent representations of its aliased version for the case when the original signal is sampled by 1/N its Nyquist frequency. A special case when N = 2 is first assumed and then will be generalised to the case when N ≥ 3.
The Fourier transform of Equation (16) across both azimuth and range dimension for the case when M = N = 2 can be expressed as the following: The received signals in Equations (18) and (19) can be expressed as a function of the impulse responses h az,t1 and h az,t2 using the principle of DPC Equations (7)-(10) as the following: where H t1 and H t2 denote the frequency channel impulse response whose EPCs are located in y t1 and y t2 , respectively, as shown in  (20) and (21) can be expressed as a function of H 12 = H 21 which corresponds to the impulse response of a SISO case in which the EPC is located midway between the two receivers as the following (noise is not included in the following equations for the sake of clarity), The azimuth spectrum of the received signal aliased as the Doppler bandwidth is assumed to be NPRF (i.e. 2PRF) as illustrated in Figure 5 which can be expressed as the following (∀n = 1, 2) [18]: e Z 1 ð f r ; f a Þ and e Z 2 ð f r ; f a Þ can be combined as explained in [18] such that the original spectrum U( f r , f a ) is recovered as expressed in the following equations for the interval I 1 which covers the first half of the Doppler spectrum: The spatial filter on the jth Doppler frequency interval when the signal is received by the nth receiver is denoted by P nj . The following can be expressed for interval I 2 (i.e. second half of the Doppler spectrum) after shifting to I 1 to allow for setting up of the linear systems of Equation (30).

F I G U R E 5
The spectrum at the nth receiver when the signal is transmitted by the mth transmitter ALSHAYA ET AL.
The reconstruction filters can be computed as the following: where Pð f a Þ ¼

| N ≥ 3
Without loss of generality, it is assumed when N = M = 3 other cases can be easily derived using the method which will be presented in the following. The Fourier transform of Equation (16) across both azimuth and range dimensions can be expressed as the following: The received signal in Equations (33)-(35) can be expressed as a function of the impulse responses h 21 , h 22 and h 23 , respectively, using the principle of DPC Equations (11)-(15) as the following: where As the constant phase terms in Equations (39)-(41) are relatively small as they are functions of the slant range R l , the following approximation can be used in order to make all 1, 2, 3): where The spectra of the transmitting waveforms X m ( f r ), ∀m ∈ [1-3] are chosen such that B( f r , f a ) occupies the whole fast time spectrum to identify the impulse response in the range dimension. Next, the received signals can be expressed as a function of H 22 which corresponds to the impulse response of a SISO case in which the EPC is located at the second receiver as the following(∀n = 1, 2, 3): where As described previously, it is expected to have aliasing in the azimuth dimension as the Doppler bandwidth of the received signals is assumed to be 3PRF while the sampling frequency in the azimuth dimension is f p = PRF. The ambiguity can be removed to perfectly reconstruct U( f r , f a ) in the same way as the one described previously for the case when M = N = 2 (i.e. starting from Equation (25)). The reconstruction filters in this case can be computed as the following: where Qð f a Þ ¼ Pð f a Þ ¼ It should be noted that the centre frequency of each Doppler interval (i.e. I 1 , I 2 and I 3 ) can be computed by the following Equation [19]: I 1 , I 2 and I 3 are then concatenated to obtain the original spectrum U( f r , f a ).

| Impulse response estimation and image formation
The image formation algorithm used is chirp scaling (CS) [20], for the sake of computational simplicity. Throughout  [14] as explained next. The recovered signal U( f r , f a ) is multiplied in the range-Doppler domain with the following phase function to adjust the all range migration to a reference range r ref : βð f a Þ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a( f a ) denotes the CS factor, k r is the modulation rate of the transmitted chirps and the azimuth frequency f a varies within the following range.
The impulse response in the range dimension is then identified using the proposed FDSI-based algorithm at every Doppler frequency [14].
The next step is to estimate the range impulse response at every Doppler frequency using our proposed FDSI algorithm in [14]. The signal in Equation (53) can be rewritten as the following: where K is the length of the waveform transmitted, w s [n t , f a ] is the noise after being multiplied with H CS,1 in the range-Doppler domain and b s [n t , f a ] is the combined full bandwidth transmitted LFM waveform with a chirp rate that varies with the Doppler frequency as given below The scaled signal can be expressed as a function of the Doppler frequency ( f a ) as the following: -691 The channel impulse response h( f a ) in Equation (60) is zero padded by (K − 1) in order to make the matrix B s ( f a ) circulant as the following: where B s,C ( f a ) is as defined in [21]. The channel impulse response can be estimated as the following: The discrete fourier transform (DFT) structure of B s,C ( f a ) expressed below can be used to estimate the impulse response in the frequency domain where F is the DFT matrix, X is the Fourier transform of the first row of the matrix B s,C ( f a ). It should be pointed out that only the first L elements of the estimated padded impulse response are kept while the other elements are forced to be zero. The impulse response estimation described above should be performed for all Doppler frequencies.
The estimated impulse response is then multiplied by the following phase function in the range frequency-Doppler domain to perform the bulk range cell migration correction: Residual phase compensation and azimuth compression are performed by multiplying the data in the range time-Doppler domain with the following: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi f Δφð The image is finally formed by calculating the inverse azimuth Fourier transform.

| Performance of range profile reconstruction
The number of transmitters/receivers separated by d = 1.5 m are assumed to be M = N = 3, and the spectra of the transmit waveforms are shown in Figure 6 in which the total transmit bandwidth is 100 MHz. The point spread function of our proposed FDSI-based algorithm is compared with the point spread function of the matched filter, as shown in Figure 7. The waveform used in the matched filter is an LFM waveform with a bandwidth of 100 MHz. It is clear from Figure 7 that the point spread function of our proposed algorithm exhibits the IRCI-free property, although the spectra of the transmitted waveforms are partially overlapped. The IRCI-free property is possible because our proposed estimation algorithm takes into account the fact that the spectrum of the combined signal is non-flat which minimises the sidelobes level of the estimated range profile, unlike the case of the matched filter (i.e. which assumes a flat spectrum). It should be noted that the point spread function shown in Figure 7 is generated under the assumption that the scatterer falls at the exact sample time in that case, our proposed algorithm generates ideally zero level F I G U R E 6 The spectra of the transmitted waveforms F I G U R E 7 The normalised point spread function of the proposed algorithm and matched filter 692sidelobes, unlike the case of the matched filter. This ideal scenario is presented in the article to demonstrate the IRCIfree property in a way similar to the IRCI-free approaches in the literature [15,16]. The point spread function of the proposed MIMO FDSI-based algorithm and matched filter is shown in Figure 8 for different signal to noise ratios, where it is clear that the matched filter case is more robust to white noise.
Next, consider the case in which there are a number of scatterers located at different ranges. The estimated impulse response is shown in Figure 9 where one can see that the sidelobes' level is not affected (i.e. IRCI-free property is maintained).

| Azimuth ambiguity removal
Consider the simulation parameters listed in Table 1 and assume that there is a single scatterer located at the swath centre.
The transmitted waveforms are assumed to be subband LFM. The azimuth and range cuts for the noiseless case are shown in Figures 10 and 11, respectively. One can see that no azimuth ambiguity is present for the case of MIMO SAR, unlike the case of the conventional SISO SAR.

| Azimuth ambiguity analysis
The azimuth ambiguity due to the interbeams overlapping in the configuration shown in Figure 1 and proposed in [14] is compared with the proposed system configuration shown in Figure 3 in which the transmitters and receivers use the same  -693 wide antenna beams. The azimuth ambiguity to signal ratio (AASR) defined in [21] will be used for the comparison.
The antenna beam patterns used in the simulation are as shown in Figure 2. The proposed configuration uses only the wide Tx beam shown in Figure 2 for both transmission and reception. It should be noted that all the side beams of the configuration shown in Figure 1 have the same azimuth ambiguity characteristics. Therefore, only one of the side beams is considered. The AASR as a function of the PRF is shown in Figure 12, where it is clear from the AASR curves that the azimuth characteristics of the system configuration with the narrow receiving beams are worse because the sidelobes of the receiving beams shown in Figure 2 are within the mainlobe of the transmit beam. In addition, the azimuth ambiguity of the system configuration with narrow receiving beams is a function of azimuth spacing among the receive beams in addition to the PRF. Therefore, the AASR cannot be reduced unless the azimuth spacing among the receive beams is increased. It is worth pointing out that digital beamforming in elevation on receive is proposed in [21] to mitigate the effect of interbeams overlapping, but this will be at the expense of the system complexity.
It is assumed in the AASR curve of the proposed configuration that the azimuth ambiguity resulted from using a sampling frequency lower than the Doppler bandwidth is removed using a set of spatial filters, as explained previously.

| Raw data simulation
The multiple subband MIMO SAR data is constructed using airborne SISO SAR data 'Gotcha' [22] in order to demonstrate the performance of our proposed algorithm in a cluttered environment. The formed image of the original SISO SAR data is shown in Figure 13 which is going to be used later for comparison.  6. The subbands are superimposed and then down-sampled by a factor of N in the azimuth dimension to construct the multiple subbands MIMO SAR data at the nth receiver.
Consider the main system parameters listed in Table 2. The formed image using a single received channel is shown in Figure 14, where it is clear that the image is aliased. The formed image using the two receiving channels is shown in Figure 15, where it is obvious that the azimuth ambiguity is removed. The formed image using the proposed algorithm is compared with the original SISO SAR data by showing the range and cross-range cuts of the corner reflector indicated by the red rectangular in Figures 13 and 15. The range and crossrange cuts are shown in Figures 16 and 17, respectively. It is clear that the quality of the formed image using the proposed algorithm is not degraded as compared with the original SISO SAR data.

| CONCLUSION
An IRCI-free multiple subband MIMO SAR algorithm to obtain HRWS imaging is presented. The proposed algorithm utilises the available bandwidth to the maximum efficiency as it does not require separating the subbands at the receiver (i.e. no need to add guard bands between adjacent subbands). It is IRCI-free; although the spectra of transmitting waveforms are partially overlapped which facilitates the use of conventional LFM waveforms, which have many desirable properties such as constant envelope and unity peak to average power. The performance of azimuth ambiguity characteristics of the proposed configuration is analysed and compared with the configuration with the narrow receiving beams. Moreover, a PRF lower than the Doppler bandwidth is used to obtain a wide swath and the resultant azimuth ambiguity is removed using a set of spatial