Three-Dimensional Microwave Imaging Using Synthetic Aperture Technique

1. Introduction

With the ability of two-dimensional (2-D) microwave imaging, synthetic aperture radar (SAR) has been an important imaging tool for civilian and military applications. The basic idea of 2-D SAR is to synthesize a linear array by moving a high-range-resolution (HRR) radar long a straight path, and obtain the additional azimuthal resolution. To exact the height information from the 2-D SAR, interferometric SAR (InSAR) technique, which requires multiple antennas or repeated flight paths, has been developed and is widely used for remote sensing applications.

However, since the interferometric SAR technique is based on the 2-D SAR images, it will be invalid, when there is more than one scatterer projected in the same pixel of the 2-D SAR image. This disadvantage makes it difficult to be used in high-precision 3-D RCS measurement and topographical survey in urban region.

To improve the ability of microwave remote sensing, some new 3-D SAR systems, such as circle SAR, elevation circular SAR, curve SAR, and linear array SAR have been developed based on the synthetic aperture technique. The basic idea of them is to produce 2-D resolution by moving the HRR radar in 2-D / 3-D space and obtain the third dimensional resolution using pulse compression technique. This chapter will discuss the principle and imaging processing technique of 3-D SAR.

In section 2, an approach to calculate the oscillatory integral has been introduced, which could simplify the analysis of 3-D SAR ambiguity function. In section 3, the ambiguity function and spatial resolution of 3-D SAR are discussed. The backprojection method and experiment data processing are presented in section 4 and 5 respectively. The multiresolution approximation techniques that can reduce the computational cost of 3-D SAR are discussed in section 6.

2. Preliminary

Compared with the traditional SAR, the echo model of 3-D SAR is more complex. To simplify the analysis of 3-D SAR, the calculation of oscillatory integral using density function will be introduced in this section.

2.1. A Simple Example

At first, let’s observe a simple example. Assume that there is a discrete function $$f(n)=(\mathrm{0,1}/\mathrm{4,2}/\mathrm{4,0,1}/\mathrm{4,1}/4)$$ , and we need to calculate the sum of $$\exp (j\cdot 2\pi \cdot f(n))$$ . Obviously, we have:

$${\displaystyle \sum _n\exp (j\cdot 2\pi \cdot f(n))}=e^{j\cdot 2\pi \cdot 0}+e^{j\cdot 2\pi \cdot 1/4}+e^{j\cdot 2\pi \cdot 2/4}+e^{j\cdot 2\pi \cdot 0}+e^{j\cdot 2\pi \cdot 1/4}+e^{j\cdot 2\pi \cdot 1/4}$$ E1

According to the commutative law of addition, eq. (1) can be rewritten as:

$${\displaystyle \sum _n\exp (j\cdot 2\pi \cdot f(n))}=2\cdot e^{j\cdot 2\pi \cdot 0}+3\cdot e^{j\cdot 2\pi \cdot 1/4}+1\cdot e^{j\cdot 2\pi \cdot 2/4}+0\cdot e^{j\cdot 2\pi \cdot 3/4}$$ E2

where, coefficients 2, 3, 1 and 0 represent the frequency (in the sense of probability, which is defined as the number of times that value occurs in the data set) of every exponential term. Thus, by introducing the concept of density function $$D(i)=\left\{\mathrm{2,3,1,0}\right\}$$ , we have:

$${\displaystyle \sum _n\exp (j\cdot 2\pi \cdot f(n))}={\displaystyle \sum _{i=0}^{3}D(i)e^{j\cdot 2\pi \cdot i/4}}$$ E3

Obviously, eq. (3) matches the definition of Discrete Fourier Transform (DFT). Denoting the DFT of $$D(i)$$ as $$\widehat{D}(k)$$ , we have:

$${\displaystyle \sum _n\exp (j\cdot 2\pi \cdot f(n))}=\widehat{D}(1)$$ E4

Similarly, we have:

$${\displaystyle \sum _n\exp (j\cdot 2\pi \cdot k\cdot f(n))}={\displaystyle \sum _{i=0}^{3}D(i)e^{j\cdot 2\pi \cdot k\cdot i/4}}=\widehat{D}(k)$$ E5

This example indicates that the exponential sum of a finite discrete-time function can be calculated using its density function. However, since the commutative law of addition holds only when $$f(n)$$ is finite, and the concept of frequency is meaningful for finite set, a more precise definition of density function using Lebesgue measure and a theorem that extends eq. (4) to the continuous-time function will be presented in next subsection.

2.2. Calculation of the oscillatory integral using density function

According to measure theory, functions are divided into four classes, simple function, Bounded Function Supported on a set of Finite Measure (BFFM), non-negative function and integrable function (the general case). For the analysis of array whose size is finite, the BFFM assumption is sufficient.

Given a BFFM $$f(t)$$ with support $$F$$ , analogous to the definition of Cumulative Density Function (CDF) in probability theory, define the Cumulative Density Function of $$f(t)$$ as:

$$C_f(y)\triangleq m(F_y)$$ E6

$$F_y\triangleq \{t:f(t)\le y;y\in ℝ\}$$ E7

where, $$F_y$$ is the subset of $$F$$ . $$m(F_y)$$ denotes the Lebesgue measure of $$F_y$$ , which describe the volume (area) of $$F_y$$ . Especially, when $$F$$ is finite set, $$m(F_y)$$ is the cardinality of subset $$F_y$$ .

Then we define the Density Function (DF) of $$f(t)$$ as the derivative of $$C_f(y)$$ :

$$D_f(y)\triangleq d{C}_f(y)/dy$$ E8

Obviously, $$D_f(y)$$ satisfies:

$$D_f(y)\ge 0$$ E9

$${\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{D}_f(y)dy}=m(F)+\infty$$ E10

Properties 1 and 2 of $$D_f(y)$$ indicate that $$D_f(y)$$ is absolutely integrable, and its Fourier transform exists.

Using the concept of density function, we can calculate an oscillatory integral using the following theorem.

Theorem 1: Given a BFFM phase function

$$f(t)$$ supported on a set $$E$$ , we have:

$${\displaystyle {\int }_E{e}^{jf(t)}dt}={\widehat{D}}_f(1)$$ E11

where, $${\widehat{D}}_f(·)$$ denotes the Fourier transform of $$D_f(y)$$

Proof:

The proof of theorem 1 includes two steps: firstly, we consider $$f(t)$$ as a simple function, and then extend the conclusion to BFFM function.

According to the definition in measure theory, a simple function is a finite sum of a group of characteristic functions:

