SAR Tomography as an Add-On to PSI : Detection of Coherent Scatterers in the Presence of Phase Instabilities

The estimation of deformation parameters using persistent scatterer interferometry (PSI) is limited to single dominant coherent scatterers. As such, it rejects layovers wherein multiple scatterers are interfering in the same range-azimuth resolution cell. Differential synthetic aperture radar (SAR) tomography can improve deformation sampling as it has the ability to resolve layovers by separating the interfering scatterers. In this way, both PSI and tomography inevitably require a means to detect coherent scatterers, i.e., to perform hypothesis testing to decide whether a given candidate scatterer is coherent. This paper reports the application of a detection strategy in the context of “tomography as an add-on to PSI”. As the performance of a detector is typically linked to the statistical description of the underlying mathematical model, we investigate how the statistics of the phase instabilities in the PSI analysis are carried forward to the subsequent tomographic analysis. While phase instabilities in PSI are generally modeled as an additive noise term in the interferometric phase model, their impact in SAR tomography manifests as a multiplicative disturbance. The detection strategy proposed in this paper allows extending the same quality considerations as used in the prior PSI processing (in terms of the dispersion of the residual phase) to the subsequent tomographic analysis. In particular, the hypothesis testing for the detection of coherent scatterers is implemented such that the expected probability of false alarm is consistent between PSI and tomography. The investigation is supported with empirical analyses on an interferometric data stack comprising 50 TerraSAR-X acquisitions in stripmap mode, over the city of Barcelona, Spain, from 2007–2012.


