# 3D Imaging Algorithm for Down-Looking MIMO Array SAR Based on Bayesian Compressive Sensing

Authors: Xiaozhen Ren; Lina Chen; Jing Yang
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

Down-looking MIMO array SAR can reconstruct 3D images of the observed area in the inferior of the platform of the SAR and has wide application prospects. In this paper, a new strategy based on Bayesian compressive sensing theory is proposed for down-looking MIMO array SAR imaging, which transforms the cross-track imaging process of down-looking MIMO array SAR into the problem of sparse signal reconstruction from noisy measurements. Due to account for additive noise encountered in the measurement process, high quality image can be achieved. Simulation results indicate that the proposed method can provide better resolution and lower sidelobes compared to the conventional method.


---

## Body

## 1. Introduction

Traditional synthetic aperture radar (SAR) is a microwave sensor which can reconstruct two-dimensional (2D) images of the observed area with weather independence and all-day operation capabilities [1, 2]. However, traditional 2D SAR works in side-looking mode and often meets with shading and lay over effects in urban and mountain areas. Compared with 2D SAR, 3D SAR has distinct advantage in estimation of forest height, 3D digital maps, complex terrain mapping, and so on. Multibaseline SAR tomography is an advanced 3D SAR imaging mode, which forms an additional synthetic aperture in the height direction. Therefore, it has resolving capability along this dimension. Unfortunately, for the current SAR tomography, it is almost impossible to avoid an uneven track distribution in repeat-pass data acquisition, which is just the main reason for the strong ambiguity in height [3, 4].Down-looking array SAR is an innovative imaging mode, which obtains range resolution by pulse compression, azimuth resolution by virtual aperture synthesis with platform movement, and cross-track resolution by a linear array antenna [5–7]. Down-looking array SAR can overcome restrictions of shading and lay over effects in side-looking SAR and also avoid the height ambiguity problem in SAR tomography caused by the uneven track distribution. However, in order to avoid the grating lobe effect in the cross-track direction, a large number of antenna elements are required, which increase the cost and complexity of the equipment. Therefore, multiple-input-multiple-output (MIMO) antenna array is often used in down-looking array SAR, which can reduce the number of real antenna elements largely for a given size of antenna array [8–10]. The common methods used for down-looking MIMO array SAR imaging are usually based on matched filter, which often suffer from low resolution and high sidelobe interference in the images [7, 8]. Moreover, the cross-track resolution is limited by the length of linear array. Hence, high resolution imaging algorithms are desired.In recent years, Bayesian compressive sensing (BCS) has caused widespread concern, showing significant advantages to sparse signal reconstruction [11]. BCS methods provide certain improvements compared with norm-based CS methods in low noise level, by exploiting the sparseness prior distribution of the image scene. In addition, the Bayesian framework takes into account the additive noise encountered when implementing compressed sampling. Therefore, some BCS based methods for SAR applications have been concerned about recently [12, 13]. As the 3D illuminated scene contains only a very small strong scattering centers compared with the total space cells, representing strong spatial sparsity in high frequency radar application, we propose a new 3D imaging algorithm for down-looking MIMO array SAR based on BCS.The rest of the paper is organized as follows. Section2 presents the geometry and principle of down-looking MIMO array SAR system. In Section 3, a new 3D imaging algorithm for down-looking MIMO array SAR is described in detail. The performance of the method is investigated in Section 4. Finally, Section 5 gives a brief conclusion.


## 2. Down-Looking MIMO Array SAR

### 2.1. Geometrical Model

The geometry of down-looking MIMO array SAR is shown in Figure1. x, y, and r denote the azimuth, cross-track, and slant range direction, respectively. The radar platform flies along the x-axis corresponding to azimuth direction, with velocity v at height H. The thinned linear antenna array, which contains M transmitting antenna elements and N receiving antenna elements, is mounted in the cross-track direction along the wings. The transmitting antenna elements are located at the tips with the distance d, the receiving antenna elements are centered at the y-axis with spacing M
d
/
2, and the distance between the transmitting and receiving antenna elements is d
/
2. The thinned linear antenna array works in the time division mode. Each time, only one transmitting antenna element transmits signal and all the receiving antenna elements receive echo simultaneously. The transmitting antenna elements work sequentially and an aperture synthesis period is acquired until all the transmitting antenna elements have worked once. According to the principle of equivalent phase center, the thinned linear array formed by the above positions and work mode is equal to a fully distributed virtual uniform linear array [5, 10]. The virtual antenna array is composed of M
N virtual elements and works in self-transmitting and receiving mode. That is to say, each virtual antenna element transmits and receives signal by itself. These virtual antenna elements are uniform distributed along the wings and centered at the y-axis. Each virtual antenna element is located at the mean position of a real single transmitting element and a real single receiving element, and the distance between individual virtual antenna elements is d
/
2. Figure 2 shows the transmitting and the receiving order of each antenna element for a down-looking MIMO array SAR.Figure 1
Geometry of down-looking MIMO array SAR.Figure 2
Transmitting and receiving order of down-looking MIMO array SAR.


### 2.2. Equivalent Phase Error Compensation

