1 Introduction

Microwave imaging technique for evaluation of hidden objects in a structure has been widely applied in the areas of breast tumor imaging,Ground Penetrating Radar (GPR),through-wall imaging and concealed weapon detection^{[1, 2, 3]}. It has advantages over other Non-Destructive Evaluation (NDE) methods regarding low cost,portability, good penetration in nonmetallic materials, good resolution,and contactless detection. The imaging algorithms for microwave imaging technique can be classified into quantitative and qualitative methods^{[1, 4]}. Quantitative methods retrieve the shapes and the electromagnetic parameters of unknown targets by solving a nonlinear inverse scattering problem. The disadvantage of quantitative methods is time consuming.

Qualitative methods usually make some assumptions and mathematical approximations for lower computational complexity^[5]. The classical method is Frequency-Wavenumber (F-K) migration, also known as Stolt migration^{[3, 6]}. F-K migration can be computed efficiently by Fast Fourier Transform (FFT). The algorithm works well for homogeneous background medium. However,in the case of layered medium,the image becomes unfocused and the hidden objects are incorrectly located. Another classical qualitative method involving Back-Projection (BP) algorithm is based on delay-and-sum beamforming. The image is reconstructed by summation of the backscattered data with proper phase compensation. The phase of the scattered signal in lossless medium is proportional to the round-trip distance in homogenous medium. For layered medium,the electromagnetic wave refracts across the discontinuity at the interface. The phase for two-layered medium can be solved numerically^{[7, 8]}. But for three or more layers,the computational complexity will increase drastically. An efficient tomographic algorithm was proposed for 2-D microwave imaging through layered media in Ref. ^[9]. In this paper, an algorithm based on Phase Shift Migration (PSM) is presented for microwave imaging of layered medium.

PSM can handle velocity variations with depth^[10]. F-K migration can be regarded as a special case of PSM. PSM originates from seismic engineering and geophysics^{[10, 11]}. It recursively extrapolates the phase shift along depth in steps and uses the result of each step as input to the next iteration.

In this paper,the phase shift for layered medium is derived without recursion using the phase matching at the interface. So the image can be reconstructed in a fairly straight forward way. We also compensate the discontinuity of the interface using Fresnel transmission coefficients. The analysis on backscattered transfer function for layered medium makes clear to understand the assumptions and mathematical approximations in PSM.

This paper is organized as follows. Section 2 gives the derivation of PSM for layered medium. Section 3 analyzes the backscattered transfer function for layered medium and achieves the identical imaging equation by making some assumptions and approximations. In Section 4,numerical and experimental results are provided to show the effectiveness of PSM for layered medium. Finally,Section 5 gives a brief conclusion. 2 PSM for Layered Medium

PSM method was first proposed by Gazdag in 1978^[10]. It is based on the concept of exploding source model^[11],which assumes that a fictitious point source explodes at a reference time of $t=0$ at the target locationand radiates EM wave to the receiver. The velocity in the medium is replaced by half of the true propagation velocity for exploding source model.

Assuming that a point object is located at $r(x,y,z)$ and an antenna is located at $r'(x',y',z')$. The antenna illuminates the target with angular frequency . The antenna synthesizes a rectangular aperture on a plane $z=z'$ parallel to xy. If we make $z'=0$,EM field $f (x,y,z,t)$ can be described by the following 3-D scalar wave equation for homogeneous medium:

$\left( {{{{\partial ^2}} \over {\partial {x^2}}} + {{{\partial ^2}} \over {\partial {y^2}}} + {{{\partial ^2}} \over {\partial {z^2}}} - {1 \over {v_m^{^2}}}{{{\partial ^2}} \over {\partial {t^2}}}} \right)f\left( {x,y,z,t} \right) = 0$

(1)

where ${v_m} = v/2$,$v$ is the propagation velocity of EM wave in the medium.

Applying 3-D Fourier transform to Eq. (1) over $x,y$,and $z$,we get

$\left( {k_x^2 + k_y^2 + {{{\partial ^2}} \over {\partial {z^2}}} - {k^2}} \right)F\left( {{k_x},{k_y},z,w} \right){\rm{ = }}0$

