UNDERICE, AND INSIDE DISTRIBUTED TARGETS

TOMOGRAPHIC DOPPLER IMAGING

REMOTE SENSING, EARTH OBSERVATION, RADAR

[0001] Abstract: The present invention relates to a method for sensing complex tomographic images of the internal matter of any object and below the Earth/sea/ice, for kilometres, from space satellite and airborne platforms. Processing the coherent vibrational Doppler information in the single Synthetic Aperture Radar (SAR) im- age, in the single-look-complex (SLC) and/or raw configuration, the sound/phonon information is exploited, and penetrating tomographic spatial high-resolution imag- ing over a depth of several kilometers from the Earth’s surface is allowed. Dense tomographic measurements of underground, undersea, under-ice and inside extended targets are extrapolated and exploited from any SLC and/or raw SAR images of any satellite and airborne with any geometry, polarization, Spotlight and/or Stripmap sensor modes, any transmission frequency, and any chirp-Doppler bandwidths. This patent allows a new observation domain of SAR sensors: Underground, undersea, under-ice and inside solid bodies.

BACKGROUND OF THE INVENTION

[0002] The invention is a processing method embedded into a software that pro- duces tomographic imaging below/inside the Earth, below/inside the ice, below /inside the sea, and inside any body made of matter which is located on the Earth/sea/ice surface. The results are complex tomographic images of the internal matter, focused in high-resolution. Data are released at any information system (GIS) standard.

The invention processes any Synthetic Aperture Radar (SAR) data at any process- ing level, and accepts any SAR satellite/air borne complex data, observed operating at any frequency, at any chirp and Doppler bandwidths, at any acquisition mode, like Scansar, Stripmap, and Spotlight, at any polarization mode. The minimum is one SAR image and also multiple SAR images at any observation geometry like incidence angle and squint orientation.

The processing system/software allows the possibility to add a new observation domain to any SAR satellite/airborne systems. This domain is anything observable by acoustic- waves below/inside the Earth, below/inside the ice, below/inside the sea and inside any subsurface and surface man-made object like buildings, bridges, dams, and any other object made of matter (natural or man-made).

The processing method is traduced into a system/software, implemented using any programming language, and compiler. It accepts as input a single or multiple SAR images in any configuration (processing level), especially single- look-complex generated at any precision in terms of radiometric resolution.

The processing system/software allows the synthesis (we mean the extrapolation) of vibrational/phononic/sounding information present on all pixels of the satellite SAR data, using the Doppler centroids displacement/shifts and abnormalities of targets due to micro-metric and any other nature of motion, and higher order motions. The tomographic analysis is performed by focusing the complex sound information, using any known technique based on fast Fourier transform, and pulse compression. The processing system /software extrapolates any nature of mechanical information that is always embedded in the Earth surface matter, due to any kind of mechanical perturbation like seismic perpetual ripple and man-made vibrations due to human activity, also if expressed at very low magnitude. The sounding data extrapolation is performed by measuring sub-pixel spatial displacement due to electromagnetic waves and matter Doppler interaction of multiple Doppler sub-apertures SAR data (refocused at lower azimuth spatial resolution).

This invention with respect to any other sounding in-situ system imaging that has been already patented, has the crucial difference that is made from space/ airborne and using any satellites for SAR Earth observation data. This invention recasts the innovative solution of airborne and Space-Sonar, that extrapolates the vibration information from electromagnetic waves-matter interaction (photon versus phonon interaction).

DESCRIPTION OF THE INVENTION

[0003] One of the main issues with SAR is that due to the poor penetrating action of electromagnetic waves within solid bodies, the ability to observe through distributed targets from space is today precluded. Indeed, imaging is only possible on targets distributed on the scene surface. This invention describes an imaging method based on the analysis of micro-motions of seismic waves and any other mechanical wave present in the underground of Earth’s surface.

