WO2015008310A1

WO2015008310A1 - Method for filtering of interferometric data acquired by synthetic aperture radar (sar)

Info

Publication number: WO2015008310A1
Application number: PCT/IT2014/000185
Authority: WO
Inventors: Gianfranco Fornaro; Antonio PACIULLO; Diego REALE; Simona VERDE
Original assignee: Consiglio Nazionale Delle Ricerche
Priority date: 2013-07-19
Filing date: 2014-07-11
Publication date: 2015-01-22
Also published as: ITRM20130426A1; EP3022582A1

Abstract

The present invention concerns a Method for filtering interferometric data acquired by Synthetic Aperture Radar (SAR), comprising the following steps: (A) acquiring a plurality of images N by means of at least one sensor (A₁,..., A_M), which runs a plurality of orbits O_n, with n=1,..., N, in which the sensor (A₁,..., A_M) carries out a plurality of detections of at least one area t _n, each detection during each orbit O_n being carried out in time references t _n, with η=1...., N, emitting a radiation having a predetermined wavelength λ; so as to obtain a multipass data set comprising a set of N images of said at least one area, each of said images being composed of a plurality of pixels P in which targets may be present even possibly interfering due to effects of geometric distortions in the presence of structures with vertical development, said images being geometrically recorded compared to a reference image, to which a reference orbit is associated; (B) determining, for each pixel P of each of said N images, a signal column vector x(P) of length N, constituted by x_n(P) signals, n = 1..., N detected by said sensor (A₁,,..., A_M) in each of said orbits O_n; (C) determining, for each pixel P of each of said images of said set of N images, a covariance matrix associated with said signal vector \(P); (D) determining the components of backscattering contributions of said signal vector x(P) by decomposing said covariance matrix, each backscattering contribution ŷ_κ having a value of the square of its norm λ_κ, with k =1,....,κ < N, and said decomposition being obtained by at least one independent components separation method, such as a method based on the identification of orthogonal components (PCA) and/or a method based on the identification of scattered components, said backscattering contributions ŷ_κ being obtained from the decomposition into eigenvectors and eigenvalues of said covariance matrix by using the following relations:(see formula 1), in which the square of the norm Ak is the k-th ordered eigenvalue and uk is the corresponding eigenvector; and (E) arranging said backscattering contributions ŷ_κ in descending order according to said square of the norm λ and selecting a number K of said significant backscattering contributions ŷ_κ, according to the higher norm with respect to at least one threshold, as a threshold set with respect to the norm λ of the main backscattering contribution ŷ,. The present invention also concern a synthetic aperture radar remote sensing system.

Description

Method for filtering of interferometric data acquired by Synthetic Aperture

