Ultra wideband ground penetrating radar imaging of heterogeneous solids

• US 5,835,054 A
• Filed: 05/09/1997
• Issued: 11/10/1998
• Est. Priority Date: 03/01/1996
• Status: Expired due to Fees

• Current Assignee: Lawrence Livermore National Security LLC (Government of the United States of America)
• Original Assignee: Regents of the University of California (University of California)
• Inventors: Warhus, John P., Mast, Jeffrey E.
• Primary Examiner(s): Sotomayor, John B.
• Application Number:
  US08/853,548
• Time in Patent Office: 550 Days
• Field of Search: 342/22, 342/25, 342/27, 342/179, 342/194, 342/195, 342/197, 342/196, 342/180
• US Class Current: 342/22
• CPC Class Codes: G01S 13/0209 Systems with very large relative bandwidth, i.e. larger than 10 %, e.g. baseband, pulse, carrier-free, ultrawideband
• CPC Class Codes: G01S 13/867 Combination of radar systems with cameras
• CPC Class Codes: G01S 13/885 for ground probing prospecting or detecting using electromagnetic waves G01V3/12
• CPC Class Codes: G01S 13/89 for mapping or imaging
• CPC Class Codes: G01S 13/90 using synthetic aperture techniques , e.g. synthetic aperture radar [SAR] techniques
• CPC Class Codes: G01S 7/292 Extracting wanted echo-signals
• CPC Class Codes: G01V 3/12 operating with electromagnetic waves operating with millimetre waves G01V8/005

Abstract

A non-invasive imaging system for analyzing engineered structures comprises pairs of ultra wideband radar transmitters and receivers in a linear array that are connected to a timing mechanism that allows a radar echo sample to be taken at a variety of delay times for each radar pulse transmission. The radar transmitters and receivers are coupled to a position determining system that provides the x,y position on a surface for each group of samples measured for a volume from the surface. The radar transmitter and receivers are moved about the surface, e.g., attached to the bumper of a truck, to collect such groups of measurements from a variety of x,y positions. Return signal amplitudes represent the relative reflectivity of objects within the volume and the delay in receiving each signal echo represents the depth at which the object lays in the volume and the propagation speeds of the intervening material layers. Successively deeper z-planes are backward propagated from one layer to the next with an adjustment for variations in the expected propagation velocities of the material layers that lie between adjacent z-planes.

Claims

1. A linear array imaging radar system, comprising:

a linear array of radar transceivers connected to provide a plurality of radar echo samples taken at a variety of delay times, and thus depths into an imaged solid, for each of a series of radar pulse transmission from each radar transceiver;

a position determining system connected to the linear array of radar transceivers and a multiplexer for them that provides an x,y position over a surface of said imaged solid for each group of samples measured for a volume from said surface, wherein the radar transmitter and receiver are moved about said surface to collect such groups of measurements from a variety of x,y positions, and wherein a plurality of return signal amplitudes represent the relative reflectivity of objects within said volume and the delay in receiving each signal echo represents the depth at which an object lays in said volume and the propagation speeds of the intervening material layers; and

means for successively backward propagating deeper z-planes from one layer to the next with an adjustment for variations in the expected propagation velocities of a variety of material layers that lie between adjacent said z-planes, wherein each said z-plane contributes to a tomographic series of slices through said volume that contribute to the visualization of said objects in said volume in three dimensions.

2. The system of claim 1, wherein the means for successively backward propagating comprises a computer-implemented process for:

placing said data in the frequency domain with a fast Fourier transform (FFT) with respect to time, the result being mathematically expressible as r(x,y,f);

scaling the frequency with media propagation velocity to operate with lambda, wherein the result may be expressed conveniently as r(x,y,λ

);

fast Fourier transforming with respect to (x,y) to yield a result expressible as R(fx,fy,λ

), where "fx" and "fy" denote spatial frequencies in the x-axis and y-axis;

backward-propagating (110), if this is the first time through, to a next deeper layer by multiplying R(fx,fy,λ

) with a spatial filter H(fx,fy,λ

,Δ

z) to yield image layer I(fx,fy,λ

,z_i+1), otherwise multiplying said spatial filter with a previous result, I(fx,fy,λ

,z_i), to backward propagate from one deep layer to the next deeper layer, said spatial filter being represented by, ##EQU3## where the factor of two inside the square root accounts for the two-way travel path of the radar pulse, the spatial spectrum of the equivalent source is,

space="preserve" listing-type="equation">S(f.sub.x,f.sub.y,z)=∫

R(f.sub.x,f.sub.y,λ

)H.sup.-1 (f.sub.x,f.sub.y, λ

)dλ

,

and where R is the received wavefield and H^-1 is the pseudo inverse, the estimate of the source distribution is given by the inverse two-dimensional Fourier transform S, wherein said process of backward propagation begins with a starting z-plane that can be equated to a hologram, and both time and phase information are recorded in such starting z-plane that represent what is illuminating it in radio frequency from deeper z-planes and projecting only to the next deeper z-plane that is then used as a basis for projecting to a still deeper z-plane, and wherein all said projected z-planes are an equal distance Δ

z apart;

setting, when still deeper z-planes are to be resolved, a next z-level and adjusting for any expected change in propagation velocity for a new media;

integrating or summing said images I(fx,fy,λ

,z) to form an output (Ifx,fy,t) and super-positioning all said information collected at various radar frequencies for improving estimates of positions of objects; and

inverse Fast Fourier transforming with respect to (fx,fy) to yield a set of final images for said z-planes, expressed as I(x,y,z).