$$\begin{array}{l}f(t)={\displaystyle \sum _{k=1}^K{a}_k{\chi }_{E_k}(t)}\\ {\chi }_{E_k}(t)=\{\begin{array}{c}\begin{array}{cc}1& t\in E_k\end{array}\\ \begin{array}{cc}0& t\notin E_k\end{array}\end{array}\end{array}$$ E12

where, $$a_k$$ is constant, $$E_k$$ denotes a measurable subset of set $$E$$ $${\chi }_{E_k}(t)$$ denotes the characteristic function of $$E_k$$ .

Case 1: rational number

Assume that $$a_k$$ are all rational number, there exists an equal-interval infinite rational number set:

$$B=\{-\infty \mathrm{,...,}\frac{-1}{Q},0,\frac{1}{Q}\mathrm{,...,}+\infty \},Q\in \mathbb{N}$$ E13

satisfying $$a_k\in B$$ in the example, $$\left\{a_k\right\}=\{\mathrm{0,1}/\mathrm{4,2}/4\}$$

$$B=\{\mathrm{0,1}/\mathrm{4,2}/\mathrm{4,3}/4\}$$ E14

Using set $$B$$ , we can construct a group of characteristic function $${\chi }_{F_i}(t)$$ , and $$f(t)$$ could be written as:

$$f(t)={\displaystyle \sum _{i=-\infty }^{+\infty }{b}_i{\chi }_{F_i}(t)}$$ E15

where, $$b_i\in B$$ $$F_i=E_k$$ , when $$b_i=a_k$$ ; otherwise, $$F_k=\Phi$$ .

According to the definition of Lebesgue integral, we have:

$${\displaystyle {\int }_E{e}^{jf(t)}dt}={\displaystyle \sum _{i=-\infty }^{+\infty }{e}^{j{b}_i}m(F_i)={\displaystyle \sum _{i=-\infty }^{+\infty }{e}^{j{b}_i}{D}_f(i)}}$$ E16

Case 2: real number

Assume that $$a_k$$ are all real number, according to the real analysis, for every real number $$a$$ , there exists a sequence $$\text{{}{a}_n\}$$ of rational numbers can approximate to it. Thus, there exists a sequence $$\text{{}{B}_Q\}$$ of rational number sets can approximate to all of the $$a_k$$ , i.e.:

$${\displaystyle {\int }_E{e}^{jf(t)}dt}=\underset{Q\to +\infty }{\lim }{\displaystyle \sum _{i=-\infty }^{+\infty }{e}^{j{b}_i^Q}{D}_f^Q(i)}$$ E17

with the increase of $$Q$$ , the interval $$1/Q$$ of $$B$$ trends toward zero, and eq. (13) can be written in integral form as:

$${\displaystyle {\int }_E{e}^{jf(t)}dt}={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{e}^{jy}{D}_f(y)dy}$$ E18

Since the Fourier transform of $$D_f(y)$$ exists, we have:

$${\displaystyle {\int }_E{e}^{jf(t)}dt}={\widehat{D}}_f(1)$$ E19

For a BFFM function $$f(t)$$ bounded by $$M$$ and supported on a set $$E$$ , there exists a sequence $$\text{{}{f}_n\}$$ of simples functions, with each $$f_n$$ bounded by $$M$$ and supported on a set $$E$$ , and such that:

$$f_n(t)\to f(t)forallt$$ E20

Thus, eq. (8) holds for all BFFM functions.

Theorem 1 provides a method to calculate the oscillatory integral without any approximation. Compared with the principle of stationary phase (PSP), this method does not need $$f(t)$$ be derivable, and holds for all BFFM function.

Using theorem 1, we can obtain the following corollary directly by rewritten eq. (12) as:

$${\displaystyle {\int }_E{e}^{j\cdot u\cdot f(t)}dt}={\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{e}^{j\cdot u\cdot y}{D}_f(y)dy}={\widehat{D}}_f(u)$$ E21

Corollary 1:

$${\displaystyle {\int }_E{e}^{j\cdot u\cdot f(t)}dt}={\widehat{D}}_f(u)$$ E22

As it will be seen in the next section, this corollary is crucial for the analysis on the ambiguity function of 3-D SAR.

3. Principle of 3-D SAR

3.1. Introduction on the typical 3-D SARs

In this subsection, a brief discussion on the typical 3-D SAR systems including circle SAR (CSAR), elevation circular SAR (E-CSAR), curve SAR, and linear array SAR (LASAR) will be proposed.

In 1999, Tsz-King Chan, Yasuo Kuga, and Akira Ishimaru proposed a novel method for radar topographical imaging which required the SAR platform to move in a circular orbit, and named it as circular SAR. In their experiment, the transmitting and receiving antennas were mounted on two separate wooden rings that were individually driven by stepping motors with an angular precision of approximately 0.02. Imaging result of a model helicopter of length 30 cm has been obtained and published (Tsz-King Chan; Kuga, Y.; Ishimaru, A., 1999). In 2001, at the Radar Division of Georgia Tech Research Institute (GTRI), a 3-D inverse synthetic aperture radar (SAR) system has been developed that performs synthetic aperture measurement via a linear motion of the radar in the elevation domain, and a circular (turntable) motion of the target in the range and cross-range domains, which was named as elevation circular SAR (E-CSAR) system. Its geometry is shown in Figure 1.

Figure 1.

Geometry of E-CSAR.

The radar motion in elevation provides target coherent radar cross section (RCS) as a function of the elevation (or depression) angle. The target’s circular motion yields the azimuthal look angle information. The imaging results of T-72 tank have been obtained and published (Bryant, M.L., Gostin, L.L., Soumekh, M. 2003). In fact, the concept “elevation circular SAR” could be extended by controlling the radar moving around the target in a helix trajectory which is shown in Figure 2. The cylindrical surface produces the 2-D resolution vertical (approximately) to the range resolution. Furthermore, to simplify the motion control, the helix trajectory could be composed by the circle motion of the radar and the rectilinear motion of the target, which is shown in Figure 3.

Figure 2.

Geometry of E-CSAR with helix trajectory

Figure 3.

Geometry of E-CSAR composing helix trajectory by circle motion and rectilinear motion

The advantage of E-CSAR is that we can obtain the RCSes of the target in different elevation angle and azimuthal angle in one observation session; its disadvantage is that the target’s size must be smaller than the diameter of the cylinder, which makes it difficult to be employed for large-size target. In fact, since the RCS of the target varies with the elevation angle and azimuthal angle, the synthetic aperture is a local region of the cylndrical surface, which is illustrated in Figure 4 (left).

Figure 4.

Approximation of a local region of the E-CSAR using a synthetic plane aperture.

Thus, we could approximate the local region as a plane and control the antenna phase centre moving in the plane by using a 2-D motion control platform or a linear array mounted on a 1-D motion control platform, which is shown in Figure 4 (right), and obtain the radar cross section (RCS) in one specific direction in one observation session. To obtain the RCSes in different elevation angle and azimuthal angle, one can just rotate the target or the platform. Compared with the E-CSAR, the size of the synthetic plane aperture is small and could be used in 3-D RCS measurement for large-size target.

Besides the E-CSAR, the curve SAR has also been researched in the radar community. In 1995, Jennifer L.H., Webb and David C. Munson, Jr. considered the problem of spotlight-mode synthetic aperture radar (SAR) imaging for an arbitrary radar path to reconstruct a 2-D image of 3-D surfaces. In 2004, Sune R. J. Axelsson researched the beam characteristics of 3-D SAR in curved or random paths in detail, and concluded that the SAR sidelobe suppression of a single circle path was worse than that of a circular antenna of similar size due to the fact that only a line boundary was used as SAR aperture. The spiral paths and random paths were discussed to improve the beam characteristic of 3-D SAR. The geometry of typical curve SAR is shown in Figure 5. The radar is mounted on a platform with curve path to synthesize a 2-D aperture.

Figure 5.

Geometry of curve SAR.

The advantage of curve SAR is that the size of the synthetic aperture could be far larger than the other 3-D SAR systems, which means high-resolution in the cross-track (x) direction; its disadvantage is that the motion control is too difficult to be implemented for the application of topographical survey in practice.

In 1996, Bassem R. Mahafza and Mitch Sajjadi proposed the concept “linear array SAR”, which mounted a linear array on a platform with rectilinear motion and synthesized a 2-D plane array. Its geometry is shown in Figure 6. In 2004, R. Giret, H. Jeul and, P. Enert conceived a millimeter-wave imaging radar onboard an UAV, and designed a 3-D millimeter-wave imaging plan. In their plan, a linear array was mounted above the ground (perpendicular to the ground plane) and vehicles passed under the system to obtain the 3-D image of the vehicles.

Figure 6.

Geometry of linear array SAR.

Compared with the curve SAR, the motion control of linear array SAR is simpler. While, to achieve high cross-track resolution, the linear array must be rather long and the number of element is large, which is difficult and expensive to be implemented. To reduce the system complexity and cost, M. Weiß and J.H.G. Ender introduced the concept “MIMO radar” into the linear array SAR. With this concept, one can synthesize a sparse linear array SAR with relatively low cost, whose equivalent geometry is shown in Figure 7.

Figure 7.

Geometry of sparse LASAR, only the colored elements active at one pulse repetition period.

The disadvantage of LASAR and sparse LASAR is that the cross-track resolution is determined by the length of the linear array. Since the length of the linear array is limited by the size of the platform, its cross-track resolution will be the bottleneck. Theoretically, the curve SAR could be considered as a kind of sparse LASAR (shown in Figure 7 in the blue dash-line).

In a word, the key problem of 3-D SAR is to vary the position of antenna phase centre (APC) in the 2-D / 3-D space. This work could be implemented mechanically (such as 2-D motion control platform and aircraft), or electrically (such as linear array). By moving HRR radar in 2-D plane using high precision motion control platform, we can build a low-cost 3-D RCS measurement device. The linear array SAR with MIMO technique might be the most feasible 3-D SAR system for the topographical survey application, thought there are still some problems, such as, the balance between the length of linear array and the cross-track resolution and the compensation of motion measurement error.

3.2. General echo model of 3-D SAR

For the traditional SAR, the echo is always considered as a function of fast-time $$\tau$$ and slow-time $$t$$ . While, for 3-D SAR, there might be more than one channel echo received at one pulse repetition period, such as LASAR and sparse LASAR, and it is not convenient to describe the 3-D SAR echo using slow-time $$t$$ . In fact, the synthetic aperture technique produces additional resolution by moving the antenna phase centre in the spatial domain, and we should pay more attention on the change of the antenna phase centre in spatial domain rather than that in the time domain. Thus, a general echo model is built in this subsection, which can describe different 3-D SAR systems.

Given a scatterer with position ${\overline{P}}_{\omega }$ , its slant range to the antenna phase centre with position $${\overline{P}}_{apc}$$ is:

$$R({\overline{P}}_{\omega }{\overline{P}}_{apc})\triangleq {\Vert {\overline{P}}_{apc}-{\overline{P}}_{\omega }\Vert }_2$$ E23

where, $${\Vert ·\Vert }_2$$ denotes the 2-norm of vector.

Given the transmitted baseband signal $$f(t)$$ , ignoring the radiation pattern, the scatterer’s echo $$D(\tau ;{\overline{P}}_{\omega };{\overline{P}}_{apc})$$ can be written as:

$$D(\tau ;{\overline{P}}_{\omega };{\overline{P}}_{apc})=\exp (j\cdot 2\pi \cdot 2R({\overline{P}}_{\omega }{\overline{P}}_{apc})/\lambda )\cdot f(\tau -2R({\overline{P}}_{\omega },{\overline{P}}_{apc})/c)$$ E24

where, $$\tau$$ denotes the fast time domain, $$\lambda$$ denotes the wave length of the carrier. The first term in eq.(21) is the Doppler term arising from the relative position changes of the antenna phase centre with respect to the target. The second term is the fast-time term which causes the range resolution. Note that, in some cases, the transmitter and receiver might be operated independently, and term $$2R({\overline{P}}_{\omega },{\overline{P}}_{apc})$$ in eq. (21) should be rewritten as $$R({\overline{P}}_{\omega },{\overline{P}}_{apc}^T)$$ + $$R({\overline{P}}_{\omega },{\overline{P}}_{apc}^R)$$ , and the analysis should be modified correspondingly.

To describe the relative position changes, we introduce the concept “antenna phase centre set” $$P$$ denoting the collection of the positions of the antenna phase centre (note that the elements in the antenna phase centre set might be repetitive.).

For traditional 2-D SAR, its antenna phase centre set could be expressed as:

$$P=\{x,y,z|x=x_0,y=v\cdot t,z=z_0;t\in T\}$$ E25

where, $$v$$ denotes the speed of the platform, $$T$$ denotes the slow-time domain, $$z_0$$ denotes the height of the platform.

For E-CSAR, we have:

$$P=\{x,y,z|x=\upsilon \cos (t),y=\upsilon \sin (t),z=v_{h}t;t\in T\}$$ E26

where, $$\upsilon$$ denotes the radius of the cylinder, $$v_h$$ denotes the speed in the vertical direction.

For curve SAR, we have:

$$P=\{x,y,z|x=x(t),y=v\cdot t,z=z_0;t\in T\}$$ E27

where, $$x(t)$$ and $$v\cdot t$$ compose the curve trajectory.

For linear array SAR, we have:

$$P=\{x,y,z|x\in X,y=v\cdot t,z=z_0;t\in T\}$$ E28

where, $$X$$ denotes the set of the x positions of the linear array, e.g., $$X=\{x|i\cdot d;i=\mathrm{0,1,...,}N-1\}$$ $$N$$ denotes the element number of the linear array, $$d$$ denotes the element interval.

For sparse LASAR, $$P$$ is a group of random positions, and we just simplify denote it as $$P$$

Using the antenna phase centre set, we can easily express the echo of 3-D SAR as:

$$D(r;{\overline{P}}_{\omega })=\{D|D(r;{\overline{P}}_{\omega };{\overline{P}}_{apc});{\overline{P}}_{apc}\in P\}$$ E29

Note that, given r, $$D$$ is a set defined under the antenna phase centre set $$P$$ rather than a number.

Thought this echo model is more abstract than the classical one, as it will be seen in the next subsection, it will simplify the analysis of 3-D SAR ambiguity function. We could build the direct relationship between the antenna phase centre set $$P$$ (describes the shape of the synthetic aperture) and its ambiguity function, which make it easy for the 3-D SAR analysis and design.

3.3. Ambiguity function

Ambiguity function (AF) is one of the crucial concepts in the radar theory. For the pulse-Doppler (PD) radar, ambiguity function is a 2-D function of time delay and Doppler frequency. For imaging radar, it describes the interaction of different scatterers in the image space, which is also called point spread function. A well-designed imaging radar should have narrow mainlobe, low peak sidelobe ratio (PSLR) and low integrated sidelobe ratio (ISLR). In this section, we will discuss the ambiguity function of 3-D SAR.

Based on the echo model built in last subsection, the ambiguity function $$\chi ({\overline{P}}_{\omega })$$ of 3-D SAR can be defined as:

$$\chi ({\overline{P}}_{\omega })\triangleq \frac{{\displaystyle \sum _p{\displaystyle \int D(\tau ;{\overline{P}}_{\omega })\cdot D_{}^{\text{*}}(\tau ;\overline{0}})d\tau }}{{\displaystyle \sum _p{{\displaystyle \int \left|D[\tau ;\overline{0}]\right|}}^{2}d\tau }}$$ E30

where, superscript $$*$$ denotes complex conjugate, $$\overline{0}$$ denotes the position of the reference point.

Since the integration with respect to the fast time $$\tau$$ in eq. (27) is the range-compression operation, eq. (27) can be rewritten as:

$$\chi ({\overline{P}}_{\omega })=\frac{1}{M}{\displaystyle \sum _p\exp (j\cdot 2\pi \cdot [R({\overline{P}}_{\omega }{\overline{P}}_{apc})-R(\overline{0},{\overline{P}}_{apc})]/\lambda )\cdot {\chi }_{}^R(r-2R(\overline{0},{\overline{P}}_{apc}))}$$ E31

where, $$r$$ denotes the range domain, $${\chi }_{}^R(r)$$ denotes the ambiguity function in the range direction, which is a sinc function.

Approximating $${\chi }_{}^R(r)$$ as the impulse function, the range AF in eq.(28) could be moved out of the summation by range migration adjustment during imaging processing, and eq. (28) could be rewritten as:

$$\chi ({\overline{P}}_{\omega })\approx \frac{1}{M}\left\{{\displaystyle \sum _p\exp (j\cdot 2\pi \cdot [R({\overline{P}}_{\omega },{\overline{P}}_{apc})-R(\overline{0}{\overline{P}}_{apc})]/\lambda )}\right\}\cdot {\chi }_{}^R(r-2R(\overline{0}{\overline{P}}_0))$$ E32

where, $${\overline{P}}_0$$ denotes the centre of the synthetic aperture.

From eq.(29), the ambiguity function of 3-D SAR is the product of two terms: the first term in the curly brace is the Doppler term, which is caused by the synthetic aperture and produces resolution in the aperture direction(s); the second term is the fast-time term which produces range resolution. We define the synthetic aperture ambiguity function as:

$$\begin{array}{l}{\chi }^p({\overline{P}}_{\omega })\triangleq \frac{1}{M}\left\{{\displaystyle \sum _p\exp (j\cdot 2\pi \cdot \Delta R({\overline{P}}_{\omega },{\overline{P}}_{apc})/\lambda )}\right\}\\ \Delta R_{}^{\omega }({\overline{P}}_{apc})\triangleq R({\overline{P}}_{\omega },{\overline{P}}_{apc})-R(\overline{0};{\overline{P}}_{apc})\end{array}$$ E33

where, $$\Delta R_{}^{\omega }({\overline{P}}_{apc})$$ denotes the difference between $$R({\overline{P}}_{\omega },{\overline{P}}_{apc})$$ and

$$R(\overline{0},{\overline{P}}_{apc})$$ E34

And the AF of 3-D SAR could be written as the product of $${\chi }_{}^R(r)$$ and

$${\chi }^p({\overline{P}}_{\omega })$$ E35

$$\chi ({\overline{P}}_{\omega })={\chi }^p({\overline{P}}_{\omega })\cdot {\chi }_{}^R(r-2R(\overline{0},{\overline{P}}_0))$$ E36

Eq (31) indicates that the AF of 3-D SAR could be divided as a range AF and a synthetic aperture AF and analyzed independently. Since the range AF is a sinc function and independent to the antenna phase centre set, $${\chi }^p({\overline{P}}_{\omega })$$ should be paid more attention to.

For further analysis, we approximate $$\Delta R_{}^{\omega }({\overline{P}}_{apc})$$ using the multivariable Taylor’s theorem, and have:

$$\Delta R_{apc}^{\omega }\approx {\overline{P}}_{apc}\cdot {\overline{P}}_{\omega }^T/R(\overline{0},{\overline{P}}_0)$$ E37

where, superscript $$T$$ denotes the transpose operator.

Rewrite $${\overline{P}}_{\omega }$$ in the spherical coordinates as $${\overline{P}}_{\omega }=\gamma \cdot \widehat{\varsigma }$$ , we have:

$$\begin{array}{l}\Delta R_{apc}^{\gamma \cdot \widehat{\varsigma }}=\theta \cdot {\overline{P}}_{apc}\cdot {\widehat{\varsigma }}^T\\ \theta =\gamma /R(\overline{0},{\overline{P}}_0)\end{array}$$ E38

where, $$\gamma$$ denotes the radius of $${\overline{P}}_{\omega }$$ $$\widehat{\varsigma }$$ denotes the direction of $${\overline{P}}_{\omega }$$ $$\theta$$ is the ratio of $$\gamma$$ to

$$R(\overline{0},{\overline{P}}_0)$$ E39

Substituting eq. (33) into eq. (30), we have:

$${\chi }_{\widehat{\varsigma }}^p(\theta )=\frac{1}{{\rm M}}\cdot {\displaystyle \sum _p\exp \left[j\cdot 2\pi \cdot \theta \cdot ({\overline{P}}_{apc}\cdot {\widehat{\varsigma }}^T)/\lambda \right]}$$ E40

Regarding $$({\overline{P}}_{apc}\cdot {\widehat{\varsigma }}^T)/\lambda$$ as a phase function, it is the projection of antenna phase centre set $$P$$ onto the $$\widehat{\varsigma }$$ direction. Using corollary 1, we have the AF in the $$\widehat{\varsigma }$$ direction:

$${\chi }_{\widehat{\varsigma }}^p(\theta )=\frac{1}{{\rm M}}\cdot {\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{e}^{j\cdot 2\pi \cdot \theta \cdot y}{D}_{\widehat{\varsigma }}(y)dy}=\frac{1}{M}\cdot {\widehat{D}}_{\widehat{\varsigma }}(\theta )$$ E41

The physical meaning of eq. (35) is shown in Figure 8. We can obtain the density function $$D_{\widehat{\varsigma }}(i)$$ by counting the number of elements whose projection on the $$\widehat{\varsigma }$$ direction is in the neighbourhood of $$i$$ , and the AF in the $$\widehat{\varsigma }$$ direction is the Fourier transform of $$D_{\widehat{\varsigma }}(i)$$ approximately.

Figure 8.

Explanation on eq. (35), the AF in the $$\widehat{\varsigma }$$ direction is the Fourier transform of $$D_{\widehat{\varsigma }}(i)$$ approximately.

Eq. (35) builds the direct relationship between the antenna phase centre set $$P$$ and the synthetic aperture ambiguity function $${\chi }_{\widehat{\varsigma }}^p(\theta )$$ . Then, the synthetic aperture ambiguity function of typical antenna phase centre sets will be discussed.

Z-shaped trajectory

As a kind of simple continuous trajectory, Z-shaped trajectory is easy to be implemented using 2-D motion control platform and has been used in our experiments, which is shown in Figure 9-a. Figure 9-b and Figure 9-c are its ambiguity functions obtained by simulation and experiment respectively, and its grating lobe is high and dense.

Figure 9.

a) Z-shaped trajectory; (b) 2-D AF by simulation; (c) 2-D AF by experiment; (d) DF vertical to the edge, that two impulses are added on a rectangle function; (e) AF corresponding to Figure (d), whose grating lobes are quite high and dense.