Introduction
Persistent scatterer interferometry (PSI) [1][2][3][4][5][6][7] is nowadays an operational geodetic technique for the monitoring of surface deformation with spaceborne synthetic aperture radar (SAR) data stacks.These stacks typically comprise several repeat-pass SAR acquisitions, spanning from months to years.PSI techniques attempt to extract the interferometric phase components correlated with the scatterer motion.The quality of the deformation estimates is tied to the precision of the interferometric phases.Temporal and geometric decorrelation, as well as uncompensated platform motion and atmosphere-induced optical path delay variations, are among the factors that cause random instabilities in phase.For these reasons, a quality control is necessary during the processing as well as when reporting the final results.
The single dominant scatterers that exhibit long-term phase stability are generally termed as persistent scatterers (PS).PSI processing approaches often use a classifier to identify a priori a set of PS candidates, e.g., the permanent scatterers [1] approach uses the dispersion index as a proxy for phase stability.The PSI approaches based on the interferometric point target analysis (IPTA) framework, as in [3,8], employ low spectral diversity [3, [9][10][11] as a proxy for phase stability in addition to the stability of the backscattering amplitude.Low dispersion index and low spectral diversity are indicative of good phase quality.The observed differential interferometric phases are fit to a phase model and the unknown parameters, such as the deformation velocity and the residual topography, are thereby estimated.The dispersion of the residue of the fit is a means to characterize the quality of the estimates.It is often used to compute the multi-interferogram complex coherence (MICC) [1,12,13] which can in turn be used as a test statistic to perform statistical detection i.e., to decide among the hypotheses whether a given PS candidate is a phase coherent single scatterer or if it comprises noise only.The statistics of the noise impact the probability of false alarm in the detection process.
An inherent limitation associated with PSI techniques is the fact that a phase-only model cannot consider multiple coherent scatterers with different complex reflectivity interfering in the same range-azimuth resolution cell.The cumulative phase response in this case is mismatched to the interferometric phase model, which is essentially based on the assumption of a single scatterer.Consequently, it may lead to erroneous estimation of the deformation parameters.Therefore, PSI processing approaches typically reject the cells that contain backscattering contributions from multiple scatterers, as for the case of layovers.
The aforementioned limitation can be alleviated by SAR tomography [14][15][16][17], which exploits both the amplitude and the phase of the received signal, thereby permitting a higher order analysis [18].It allows 3-D reconstruction of the scene reflectivity-a feature that renders it possible to resolve the layover problem [19][20][21][22].Additionally, differential SAR tomographic methods [23][24][25] allow a joint spatio-temporal inversion of the coherent scatterers in layover, i.e., the position along the elevation axis as well as the deformation velocity of the interfering scatterers are simultaneously estimated.Therefore, differential SAR tomography has been proposed as an add-on to PSI techniques to improve deformation sampling by resolving the scatterers in layover that are rejected in the PSI processing [26][27][28][29].Inevitably, a detection strategy is again required to classify whether the detection of one or more scatterers in the same resolution cell is true or false.In this context, it is pertinent to carry forward the same quality criteria as used in the prior PSI analysis so that the combined use of PSI and tomography holds compatibility.
The prevailing detection mechanisms for SAR tomography, such as the generalized likelihood ratio tests in [13, 24,30], consider an additive noise model for the received complex signal vector.The source of the noise is attributed to the clutter in the resolution cell.However, the instabilities in the observed interferometric phases, albeit considered additive in the phase-only model, naturally represent themselves as multiplicative noise in the tomographic signal model.Therefore, in order to carry the impact of the phase instabilities from an interferometric to tomographic analysis, the detection strategy employed for hypothesis testing needs to account for the phase instabilities as a multiplicative disturbance in tomographic inversion.
Keeping in view the aforementioned concerns, this paper describes a strategy for the detection of single and double scatterers with SAR tomography whereby the hypothesis testing is directly linked to the MICC-based test statistic for PS detection in the prior PSI processing.As a whole, this paper is a follow-up to the earlier works in [12,27,31].Section 2 presents the mathematical models typically used for SAR interferometry and tomography, as well as the associated detection mechanisms.Section 3 presents the processing methodology adopted in the paper.The data stack for empirical analysis is introduced in Section 4. The results obtained are presented in Section 5, followed by a discussion in Section 6.

Models
We consider the availability of a coregistered, single-reference interferometric SAR data stack comprising M layers of repeat-pass interferograms.For a given range-azimuth resolution cell in an interferometric layer, we denote the received single-look complex (SLC) signal as y m = z m exp (−jϕ m ), where z m = |y m | ∈ R is the amplitude of the received signal, and ϕ m is the observed interferometric phase.The subscript m, where m ∈ {0, 1, . . . ,M − 1}, is used to indicate a specific layer in the interferometric stack.In the following text, an underlined symbol represents a quantity that has been modeled as stochastic, or when the distinction between observables versus observations is emphasized.Bold symbols represent vectors, or matrices when capitalized.

Interferometric Phase Model
The interferometric phase observable, ϕ m , is generally modeled as a sum of several phase contributions [32,33]: where ϕ disp is the phase change due to the linear displacement of the target as a function of time within the resolution cell: λ is the wavelength, v is the deformation velocity in the line of sight (LOS), and t m is the temporal baseline for the m th interferogram.ϕ geo is the phase variation due to sensor-to-target geometry.Neglecting higher order terms [16,34], where b ⊥ m and b m are the orthogonal and parallel components of the spatial baseline for the mth interferogram, respectively.ρ 0 is the range distance from the sensor to the target location for the reference acquisition.s represents the elevation, i.e., the position of the target in the axis perpendicular to the LOS.In case of thermal expansion, the additional phase variations are linearly modeled as follows [27,35]: where T m is the temperature change (with respect to the temperature for the reference layer), and η is the phase-to-temperature sensitivity.The term 2π p, where p ∈ Z, is added to account for phase wrapping.The phase variations ϕ atm m are due to the optical path length variations while propagation through the atmosphere.They are modeled as stochastic variables due to the temporally varying nature of atmospheric refractivity [36][37][38][39].The phase decorrelation term, ϕ decor m is, by definition, a random quantity, which is typically modeled as an additive phase noise.The parameters s, v and η are treated as deterministic unknowns in this work.
The interferometric phase model in Equation ( 1) is implicitly assuming the presence of a single coherent scatterer in the resolution cell.In case of multiple coherent scatterers in the same resolution cell, it is not possible to write the interferometric phase, ϕ m as a sum of the aforementioned sources of phase variations, independently of the reflectivity of the individual scatterers.

PSI: Model of Observation Equations
While several approaches to parameter estimation with PSI have been proposed over time, as in [1][2][3][4][5][6], the functional model of interferometric phase observation equations common to these approaches is as follows [33]: where ϕ is the M × 1 vector of interferometric phase observables, A is the design matrix, and p is the vector of the aforementioned unknown parameters.w is the M × 1 vector of phase residuals which collectively represent the phase instabilities owing to decorrelation, uncompensated atmospheric phases and model imperfections.The residuals in each layer are assumed to be zero-mean and independent random variables: E {w} = 0; and D {w} = E ww H = Q ww is the covariance matrix for the residuals.If it can be assumed there are no phase unwrapping issues, and the data stack can be phase calibrated by compensating for the atmospheric phase with external data-although both assumptions are simplistic-then the remaining unknowns are s, v and η.The design matrix is then constituted by the coefficients of these parameters (from Equations (2-4)) [1,33].Under Gauss-Markov conditions, the best linear unbiased estimate of the parameter vector using weighted least squares is given as [33]: The covariance matrix of the estimated parameter vector, Q p p = D p is as follows: The quality of the estimates is, therefore, dependent on the dispersion of the residuals.The vector of the estimated phase residuals is as shown below:

PSI: Statistics for PS Detection
For each PS candidate, we distinguish between the following two hypotheses: H 0 -the null hypothesis.The range-azimuth resolution cell does not contain any coherent scatterer and comprises merely clutter; H 1 -the alternative hypothesis.The cell contains a phase coherent single scatterer, i.e., a PS.
In the presence of a coherent scatterer whose phase response is well-matched to the model in Equation (1), the phase residuals are expected to have a low dispersion around the expected value of zero.Contrarily, in the absence of a coherent scatterer, the observed phase and the residuals are expected to have a wider dispersion.With these considerations, we assume that the phase residuals generally follow a von Mises (circular normal) distribution.The probability density function (PDF) is given by [40]: where the support of the distribution is any 2π interval.The parameter µ = E {w} represents the 'preferred direction', which we consider to be zero under both H 0 and H 1 .The support is then the interval [−π, π) and the distribution is symmetric about zero.The parameter κ ≥ 0 is a measure of 'concentration' of the distribution around the mean value, i.e., κ −1 behaves analogously to the dispersion of a linear random variable.I o (κ) is the modified Bessel function of the first kind and order zero.Under H 1 , we consider the residuals to exhibit a higher concentration around µ.
NB: The term circular distribution as used in this paper refers to a directional distribution with support on the circumference of unit circle [40].

Test Statistic
A commonly used statistic to test among the two hypotheses is the ensemble coherence, as defined below [5,32]: An unbiased estimator of the coherence, given M interferometric layers, is the multi-interferogram complex coherence (MICC) [12,32]: where X and Y are the sum of cosine and sine terms in the expression, respectively, and the length of the resultant, R = The overscore indicates sample mean.Hereafter, we refer to MICC simply as the sample coherence.The phase residuals ŵm are assumed to be independent and identically distributed (i.i.d.) random variables.
In the context of interferometry, we typically use the coherence values normalized between 0 and 1, i.e., | γ|, instead of the resultant R = M| γ|.However, in the directional statistics literature, the use of the term R is more common.Here, we state both to facilitate cross-referencing with the literature.The sample mean direction μ, computed with sample coherence for any random sample (w 1 , w 2 , . . ., w m ) from a von Mises population, is the maximum likelihood estimator of the preferred direction µ when R is well-defined [41,42].This property is characteristic of von Mises populations on a circle, analogous to a similar property holding for Gaussian distribution on a real line whose location parameter is estimated with maximum likelihood by the sample mean [40,43].

Statistics under H 0
The statistics of the sample coherence depend on the distribution of the phase residuals.With reference to Equation (11), the phase residuals can be considered as angles subtended by phasors of unit length.Under H 0 , when the phasors have no preferred direction, we consider the limiting case of von Mises distribution when κ → 0 [40]: where U (w) is the circular uniform distribution.In this case, E {x} = E {y} = 0; therefore, E γ = 0.
The second order moments are E x 2 ; H 0 = E y 2 ; H 0 = 1 2 .The terms x and y are not independent (as x 2 + y 2 ≡ 1), but they are uncorrelated as E {x • y} = 0 [12,40].The variance of the addends in the Equation ( 12) is finite.Therefore, under the assumption of a large sample size, multivariate central limit theorem holds, and we consider the joint distribution of ( X, Ȳ) to be converging to a Gaussian distribution, N 2 0, Σ γ where R is then approximately Rayleigh-distributed, and its PDF is as follows [40]: where 0 ≤ r ≤ M. Referring to [12,40], the probability of false alarm can be computed as the upper tail of the Rayleigh distribution, as follows: It can be equivalently expressed as where T γ is the detection threshold such that 0 ≤ T γ ≤ 1.

Statistics under H 1
In case of H 1 , the probability distribution of R is given by [40], J 0 is the Bessel function of the first kind and zero-order.A closed form expression for the PDF is not available.We again assume a large sample size and invoke the multivariate central limit theorem.It allows us to consider the joint distribution of X and Ȳ to be asymptotically normal, and expectation and the variance of the sample coherence can be approximated as follows [44]: var where ν j = E {cos (jw)}: For sufficiently large κ, the von Mises distribution for the phase residuals can be approximated by a linear normal distribution with σ 2 w = κ −1 [40].The coherence in this case is given by [12,31]: For a discussion on the details about the corresponding probability of detection, interested readers are referred to earlier works in the literature [12,13].
Since exact closed-form expressions for the PDF of | γ| are not available, we resort to numerical methods to compare the estimate of the coherence magnitude for the general case of κ > 0 against the estimate in case of the aforementioned linear normal approximation.For selected values of κ between [1, 10], we perform 10 5 Monte Carlo simulations of the residual phase vector, w (comprising M instances of von Mises distributed random variables), and compute the coherence magnitude.
The results are shown in Figure 1 for three different values of M. The estimate under the normal approximation (Equation ( 24)) is also shown.It can be seen that the normal approximation for the limiting case tends to overestimate the coherence magnitude.The overestimation decreases for increasing values of κ.For κ > 3, the difference between the coherence estimate under the assumption of von Mises distribution and the normal approximation is less than 5% on average.With increasing number of acquisitions, the variance in the estimation of the coherence magnitude decreases (in agreement with Equation ( 21)).

SAR Tomography: Mathematical Model
In the absence of noise, for a given range-azimuth resolution cell, the mathematical model for SAR tomography (3-D SAR) can be written as [16,19,21,26,45]: where α is the complex reflectivity and I s is the support of s.This model assumes there has been no displacement in the line of sight during the observation time period.Differential SAR tomography [23,25] with extended phases models [25,27,46] allows modeling linear displacement as well as seasonal or temperature-induced motion: where ψ m is the sum of the deterministically modeled phase components as a function of the unknown parameters, i.e., It is assumed that the phase terms (and hence the spatial and temporal baselines, and temperature changes) are mutually independent of each other.A general mathematical model for SAR tomography can be defined as follows [27,47]: where P represents the support of the parameter vector (i.e., the parameter space), and p ∈ P. It is analogous to a multi-dimensional Fourier transform [48].In case the resolution cell contains a single point source with dirac delta response, α(p) = τ 1 δ(p − p 1 ), with τ 1 ∈ C, Equation (28) reduces to the following: For the general case of Q point sources in the presence of clutter, the tomographic model is further extended as follows: where n m represents additive noise which is typically modeled as zero-mean complex Gaussian (with symmetric variances for the real and imaginary parts).We assume the noise samples are i.i.d.across the stack, i.e., D {n} = σ 2 n I M , with σ 2 n > 0. d m represents the coherent sum of the deterministic components in the signal vector.τ q is the reflectivity, and ψ m (p q ) is the modeled phase for the qth scatterer.

SAR Tomography: Model Inversion and Parameter Estimation
We use single-look beamforming for the inversion of the general tomographic model to estimate the unknown scatterer reflectivity as a function of the parameter vector p for a given range-azimuth resolution cell as follows [13,16]: where ., .represents the inner product, a (p) is the steering vector as a function of p, and y is the vector comprising the SLC observations: The steering vector is structured as follows: For the estimation of the unknown parameters, we use the estimated absolute reflectivity as the objective function in the following maximization: As more than one coherent scatterer may be present in the same resolution cell, successive maxima after the global maximum may indicate the presence of more scatterers.Assuming a maximum of two scatterers, an estimate of the parameter vector for the second scatterer is obtained as follows: where δp indicates the Rayleigh resolution for the tomographic profile along each of the unknown parameters.Equation (32) implies that noise in the SLC vector will cause errors in the reconstructed target reflectivity.As a consequence, errors will propagate in the estimation of the parameters using the aforementioned maximizations.Therefore, a scatterer detection strategy is needed to classify whether a given resolution cell contains one or more phase coherent scatterers, or is merely clutter.

SAR Tomography: Statistics for Scatterer Detection
A commonly used test statistic for coherent scatterer detection in the context of tomography is the absolute value of the estimated reflectivity, | α|.The same hypotheses are carried forward as introduced in Section 2.3, except for the change that now we consider them for multiple coherent scatterer candidates for each pixel.We consider a maximum of two candidates per pixel.In case only one of the candidates fulfills H 1 , we call the pixel a single scatterer.In case both the candidates fulfill H 1 , the pixel is called a double scatterer.

Statistics under H 0
In case the received signal is merely clutter, the received signal vector y = n.Using Equation ( 32), where the third equality follows from rotational invariance of the Gaussian distributed samples, and, therefore, the inconsequential difference between ǹ and n will be dropped.Since ϕ m = ∠n m under H 0 , the observed interferometric phase (and the residual phase in this case) follows a uniform distribution [12].Along similar lines as in Section 2.3, the joint distribution of the real and imaginary parts of α is a zero-mean Gaussian with the following covariance matrix: The PDF of | α| in this case is Rayleigh, and the right tail probability to compute the probability of false alarm is as follows: where y 2 2 = ∑ m z 2 m is the squared L2-norm of the observed signal vector.

Statistics under H 1
In general, the received signal contains clutter besides the possibility of backscattering contribution from point-like sources.We assume that, under H 1 , the deterministic backscatter from the point sources is dominant over the clutter, i.e., |d m | |n m | ∀ m.This assumption allows us to consider that the observed phase owes primarily to the vector sum of the backscatter from point-like sources (and not the clutter).Using Equations ( 30) and (32), the expression for the estimated reflectivity can then be stated as follows: Formally, the origin of phase instability, ŵm in Equation ( 43) is not the clutter, rather it is phase disturbances such as uncompensated atmospheric phase delay variations or residual motion [31], or phase model imperfections.Using Equations ( 10) and ( 39), From Equation (44), it is clear that phase instability is disturbing tomographic reconstruction in a multiplicative sense.The ensemble coherence has a direct impact on the expected value of the retrieved reflectivity profile, and thereby on the hypothesis testing.Closed-form expression s for the PDF of | α| are not available when the residuals are assumed to follow a von Mises distribution with κ > 0. A Rician approximation can be taken, as suggested in [31], when the residuals can be considered to be normally distributed (i.e., the limiting case when κ → ∞).The probability of detection, f D for a fixed false alarm rate can then be studied as the area under the upper tail of the Rician distribution [49].
Nonetheless, we resort to Monte Carlo simulation to study the probability of detection numerically in terms of the inverse coefficient of variation (iCV) as defined below for the text statistic α: This definition has been referred to as the signal-to-noise ratio (SNR) in [31].Although, in the field of signal processing iCV is often referred to as the SNR, we avoid referring it so.In our context, formally the denominator in Equation ( 46) is not representing the noise power, neither additive (σ 2 n ) nor multiplicative (σ 2 w ), but rather the dispersion of the test statistic.Considering n ≈ 0, and dropping the dependence on p to simplify notation, Using the assumption that the residual phases are i.i.d.random variables, the covariance term in Equation ( 47) simplifies as follows: cov e j ŵl , e −j ŵk = (1 where δ [ ] is the unit sample function.Using this result, Equation (47) reduces to the following [31]: Therefore, Since y 2 ≤ y 1 ≤ √ M y 2 [50], we reach the following bounds on the iCV for a given level of coherence: iCV is a function of the ensemble coherence as well as on the ratio of the L1 to L2 norm of the signal vector.While the coherence is in turn a function of the concentration of the phase residuals (as shown in Figure 1), the L1-L2 ratio is influenced by the (1) number of acquisitions and (2) the number of point-like scatterers in the same resolution cell.Figure 2 shows the variation of the empirically estimated iCV against the concentration parameter for different numbers of scatterers, for M = 50 acquisitions as an example.In addition, 10 5 realizations of the phase residue are generated under a von Mises distribution for each value of κ selected between (0, 20].The dashed lines in Figure 2 highlight the upper and lower bounds on iCV.The upper bound is reached theoretically when Therefore, the greater the number of acquisitions, the higher is the achievable iCV.At a given concentration of phase residuals, the iCV decreases for an increasing number of scatterers.The iCV estimates for Q = 1 converge at the upper bound.The impact of the number of scatterers on the iCV is further discussed in Appendix A. Concentration parameter, κ iCV (dB) Monte Carlo runs with M = 50 acquisitions iCV against concentration parameter Figure 3a is a plot of the numerically estimated f D against the iCV.The detection thresholds are set to ensure a fixed level of probability of false alarm, f F ∈ 10 −2 , 10 −3 , 10 −4 given M = 50 acquisitions.Lower levels of f F provide higher f D , indicating the trade-off typically observed for statistical detectors [49].At the same time, we observe a slight dependency of f D on the number of scatterers.Even for a fixed level of iCV, the f D is lower for a higher number of scatterers.However, this dependency diminishes as the level of the false alarm is relaxed.
Figure 3b shows f D against iCV while fixing f F at 10 −3 for single and double scatterers, for different number of acquisitions in the stack, M ∈ {25, 35, 50, 75}.We observe slight dependency of f D on M, though it tends to diminish as the number of acquisitions in the stack grow larger.The aforementioned simulations have been performed in the absence of clutter.We repeat them next with varying levels of clutter, expressed in terms of the signal-to-clutter ratio (SCR): d 2  2 /σ 2 n .Samples to simulate clutter are generated as instances of zero-mean Gaussian noise with variance σ 2 n .Figure 4a shows the iCV observed for the case of single and double scatterer for three different, but fixed, levels of SCR ∈ {6, 3, 0} dB.As expected, the iCV decreases with decreasing SCR.The case of SCR = 0 dB, i.e., when the intensity of the deterministic backscatter from point scatterers equals that of the clutter, contradicts the assumption used in deriving Equation (43).Nonetheless, we perform the simulation as a worst-case analysis.Figure 4b shows f D against the iCV for this case.The plots shown are nearly identical to those shown in Figure 3a.This is an auspicious finding as it implies that, for fixed levels of iCV, f D can be characterized nearly independently of the origin (additive or multiplicative) and level of noise.
The simulation results in Figures 3 and 4 collectively imply that when the number of acquisitions are sufficiently large and the false alarm setting is not too strict, the empirically estimated iCV can be considered to fully characterize the f D , even in the presence of clutter.

Methods
This section presents the overall methodology adopted for the interferometric and tomographic processing of a real interferometric data stack.The models discussed in the previous section form the basis of this methodology.The data undergoes several preprocessing steps.A reference scene is selected, and a multilooked intensity image of the reference scene is used to geocode and coregister all the acquisitions in the stack.An external digital elevation model (DEM) is used in the process [51,52].A suitable reference point is selected to compute double-differenced interferograms.

Interferometric Processing with IPTA
We use the IPTA [3,8] framework for the PSI processing, whereby parameter estimation and phase calibration of the data stack are performed side by side using an iterative approach to least squares regression.An initial list of PS candidates is prepared on the basis of high temporal stability of the backscattering and low spectral diversity.The phase model assumed is as given in Equation (1).Point differential interferograms are obtained by subtracting the topographic phase computed using the DEM.A multiple linear regression is used for each candidate to obtain an initial estimate of s and v, as well as the phase unwrapping integer, p.The quality of the estimates is assessed in terms of the root-mean-square (RMS) phase deviation, σw of the residual phase.At the initial stage, atmospheric phases in each interferometric layer have not be corrected, and the possible temperature-induced phase variations of candidates on structures experiencing thermal expansion have also not been accounted for.Therefore, the residual phase typically exhibits a high dispersion.The PS candidates for which σw is higher than a pre-selected threshold, σ c are masked out.The residue of the remaining candidates is analyzed further.Assuming the atmospheric phase screen (APS) to be spatially low-frequency and temporally uncorrelated, we estimate it by spatial filtering and unwrapping of the phase residue in the neighborhood of the candidates that satisfied the quality criterion.The estimated APS is subtracted and point-differential interferograms are re-computed for the full list of PS candidates, and this time the phases related to the initial estimates of residual height, linear deformation and the atmospheric phase are subtracted as well.The resulting point differential interferograms are unwrapped and the regression is iterated.It is expected that the quality of the candidates would improve since an estimate of the atmospheric phase has been subtracted prior to the regression.σ w is computed again for all the candidates, and compared against σ c to mask out those with relatively low quality.For the retained candidates, the newly estimated regression coefficients (residual height and deformation velocity) act as 'corrections' on the previous estimates.The new phase residue is added to the previous estimate of the atmospheric phase, re-filtered and unwrapped to give a new estimate of the atmospheric phase.The process is iterated several times.In this way, there is progressive improvement in the quality of the estimates in consecutive iterations.For more details on various time-series processing strategies using the IPTA framework, the interested readers are referred to earlier works [3,8,9,53].
For the candidates that are potentially undergoing thermal expansion, another regression-based routine is used that models it assuming that the corresponding phase variations are linearly dependent on the temperature changes [54][55][56].The estimated regression coefficient is the phase-to-temperature sensitivity, η.Further details are available in the earlier work in [35].
After several iterations, the APS is well isolated and we obtain iteratively-refined estimates of the parameter vector p for the PS candidates that satisfy the quality criterion.Assuming that these PS are of sufficiently good quality that the limiting case of von Mises distribution for the phase residuals being approximated by a linear normal distribution is justified, we compute the sample coherence threshold corresponding to σ c using Equation ( 24) as follows: In turn, the corresponding probability of false alarm (theoretical) is computed using Equation ( 18).It is important to mention here that the aforementioned assumption is not mandatory to choose the threshold; in fact, a threshold can be set directly on the coherence (as typically done for interferometric processing) [2,12,57].In our context where we perform PSI processing with the IPTA toolbox (which allows quality assessment in terms of the residual phase statistics), the relation in Equation (52) provides a means to compute the coherence threshold corresponding to the quality criteria in our PSI processing.

Single-Look Differential SAR Tomography with Extended Phase Model
Prior to tomographic inversion, the interferometric data stack requires a precise phase calibration.For the pixels containing PS, we already have an estimate of the atmospheric phases from the PSI processing.Given a sufficient distribution of the PS over the imaged scene or the region of interest, we interpolate these phases over the surrounding pixels that may or may not have been PS candidates.Single-look differential tomographic inversion is applied for each pixel.The extended phase model, given in Equation (27), is used to set up the steering vectors.The reflectivity profile, α (p) is estimated as a function of the unknown parameters.Scatterer localization and parameter estimation for a maximum of two scatterers in each resolution pixel is performed, as stated in Equations ( 35) and (36).The amplitude of the estimated reflectivity is compared against a threshold for each potential scatterer to accept or reject the null hypothesis.
We propose to set the detection threshold in such a way that the desired probability of false alarm from PSI processing is carried forward for the detection of coherent scatterers at this stage.Equating the Equations ( 18) and ( 40), we set the detection threshold for tomography as follows: In turn, the decision between H 1 and H 0 is made for each candidate as follows: |α| (p) In this way, the same quality criterion that is used for setting the threshold T γ c in the PSI processing also determines the threshold for scatterer detection in tomographic processing.Hence, a consistency is achieved for the synergistic use of tomography as an add-on to PSI.
It is to be noted that Equation ( 53) is independent of how the threshold T γ c for PSI processing was selected, whether as a direct choice on the coherence, or using the standard deviation of the residual phase according to Equation ( 52) under the assumption of linear normal distribution of the residual phases for the PS.Therefore, this assumption is not a limiting factor for the application of the proposed detection strategy in general.

Data
The interferometric data stack used in the work comprises 50 TerraSAR-X stripmap acquisitions over the city of Barcelona, Spain in repeated passes.It is the same stack as used in our earlier work in [27].The temporal span of the acquisitions extends from 2007 to 2012.The images have been oversampled by a factor of 2 to allow for more accurate coregistration.The resolution in range and azimuth is 1.2 m and 3.3 m, respectively.The orthogonal component of the total spatial baseline is 503 m, which provides resolution in elevation axis of ∼19 m.The distribution of the spatial and temporal baselines, as shown in Figure 5a, is highly non-uniform.The corresponding 2-D point spread function (PSF) is shown in Figure 5b.The PSF represents the impulse response of the tomographic system for the given distribution of the baselines, for an ideal point scatterer at zero elevation and with no deformation.The footprint of the acquisitions in map coordinates is shown in Figure 5c.Apart from a dense urban stretch, some part of the viewed scene extends over the Balaeric sea.

Results on Real Data
This section presents the results obtained on the real interferometric data stack introduced in the previous section.

Interferometric Processing
An initial list of PS candidates was prepared on the basis of low spectral diversity and high stability of the backscattering amplitude that is characteristic of single dominant scatterers [3].There was no candidate in unexpected areas, such as the water surface or radar shadows.After several iterations of the least squares regression within the IPTA framework, as outlined in Section 3.1, a subset of the initial candidates is retained such that σ w ≤ σ c = 1.1 rad for each candidate.Figure 6 shows these candidates from the last iteration.These are 936,649 in number, and spread over an area of nearly 4 km 2 .In sub-figure a-c, the color coding represents the estimated parameters, namely residual height, deformation velocity in the LOS and phase-to-temperature sensitivity, respectively.The sample coherence for these candidates is as shown in sub-figure d.
Corresponding to σ c = 1.1 rad, the coherence threshold T γ c = 0.55 according to Equation ( 52), and the theoretical probability of false alarm according to Equation ( 18) is 3.3 × 10 −7 .As stated in Section 3.1, the use of Equation ( 52) to convert a threshold in terms of residual phase standard deviation to corresponding threshold on coherence requires the assumption that the von Mises distribution can be approximated by a linear normal distribution (for the case of ideal, noise-free PS, with κ → ∞).In order to assess the suitability of this assumption, we require estimates of the concentration parameter.Using Equations ( 20) and ( 22), E {| γ|} I 1 (κ) I 0 (κ) ; to estimate κ, this expression needs to be inverted, which involves inversion of the ratio of modified Bessel functions (first kind) of first and zero order.We do not have a closed-form expression for such an inverse relation; we use the following piece-wise defined approximation [43,58]: The concentration parameter is estimated for each PS, and a histogram of the parameters is shown as an inset in Figure 6d.The mean and the median values are 5.4 and 4.1, respectively.In existing literature in the field of directional statistics, we can find precedence where concentration parameters greater than 2 are considered reasonable to approximate von Mises distribution as wrapped normal distribution (i.e., linear normal distribution wrapped between −π to π rad) [58].

Tomographic Processing and Empirical Analysis of False Alarms
The APS isolated in the IPTA-based PSI processing is extrapolated over the scene and compensated for over the entire scene in each layer of the interferometric stack.In this way, each pixel is considered to be phase calibrated so that tomographic inversion can be applied next.Given that the city of Barcelona has several high-rise buildings, the elevation extent, I s is set as [−60, 300] m.The parameter space for the deformation parameters is as follows: I v ∈ [−10, 10] mm/yr and I η ∈ [−1, 1] rad/K.The discretization in each dimension is 1/2.5 times the Rayleigh resolution, followed by a local refinement of the estimated reflectivity around the two candidate peaks at one-tenth the resolution.Using Equations ( 52) and (53), and keeping σ c = 1.1 rad, we threshold the reflectivity of the two candidates to perform the detection process.The point cloud of single scatterers thus detected is shown in Figure 7. 1454 false alarms occur over the water surface.52) and ( 53)).The color coding represents the estimated height.Some false alarms can be seen over the water surface , as highlighted in the inset.
A significant portion of the viewed sea extends over the sea, which is favorable in our context as it can be used as a test bed to conduct an empirical analysis of the false alarm rate.We perform sample coherence-based detection, as well as tomographic inversion and detection, for the range-azimuth pixels over the sea and observe the variation of the false alarm rate.These pixels constitute 1.4 million independent resolution cells.The results are shown in Figure 8.The solid lines in the figure represent different cases of tomographic inversion and detection: (1) [α (s, v, η); 3-D inv.]: 3-D inversion and detection on the reflectivity, α retrieved as a function of elevation (s), deformation (v) as well as thermal expansion (η) where the support in each dimension is as for the results shown in Figure 7, (2) [α (s, v); 2-D inv.]: 2-D inversion i.e., thermal expansion is not considered, (3) [α (s); 1-D inv.]: 1-D inversion, whereby the reflectivity is retrieved only along the elevation profile, (4) [α (s); reduc.supp.]:1-D inversion with the elevation support reduced to [−25, 50] m, and (5) [α (no fitting)]: 1-D inversion without the maximizations to detect peaks in the reflectivity, i.e., no fitting is performed in the parameter space to estimate the unknown elevation and deformation parameters.(6) [ γ (no fitting)]: The dot-dashed line represents the PSI case whereby the thresholds are applied on the sample coherence without any parameter fitting.(7) [ γ (analytical)]: The black curve with diamond symbols shows the probability of false alarm (theoretical) according to the Equation (18).The bottom x-axis in the figure shows the detection thresholds, T γ and T α (normalized between 0 and 1 as per Equation ( 53)), while the top x-axis shows the equivalent standard deviation of the residual phase according to Equation (52).The area shaded in gray indicates the region in the figure where the results may not be sufficiently accurate due to a limited number of independent range-azimuth resolution cells over the water surface.Given we have only 1.4 million of these cells, and assuming the test statistics are normally distributed over the scene, we can estimate a probability of false alarm no less than 1.1 × 10 −3 with a relative absolute error of 5% for 95% of the time [49].(no fitting) fitting) (analytical, eq.18)  Figure 9 shows the point cloud of single scatterers obtained with tomographic inversion and detection with σ c = 1.0 rad.In comparison with Figure 7, we can see a reduction in the false alarms.Now, we observe only 194 false alarms.3-D tomographic inversion has been applied with the same support in each dimension as for the results shown in Figure 7.We have estimates of height, deformation velocity as well as the phase-to-temperature sensitivity.Figure 10 shows the point cloud of double scatterers obtained with the same threshold.They are separated as lower and upper scatterers, according to the estimated height for each of the two scatterers in layover.2.14 × 10 6 single scatterers and 1.01 × 10 4 double scatters (lower + upper) have been detected.The inset in the sub-figures in Figure 10 shows a commercial complex, namely Diagonal Mar, in focus.The red polygon encloses a high-rise building, which is partly in layover with the roof of a nearby building.These are the same test sites as in our earlier work in [27].7 where a more relaxed detection threshold (corresponding to σ c = 1.1 rad) is used, fewer false alarms are observed here, as highlighted in the inset.