(2)

When upgoing EM wave is considered only,the solution of Eq. (2) can be expressed as:

$F\left( {{k_x},{k_y},z,w} \right){\rm{ = }}F\left( {{k_x},{k_y},z{\rm{ = }}0,w} \right){{\rm{e}}^{jk{z^z}}}$

(3)

where ${k_x},{k_y}$,and ${k_z}$ are the wavenumbers along $x$, $y$,and $z$ directions,respectively. Note that in Eq. (3)${k_z} = \sqrt {{k^2} - k_x^2 - k_y^2} $ ,where ${k^2}=4{\left( {w/v} \right)^2}$. $k$ is the magnitude of the wavenumber vector $k$.

Doing 3-D inverse Fourier transform to $F({k_x},{k_y},z,w)$ over $x,y$,and $t$,we get:

$\eqalign{ & f\left( {x,y,z,t} \right) = \int\!\!\!\int\!\!\!\int F \left( {{k_x},{k_y},z{\rm{ = }}0,w} \right) \cr & \qquad \qquad \qquad \cdot {{\rm{e}}^{jk{z^z}}}{{\rm{e}}^{jk{x^x}}}{{\rm{e}}^{jk{y^y}}}{{\rm{e}}^{j{w^t}}}d|{k_x}d{k_y}dw \cr} $

(4)

According to the imaging principle^[11],the wavefront shape at $t=0$ corresponds to the reflector shape. Let $t=0$,we have the reflectivity at $r(x,y,z)$:

$\eqalign{ & I\left( {x,y,z} \right) = f\left( {x,y,z,t = 0} \right) \cr & \, \qquad \qquad = \int {FT_{xy}^{ - 1}\left\{ {F\left( {{k_x},{k_y},z{\rm{ = }}0,w} \right){{\rm{e}}^{jk{z^z}}}} \right\}} dw \cr} $

(5)

$F({k_x},{k_y},z=0,w)$ can be obtained by applying 2-D spatial Fourier transform to $f (x,y,z=0,w)$,which represents the collected data on the antenna plane. The data can be collected with a time-domain system as well. In this case,$f (x,y,z=0,w)$ can be simply obtained by applying Fourier transform to $f (x,y,z=0,t)$ overtime.

The multiplication by ${{\rm{e}}^{jk{z^z}}}$ in Eq. (5) represents phase shift factor. It means that the extrapolation along $z$ direction in the frequencywavenumber domain is a phase shift operation.

Consider a point object embedded in an Nlayer medium,as shown in Fig. 1. The electromagnetic properties vary with depth only,$i.e$. the z-direction. For layered medium,the $z$ component of the wavenumber vector is not constant:

${k_{iz}} = \sqrt {4{{\left( {w/{v_i}} \right)}^2} - k_x^2 - k_y^2} ,i = 1,2, \cdot \cdot \cdot ,N$

(6)

We can generally define the phase shift factor f or the point target in layer $n$ as:

${P_{\left( z \right)}} = {{\rm{e}}^{j{k_1}{z^{d1}}}}\left( {\prod\limits_{i = 1}^{n - 1} {{{\rm{e}}^{j{k_1}z\left( {{d_i} - {d_{i - 1} \ \ }} \right)}}} } \right){{\rm{e}}^{j{k_n}{z^{\left( {z - {d_{n - 1} \ \ }} \right)}}}}$

(7)

So the reflectivity at $r(x,y,z)$ in region $n$ can be written as:

$I\left( {x,y,z} \right) = \int {FT_{xy}^{ - 1}\left[ {F\left( {{k_x},{k_y},z{\rm{ = }}0,w} \right){{\rm{e}}^{j{k_1}{z^{d1}}}} \\ \qquad \qquad \quad \cdot \left( {\prod\limits_{i = 1}^{n - 1} {{{\rm{e}}^{j{k_i}z\left( {{d_i} - {d_{i - 1} \ \ }} \right)}}} } \right){{\rm{e}}^{j{k_n}{z^{\left( {z - {d_{n - 1 \ \ } }} \right)}}}}} \ \ \right]} dw$