This phenomenon could be explained using eq. (35). Observing Figure 9-a, we find that when the direction is vertical to the edge of the triangle function, all of the elements on one edge are projected on the same point, and there are two impulses are added on the density function (Figure 9-d). Consequently, the sidelobe of the directional AF (which is the Fourier transform of Figure 9-d and shown in Figure 6.e) is high, whose PSLR and ISLR are 7.60dB and -3.28 dB respectively.

Figure 10.

a) DF vertical to the edge with short period; (b) 2-D AF by simulation; (c) AF corresponding to Figure (a), whose grating lobes are sparse

High ISLR and PSLR mean that the sidelobe of strong scatterer will submerge the weak scatterer, and cause measurement error in the 3-D RCS measurement application. This problem could be solved by increasing the periodicity of the Z-shaped trajectory. With the increase of periodicity, the impulses in the density function increase correspondingly (Figure 10 a), and the grating lobe becomes sparse (Figure 10. b and c). Just as the grating lobe problem in the theory of antenna array, when the period of the Z-shaped trajectory is less than 1/2, the grating lobe will be eliminated completely.

Dense square array

For LASAR, its synthetic aperture is a full-element square array, which is shown in Figure 11.a. Since the number of elements is equal for different x, its density function in the x direction is rectangle function, which is shown in Figure 11.b. The directional AF in the x direction is the Fourier transform of the rectangle function, which is a sinc function and shown in Figure 11.e. the black solid line is obtained by numerical simulation, the red dot line is the Fourier transform of Figure 11.b. The peak sidelobe ratio (PSLR) and integrated sidelobe ratio (ISLR) are -11.80 dB and -8.41dB respectively, which is near to the sinc function (-13.30 dB and -10.16 dB respectively).