The only way to penetrate matter for kilometers is to use vibrations: we observed this communication channel to be present within all SAR images. So as we discovered this embedded information entropy, we are patenting the extrapolation method, using the Doppler sub-aperture strategy and measuring the displacement of any pixels through any multiple Doppler sub-aperture SAR image.

The SAR synthesizes the electromagnetic image through a "side looking” acquisi- tion, according to the observation geometry shown in Figure 0.1, where:

• r is the zero-Doppler distance (constant);

• R is the slant-range;

• R ₀ is the reference range at t = 0;

• d _a is the physical antenna aperture length;

• V is the platform velocity;

• d is the distance between two range acquisitions;

• G _sa is the total synthetic aperture length;

• t is the acquisition time variable;

• T is the observation duration;

• t = 0 and t = T are the start and stop time acquisition respectively;

• L = is the azimuth electromagnetic footprint width;

• θ is the incidence angle of the electromagnetic radiation pattern.

The invention estimates the micro-motion present on any target located on the Earth surface processing a single or multiple SAR images in the single-look-complex or any other configuration, by any specific method of pixel-tracking. The seismic ripple or any other kind of micro-motion is estimated through the Doppler sub-apertures anal- ysis of SAR data. Multiple Doppler sub-apertures, SAR images with lower azimuth resolution, are generated to estimate the vibrational trend of all pixels constituting SAR data. The infra-chromatic displacement is calculated through any technique that measures sub-pixel complex displacement. The displacement is generated by the micro-motion and any higher-order variation of micromotion. The temporal trend of this micro-motion is calculated by evaluating the orbital linear evolution of the Doppler spectrum present in the single image, according to the strategy depicted in Figure 1.

[0004] Here we explain how mechanical vibrations extrapolation is made. The total Doppler band of the SAR data are , while the total chirp band is . In the image, two squares are observed, one black and one blue. The black square has a base equal to where the parameter represents the Doppler band that is accepted to be lost. We found the best value of This means that the master image is focused using the entire chirp band with amplitude and half of the Doppler band The slave image, again estimated from the spectrum of the same SAR image made by acquiring the entire Doppler band, slightly shifted along the azimuth frequency by which represents the precise vibrational frequency one wishes to observe. The higher this parameter, the lower the mechanical frequency observed. Once the master and slave images have been refocused in accordance with the diagram depicted in Figure 1, these bands are shifted continuously until they reach position (2) and finally position (3). The orbital shift determines the temporal trend of vibrations. Each small displacement of the two bands (master and slave), rigidly held at a distance B _shift will have an expression in the displacement calculation in the image domain, where a continuous variation of the position of the green square (slave) of Fig. 0.3 will be observed. This estimated temporal trend of spatial shifts is just the estimated micro-vibration on the pixel considered. The above-mentioned procedure is repeated on all pixels constituting the SAR data. Fig. 0.4 (a), (b), (c), and (d) represent the vibration acquisition scheme in more detail. In Fig. 0.4 (a), surface vibrations are modeled by means of a system made by several springs held by masses. These facilities oscillate perpetually, providing their own oscillatory contribution to each pixel in the image

(in the present case we consider a strip of adjacent pixels representing a tomographic line of investigation). These springs generate a vibration that is visible in Fig. 0.4 (b). Fig. 0.4 (c) depicts the spring in detail that can be extended dynamically from L towards L ₀ . The last sub-picture represents the satellite acquisition scheme along the orbit and we can see how the strategy of partial focuses employing lower Doppler bands allows the system to observe mechanical vibrations. The parameter represents the spatial orthogonal baseline that is used to focus the tomograms. We suppose now the spring being perturbed by an impulse force.

According to this perturbation the rope begins to vibrate describing an harmonic motion (in this context we are not considering any form of friction). The resulting perturbation moves the rope through the space-time in the form of a sinusoidal func- tion. The seismic wave will then reach a constraint end that will cause it to reflect in the opposite direction. The reflected wave will then reach the opposite constraint that will make it reflect back in the original direction and return to the initial loca- tion, maintaining the same frequency and amplitude. According to Classical Physics principles, the rebounding wave is superimposed on the arriving wave, and the in- terference of two sine waves with the same amplitude and frequency propagating in opposite directions leads to the generation of an ideal and perpetual standing wave on the spring. Each vibrational channel is now considered when the spring is able to oscillate into the three-dimensional space, according to specific perturbation nature. When the Earth vibrates, it happens that the length of the spring must also fluctuate.

