Combination of Super-Resolution PSI and Traditional PSI by Identification of Homogeneous Areas

The performance of Persistent Scatterer Interferometry (PSI) depends heavily on Persistent Scatterer (PS) density. In order to increase PS density, we can apply Super-Resolution reprocessing algorithms in PSI. Involving the reprocessing algorithms and the peak-detection-based Persistent Scatterer Candidate points (PSCs) selection method, the full PSI chain is referred to as Super-Resolution PSI (SR-PSI). The implementation of the Super-Resolution reprocessing algorithm, however, is computationally intensive, which makes SR-PSI time-consuming. In this work, we propose to improve the efficiency by constraining the Capon-based reprocessing to the non-homogeneous areas (e.g., urban areas). We notice that the Capon algorithm performs similarly as the Fourier-based algorithm for homogeneous regions (e.g., grassland), thus we can use Single Look Complex (SLC) images for these areas. With the Coefficient of Variation (CV) as the index, we divide the full image into two classes: homogeneous areas, for which we select PSCs from the original stack, and non-homogeneous areas, for which we extract PSCs from the Capon-based reprocessed images. Then we combine the PSCs of both cases for further PSI processing. We applied the combination method to a stack of TerraSAR-X data. The results show that the proposed approach is more computationally efficient than the original SR-PSI with the effectiveness uncompromised, especially for applications aiming at the urban deformation.


I. INTRODUCTION
Since it was first proposed by [1], [2] in 2000, Persistent Scatterer Interferometry (PSI) has been intensively used for the estimation of relative surface deformation in various areas, e.g., landslide monitoring [3], building deformation assessment [4], with accuracy down to millimeter scale. Since higher PS density enables us to retrieve more details concerning the deformation phenomena, the performance of PSI depends heavily on the PS density. The PS density is roughly proportional to the bandwidth of the SAR system (the wider the bandwidth, the finer the resolution) [5], thus data from higher-resolution SAR systems such as TerraSAR-X [6], [7] and Cosmo-Skymed [8], [9] provides higher PS density. A more feasible and economical approach, though, is through The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . the application of super-resolution reprocessing algorithms [10]. The corresponding processing chain, which is called Super-Resolution PSI (SR-PSI), is also proposed by the reference [10]. Compared with the regular PSI chain, the main change of SR-PSI is that it involves the Super-Resolution reprocessing algorithms and applies the peak-detection-based Persistent Scatterer Candidate (PSC) selection method.
The idea of SR-PSI comes from the applications of modern spectral estimation algorithms in improving the resolution and suppressing sidelobe levels [11]- [13]. A range of spectral estimation algorithms are available, and the Capon algorithm seems to be the most promising one concerning PSI [10]. Thus we only focus on the Capon algorithm in SR-PSI. The Capon algorithm, however, is computation-consuming, which may make it difficult to deploy in practice, especially when a large number of SAR images are involved. One attribute of the Capon algorithm is that the performance is similar as the original SLC for homogeneous regions (agricultural area, for example) [13]. Therefore, we are able to identify and filter out these areas in the Capon-based reprocessing, to save computation. First, we segment the full image into overlapped square-shaped chip images, which are classified into two groups: homogeneous areas and nonhomogeneous areas. For the homogeneous areas, we select PSCs from the original images. For the non-homogeneous ones, we run the Capon-based reprocessing algorithm on the chip images, followed by a peak-detection based PSC selection approach [10]. We combine PSCs in both cases for further PSI processing, mainly including PSC network construction, temporal ambiguity resolution, spatial ambiguity resolution and atmosphere phase estimation. Finally, we obtain the time-series deformation, the deformation velocity and the height of each Persistent Scatterer (PS).
A range of literature, e.g., [14]- [18], have studied the issue of identifying the homogeneous areas. Previous research has proven that the simple Coefficient of Variation (CV), which is defined as the division between the standard deviation and mean of the amplitudes, can be used as an effective index of homogeneous areas [14], [15]. If we set a specific threshold of CV, the images with lower CV will be labelled as homogeneous areas. In order to calculate the threshold, we shall analyze the distribution of the squared mean amplitude and the distribution of the standard deviation of the amplitudes, respectively. Then both their ranges and the threshold can be readily derived.
We obtained 43 TerraSAR-X images around Rotterdam, the Netherlands as test data. We then applied the proposed approach to the data to evaluate the proposed workflow. Section II presents the combination of SR-PSI and traditional PSI, using CV as the index of homogeneous areas. Section III discusses the results of the experiment which applied the proposed method to TerraSAR-X data. Section IV draws the main conclusions.