Then, let’s observe the density function in the diagonal direction. From Figure 11.a, it is obvious that the number of elements increases linearly from one endpoint of the diagonal to the centre, and decreases linearly from the centre to the other endpoint. In consequence, the density function in the diagonal direction is a triangle function, which is shown in Figure 11c. Figure 11 f is the AF in the diagonal direction. The black solid line is obtained by numerical simulation, the red dot line is the Fourier transform of Figure 11.c. The PSLR and ISLR are -22.86 dB and -20.71dB respectively, which is near to the Fourier transform of triangular window (-26.82dB and -22.02 dB respectively)

Figure 11.d is the 2-D AF of dense square array. We find that it has star-shaped AF.

Figure 11.

a) Dense square array; (b) DF in the x direction, which is a rectangle function; (c) DF in the diagonal direction, which is a triangle function; (d) 2-D AF of dense square array, whose sidelobes are distributed in the x and y directions mainly; (e) AF in the x direction, which is a sinc function; (f) AF in the diagonal direction, which is the Fourier transform of triangle function.

Random sampling

For sparse LASAR, its synthetic aperture is a random sampling in the full-element square array, which is shown in Figure 12.a

Figure 12.

a) Random sampling square array; (b) DF in the x direction, which is a rectangle function with noise; (c) DF in the diagonal direction, which is a triangle function with noise; (d) 2-D AF; (e) AF in the x direction, which is similar to the Fourier transform of Figure b; (f) AF in the diagonal direction, which is similar to the Fourier transform of Figure c.