This phenomenon causes oscillations in the tension domain of the spring. It is clear that these oscillations (i.e. the longitudinal ones) propagate through a frequency approximately twice as high as the frequency value of the transverse vibrations. The coupling between the transverse and longitudinal oscillations of the spring can es- sentially be modeled through non-linear phenomena. The estimated complex vibra- tion information, observed along the tomographic view-direction is focused along the height (or depth) dimension. This procedure is used to develop high-resolution tomographic underground imaging.

[0005] Here we explain a detailed mathematical formulation of the Doppler Sub- Apertures Model of SAR data. The SAR data belonging to the electromagnetic image are formed through the focusing process that involves the application of a two- dimensional matched filter acting in the range direction and in the azimuth direction. The SLC signal resulting from compression is given by :

Equation (1) represents the focused SAR signal generated by the back-scattered elec- tromagnetic energy of a point target supposed to be stationary. The parameter is showed in Fig. 0.4 (d), and represents the total spatial synthetic aperture of the SAR acquisition (the total orbital distance in the case or space-borne acquisition). The parameter λ is the SAR radio-frequency wavelength. The terms , and are the total chirp and Doppler bandwidths respectively. The total synthetic aperture is equal to L _sa = 2Nd and the azimuth resolution In (1) the parameter identifies the position in range where the maximum of the sinc function is positioned, while in azimuth it is centered around "zero”.

In the case where the peak of the sine function has a nonzero coordinate along the azimuth dimension, equation (1) can be recast as:

In (2) , and are the slant-range, azimuth position of the single-target beam center, into the image coordinates. In this context the DFT is equal to: which recast to a rectangular shape spectrum. In (3) the parameter DFT2 is the two-dimensional digital Fourier transform. The rectangular-shaped spectrum of any SAR image is schematized in Fig. 0.2.

Here we test the above-mentioned strategy that employs Doppler sub-apertures, with the aim to measure target motion. Figure 0.2 represents the used bandwidth allocation strategy. From the single SAR image we calculate the DFT2 which, ac- cording to (3), has a rectangular shape. As it can be seen from Figure 0.2, is the total Doppler band synthesized with the SAR observation, while is the bandwidth not processed from the matched- filter bound- aries, to obtain a sufficient sensitivity to estimate target motions. In this context formula (3) is the focused SAR spectrum, at maximum resolution, thus exploiting the whole band , in accordance with the frequency allocation strategy shown in Figure 0.2, the following range-Doppler sub-apertures large-matrix is con- structed for the master multi-dimensional information: and for the slave, the following large-matrix is presented:

In (4), and (5) the term N _D is the number of Doppler sub-aperture refocused SAR images, extrapolated from the single full-bandwidth (chirp and Doppler) SAR data. The explanation of the chirp-Doppler sub-aperture strategy, represented in Figure 0.2 is the following: master and slave sub-bands are generated by focusing the SAR image, where the matched-filter is set to exploit a range-azimuth bandwidth equal to - The not-processed azimuth bandwidths are divided into N _D equally-distributed bandwidths steps, where N _D represents the sampling-rate of the mechanical wave existing on the Earth that we are observing digitally. At this point N _D rigid shifts of the master-slave system are made along the azimuth bandwidth domain, so as to populate the entire row of equations (4) and (5). The process is repeated N _D times for each shift in azimuth, in fact, Figures 0.2 (1), (2), and (3) represent the azimuth frequency variation strategy when the Doppler bandwidth is located at a distance of N _D . At each Doppler frequency shift every element of (4) and (5) is populated.

Doppler Sub- Aperture Strategy

