Using FDFD technique in two-dimensional TE analysis for modeling clutter in wall penetrating radar.
Previous work on EM field modeling for WPR commonly used time domain-based techniques to determine the true time-varying fields that are backscattered from the interaction of an excitation with the full 2D dielectric map of a scene. Several time domain approaches to resolving the inner layout of a building use ray tracing techniques with UWB radar imaging . In WPR, depending on the application (counter-terrorism, urban warfare, etc.), several aspects of the scene can be important. For typical WPR applications, a sensor transmits from a remote location into a dielectric barrier (clutter) of unknown internal configuration containing an unknown number of targets of interest (TOIs). For monostatic radar, the scattered field received at the sensor contains the direct scattered energy due to clutter plus the TOIs plus the interactions between them. Precise TOI localization requires isolating the scattered energy due to the TOIs alone.
The focus of this paper is around generating an accurate model of the clutter in the frequency domain. Due to the geometric complexity and inhomogeneity involved in a typical scene, precise model formulation containing minor detail of the internal structure becomes difficult. Instead, emphasis is placed on characteristics that contribute most to the scattered field. These include the overall dimensions and layout of the infrastructure, with focus on the orientation of the internal walls, doorway locations, and relatively large, static internal objects. TOIs typically have a significantly smaller radar cross-section (RCS) profile compared to the surrounding clutter and therefore contribute significantly less to the scattered field of the total scene. Human TOIs, which typically move throughout the scene, can have varying RCS as a function of time. This results in a reduced target-to-clutter ratio when integrating the full spatial scene over several points in time.
The performance of resolving small internal TOIs relies heavily on the accuracy of the clutter model. For typical construction, the clutter map will consist of several layers of walls. While partial clutter may be constructed from external observations, the internal configuration would remain hidden from any remote view point. Outside of the field scattered directly back from the outer edge of the external wall, closest to the sensor, the greatest contribution to the scattered field comes from the direct scattering and interaction between the internal and remaining external walls orientated perpendicular to the Poynting vector of the incident field. For a fixed sensor location, measuring the frequency response of a given clutter map, at several discrete frequencies, provides useful information regarding the internal dimensions and orientation of the dielectric scene. For a horizontally polarized electric field ([TE.sub.z] mode), the frequency can be tuned to look for objects of certain horizontal extent in typical construction, such as headers above doorways, while being desensitized to vertical features such as rebar and electrical wiring aligned with the studs.
2. Finite Difference
Because of the complicated geometries involved with a clutter model, analytic field solutions to Maxwell's equations are difficult to solve accurately. Finite difference computational electromagnetic techniques offer a means to accurately solve for the field solutions within these complicated environments.
2.1. FDTD. Time domain-based modeling techniques such as finite difference time domain (FDTD) have been shown to provide the necessary field analysis for determining the scattered field response, for a given excitation [2-4]. FDTD is fast and accurate but also can be quite computationally intensive for a large 2D space. While FDTD can provide a good time-dependent wave analysis through the walls of the structure (short pulse UWB scattering effects) 1], FDTD suffers from the time constant required to reach a steady-state response for a resonant structure. For 2D FDTD, good accuracy (accuracy [approximately equal to] [[DELTA].sup.2]) is maintained for a uniform spatial increment [DELTA] < [lambda]/10 in both dimensions. In order to maintain stability, the Courant condition also requires that the time increments be limited to [DELTA]t [less than or equal to] [DELTA]/([square root of 2]V), where V is the wave velocity (V = c/[square root of [epsilon].sub.r]. The dimensions involved in most WPR applications require a large computational space [I, J]. Along with this is that the geometries involved, with typical external and internal wall configurations, establish resonant structures. In the case of FDTD, a significant number of iterations may be required to settle transients associated with these resonances, and reach steady-state conditions, while maintaining stability. The CPU timing associated with 2D FDTD is proportional to 1/[[DELTA].sup.3].
2.2. FDFD. Finite difference frequency domain (FDFD) offers some advantages over FDTD, however is significantly more computationally intensive considering the dimensions involved with WPR. FDFD solves for the fields at all grid locations all at once by solving M = I * J simultaneous equations. The result is already the steady-state frequency response of the dielectric scene for a single frequency excitation. To maintain the equivalent accuracy of FDTD, the 2D uniform spatial increment is limited to [DELTA] < [lambda]/20. For a 1 GHz excitation in free space ([DELTA] < 15 mm) and a structure that extends 15 meters in both dimensions, it requires solving M > 1 * [10.sup.6] simultaneous equations. For 2D FDFD, the CPU timing is proportional to M/([[DELTA].sup.4] ln([DELTA])). Recognizing that the number of grid points that do not contain background dielectric is sparse within the full 2D space, the computational loading for FDFD is significantly reduced. Combining this with Moore's law, which states the number of transistors on integrated circuits doubles every two years, supports FDFD as a practical choice for real-world implementation. With the ability to assign each cell a unique permittivity and permeability value, FDFD is convenient for handling media with nonuniform, frequency dependent properties, which require special methods in the time domain .
3.1. Dielectric Space. For this paper, a dielectric scene, consisting of five small (~5 cm diameter) metallic cylinders and a human elliptic sphere (TOIs) within a rectangular structure, is constructed over a subdivided free space background (Figure 1). The ultimate goal of resolving the location of the TOIs requires accurately removing the effects of the clutter, which is the focus of this paper. In this case, all structure-related dielectric (walls, headers, rebar, etc.) is classified as clutter. The external wall is modeled using 0.5 meter thick dielectric with an effective dielectric constant ([[epsilon].sub.r_ext]) of 4. To accurately represent typical concrete or cinder block wall construction, pairs of 2.5 cm thick vertical rebar are placed within the wall at an interval of 1 meter around the perimeter with 4 larger 10 cm beams at each corner of the external structure. The internal walls are 14 cm thick dielectric ([e.sub.r_inner] = 4), containing metallic door headers which represent a transition point between rooms. An actual scene may or may not have a closed door below each header at these locations; however each transition point will most likely contain a header, providing a consistent feature in determining the internal layout of the structure and making it worthy of focus in this analysis. The physical extent of headers for standard doorways is predictable and provides a significant scattered field to a transverse electric polarized incident field. In this case, the headers were set to 1 meter in length.
3.2. FDFD Analysis. This paper proposes using 2-dimensional finite difference frequency domain transverse electric mode (2D-FDFD [TE.sub.z]) analyses  to solve for the total scattered fields from the full scene, using an 8-layer perfectly matched layer (PML) along the edge of the computational 2D space [7, 8]. If 2D space is represented as an I x J grid, FDFD determines the steady-state field response at all grid locations by simultaneous solving of a set of M = I * J equations. With FDFD, a tradeoff between computational accuracy and time is used to determine the grid size for a given excitation frequency [f.sub.s], with good accuracy ([approximately equal to] [[DELTA].sup.2]) maintained for grid spacing [DELTA] < [lambda]/20. A [TE.sub.Z] ([E.sub.x], [E.sub.y], [H.sub.z]) mode analysis was chosen to mitigate the dispersive effects through the external walls containing vertically aligned rebar. For [TE.sub.z], the strongest scattered fields will originate at locations that have a linear dielectric boundary with significant extent in the (x, y) plane. This will be most sensitive to wall-air boundaries especially where the dielectric has the greatest value within the wall, such as headers above door openings or other horizontal supports within the internal structure.
With the dielectric mapping of the full scene given by k(x,y), the background wave number [k.sub.b](x,y) consists of free space where the regions containing nonfree space are expressed as [k.sub.b](x,y) + [DELTA]k(x,y). Assuming a source free space, the Helmholtz equation states the following for the incident ([H.sup.i.sub.z]) and total ([H.sup.i.sub.z] + [H.sup.s.sub.z]) fields:
[nabla] x (1/[k.sup.2.sub.b](x,y))[nabla] x [H.sup.i.sub.z] - [H.sup.i.sub.z] = 0 (1)
[nabla] x (1/[k.sup.2](x,y)) [nabla] x ([H.sup.i.sub.z] + [H.sup.s.sub.z]) - ([H.sup.i.sub.z] + [H.sup.s.sub.z]) = 0 (2)
Using (1), we can rewrite (2), solving for the scattered fields,
[nabla] (1/[k.sup.2](x,y))[nabla] x [H.sup.s.sub.z] - [H.sup.s.sub.z] = - [nabla] x D(x,y) [nabla] x [H.supi.sub.z] (3)
with D(x,y) = (1/[k.sup.2](x,y) - 1/[k.sup.2.sub.b](x,y)).
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)
If we assume uniform grid spacing in both [increment of x] = [increment of y] = [DELTA], the wave number for each grid location is given as [k.sup.2](i, j) = [[omega].sup.2][[mu].sub.o][[epsilon].sub.o[[epsilon].sub.r](i[DELTA], j[DELTA]). For 2nd order derivation, the subgrid values for the wave number are needed and expressed as [k.sup.2](i+ 1/2, j) = ([k.sup.2](i+ 1, j) + [k.sup.2](i, j))/2. The finite difference equation of (4) represented in grid space is given by
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)
3.2.1. Incident Field Generation. The incident field [H.sup.i.sub.z] is found at all grid locations for a single frequency (CW at [f.sub.s]), using free space for the background dielectric ([k.sub.b] = [k.sub.0]) at all grid locations. The finite difference equation of (1), represented in grid space, is then given by
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)
For an I x J computational space with M = I * J total grid locations, FDFD determines [H.sup.i.sub.z] at each location by solving M simultaneous equations through the following relationship:
[A] [[H.sup.i.sub.z]] = [S], (7)
where [A] is an M x M sparse operator matrix with each row containing (4 - [k.sup.2.sub.o][[DELTA].sup.2]) along the diagonal and 1 at the corresponding (i + 1, j), (i - 1, j), (i, j + 1), and (i, j - 1) locations. [[H.sup.i.sub.z]] is an Mx1 matrix containing the field solution at each grid location and [S] is an Mx1 source matrix of zeros except for at the grid locations of the excitations. The row of A corresponding to an excitation location replaces the ones with zeros as well as the (4-[k.sup.2.sub.o][[DELTA].sup.2]) with a value of 1. As mentioned earlier, this space is surrounded by an 8-layer PML to handle the wave conditions along the boundary.
A linear array of equally spaced current source excitations, orientated parallel to the external wall interface, is used to excite free space. For a single [??] directed current source excitation located at the half interval [[i.sub.s], [j.sub.s] + 1/2], the energy of the adjacent grid X indices is set to nonzero values with opposing polarity. Consider
S [I x ([J.sub.s] - 1) + [i.sub.s] = [H.sub.js] [??],
S[I x [j.sub.s] + [i.sub.s] = -[H.sub.js][??] (8)
Likewise, the [??] directed current source at the half interval [[i.sub.s] + 1/2, [j.sub.s]] is generated by setting the energy of the adjacent grid Y indices equivalently
S [I x ([J.sub.s] - 1) + [i.sub.s] + 1] = [H.sub.js] [??],
S [I x ([J.sub.s] - 1) + [i.sub.s]] = [H.sub.js] [??](9)
Figure 2 shows the incident field generated from a left side vertically aligned linear array using 23 uniformly spaced [??] directed current source excitations. Using the established dielectric scene of Figure 1 and the numerically generated incident field of Figure 2, the values [H.sup.i.sub.z] are used in (5) to determine the scattered field at each grid location, as shown in Figure 3. The total field at each grid location is then found from adding the scattered field to the incident field, [H.sup.T.sub.z] = [H.sup.s.sub.z] + [H.sup.i.sub.z], as shown in Figure 4.
4. Numerical Results
4.1. Clutter Removal: Direct Subtraction. The dominant sources of the scattered field are primarily the external and internal wall dielectrics and their interaction with one another, as seen in Figure 3. TOI localization through imaging relies on the ability to isolate the scattered field due to the TOIs alone. One approach to removing the effects of the clutter elements is to directly subtract the scattered fields of the known clutter, determined through FDFD modeling, from the frequency response of the full scene. There may be several layers to the clutter where each layer is removed at a different stage of analysis. For example, an initial clutter model could be the external wall closest to the sensor, with the internal wall and hidden external wall configurations unknown. Several observed factors maybe used in constructing the initial model. For instance, the visible material, outer dimensions, and some architectural rules of thumb, such as the use of rebar and other structural reinforcing material in buildings of certain height and extent, can all be used as factors for the initial model. In this case, there are multiple layers of clutter. Therefore several stages of clutter removal are required, where each stage provides additional maturity to the overall model.
Using the external walls as the initial clutter model (Figure 5), the remaining scene of Figure 1 consists of the internal walls and TOIs and is shown in Figure 6. Figure 7 shows the scattered fields ([H.sup.s.sub.z-Ext-Walls]) due to the dielectric scene of Figure 5 (external walls alone) and Figure 8 shows the scattered fields ([H.sup.s.sub.z-CS2+Tois]) due to the dielectric scene of Figure 6 (clutter stage2 + TOIs). Subtracting the scattered fields due to the external walls alone (Figure 7) from the scattered fields due to the full scene (Figure 3) provides the difference field of Figure 9. Consider
[H.sup.s.sub.z_Difference] = [H.sup.s.sub.z-Full] - [H.sup.s.sub.z-Ext-Walls] (10)
Comparing Figures 8 and 9, the fields within the externals walls are consistent for the most part, with the exception of the various modes generated by [H.sup.s.sub.z-CS2+TOIs] interacting with the external walls. The fields that would be directly measured by a remote sensor, outside of the external walls, vary significantly between the two figures. This is due primarily to the transfer function through the external walls as well as other second order clutter effects (external to external, internal to internal, and external to internal clutter interactions). This provided motivation to be able to fully characterize both stages of clutter in order to be able to provide an accurate model of the full clutter scene, due to the strong clutter-to-clutter interactions of Figure 9. For image reconstruction techniques that use subspace analysis, such as MUSIC and SVD, the ability to isolate the singular values that belong to the targets of interest while minimizing the effects of clutter relies heavily on how accurately the clutter can be identified and therefore removed from the total scattered field, especially in regions of strong clutter interaction. With the full clutter map of Figure 10 identified, the scattered field due to this model ([H.sup.s.sub.z-Clutter]) is removed from [H.sup.s.sub.z-full] to produce the difference field of Figure 12. Comparing this difference field with the modeled scattered field due to the TOIs alone (Figure 11) shows good comparison of fields in the outer region. Removing the full clutter map would therefore provide a remote sensor with the best opportunity to locate the TOIs, in the presence of strong clutter, through imaging.
4.2. Construction of the Full Clutter Map. As mentioned, the dominant scattering in [H.sup.s.sub.z-full] is due to the external and internal wall dielectric configuration. Figure 3 shows that the dominant scattering occurs at dielectric boundaries running perpendicular to the poynting vector of the incident field. In more practical cases, observation of the external wall closest to the sensor will be the easiest, especially at a remote location. Therefore the section of the external structure closest to receiver is used as the initial clutter model (Figure 13) to generate the difference fields ([H.sup.s.sub.z_Difference]) through direct subtraction as described in the previous section. Using FMCW, the received data is processed in the frequency domain. The modeled frequency response of the clutter can be subtracted from the received frequency response of the scene to precondition the data before applying a fast fourier transform (FFT) process for determining the range profile. For the experiment in this paper, we measure the received energy across a 23-element array at each discrete frequency, across a stepped frequency ramp of bandwidth BW, at a frequency step interval of [DELTA]f. A range profile is determined through a 1-D FFT across the stepped frequency range with a maximum unambiguous range of [R.sub.unamb] and range resolution per FFT bin of [R.sub.res]. Consider
[R.sub.unamb] = c/2[DELTA]f (11)
[R.sub.res] = c/2BW (12)
The energy of each range bin is assigned a perpendicular distance from the receiver according to (12), for each array element. Determining the difference field across the array elements for a given receiver location is the first step in imaging the clutter. [H.sup.s.sub.z_Difference] is the scattered field difference that results from subtracting the scattered fields due to the dielectric maps of Figures 1 and 13. Setting the frequency range to 600-1800 MHz ([R.sub.res] [approximately equal to] 0.125 m) and a frequency step ([DELTA]f) to 2 MHz ([R.sub.unamb] [approximately equal to] 75 m), a 2D image is updated using the logarithmic magnitude of the 1-D FFT bin array for each receive array element at a perpendicular range of [R.sub.bin] = bin * [R.sub.res].
Resolving the different internal dielectric interface configurations requires excitations from all directions along the outer perimeter of the structure. For example, vertically aligned walls in Figure 1 are dominate scatterers for left and right side transmit-receive (TR) locations, as well as horizontally aligned walls being dominate for top and bottom TR locations. The frequency response of the full scene (Figure 1) is recorded at the receive elements by sweeping the TR locations along the outer edge of the left external wall across the full 1200 MHz bandwidth. The left side sweep is repeated, this time capturing the frequency response of Figure 13(a) as the model of the known clutter. A difference field is generated as a result of subtracting the frequency response from the two dielectric scenes at the receive elements for each receiver location. An FFT of the result provides down range information from the receiver which is mapped onto a 2D image with 7 mm pixel resolution. Since the FFT range bins have 125 mm resolution, the image is greatly oversampled. For each TR location, a 2D blurring function (spatial low pass filter) is applied to the FFT data to spread the energy over several adjacent pixels. The results are integrated across the left edge sweep and shown in the image of Figure 14. With the effects of the left external wall (clutter map) subtracted, the FFT energy peaks at the locations of the vertically aligned inner wall and right side external wall, now the dominate scatterers for the left side TR locations. There is enough information from the left edge sweep analysis to locate the vertical inner wall as well as the opposing external wall; however horizontal features are difficult to observe when probing the scene from this aspect. The results hold for the right edge scan shown in Figure 17. Horizontally aligned clutter features are more obvious from the top and bottom edge scans shown in Figures 15 and 16, respectively.
Careful inspection of Figure 9 shows significant scattering from the lower left header guided back to the receiver along the inner wall dielectric. Figure 14 shows a peak in FFT energy where the internal wall intersects the left external wall. This represents the focused energy from the lower left header. Figures 15 and 16 show significant FFT energy at a linear extent from that intersection, caused by the interaction of the horizontally aligned header with the incident fields from top and bottom TR sweep locations, respectively. These features provide a high confidence factor that an internal wall resides at this location. This can be observed for the other 3 equivalent intersections points around the structure, combining the results of Figures 14-17.
The integrated result of sweeping the TR location along all four sides of the building is shown in Figure 18. Using these images, a simple algorithm could be implemented that looks for peak energy above a threshold, with linear features matching the characteristics of the inner walls (matched filter).
Furthermore, within those regions, an additional threshold is established to determine the location of the headers based on typical header sizes and locations. Implementing such an algorithm on Figures 14-17 produced the net result of Figure 19, which averages the results of applying the header matched filter (HMF) algorithm to each image. The result of Figure 19 correlates well with the actual header locations of Figure 10. This simple algorithm is adequate for identifying the domain scattering features of the internal clutter. A more advanced algorithm could provide a more detailed clutter map for higher order classification. TOI localization will ultimately rely on the accuracy of the clutter model.
In this paper, FDFD was used to formulate the scattered fields from a 2D dielectric scene of a typical WPR application. The goal of characterizing dominate scattering features for an arbitrarily configured clutter map started with an initial model of the external wall, closest to the sensor, from visual observations. Based on typical construction, vertical rebar was included in the external wall model to determine its effect on TEZ mode propagation. It was shown that FDFD provides enough accuracy to model and ultimately remove the effects of clutter, through direct subtraction, at the cost of being computationally expensive for the large structures typical in WPR application. FDFD is geometrically versatile and has its advantages over FDTD for complex shaped and inhomogeneous structures. FDFD is a convenient technique, and with the advances in available computer resources, it becomes a practical and robust means for solving the 2D field solutions for a variety of WPR environments.
This paper focuses on the construction of the clutter map in 2D space. In practical application, the clutter itself represents the dominate scattering with internal targets of interest, contributing a relatively small amount to the scattered field of the full dielectric scene. This allows the dominate features of the clutter to be characterized, even in the presence of TOIs. Starting with an initial clutter model, the frequency response, as observed from the sensor, is generated and stored for a range of discrete frequencies. The data can be accessed and directly subtracted from the frequency response of the received signal in FMCW radar to remove the effects of the known clutter. An FFT of the result provides range data which, when combined with location of the sensor, provides a 2D mapping of the scattered scene. Repeating this process for a number of receiver locations along the outer edges of the external structure reveals dominate scattering features of the internal structure. A specific interest for this paper was to identify doorway headers of known extent. This provides useful information concerning the internal configuration of the structure from a remote location. Headers typically have fixed horizontal dimensions which can be sensitive to certain frequencies for TEZ mode wave propagation. A simple threshold-based algorithm could be implemented to resolve header locations and help to map the internal wall configuration throughout the structure. It was shown that the integrated results improve significantly as the scene is observed from several locations around the exterior of the structure.
Future work in this area will include taking the next steps to ultimately resolve the locations of the TOIs within a known clutter environment. This paper acts to establish a process for generating the clutter model which will lead to imaging of the 2D scene and include the effects of the TOIs alone. This analysis will include the impact that an incorrect clutter model (dimensions and dielectric properties) will have on the accuracy of TOI localization.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
 C. Le, T. Dogaru, L. Nguyen, and M. A. Ressler, "Ultrawideband (UWB) radar imaging of building interior: measurements and predictions," IEEE Transactions on Geoscience and Remote Sensing, vol. 47, no. 5, pp. 1409-1420, 2009.
 X. Gu and Y. Zhang, "Autofocus imaging simulation for through-wall radar by using FDTD with unknown wall characteristics," in Proceedings of the Asia-Pacific Microwave Conference (APMC '10), pp. 1657-1660, December 2010.
 B. Boudamouz, P. Millot, and C. Pichot, "Through the wall radar imaging with MIMO beamforming processing," in Proceedings of the 3rd Microwaves, Radar and Remote Sensing Symposium (MRRS '11), pp. 251-254, August 2011.
 T Dogaru and C. Le, Through-the-Wall Radar Simulations for Complex Room Imaging, Army Research Laboratory, 2010.
 L. Kuzu, V. Demir, A. Z. Elsherbeni, and E. Arvas, "Electromagnetic scattering from arbitrarily shaped chiral objects using the finite difference frequency domain method," Progress in Electromagnetics Research, vol. 67, pp. 1-24, 2007.
 C. Rappaport, A. Morgenthaler, E. Bishop, Q. Dong, and M. Kilmer, "Finite difference frequency domain (FDFD) modeling of two dimensional TE wave propagation and scattering," in URSI Symposium Digest, pp. 1134-1136, Pisa, Italy, May 2004.
 J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of Computational Physics, vol. 114, no. 2, pp. 185-200, 1994.
 E. A. Marengo, C. M. Rappaport, and E. L. Miller, "Optimum PML ABC conductivity profile in FDFD," IEEE Transactions on Magnetics, vol. 35, no. 3, pp. 1506-1509, 1999.
David Insana (1) and Carey M. Rappaport (2)
(1) Northeastern University, Boston, MA 02115, USA
(2) CenSSIS, Northeastern University, Boston, MA 02115, USA
Correspondence should be addressed to David Insana; email@example.com
Received 19 July 2013; Revised 28 November 2013; Accepted 12 December 2013;
Published 19 January 2014
Academic Editor: Danilo Erricolo
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Author:||Insana, David; Rappaport, Carey M.|
|Publication:||International Journal of Antennas and Propagation|
|Date:||Jan 1, 2014|
|Previous Article:||A 2.4 GHz cross rhombic antenna for a cube satellite application.|
|Next Article:||Indoor massive MIMO channel modelling using ray-launching simulation.|