Discussion
This section provides an itemized discussion of the results presented in the previous section.

Interferometric Processing
The PSI solution, as shown in Figure 6, provides a good coverage over the viewed scene, which is typical with high resolution X-band interferometric imagery over urban areas such as Barcelona city [59,60].The PS heights fit reasonably with actual 3-D structures, as shown for selected buildings in our earlier work in [26,27].The PSI solution reveals deformation along the shoreline, which was partly observed in [59] as well.Several PS on high-rise buildings show temperature-dependent phase variations, which can be attributed to thermal expansion of the structures [30,35,46,61,62].The observed coherence is high, and the estimated concentration parameters are all non-zero.With reference to Figure 1, the fact that the mean and the median value of κ are greater than 3 substantiates the assumption of linear normal statistics for the PS (since the approximation of von Misses as linear normal distribution is accurate to within 5% error on average).
Interestingly, we do not observe false alarms over the sea patch in the scene.This is due to the fact that we have used high stability of the backscattering amplitude and low spectral diversity as pre-classifiers to set up the initial PS candidate list.These classifiers are proxies for temporally coherent, single dominant scattering; therefore, they already preclude PS candidates from appearing on the water surface.Hence, no PSI solution has been sought (no regression fitting) on the pixels over the sea patch.In the context of tomography, these pre-classifiers cannot be used since they would tend to reject double scatterers as well.