The decomposition of the SAR data into Doppler sub apertures which is formalized in this subsection, is performed starting from the spectral representation of the fo- cused SAR data. To this end, notice that the generic i— th azimuth sub-aperture 2-dimensional DFT of (2) is given by:

From the last equation, it turns out that a single point stationary target has a two- dimensional rectangular nature with total length proportional to the range-azimuth bandwidths respectively. The phase term exp is due to the sine function dislocation in range and azimuth when the SLC SAR data are considered. In the SAR, the movement of a point target with velocity in both range and azimuth direction is immediately warned by the focusing process, resulting in the following anomalies recasting on the focused SAR image domain: • azimuth displacement in the presence of target constant range velocity;

• azimuth smearing in the presence of target azimuth velocity or target range accelerations;

• range-walking phenomenon, visible as range defocusing, in the presence of tar- get range speed, backscattered energy can be detected over one or more range resolution cells.

In practical cases, the backscattered energy from moving targets is distributed over several range-azimuth resolution cells. As a matter of fact, considering the point-like target T ₁ (of Figure 0.1) that is moving with velocity v _t whose range-azimuth and acceleration components are {v _r , v _a }, and {a _r , a _a }, respectively, then it is possible to highlight

Considering the following Taylor expansion: and that and can be written in the following form: The term can be neglected and by approximating equation (10) can be written Like:

Recasting (11) in terms of x = Vt, we obtain where:

• (due to range velocity);

• (due to range acceleration);

• (due to azimuth velocity).

Thus, the above terms modify the received signal, as shown in , and should be taken into account in equation (6).

[0006] Here we explain how complex vibrational raw data are focused to obtain a tomogram. Data are composed by a temporal expression of . Figure 0.4 (c, d) illustrates the oscillating model in the Euclidean space-time coordinates (x,y,z,t), where the satellite motion has been purified from any orbital distortions, so that the geometric parameters used to perform the tomographic focusing can be rigorously understood. From Figure 0.4 (c), L is the length of the spring when it is at its maximum tension, while L ₀ is its length when no mass is present. Finally, the spring has been considered to have an elastic constant equal to ξ. The vibrational force applied to the mass m1 of Figure 0.4 (c) is equal to : If (13) is expanded in the following series: where a precise approximation of (14) is the following cubic restoring force:

Considering (15), the non-linearity dominates when L ≈ L ₀ . If we define: and

Considering (15) we have:

If we consider damping and forcing (18) is modified as: where f(ω t) is the forcing term and A is the damping coefficient. If non-linearity of (19) is sufficiently low, it can be reduced into the following two-degree-of-freedom linear harmonic oscillator:

In (20) {a, b} are the instantaneous shifts estimated by the coregistrator. The harmonic oscillator (20) contains an approximation of the displacement parameters estimated by (12). According to Figure 0.4 (d) the vector representation of k samples of the time-domain function (20) consisting in the following multi-frequency data input is considered:

The steering matrix A(z) = [a(z _min ), . . . , a(z _MAx )], ∈ C ^kxF contains the phase in- formation of the Doppler frequency variation of the sub-aperture strategy, associated to a source located at the elevation position z ∈ {z _min , Z _MAX } , where is the i— th orthogonal baseline which is vis- ible in Figure 0.4 (d), and is the i-th slant-range distance. The standard sonic tomographic model is given by the following relation:

Y = A(K _z , z)h(z). (23) where in (23) h(z) ∈ C ^1xF , inverting (23) the following tomographic solution is finally found:

In (24) the steering matrix A(K _Z , z) represents the best approximation of a matrix operator performing the Digital Fourier Transform (DFT) of Y. The tomographic image h(z), which represents the spectrum of A(K _z , z), is obtained by doing pulse compression.