A. COEFFICIENT OF VARIATION AS THE INDEX OF HOMOGENEOUS AREAS
Homogeneous areas are the regions with consistent spatial statistical characteristics. Usually, homogeneous areas are comprised by distributed scatterers, which can be explained in Figure 1. In Figure 1, the reflected echo is the sum of echoes from a large number of small discrete scatterers illuminated by the beam: where A φ indicate the amplitude and the phase respectively, N is the number of discrete scatterers. The in-phase component Acosφ and orthogonal component Asinφ can be modeled as a normal distribution with a mean value of 0 and a standard deviation of σ [19], where σ indicates the characteristic of distributed scatterers. The CV has been used as an effective indicator of scene heterogeneity [14], [15], which is defined as where σ A and µ A are the standard deviation and mean of the amplitudes. It is easy to conclude that the chip images with lower CV can be recognized as homogeneous areas: In the following subsection, we shall analyze the distribution of σ A and µ A , followed by the derivation of the range. Then the threshold C 0 can be readily calculated.
The intensity I of homogeneous areas is the sum of two squared independent normal distribution, thus the probability density function of I can be modeled by an exponential distribution. Furthermore, I can be modeled by: with E(I ) = 2σ 2 , Var(I ) = 4σ 4 . In the case of L-look processing, i.e., averaging L independent observations for the each individual sample is performed, the averaged intensity I is gamma distributed ( [19], [20]): where (·) is the gamma function, (L) = (L − 1)!. The mean and variance of intensity I can be expressed by E(I ) = 2σ 2 , Var(I ) = 4σ 4 /L. Then the amplitude is characterized by the square root of the gamma distribution. Using the variable relation, where A denotes the amplitude, and inserting (5) into (6), yields the distribution of A: It's a Nakagami distribution [21]. We calculated the probability density (PSD) function with different L and σ , which is presented in Figure 2. We can see that the PSD of the amplitude is related to the number of looks and the standard deviation σ . The mean and variance of A can be VOLUME 8, 2020  expressed by A typical size of the chip images is 16 × 16 pixels. Hence, suppose the number of samples contained in the individual chip images is N , which usually satisfy N 30. CV can be expressed by Likely, we can obtain where E(A), Var(A) can be expressed by (8). The distributions of N i=1 A 2 i andĀ are given by (10) and (11), respectively. With significance level set as 0.01, the range of Thus, with significance level of 0.01, FIGURE 5. Squared mean intensity image of the original TerraSAR-X stack. We can observe several kinds of terrain, including urban areas, water, and grassland, from the varied scattering characteristics. The part marked by the dashed line is shown in detail in Figure 11.
the threshold can be set as where E(A), Var(A) are given by (8). From (8) and (12), we can observe that the threshold is related to the number of looks and the number of samples. In the case of L = 1, N = 256, E(A) = σ √ π/2, Var(A) = σ 2 (4 − π)/2, and C 0 = 0.89. Similarly, we can calculate C 0 as 0.79 and 0.75 for significance level 0.05 and 0.1, respectively.