Radar (SAR

***_**

The present invention relates to a method for filtering of interferometric data acquired by Synthetic Aperture Radar (SAR).

More specifically, the invention concerns a method for filtering data from Synthetic Aperture Radar (SAR) interferometry, also called CAESAR (Component Analysis and Extraction Sinthetic Aperture Radar) acquired on the same area with angular and possibly temporal diversity, which, jointly using SAR tomography techniques (see F. Lombardini, PI, 2007 A 12 EP/08709820, 8 in Feb. 2007; G. Fornaro, F. Serafino, F. Soldovieri, European patent application EP2017647A1 ) and the principal component analysis (principal Component Analysis - PCA) techniques, it is capable to identify and to separate in each pixel of the image obtained from a SAR radar the significant backscatter mechanisms (backscattering), possibly interfering with each other,

As it is well known, currently synthetic aperture radar (SAR) are widely used in the context of remote sensing. In fact, the advantages of this detection technology are known, which allows to obtain microwave high-resolution images by real antennas having dimensions equippable on board of satellites or aircraft.

In fact, microwave interactions with water molecules are very limited in terms of signal attenuation. Therefore, producing images both day and night, and in all weather conditions, is possible by means of this technology.

One of the most important characteristics of SAR radar is that of being a coherent sensor. In fact, images amplitude is related to the targets capability to backscatter incident radiation, while phase is sensitive, on wavelength scale (centimeters), to the distance of the object from the radar.

Currently, SAR sensing technology is a powerful tool for continuous monitoring of dynamic processes on the Earth's surface.

In the last decade, the remote sensing field has seen a rapid growth and the programs of the major international space agencies devoted to put into orbit SAR sensors increases (e.g. COSMOSKYMED and COSMO/SKYMED II generation for Italy, TerraSAR-X and Tandem-X for Germany, Sentinel-I for the European Space Agency), Regularly acquired data archives (multipass), in subsequent times (multitempora!) on repeated orbits (multi-angle) are available, which are processed with coherent combination (interferometric) complex techniques. In general, these techniques allow an accurate localization of the targets, and then the measurement of the soil topography and the ability to monitor the movements with centimetric/millirnetric accuracies.

Differential interferometry (DlnSAR) techniques are specifically used for monitoring the soil surface deformation. Said DlnSAR techniques are based mainly on the use of the signal phase backscattered from the scene or area illuminated by the sensor. Currently, measurements of historical series are provided with processing techniques oriented to the observation of targets distributed or with point techniques, compared to the spatial resolution of the radar system.

Algorithms belong to the first class that, in order to limit the changes effects (decorrelations) of radar targets response with respect to view angle (angular or spatial decorrelation) and time (temporal decorrelation) variations, restrict the analysis to a selection of interferograms obtained by strict constraints on the distance (baseline) and on the acquisitions time in the construction of interferometric pairs. For example, the technique Small BAseline Subset - SBAS, described in Berardino et al. "A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms", IEEE Trans. Geosci. Remote Sens., Vol. 40 (1 1 ), pp. 2375-2383 2002, assuming a diffusion mechanism (scattering) spatially distributed, contrasts decorrelation phenomena both limiting the lines bases or baselines, as well as by using average spatial operations (mu!tilook). This latter operation is carried out at the expense of a degradation of the final product resolution. The SBAS technique allows to analyze low-resolution deformations over broad areas but turns out to be less appropriate for the reconstruction and analysis of single structures deformation at the highest available resolution of the radar.

In a dual way, techniques oriented to the monitoring of the deformations of localized targets corresponding to strong scatterers, known in the literature with the name of Persistent or Permanent scatterers (PS), are based on the development of SAR data at maximum spatial resolution available (for example persistent Scatterer Interferometry - PSI in A. Ferretti, C. Prati, and F. Rocca, "Nonlinear subsidence rate estimation using the permanent scatterers in differential SAR interferometry^" !EEE Trans. Geosci. Remote Sens., vol. 38 (5), pp. 2202- 2212, 2000).

Said techniques assume that the electromagnetic response of a given object in a given resolution cell is correlated in all available acquisitions, or better t remains consistent (persists) throughout the observation period and in case of variation of the angle diversity of the sensor observation. This assumption is valid only for a limited number of objects (i.e. the Persistent Scatterers), on which applying specific algorithms (see also the European patent no. EP-1 183551 and the Italian patent application MO2007A000363) is possible.

Unlike interferometry processing techniques oriented to the distributed targets observation (such as SBAS), the Persistent Scatterers Interferometry techniques or PSI do not bind angular or temporal differences in the nterferograms generation. The main drawbacks of this technique are associated to the use of only SAR data phase and to the limitation of the analysis of only point targets, for which it is assumed a scattering of the "localized and dominant" type, which is maintained also correlated on high baseline.

In real cases, a vast majority of pixels is still influenced by the effects of decorrelation. The degree of correlation of a pixel may also vary significantly from image to image and then from interferogram to interferogram.

Recently, a technique has been developed, known as SQUEESAR (see international patent applications no. WO201 1/003836 A1 and EP2010/059494, A. Ferretti, A. Fumagalli, C. Novali, C. Prati, F. Rocca, and A. Rucci, "A new algorithm for processing interferometric data-stacks: SqueeSAR ^", IEEE Trans. Geosci. Remote Sens., vol. 49 (9), pp. 3460- 3470, Sep. 201 1), which extends the PSI technique, allowing the monitoring of distributed scatterers, typical of rural areas. This technique, using the coherence measurements estimated by SAR data, uses an iterative procedure to estimate (for each possible pair of interferometric data set obtained from a multiangular/multitemporal SAR radar) the scattering mechanism phase equivalent to that of a persistent scatterer.

A limit of all the analyzed techniques is not to consider the presence of multiple scattering mechanisms in the individual pixels of acquired radar images, in fact, because of radar view perspective distortions, the prevalent vertical development complex scenarios, such as those urbanized, frequently present overlapping phenomena (layover) of the responses from different targets on the soii (e.g., the ground and the buildings wails and roofs, bridges and the structural components of the soil, etc.).

Therefore, due to the above, DlnSAR data processing technology does not allow to separate the contributions linked to distinct scatterers that interfere in the same cell of the SAR image resolution, i.e. to solve the layover phenomenon.

Differential SAR Tomography (see Fornaro et al. Tour-dimensional SAR imaging for height estimation and monitoring of single and double scatterers" IEEE Trans. Geosci. Remote Sens., Vol. 47 (1 ), p. 224-237 , 2009), exploiting the difference in view angle, uses a processing technique aimed at the realization of (synthetic) antennas also in the direction orthogonal to flight and radar view direction. This technique allows the creation (through a digital processing) of a "radar scanner" with a narrow lighting beam (the beam) which can be (digitally) oriented to carry out high accuracy investigation of individual structures (e.g. buildings). Unlike SBAS and PSI techniques, tomography s based on the development of the entire complex information, considering, then, both the module and the phase of the SAR data. Scientific literature studies (A. De Maio, G. Fornaro, A. Pauciullo, "Detection of Single Scatterers in Multidimensional SAR Imaging", IEEE Trans. Geosci. Remote Sens., Vol. 47 (7), pp. 2284 - 2997, July 2009) have shown that, with reference to individual persistent targets, this technique allows to improve the performance compared to the PSI algorithms type in terms of detection probability and estimation accuracy of the parameters of the location and displacement of the persistent scatterers, just because it uses the module and phase of the radar images. In addition to the above, SAR tomography offers the advantage of allowing the separation of any contributions interfering in the same resolution ceil due to layover phenomenon.

As for PSI technique, SAR tomography requires, however, at a preliminary stage, a compensation operation of the delays associated with the radiation propagation in the atmosphere. These delays can be estimated and compensated in the data using the traditional techniques of differential interferometry operating at low resolution (e.g., SBAS, SqueeSAR), which, by adopting multilooking operations, improve the spatial coverage.

Prior art also includes the article of Fabrizio Lombardmi et a!. "Superresoiution Differential Tomography: Experiments on Identification of Multiple Scatterers in Spaceborne SAR Data", IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 4, 1 April 2012, pages 1 1 17-1129, which describes a processing method for data acquired by Adaptive Differential Tomography synthetic aperture radar (SAR) for scatterers super-resolution.

SAR tomography is a well established multibaseline-muititemporal technique, that allows to solve more moving scatterers placed at different heights in the same pixel of the SAR image. SAR tomography is applied to SAR data, after compensation of propagation effect in the atmosphere (generally referred to as the calibration phase), in this document, already calibrated data and a model of linear displacement over time are assumed available.

This solution, however, does not allow to identify and separate backscatter dominant mechanisms (main components) in each pixel of SAR image, possibly interfering in the same resolution cell upstream of the processing chain (for example at the level of interferogram generation) and in a phase, which is intended, in fact, to estimate the atmospheric delay and deformation, even nonlinear in a gross scale.

This document, being a tomographic technique, uses a specific model for the scatterer response, known as the steering vector, parameterized with respect to the target elevation and deformation average speed. As a result, in the described technique, a structured separation of multiple components is performed.

In light of the above, it is, therefore, object of the present invention to propose a method for processing data collected by synthetic aperture radar (Synthetic Aperture Radar - SAR) that, operating with multilook data, allows to identify and to separate, in each pixel common to all images captured by SAR radar, dominant backscattering mechanisms (main components), even possibly interfering in the observed scene, jointly exploiting the PCA technique and SAR Tomography.

A further object of the present invention consists, therefore, in developing a technique capable of filtering the backscattered dominant contribution, whether punctiform or distributed, with respect to the resolution ceil, thus counteracting the effects related to temporal and angular decorrelation phenomena and improving, then the quality of the interferometric products. A further scope of the present invention is, using principal components analysis technique (Principal Component Analysis - PCA) in combination with tomography, that of identifying, distinguishing, separating multiple scattering mechanisms, interfering in the same pixel, thus allowing the selection of interferings according to power or altitude ordering criteria, and, thereby, an efficient and effective method to solve the problem of undesired layover, especially prevalent in urban scenarios.

It is also object of the present invention to filter out the contributions associated with the dominant target, possibly even interfering with secondary and higher order contributions, in stacks of interferometric SAR acquisitions.

It is therefore object of the present invention a method for filtering interferometric data acquired by Synthetic Aperture Radar (SAR), comprising the following steps: (A) acquiring a plurality of images N by means of at least one sensor (Ai , AM), which runs a plurality of orbits 0_R, with n=1 , ... , N , in which the sensor (A-i , ... , A_M) carries out a plurality of detections of at least one area t_n , each detection during each orbit O_N being carried out in time references t_n , with n- , , . , , Ν , emitting a radiation having a predetermined wavelength λ , so as to obtain a multipass data set comprising a set of /V images of said at least one area, each of said images being composed of a plurality of pixels P in which targets may be present even possibly interfering due to effects of geometric distortions in the presence of structures with vertical development, said images being geometrically recorded compared to a reference image, to which a reference orbit is associated; (B) determining, for each pixel P of each of said N images, a signal column vector x(P) of length Λ^? , constituted by x_n (P) signals, n - \,..., N detected by said sensor (A-i , . . . , A_M) in each of said orbits O_r.; (C) determining, for each pixel P of each of said images of said set of N images, a covariance matrix associated with said signal vector \(P) ; (D) determining the components of backscattering contributions "^!k of said signal vector \( P) by decomposing said covariance matrix, each backscattering contribution ^^k having a value of the square of its norm Α,. , with k I A . \ , and said decomposition being obtained by at least one independent components separation method, such as a method based on the identification of orthogonal components (PCA) and/or a method based on the identification of scattered components, said backscattering contributions ^^k being obtained from the decomposition into eigenvectors and eigenvalues of said covariance matrix by using the following relations: in which the square of the norm A_k is the k-th ordered eigenvalue and u , is the corresponding eigenvector; and (E) arranging said backscattering contributions in descending order according to said square of the norm /;, and selecting a number K of said significant backscattering contributions ^γ* , according to the higher norm with respect to at least one threshold, as a threshold set with respect to the norm A, of the main backscattering contribution γ, .

Always according to the invention, said step (C) could comprise the step of determining, for each pixel /' of said set of N images, a set of pixels w^'m(P) constituted by said pixel P and a plurality of pixels Q spatially close to it, selected using at least one criteria of statistical similarity between the plurality of signal vectors

and the signal vector %(P) , and providing the estimate of the covariance matrix of said vector

\{ P ) as: C ^ *( χ^■'^''' (^'ø where N_p is the number of pixels constituting said set of pixels win(P) .

Advantageously according to the invention, said step (A) could comprise the following substeps: (A.1 ) focusing in azimuth and range variables, said detected multipass data set; (A.2) aligning at spatial level, with respect to said azimuth and range variables, A⁷ - 1 of said images with respect to said reference image; (A.3) determining, for each pixel P of each image, the geometrical distances of each of said at least one sensor (A^ , , . , A_M) in each of said orbits by means of a reference altimeter digital map; and (A.4) subtracting, for each pixel P of each image of said set of images, the phase corresponding to the calculated distances of said sensor (A , , . . .AM) in each of said orbits O_N .

Further according to the invention, said methos could further comprise the following step: (F) ordering said possible targets that interfere in the pixel P according to the elevation (s).

Always according to the invention, said step (F) could comprise the following substeps: (F.1 ) determining, for each estimates pair of said backscatter contributions _k and y_m with k 1,..., A^" and m - \,..., K , the estimate of the difference of elevation values A associated with the corresponding targets; and (F.2) ordering said targets according to the sign of said estimation of the difference of elevation values Av .

Always according to the invention, said step (F.1 ) could comprise the following substeps: (F. .1 ) determining the baselines b, _n , n = 1 , N, for each of said time references t„ , for said at least one sensor (Ai , A_M) on each of said orbits O_n, said baselines h , _n being the parallax distances of each single acquisition with respect to said reference orbit of said step (A); (F.1 .2) calculating the column vector ς of length N constituted by the spatial frequencies _n , n = 1 , N , corresponding to acquisitions and obtained as ζ_Λ = 2b _r,/(Ar) (F.1 .3) calculating the column vector κ of length N constituted by the temporal frequencies k_n , n = 1 , N, corresponding to the acquisitions and obtained as k_n = 2 t ; (F.1 .4) calculating the direction vector aO, v) , as: aO, v) = exp(j2ngs + jlm s) where (s, v) is an arbitrary point of the elevation plane ( s ) and speed ( v ); (F.1 .5) calculating the beat vector Au_km between the eigenvectors u . and u .„ corresponding to the two backscatfering contributions at issue obtained as Διι^, - II , o u'. , where the symbol indicates the Hadamard product and the symbol "*" the complex conjugation operator; (F.1 .6) calculating the Euclidean norm scalar product between said direction vector and said beat vector; (F.1 .7) providing an elevation search interval, preferably symmetrica! with respect to zero [-½ ,½ ] ; (F.1.8) providing a speed search interval, preferably symmetrical with respect to zero [-v_A, . v_M ] ;

(F.1 .9) maximizing the module of said scalar product with respect to the elevation and speed variables pair (s, v) , searching for the solution in the domain defined by the Cartesian product [~ί_Μ,¾ ] [™¾, _Μ] by a discretization process; and (F.1 .10) taking the argument on the variable Λ· of said maximization, as estimation AS„

Still according to the invention, said methos could comprise, after said step (E), the following steps: (G) generating, for each pixel P , an appropriate combination of N images, obtained by selecting one of the components obtained by said steps (A) to (E), and applying at least one multipass data differential interferometric processing technique, so as to separate the signal associated with the surface deformations from the signals associated with the atmospheric delays.

Further according to the invention, said said step (G) could comprise the following steps: (H.1 ) generating for each pixel P an appropriate combination of images, obtained by selecting one of the components obtained by said steps (A) to (E); (H.2) calibrating the data of said step (A.4) by the low resolution atmosphere and deformation signals estimated in step (G); and (H.3) estimating the topography and the deformations of said area by applying said steps from (F. L I ) to (F.1 .10), replacing said vectors of said step (F.1 .5) with the calibrated data according to said step (H .2).

Advantageously according to the invention, said method could comprise a plurality of sensors (A< , ... , A_M), and in that said at least one sensor can be constituted by one or more antennas.

Always according to the invention, each of said sensors (Ai AM) could acquire said data simultaneously with the other in each of said orbits O_n, for n = 1 , N, so that t_n} =... = t_rM .

It is further object of the present invention a synthetic aperture radar remote sensing system, comprising at least one sensor (Ai , . . . , AM), adapted to emit a radiation towards at least one area and to receive the return radiation, at least one transport or moving means, such as a satellite or an airplane, on which said at least one sensor is placed, said transport means performing a plurality of passages on said area according to a plurality of trajectories, data storage means, connected with said at least one sensor (A^ , ... , AM), to store the data it has collected, and means for processing the data detected by said at least one sensor (Ai , ... , AM), connected with said data storage means, characterized in that said processing means run on said data detected by said at least one sensor (Ai , ... , AM) the method for processing the data detected by synthetic aperture radar according to the method as defined in claims 1 -1 1 .

Always according to the invention, said system could comprise a plurality of sensors (A\ , ... , AM), placed on a corresponding plurality of transport means.

Still according to the invention, said sensor could be an antenna (Ai , ... , AM) .

The present invention will be now described, for illustrative but not limitative purposes, according to its preferred embodiments, with particular reference to the figures of the enclosed drawings, wherein:

figure 1 shows the geometry of the acquisition system of the system for filtering interferograms obtained from data acquired by Synthetic Aperture Radar (SAR) according to the present invention;

figure 2 shows a synthetic block diagram of the method for the filtering interferograms according to the present invention; figure 3 shows the amplitude of a radar image corresponding to the processed scene, relating to an area with vegetation;

figures 4a and 4b show a comparison between the interferograms related to the area in figure 3;

figure 5 shows the amplitude of a radar image relating to an area characterized by vertical structures; and

figures 6a-6c show a comparison between interferograms relating to the area in figure 5.

!n the various figures, similar parts will be indicated by the same reference numbers.

Referring to figure 1 , the geometry of the acquisition system is shown in a section orthogonal with respect to the flight trajectory of the sensor (for constant azimuth) i.e. an antenna A, for which, for any fixed range r, and assuming to neglect diffraction effects and mutual interaction between the targets, the relationship between the distribution of the backscattering coefficient in elevation s, which will be called y(s) , and the data x„ acquired from one or more antennas, in which x„ is the signal for a pixel on the n-th image, after an appropriate amplitude and phase geometric calibration pre-processing (see G. Fornaro, F. Serafino, F . Soldovieri, European patent application EP20 7647A1 ) is a Fourier transform type:

where frequencies ζ_η , n = 1 , N, are equal to _n - 2b_in /(Ar) , λ is the wavelength of the incident radiation and the quantities b_{l n} , n = 1 , N, are the parallax displacements , i.e. the components of the distance vector of the acquisitions from a reference orbit, orthogonal to the view direction of the scene from the reference orbit. These components are commonly referred to as base lines or baselines.

The operator that links the data and the unknowns defined by equation [1] is of linear type and semi-discrete and can be reversed with different techniques, in order to obtain a 3D reconstruction (tomographic approach) of the backscattering profile. The most commonly used technique, as mentioned in the preamble, is referred to as "beam forming" (Beam-Forming) and it is based on the application of the operator added (conjugate transpose of the matrix operator resulting from the discretization of equation [1]).

Where it is desired to reduce the effects of the disturbs and of the noise, 3D reconstruction should be conducted with simuitaneous acquisitions obtained by alignments of antennas.

Examples of multiple antennas simultaneous acquisitions are provided by muitistatic systems on-board of aerial platform and by "Tandem-X" system on the satellite platform.

Monostatic satellite systems, i.e. constituted by a single antenna, can "synthesize" an alignment of antennas due to the characteristic of their platform to repeat its orbit. By real data provided by antenna alignments synthesized in following steps, it is known and it has been shown that, on sufficiently stable targets (typically anthropic structures), it is possible to recover the three-dimensionality of the investigated scene.

If acquisitions are synthesized in later times, it is possible to combine tomographic and differential approaches. In particular, the problem of locating and monitoring targets from multipass and multiview SAR data can be seen as a 4D imaging problem (elevation and speed, in addition to azimuth and range), which consists, then, in the linear inversion of a two-dimensional Fourier transform operator i.e. 2D (or higher order, where it is wished to take into account further components of mouldable deformation, such as those thermal) for each pixel, the azimuth and range.

The known concept of 3D imaging by SAR multiview data is, therefore, extended to the temporal dimension, and the variable speed (scaled by an appropriate factor inversely proportional to the wavelength) assuming the meaning of spectral frequency that describes the movement of the target to the soil. In this way, it is possible to separate the

movements, even nonlinear, of possibly interfering target in the same pixels azimuth and range at different elevations (or altitudes).

In the most common case of a single detection antenna, indicating with t, , n = 1 N, the instants in which baselines b_ln, n ~ 1 , ... , N, are synthesized, where the time reference is fixed with respect to antenna (synthetic), to which also base lines relate overlaying the harmonics of the speed corresponding to deformation mechanisms possibly nonlinear as described in (G. Fornaro, F. Serafino. F. Soldovieri, RM2007A000399 - EP2017647A1 , July 19, 2007), it is possible to write the following linear relationship (n = 1 ,... ,N): ^x» " j i^{? , y}-

fin κ_η v]ds dv [2] where the frequencies κ„ are given by the ratio 2 t_n(X , which highlights a double Fourier transformation relationship.

Equation [2] represents the relationship that exists between the data x_H , namely between the vector of detected signals of each pixel of the n-th image, after the pre-treatment, and a distribution function y(s, v) of targets in elevation, which we will call distribution function of the backscattering coefficient each one in independent motion and not necessarily uniform, for a system that acquires multitemporal and muitiview data with distributions of times and base lines respectively equal to t_n and h . _n , n = 1 .....N.

The linear and semi-discrete operator of the equation [2] can be completely discretized by considering a grid of points, i = 1 , M_s and j = 1 , M_V! distributed over the domain of interest (integration domain) of the eievation/speed plane. In this way the following matrix problem is obtained:

x - Ay [3] where the vector M=M_sxM_v-dimensiona! γ , of the unknowns Y(s_t , v . ) , i=1 , ....M_s and j=1 ,.. . ,M_v, and the N-dimensionai data x... vector x , n-1 , N, are related by the linear operator represented by N*M dimensional matrix A , which M-th column is given by the direction vector:

a,.. = 'ά($ ,ν,.) = exp /2;r( , + ν )] , [4]

where ς and κ are the N-dimensional vectors of the spatial and temporal frequencies, respectively, m = 1 , M, i = 1 , _s and j = 1 , .... M_y.

The Beam-Forming inversion technique returns the following solution of the matrix problem of the equations system [3]:

where H is the conjugation and transposition operator of the vectors, which represents the projection (scalar product in Euclidean norm) of the data along each direction vector defined in equation [4]. In other words, consistent with the direct mode! of equation [3], the inversion technique assumes that any scattering mechanisms are present along directions of the C^N vector space (with C complex numbers field) structured in accordance with the equation [4]. In particular, the dominant scattering mechanism corresponds to the direction a.,, , that maximizes the module of equation [5] and the relative position in the domain of interest of the elevation/speed plane is given by:

(s, V) , v , )| = arg max a" (xx^h )a„, = arg max P- (s_i ,v . ) [6]