(8)

Eq. (8) does not consider the discontinuity between layers. One way to consider the discontinuity between layers is to use Fresnel transmission coefficients:

$I\left( {x,y,z} \right) = \int {FT} _{xy}^{ - 1}\left[ {F\left( {{k_x},{k_y},z = 0,w} \right) \\ \qquad \qquad \quad \cdot {{{e^{j{k_{1z}}^{_{^{d1}}}}}\left( {\prod\limits_{i = 1}^{n - 1} \ {{e^{j{k_{_{_{iz\left( {{d_{_i}} - {d_{_i}} - 1} \ \ \right)}}}}}}} } \ \ \ \ \right) \ {{\rm{e}}^{j{k_{nz\left( {z - {d_n} - 1} \right) \ \ }}}}} \over {\prod\limits_{i = 1}^{n - 1} {{T_{i,i + 1}} \ \ {T_{i + 1,i}}} }}} \ \ \right]dw$

(9)

where ${T_{i,i + 1}} =\left( {{\rm{2}} {u_i}{\rm{ + 1}}{k_{iz}}} \right)/\left( {{u_i}{\rm{ + 1}}|{k_{iz}} + {u_i}{k_{i + 1,z}}} \right)$ , ${T_{i + 1,i}}=2{u_i}{k_{i + 1,z}}/\left( {{u_i}{k_{i + 1,z}} + {u_{i + 1}}{k_{i,z}}} \right)$.

Generally,when the antenna plane is $z=z'$, the reflectivity at $r(x,y,z)$ can be expressed as:

$I\left( {x,y,z} \right) = \int {FT} _{xy}^{ - 1}\left[ {F\left( {{k_x},{k_y},z = z',w} \right) \\ \qquad \qquad \quad \cdot {{{e^{j{k_{1z}}^{\left( {_{^{d1}} - z'} \right)}}}\left( {\prod\limits_{i = 1}^{n - 1} {{e^{j{k_{^{_{_{_{iz\left( {{d_{_i}} - {d_{_i}} - 1} \ \ \right)}}}}}}}}} } \ \ \ \ \right){{\rm{e}}^{j{k_{nz\left( {z - {d_n} - 1} \ \ \right)}}}}} \over {\prod\limits_{i = 1}^{n - 1} {{T_{i,i + 1}} \ \ {T_{i + 1,i}}} }}} \ \ \right]dw$

(10)

3 Analysis on PSM with Backscattered Transfer Function

The backscattered transfer function should be figured out to calculate the scattering field. Scattering problem can be typically formulated as the superposition of the scattering from a set of points,assuming that no interaction occurs between points. This leads to the following model:

$f\left( {r',w} \right) = \mathop{\int\!\!\!\int\!\!\!\int}\limits_{\kern-5.5pt V} I \left( r \right){G_{GT}}\left( {r,r',w} \right)dr$

(11)

where $f (r',w)$ is the collected data,$I(r)$ is the reflectivity of the point target. $GRT(r,r',w)$ indicates the backscattered transfer function,which accounts for the two-way propagation process:

$G_{RT}(r,r',w) = G_1(r,r',w)G_2(r,r',w)$

(12)

$G_1(r,r',w)$ is the Green’s function for the wave propagation process from the transmitting antenna to the point target. $G_2(r,r',w)$ is the Green’s function for the wave propagation process from the point target to the receiving antenna. The reflectivity function can be reconstructed as:

$I\left( r \right) = \int {dw\int\!\!\!\int\limits_{S'} f } \left( {r',w} \right)G_{_{GT}}^{ - 1}\left( {r,r',w} \right)dr'$

(13)

The Green’s functions for layered medium can be written as^{[12, 13]}:

$\eqalign{ & {G_1}\left( {r,r',w} \right) = {{ - j} \over {8{\pi ^2}}}\int\!\!\!\int {{{{e^{ - j{k_x}\left( {x - x'} \right) - j{k_y}\left( {y - y'} \right)}}} \over {{k_{1z}}}}} \cr & \qquad \qquad \qquad \cdot {F_ - }\left( {z,z'} \right)d{k_x}d{k_y} \cr} $

(14)

$\eqalign{ & {G_2}\left( {r',r,w} \right) = {{ - j} \over {8{\pi ^2}}}\int\!\!\!\int {{{{e^{ - j{k_x}\left( {x - x'} \right) - j{k_y}\left( {y - y'} \right)}}} \over {{k_{1z}}}}} \cr & \qquad \qquad \qquad \cdot {F_ + }\left( {z,z'} \right)d{k_x}d{k_y} \cr} $

(15)

Neglecting the upgoing wave below region n, ${F_ - }(z,z')$ and ${F_ + }(z',z)$ can be simplified as:

${F_ - }\left( {z,z'} \right) = {{\tilde T}_{1n}}{e^{^{ - j{k_{1z\left( {{d_1} - z'} \right)}}}}}{e^{^{ - j{k_{nz\left( {z - {d_n} - 1} \right)}}}}}$

(16)

${F_ + }\left( {z,z'} \right) = {{\tilde T}_{n1}}{e^{^{ - j{k_{1z\left( {{d_1} - z'} \right)}}}}}{e^{^{ - j{k_{nz\left( {z - {d_n} - 1} \right)}}}}}$

(17)

where ${k_{iz}} = \sqrt {{{\left( {w/{v_i}} \right)}^2} - k_x^2 - k_y^2} $ . ${{\tilde T}_{1n}}$ is the generalized transmission coefficient for the layered medium from layer 1 to layer $n$. ${{\tilde T}_{n1}}$ is the generalized transmission coefficient for the layered medium from layer n to layer 1. ${{\tilde T}_{1n}}$ and ${{\tilde T}_{n1}}$ are defined a s^[12]:

${\tilde T_{1n}} = \frac{{\prod\limits_{i = 1}^{n - 1} {{e^{ - j{k_{iz}}\left( {{d_i} - {d_{i - 1}}} \right)}}} {T_{i,i + 1}}}}{{1 - {R_{i + 1,i}}{{\tilde R}_{i + 1,i + 2}}{e^{ - j2{k_{i + 1,z}}\left( {{d_{i + 1}} - {d_i}} \right)}}}}$

(18)

${\tilde T_{n1}} = \frac{{\prod\limits_{i = 1}^{n - 1} {{e^{ - j{k_{iz}}\left( {{d_i} - {d_{i - 1}}} \right)}}{T_{i + 1,i}}} }}{{1 - {R_{i,i + 1}}{{\tilde R}_{i,i - 1}}{e^{ - j2{k_{iz}}\left( {{d_i} - {d_{i - 1}}} \right)}}}}$

(19)

The denominator of ${{\tilde T}_{1n}}$ and ${{\tilde T}_{n1}}$ can be physically interpreted as a consequence of the multiple reflections of the layered medium. For the sake of simplicity,the multiple reflections are not considered. Ignoring the spreading term and the constant term,the backscattered transfer function can be simplified as:

$\eqalign{ & {G_{RT}}(r,r',w) \cr & \qquad = \left( {\prod\limits_{i = 1}^{n - 1} {{T_{i,i + 1}}{T_{i + 1,i}}} } \right) \cr & \qquad \quad \cdot {\left( {\int\!\!\!\int {{e^{ - j{k_x}\left( {x - x'} \right) - j{k_y}\left( {y - y'} \right)}}\left( {\prod\limits_{i = 1}^{n - 1} {{e^{ - j{k_{iz}}\left( {{d_i} - {d_i} - 1} \right)}}} } \right) \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \\ \cdot } {e^{^{ - j{k_{1z\left( {{d_1} - z'} \right)}}}}}{e^{^{ - j{k_{nz\left( {z - {d_n} - 1} \right)}}}}}d{k_x}d{k_y}} \right)^2} \cr} $

(20)

When the target is in the far field of the antenna, the integral in Eq. (20) can be evaluated asymptotically with the method of stationary phase^{[12, 14, 15]}. The backscattered transfer function can be approximated by:

$\eqalign{ & {G_{RT}}(r,r',w) \cr & \qquad = \left( {\prod\limits_{i = 1}^{n - 1} {{T_{i,i + 1}}{T_{i + 1,i}}} } \right) \cr & \qquad \quad \cdot \int\!\!\!\int {{e^{ - j{k_x}\left( {x - x'} \right) - j{k_y}\left( {y - y'} \right)}}\left( {\prod\limits_{i = 1}^{n - 1} {{e^{ - j{k_{iz}}\left( {{d_i} - {d_i} - 1} \right)}}} } \right)} \cr & \qquad \quad \cdot {e^{^{ - j{k_{1z\left( {{d_1} - z'} \right)}}}}}{e^{^{ - j{k_{nz\left( {z - {d_n} - 1} \right)}}}}}d{k_x}d{k_y} \cr} $

(21)

where $k_{iz}$ is expressed as Eq. (6). So the reflectivity image can be reconstructed as:

$I\left( {x,y,z} \right) = \int {FT} _{xy}^{ - 1}\left[ {F{T_{xy}}\left[ {f\left( {x',y',z',w} \right)} \right] \\ \qquad \qquad \quad \cdot {{{e^{j{k_{1z}}^{\left( {_{^{d1}} - z'} \right)}}}\left( {\prod\limits_{i = 1}^{n - 1} {{e^{j{k_{^{_{_{_{iz\left( {{d_{_i}} - {d_{_i}} - 1} \ \right)}}}}}}}}} } \ \ \ \right){{\rm{e}}^{j{k_{nz\left( {z - {d_n} - 1} \right)}}}}} \over {\prod\limits_{i = 1}^{n - 1} {{T_{i,i + 1}} \ \ \ {T_{i + 1,i}}} }}} \ \ \ \right]dw$

(22)

which is identical to PSM for layered medium given by Eq. (10).

The PSM for layered medium can be viewed as an approximate inversion of the inverse scattering problem. The assumption of the proposed algorithm can be summarized as follows. (1) Under first-order Born approximation,multiple scattering within the target is ignored. (2) For the object buried in layer n,the upgoing wave below region $n$ is ignored. That means the generalized reflection coefficient for the layered medium below region $n$ is set to zero. (3) The multiple reflections in layered medium are not considered. (4) The propagation decay is ignored. (5) The imaging area is in the far field of the antenna. 4 Results

In order to show the high computational efficiency and accurate image formation of the algorithm, both numerical simulation and experimental measurement are demonstrated in this section. Matlab is used to perform the simulation and imaging on a 64 bit PC with 16 GB RAM and four-core 2.5 GHz CPU. 4.1 Numerical results

Firstly,we present a numerical example for the imaging of the targets buried in a layered medium shown in Fig. 2. The antenna is placed on the plane $z=0$ with a standoff distance of 2 cm. The area is scanned along x direction from $x=-15$ cm to 15 cm with a step size of 0.5 cm. The first layer is wood with dielectric constant $ε_1=2$. The second layered is concrete with dielectric constant $ε_2=6$. Thicknesses of wood and concrete are chosen to be 5 cm and 7.5 cm,respectively. The buried objects are metallic cylinders with 1.5 cm radius located at (-5 cm,5 cm) and (5 cm,11 cm). The measured data is generated in time domain using 2-D Finite Difference Time Domain (FDTD) method and is transformed to frequency domain by Fourier transform. The frequency range covers from 3 GHz to 8 GHz with a step size of 25 MHz.

The reflected signal from interfaces of layers is relatively stronger compared with the scattered signal from targets. To clearly show the image of targets,we can estimate the background by averaging collected data over the scanning area and subtract the average from the collected data.

$f'\left( {x,w} \right) = f\left( {x,w} \right) - \sum\limits_x {{{f\left( {x,w} \right)} \over {{N_x}}}} $

(23)

The collected data was processed by the proposed algorithm and $F-K$ migration for comparison. For $F-K$ migration,the effective dielectric constant can be approximately obtained by Root Mean Square (RMS)^[16]:

$\varepsilon_{RMS} = {\left( {{{\sum\limits_{i = 1}^N {\varepsilon _i^2\Delta {t_i}} } \over {\sum\limits_{i = 1}^N {\Delta {t_i}} }}} \right)^{{1 \over 2}}}$

(24)

where $Δt_i$ is the round-trip delay of the $i-th$ layer. The RMS dielectric constant is calculated to be 4.6.

Fig. 3(a) shows the reconstructed image using the proposed algorithm. The pixel size of the image is 256×256. It takes only 1 s to get the image. However,more than 2 min is needed to reconstruct the image processed by $BP$ algorithm. The horizontal lines are added to show the interfaces of layers. The targets are correctly located. In Fig. 3(a),the shadows at the interfaces of layered medium are introduced by the interaction between the objects and the layered medium. The image processed by F-K migration with RMS dielectric constant is shown in Fig. 3(b). As can be seen from Fig. 3(b),the targets are displaced away from the target locations. The image entropy is employed to quantitatively evaluate the focus quality of reconstructed images. The better focus corresponds to smaller entropy. The entropy is defined as^[17]:

$\xi = {{{{\left( {\sum {\sum {{{\left| {{x_{ij}}} \right|}^2}} } } \right)}^2}} \over {\sum {\sum {{{\left| {{x_{ij}}} \right|}^4}} } }}$

(25)

where ${\left| {{x_{ij}}} \right|}$ denotes the modulus at the grid cell $(i,j)$ of the reconstructed image. The values of image entropy of Fig. 3(a) and Fig. 3(b) are 880 and 9689,respectively. So Fig. 3(a) is more focused than Fig. 3(b). 4.2 Experimental results

The proposed imaging algorithm has been validated by experimentally for detecting an object embedded in a three-layered medium. The experimental setup consists of Agilent N5230A (10 MHz~20 GHz)VNA and AEL H-1498 doubleridge horn antenna mounted on a 1-D scanner. The scanner moves the antenna parallel to the layered medium with a standoff distance of 3 cm. A 30.48 cm aperture is synthesized with 25 measurements spaced 1.27 cm in $x$ direction. The aperture center is set to $x=0$,$z=0$,where $z=0$ corresponds to the measurement plane. Data is collected at 1001 points in the frequency range from 2 GHz to 12 GHz with a step of 10 MHz. The complex S-parameter S11 is used as the radar response for layered medium imaging. The scenario is shown in Fig. 4(a). The three layered medium is comprised by two wood blocks and concrete block. The thicknesses of wood block and concrete block are 1.85 cm and 4.46 cm,respectively. Two metallic disks with 1.3 cm radius are placed at the interface between the second layer and third layer,as shown in Fig. 4(b). The locations of the embedded object are (-2 cm,9.3 cm) and (3 cm,9.3 cm).

The dielectric constants of wood block and concrete block are estimated as 2.12 and 7.25 by time domain reflectometry^[17]. Fig. 5 shows the resulting image of embedded targets processed by 6 Journal of Radars Vol. x the proposed algorithm. The pixel size of the image is 256×256. The processing time is about 2.5 second. The indications of metallic disks in the image are focused and properly located. The interaction of the two objects leads to the ghost located at the center of the object locations. Faint artifacts at the left and right sides exist because the concrete block is not homogeneous.

5 Conclusion

In this paper,phase shift migration for imaging of layered medium is proposed. Based on the analysis of backscattered transfer function,phase shift migration for layered medium can be viewed as an approximation of inverse problem. The proposed algorithm for layered medium is efficient due to the following reasons. Firstly,the coherent summation in BP algorithm is computed efficiently with FFT. Secondly,the solving of nonlinear equation in order to find the refraction points in BP algorithm is avoided. Only the phase shift factor is updated along with depth. Thirdly,instead of the pixel-by-pixel reconstruction in BP algorithm,the proposed algorithm efficiently reconstructs the image with IFFT. Numerical and experimental results are presented to show the effectiveness of the proposed algorithm for real-time imaging of layered medium.

Acknowledgement The authors would like to thank ElectroScience Laboratory,the Ohio State University for providing the experimental data.