Tomographic Processing and Empirical Analysis of False Alarms
We applied tomographic processing over the entire scene, regardless of any surface classification.The point cloud shown in Figure 7 is obtained using the same cut-off phase standard deviation, σ c = 1.1 rad, as for the iterative least squares based PSI processing.Nevertheless, several false alarms are visible over the sea patch.A simple mask (based on SAR multi-look intensity with spatial constraints for example) could have allowed us to remove the sea patch from the processing, but we choose to show these false alarms to highlight that similar false alarms may arise (due to noise) within the urban stretch as well though they may remain unnoticed.
Figure 8, which shows the results of a false alarm analysis exclusively conducted over the sea patch, reveals that the false alarm rates can typically be higher in practice in comparison with the theoretical probability of false alarm (as the area under the upper tail of Rayleigh distribution).The maximizations (Equations ( 35) and ( 36)) allow degrees of freedom to fit the data; when the noise is fit incorrectly with the data model, it may lead to a false alarm.The false alarm rate can be seen to decrease from 3-D to 2-D inversion, as reducing the dimensionality reduces the degrees of freedom to fit the data.Similar reduction in false alarms is observed when moving from 2-D to 1-D inversion, or when we reduce the support of the elevation in case of 1-D inversion.These findings imply that in case some a priori information is available-e.g., if significant thermal expansion is not expected (as is usually the case for buildings of low height [27]), or if the support of deformation velocity can be reduced on the basis of local leveling measurements, or if the support for height corrections can be reduced given a digital surface model is available-then a reduction in false alarm rate can be achieved in practice.
Figure 8 also shows the case where no parameter fitting is performed, for both tomography as well as sample coherence based detection.The latter case, i.e., [ γ (no fitting)], matches closely with the theoretical relationship in Equation (18), indicating that the area under the upper tail of the PDF of | γ| approaches that of a Rayleigh distribution.However, in the former case, i.e., [α (no fitting)], it can be observed that the estimated false alarm rate is slightly lower than the probability of false alarm according to the analytical expression for MICC-based detection, in turn implying deviation from the statistics of a Rayleighian process.It can be explained following the findings in an earlier work in [13].In this work, a generalized likelihood ratio test (GLRT) was compared against MICC for scatterer detection in the presence of additive noise with Gaussian statistics.It is to be noted that in our case the false alarm analysis is conducted on cells over the water surface; therefore, the origin of noise in the observed SLC values lies in the backscattering characteristics (rather than phase mis-calibration).In this particular context, an additive noise model is appropriate, and, consequently, the detection for a scatterer under Equation (54) in our work becomes identical to the GLRT in [13].It was found in [13] that the GLRT provides a lower probability of false alarm compared to MICC (as we observed).For a discussion on the performance analysis of radar detectors where the actual PDF of the amplitude of complex-valued noise/clutter deviates from Rayleigh statistics, interested readers are referred to [63][64][65][66].
Figures 9 and 10 show the single and double scatterers, respectively, detected with σ c = 1.0 rad.As expected, we observe fewer false alarms, and at the same time fewer scatterers are detected.Double scatterers constitute <1% of the total scatterers detected over the scene.The gain in deformation sampling due to double scatterer detections [27], relative to the PSI solution, are around 2% for Diagonal Mar complex and 4% for the selected building marked in red, respectively.If the threshold is relaxed to σ c = 1.1 rad, the gain improves to 6.4% for Diagonal Mar and 17% for the individual building.The interferometric data stack and the test sites in this work are the same as in our earlier work in [27].The detection strategies are, however, different.The sequential GLRT with cancellation (SGLRTC), as proposed in [24], was used for hypothesis testing in the earlier work.The quality of the detected scatterers was empirically evaluated only after the detection, and in turn compared with the quality of the PS (obtained independently in the prior PSI processing).In other words, the detection threshold for hypothesis testing had to be adjusted a posteriori to achieve comparable quality.The results thus obtained in [27] show a gain in deformation sampling of around 2.5% for Diagonal Mar complex and 10% for the selected building.On the other hand, the detection strategy proposed in this work allows the use of quality criterion during the hypothesis testing itself.Nonetheless, it needs to be noted that the SGLRTC and the proposed strategy are not directly comparable.SGLRTC explicitly assumes an additive noise model for SAR tomography, thus it cannot formally address multiplicative noise arising due to phase instabilities such as atmospheric disturbances.Moreover, it is a subspace method where the first scatterer is canceled out before a second scatterer is searched for [24].Therefore, the test statistics (and the corresponding threshold settings) for double scatterer detection under the proposed detection strategy are not the same as in SGLRTC.