Consider the data acquisition shown in Figure1. At the slow time t
m, the position of the ith transmitting antenna element is given by (x, y
T
i, H), where x
=
v
t
m is the azimuth position, and y
T
i is the cross-track position of the ith transmitting antenna element. The position of the jth receiving antenna element is given by (x, y
R
j, H), where y
R
j is the cross-track position of the jth receiving antenna element. For a point scatterer P positioned at (x
p, y
p, z
p), the transmitting and receiving paths R
T
i and R
R
j are given by
(1)
R
T
i
=
(
x
-
x
p
)
2
+
(
y
T
i
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
T
i
-
y
p
)
2
+
R
P
2
,
R
R
j
=
(
x
-
x
p
)
2
+
(
y
R
j
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
R
j
-
y
p
)
2
+
R
p
2
,
where
(2)
R
p
=
(
x
-
x
p
)
2
+
(
H
-
z
p
)
2
.
Then the complete travelling path of the wave from the ith transmitting antenna element to point scatterer P to the jth receiving antenna element is given by
(3)
R
i
,
j
=
R
T
i
+
R
R
j
=
(
y
T
i
-
y
p
)
2
+
R
P
2
+
(
y
R
j
-
y
p
)
2
+
R
P
2
.
From (2) we can get that the complete travelling path contains two square roots, which will lead to complex computation for the following imaging process. According to the principle of equivalent phase center, the above complete travelling path can be equal to the dual echo paths from the virtual antenna element located at y
=
(
y
T
i
+
y
R
j
)
/
2 to point scatterer P [5]. And the equivalent echo path can be written as
(4)
R
=
(
x
-
x
p
)
2
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
-
y
p
)
2
+
R
p
2
.
Then, the phase difference between the virtual antenna element and the reality antenna element is given by
(5)
Δ
φ
=
2
π
λ
(
R
i
,
j
-
2
R
)
=
2
π
λ
(
(
y
T
i
-
y
p
)
2
+
R
P
2
+
(
y
R
j
-
y
p
)
2
+
R
P
2
m
m
m
m
-
2
(
y
-
y
p
)
2
+
R
p
2
)
≈
π
(
y
T
i
-
y
R
j
)
2
2
λ
R
p
.When the scatterers are located at the far field, the phase differenceΔ
φ is approximate to zero. Otherwise, the phase compensation should be implemented before the image processing. From (5) it can be seen that the compensated phase varies with R
p; in practice we use the center of view field as reference point to compensate the whole view field.Furthermore, the down-looking MIMO array SAR works in the time division mode and the virtual antenna elements obtained from different pulse repetition period are not in a straight line with the movement of the platform. Therefore, in order to obtain a fully distributed virtual uniform linear array, the motion compensation should be implemented. The compensated phase caused by the overtake or lag phases owing to the position difference of the antenna elements is given by(6)
Δ
φ
l
,
k
=
4
π
λ
(
(
l
-
L
2
)
1
/
2
(
(
M
v
PRF
(
l
-
L
2
)
+
Δ
x
k
-
x
p
)
2
m
m
m
m
m
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
M
v
P
R
F
(
l
-
L
2
)
)
1
/
2
m
m
m
m
-
(
(
M
v
PRF
(
l
-
L
2
)
-
x
p
)
2
m
m
m
m
m
m
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
M
v
P
R
F
(
l
-
L
2
)
)
1
/
2
)
,
m
m
m
m
m
m
m
l
=
1,2
,
…
,
L
;
k
=
1,2
,
…
,
M
,
where L is the azimuth sample number and Δ
x
k is the move spacing between the kth and the first transmitting antenna element in the azimuth direction. v is the velocity of the platform and PRF is the pulse repetition frequency.After the process above, the collected data of down-looking MIMO array SAR can be regarded as received by the fully distributed virtual linear array.


## 2.1. Geometrical Model

The geometry of down-looking MIMO array SAR is shown in Figure1. x, y, and r denote the azimuth, cross-track, and slant range direction, respectively. The radar platform flies along the x-axis corresponding to azimuth direction, with velocity v at height H. The thinned linear antenna array, which contains M transmitting antenna elements and N receiving antenna elements, is mounted in the cross-track direction along the wings. The transmitting antenna elements are located at the tips with the distance d, the receiving antenna elements are centered at the y-axis with spacing M
d
/
2, and the distance between the transmitting and receiving antenna elements is d
/
2. The thinned linear antenna array works in the time division mode. Each time, only one transmitting antenna element transmits signal and all the receiving antenna elements receive echo simultaneously. The transmitting antenna elements work sequentially and an aperture synthesis period is acquired until all the transmitting antenna elements have worked once. According to the principle of equivalent phase center, the thinned linear array formed by the above positions and work mode is equal to a fully distributed virtual uniform linear array [5, 10]. The virtual antenna array is composed of M
N virtual elements and works in self-transmitting and receiving mode. That is to say, each virtual antenna element transmits and receives signal by itself. These virtual antenna elements are uniform distributed along the wings and centered at the y-axis. Each virtual antenna element is located at the mean position of a real single transmitting element and a real single receiving element, and the distance between individual virtual antenna elements is d
/
2. Figure 2 shows the transmitting and the receiving order of each antenna element for a down-looking MIMO array SAR.Figure 1
Geometry of down-looking MIMO array SAR.Figure 2
Transmitting and receiving order of down-looking MIMO array SAR.


## 2.2. Equivalent Phase Error Compensation