Figure 12. b and c are the density functions in the x and diagonal directions. Comparing with their counterparts of dense square array, we find that its density functions could be considered as the density functions of dense square array modulated by a noise. As a result, its PSLR in the x and diagonal directions (Figure 12.e and 12.f) are similar to those of dense square array. However, since the noise modulated on the density functions increases the high-frequency components, which are corresponded to the far-area sidelobes, its ISLR is higher than that of dense square array.

Figure 12.d is the 2-D AF of uniform distribution sparse array. Comparing with Figure 12.d, we find that its far-area sidelobe is higher than that of dense square array (Note the color bar). Energy leak is the main disadvantage of sparse array. One can improve the PSLR and ISLR by increasing the random sampling number, since its mainlobe energy is proportional to the square of the sampling number, and the sidelobe energy is proportional to the sampling number(Gauss distribution).

3.4. Spatial Resolution

The range resolution of 3-D SAR is produced by the pulse compression technique, and we can obtain the resolution formula directly as:

$${\rho }_R=c/(2B)$$ E42

where, $$c$$ denotes the speed of light, $$B$$ denotes the signal bandwidth.

The other two dimensional resolutions are produced using the array theory, and we can write the resolution formula as:

$$\rho =\lambda /(2\theta )$$ E43

where, $$\lambda$$ denotes the wave length, $$\theta$$ denotes the aperture angle. Remark that the “2” in eq. (37) indicates that the transmitter and receiver moving cooperatively. If the transmitter or receiver is fixed, the resolution formula should approximately be $$\rho =\lambda /\theta$$ .

Note that, the aperture angle $$\theta$$ is influenced by the size of the array, the beam angle of the T/R antenna and the scatterer angle (the angle in which the RCS could be considered as constant.), which are shown in Figure 13.

Figure 13.

Influenced of the array size, aperture angle and scatterer angle on the resolution

The resolution is restricted by the worst factor, i.e.:

$$\rho =\max (\frac{\lambda }{2{\theta }_{ape}},\frac{\lambda R}{2L},\frac{\lambda }{2{\theta }_{sat}})$$ E44

The first two factors could be optimized in the design of 3-D SAR system; the last factor arises from the scattering mechanism and is difficult to be reduced. It also means that we can not improve the resolution of 3-D SAR unlimitedly.

4. Backprojection Method

Backprojection (BP) algorithm is a 3-D SAR imaging algorithm based on the time domain correlation (TDC) technique, which coherently adds the data at the fast-time bin that corresponds to the location of a point for all synthetic aperture locations. The BP algorithm can be considered as the implementation of the definition of ambiguity function, and has been used in E-CSAR data processing.

The input of the BP operator is the raw data, antenna phase centre set and the scatterer’s position; the output is the RCS of the scatterer.

Let $$D_{\text{II}}$$ , $$P$$ and $${\overline{P}}_{\omega }$$ be the raw data after range compression, antenna phase centre set and the scatterer’s position, the BP operator can be expressed as:

$$C[D_{\text{II}},P,{\overline{P}}_{\omega },]\to {\sigma }_{\omega }$$ E45

The implementation of the BP operator $$C[·]$$ is presented in eq. (40):