Conclusions
In the context of SAR tomography as an add-on to PSI to potentially improve deformation coverage, following the directions set in earlier works in [12,27,31], this paper reports the application of a detection strategy that allows for extending the same quality considerations to tomography as used in the prior PSI processing.In interferometric processing, the quality is typically assessed on the basis of the residual phase, either in terms of the phase dispersion (phase standard deviation) or the ensemble coherence computed using the residue of the fit.In both cases, under the proposed detection strategy, the quality parameters can be used to set up the threshold for hypothesis testing of coherent scatter candidates following tomographic inversion.Moreover, the theoretical probability of false alarm remains the same between the PSI and tomography.The paper also highlighted that while the instabilities in phase are typically modeled as additive noise, their impact on tomography is multiplicative in nature.The experiments performed in this work with simulated data consider both multiplicative noise as well as additive disturbances (clutter) in the tomographic model.It is shown that the inverse coefficient of variation is a suitable parameter to assess the probability of detection, irrespective of the origin of noise.The proposed detection strategy is also tested on real data.An assessment of the variation of the observed false alarm rates against the thresholds set according to the proposed detection strategy has been conducted.An interferometric data stack comprising 50 Terra-SAR-X acquisitions over the city of Barcelona, Spain is used.Single-look beamforming for 1/2/3-D tomographic inversion, depending on whether the phase model used considers only the scatterer height, or height plus deformation velocity, or additionally thermal expansion, is performed.The results show that higher dimensionality and larger support sizes in each dimension lead to higher false alarm rates due to larger parameter space that may incorrectly fit noise to the data model.These results also suggest that in case a priori information can reduce the dimensionality and/or support sizes, it should be adopted by the user to reduce the false alarm rate in practice.For the case of 3-D tomographic inversion, with detection thresholds set in accordance with residual phase standard deviation below 1.1 rad for the prior PSI processing, the empirically estimated false alarm rate is <1.1 × 10 −3 .The gain in deformation sampling (due to layover resolutions) is 17% for a selected high-rise building.For a commercial complex in Diagonal Mar locality, it is 6.4%.As a whole, the number of double scatterers detected in the urban scene are <1% of the total detected scatterers.These results show that, for urban areas like Barcelona, when using interferometric data stacks comprising the typical stripmap products, the application of SAR tomography as an add-on to PSI is mainly useful for a detailed analysis of selected urban zones or individual buildings in layover.