Consider the data acquisition shown in Figure1. At the slow time t
m, the position of the ith transmitting antenna element is given by (x, y
T
i, H), where x
=
v
t
m is the azimuth position, and y
T
i is the cross-track position of the ith transmitting antenna element. The position of the jth receiving antenna element is given by (x, y
R
j, H), where y
R
j is the cross-track position of the jth receiving antenna element. For a point scatterer P positioned at (x
p, y
p, z
p), the transmitting and receiving paths R
T
i and R
R
j are given by
(1)
R
T
i
=
(
x
-
x
p
)
2
+
(
y
T
i
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
T
i
-
y
p
)
2
+
R
P
2
,
R
R
j
=
(
x
-
x
p
)
2
+
(
y
R
j
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
R
j
-
y
p
)
2
+
R
p
2
,
where
(2)
R
p
=
(
x
-
x
p
)
2
+
(
H
-
z
p
)
2
.
Then the complete travelling path of the wave from the ith transmitting antenna element to point scatterer P to the jth receiving antenna element is given by
(3)
R
i
,
j
=
R
T
i
+
R
R
j
=
(
y
T
i
-
y
p
)
2
+
R
P
2
+
(
y
R
j
-
y
p
)
2
+
R
P
2
.
From (2) we can get that the complete travelling path contains two square roots, which will lead to complex computation for the following imaging process. According to the principle of equivalent phase center, the above complete travelling path can be equal to the dual echo paths from the virtual antenna element located at y
=
(
y
T
i
+
y
R
j
)
/
2 to point scatterer P [5]. And the equivalent echo path can be written as
(4)
R
=
(
x
-
x
p
)
2
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
=
(
y
-
y
p
)
2
+
R
p
2
.
Then, the phase difference between the virtual antenna element and the reality antenna element is given by
(5)
Δ
φ
=
2
π
λ
(
R
i
,
j
-
2
R
)
=
2
π
λ
(
(
y
T
i
-
y
p
)
2
+
R
P
2
+
(
y
R
j
-
y
p
)
2
+
R
P
2
m
m
m
m
-
2
(
y
-
y
p
)
2
+
R
p
2
)
≈
π
(
y
T
i
-
y
R
j
)
2
2
λ
R
p
.When the scatterers are located at the far field, the phase differenceΔ
φ is approximate to zero. Otherwise, the phase compensation should be implemented before the image processing. From (5) it can be seen that the compensated phase varies with R
p; in practice we use the center of view field as reference point to compensate the whole view field.Furthermore, the down-looking MIMO array SAR works in the time division mode and the virtual antenna elements obtained from different pulse repetition period are not in a straight line with the movement of the platform. Therefore, in order to obtain a fully distributed virtual uniform linear array, the motion compensation should be implemented. The compensated phase caused by the overtake or lag phases owing to the position difference of the antenna elements is given by(6)
Δ
φ
l
,
k
=
4
π
λ
(
(
l
-
L
2
)
1
/
2
(
(
M
v
PRF
(
l
-
L
2
)
+
Δ
x
k
-
x
p
)
2
m
m
m
m
m
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
M
v
P
R
F
(
l
-
L
2
)
)
1
/
2
m
m
m
m
-
(
(
M
v
PRF
(
l
-
L
2
)
-
x
p
)
2
m
m
m
m
m
m
+
(
y
-
y
p
)
2
+
(
H
-
z
p
)
2
M
v
P
R
F
(
l
-
L
2
)
)
1
/
2
)
,
m
m
m
m
m
m
m
l
=
1,2
,
…
,
L
;
k
=
1,2
,
…
,
M
,
where L is the azimuth sample number and Δ
x
k is the move spacing between the kth and the first transmitting antenna element in the azimuth direction. v is the velocity of the platform and PRF is the pulse repetition frequency.After the process above, the collected data of down-looking MIMO array SAR can be regarded as received by the fully distributed virtual linear array.


## 3. Three-Dimensional Imaging Algorithm for Down-Looking MIMO Array SAR

Based on the principle of equivalent phase center, the thinned linear array formed by the time division mode can be equal to a virtual linear array, and each virtual antenna element transmits and receives signal by itself. The linear frequency modulated pulse signal transmitted by themth virtual antenna element is given by
(7)
s
(
t
^
)
=
rect
(
t
^
T
r
)
exp
⁡
[
j
2
π
(
f
c
t
^
+
1
2
γ
t
^
2
)
]
,
where t
^ denotes the fast time, T
r denotes the pulse width, f
c is the carrier frequency, γ is the chirp rate, and rect
(
·
) is the unit rectangular function.For an arbitrary point scatterer positioned atP
(
x
p
,
y
p
,
z
p
), the echo signal received by the mth virtual antenna element can be expressed as
(8)
s
1
(
t
^
,
τ
,
y
m
)
=
rect
(
t
^
-
2
R
/
c
T
r
)
rect
(
τ
-
τ
p
T
a
)
×
exp
⁡
[
j
π
γ
(
t
^
-
2
R
c
)
2
-
j
4
π
λ
R
]
,
where τ denotes the azimuth time, τ
p
=
x
p
/
v is the azimuth time of the point scatterer P, c is the light velocity, T
a is the observing duration, λ is the wavelength, and R is the instantaneous distance between the mth virtual antenna element and the point scatterer
(9)
R
=
(
v
τ
-
x
p
)
2
+
(
y
m
-
y
p
)
2
+
(
H
-
z
p
)
2
≈
R
p
+
(
v
τ
-
x
p
)
2
2
R
p
+
y
m
2
-
2
y
m
y
p
2
R
p
,
where R
p
=
(
H
-
z
p
)
2
+
y
p
2 and y
m is the cross-track position of the mth virtual antenna element.


### 3.1. Range Compression