5 where m→ (/, j) indicates the index corresponding to the copy (i, j), and

is the power of the reconstruction of the function at the backscatter point

The strong limitation of this technique (and all of the imaging techniques that have as objective the inversion of the linear problem according to the equation [3]) lies in assuming that the data acquired, following an appropriate pre-processing, follow exactly the model is according to the equation [3].

As it will be described in detail below, any phase errors on the data, due for example to a not perfect compensation of the atmospheric effects, determine a deviation from the model described in equation [3], which basically implies a rotation of the direction vectors of equation [4]. This

20 rotation can not, of course, be taken into account in the reconstruction process, with the consequent performance loss of the inversion technique.

The present invention concerns the application of a separation technique of different scattering mechanisms in the remote sensed data, that is robust with respect to possible deviations of the same data from the 25 model expressed in equation [3].

The block diagram of the solution is shown in figure 2, in which it is highlighted the possibility, by multilook spatial operations, of estimating the covariance matrix, to identify, separate and select, starting from a stack or a set of acquisitions by multiangular/multitemporal SAR radar, one

3 0 or more sets or stacks corresponding to one or more dominant scattering mechanisms interfering in the same pixels due to layover effects on structures with vertical development, which does not require special phase calibration operations as those due to weather delays.

In formal terms, a first step to release the scattering mechanisms 35 from the particular structure of the direction vectors is to avoid the structure equation [4] assumption, seeking the direction associated with all scattering mechanisms between all vectors a C^N with fixed norm (for example, all the direction vectors structured according to the equation [4] have generally norm equal to N ). The elimination of such a structure can not be made in equation [6], because it would lead to the trivial solution a oc x .