Figure 1 .
Figure 1.Estimates of the coherence magnitude obtained with 10 5 Monte Carlo iterations assuming the residual phases have a von Mises distribution with concentration parameter, κ.Each solid line indicates the estimates for a specific number of acquisitions, M in the data stack.The vertical bars represent ± 1-σ from the mean.The dashed line shows the coherence magnitude under the assumption that the residual phases follow a linear normal distribution, cf.Equation (24) (assuming σ 2 w = κ −1 ).

Figure 2 .
Figure 2. Empirically estimated inverse coefficient of variation (iCV) of the test statistic α against concentration parameter for von Mises distributed phase residuals, for different number of scatterers, Q in the same resolution cell.The dashed lines enclosing the gray region indicate the theoretical bounds on the iCV (cf.Equation (51)) , where ν 1 = I 1 (κ)I 0 (κ) (Equation (22)).

Figure 4 .
Figure 4. Numerical analysis of the inverse coefficient of variation (iCV) of the test statistic α when point scatterers are embedded in different clutter levels, for M = 50 acquisitions.(a) iCV against concentration of the phase residuals for different levels of signal-to-clutter ratio (SCR) and number of scatterers, Q ∈ {1, 2}; (b) probability of detection against iCV for fixed levels of false alarm, f F ∈ 10 −2 , 10 −3 , 10 −4 , and Q ∈ {1, 2, 3}.