The tomographic resolution is equal to where λ is the sound wavelength over the earth, R is the slant range, and A is the orbit aperture considered in the tomographic synthesis, in other words A is proportional to the Doppler bandwidth used to synthesize the sub-apertures. The maximum tomographic resolution obtain- able using this SLC data, synthesized at 22 kHz, is as follows. Considering an average speed of propagation of the seismic waves of about and the maximum observable frequency of investigation (approximately 22000 Hz), the wavelength of these vibrations is equal to about . Considering the above pa- rameters, extending the tomography to the maximum orbital aperture equal to half the total length of the orbit, therefore about 75000 m, with R = 650000m the tomo- graphic resolution is equal to . This is the tomographic resolution set to calculate all the results estimated by the invention.

[0007] Here the computational steps required to perform seismic tomography are explained. This subsection is proposed to explain well the computational scheme to ensure the repeatability of the experiments. The computational scheme, consisting of 11 interconnected blocks, is depicted in Figure 0.5. Block number 1 contains the SAR image in SLC format, processed to obtain tomography, while block number 2 represents the two-dimensional DFT (DFT2) operator. The input of computational stages 3 and 4 represent the copy of the DFT2 of 1, and therefore contains the same data, having a common source.

Computational stage 3 is programmed to generate the sub-aperture represented by the black squared Doppler spectrum visible in Figure 0.2, while block 4 is programmed to generate the blue squared spectrum, also visible in Figure 0.2. Computational blocks 5 and 6 perform the inverse DFT2 (IDFT2) which is used to return to a lower azimuth resolution SLC SAR image. Computational block 7 performs pixel- tracking for all those pixels for which tomography needs to be trained. The result of pixel-tracking is the generation of complex vectors that will form block 8 (non- computational) which represents the raw tomographic complex data that must be focused in elevation or depth but in the orthogonal dimension of the slant range. The focus of the raw tomographic signal is made by computational block 9, representing the DFT mathematical operator. Block 10 (also this is not computational) represents the focused tomographic image. Finally, block 11 performs the geocoding of the tomogram, using a three-dimensional geographic coordinates reference system.

EXPERIMENTAL RESULTS

[0008] For a detailed and practical understanding of the present invention, some preferred sample cases, provided purely as explanatory and non-limiting examples, will now be illustrated with reference to the annexed drawings (not to scale), wherein:

Fig. 0.6 is a schematic illustration of the acquisition geometry currently used for acquiring SAR images of one single area of the surface of Earth in order to reconstruct the tomographic structure of the internal matter of a volcano;

Fig. 0.7 is a schematic illustration of the acquisition geometry currently used for acquiring SAR images of one single area of the surface of Earth in order to reconstruct the tomographic structure of the internal matter. This figure shows a void space located below the Earth’s surface and a pyramid-shaped superstructure extending above the Earth’s surface;

Fig. 0.8 (a,b) represents the tomographic reconstruction of what was observed by the SARI satellite in Fig. 0.8. Fig. 0.8 (a) is the magnitude of the complex raw vibrational data representing the sound information in the time domain, while Fig. 0.8 (b) represents the compressed (focused) vibrational information (in magnitude), thus it is the tomogram, that is the matter inside, in this case below ground. For both figures, the axes are not to scale and represent range and depth distance, in any metric reference.

Fig. 0.8 (a,b) represents the tomographic reconstruction of what was observed by the SAR2 satellite in Fig. 0.8. Fig. 0.8 (a) is the magnitude of the complex raw vibrational data representing the sound information in the time domain, while Fig. 0.5 (b) represents the compressed (focused) vibration information (in magnitude) , thus it is the tomogram, that is the matter inside, in this case the pyramid which is built above the ground. For both figures, the axes are not to scale and represent range and depth distance, in any metric reference.

Bibliography

[1] J. C. Curlander and R. N. McDonough, Synthetic aperture radar: systems and signal processing. 1991. 7

[2] R. K. Raney, “Synthetic aperture imaging radar and moving targets,” IEEE Transactions on Aerospace and Electronic Systems, vol. AES-7, no. 3, pp. 499- 505, 1971. 11, 12

[3] N. B. Tufillaro, “Nonlinear and chaotic string vibrations,” American Journal of Physics, vol. 57, no. 5, pp. 408-414, 1989. 12