$$C[D_{\text{II}},P,{\overline{P}}_{uvw},]\triangleq {\displaystyle \sum _p{D}_{\text{II}}\cdot \exp (-j\cdot 2\pi \cdot 2R({\overline{P}}_{uvw},{\overline{P}}_{apc})/\lambda )\cdot {\chi }_{}^R(r-2R({\overline{P}}_{uvw},{\overline{P}}_{apc})}$$ E46

where, $${\overline{P}}_{uvw}$$ denotes the pixel in the image space.

From eq. (40), we know that just like the 2-D BP algorithm, the 3-D BP algorithm can roughly be divided into four steps: range-compression, interpolation, resampling and coherent summation, whose block diagram is shown in Figure 14. Processing the 3-D image region one pixel by one pixel, we can obtain the 3-D RCS distribution finally.

Figure 14.

Block diagram of 3D BP algorithm.

From eq. (40), the computational cost $${\Xi }_C$$ of the single-scatterer compression operator $$C[·]$$ can easily be calculated as:

$${\Xi }_C=M\cdot ({\Xi }_{int}+{\Xi }_{coh})$$ E47

where, $$M$$ denotes the total element number of antenna phase centre set, $${\Xi }_{int}$$ and $${\Xi }_{coh}$$ denote the computational costs of the interpolation operation and the coherent summation operation respectively.

Ignoring the computational cost of the range-compression operation, for a 3-D image region with size L×W×H(pixel³), the total computational cost of 3-D BP algorithm is:

$${\Xi }_{BP}^{}=\text{ L}\cdot \text{W}\cdot \text{H}\cdot {\Xi }_C$$ E48

Compared with 2-D BP, there are two factors that cause the computational cost of 3-D BP algorithm is far larger than that of 2-D BP: the increase of the acquired data and the extension of image region. The former factor can partially be solved using the sparse array technique. The latter factor can be solved using the multiresolution approximation technique, which will be introduced in the following section.

5. Imaging Processing of 3-D SAR

To verify the feasibility of “one-active” LASAR, a series of experiments have been carried out. The typical experiment plan is shown in Figure 15, and includes three parts: “one-active” LASAR, reference points and scene area.

Figure 15.

Experiment plan of 3-D SAR.

The “one-active” LASAR consists of two parts: radar module and motion control module. The radar module is used to transmit the LFM signal and receive the echo from the observation scene. The whole system works on the X-band with signal bandwidth about 120MHz and pulse repetition frequency 20Hz. The motion control module is used to control the transmitter and the receiver moving in a 2-D plane, and synthesize a virtual 2-D antenna array. The motion control module consists of a set of high-precision transfer device with effective length 2mⅹ2m and two high-precision motors, which can compose any continuous 2-D curve. The Z-shaped trajectory with period 1/5, 1/20 and 1/30 are used in the experiment.

Figure 16 is the experiment data of “one-active” LASAR after range-compression. The top black line is the echo of reference point after range compression; the bottom stripped area is the echo of scene area.

Figure 16.

Experiment data after range-compression

Figure 17.

Time-frequency spectrum of single scatterer

To analyze the signal characteristic of “one-active” LASAR, we select one row data in one range-bin, and obtain its time-frequency spectrum by Short-Time Fourier transformation (STFT). Figure 17 is the typical time-frequency spectrum of single scatterer. We find that the time-frequency characteristic of 3-D SAR is more complex than that of traditional SAR (chirp signal), and contains more information on the target.

Figure 18.

a photo of the whole scene

Figure 19.

a photo of fence and lamp Fig. 19-b imaging result of the fence and lamp

Figure 20.

b, c and d, imaging results of the whole scene in the top view, side view and 3-D view respectively.

Figure 21.

a photo of metal spheres Fig. 20-b imaging result of spheres (side view) Fig. 20-c imaging result of spheres (top view)

Figure 18-a is the photo of observation scene, Figure 18-b, c and d show its imaging result in the top view, side view and 3-D view. And we find that the imaging result can depict the main features of the observation scene.

Figure 19-a is the photo of a local region of the observation scene, which contains a metal fence and a street lamp with metal lampstandard. Figure 19-b is its imaging result. Obviously, the imaging result can be divided into three parts, which correspond to the metal fence, corner of the metal fence (far area in the photo) and the metal lampstandard soundly. Especially, according to Figure 18-c, we can read the height of the lamp is about 14m, which matches the real height (14.5m) correctly.

Figure 20-a is the photo of another local region of the observation scene, which contains two copper spheres (a, c) and a stainless steel sphere (b). since its RCS is too low, the copper sphere c is not shown in Figure 20-b and c (could be found in Figure 18-b). Figure 20-b is the imaging result (side view) of spheres a and b. From it, we can read the relative height of spheres a and b is 1.2m, which matches the measurement value (1.09m) soundly. Figure 20-c is the imaging result (top view) of spheres a and b. From it, we can read the relative distance of spheres a and b is 15.8 m, which matches the measurement value (14.96m) soundly.

The above experiment results demonstrate the ability of 3-D SAR in the application of 3-D RCS measurement.

6. Multiresolution Approximation Techniques

Unlike the 2-D SAR, a great deal of 3-D image region contains no scatterer (such as atmosphere) or is shadowed by the other scatterers, and it is not necessary to compress all the pixels in the image region. Based on this idea, multiresolution approximation techniques could be employed to reduce the computational cost of 3-D SAR.

6.1. Scattering Model

Much excellent work has been done on the modeling radar backscatter for both naturally occurring terrain and man-made objects. One of the most popular models is the three-component scattering model developed by Anthony Freeman and Stephen L. Durden. In this model, the scattering mechanism of target is divided into three components, including the rough surface scattering, double-bounce scattering and canopy scattering.

The rough surface scattering component assumes that the backscatter is reciprocal, such as road and bare soil, whose mechanism is illustrated in Figure 21 (left).

The double-bounce scattering component is modeled by scattering from a dihedral corner reflector, where the reflector surface can be made of different dielectric materials, such as a ground-trunk interaction, whose mechanism is illustrated in Figure 21 (middle).

The canopy (volume) scattering component assumes that the radar echo is from a cloud of randomly oriented, very thin, cylinder-like scatterers, such as the forest canopy.

Figure 22.

Scattering mechanisms of three-component model: rough surface scattering (left), double-bounce scattering (middle) and tree scattering model (right) containing both canopy scattering and double-bounce scattering.

In the application of 3-D RCS measurement, the scatterers always concentrate in the local region of the 3-D image space which are shown in Figure 18. In the application of topographical survey, since the thickness of the scattering layer is far smaller than the height of the imaging region, it is convenient to consider the scattering layer as a surface in a low height-resolution level. And one can obtain the scattering layer by searching the neighborhood of the scattering surface.

6.2. Multiresolution approximation

In this subsection, we will introduce the basic concept on the multiresolution wavelet approximation, which is necessary for the design of 3-D LASAR imaging method via multiresolution approximation.

Multiresolution approximation

The multiresolution approximation of $f(t)$ is defined as the orthogonal projection $$P_{{\text{V}}_{\text{i}}}[f]$$ on a multiresolution approximation subspace of $L^2(ℝ)$ . The multiresolution approximation of $L^2(ℝ)$ is a sequence ${\left\{V_j\right\}}_{j\in \mathbb{Z}}$ of closed subspaces of $L^2(ℝ)$ that obeys the following 5 properties:

$$\begin{array}{l}\begin{array}{cc}\forall (j,k)\in {\mathbb{Z}}^2,& f(t)\in V_j\end{array}\iff f(t-2^{j}k)\in V_j\\ \begin{array}{cc}\forall j\in \mathbb{Z},& V_j\end{array}\subset V_{j+1}\\ \begin{array}{cc}\forall j\in \mathbb{Z},& f(t)\in V_j\end{array}\iff f(2t)\in V_{j+1}\\ \underset{j\to -\infty }{\lim }{V}_j=\underset{j=-\infty }{\overset{+\infty }{\cap }}{V}_j=\left\{0\right\}\\ \underset{j\to +\infty }{\lim }{V}_j=\text{closure}\left(\underset{j=-\infty }{\overset{+\infty }{\cup }}{V}_j\right)=L^2({ℝ}^2)\end{array}$$ E49

$$\begin{array}{cc}\forall j\in \mathbb{Z},& V_j\end{array}\subset V_{j+1}$$ E50

where, $j$ denotes the approximation level, and the resolution at level $j$ is

$$2^{-j}$$ E51

Property (43-1) means that $V_j$ is invariant by any translation proportional to the scale $2^j$ . The inclusion (43-2) is a causality property which proves that an approximation at a resolution $2^{j+1}$ contains all the information to compute an approximation at a coarser resolution $2^j$ . Recursive eq. (43-3) specifies the relationship between approximation subspaces. The property (43-4) implies that we lost all the details of $f$ when the level goes to $-\infty$ ; on the other hand, when the level goes $+\infty$ , property (43-5) imposes that the signal approximation converges to the original signal.

According to the approximation theory, the basis of $V_j$ can be generated by dilating and translating a scaling function $$\varphi (t)$$ :

$${\varphi }_{j,n}(t)=\frac{1}{\sqrt{2j}}\varphi (\frac{t-n}{2^j})$$ E52

Thus, the multiresolution approximation $${\tilde{f}}_j(t)$$ of $f(t)$ at level $j$ can be calculated as:

$${\tilde{f}}_j(t)={\displaystyle \sum _{n=-\infty }^{+\infty }f,{\varphi }_{j,n}(t)\cdot {\varphi }_{j,n}(t)}$$ E53

where, $·$ denotes the inner product.

Conjugate mirror filter

To ensure that the multiresolution approximation can be conducted recursively, it is necessary to analyze the relationship between the approximations of $$f(t)$$ at level i and i+1.

The multiresolution causality property (42-2) imposes that $V_j\subset V_{j+1}$ . Since $${\varphi }_{j+\mathrm{1,}n}(t)$$ is an orthonormal basis of $V_{j+1}$ , we can decompose $${\varphi }_{j\mathrm{,0}}(t)$$ as:

$${\varphi }_{j\mathrm{,0}}(t)={\displaystyle \sum _{n=-\infty }^{+\infty }h[n]\cdot {\varphi }_{j+\mathrm{1,}n}(t)}$$ E54

With:

$$h(n)={\varphi }_{j\mathrm{,0}}(t),{\varphi }_{j+\mathrm{1,}n}(t)$$ E55

where, eq. (47) is called two-scale relation, $$h(n)$$ denotes the conjugate mirror filter corresponding to the scaling function $$\varphi (t)$$ , provide that:

$$\widehat{\varphi }(\omega )={\displaystyle \prod _{p=1}^{+\infty }\frac{\widehat{h}(2^{-p}\omega )}{\sqrt{2}}\widehat{\varphi }(0)}$$ E56

According to the fast orthogonal wavelet transform developed by Stephane. G. Mallat, the reconstruction of $${\tilde{f}}_{j+1}(t)$$ can be implemented by a two-channel multirate filter bank:

$$\begin{array}{l}{\tilde{f}}_{j+1}(n)={\stackrel{\vee }{f}}_j(n)\ast h(n)+{\stackrel{\vee }{d}}_j(n)\ast g(n)\\ \begin{array}{cc}{\stackrel{\vee }{f}}_j(n)=\{\begin{array}{c}\begin{array}{cc}{\tilde{f}}_j(n)& n=2p\end{array}\\ \begin{array}{cc}0& n=2p+1\end{array}\end{array}& p=\mathrm{0,1,...}\end{array}\\ \begin{array}{cc}{\stackrel{\vee }{d}}_j(n)=\{\begin{array}{c}\begin{array}{cc}{\tilde{d}}_j(n)& n=2p\end{array}\\ \begin{array}{cc}0& n=2p+1\end{array}\end{array}& p=\mathrm{0,1}\end{array}\end{array}$$ E57

where, $${\stackrel{\vee }{d}}_j(n)$$ denotes detail coefficients, $$g(n)$$ denotes the high-pass filter corresponding to $$h(n)$$ , satisfying:

$$\widehat{g}(\frac{\omega }{2})=e^{-j\omega }{\widehat{h}}^{*}(\omega +\pi )$$ E58

Figure 23.

Diagrams on wavelet interpolation, the left top, right and left bottom are the diagrams of 1-D fast wavelet transform, 2-D fast wavelet transform and 2-D wavelet interpolation respectively

The diagram of 1-D wavelet reconstruction is shown in Figure 22 (left top), where, upsample operator “↑2” inserts zeros at odd-indexed elements.

2-D multiresolution surface approximation

The properties of two-dimensional wavelet are essentially the same as in one dimension. A separable two-dimensional wavelet transform can be factored into one-dimensional wavelet transforms along the rows and columns.

Assume that the conjugate mirror filters corresponding to the one-dimensional wavelet transform are denoted as $$h(n)$$ and $$g(n)$$ . According to the fast two-dimensional wavelet transform, the reconstruction of a two-dimensional function $${\tilde{f}}_{j+1}(n_1,n_2)$$ can be implemented by the following equation:

$$\begin{array}{l}{\tilde{f}}_{j+1}(n_1,n_2)={\stackrel{\vee }{f}}_j(n_1,n_2)\ast [h(n_1)h(n_2)]+\stackrel{\vee }{d_j^1}(n_1,n_2)\ast [h(n_1)g(n_2)]\\ +\stackrel{\vee }{d_j^2}(n_1,n_2)\ast [g(n_1)h(n_2)]+\stackrel{\vee }{d_j^3}(n_1,n_2)\ast [g(n_1)g(n_2)]\end{array}$$ E59

where, $$\stackrel{\vee }{d_j^1}(n_1,n_2)$$ $$\stackrel{\vee }{d_j^2}(n_1,n_2)$$ and $$\stackrel{\vee }{d_j^3}(n_1,n_2)$$ denotes

the detail coefficients matrices. And the corresponding diagram is shown in Figure 22 (right).

In the application of surface prediction, since there is no information on the detail coefficients matrices, we just consider that all of them are zero matrices, and the reconstruction formula can be simplified as:

$${\tilde{f}}_{j+1}(n_1,n_2)={\stackrel{\vee }{f}}_j(n_1,n_2)\ast [h(n_1)h(n_2)]$$ E60

And the corresponding diagram is shown in Figure 22 (left bottom).

Typical conjugate mirror filters

According to the discussion in above section, we conclude that the 2-D surface multiresolution prediction is specified by the conjugate mirror filter $$h(n)$$ . In this subsection, we present some typical conjugate mirror filter as follows:

Shannon conjugate mirror filter:

$$\widehat{h}(\omega )=\{\begin{array}{cc}\begin{array}{c}\sqrt{2}\\ 0\end{array}& \begin{array}{c}\omega \in [-\pi /\mathrm{2,}\pi /2]\\ others\end{array}\end{array}$$ E61

Meyer conjugate mirror filter:

$$\widehat{h}(\omega )=\{\begin{array}{c}\begin{array}{cc}\sqrt{2}& \omega \in [-\pi /\mathrm{3,}\pi /3]\end{array}\\ \begin{array}{cc}0& \omega \in [-\pi ,-2\pi /3]\cup [-2\pi /\mathrm{3,}\pi ]\end{array}\end{array}$$ E62

DB conjugate mirror filter:

$${\left|\widehat{h}(\omega )\right|}^2=2{\cos }^{2p}(\frac{\omega }{2})P({\sin }^2(\frac{\omega }{2}))$$ E63

where, polynomial $P(y)$ satisfies:

$${(1-y)}^{p}P(y)+y^{p}P(1-y)=1$$ E64

Spline biorthogonal conjugate mirror filter

$$\widehat{h}(\omega )=\sqrt{2}\exp (\frac{-j\epsilon \omega }{2}){\cos }^p(\frac{\omega }{2})$$ E65

where, $$\epsilon =0$$ if $p$ is even and $$\epsilon =1$$ for $p$ is odd.

For the different applications, the conjugate mirror filter might affect the estimation error of the prediction operator, it is sensible to select the conjugate mirror filter according to the application.

6.3. Surface Tracing Technique

In the application of topographical survey, the scatterers combine a surface in the 3-D image space, we can trace the surface and focus those scatterers near it via specific searching method. Thus, the 3-D SAR imaging processing can be reduced to a 2-D imaging problem, and the computational cost will be reduced greatly. Those methods based on the above idea are named as surface-tracing-based 3-D imaging method (STB 3-D imaging method).

Principle and steps

Figure 24.

Principle of STB 3-D SAR imaging technique.

The principle of STB 3-D imaging method is illustrated in Figure 23. Assume that the line of radar sight (i.e. LOS) is parallel to the elevation (z) direction, and the 3-D scattering surface $$\Upsilon$$ in the observation scene can be expressed as following:

$$\begin{array}{cc}\Upsilon =\{(x,y,z)/z=h(x,y)& (x,y)\in \Omega \subset {ℝ}^2\end{array}\}$$ E66

where, $$(x,y)$$ is called note, $$\Omega$$ is the note set, $$h$$ denotes the elevation function.

For a given note $$(x_0,y_0)$$ , the RCS distribution in z direction is a quasi-impulse function which is shown in Figure 23 (right). The index of the maximum is the elevation of note $$(x_0,y_0)$$ . Searching the maximum of every note $$(x,y)$$ in $$\Omega$$ , one can reconstruct the DEM of the scene eventually. Generally, the steps of STB 3-D imaging method are stated as follows:

Step-1 Initiation:

Compress in a subset $${\Omega }_0$$ of the note set $$\Omega$$ by 3-D BP imaging method, find the maximum and the associated index of every note, and obtain the initial subsurface $${\Upsilon }_0$$ ;

$$\begin{array}{cc}{\Upsilon }_0=\{(x,y,z)/z=h(x,y)& (x,y)\in {\Omega }_0\subset \Omega \end{array}\}$$ E67

Step-2 Prediction:

Expand the note set $${\Omega }_0$$ to $${\Omega }_1$$ ( $${\Omega }_0\subset {\Omega }_1\subseteq \Omega$$ predicate the surface on $${\Omega }_1$$ using the known surface $${\Upsilon }_0$$ by a surface prediction operator $$P[·]$$ , and obtain the estimation of the surface on $${\Omega }_1$$ , denoted as $${\widehat{\Upsilon }}_1$$

$$P[{\Upsilon }_0]\to {\widehat{\Upsilon }}_1$$ E68

Step-3 Searching:

Search in the neighborhood of $${\widehat{\Upsilon }}_1$$ , and obtain the actual surface $${\Upsilon }_1$$

$$S[{\widehat{\Upsilon }}_1]\to {\Upsilon }_1$$ E69

Step-4 : Recursion:

Replace $${\Omega }_0$$ and $${\Upsilon }_0$$ in step 2 by $${\Omega }_1$$ and $${\Upsilon }_1$$ , and repeat the step 2-4 until obtaining the 3-D surface $$\Upsilon$$

where, the most important steps are predication and searching, which will be discussed the in the rest of this section.

Surface prediction operator

The aim of prediction operator $$P[·]$$ is to estimate the likely elevation of the surface using the known subsurface, which is necessary for the searching operator.

The input of the prediction operator is the known subsurface on a note set $${\Omega }_i$$ ; the output is the estimation of the subsurface on the note set $${\Omega }_{i+1}$$

$${\Omega }_i\subset {\Omega }_{i+1}$$ E70

Let $${\Upsilon }_i$$ and $${\Upsilon }_{i+1}$$ be the subsurfaces on the note sets $${\Omega }_i$$ and $${\Omega }_{i+1}$$ respectively, the prediction operator $$P[·]$$ can be expressed as:

$$P[{\Upsilon }_i]\to {\widehat{\Upsilon }}_{i+1}$$ E71

And the prediction error can be defined as:

$$e=z-\widehat{z}$$ E72

where, $$\widehat{z}$$ and $$z$$ denote the predicated elevation and the actual elevation at given note respectively.

Mathematically, the prediction operator can be implemented using multivariate interpolation technique. According to the interpolation method, the prediction operator includes the polynomial prediction operator, ridge prediction operator, spline prediction operator and multiresolution (wavelet) prediction operator, etc.

In the viewpoint of predication strategy, the prediction operator can roughly be classified as two classes: local prediction operator and multiresolution prediction operator. The former operator starts from the local region of the note set and expands the known region from the edge, which is shown in Figure 24-a. The latter one starts from a coarse resolution scene and improves the resolution of the scene recursively, which is shown in Figure 24-b.

Figure 25.

Local prediction (left) and multiresolution prediction (right)

In the case of high-resolution 3D SAR imaging, since the elevations of the neighbor notes are affected by the fluctuation of ground greatly, the prediction error of the local prediction operator is larger than that of the multiresolution prediction operator. On the other hand, for an N×N scene, the recursion times of the local prediction operator are N, and the recursion times of the multiresolution prediction operator are ${\log }_2(\text{N})$ . Less recursion times always mean less interpolation operation and less computational cost.

Searching operator

The searching operator $$S[·]$$ is to find out the maximum and the associated index at every note.

The input of the searching operator includes the raw data, note, estimated elevation, threshold; the output includes the actual elevation and the associated RCS.

Let $$D$$ , $(x_0,y_0)$ $\widehat{Z}$ $\Theta$ $Z$ and $\sigma$ be the raw data, note, estimated elevation, threshold, actual elevation and the associated RCS respectively, the searching operator can be expressed as:

$$S[D,(x_0,y_0),\widehat{Z},\Theta ]\to [Z,\sigma ]$$ E73

The prediction operator can be implemented using the following equation:

$$\{\begin{array}{c}\sigma (x_0,y_0,z_{\max })\ge \sigma (x_0,y_0,z_{\max }-1)\\ \sigma (x_0,y_0,z_{\max })\ge \sigma (x_0,y_0,z_{\max }+1)\\ \sigma (x_0,y_0,z_{\max })\Theta \end{array}$$ E74

It indicates that the maximum is larger than the adjacent pixels and the detection threshold.

The detection threshold $$\Theta$$ can be selected based on the constant false alarm rate (CFAR) criteria. Assume that the RCS and the noise both obey the normal distribution, $$\Theta$$ can be calculated as:

$$\Theta ={\nu }_0\cdot erf{c}^{-1}(P_{false})+{\mu }_0$$ E75

where, $${\mu }_0$$ and $${\nu }_0$$ denote the mean and variance of the signal respectively, and can be obtained in the step of initiation, $$erf{c}^{-1}(·)$$ denotes the inverse complementary error function.

Numerical results

In this subsection, some numerical experiments are conducted to demonstrate the procedures of the STB 3-D BP algorithm.

Figure 26.

a is the terrain used in the numerical experiments. Figure 25-b, c, d and e are the imaging results of the STB 3-D BP algorithm with 1^st, 2^nd, 3^rd, and 4^th iterations respectively. From them, we find that the STB 3-D BP algorithm can reconstruct the DEM of the scene correctly. By further analysis (Shi Jun, Zhang Xiaoling, Jianyu Yang;, Wang Yinbo, 2008), we find that, in general case, the total times to call the compress operator of the STB 3-D imaging algorithm is a dozen of times larger than that of 2-D BP algorithm generally, and is far smaller than that of the 3-D BP algorithm that is about several hundreds times larger than 2-D BP (determined by the height of the 3-D imaging scene).

6.4. Subaperture Approximation Technique

The STB based multiresolution approximation technique can only be used for the topographical survey application, which can reduce the size of image space. In fact, the multiresolution approximation technique can also reduce the number of antenna elements that need to be processed, since only a subaperture is necessary to obtain a low-resolution image. Based on this idea, we can obtain a low-resolution image using a subaperture, detect the interested regions with scatterers, and process the interested regions with a larger subaperture iteratively until obtain the fine-resolution image. This kind of multiresolution approximation technique is named as subaperture approximation technique, which can be used in both 3-D RCS measurement and topographical survey applications. Limited by the length of this chapter, this topic will be discussed hereafter.

7. Summary

Using synthetic aperture technique, we can obtain the 3-D RCS distribution of target. The precondition of synthetic aperture technique is that the target’s RCS does not vary with time in one aperture. The synthetic aperture can be implemented mechanically (such as CSAR, E-CSAR and curve SAR), or electrically (such as linear array SAR). By moving HRR radar in 2-D plane using high precision motion control platform, we can build a low-cost 3-D RCS measurement device. The linear array SAR with MIMO technique might be the most feasible 3-D SAR system for the topographical survey application, though there are still some problems, such as, the balance between the length of linear array and the cross-track resolution and the compensation of motion measurement error.

Ambiguous function (AF) of 3-D SAR is the product of the range AF and the synthetic aperture AF, which can be analyzed independently. The range AF is a sinc function without any window function; the synthetic aperture AF could be analyzed using the theory of array antenna. The resolution in the synthetic aperture direction(s) is restricted by the size of the array, the beam angle of the T/R antenna and the scatterer angle. Since the RCS varies in different elevation angle and azimuthal angle, we can not improve the resolution of 3-D SAR unlimitedly.

Backprojection method can be employed in 3-D SAR imaging processing. Its disadvantage is the high computional cost. Unlike the 2-D microwave image, the scatterers always concentrate in the local regions of 3-D image space. Based on this feature, the multiresolution approximation technique could be employed in imaging processing, which can reduce the computational cost greatly.