Figure 5 .
Figure 5. Data characteristics.(a) distribution of spatial (orthogonal component) and temporal baselines; (b) 2-D point spread function (PSF); (c) footprint of the reference acquisition over Spain.

Figure 6 .
Figure 6.PSI solution obtained with iterative least-squares regression-based processing using the interferometric point target analysis (IPTA) toolbox.The colored dots are the PSs identified in the PSI processing.(a) estimated height, relative to the WGS-84 reference ellipsoid; (b) deformation velocity in the line-of-sight; (c) phase-to-temperature sensitivity; (d) sample coherence, and histogram of the estimated concentration parameter (shown as inset).

Figure 7 .
Figure 7. Point cloud of single scatterers obtained with differential SAR tomography.The detection threshold is set corresponding to σ c = 1.1 rad under the proposed detection scheme (see Equations (52) and (53)).The color coding represents the estimated height.Some false alarms can be seen over the water surface , as highlighted in the inset.

Figure 8 .
Figure 8. False alarm rate observed over the sea patch in the viewed scene at different detection thresholds.The colored solid lines represent the case of 3/2/1-D tomographic inversion.The detection is performed on the retrieved reflectivity, |α| according to Equation (53).The dot-dashed lines shows the case of PSI whereby the detection is performed on the sample coherence, | γ| without fitting any phase model to the observed interferometric phases.

Figure 9 .
Figure 9. Point cloud of single scatterers obtained with differential SAR tomography.The detection threshold is set corresponding to σ c = 1.0 rad under the proposed detection scheme, see Equations (52) and (53).(Top) Estimated height, relative to the WGS-84 reference ellipsoid.(Middle): Deformation velocity in the line-of-sight.(Bottom) Phase-to-temperature sensitivity.In comparison with Figure 7 where a more relaxed detection threshold (corresponding to σ c = 1.1 rad) is used, fewer false alarms are observed here, as highlighted in the inset.

Figure 10 .
Figure 10.Point cloud of double scatterers obtained with differential SAR tomography.The detection threshold is set corresponding to σ c = 1.0 rad under the proposed detection scheme.(Top) Estimated height, relative to the WGS-84 reference ellipsoid.(Middle) Deformation velocity in the line-of-sight.(Bottom)Phase-to-temperature sensitivity.The left column shows the lower layer and the right column shows the upper layer of the double scatterers, respectively.The inset focuses on a commercial complex (Diagonal Mar).The red polygon encloses a single building, part of which is in layover with a nearby building of shorter height.