Point coherence estimation in interferometry

Abstract

The present invention concerns a computer program product comprising software code portions that are loadable / storable on, and executable by, electronic processing means, and such that to cause, when loaded / stored, the electronic processing means to become configured to carry out a point coherence estimation (PCE) method for points in interferometric images that image a common scene, each at a respective acquisition time. The PCE method includes selecting pairs of points selected in the interferometric images, wherein, for each selected pair of points, a respective equation relates a relative interferometric coherence between the two points of said selected pair with interferometric coherences of said two points, whereby a system of equations is derived from the respective equations of the selected pairs of points and is solved to obtain the respective interferometric coherence of each point.

Description

[0001]“POINT COHERENCE ESTIMATION IN INTERFEROMETRY” CROSS-REFERENCE TO RELATED APPLICATIONS This Patent Application claims priority from European Patent Application No. 23176143.8 filed on May 30, 2023, the entire disclosure of which is incorporated herein by reference. TECHNICAL FIELD OF THE INVENTION The present invention relates, in general, to synthetic aperture radar (SAR) interferometry (also concisely called InSAR) and, in particular, to point coherence estimation in SAR interferometry. However, it can be relevant also to other fields where interferometric images or, more in general, phase signals are involved, such as, for example, laser and sonar interferometry, nuclear magnetic resonance, etc. STATE OF THE ART Nowadays, InSAR is a powerful technology (cf. M. Crosetto et al., “Persistent Scatterer Interferometry: A review”, ISPRS Journal of Photogrammetry and Remote Sensing, vol. 115, pp. 78-89, 2016, and D. Ho Tong Minh et al., “Radar Interferometry: 20 Years of Development in Time Series Techniques and Future Perspectives”, Remote Sens., vol. 12, pp. 1364, 2020 – hereinafter, for the sake of conciseness, denoted as Ref1 and Ref2, respectively) for monitoring from satellite very large areas (even entire continents – cf. M. Costantini et al., “EGMS: Europe-Wide Ground Motion Monitoring based on Full Resolution InSAR Processing of All Sentinel-1 Acquisitions”, IEEE International Geoscience and Remote Sensing Symposium - IGARSS, Kuala Lumpur, Malaysia, pp. 5093-5096, 2022 - hereinafter, for the sake of conciseness, denoted as Ref3) and for detecting motions (typically due to subsidence, landslides, earthquakes, and volcanic phenomena) of imaged scene (e.g., a ground surface, including buildings and infrastructures), with millimetric precision (and sub-metric spatial localization, making it possible to distinguish different parts of landslide bodies, buildings or infrastructures). In fact, thanks to the coherent nature of the SAR emitted radiation, the phase difference of images acquired at different times provides precise information about the possible displacements and deformations of the observed scene occurred between the acquisitions (and also about the three-dimensional (3D) shape of the scene if the acquisitions are taken with slightly different lines of sight). InSAR is a complex technology, and one of its key steps is the identification, among billions of pixels in an InSAR image stack, of the points returning coherent backscattering signals at different acquisition times (these points typically correspond to man-made structures, rocks, or bare soil), whose radiometric and geometric characteristics do not change over time. The concept of persistent scatterer (PS), characterized by the property of backscattering a coherent signal over long series of acquisitions, brought a breakthrough to InSAR (in this respect, reference can be made to A. Ferretti et al., “Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry” IEEE Transactions on Geoscience and Remote Sensing, vol. 38, no. 5, pp. 2202-2212, Sept. 2000, and A. Ferretti et al., “Permanent scatterers in SAR interferometry” IEEE Transactions on Geoscience and Remote Sensing, vol. 39, no.1, pp. 8-20, Jan. 2001 - hereinafter, for the sake of conciseness, denoted as Ref4 and Ref5, respectively), while the alternative approach of working with small baseline (SBAS) interferometric pairs was brought to maturation (cf. P. Berardino et al., “A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms” IEEE Transactions on Geoscience and Remote Sensing, vol. 40, no. 11, pp. 2375– 2383, 2002 - hereinafter, for the sake of conciseness, denoted as Ref6). Since then, different advanced multi- interferogram InSAR techniques have been proposed. including techniques aimed at retrieving weak coherent signals from “smoother” surfaces characterized by distributed scattering properties, sometimes called distributed scatterers (DSs). The different techniques are based on various combinations of statistics of the stack image amplitudes (e.g., amplitude dispersion, signal-to-clutter ratio) and / or phases (e.g., exploiting interferometric coherence) in the spatial and / or temporal domain. In this respect, reference can be made to: Ref7) N. Adam et al., “Wide area persistent scatterer interferometry: Algorithms and examples”, Proc. of Fringe 2011, ESA SP (2011): 1-5; Ref8) K. Goel et al., “Single-network wide-area persistent scatterer interferometry: Algorithms with application to Sentinel-1 inSAR data”, American Geophysical Union (AGU) Fall Meeting, 14-18 Dec. 2015, San Francisco, US; Ref9) R. Lanari et al., “A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms”, IEEE Transactions on Geoscience and Remote Sensing, vol. 42, no. 7, pp. 1377–1386, 2004; Ref10) A. Hooper et al., “A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers”, Geophysical research letters, vol. 31, no. 23, 2004; Ref11) A. Hooper et al., “Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to volcan Alcedo, Galapagos”, Journal of Geophysical Research: Solid Earth, vol. 112, no. B7, 2007; Ref12) B.M. Kampes, “Radar Interferometry – Persistent Scatterer Technique” Springer, 2006, ISBN-10 1-4020-4576-X (HB); Ref13) P. Blanco-Sanchez et al., “The coherent pixels technique (CPT): An advanced DInSAR technique for nonlinear deformation monitoring”, Pure and Applied Geophysics, vol. 165, no. 6, pp. 1167–1193, 2008; Ref14) F. Zhao and J. J. Mallorqui, “A Temporal Phase Coherence Estimation Algorithm and Its Application on DInSAR Pixel Selection”, IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 11, pp. 8350-8361, Nov. 2019; Ref15) M. Costantini et al., “A New Method for Identification and Analysis of Persistent Scatterers in Series of SAR Images”, IEEE International Geoscience and Remote Sensing Symposium - IGARSS, Boston, MA, USA, pp. II- 449-II-452, 2008; Ref16) M. Costantini et al., “Enhanced PSP SAR interferometry for analysis of weak scatterers and high definition monitoring of deformations over structures and natural terrains” IEEE International Geoscience and Remote Sensing Symposium - IGARSS, Melbourne, VIC, Australia, pp. 876-879, 2013; Ref17) M. Costantini et al., “Persistent Scatterer Pair Interferometry: Approach and Application to COSMO-SkyMed SAR Data”, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, manuscript ID JSTARS- 2014-00117; Ref18) G. Zeni et al., “Long-term deformation analysis of historical buildings through the advanced SBAS-DInSAR technique: The case study of the city of Rome, Italy”, Journal of Geophysics and Engineering, 8(3), S1-S12, 2011; Ref19) O. Mora et al., “A. Linear and Nonlinear Terrain Deformation Maps from a Reduced Set of Interferometric SAR Images”, IEEE Transactions on Geoscience and Remote Sensing, vol. 41, pp. 2243–2253, 2003; Ref20) A. Ferretti et al., “A New Algorithm for Processing Interferometric Data-Stacks: SqueeSAR”, IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 9, pp. 3460-3470, September 2011; Ref21) E.A. Hetland et al. “Multiscale InSAR Time Series (MInTS) Analysis of Surface Deformation”, J. Geophys. Res. Solid Earth, vol. 117, pp. 8731, 2012; Ref22) K. Goel, N. Adam, “A Distributed Scatterer Interferometry Approach for Precision Monitoring of Known Surface Deformation Phenomena”, IEEE Transactions on Geoscience and Remote Sensing, vol. 52, pp. 5454–5468, 2014; Ref23) G. Fornaro et al., “CAESAR: An Approach Based on Covariance Matrix Decomposition to Improve Multibaseline– Multitemporal Interferometric SAR Processing” IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 4, pp. 2050-2065, April 2015; Ref24) H. Ansari et al., “Sequential Estimator: Toward Efficient InSAR Time Series Analysis”, IEEE Transactions on Geoscience and Remote Sensing, vol. 55, no. 10, pp. 5637- 5652, October 2017; Ref25) H. Ansari et al., “Efficient Phase Estimation for Interferogram Stacks”, IEEE Transactions on Geoscience and Remote Sensing, vol. 56, pp. 4109–4125, 2018. OBJECT AND SUMMARY OF THE INVENTION Object of the present invention is that of providing an accurate, robust and computationally efficient methodology for estimating coherence of points (where the term point can also refer to a set of adjacent or close pixels, possibly grouped together to form a new pixel, obtained by pre- processing) in interferometric SAR images. This and other objects are achieved by the present invention in that it relates to a computer program product, as defined in the appended claims. In particular, the computer program product according to the present invention comprises software code portions that are: • loadable / storable on, and executable by, electronic processing means; and • such that to cause, when loaded / stored, the electronic processing means to become configured to carry out a point coherence estimation (PCE) method for points in interferometric images that image a common scene, each at a respective acquisition time; wherein: • the PCE method includes selecting pairs of points in the interferometric images; • for each selected pair of points, a respective equation relates a relative interferometric coherence between the two points of said selected pair with interferometric coherences of said two points; and • a system of equations is derived from the respective equations of the selected pairs of points and is solved to obtain the respective interferometric coherence of each point. BRIEF DESCRIPTION OF THE DRAWINGS For a better comprehension of the present invention, preferred embodiments, which are intended purely by way of non-limiting, non-binding examples, will be described hereinafter with reference to the attached drawings (all not to scale) that are related to results of experimental tests performed by the Applicant in order to evaluate performance of the point interferometric coherence estimation method according to the present invention. DESCRIPTION OF EMBODIMENTS OF THE INVENTION General Description The following description is presented to enable a person skilled in the art to comprehend, make and use the invention. Various modifications to the embodiments will be readily apparent to those skilled in the art, without departing from the scope of the present invention as claimed. Thence, the present invention is not intended to be limited to the embodiments shown and described but it is to be accorded the widest scope of protection consistent with the features defined in the appended claims. The present invention is implemented by means of a software program product, loadable / storable in a memory of an electronic processor and executable by the latter, and comprising software code portions for implementing, when the software program product is run on the electronic processor, the point coherence estimation (PCE) method according to the present invention described hereinafter. Throughout this document, for the sake of description simplicity and without losing generality, the term “point(s)” will be used as synonym of the term “pixel(s)” or, sometimes, could refer to set(s) of adjacent or close pixels (possibly grouped together to form a new pixel) obtained by pre-processing such as, for example, multilook or distributed scattering techniques. In addition, the term “arc(s)” will be used as synonym of the term “pair(s) of points”. When the terms “point(s)” and “arc(s)” will be used in the following, the correct meanings thereof will be doubtless understandable on the basis of the respective context in which these terms will be used. Also, the terms image and acquisition, as well the terms set, stack and series of images or acquisitions, will be used as synonyms. In fact, without losing generality, in the description reference is made to stacks of images obtained by temporal series of acquisitions. Therefore, the terms interferometric coherence, temporal interferometric coherence, and sample interferometric coherence (sometimes omitting the word interferometric) will be used as synonyms, to indicate either the intrinsic property of interferometric coherence or its estimated value. The correct meanings will be doubtless understandable based on the respective context. Moreover, the term noise is used to mean disturbances affecting measurements in a manner not characterized by systematicity, e.g. spatially and temporally uncorrelated disturbances. The present invention concerns a novel algorithm, named point coherence estimation (PCE), for identifying coherent points in a stack of interferometric SAR (InSAR) images. The interferometric coherence (related to phase noise) of each point is derived from the interferometric coherences between pairs of points, which can be estimated directly, through an effective and clean procedure, without the need for critical assumptions or approximations or articulated procedures, such as amplitude or phase calibrations, pre-selections of candidate coherent points, spatial or temporal averages and interpolations, hypotheses about the probability of distribution of the interferometric phase, etc. According to the present invention, pairs of neighboring points are considered (typically within a few tens or hundreds of meters) so that the systematic phase contributions slowly variable with position (such as atmospheric propagation delays, orbital artifacts in satellite observations, etc.) are similar in the two points of each considered pair (at least for the majority of the considered pairs), i.e., their differences are negligible w.r.t. a phase cycle. Thanks to this property, the so-called temporal interferometric coherence of each considered pair of points (or arc) is well defined and maximized estimating the components of the phase differences that varies systematically at the different acquisitions, e.g. according to a model describing the elevation and motion differences of the two points as functions of the times and geometries of the stack of acquisitions. As a result, only noise (e.g., temporal, spectral, geometric decorrelation, thermal noise), i.e., the components of the phase differences that are statistically independent spatially and temporally (i.e., in the different selected points and in the different acquisition times), remain and determine the temporal interferometric coherences of the arcs. In the hypothesis that the phase noises in two neighboring points are statistically independent (which might require to exclude pairing nearest neighboring pixels if the images are oversampled), it can be easily demonstrated that the expected value of the temporal interferometric coherence for an arc is equal to the product of the temporal interferometric coherence expected values for the two points forming the arc. So, an overdetermined system of equations involving all the considered pair of points can be written. For completeness, inequalities can be added expressing the property that the absolute value of the interferometric coherence expected values for the two points forming the arc must be less or equal to one. Such overdetermined system of equations (and optionally inequalities) can be solved to obtain the interferometric coherence of each point from the temporal interferometric coherences of the arcs, which can be directly estimated. In particular, applying logarithm to the obtained equations (and optional inequalities), an overdetermined system of linear equations (and optionally inequalities) involving all the considered pair of points can be written. Different efficient solvers exist to solve overdetermined systems of linear equations (and optionally inequalities), which often find the solution corresponding to minimize, according to the L1 or L2 norm, the residuals of the equations. A reliable and consistent estimate of the temporal interferometric coherence of each single point is then obtained, based on which the points with good coherence to be further processed according to InSAR procedures can be selected. It can be noted that the selected coherent points can also be called persistent scatterers (PSs), to be understood in a broad sense as scatterers that keep their interferometric coherence over time, despite the type of scattering mechanism, which can be point-like or distributed. Moreover, it is worth reminding that the PCE method can be applied to full-resolution data as well as to data with resolution degraded for a previous processing such as a multi-look or distributed scattering processing, always paying attention to the fact that, depending on the chosen sampling, the nearest neighboring pixels might be correlated and should not paired when applying the present method. It is important to note that the proposed method does not need any assumption about the probability distribution of the phase noise. However, it is interesting to observe that, when considering a Gaussian probability distribution, each equation of the above overdetermined system expresses the well-known property that the variance of each phase noise difference between a pair of points is the sum of the phase noise variances of the two points if the phase noises in the two points are uncorrelated. Moreover, this well-known property can be used also to formulate an alternative embodiment of the present invention. In fact, for given statistical distributions of the phase noise difference between pairs of points, it is possible to derive (either analytically or numerically) for each arc the phase noise variance from the interferometric coherence, and then solve the overdetermined linear system of equations expressing the mentioned sum of variances property to obtain the phase noise variance (and the related coherence) of each considered point. It is also important to note that, according to the present invention, the information brought from any considered arc is exploited, regardless of whether the arc connects coherent or noisy points. Thanks to this property, the PCE method according to the present invention is quite stable, and applying the algorithm to the whole set or to a subset of the points produces similar results. In fact, it is possible, and it can be convenient in some cases, to iteratively apply the method to previously selected points. In addition, this stability makes it possible also to apply the present technique to a set of preselected coherent point, or PS, candidates, obtained for example by the AD and SCR methods with very relaxed thresholds, or by other quick techniques. However, there is no need for preselection of candidate coherent points or for any iterative selection of coherent points. Such procedures might be applied to further reduce the computational time, which in any case is absolutely affordable. In fact, the present algorithm requires the calculation of the temporal interferometric coherence of a set of arcs, which is common to most InSAR methods, and the solution of an overdetermined system of linear equations (and possibly inequalities), for which very efficient solvers can be used. Mathematical formulation Let us consider a stack of co-registered complex SAR images, and let ^^^^, ^^denote their phases, with ^^ ∈ ^^ denoting the image acquisition times, and ^^ ∈ ^^ the image stack points (or a subset of them). Let us denote with ^^ =(^^, ^^′)∈ ^^ a set of arcs connecting pairs of spatially close points ^^ ∈ ^^ and ^^′∈ ^^, and such that all considered points result connected by arcs and the system of equations defined below in (10) result determined and generally overdetermined. The paired points should be spatially close enough that differences of systematic phase contributions between the two points are negligible compared to a phase cycle, at least for most of the considered pairs; in satellite InSAR applications this means that the distance on the ground between the connected points should be smaller than a few tens or hundreds meters, so that the systematic phase contributions slowly variable with position, such as atmospheric delays and orbital artifacts, are similar in the two points. Let us denote with ^^ =(^^, ^^′)∈ ^^ a set of arcs connecting pairs of acquisition dates ^^ ∈ ^^, ^^′∈ ^^ (in the unique master configuration, which is the typical choice in case of full- resolution data, ^^ includes the pairs between one acquisition date chosen as reference and all the other acquisition dates; however, other configurations can be considered, in particular in the case of data not at full- resolution). Let us define i.e., the double difference, modulo 2 ^^, of the SAR image phases between two acquisition times ^^ =(^^, ^^′)∈ ^^ and two close points ^^ =(^^, ^^′)∈ ^^. Since the components of the phase that are slowly variable in space (such those depending on atmosphere propagation delay or inaccurate orbital knowledge) are almost identical in points that are sufficiently close in the scene (ideally, within a few tens or hundred meters) and so their difference is negligible w.r.t. a phase cycle, it is well known (cf. Ref4) that the phase difference can be modeled as where ^^^^is a function modeling the phase difference signal as a function of ^^ and depends on a few parameters ^^^^determined as ^^^^= arg where, in the simplest case, the weights ^^^^, ^^are all identical and serve as a normalization factor, i.e., with ^^ being the cardinality of ^^ (i.e., the number of independent pairs in the time domain). More in general, the weights ^^^^, ^^can depend on the amplitude values or other parameters of the acquired SAR images, or they can be chosen according to different criteria. In the simplest case, the function ^^^^depends only on the relative mean velocity v^^and the relative height ℎ^^of the targets on the ground corresponding to the arc ^^, i.e., where ^^ is the wavelength, ^^ is the range between the sensor and the target in the master acquisition, ^^ is the incidence angle, ^^^^is the perpendicular-to-the-line-of-sight baseline. However, the function ^^^^can be more complex and, in fact, typically includes other terms describing motion acceleration and higher order motion terms, as well as sinusoidal terms describing periodical motions which could be correlated to temperature or other phenomena. Since the components of the interferometric phase that are slowly variable in space basically cancel out in the difference between close points, whereas the components of the phase differences that varies systematically at the different acquisitions are described by the estimated model (5) or by a more complex one, only the interferogram phase noises (e.g., thermal noise and noises due to temporal, spectral and geometric decorrelation) ^^^^, ^^and ^^^^′, ^^, in the two points contribute significantly to ^^^^, ^^, at least for the majority of arcs We can consider ^^^^, ^^and ^^^^, ^^as samples, at the time pairs ^^ =(^^, ^^′)∈ ^^, of random variables ^^^^and ^^^^defined on each arc ^^ =(^^, ^^′)∈ ^^ and point ^^ ∈ ^^ (or ^^′∈ ^^), respectively. The so-called temporal interferometric coherence (cf. Ref4) associated to each arc ^^ =(^^, ^^′)∈ ^^ can be operationally estimated according to the sample coherence operator with ^^^^and ^^^^, ^^defined by (3) and (4) respectively. The indices|^^^^(^^^^)|or|^^^^(^^^^)|2are typically used in InSAR, but other quantities different from (7) could be used to measure the statistical similarity of the SAR signals in two points, without affecting the validity of the method according to the present invention. However, without loss of generality, reference will be made in the following to the sample interferometric coherence operator defined in (7) and to the interferometric coherence operator defined in (8). The approach according to the present invention is based on a novel technique that makes it possible to estimate the interferometric coherence (related to the noise) associated with each single point of the full-resolution interferometric data stack from temporal interferometric coherences associated to the arcs (7), in a clean and simple way, differently from methods using articulated combinations of coherence values on pairs of points (cf. Ref15, Ref16, Ref17), or needing spatial averaging of neighboring pixels (cf. Ref10, Ref11, Ref14). However, as already noted, the present method can be applied also to data with degraded resolution, for example considering points obtained after a multi-look or distributed scattering processing. Let us establish some further notation. Given a random variable ^^, let ^^(^^)be the expected value of ^^. Moreover, let us define the interferometric coherence operator ^^(^^)for a random variable ^^ as ^^( ^^) = ^^( ^^^^ ^^).(8)The sample temporal coherence γN(^^^^)defined in (7) can be considered an estimator of ^^(^^^^). Moreover, let us define the log coherence operator Γ(^^)as Γ( ^^) = − ln| ^^( ^^)|2.(9)Continuing the discussion of (6), let us assume that for any arc ^^ =(^^, ^^′)∈ ^^ the phase noises ^^^^and ^^^^′are independent variables. The assumption of statistical independence is typically well verified unless the considered images are oversampled, in which case the set of arcs ^^ =(^^, ^^′)∈ ^^ connecting pairs of spatially close points ^^ ∈ ^^ and ^^′∈ ^^ can be defined excluding pairs of very close points for which the statistical independence assumption is not valid. Using the well-known property that the expectation value of the product of independent variables is the product of their expectation values, it follows from (6) and (8) that the interferometric coherences ^^(^^^^)associated with any arc ^^ =(^^, ^^′)∈ ^^ is related to the interferometric coherences from which it follows and, taking the logarithm of both terms of (11), i.e. substituting (9): Note that (10) defines a system of equations for the complex unknown variables ^^ ∈ ^^, whereas (11) is the 2 corresponding systems for the variables | ^^( ^^^^)| , and (12) is the linear system obtained for the variable Γ( ^^^^). The right hand terms of systems (10), (11) or (12) are quantities that can be operationally estimated from (7). An unbiased estimation for the right terms of (11) and (12) is derived below in (15). Also observe that, for how the pairs of points ^^ =(^^, ^^′)∈ ^^ are built, system (10) and the derived systems (11) and (12) are or can be made determined and generally overdetermined (in the worst case redefining the set of considered points ^^ ∈ ^^ as those for which such systems can be made determined or overdetermined). System overdetermination makes for the robustness of the solution. It is worth saying that different weights may be given to the equations of the systems above, which may be useful in some cases, for example to consider that the estimates of the temporal coherences in the right-hand sides can be biased for low coherence values. Moreover, it can be straightforwardly demonstrated from definition (8) that 0 ≤ 2 | ^^( ^^^^)| ≤ 1, and then, from (9), that Γ( ^^^^) ≥ 0. These inequalities can be enforced as constraints to make the solutions of systems (11) and (12) more robust, although the solution obtained without inequality constraints is good enough in most cases. Finally, note that system (12) is a linear system (either unconstrained or with inequality constraints, and generally overdetermined), and can be solved very efficiently. Many efficient solvers find the solution corresponding to minimize, according to the L1 or L2 norm, the residuals of the equations. The solution of any of the systems (10), (11) or (12) provides a reliable and consistent estimate of the phase 2 coherence | ^^( ^^^^)| for each point ^^ ∈ ^^, based on which the points with good interferometric coherence can be selected and used for the next steps of InSAR processing and applications. It is important to note that possible small deviations from the assumptions made above obviously do not invalidate the method and do not critically affect its performance. Moreover, it can be worth noting an interesting property. Assume that (6) is not exactly true (e.g., because the components of the interferometric phase that are slowly variable in space do not completely cancel out in the difference between close points, or the components of the phase differences that varies systematically at the different acquisitions are not perfectly described by the estimated model (5) or by a more complex one), and let ^^^^, ^^be the residual of equation (6), i.e. the deviation from equality in (6). If ^^^^, ^^is statistically independent from the noise terms ^^^^, ^^and ^^^^′, ^^(which is typically true by definition of noise), then ^^^^, ^^will produce an additional multiplicative term in the left-hand side of (10) and (11), and an additional additive term in the left-hand side of (12). For a limited cardinality ^^ of ^^, the temporal sample coherence|γN(^^^^)|2defined according to (7) is a biased estimator of|^^(^^^^)|2, ^^ =(^^, ^^′)∈ ^^. In fact, its expected value can be easily calculated (see demonstration of (22)and (23) belowErrore. L'origine riferimento non è stata trovata.) and results to be: The estimator is however consistent as, from (11)and (13): From equation (13) it results that an unbiased estimation to be used for the right-hand term of system (11) and, after taking the logarithm, system (12), is: It is worth noting that the approach according to the present invention is independent of the probability distribution of the phase noise. However, it is interesting to note that given the probability distribution for the deviations ^^^^, ^^ =(^^, ^^′)∈ ^^, under certain conditions a relation can be determined analytically or estimated numerically between the variance ^^2(^^^^)and the interferometric coherence γ(^^^^). For example, if the variables ^^^^, ^^ =(^^, ^^′)∈ ^^ have a zero-mean Gaussian probability distribution, from the definition (8) of the coherence operator ^^(^^)it can be easily derived (or directly seen using the known characteristic function of a Gaussian distribution) that Then, reminding that for any arc ^^ =(^^, ^^′)∈ ^^ the phase noises ^^^^and ^^^^′are assumed statistically independent, and the well-known property that the variance of the sum of uncorrelated variables is the sum of the variances of the variables, the following equations can be written for ^^ =(^^, ^^′)∈ ^^: which is a linear system of equations in general analogous to (12), and exactly identical to (12) when the variables ^^^^have zero-mean Gaussian probability distribution. The mentioned procedure thus constitutes an alternative way to obtain a system of equations that makes it possible to determine the variance and the related interferometric coherence of points from the interferometric coherence of arcs, which can be directly estimated. From the foregoing it is evident that the present invention provides a novel algorithm for point coherence estimation (PCE), and then for selection of coherent points, or PSs, in SAR interferometry. Experimental results In all the experimental tests performed, the method has proved to be very effective, providing for each single full- resolution point a reliable measure of the temporal interferometric coherence (which is under certain conditions related to the phase noise variance), therefore making it possible to detect a very large number of coherent points, or PSs, with very few false detections. Hereinafter, some tests performed on two stacks of Sentinel-1 interferometric SAR images acquired over a pre- alpine area in Piemonte, Italy (86 acquisitions from January 2020 till October 2022), and over an area between Sicily and Calabria, Italy, including the Etna volcano (120 acquisitions from January 2020 till December 2022), respectively, will be described. The analyzed areas are affected by different kinds of displacement phenomena associated to natural and anthropic activities, and include different types of land cover, among which continuous and discontinuous urban fabric, transport infrastructures, bare soil, agricultural fields, mountains, and a big volcano. The quality of the obtained results can be clearly appreciated by visual inspection of the selected coherent points, or PSs, in very high-resolution optical images of the ground. In this respect, reference can be made to Figures 1 and 2, wherein: • Figure 1 includes images related to interferometric analysis of a stack of 86 Sentinel-1 SAR images (track 015) acquired from January 2020 till October 2022 in a pre-alpine area, Piemonte, Italy; in the first two image rows, the PSs identified with the present method are compared with the PSs selected by the joint use (logical OR) of the signal-to- clutter ratio (SCR) and the amplitude dispersion (AD) methods, in three areas with different types of ground cover, recognizable from the background very high resolution optical image; the present algorithm largely outperforms the classical, although basic, SCR and AD methods, both in density and coverage of the identified PSs; the thresholds for PS identification were chosen in order to have about the same false detection rates for the three methods; in the third image row, the corresponding mean velocities obtained for the PSs identified with the present algorithm are shown; • Figure 2 shows mean velocities obtained for the coherent points, or PSs, identified with the present algorithm on a stack of 120 Sentinel-1 interferometric SAR images (track 044) acquired from January 2020 till December 2022 in Sicily and Calabria, Italy; from left to right, the overall view of the processed area, the Etna volcano, characterized by strong motion patterns that follow fault lines, and the famous town of Taormina; the present method is able to identify PSs wherever interferometric coherence is expected, including bare soil and areas characterized by strong and complex displacement phenomena. More in detail, some comparisons of the method according to the present invention are shown in Figure 1 with respect to the PSs identified by the classical, although basic, methods of the amplitude dispersion (AD) and / or signal-to- clutter ratio (SCR). Choosing thresholds such that the three methods approximately have the same PS false detection rate, the present method provides significant improvements not only to the PS density, but also to the PS coverage of the ground, i.e., more areas and objects are covered by PS measurements. The selected coherence threshold in the reported examples was 0.66, which corresponds to a phase noise standard deviation of 0.91 radians if a Gaussian probability distribution is assumed for the phase noise. In order to better clarify the difference between the three considered PS identification methods, two-dimensional (2D) histograms relating AD, SCR and the temporal interferometric coherence estimated by the present method were computed. In this connection, reference can be made to Figure 3 which shows histograms comparing the three different persistent scatterer selection methods. The first two images show the 2D histograms, based on all processed pixels relative to the stack of Figure 1, of the temporal interferometric coherence (estimated with the present algorithm) and the AD (left) or the SCR (center). The different gray levels indicate histogram frequency. The dashed lines indicate the thresholds used with the different methods for the PS selection shown in Figure 1. The right image shows the 2D histogram, based on the points identified as PSs by the present algorithm, of the SCR and AD. The rectangle in the third image delimits the portion of these PSs that would be not detected by considering AD or SCR. The analysis of the three histograms in Figure 3 shows that the algorithm according to the present invention is able to identify also PSs characterized by low SCR and high AD, confirming the effectiveness of the method. FURTHER FORMULATIONS Factorization of the arc coherence The residual interferometric phase ^^^^at any point ^^, after the removal of the signal model, can be assumed a zero- mean random variable with a probability density function . The interferometric coherence associated to the residual interferometric phase ^^^^of a point ^^ ∈ ^^ is defined as the expectation of the phasor ^^^^ ^^ ^^: Let ^^^^ ^^( ^^^^) be the probability density function associated to the interferometric phase residual difference ^^^^= ^^^^− ^^^^′on the arc between two points ^^ = ( ^^, ^^′) ∈ ^^, after removal of the signal model. The interferometric coherence associated to the arc ^^ = ( ^^, ^^′) ∈ ^^ is defined as: (19) Given an arc ^^ =(^^, ^^′)∈ ^^, the random variables ^^^^and ^^^^′can be assumed mutually independent and the probability function ^^^^ ^^(^^^^)can be considered as the following convolution between the probability functions of the interferometric phase residuals in the points ^^ and ^^′: where ∗ denotes the convolution operation. Inserting such convolution equation in the (19) and performing the change of variable ^^^^′ = ^^^^− ^^^^, the factorization property can be derived as follows: This result can be considered a particular case of the well-known factorization of characteristic function of sum of independent random variables. Moreover, it can be directly derived (very simply, but maybe with some lack of notation rigor) using the well-known property that the expected value of the product of independent random variables is equal to the products of the expected values: Expectation of the arc coherence estimator 1 Formula (7), when ^^^^, ^^= ^^, with ^^ the cardinality of ^^, is an estimator of the arc interferometric coherence ^^( ^^^^) and ^^ = ( ^^, ^^′) ∈ ^^. In the following, the expectation of this estimator is calculated as a function of the number of the interferograms used in the estimation, ^^, and of the arc coherence ^^( ^^^^) or of the corresponding point coherences and ^^( ^^^^′). If independence of the residual phases ^^^^and ^^^^′in the erent points ^^ and ^^′is considered: ^ ^ If independence of the residual phases in the point ^^ at erent times ^^ and ^^′is assumed: ^ ^ Using the factorization property (11) of the arc coherence, the following equation can be also written: Pointwise temporal coherence estimator For each point ^^ ∈ ^^, let define the pointwise temporal coherence estimator, calculated on ^^ interferograms, as: ^ (24) where ^^^^, ^^are the residual interferometric phases of the point ^^ at ^^ = ( ^^, ^^′) ∈ ^^. Bias of the arc coherence estimator The following equation characterizes the bias in function of N and of the coherences in the points ^^ and ^^′. ^^ ^^ ^^ ^^ As expected, the arc coherence estimator is biased and the bias decreases for large values of N and for high values of the coherences in the points ^^ and ^^′. Factorization of the arc coherence estimator In general, the arc coherence estimator can be expressed through the pointwise coherence estimators as in the following equation: Further description of an Embodiment of the Invention In the following an embodiment of the invention, corresponding to the case mentioned in the discussion of (17) above, is further discussed. Given the probability distribution of the phase noise, under certain conditions the relation between the variance of the phase noise and the corresponding temporal coherence can be either calculated analytically or retrieved by simulation. A lot of works have been published about phase statistics, particularly for the case of distributed scattering. In this respect reference can be made to: Ref26) J.-S. Lee et al., “Intensity and phase statistics of multilook polarimetric and interferometric SAR imagery”, IEEE Transactions on Geoscience and Remote Sensing, vol. 32, no. 5, pp. 1017–1028, 1994; Ref27) R. Touzi and A. Lopes, “Statistics of the stokes parameters and of the complex coherence parameters in one- look and multilook speckle fields” IEEE Transactions on Geoscience and Remote Sensing, vol. 34, no. 2, pp. 519–531, 1996; Ref28) R. Bamler and P. Hartl, “Synthetic aperture radar interferometry”, Inverse problems, vol. 14, no. 4, p. R1, 1998; Ref29) P. A. Rosen et al., “Synthetic aperture radar interferometry” Proc. IEEE, vol. 88, no. 3, pp. 333-382, March 2000; Ref30) R. F. Hanssen, “Radar interferometry: data interpretation and error analysis”, Springer Science & Business Media, 2001, vol. 2; Ref31) A. Pepe, “Theory and statistical description of the enhanced multi-temporal InSAR (E-MTInSAR) noise- filtering algorithm”, Remote Sensing, 11(3), 363, 2019. However, specific assumptions about the scattering mechanisms are not needed, but plausible distributions for the phase noise can be directly assumed (e.g., Gaussian, uniform over a small interval, etc.). As a matter of fact, the choice of the phase noise probability distribution is not critical, since the relation between phase variance and temporal coherence differs significantly only for low coherence values, and the method is quite insensitive to it. From the phase noise variances of the pairs of neighboring points, in the hypothesis that there is no statistical correlation between the phase noises of the paired points, the phase noise variance of each point can be obtained by solving a minimization problem involving all the considered pairs of points. Lastly, it is also possible to convert back the estimated phase noise variance of each single point to the corresponding temporal coherence value, using the inverse method to that considered for passing from coherence to variance, which further reduces the impact of the phase noise probability distribution chosen to derive such relations. Finally, the PSs can be identified among the points of the full-resolution interferometric data stack based on the retrieved phase noise variance, or the corresponding temporal coherence. ***** In conclusion, it is clear that numerous modifications and variants can be made to the present invention, all falling within the scope of the invention, as defined in the appended claims. Moreover, it is important to note that the present invention can be advantageously exploited also with interferometric images different than SAR images or, more in general, with other data sets involving phase signals.

Claims

CLAIMS 1. Computer program product comprising software code portions that are: • loadable / storable on, and executable by, electronic processing means; and • such that to cause, when loaded / stored, the electronic processing means to become configured to carry out a point coherence estimation method for points in interferometric images that image a common scene, each at a respective acquisition time; wherein: • the point coherence estimation method includes selecting pairs of points in the interferometric images; • for each selected pair of points, a respective equation relates a relative interferometric coherence between the two points of said selected pair with interferometric coherences of said two points; and • a system of equations is derived from the respective equations of the selected pairs of points and is solved to obtain the respective interferometric coherence of each point.

2. The computer program product of claim 1, wherein the point coherence estimation method includes selecting the points in the interferometric images such that the system of equations is at least determined or, preferably, overdetermined.

3. The computer program product according to claim 1 or 2, wherein most of the pairs of points selected in the interferometric images are formed by two points that are: • spatially close enough that differences of systematic phase contributions between the two points are negligible compared to a phase cycle, and • spatially far enough for phase noise to be statistically independent at the two points.

4. The computer program product according to any claim 1-3, wherein, for each selected pair of points, therespective equation relates the relative interferometric coherence between the two points of said selected pair with the product of the two interferometric coherences of the two points of said selected pair; wherein the relative interferometric coherence between the two points of each selected pair is computed by a given mathematical procedure, while the single interferometric coherences of the points are computed by solving the system of equations.

5. The computer program product of claim 4, wherein the given mathematical procedure includes: • a known procedure for calculating the relative interferometric coherence between a selected pair of points in the interferometric images, or • said known procedure and further applying a linear function with coefficients depending on number of samples to the relative interferometric coherence obtained with said known procedure thereby obtaining an unbiased coherence estimate.

6. The computer program product according to claim 4 or 5, wherein, for each selected pair of points, a logarithmic function is applied to terms of the respective equation thereby obtaining a respective new linear equation that relates the logarithm of the relative interferometric coherence between the two points of said selected pair with the sum of the logarithms of the two interferometric coherences of the two points of said selected pair, whereby the derived system of equations is a system of linear equations.

7. The computer program product according to any preceding claim, wherein the point coherence estimation method comprises including in the derived system of equations to be solved also one or more predefined inequality constraints related to a priori knowledge on admissible values for unknown variables of the system.

8. The computer program product according to claim 1 or 2, wherein, for each selected pair of points, the respectiveequation relates the sum of variances of phase noises of the two points of said selected pair with a given variance of phase noise difference between the two considered points; wherein said given variance is estimated analytically or numerically based on a given probability distribution from the relative interferometric coherence between the two points of the selected pair, which in turn is computed by a given mathematical procedure; whereby the derived system of equations is a system of linear equations.

9. The computer program product of claim 8, wherein most of the pairs of points selected in the interferometric images are formed by two points that are: • spatially close enough that differences of systematic phase contributions between the two points are negligible compared to a phase cycle, and • spatially far enough for phase noise to be statistically independent or uncorrelated at the two points.

10. The computer program product according to claim 8 or 9, wherein the given mathematical procedure includes: • a known procedure for calculating the relative interferometric coherence between a pair of points in the interferometric images, or • said known procedure and further applying a linear function with coefficients depending on number of samples to the relative interferometric coherence obtained with said known procedure thereby obtaining an unbiased coherence estimate.

11. The computer program product according to any claim 8-10, wherein the given probability distribution is a zero- mean Gaussian probability distribution.

12. The computer program product according to any preceding claim, wherein the point coherence estimation method further includes select, based on estimated interferometric coherences and / or phase noise variances, points to be further processed in subsequent steps ofinterferometric processing and / or interferometry applications.

13. Processor storing the computer program product as claimed in any preceding claim.

Images