B. COMBINATION OF SR-PSI AND TRADITIONAL PSI WITH CV-INDEX
A detailed discussion of SR-PSI has been presented in [10]. By applying the Capon-based reprocessing algorithm, it obtains the stack of super-resolution images. The output of the Capon algorithm can be given bŷ where H represents the Hermitian operator, a is the 2-D Fourier matrix, L 1 , L 2 is the number of snapshots,R is the sample covariance matrix, which is given bŷ and g can be expressed by where z l 1 ,l 2 represents the vectorization of sub-matrices of the 2-D data matrix. (13), (14), (15) define the computation of the complexvalued estimation for each frequency pair (ω 1 , ω 2 ). While a direct implementation of the Capon estimator is computationally demanding, several computationally efficient algorithms have been proposed to lighten the computational burden [22]- [24]. These available efficient schemes generally exploit the structure of the involved matrix, followed by fast calculation through means of the 2-D fast Fourier transform (FFT). Nevertheless, even if we take the existence of these algorithms into account, Capon based re-focusing remains computationally intensive.
We propose to combine the SR-PSI and the traditional PSI by identifying the homogeneous areas, which is based on two nationalities: 1) The performance of the Capon algorithm is similar to the Fourier transform for homogeneous areas [13].
2) The corresponding objects of PS are usually artificial features, which are not distributed in homogeneous areas. Thus we cannot select PSCs or can only select sparse PSCs from homogeneous areas.
Therefore, we can filter out these homogeneous areas in the Capon-based reprocessing.
Since the Capon algorithm is memory-intensive, we usually employ a chipping and mosaicking strategy in the Caponbased reprocessing [11], [13]. It is demonstrated in Figure 3, where an overlapping rate of 50% is typically used to mitigate the edge effects [13]. We adapt the chipping and mosaicking processing by classifying the individual chip images into homogeneous area and non-homogeneous areas. Since the traditional PSC selection method cannot be applied directly to the super-resolution reprocessed images, we can use the peakdetection-based PSC selection method proposed by [10]. On the non-homogeneous chip images, we run the Caponbased reprocessing algorithm, followed by a peak-detection based PSC selection [10]. For the homogeneous chip images, we can select PSCs from the original SLC (or SLC windowed by some functions to taper the sidelobe levels), by setting an upper limit on the normalized amplitude dispersion [2].
Additionally, we need to include the neighboring chip images of the non-homogeneous areas in the actual Capon-based reprocessing.
We present the full workflow in Figure 4, which includes a regular InSAR process, PSC selection, and other PSI processing. The regular InSAR process includes master image selection, image coregistration and reference phase subtraction. In the last step, the extracted PSCs are exported for PSC network construction, temporal ambiguity resolution, spatial ambiguity resolution and atmosphere phase estimation [2], [25]. Finally, we can obtain the time-series deformation, the deformation velocity and the height of each PS. Compared with the workflow of the original SR-PSI in [10], the main improvement is that we constrain the Capon-based reprocessing to the non-homogeneous areas.

A. DATA DESCRIPTION
To evaluate the performance of the proposed method, we chose a test site around Rotterdam, the Netherlands. We obtained 43 TerraSAR-X images, corresponding to 43 epochs ranging from January 15, 2016 to May 22, 2018. Choosing the image of January 1, 2017 as the master image, we processed these images with the Delft Object-oriented Radar Interferometric Software (DORIS) [26] and acquired a stack of interferometric images with size of 2048 × 2048 pixels. Figure 5 shows the squared temporal mean of the intensity, from which we can observe several kinds of terrain, including the urban area, water, and grassland.

B. CV AS THE INDEX OF HOMOGENEOUS AREAS
With the size of the chip image set to 16 × 16, we divided the full image into 128 × 128 chip images. We calculated the temporal average CV for each individual chip image, which is shown in Figure 6. We can see that most   of the grassland and water are indicated by lower average CV. We can also observe high average CV in the urban areas.
Since homogeneous chip images usually has a lower image entropy [27], we use image entropy to cross validate the effectiveness of CV. For any given vector − → Y , the image entropy is defined as where .
We calculated the image entropy. Subsequently, we drew the scatter plot of the chip images using the CV and the image entropy as the x-axis coordinate and y-axis coordinate, respectively, as shown in Figure 7. From the scatter plot, we can observe the CV and the image entropy are highly negative correlated, which is further confirmed by their correlation coefficient, -0.98. We calculated the linear function between the CV and the image entropy through the least square approach. The relationship is expressed by: where S is the image entropy and C represents the CV. Both Figure 7 and (18) show that the chip images with lower CV have high image entropy, which suggests that CV can be used as an effective index of homogeneous areas.