Consider, then, the quantity y_a = a^'^ , which represents the projection of the data vector x on the generic vector a e C^N, which is assumed subject to the condition ||a|| ~ i . Consider, in addition, the average power of such projection: P (aa ^H . C_x ) [8]

where C ^ E(XX^H ) is the data covariance matrix, assumed zero mean,

E(.) is the statistical average operator, trace{} is the trace operator of a matrix and (A, B) is the scalar product operator between matrices in Frobenius norm.

Releasing the scattering mechanisms from the particular structure of the vectors in equation [4], it is possible to search the directions associated with the scattering mechanisms among all vectors a e C^N with fixed norm such that the scalar product between the matrix associated with the outer product (corresponding to all interferometric pairs) aa ^,: , and the estimated covariance matrix is above the threshold,

A direct way to get the directions with main components is obtained by the spectral decomposition of the covariance matrix as follows:

V where _k and u . are the k-th eigenvalue and the corresponding eigenvector (with unit norm) of the above matrix.

The eigenvectors set {uj^, , therefore, constitute an orthonormal basis of the vector space C^N. Also, being €.. a Hermitian matrix ( C_T =€_<' ), its eigenvalues are all real and non-negative; for convenience they are assumed ordered in descending way: A≥ A₂≥ ...≥ A_N .

It is shown that, varying a e C^N among unitary norm vectors, the average power in equation [8] is maximum in correspondence of the eigenvector u_{ of C. associated to ther maximum eigenvalue λ and is really equal to A, .

The dominant scattering mechanism, defined as the one to which corresponds the maximum average energy, is thus associated with the eigenvalue-eigenvector pair (Λ,, , ΙΙ, ) , which lies along the direction u, and has average energy Λ, , i.e. γ_χ - _Λ /ζα_; .

Moreover, varying a e C.^v among unitary norm vectors, where C,'^v is the vectorial subspace orthogonal to u, , the average power in equation [8] is maximum in correspondence of the eigenvector u, of C, associated to the second maximum eigenvalue λ₂ and is precisely equal to λ₂ . The second scattering mechanism, defined as the one which corresponds to the maximum average power in the subspace orthogonal to the direction of the dominant mechanism, it is therefore associated to the eigenvalue- eigenvector pair (A,_j , u₂ ) , i.e. lies along the direction u , , and has average energy ? , i.e. f ^~ _v -.· " . .

In general, the k-th scattering mechanism, defined as the one which corresponds to the highest value of the average power in the subspace orthogonal to the subspace identified by the directions of the first k-1 scattering mechanisms, lies along the direction u . and has average energy A_k , i.e. v, =V^ u* -

Note that, differently from tomographic reconstruction algorithms, the selection of scattering mechanisms using covariance matrix spectral decomposition does not require any assumption about the structure of the directions to be identified. This feature allows to apply the above procedure for the selection of scattering mechanisms directly on data focused in range and azimuth, but not yet calibrated, i.e. before carrying out any compensation of undesired phase contributions due to atmospheric effects and possibie deformation of the soil. In fact, the not calibrated data can be expressed as y = ψ ^ χ , where ψ is a complex phase vector (i.e. unitary norm elements), which takes into account of the abobve mentioned undesired contributions and the symbol " ° " indicates the Hadamard product (i.e element by element) between matrices. The corresponding covariance matrix is then:

C_x = ψψ^Λ o - ^ _ι.ψψ o u_A.u -∑l_k (ψ o , )(ψ o u_k )^!i [10]

The above equation shows that a not-perfect phase calibration of the data determines a relative rotation of the eigenvectors of the covariance matrix (directions of the scattering mechanisms), but leaves unaffected the eigenvalues (power scattering mechanisms). Therefore, the application of this procedure, applied to a not-calibrated data, still allows to identify the direction to which the maximum power value is associated considering the following orthogonal subspaces.

In real cases, the covariance matrix of the data has to be estimated from the data itself. For this purpose, the present invention implements the above procedure, where the matrix C , is replaced by an estimate. In particular, called P the processing pixel and indicating explicitly the dependence of the data x( P) , it is considered the estimate:

< ' ∑x(0x(0 [1 1 ] where winfPj is a set, a window, of pixels, that contains the pixel P with cardinality N_P. on which the collected data x(0 (look) are statistically similar to the processed data

. If the looks x(Q) are independent, the estimate of the covariance matrix has rank greater than one, then separating different scattering mechanisms is still possible.

In the following examples, the selection of the pixels statistically similar to the determination of the covariance matrix has been performed using a simple statistical Kolmogorov-Smirnov test, but solutions or different systems are possible, as those based on the non-local filtering methods, used for the reduction of speckle noise in the amplitude images.

Unlike tomographic technique that operates a hypothesis of structuring the response according to a parametric model of the response of the target (typically with the speed and topography parameters), which is valid only downstream of the compensation (phase calibration) of the mismatch contributions of the model (in primis propagation delays in the atmosphere and non-linear deformations), the aspect of the separation of possible interfering contributions in the proposed technique occurs upstream, directly on the focused data and downstream of the images single alignment operation (registration), without assuming a specific structure of the response.

In the proposed technique according to the present invention, this is possible due to a specific use of the eigenvalues and eigenvectors extracted from the covariance matrix of the focused and recorded interferometric acquisitions stack. In fact, while in the classical interferometry techniques the mediated interferograms, which constitute the elements of the covariance matrix, are treated separately according to a processing chain that leads to the estimation of deformation and tropospheric delays, the proposed technique according to the invention for each pair of eigenvalue

(u_k ), obtained by processing the covariance matrix, generates a new acquisitions stack ( ¾u* )_< which corresponds to the contribution of the individual scattering mechanism. From extracted stack interferograms filtered from noise and other possible contributions from significant interfering scattering can then be generated.

In other words, the method according to the invention decomposes the covariance matrix of the original data that contains the interferograms corresponding to the starting acquisitions stack, in matrixes (dyads) containing the interferograms corresponding to the contributions of the dominant, secondary, etc., scattering.

Figures 3-6 show the results obtained through the present invention.

In particular, figure 3 shows the amplitude of a radar image corresponding to the processed scene, relative to an area with vegetation, and thus affected by temporal decorrelation phenomena, in order to facilitate the visualization of figures 4a and 4b.

Figures 4a and 4b show a comparison between interferograms for the area in figure 3. In particular, figure 4a shows a full-resolution original interferogram, in which the phase signal noise is due to decorrelation phenomena. Instead, figure 4b shows an interferogram reconstructed by the present invention, in which dominant backscatter contribution has been filtered. The covariance matrix has been estimated using a uniform spatial average in a set of win(/^J) pixels. It is seen the overall quality improvement of the phase signal.

Figure 5, in order to facilitate the i terpretation of the figures 6a~6c, shows the amplitude of a radar image of the processed scene relative to an urban area characterized by a dense presence of buildings and, more generally, by vertical structures, in which the layover phenomenon is strongly present. Figures 6a-6c show a comparison between interferograms of the area in figure 5.

In particular, figure 6a shows an interferogram averaged spatially adaptively without the application of the interferograms filtering method according to the present invention. Figure 6b shows an interferogram obtained by the decomposition by the interferograms filtering method according to the present invention relative to the backscatter contribution of the soil, i.e. the phase signal constructed by selecting on layover mechanisms the scatterers located at low altitudes and which shows a mitigation of the contributions related to the buildings topography. The distinction between two general scattering mechanisms k and I in lows and highs, is made by detecting the sign (in elevation, s) of the peak corresponding to the reconstruction obtained by applying to the vector data u _k c u (on which the effects of atmospheric delay are intrinsically deleted in the product operation) the known tomographic technique in [2] with vectors structured according to [4]. Finally, figure 6c shows an interferogram obtained by the decomposition by the interferograms filtering method according to the present invention relative to the backscattering contribution from buildings, or the phase signal constructed by selecting on layover mechanisms the scatterers located at high altitudes and showing an over-emphasis of the contributions related to buildings topography.

An advantage of the interferometric data filtering method acquired by Synthetic Aperture Radar according to the invention is allowing to generate filtered interferograms, in which topographical contribution can also be emphasized or de-emphasized,

A further advantage according to the invention, is to perform an "unstructured" separation of the scattering components: it does not use, in fact, a specific model for the target response but it estimates its structure directly from the data using the Principal Components Analysis (PCA). This analysis involves the evaluation of the data covariance matrix, from which the eigenvectors are extracted, automatically identifying the structure of scattering basing, therefore, only on the measures of eigen and mutual power on the array (spatial-ternpora!) corresponding to acquisitions. Moreover, it should be noted that the not-structured separation according to the invention can also be applied to calibrated data as an alternative to traditional SAR Tomography approaches for reconstructing and monitoring of the observed scene in detail scale. In this case, the invention allows to obtain high density of measuring points, at the expense of a slight loss in spatial resolution.

The present invention has been described for illustrative but not limitative purposes, according to its preferred embodiments, but it is to be understood that modifications and/or changes can be introduced by those skilled in the art without departing from the relevant scope as defined in the enclosed claims.

Claims

1 . Method for filtering interferometric data acquired by Synthetic Aperture Radar (SAR), comprising the following steps:

(A) acquiring a plurality of images N by means of at least one sensor (Ai AM), which runs a plurality of orbits O_n, with n=1 , , , , , J¥ , in which the sensor (Ai , ... , A ) carries out a plurality of detections of at least one area t„ , each detection during each orbit O_n being carried out in time references t„ , with η=1 , ... , Λ' , emitting a radiation having a predetermined wavelength λ , so as to obtain a multipass data set comprising a set of N images of said at least one area, each of said images being composed of a plurality of pixels P in which targets may be present even possibly interfering due to effects of geometric distortions in the presence of structures with vertical development, said images being geometrically recorded compared to a reference image, to which a reference orbit is associated;

(B) determining, for each pixel P of each of said N images, a signal column vector %(P) of length N , constituted by x P) signals, n = \,..., N detected by said sensor (Ai , ... , A_M) in each of said orbits O_n;

(C) determining, for each pixel P of each of said images of said set of N images, a covariance matrix associated with said signal vector x( i ;

(D) determining the components of backscattering contributions i of said signal vector x(P) by decomposing said covariance matrix, each backscattering contribution ^γ¾ having a value of the square of its norm A_t , with k = ι,.,., κ < N , and said decomposition being obtained by at least one independent components separation method, such as a method based on the identification of orthogonal components (PCA) and/or a method based on the identification of scattered components, said backscattering contributions ^γ·^· being obtained from the decomposition into eigenvectors and eigenvalues of said covariance matrix by using the following relations: y_k = ¾ u_t , k = \,..., K ,

in which the square of the norm X_k is the k-th ordered eigenvalue and u, is the corresponding eigenvector; and

(E) arranging said backscattering contributions ^^k in descending order according to said square of the norm l_k and selecting a number K of said significant backscattering contributions ^γ* , according to the higher norm with respect to at least one threshold, as a threshold set with respect to the norm A, of the main backscattering contribution γ, .

2, Method according to claim 1 , characterized in that said step (C) comprises the step of determining, for each pixel /' of said set of N images, a set of pixels win( P) constituted by said pixel P and a plurality of pixels O spatially close to it, selected using at least one criteria of statistical similarity between the plurality of signal vectors x(0 and the signal vector x(P) , and providing the estimate of the covariance matrix of said vector x(P) as:

where N_t> is the number of pixels constituting said set of pixels win( ) .

3, Method according to anyone of the preceding claims, characterized in that said step (A) comprises the following substeps:

(A.1 ) focusing in azimuth and range variables, said detected multipass data set;

(A.2) aligning at spatial level, with respect to said azimuth and range variables, N - 1 of said images with respect to said reference image;

(A.3) determining, for each pixel P of each image, the geometrical distances of each of said at least one sensor (A-i , ... , Ay) in each of said orbits by means of a reference altimeter digital map; and

(A.4) subtracting, for each pixel /' of each image of said set of images, the phase corresponding to the calculated distances of said sensor (Α·, , ... ,A_M) in each of said orbits O_n.

4, Method according to anyone of the preceding claims, characterized in that it further comprises the following step:

(F) ordering said possible targets that interfere in the pixel /' according to the elevation (s).

5, Method according to claim 4, characterized in that said step (F) comprises the following substeps:

(F.1 ) determining, for each estimates pair of said backscatter contributions γ,. and _m with k - 1 K and m - I,..., K , the estimate of the difference of elevation values As associated with the corresponding targets; and

(F.2) ordering said targets according to the sign of said estimation of the difference of elevation values As .

6, Method according to claim 5, characterized in that said step (F.1 ) comprises the following substeps:

(F.1 .1) determining the baselines b_{l n} , n = 1 , N, for each of said time references t_n , for said at least one sensor (A-i , AM) on each of said orbits O_n, said baselines b _n being the parallax distances of each single acquisition with respect to said reference orbit of said step (A);

(F.1 .2) calculating the column vector ς of length N constituted by the spatial frequencies ζ.. . n - 1 , N, corresponding to acquisitions and obtained as ζ_η ~ 2b_L (Xr);

(F.1 .3) calculating the column vector κ of length N constituted by the temporal frequencies k„ , n = 1 , N, corresponding to the acquisitions and obtained as k_n ~ 2 t X ;

(F. I . ) calculating the direction vector a(.s, v) , as:

a(s, v) = exp(/2/r¾s + jlma)

where {s, v) is an arbitrary point of the elevation plane ( s ) and speed ( v );

(F. I .5) calculating the beat vector Au.„, between the eigenvectors u and u.„ corresponding to the two backscattering contributions at issue obtained as Au_hll = u, « u^* _; , where the symbol " ° " indicates the Hadamard product and the symbol "^*" the complex conjugation operator;

(F.1.

6) calculating the Euclidean norm scalar product between said direction vector and said beat vector;

(F.1.7) providing an elevation search interval, preferably symmetrical with respect to zero [-¾, ¾ ] :

(F.1.8) providing a speed search interval, preferably symmetrical with respect to zero [-v_w , v_M ] ;

(F.1 .9) maximizing the module of said scalar product with respect to the elevation and speed variables pair (s, v) , searching for the solution in the domain defined by the Cartesian product [-¾ , s_M ] ^χ [~~v_M , v_M j by a discretization process; and

(F.1 .10) taking the argument on the variable s of said maximization, as estimation Δ? .

7. Method according to anyone of the preceding claims, characterized in that it comprises, after said step (E), the following steps:

(G) generating, for each pixel P , an appropriate combination of N images, obtained by selecting one of the components obtained by said steps (A) to (E), and applying at least one multipass data differential interferometric processing technique, so as to separate the signal associated with the surface deformations from the signals associated with the atmospheric delays.

8. Method according to claim 7, characterized in that it comprises, after said step (G), the following steps:

(H. 1 ) generating for each pixel P an appropriate combination of images, obtained by selecting one of the components obtained by said steps (A) to (E);

(H.2) calibrating the data of said step (A.4) by the low resolution atmosphere and deformation signals estimated in step (G); and

(H.3) estimating the topography and the deformations of said area by applying said steps from (F.1 . 1 ) to (F. 1 . 1 0) , replacing said vectors of said step (F.1 .5) with the calibrated data according to said step (H.2).

9. Method according to anyone of the preceding claims, characterized in that it comprises a plurality of sensors (A i , . .. , AM), and in that said at least one sensor can be constituted by one or more antennas.

10. Method according to anyone of the preceding claims, characterized in that each of said sensors (Ai ..... A_M) can acquire said data simultaneously with the other in each of said orbits O_n, for n - 1 , N, so that (..,, = ... = /„,_,, .

1 1 . Synthetic aperture radar remote sensing system, comprising at least one sensor (Ai , ... , A_M), adapted to emit a radiation towards at least one area and to receive the return radiation,

at least one transport or moving means, such as a satellite or an airplane, on which said at least one sensor is placed, said transport means performing a plurality of passages on said area according to a plurality of trajectories,

data storage means, connected with said at least one sensor (Ai

A_M), to store the data it has collected, and

means for processing the data detected by said at least one sensor (AT , . . . , A_M), connected with said data storage means, characterized in that said processing means run on said data detected by said at least one sensor (A^ ... , AM) the method for processing the data detected by synthetic aperture radar according to the method as defined in claims I 'l l .

1 2. System according to claim 1 1 , characterized in that it comprises a plurality of sensors (Ai , ... , AM), placed on a corresponding plurality of transport means.

1 3. System according to anyone of claims 1 1 or 12, characterized in thai said sensor is an antenna (Ai,..., A^).