Transform the signal expressed in (8) into the range frequency domain
(10)
S
1
(
f
r
,
τ
,
y
m
)
=
rect
(
f
r
T
r
γ
)
rect
(
τ
-
τ
p
T
a
)
×
exp
⁡
(
-
j
π
f
r
2
γ
)
exp
⁡
(
-
j
4
π
f
r
c
R
)
exp
⁡
(
-
j
4
π
λ
R
)
,
where f
r is range frequency.The matched filter function for the range compression is given by(11)
H
1
(
f
r
)
=
rect
(
f
r
T
r
γ
)
exp
⁡
(
j
π
f
r
2
γ
)
.After the range compression, the received signal can be written as(12)
s
2
(
t
^
,
τ
,
y
m
)
=
rect
(
τ
-
τ
p
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
)
]
×
exp
⁡
(
-
j
4
π
R
λ
)
.


### 3.2. Azimuth Compression

From (9) and (12) we get that the distance between the antenna element and the scatterer varies with the azimuth position of the radar platform, which leads to the coupling of the envelope in the range-azimuth plane. Therefore, the range migration in the azimuth direction should be removed first before the azimuth compression. And the correction function for the range migration can be derived from (9)
(13)
H
2
(
f
r
,
τ
)
=
exp
⁡
(
j
2
π
(
v
τ
-
x
p
)
2
c
R
p
f
r
)
.Transform the signal expressed in (12) to the range frequency domain and correct the range migration using H
2
(
f
r
,
τ
). Then, after the range inverse Fourier transform, the signal becomes
(14)
s
3
(
t
^
,
τ
,
y
m
)
=
rect
(
τ
-
τ
p
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
×
exp
⁡
(
-
j
4
π
R
p
λ
)
·
exp
⁡
[
-
j
2
π
v
2
(
τ
-
τ
p
)
2
λ
R
p
]
×
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
,
where R
c
=
R
p
+
(
y
m
2
-
2
y
m
y
p
)
/
2
R
p.An azimuth Fourier transform is then performed on each range gate to transform the data into the range time-azimuth frequency domain, and the signal becomes(15)
S
3
(
t
^
,
f
a
,
y
m
)
=
rect
(
f
a
γ
a
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
×
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
·
exp
⁡
(
-
j
4
π
R
p
λ
)
×
exp
⁡
[
j
π
f
a
2
γ
a
]
exp
⁡
(
-
j
2
π
τ
p
f
a
)
,
where γ
a
=
2
v
2
/
λ
R
p. Then a multiplication of the signal with the azimuth matched filter function H
3
(
f
a
) is performed, where
(16)
H
3
(
f
a
)
=
exp
⁡
(
-
j
π
f
a
2
γ
a
)
.
Then, we get the azimuth compressed signal by performing inverse Fourier transform in the azimuth direction
(17)
s
4
(
t
^
,
τ
,
y
m
)
=
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
sinc
[
T
a
γ
a
(
τ
-
τ
p
)
]
·
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
exp
⁡
(
-
j
4
π
R
p
λ
)
.


### 3.3. Cross-Track Compression

From (17) it can be seen that the distance between the antenna element and the scatterer also varies with the cross-track position of the antenna element. Therefore, the range migration in the cross-track direction should be removed before the cross-track compression. And the amount of the range migration to be corrected can be given by (9)
(18)
Δ
R
c
=
y
m
2
-
2
y
m
y
p
2
R
p
.Transform the signal expressed in (17) into the range frequency domain, and correct the range migration in the range frequency domain. Then the signal becomes
(19)
s
4
(
t
^
,
τ
,
y
m
)
=
sinc
[
T
r
γ
(
t
^
-
2
R
p
c
)
]
sinc
[
T
a
γ
a
(
τ
-
τ
p
)
]
·
exp
⁡
(
-
j
2
π
y
m
2
λ
R
p
)
exp
⁡
(
j
4
π
y
m
y
p
λ
R
p
)
exp
⁡
(
-
j
4
π
R
p
λ
)
.
The first phase term in (19) represents a quadratic distortion, which can be compensated by
(20)
H
4
(
y
m
)
=
exp
⁡
(
j
2
π
y
m
2
λ
R
p
)
.In high frequency radar application, the interest scene can be modeled by a limited number of strong scattering centers reflecting impinging electromagnetic waves isotropically to all receivers, representing strong spatial sparsity [14]. Therefore, the cross-track imaging process of down-looking MIMO array SAR can be transformed into the problem of sparse signal reconstruction from noisy measurements. After the 2D imaging process in the range and azimuth directions, the signal in the range-azimuth cell corresponding to (
t
^
=
2
R
p
/
c
,
τ
=
τ
p
) by neglecting the constant phase term of (19) can be written as
(21)
s
5
(
t
^
,
τ
,
y
m
)
=
∑
p
=
1
P
σ
p
(
t
^
,
τ
)
exp
⁡
(
j
2
π
f
p
y
m
)
,
where σ
p
(
t
^
,
τ
) and f
p
=
2
y
p
/
λ
R
p denote the backward scattering coefficient and the frequency of the pth point scatterer, respectively.For numerical analysis, (21) can be described by discrete system model
(22)
s
=
Φ
σ
,
where s
=
[
s
5
(
t
^
,
τ
,
y
1
)
,
s
5
(
t
^
,
τ
,
y
2
)
,
…
,
s
5
(
t
^
,
τ
,
y
M
N
)
]
T is the signal vector corresponding to the M
N virtual antenna elements, σ
=
[
σ
1
(
t
^
,
τ
)
,
σ
2
(
t
^
,
τ
)
,
…
,
σ
P
(
t
^
,
τ
)
]
Tis the complex-valued scatter coefficient vector in the cross-track direction, and the matrix Φ can be constructed as
(23)
Φ
=
[
ϕ
1
,
ϕ
2
,
…
,
ϕ
p
,
…
,
ϕ
P
]
,
where
(24)
ϕ
p
=
[
exp
⁡
(
j
2
π
f
p
y
1
)
,
m
m
exp
⁡
(
j
2
π
f
p
y
2
)
,
…
,
exp
⁡
(
j
2
π
f
p
y
M
N
)
]
T
,
where (
·
)
T represents the transpose operation.In the more realistic case some noise is added on the measurements(25)
s
=
Φ
σ
+
n
with n a complex Gaussian vector with zero mean and power σ
n
2.Then, the Bayesian compressive sensing method is employed to estimate theσ. From the aspect of denoising, Laplace distribution is often used as the sparseness prior [15]. Hence, the probability distribution function of σ can be given by
(26)
p
(
σ
∣
σ
δ
)
=
∏
i
=
1
P
p
(
σ
i
∣
σ
δ
)
=
∏
i
=
1
P
1
2
σ
δ
exp
⁡
(
-
2
σ
i
σ
δ
)
,
where σ
δ is the scale parameter of the Laplace distribution.And the probability density function of Gaussian noisen is given by
(27)
p
(
n
∣
σ
n
2
)
=
(
1
2
π
σ
n
2
)
M
N
exp
⁡
(
-
∥
n
∥
2
2
2
σ
n
2
)
,
where ∥
·
∥
2
2 denotes the l
2 norm.Therefore, the likelihood function of the received datas is
(28)
p
(
s
∣
σ
,
σ
n
2
)
=
(
1
2
π
σ
n
2
)
M
N
exp
⁡
{
-
∥
s
-
Φ
σ
∥
2
2
2
σ
n
2
}
.Based on Bayesian theory, the maximum a posteriori (MAP) estimator is used to estimateσ as
(29)
σ
^
=
argmax
[
p
(
σ
∣
s
)
]
=
argmax
[
p
(
s
∣
σ
,
σ
n
2
)
·
p
(
σ
∣
σ
δ
)
]
.Apparently, the MAP solution ofσ can be estimated by maximizing the log posterior of σ as
(30)
σ
^
=
argmax
J
(
σ
)
=
argmax
[
ln
⁡
p
(
s
∣
σ
,
σ
n
2
)
+
ln
⁡
p
(
σ
∣
σ
δ
)
]
=
argmax
[
-
1
2
σ
n
2
∥
s
-
Φ
σ
∥
2
2
-
2
σ
δ
∑
i
=
1
P
σ
i
]
=
argmin
[
∥
s
-
Φ
σ
∥
2
2
+
μ
∥
σ
∥
1
]
,
where μ
=
2
2
σ
n
2
/
σ
δ.Then, the quasi-Newton iterative method with Hessian update scheme is used to solve the optimization problem [16]. The gradient of the objective function J
(
σ
)with respect to σ is given by
(31)
∇
J
(
σ
)
=
H
~
(
σ
)
σ
-
2
Φ
H
s
,
where
(32)
H
~
(
σ
)
=
2
Φ
H
Φ
+
μ
Λ
(
σ
)
,
Λ
(
σ
)
=
diag
⁡
{
(
|
(
σ
)
i
|
2
+
ε
)
-
1
/
2
}
.HereH
~
(
σ
) is used as an approximation to the Hessian, and σ can be obtained from the following quasi-Newton iteration:
(33)
σ
^
(
n
+
1
)
=
σ
^
(
n
)
-
γ
[
H
~
(
σ
^
(
n
)
)
]
-
1
∇
J
(
σ
^
(
n
)
)
,
where γ is the step size. After substituting (31) into (33), we obtain the following iterative algorithm:
(34)
H
~
(
σ
^
(
n
)
)
σ
^
(
n
+
1
)
=
(
1
-
γ
)
H
~
(
σ
^
(
n
)
)
σ
^
(
n
)
+
2
γ
Φ
H
s
.
The iteration is stopped when ∥
σ
^
(
n
+
1
)
-
σ
^
(
n
)
∥
2
2
<
ζ, where ζ is a small positive constant.Then, the 3D image of down-looking MIMO array SAR can be obtained until all the range-azimuth cells have been processed using the same procedure.In conclusion, the 3D imaging processing flow of the proposed method for down-looking MIMO array SAR can be shown in Figure3.Figure 3
The flowchart of proposed method for down-looking MIMO array SAR imaging.


### 3.4. Estimation of the Parametersσ
n
2 and σ
δ

It is clear that the imaging performance of the proposed method is related to the correct selection of the noise levelσ
n
2. Therefore, the noise level must be estimated accurately. Because the noise always distributes evenly and there are many range cells containing noise only in down-looking SAR imaging, the noise level estimation is available by setting an energy-based threshold to select the noise cells [17]. The threshold is given by
(35)
thres
=
E
m
+
[
∑
i
=
1
P
Q
(
E
i
-
E
m
)
2
P
Q
]
1
/
2
,
where E
i is the energy of the ith range cell and E
m denotes the mean energy of all range cells. If the energy of a range cell is below the threshold, this range cell is selected as noise cell. Then, the energy of all the selected cells can be used as the noise level σ
n
2.In addition, the maximum likelihood estimate method is used to choose the objective statistical parameterσ
δ. According to (26), the log-likelihood function for the Laplace distribution can be written as
(36)
ℓ
(
σ
δ
∣
σ
1
,
σ
2
,
…
,
σ
P
)
=
-
P
ln
⁡
2
σ
δ
-
∑
i
=
1
P
2
σ
i
σ
δ
.
Maximizing the log likelihood function (36) with respect to σ
δ produces the following equation:
(37)
σ
^
δ
=
1
P
∑
i
=
1
P
2
σ
i
.
From (37) it is clear that we can average the estimates of all pixel values to obtain the estimation of the statistical parameter σ
δ.


## 3.1. Range Compression

Transform the signal expressed in (8) into the range frequency domain
(10)
S
1
(
f
r
,
τ
,
y
m
)
=
rect
(
f
r
T
r
γ
)
rect
(
τ
-
τ
p
T
a
)
×
exp
⁡
(
-
j
π
f
r
2
γ
)
exp
⁡
(
-
j
4
π
f
r
c
R
)
exp
⁡
(
-
j
4
π
λ
R
)
,
where f
r is range frequency.The matched filter function for the range compression is given by(11)
H
1
(
f
r
)
=
rect
(
f
r
T
r
γ
)
exp
⁡
(
j
π
f
r
2
γ
)
.After the range compression, the received signal can be written as(12)
s
2
(
t
^
,
τ
,
y
m
)
=
rect
(
τ
-
τ
p
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
)
]
×
exp
⁡
(
-
j
4
π
R
λ
)
.


## 3.2. Azimuth Compression

From (9) and (12) we get that the distance between the antenna element and the scatterer varies with the azimuth position of the radar platform, which leads to the coupling of the envelope in the range-azimuth plane. Therefore, the range migration in the azimuth direction should be removed first before the azimuth compression. And the correction function for the range migration can be derived from (9)
(13)
H
2
(
f
r
,
τ
)
=
exp
⁡
(
j
2
π
(
v
τ
-
x
p
)
2
c
R
p
f
r
)
.Transform the signal expressed in (12) to the range frequency domain and correct the range migration using H
2
(
f
r
,
τ
). Then, after the range inverse Fourier transform, the signal becomes
(14)
s
3
(
t
^
,
τ
,
y
m
)
=
rect
(
τ
-
τ
p
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
×
exp
⁡
(
-
j
4
π
R
p
λ
)
·
exp
⁡
[
-
j
2
π
v
2
(
τ
-
τ
p
)
2
λ
R
p
]
×
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
,
where R
c
=
R
p
+
(
y
m
2
-
2
y
m
y
p
)
/
2
R
p.An azimuth Fourier transform is then performed on each range gate to transform the data into the range time-azimuth frequency domain, and the signal becomes(15)
S
3
(
t
^
,
f
a
,
y
m
)
=
rect
(
f
a
γ
a
T
a
)
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
×
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
·
exp
⁡
(
-
j
4
π
R
p
λ
)
×
exp
⁡
[
j
π
f
a
2
γ
a
]
exp
⁡
(
-
j
2
π
τ
p
f
a
)
,
where γ
a
=
2
v
2
/
λ
R
p. Then a multiplication of the signal with the azimuth matched filter function H
3
(
f
a
) is performed, where
(16)
H
3
(
f
a
)
=
exp
⁡
(
-
j
π
f
a
2
γ
a
)
.
Then, we get the azimuth compressed signal by performing inverse Fourier transform in the azimuth direction
(17)
s
4
(
t
^
,
τ
,
y
m
)
=
sinc
[
T
r
γ
(
t
^
-
2
R
c
c
)
]
sinc
[
T
a
γ
a
(
τ
-
τ
p
)
]
·
exp
⁡
[
-
j
2
π
(
y
m
2
-
2
y
m
y
p
)
λ
R
p
]
exp
⁡
(
-
j
4
π
R
p
λ
)
.


## 3.3. Cross-Track Compression

From (17) it can be seen that the distance between the antenna element and the scatterer also varies with the cross-track position of the antenna element. Therefore, the range migration in the cross-track direction should be removed before the cross-track compression. And the amount of the range migration to be corrected can be given by (9)
(18)
Δ
R
c
=
y
m
2
-
2
y
m
y
p
2
R
p
.Transform the signal expressed in (17) into the range frequency domain, and correct the range migration in the range frequency domain. Then the signal becomes
(19)
s
4
(
t
^
,
τ
,
y
m
)
=
sinc
[
T
r
γ
(
t
^
-
2
R
p
c
)
]
sinc
[
T
a
γ
a
(
τ
-
τ
p
)
]
·
exp
⁡
(
-
j
2
π
y
m
2
λ
R
p
)
exp
⁡
(
j
4
π
y
m
y
p
λ
R
p
)
exp
⁡
(
-
j
4
π
R
p
λ
)
.
The first phase term in (19) represents a quadratic distortion, which can be compensated by
(20)
H
4
(
y
m
)
=
exp
⁡
(
j
2
π
y
m
2
λ
R
p
)
.In high frequency radar application, the interest scene can be modeled by a limited number of strong scattering centers reflecting impinging electromagnetic waves isotropically to all receivers, representing strong spatial sparsity [14]. Therefore, the cross-track imaging process of down-looking MIMO array SAR can be transformed into the problem of sparse signal reconstruction from noisy measurements. After the 2D imaging process in the range and azimuth directions, the signal in the range-azimuth cell corresponding to (
t
^
=
2
R
p
/
c
,
τ
=
τ
p
) by neglecting the constant phase term of (19) can be written as
(21)
s
5
(
t
^
,
τ
,
y
m
)
=
∑
p
=
1
P
σ
p
(
t
^
,
τ
)
exp
⁡
(
j
2
π
f
p
y
m
)
,
where σ
p
(
t
^
,
τ
) and f
p
=
2
y
p
/
λ
R
p denote the backward scattering coefficient and the frequency of the pth point scatterer, respectively.For numerical analysis, (21) can be described by discrete system model
(22)
s
=
Φ
σ
,
where s
=
[
s
5
(
t
^
,
τ
,
y
1
)
,
s
5
(
t
^
,
τ
,
y
2
)
,
…
,
s
5
(
t
^
,
τ
,
y
M
N
)
]
T is the signal vector corresponding to the M
N virtual antenna elements, σ
=
[
σ
1
(
t
^
,
τ
)
,
σ
2
(
t
^
,
τ
)
,
…
,
σ
P
(
t
^
,
τ
)
]
Tis the complex-valued scatter coefficient vector in the cross-track direction, and the matrix Φ can be constructed as
(23)
Φ
=
[
ϕ
1
,
ϕ
2
,
…
,
ϕ
p
,
…
,
ϕ
P
]
,
where
(24)
ϕ
p
=
[
exp
⁡
(
j
2
π
f
p
y
1
)
,
m
m
exp
⁡
(
j
2
π
f
p
y
2
)
,
…
,
exp
⁡
(
j
2
π
f
p
y
M
N
)
]
T
,
where (
·
)
T represents the transpose operation.In the more realistic case some noise is added on the measurements(25)
s
=
Φ
σ
+
n
with n a complex Gaussian vector with zero mean and power σ
n
2.Then, the Bayesian compressive sensing method is employed to estimate theσ. From the aspect of denoising, Laplace distribution is often used as the sparseness prior [15]. Hence, the probability distribution function of σ can be given by
(26)
p
(
σ
∣
σ
δ
)
=
∏
i
=
1
P
p
(
σ
i
∣
σ
δ
)
=
∏
i
=
1
P
1
2
σ
δ
exp
⁡
(
-
2
σ
i
σ
δ
)
,
where σ
δ is the scale parameter of the Laplace distribution.And the probability density function of Gaussian noisen is given by
(27)
p
(
n
∣
σ
n
2
)
=
(
1
2
π
σ
n
2
)
M
N
exp
⁡
(
-
∥
n
∥
2
2
2
σ
n
2
)
,
where ∥
·
∥
2
2 denotes the l
2 norm.Therefore, the likelihood function of the received datas is
(28)
p
(
s
∣
σ
,
σ
n
2
)
=
(
1
2
π
σ
n
2
)
M
N
exp
⁡
{
-
∥
s
-
Φ
σ
∥
2
2
2
σ
n
2
}
.Based on Bayesian theory, the maximum a posteriori (MAP) estimator is used to estimateσ as
(29)
σ
^
=
argmax
[
p
(
σ
∣
s
)
]
=
argmax
[
p
(
s
∣
σ
,
σ
n
2
)
·
p
(
σ
∣
σ
δ
)
]
.Apparently, the MAP solution ofσ can be estimated by maximizing the log posterior of σ as
(30)
σ
^
=
argmax
J
(
σ
)
=
argmax
[
ln
⁡
p
(
s
∣
σ
,
σ
n
2
)
+
ln
⁡
p
(
σ
∣
σ
δ
)
]
=
argmax
[
-
1
2
σ
n
2
∥
s
-
Φ
σ
∥
2
2
-
2
σ
δ
∑
i
=
1
P
σ
i
]
=
argmin
[
∥
s
-
Φ
σ
∥
2
2
+
μ
∥
σ
∥
1
]
,
where μ
=
2
2
σ
n
2
/
σ
δ.Then, the quasi-Newton iterative method with Hessian update scheme is used to solve the optimization problem [16]. The gradient of the objective function J
(
σ
)with respect to σ is given by
(31)
∇
J
(
σ
)
=
H
~
(
σ
)
σ
-
2
Φ
H
s
,
where
(32)
H
~
(
σ
)
=
2
Φ
H
Φ
+
μ
Λ
(
σ
)
,
Λ
(
σ
)
=
diag
⁡
{
(
|
(
σ
)
i
|
2
+
ε
)
-
1
/
2
}
.HereH
~
(
σ
) is used as an approximation to the Hessian, and σ can be obtained from the following quasi-Newton iteration:
(33)
σ
^
(
n
+
1
)
=
σ
^
(
n
)
-
γ
[
H
~
(
σ
^
(
n
)
)
]
-
1
∇
J
(
σ
^
(
n
)
)
,
where γ is the step size. After substituting (31) into (33), we obtain the following iterative algorithm:
(34)
H
~
(
σ
^
(
n
)
)
σ
^
(
n
+
1
)
=
(
1
-
γ
)
H
~
(
σ
^
(
n
)
)
σ
^
(
n
)
+
2
γ
Φ
H
s
.
The iteration is stopped when ∥
σ
^
(
n
+
1
)
-
σ
^
(
n
)
∥
2
2
<
ζ, where ζ is a small positive constant.Then, the 3D image of down-looking MIMO array SAR can be obtained until all the range-azimuth cells have been processed using the same procedure.In conclusion, the 3D imaging processing flow of the proposed method for down-looking MIMO array SAR can be shown in Figure3.Figure 3
The flowchart of proposed method for down-looking MIMO array SAR imaging.


## 3.4. Estimation of the Parametersσ
n
2 and σ
δ

It is clear that the imaging performance of the proposed method is related to the correct selection of the noise levelσ
n
2. Therefore, the noise level must be estimated accurately. Because the noise always distributes evenly and there are many range cells containing noise only in down-looking SAR imaging, the noise level estimation is available by setting an energy-based threshold to select the noise cells [17]. The threshold is given by
(35)
thres
=
E
m
+
[
∑
i
=
1
P
Q
(
E
i
-
E
m
)
2
P
Q
]
1
/
2
,
where E
i is the energy of the ith range cell and E
m denotes the mean energy of all range cells. If the energy of a range cell is below the threshold, this range cell is selected as noise cell. Then, the energy of all the selected cells can be used as the noise level σ
n
2.In addition, the maximum likelihood estimate method is used to choose the objective statistical parameterσ
δ. According to (26), the log-likelihood function for the Laplace distribution can be written as
(36)
ℓ
(
σ
δ
∣
σ
1
,
σ
2
,
…
,
σ
P
)
=
-
P
ln
⁡
2
σ
δ
-
∑
i
=
1
P
2
σ
i
σ
δ
.
Maximizing the log likelihood function (36) with respect to σ
δ produces the following equation:
(37)
σ
^
δ
=
1
P
∑
i
=
1
P
2
σ
i
.
From (37) it is clear that we can average the estimates of all pixel values to obtain the estimation of the statistical parameter σ
δ.


## 4. Simulation Results

In this section, point target simulation is carried out to verify the validity of the proposed imaging algorithm. The main parameters used for simulation are listed in Table1.Table 1
Parameters used for simulation.
Parameter
Value
Carrier frequency
37.5 GHz
Pulse bandwidth
300 MHz
Pulse repetition frequency
1024 Hz
Chirp duration
1.0μs
Radar height
500 m
Radar velocity
50 m/s
Number of transmitting antenna elements
4
Number of receiving antenna elements
32
Azimuth resolution
0.4 m
Range resolution
0.5 m
Cross-track resolution
0.4 mSuppose that there are five point targets located at the scene with the azimuth-range-cross track values equal to (0, 485, 0), (8, 495, 20), (8, 495, −20), (−8, 495, 20), and (−8, 495, −20), respectively. The distributions of the five point targets are shown in Figure4(a). After raw data generation and 3D imaging processing by using the proposed algorithm, the surfaces of the final 3D image are plotted at −20 dB in Figure 4(b). As expected, the image is reconstructed in 3D space, and the whole space structure is very consistent with the real situation in Figure 4(a). Figure 5 shows three selected sections of the final 3D image of down-looking MIMO array SAR. Figure 5(a) shows the 2D image of the selected range-cross track section corresponding to azimuth position 8 m. Figure 5(b) shows the 2D image of the selected azimuth-cross track section corresponding to range position 495 m. Figure 5(c) shows the 2D image of the selected azimuth-range section corresponding to cross-track position −20 m. The below imaging results show that the point scatterers are well focused in three directions, confirming the validity of the proposed algorithm.Real spatial position and final 3D image.
(a)
Real spatial position
(b)
Final 3D image2D image of selected sections.
(a)
Selected range-cross track section
(b)
Selected azimuth-cross track section
(c)
Selected azimuth-range sectionIn order to analyze the performance of the proposed method, the imaging result obtained by Fourier transform is given for comparison. Supposed that there are two targets located at the azimuth-range-cross track unit of (8, 0.18, 495) and (8, −0.18, 495). Figure6(a) shows the range-cross track reconstruction result obtained by Fourier transform, and Figure 6(b) shows the range-cross track reconstruction result obtained by the proposed method. As the cross-track distance of the two targets is 0.36 m, which is less than the cross-track resolution of 0.4 m, the two targets can not be distinguished in the image obtained by Fourier based method. However, the proposed method can improve the spatial resolution and distinguish the two targets clearly.Comparison of the range-cross track imaging results of down-looking MIMO array SAR.
(a)
Imaging result obtained by Fourier transform
(b)
Imaging result obtained by proposed methodIn the following experiment, we take into account two targets located at the same azimuth and range positions and with the cross-track values equal to −2 m and 2 m, respectively. Moreover, additional Gaussian distributed complex noise is added to generate measurements with SNR of 5 dB. Figures7(a) and 7(b) show the cross-track distribution of the two targets obtained by Fourier transform and the proposed method, respectively. By comparing the imaging results, it can be seen that the proposed method is more robust to noise and it recovers target image together with suppressing the noise components.Comparison of the cross-track distribution obtained by Fourier transform and the proposed method.
(a)
Cross-track distribution obtained by Fourier transform
(b)
Cross-track distribution obtained by proposed method


## 5. Conclusions

Down-looking MIMO array SAR can reconstruct 3D images of the observed area and overcome restrictions of shading and lay over effects in side-looking SAR. Therefore, down-looking MIMO array SAR has challenging potential for 3D digital maps, complex terrain mapping, and so on. However, the cross-track resolution of down-looking MIMO array SAR is limited by the length of linear array. In this paper, a novel 3D imaging strategy is proposed for down-looking MIMO array SAR. Exploiting the spatial sparsity of the interest scene, we transform the cross-track imaging process of down-looking MIMO array SAR into the problem of sparse signal reconstruction from noisy measurements. Raw data of down-looking MIMO array SAR in Ka-band is simulated and the 3D image is achieved. The results of the simulated data confirm the effectiveness of the proposed method.



---

*Source: 612326-2014-05-11.xml*