C. COMBINATION OF SR-PSI AND TRADITIONAL PSI USING CV-INDEX
To exploit the relationship between the PSC number and the average CV, we reprocessed the images by the Caponbased reprocessing algorithm and calculated the PSC number for individual chip images. With the threshold of normalized amplitude dispersion set to 0.25, we can select 9480 PSCs from the Capon-based reprocessed stack compared to 5557 PSCs from the original stack, which attains a 71% increment. We refer to the chip images where the PSC number increase as the Capon-effective chip images. Conversely, we refer to other chip images as Caponineffective chip images. The histograms of the average CV are presented in Figure 8. From the left subplot, we note that the average CV of the Capon ineffective chip images mostly fall into lower value section, which indicates that most of them are homogeneous areas. In the right subplot, the average CV of the Capon effective chip images are relatively higher. It is easy to conclude that choosing the threshold of CV implies balancing the overall computation efficiency against the effectiveness of the Capon algorithm. We refer to the number of increased PSC for each chip image as PSC increment. Subsequently, we calculated the PSC increment and the number of the reprocessed images using different CV thresholds, as summarized in Table 1. We use the PSC increment in the case of SR-PSI as the reference. The minimum average CV of the Capon-effective chip images is 0.56. By setting the threshold to 0.56, we can maintain the full effectiveness of SR-PSI. In this case, we need to reprocess 83% of the chip images. With the threshold of CV set to 0.75 (in the case of significance level set as 0.1), we need to reprocess approximately 52% of the chip images while maintaining  97% of the effectiveness. As indicated by Table 1, a higher threshold of CV may lead to the degradation of the PSC increment while improving the efficiency. We should choose the threshold of CV according to the specific applications. We will make back to this in the final discussion. For the rest of this section, we set the CV threshold to 0.75. Moreover, we can observe that some chip images which are mixed by grassland and water also have high CV in Figure 6. However, the effect is negligible due to the relatively minor amount.

D. DEFORMATION ESTIMATION
We calculated the deformation results of the proposed combination. We did a full PSI time-series analysis to the detected PSCs using DePSI [28], obtaining the estimates of the linear deformation rate. To retain as many PSs as possible, we relax the normalized dispersion threshold to 0.45. For the final PS selection, we use the ensemble coherence as a quality estimator, setting the lower limit for selection to 0.8. We obtained 15676 PSs from the original stack and 21759 PSs using the combination method, which are shown in Figure 9. We can observe that most of PSs are distributed in the urban area in both cases, which accords with the nature of PS. The amount of PSs increases by approximately 39%.
The reference paper [10] has validated that the quality of the PSs extracted from the Capon-based reprocessed images are maintained. The main assumption is that the behaviour of the common PSs extracted in both cases should be consistent if the quality is maintained compared with the regular PSI. We found 13950 common PSs through a point matching process [10]. We then compared the deformation velocities, standard deviation (STD) of the residuals between the deformation model and the deformation time series and the heights of the common PSs. The joint 2-D histograms are shown in Figure 10. In the Figure 10, we can observe that most of the common PSs are concentrated along the reference line, which suggests that the quality level is maintained in the application of the combination method.
To visualize the details, we show the deformation velocities of the PSs in the region indicated by the dashed line in Figure 5. In Figure 11, we found 1771 PSs from the original stack and 2562 PSs using the combination method, which is an increment of 45%. Moreover, we can also see some PSs with new deformation velocities, which implies the advantages of the Capon-based reprocessing algorithm.

IV. CONCLUSION
In this paper, we propose to improve the efficiency of SR-PSI by the identification of homogeneous areas. In achieving the results shown in this paper, we used a CV threshold of 0.75, finding that it provided a good balance between the effectiveness of the method, measured in terms of the fraction of all possible newborn PSs that were recovered, and the computational cost. The optimum point is, of course, dependent on the application and on the available computational resources. It is nevertheless interesting to highlight different scenarios that would lead to different re-processing strategies: • For applications targeting the building environment and emphasizing local deformations, we may, indeed, choose a CV threshold aiming at recovering most of potential newborn PSs, knowing that most of these new PSs will appear in areas that already show a high PS density.
• If we are studying larger scale deformation signals, for example, subsidence bowls of some nature, we may actually be starved for PSs in low-PS density areas, which typically have a low associated CV index. In such a case, where a handful additional PS may have a very significant impact, we could choose to lower the CV threshold, but at the same time limit the reprocessing to chip-images with a low number of PSs. The approach proposed in this work provides a path for an efficient use of super-resolution algorithms in combination with PSI analyses. For the particular example shown, the improvement in efficiency was limited as large portions of the TerraSAR-X images used correspond to the built environment. The gain of efficiency can be expected to be much higher when dealing with wider and more sparsely urbanized areas. He is also a Lead Investigator for the STEREOID mission candidate. His current research interests include (In)SAR time series analysis, retrieval from ocean surface currents from radar data, and the development of distributed multistatic radar concepts.
SHAONING LI received the master's degree in geographic information system from the Wuhan University of Science and Technology, in 2013, and the Ph.D. degree from the State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, China, in 2017. He was engaged in scientific research with the Postdoctoral Research Station, Wuhan University, from 2017 to 2019. His research interest includes the remote sensing data processing. VOLUME 8, 2020