3. The system of claim 1, wherein the positioning determining system comprises:

a position transmitter that periodically transmits a radio frequency pulse that is received and measured by at least three position receivers (16-18), wherein the free-space time-of-flight to each receivers (16-18) is compared and the position of an acquisition unit (12) is solved for by triangulation and said positions determined contribute to the solution of the x,y position of said data acquisition unit (12).

4. The system of claim 1, wherein the positioning determining system comprises:

a navigation receiver that determines its position from a constellation of orbiting satellites and said position is connected to an image solution unit with each corresponding collection of radar imaging data.

5. The system of claim 4, wherein the positioning determining system further comprises:

a real-time kinematic (RTK) on-the-fly (OTF) position and navigation method for sub-centimeter accurate position solutions and a stationary reference station with a known position obtained by an independent high-accuracy survey, wherein signals are received from said constellation and RTK solutions are initialized by resolving the carrier integer ambiguities to each satellite in the constellation in said navigation receiver by using phase information from each carrier and timing and code ranging information to compute the apparent sub-centimeter level position of the reference station, and the known position is subtracted from the apparent position solution, and the correction and RTK-initialization data are then radio-linked constantly to said navigation receiver, and wherein real-time sub-centimeter position solutions are provided for said data acquisition unit by said satellite navigation receiver to an image reconstruction unit.

6. The system of claim 1, wherein the radar transmitter and receiver connected to a timing mechanism include timing pulses that are transmitted by a position transmitter or said radar transmitter and are time multiplexed.

7. A computer process for extracting tomographic images of a solid from a synthetic aperture radar swept through various points in a plane of observation, the method comprising the steps of:

collecting data from time delay information (t) provided by a radar signal penetration of a solid from a plane of observation and matching said data with the corresponding positions (x,y) of each observation, wherein said data set may be referred to in the convenient mathematical form as r(x,y,t);

placing said data in the frequency domain with a fast Fourier transform (FFT) with respect to time, the result being mathematically expressible as r(x,y,f);

scaling the frequency with media propagation velocity to operate with lambda, wherein the result may be expressed conveniently as r(x,y,λ

);

fast Fourier transforming with respect to (x,y) to yield a result expressible as R(fx,fy,λ

), where "fx" and "fy" denote spatial frequencies in the x-axis and y-axis;

backward-propagating, if this is the first time through, to a next deeper layer by multiplying R(fx,fy,λ

) with a spatial filter H(fx,fy,λ

,Δ

z) to yield image layer I(fx,fy,λ

,z_i+1), otherwise multiplying said spatial filter with a previous result, I(fx,fy,λ

,z_i), to backward propagate from one deep layer to the next deeper layer, said spatial filter being represented by, ##EQU4## where the factor of two inside the square root accounts for the two-way travel path of the radar pulse, the spatial spectrum of the equivalent source is,

space="preserve" listing-type="equation">S(f.sub.x,f.sub.y,z)=∫

R(f.sub.x,f.sub.y,λ

)H.sup.-1 (f.sub.x,f.sub.y, λ

)dλ

,

and where R is the received wavefield and H^-1 is the pseudo inverse, the estimate of the source distribution is given by the inverse two-dimensional Fourier transform S, wherein said process of backward propagation begins with a starting z-plane that can be equated to a hologram, and both time and phase information are recorded in such starting z-plane that represent what is illuminating it in radio frequency from deeper z-planes and projecting only to the next deeper z-plane that is then used as a basis for projecting to a still deeper z-plane, and wherein all said projected z-planes are an equal distance Δ

z apart;

setting, when still deeper z-planes are to be resolved, a next z-level and adjusting for any expected change in propagation velocity for a new media;

integrating or summing said images I(fx,fy,λ

,z) to form an output I(fx,fy,t) and super-positioning all said information collected at various radar frequencies for improving estimates of positions of objects; and

inverse Fast Fourier transforming with respect to (fx,fy) to yield a set of final images for said z-planes, expressed as I(x,y,z).

8. The method of claim 7, wherein:

the step of collecting data from time delay information (t) provided by said radar signal penetration the scatter from buried objects is assumed to follow Born'"'"'s approximation for weak scattering and the material of any intervening layers is assumed to be homogeneous; and

wherein successively deeper layers are assumed also to be homogeneous with any difference between layers being accounted for by adjusting the propagation velocity "v" in the step of backward-propagating (110), such that the approximation linearizes the scattering problem and relates the transmitted wavefield "u", at a given wavelength "λ

", the object distribution "o(x,y,z)" and the received wavefield "r" with,

space="preserve" listing-type="equation">r(x,y;

λ

)=∫

o(x,y,z'"'"';

λ

)u(x,y,z'"'"';

λ

)!**g(x,y,z'"'"';

.lambda.)dz'"'"',

where "g" is the Green'"'"'s function and "**" indicates convolution, and the Green'"'"'s function for homogeneous background material is, ##EQU5## where r =√

x² +y² +z² , and k =2π

/λ

, and with monostatic illumination with an isotropic radiator, the received wavefield becomes,

space="preserve" listing-type="equation">r(x,y;

λ

)=∫

o(x,y,z'"'"')**g.sup.2 (x,y,z'"'"';

λ

)dz'"'"',

such that the dependency on g² is removed by taking the derivative with respect to k,

space="preserve" listing-type="equation">r(x,y,λ

)=∂

/∂

kr(x,y,λ

)=C∫

o(x,y,z'"'"')**g (x,y,z'"'"';

λ

/2)dz'"'"',

and where C is a constant.