High-Quality Pixel Selection Applied for Natural Scenes in GB-SAR Interferometry

Phase analysis based on high-quality pixel (HQP) is crucial to ensure the measurement accuracy of ground-based SAR (GB-SAR). The amplitude dispersion (ADI) criterion has been widely applied to identify pixels with high amplitude stability, i.e., permanent scatterers (PSs), which typically are point-wise scatterers such as stones or man-made structures. However, the PS number in natural scenes is few and limits the GB-SAR applications. This paper proposes an improved method to take HQP selection applied for natural scenes in GB-SAR interferometry. In order to increase the spatial density of HQP for phase measurement, three types of HQPs including PS, quasi-permanent scatter (QPS), and distributed scatter (DS), are selected with different criteria. The ADI method is firstly utilized to take PS selection. To select those pixels with high phase stability but moderate amplitude stability, the temporal phase coherence (TPC) is defined. Those pixels with moderate ADI values and high TPC are selected as QPSs. Then the feasibility of the DS technique is explored. To validate the feasibility of the proposed method, 2370 GB-SAR images of a natural slope are processed. Experimental results prove that the HQP number could be significantly increased while slightly sacrificing phase quality.


Introduction
Landslide is one of the most frequent and harmful geological hazards all over the world. Its monitoring and early warning plays an important part in disaster prevention and reduction. Surface deformations normally occur before the macro failure of natural and engineered slopes. Deformation measurement is of great significance to forecast and early warn landslide hazards [1,2]. Spaceborne synthetic aperture radar (SAR) is a powerful and effective remote-sensing technique for surveying large areas (km 2 ). Main limitations are due to its long repetition period in the order of several days or weeks [3,4]. Ground-based SAR (GB-SAR) is a good complementary tool of spaceborne SAR suitable for monitoring deformations of small areas (from a single building to a large slope), with the advantages of all-day, all-weather, high-accuracy, and fast image acquisition [5,6].
With two radar images acquired at different times and from the same position, their interferometric phases could be explored to measure deformations based on the differential interferometry technique [7]. Interferometric phase is the most crucial element that

High-Quality Pixel
HQPs for natural scenes refer to those scatterers with high phase stability. Based on the PS concept, the idea is to firstly identify those single-point scatterers with moderate amplitude stability and high phase stability, i.e., QPSs. Then the feasibility of the DS technique is explored in GB-SAR interferometry. Although the temporal coherence of those natural scatterers is typically low, the number of pixels sharing the same statistical behavior can be large enough to make it possible for a few of them to exceed the coherence threshold. Together with the PS, QPS, and DS, the spatial density of the HQP utilized for phase analysis could be significantly increased.

Permanent Scatterer
PSs refer to those point-wise scatterers which are almost not affected by the spatial and temporal decorrelation over a long time period. The ADI criterion is commonly used by analyzing the amplitude series of a pixel and thresholding on the ADI. The ADI D A of a pixel can be expressed as where σ A and m A are the standard deviation and mean value of a pixel's amplitude series, respectively. A same numerical simulation with the published research is repeated to show the ADI criterion [10]. For a deterministic point-wise scatterer, its signal model is built as z k = g + n k (k = 1, · · · , K), where g is a pixel's reflectivity. To better analyze the relationship between the dispersion index D A and the phase standard deviation σ V , g is fixed to 1. n is the noise term, and tends to be a complex circular Gaussian distribution with power σ 2 n for both real and imaginary components. The deviation σ n gradually increases from 0.05 to 0.8. For each value of σ n , 5000 estimates of the ADI values of the amplitude series A = |z| are completed. Thirty (K = 30) data are supposed to be available.
The mean values (solid line) and the dispersion (error bars) of the estimatedD A are shown in Figure 1a, together with the values of the phase standard deviation (dotted line) σ V . It can be noted that the ADI D A is a good measure of the phase stability for pixels with high SNR (signal-to-noise ratio). Those pixels exhibiting D A value under a given threshold D T (typically lower than 0.25) can be selected as PSs.
by analyzing the amplitude series of a pixel and thresholding on the ADI. The ADI A D of a pixel can be expressed as where  A and A m are the standard deviation and mean value of a pixel's amplitude series, respectively.
A same numerical simulation with the published research is repeated to show the ADI criterion [10]. For a deterministic point-wise scatterer, its signal model is built as  k k z g n     1, , k K , where g is a pixel's reflectivity. To better analyze the relationship between the dispersion index A D and the phase standard deviation  V , g is fixed to 1. n is the noise term, and tends to be a complex circular Gaussian distribution with power  2 n for both real and imaginary components. The deviation  n gradually increases from 0.05 to 0.8. For each value of  n , 5000 estimates of the ADI values of the amplitude series  Az are completed. Thirty    30 K data are supposed to be available. The mean values (solid line) and the dispersion (error bars) of the estimated ˆA D are shown in Figure 1a, together with the values of the phase standard deviation (dotted line) V  . It can be noted that the ADI A D is a good measure of the phase stability for pixels with high SNR (signal-to-noise ratio). Those pixels exhibiting A D value under a given threshold T D (typically lower than 0.25) can be selected as PSs.
(a) (b)   varies from 0.05 rad to 0.33 rad (95% confidence interval). Since the ADI is a criterion in a statistical sense, some PSs with relatively large phase std   For rocky scenes with lots of deterministic point-wise scatterers, the 'true' PS number in the Set 1 S is enough for phase analysis. Therefore, based on the ADI criterion, another phase criterion could be further utilized together to reduce and refine the set 1 2 S +S to 1 S . However, for natural scenes, even if only the ADI criterion is utilized, the PS number in 1 2 S +S is still few. To increase the point number suitable for differential interferometry, those points with moderate amplitude stability and high phase stability, as shown with the yellow points in the set 3 S , could be utilized. This type of point is defined as QPS.  Figure 1b shows the scatter map of the dispersion index D A with the phase std σ V . For points located under the red line of D A = 0.25, their phase std σ V varies from 0.05 rad to 0.33 rad (95% confidence interval). Since the ADI is a criterion in a statistical sense, some PSs with relatively large phase std σ V > 0.25 are inevitably selected, as shown with the cyan points in the set S 2 . Those brown points in the set S 1 with D A ≤ 0.25 and σ V ≤ 0.25 are 'true' PSs which are stable in both amplitude and phase.
For rocky scenes with lots of deterministic point-wise scatterers, the 'true' PS number in the Set S 1 is enough for phase analysis. Therefore, based on the ADI criterion, another phase criterion could be further utilized together to reduce and refine the set S 1 +S 2 to S 1 . However, for natural scenes, even if only the ADI criterion is utilized, the PS number in S 1 +S 2 is still few. To increase the point number suitable for differential interferometry, those points with moderate amplitude stability and high phase stability, as shown with the yellow points in the set S 3 , could be utilized. This type of point is defined as QPS.

Quasi-Permanent Scatterer
QPS refers to a pixel with moderate amplitude stability and high phase stability, i.e., a moderate ADI D A and a small phase std σ V . The interferometric phase ϕ Int of a pixel in GB-SAR interferometry could be modeled as: where ϕ defo is the deformation phase component, ϕ atm is the phase component due to the atmospheric effects during image acquisitions, ϕ geom is a geometrical phase term due to In order to evaluate the phase stability of a pixel, the TPC γ TPC is defined [14]. N IN denotes the interferogram number, or the length of a pixel's phase sequence. The value range of γ TPC is from 0 to 1, and a larger value means a better stability of ϕ noi . A key problem is how to accurately estimate the noise term from the interferometric phase. The basic solution is to design a filtering framework based on the different temporal and spatial properties of those four phase terms in (2).
A further analysis of the numerical simulation in the Section 2.1 is made. Take interferometry process for each data sequence, and the interferometric phase ϕ Int of any two neighboring data can be regarded as the noise phase ϕ noi . Based on (3), the TPC γ TPC of each point is calculated. The mean values (solid line) and the dispersion (error bars) of γ TPC , together with the values of the phase std σ V (dotted line), are shown in Figure 2a. Contrast with Figure 1a, γ TPC decreases with the increasement of the noise std σ n . Figure 2b shows the scatter map of the TPC γ TPC with the phase std σ V . For those pixels with phase std σ V lower than 0.25 rad, i.e., left part of the green dotted line, 99.99% of them are with the TPC γ TPC larger than 0.91. To some extent, the TPC γ TPC is an equivalent criterion with the ADI D A . a moderate ADI A D and a small phase std V  . The interferometric phase Int  of a pixel in GB-SAR interferometry could be modeled as: where defo  is the deformation phase component, atm  is the phase component due to the atmospheric effects during image acquisitions, geom  is a geometrical phase term due to the aperture variation or repositioning error. noi  is the noise term. W[ ] denotes the phase wrapped function.
In order to evaluate the phase stability of a pixel, the TPC TPC  is defined [14].
IN N denotes the interferogram number, or the length of a pixel's phase sequence. The value range of TPC  is from 0 to 1, and a larger value means a better stability of noi  . A key problem is how to accurately estimate the noise term from the interferometric phase. The basic solution is to design a filtering framework based on the different temporal and spatial properties of those four phase terms in (2).
A further analysis of the numerical simulation in the Section 2.1 is made. Take interferometry process for each data sequence, and the interferometric phase Int  of any two neighboring data can be regarded as the noise phase noi  . Based on ( To better illuminate the relationship between the TPC TPC  and the ADI A D , Figure  3a  To better illuminate the relationship between the TPC γ TPC and the ADI D A , Figure 3a shows the scatter map of γ TPC with D A . By setting two proper thresholds of γ T and D T , those pixels with γ TPC ≥ γ T and D A > D T in the set S 3 could be selected as QPSs. For those pixels with the ADI D A lower than 0.25, and belonging to S 1 + S 2 , their TPC values are between 0.91 and 0.99 (90% confidence interval). The mean values (solid line) and the dispersion (error bars) of γ TPC with D A , are shown in Figure 3b. When D A equals to 0.25, the variation range of γ TPC is about from 0.91 to 0.96. Moreover, for those pixels with the TPC γ TPC larger than 0.91, and belonging to the set S 1 + S 3 , 99.99% of them are with the ADI D A lower than 0.45. Therefore, those pixels with the TPC γ TPC larger than 0.91, but with the ADI D A between 0.25 and 0.45, could be selected the QPSs. dispersion (error bars) of TPC γ with A D , are shown in Figure 3b. When A D equals to 0.25, the variation range of TPC γ is about from 0.91 to 0.96. Moreover, for those pixels with the TPC TPC γ larger than 0.91, and belonging to the set S1 + S3, 99.99% of them are with the ADI A D lower than 0.45. Therefore, those pixels with the TPC TPC γ larger than 0.91, but with the ADI A D between 0.25 and 0.45, could be selected the QPSs.

Distributed Scatterer
DSs refer to those scatterers which are affected by temporal, geometrical, and volumetric decorrelation, but their phase quality could be improved with their neighboring homogeneous pixels [15].
Ferretti et al. proposed an algorithm named by DespecKS to identify a pair of statistically homogeneous pixels [11]. To evaluate whether two pixels P 1 and P 2 are homogeneous, the two-sample KS test is utilized. By measuring the maximum difference K D between the cumulative probability distribution functions where K is the sequence length, and , ,..., K x x x 1 2 is the sorted amplitude sequence of a pixel. Figure 4a shows two curves whose amplitudes randomly vary and cannot determine their similarity. Figure 4b shows their sorted amplitude curves. With KS test, both curves satisfy the same normal distribution, and they can be regarded as a pair of homogeneous pixels.

Distributed Scatterer
DSs refer to those scatterers which are affected by temporal, geometrical, and volumetric decorrelation, but their phase quality could be improved with their neighboring homogeneous pixels [15].
Ferretti et al. proposed an algorithm named by DespecKS to identify a pair of statistically homogeneous pixels [11]. To evaluate whether two pixels P 1 and P 2 are homogeneous, the two-sample KS test is utilized. By measuring the maximum difference D K between the cumulative probability distribution functions S P 1 K (x) and S P 2 K (x) of their amplitude sequences, two pixels are homogeneous if D K is lower than a proper threshold.
where K is the sequence length, and x 1 , x 2 , . . . , x K is the sorted amplitude sequence of a pixel. Figure 4a shows two curves whose amplitudes randomly vary and cannot determine their similarity. Figure 4b shows their sorted amplitude curves. With KS test, both curves satisfy the same normal distribution, and they can be regarded as a pair of homogeneous pixels.  Take an image-pixel P as an example, and then define an estimation window centered on P. Apply two-sample KS test to count the number of its homogeneous pixels. Setting a proper number threshold, P can be regarded as a DS if its homogeneous pixels are enough. For a DS, its phase quality could be enhanced with its surrounding homogeneous pixels. It should be noted that only those DSs whose phase qualities are good enough could be utilized for phase analysis, and a further optimization is essential.

Pixel Selection
To select enough HQP for natural scenes, the key problem is how to accurately select QPS, and explore the feasibility of DS selection in GB-SAR interferometry. Take an image-pixel P as an example, and then define an estimation window centered on P. Apply two-sample KS test to count the number of its homogeneous pixels. Setting a proper number threshold, P can be regarded as a DS if its homogeneous pixels are enough. For a DS, its phase quality could be enhanced with its surrounding homogeneous pixels. It should be noted that only those DSs whose phase qualities are good enough could be utilized for phase analysis, and a further optimization is essential.

Pixel Selection
To select enough HQP for natural scenes, the key problem is how to accurately select QPS, and explore the feasibility of DS selection in GB-SAR interferometry.

PS Selection
Given N (no less than 20) GB-SAR images, two different ADI thresholds D T1 (0.1 to 0.25) and D T2 (0.25 to 0.45) are firstly utilized to take initial selection. With a small ADI threshold to select PS, the possibility of those PSs with large phase errors could be reduced effectively. Pixels with ADI value D A ≤ D T1 are set as the reference PSs. Pixels with ADI value D T1 < D A ≤ D T2 are set as the candidate set of the QPS (CS QPS ). Based on those reference PSs, the candidate set should be refined.

QPS Selection
With N GB-SAR images, every two consecutive images could be explored to construct an interferogram. Then the length of a pixel's interferometric phase sequence, or the interferogram number is N IN = N − 1. To calculate a pixel's TPC, the noise phase ϕ noi in (2) should be estimated. The first three terms ϕ defo , ϕ atm , and ϕ geom in (2) are highly spatially correlated, while ϕ noi is spatially uncorrelated. Therefore, a spatial filter could be utilized to filter out the spatial phase components and estimate ϕ noi .
The reference PSs distribute non-uniformly on the interferogram. They can be divided into a certain number of clusters according to their spatial locations with a K-means clustering algorithm, which is a method commonly utilized to automatically partition a point set into K clusters [16]. As shown in Figure 5, those scatter points are partitioned into 5 clusters which are marked with different colors. The sign 'x' denotes each cluster center. For those PSs belonging to a same cluster, calculate their mean location r and mean phase C lu ϕ . Since the interferometric phase Int ϕ of each pixel is wrapped, the circular period mean filter is utilized, as shown below.  For those PSs belonging to a same cluster, calculate their mean location r and mean phase ϕ Clu . Since the interferometric phase ϕ Int of each pixel is wrapped, the circular period mean filter is utilized, as shown below.
For a cluster, N Clu is the total pixel number belonging to it, and ϕ Clu(i) Int denotes the interferometric phase of its ith pixel. d Clu is the mean value of the complex phase of all the pixels belonging to this cluster. arg[ ] is the function to get the phase angle of a complex value.
For every pixel belonging to CS QPS , its spatial phase component ϕ Spa should be estimated. A 2-D spatial interpolation could be utilized, such as linear or Kriging interpolation. It should be noted that the complex phase exp(jϕ Clu ) rather than the wrapped phase ϕ Clu should be used for interpolation. Subtracting ϕ Spa from ϕ Int , the residual phase is ϕ noi [17].
This method is commonly utilized to compensate the space-variant phase in GB-SAR interferometry, if only those montionless PSs with no deformations are used for K-means cluster. Research has proved that the geometrical term ϕ geom caused by repeat-track error in GB-SAR interferometry could be similarly modeled with the atmospheric phase. Undoubtedly, ϕ geom can be well compensated with the interpolation method. Although the deformation term ϕ defo cannot be compensated thoroughly, its residual component during a short period for a slow-moving scene could be neglected.
Based on (3), the TPC of a pixel could be estimated with (6). ϕ (n) Int is the nth value of a pixel's interferometric phase sequence. ϕ (n) Spa is its estimated spatial phase component.
Setting a proper TPC threshold γ T , those pixels with γ TPC no smaller than γ T are selected as QPSs. Other pixels with γ TPC smaller than γ T are then selected to be the candidate set of the DS (CS DS ).

DS Selection
For a pixel belonging to CS DS , it is used as a window center, and then determined its homogeneous pixel according to whether their amplitude curves satisfy a same distribution. If the number of its homogeneous pixels is large enough, it can be selected as a DS. The key problem is how to improve the phase quality of these DSs, since their coherence is low and cannot be directly utilized for phase analysis.
The coherence matrix Γ(P) of a DS P could be estimated as, where Ω is the homogeneous pixel set of a DS P, and N HP denotes the number of its homogeneous pixels. d= exp(jϕ) is the complex vector of a pixel belonging to Ω, and ϕ = [ϕ 1 , ϕ 2 , · · · , ϕ N ] H is its phase vector. H denotes the complex conjugate. Based on the coherence matrix Γ(P), a method of maximum likelihood estimation could be utilized to extract the optimal phase λ = [ϑ 1 , ϑ 2 , · · · , ϑ N ] T of the DS P, as shown in (7). Λ is a column vector and equals to exp(jλ). · denotes the Hadamard product. A newton iteration method can be used to solve this optimization problem.
Once the optimal solution λ has been obtained, the quality of the estimated phase values ϑ i should be assessed. Another type of coherence γ DS could be utilized to be a 'goodness of fit' measure. Select every DS with γ DS higher than the TPC threshold γ T , and then substitute the original phase values with its optimized values.

HQP Selection
The processing flow developed for the HQP selection is shown in Figure 6, and it can be presented schematically as follows: (1) Two different ADI thresholds D T1 and D T2 are utilized for initial selection. Pixels with ADI value D A ≤ D T1 are set as the reference PSs. Pixels with D T1 < D A ≤ D T2 are set as CS QPS . (2) Based on the K-means clustering and spatial filtering, γ TPC of every pixel belonging to CS QPS is calculated. Setting a proper γ T , those pixels with γ TPC no smaller than γ T are selected as QPSs, and other pixels belong to CS DS . (3) For every pixel belonging to CS DS , its homogeneous pixels are determined. Define those pixels whose homogeneous pixels are enough as DSs. For each DS, phase optimization with its homogeneous pixels is made. Those DSs with coherence γ DS large than γ T are reserved. (4) Together with the PS, QPS, and DS as the HQP, the traditional PSInSAR algorithm could be utilized for differential interferometry.

Experimental Information
From 10:39 on 20 April 2019 to 8:08 on 30 April 2019 (China Standard Time), a GB-SAR system was utilized to monitor a vegetational scene (N30 • 50 39.83 , E108 • 25 3.58 ) located in Wanzhou district, Chongqing City, China. The plane form of this slope is a chair-like shape, as shown in Figure 7. The slope is about 80 m to 108 m long, 120 m to 250 m wide, with an average thickness of 13.33 m, and a volume of about 298,000 m 3 . It is a typical medium-scale middle-layer soil landslide. A village is located at the slope foot and on the slope. Those residential buildings are strong scatterers, i.e., PSs in the GB-SAR images. Figure 8 shows the photo of the GB-SAR. The instrument was developed by Beijing Institute of Technology, and provided by Suzhou Institute of Technology Reco Sensing Technology Co., Ltd. As a developed type of high-accuracy deformation measurement system, GB-SAR has been widely utilized to monitor surface deformations of natural or geological slopes.
SAR system was utilized to monitor a vegetational scene (°N30 50' 39.83'' , °E108 25'3.58'' ) located in Wanzhou district, Chongqing City, China. The plane form of this slope is a chair-like shape, as shown in Figure 7. The slope is about 80 m to 108 m long, 120 m to 250 m wide, with an average thickness of 13.33 m, and a volume of about 298,000 m 3 . It is a typical medium-scale middle-layer soil landslide. A village is located at the slope foot and on the slope. Those residential buildings are strong scatterers, i.e., PSs in the GB-SAR images.    Table 1 show the system parameters. The system works at Ku band with a wavelength of 1.85 cm. The transmitted signal is a frequency-modulated continuous wave. With the transmitting and receiving antennas moving along a high-precision mechanical track, a large aperture of 1.8 m could be synthesized. For a GB-SAR system with a fixedlength track, the Ku-band wave is conducive to a significant improvement in its azimuth resolution. In addition, it benefits increasing the measurement accuracy under the condition of a same level of noise phase, contrast with a lower-band wave. Although the Kuband electromagnetic wave cannot penetrate vegetations, those bare soils or ground surfaces make it feasible to select QPS and DS [8]. The experiment lasted for about 10 days, a total of 2370 images were acquired during the monitoring period. The average image acquisition time was about 6 min. Figure 9a shows a power image of the monitoring area at the polar coordinate, which was dB processed with respect to the average noise level, and could be regarded as a signal-to-noise ratio (SNR) map. Take the first = 30 N images as an example. Based on (1), the ADI value of every pixel could be calculated with these N amplitude images. Figure 9b shows the ADI map. Pixels with ADI values lower than 0.25 only occupy about 2.3%.  Table 1 show the system parameters. The system works at Ku band with a wavelength of 1.85 cm. The transmitted signal is a frequency-modulated continuous wave. With the transmitting and receiving antennas moving along a high-precision mechanical track, a large aperture of 1.8 m could be synthesized. For a GB-SAR system with a fixed-length track, the Ku-band wave is conducive to a significant improvement in its azimuth resolution. In addition, it benefits increasing the measurement accuracy under the condition of a same level of noise phase, contrast with a lower-band wave. Although the Ku-band electromagnetic wave cannot penetrate vegetations, those bare soils or ground surfaces make it feasible to select QPS and DS [8]. The experiment lasted for about 10 days, a total of 2370 images were acquired during the monitoring period. The average image acquisition time was about 6 min. Figure 9a shows a power image of the monitoring area at the polar coordinate, which was dB processed with respect to the average noise level, and could be regarded as a signal-tonoise ratio (SNR) map. Take the first N = 30 images as an example. Based on (1), the ADI value of every pixel could be calculated with these N amplitude images. Figure 9b shows the ADI map. Pixels with ADI values lower than 0.25 only occupy about 2.3%.
The experiment lasted for about 10 days, a total of 2370 images were acquired during the monitoring period. The average image acquisition time was about 6 min. Figure 9a shows a power image of the monitoring area at the polar coordinate, which was dB processed with respect to the average noise level, and could be regarded as a signal-to-noise ratio (SNR) map. Take the first = 30 N images as an example. Based on (1), the ADI value of every pixel could be calculated with these N amplitude images. Figure 9b shows the ADI map. Pixels with ADI values lower than 0.25 only occupy about 2.3%.  With the 1st image as the master image, and the 2nd and 30th images as the slave images, two interferograms IM 1-2 and IM 1-30 could be acquired. Figure 10a shows the phase interferogram IM 1-2 , whose temporal baseline is just several minutes. As a result of the atmospheric disturbance, for pixels located on the natural slope, their phases vary smoothly. The atmospheric refractivity affects the propagation speed of the radar wave through the troposphere. Different transmission delays could be inevitably caused due to the variant atmospheric condition, and it could even reach centimeter level at one kilometer. Atmospheric phase is one major error source in GB-SAR measurement, and must be well compensated [18]. Figure 10b shows the interferogram IM 1-30 , whose temporal baseline reaches about three hours. Its phase variation is much more noisy in contrast with that of IM 1-2 , since this vegetational scene is affected by the temporal decorrelation over a time period of several hours. With the 1st image as the master image, and the 2nd and 30th images as the slave images, two interferograms 1-2 IM and  IM could be acquired. Figure 10a shows the phase interferogram 1-2 IM , whose temporal baseline is just several minutes. As a result of the atmospheric disturbance, for pixels located on the natural slope, their phases vary smoothly. The atmospheric refractivity affects the propagation speed of the radar wave through the troposphere. Different transmission delays could be inevitably caused due to the variant atmospheric condition, and it could even reach centimeter level at one kilometer. Atmospheric phase is one major error source in GB-SAR measurement, and must be well compensated [18]. Figure 10b shows the interferogram  IM , whose temporal baseline reaches about three hours. Its phase variation is much more noisy in contrast with that of [1][2] IM , since this vegetational scene is affected by the temporal decorrelation over a time period of several hours.    Figure 11b shows the phase map of these pixels, which distribute on most parts of the whole map. A regular phase variation caused by the atmospheric disturbance can be observed. In addition, lots of noisy pixels exist, since the threshold T2 D is larger than the typical threshold T1 D utilized for PS selection. This phenomenon illuminates the feasibility of further selecting pixels with high phase quality.

Experimental Results
For the first 30 images acquired from 10:40 to 13:34 on 20 April 2019 (China Standard Time), by setting an ADI threshold of D T1 = 0.25, 25,575 reference PSs could be selected. Figure 11a shows the phase map of the reference PSs based on the interferogram IM 1-2 , which sparsely distribute on the whole map. By setting an ADI threshold of D T2 = 0.45, 159,753 pixels with ADI values D T1 < D A ≤ D T2 are selected, and they belong to the set CS QPS . Figure 11b shows the phase map of these pixels, which distribute on most parts of the whole map. A regular phase variation caused by the atmospheric disturbance can be observed. In addition, lots of noisy pixels exist, since the threshold D T2 is larger than the typical threshold D T1 utilized for PS selection. This phenomenon illuminates the feasibility of further selecting pixels with high phase quality.  Figure 11b shows the phase map of these pixels, which distribute on most parts of the whole map. A regular phase variation caused by the atmospheric disturbance can be observed. In addition, lots of noisy pixels exist, since the threshold T2 D is larger than the typical threshold T1 D utilized for PS selection. This phenomenon illuminates the feasibility of further selecting pixels with high phase quality. These reference PSs are divided into 70 clusters, as shown in Figure 12a. Those black squares denote the cluster centers, and every polygon marks the boundary of every cluster. With an inverse distance weighting interpolation, the spatial phase components of the pixels belonging to CS QPS are estimated, as shown in Figure 12b. It is worth noting that the estimated phases vary regularly both along the range and azimuth direction. For the phase components along the range direction, they are caused by the atmospheric disturbance at different times. While for the components along the azimuth direction, published researches prove that they are caused by the repeat-track error when radar moves along the mechanical track [19]. These reference PSs are divided into 70 clusters, as shown in Figure 12a. Those black squares denote the cluster centers, and every polygon marks the boundary of every cluster. With an inverse distance weighting interpolation, the spatial phase components of the pixels belonging to QPS CS are estimated, as shown in Figure 12b. It is worth noting that the estimated phases vary regularly both along the range and azimuth direction. For the phase components along the range direction, they are caused by the atmospheric disturbance at different times. While for the components along the azimuth direction, published researches prove that they are caused by the repeat-track error when radar moves along the mechanical track [19]. After compensating the spatial phase components, the TPC of those pixels belonging to QPS CS is calculated based on (6), as shown in Figure 13a. According to the analysis conclusion in Section 2.2, a TPC threshold of . T =0 91 γ equivalent to the ADI threshold of .  Figure 13c shows the probability density curves of TPC γ of the reference PSs and QPS CS . The shapes of both curves are much different. This illuminates that with a lower ADI threshold, pixels with high amplitude stability are more possible to be also with high phase stability. Figure 13d shows the scatter map of TPC . Although the phase quality of this type of PSs is doubtful, there is low necessity to further refine them for the natural scenes. The QPS belonging to the set 3 S makes up about 12.75% of the set QPS CS . Therefore, for pixels with high ADI values, it is rather hard to squeeze out enough pixels with high phase stability. After compensating the spatial phase components, the TPC of those pixels belonging to CS QPS is calculated based on (6), as shown in Figure 13a. According to the analysis conclusion in Section 2.2, a TPC threshold of γ T = 0.91 equivalent to the ADI threshold of D T = 0.25 is set. A total number of 20,364 QPSs can be selected, as shown in Figure 13b. Contrast with the PS number, the QPS number occupies about 79.6%. Pixels with the TPC γ TPC lower than the threshold γ T = 0.91 belong to the set CS DS . Figure 13c shows the probability density curves of γ TPC of the reference PSs and CS QPS . The shapes of both curves are much different. This illuminates that with a lower ADI threshold, pixels with high amplitude stability are more possible to be also with high phase stability. Figure 13d shows the scatter map of γ TPC with D A for the reference PSs and CS QPS . Its distribution trend is much similar with the simulation results shown in Figure 3a. The PS number belonging to the set S 2 , i.e., γ TPC < 0.91 and D A ≤ 0.25, is 2251 and makes up 8.8%. Although the phase quality of this type of PSs is doubtful, there is low necessity to further refine them for the natural scenes. The QPS belonging to the set S 3 makes up about 12.75% of the set CS QPS . Therefore, for pixels with high ADI values, it is rather hard to squeeze out enough pixels with high phase stability.
curves are much different. This illuminates that with a lower ADI threshold, pixels with high amplitude stability are more possible to be also with high phase stability. Figure 13d shows the scatter map of TPC  with A D for the reference PSs and QPS CS . Its distribution trend is much similar with the simulation results shown in Figure 3a. The PS number belonging to the set 2 S , i.e., . TPC 0 91   and A 0.25 D  , is 2251 and makes up 8.8%. Although the phase quality of this type of PSs is doubtful, there is low necessity to further refine them for the natural scenes. The QPS belonging to the set 3 S makes up about 12.75% of the set QPS CS . Therefore, for pixels with high ADI values, it is rather hard to squeeze out enough pixels with high phase stability. After the QPS selection, the remaining pixels are belonging to the set DS CS . For every pixel, it is used as the center of a 57  rectangular window, where 5 is the pixel number along the azimuth direction, and 7 is along the range number. The statistically homogeneous pixel of every pixel is determined based on the two-sample KS test. Figure 14a shows the homogeneous pixel map of DS CS . By setting a number threshold of T 10 N  , 111,406 DSs could be selected. Based on (7)-(9), DS  of these DSs are calculated and shown in Figure 14b. With the coherence threshold . T =0 91  , 4260 DSs are reserved. Although the DS number is only 16.7% of the PS number, it still benefits increasing the pixel number utilized for differential interferometry. Figure 14c shows the phase map of the selected DSs, which sparsely distribute on the map. Above all, a total of 50,199 HQPs are selected, including 25,575 PSs, 20,364 QPSs, and 4260 DSs. Figure 14d shows the phase map of the HQPs. Contrast with the PSs selected only with the ADI criterion, the HQPs S1 S3 S2 . After the QPS selection, the remaining pixels are belonging to the set CS DS . For every pixel, it is used as the center of a 5 × 7 rectangular window, where 5 is the pixel number along the azimuth direction, and 7 is along the range number. The statistically homogeneous pixel of every pixel is determined based on the two-sample KS test. Figure 14a shows the homogeneous pixel map of CS DS . By setting a number threshold of N T = 10, 111,406 DSs could be selected. Based on (7)-(9), γ DS of these DSs are calculated and shown in Figure 14b. With the coherence threshold γ T = 0.91, 4260 DSs are reserved. Although the DS number is only 16.7% of the PS number, it still benefits increasing the pixel number utilized for differential interferometry. Figure 14c shows the phase map of the selected DSs, which sparsely distribute on the map. Above all, a total of 50,199 HQPs are selected, including 25,575 PSs, 20,364 QPSs, and 4260 DSs. Figure 14d shows the phase map of the HQPs. Contrast with the PSs selected only with the ADI criterion, the HQPs increase by 96.3%. The spatial density of the pixels for deformation analysis has been significantly increased. Other than the increase in the HQP number, analysis on the phase quality should be further made. Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 18 could be formed by every two consecutive images. Based on the PS-InSAR algorithm, these interferograms are processed, including phase unwrapping, atmospheric compensation, and deformation measurement. Figure 15a shows the cumulative deformation map. Since these images were acquired during three hours, it was impossible for this slow-moving landslide to displace for several millimeters. Then the standard deviation of every HQP's deformation is calculated, and utilized to quantitively evaluate the phase quality of the PSs, QPSs, and DSs. Figure 15b shows the probability density curves of these three types of pixels. Most of the PSs are with deformation deviations around 0.3 mm, while for most of the QPSs and DSs, their deformation deviations are both around 0.42 mm. It is possibly caused by that a relatively low coherence threshold of .  For the first 30 images, 29 interferograms IM 1-2 , IM 2-3 , . . . , and IM 29-30 could be formed by every two consecutive images. Based on the PS-InSAR algorithm, these interferograms are processed, including phase unwrapping, atmospheric compensation, and deformation measurement. Figure 15a shows the cumulative deformation map. Since these images were acquired during three hours, it was impossible for this slow-moving landslide to displace for several millimeters. Then the standard deviation of every HQP's deformation is calculated, and utilized to quantitively evaluate the phase quality of the PSs, QPSs, and DSs. Figure 15b shows the probability density curves of these three types of pixels. Most of the PSs are with deformation deviations around 0.3 mm, while for most of the QPSs and DSs, their deformation deviations are both around 0.42 mm. It is possibly caused by that a relatively low coherence threshold of γ T = 0.91 is utilized, compared with the ADI threshold of D T1 = 0.25.
Every N = 30 images of all the 2370 images are divided into an image group. For each group, the HQP selection is taken with the same processing for the 1st group. Figure 16a shows the number curves of different types of HQPs for all the 79 groups. It should be noted that the number fluctuations of PS, QPS, and DS for time-series image groups are normal phenomena in GB-SAR applications, which is mainly caused by daily temperature variation, and could be solved by an adaptive selection threshold [7]. The PS numbers of different image groups change between 6355 and 46,920, and its average number is 17,702. In addition, the average numbers of QPS and DS are 10,052 and 3389, and their increase ratios contrast with the average PS number are 56.8% and 19.1%, respectively. Therefore, with the improved method to select the QPS and DS, the HQP number could be significantly increased. For the 12th image group, its HQP number is the most, and reaches 78,935, as shown in Figure 16b. The HQP density is much larger than that of the 1st group.  Figure  16a shows the number curves of different types of HQPs for all the 79 groups. It should be noted that the number fluctuations of PS, QPS, and DS for time-series image groups are normal phenomena in GB-SAR applications, which is mainly caused by daily temperature variation, and could be solved by an adaptive selection threshold [7]. The PS numbers of different image groups change between 6355 and 46,920, and its average number is 17,702. In addition, the average numbers of QPS and DS are 10,052 and 3389, and their increase ratios contrast with the average PS number are 56.8% and 19.1%, respectively. Therefore, with the improved method to select the QPS and DS, the HQP number could be significantly increased. For the 12th image group, its HQP number is the most, and reaches 78,935, as shown in Figure 16b. The HQP density is much larger than that of the 1st group. It is worth noting that the DS numbers in these groups are no more than 6958, and the average DS number is only 19.1% of PS. However, the procedures to take DS selection include homogeneous pixel calculation, phase optimization, and coherence thresholding, which need much computation. Although the DS technique in GB-SAR interferometry is feasible, the DS density is rather small. A lower coherence threshold could be set to select more DSs, but larger deformation deviations are unavoidable. Therefore, the necessity of DS selection in GB-SAR applications is low, unless when the PS number is rather few for a natural slope. Figure 17 shows the cumulative deformation map during the monitoring period of 10 days. Since the monitored slope is a creep landslide, no large deformation area is found  images of all the 2370 images are divided into an image group. For each group, the HQP selection is taken with the same processing for the 1st group. Figure  16a shows the number curves of different types of HQPs for all the 79 groups. It should be noted that the number fluctuations of PS, QPS, and DS for time-series image groups are normal phenomena in GB-SAR applications, which is mainly caused by daily temperature variation, and could be solved by an adaptive selection threshold [7]. The PS numbers of different image groups change between 6355 and 46,920, and its average number is 17,702. In addition, the average numbers of QPS and DS are 10,052 and 3389, and their increase ratios contrast with the average PS number are 56.8% and 19.1%, respectively. Therefore, with the improved method to select the QPS and DS, the HQP number could be significantly increased. For the 12th image group, its HQP number is the most, and reaches 78,935, as shown in Figure 16b. The HQP density is much larger than that of the 1st group. It is worth noting that the DS numbers in these groups are no more than 6958, and the average DS number is only 19.1% of PS. However, the procedures to take DS selection include homogeneous pixel calculation, phase optimization, and coherence thresholding, which need much computation. Although the DS technique in GB-SAR interferometry is feasible, the DS density is rather small. A lower coherence threshold could be set to select more DSs, but larger deformation deviations are unavoidable. Therefore, the necessity of DS selection in GB-SAR applications is low, unless when the PS number is rather few for a natural slope. Figure 17 shows the cumulative deformation map during the monitoring period of 10 days. Since the monitored slope is a creep landslide, no large deformation area is found It is worth noting that the DS numbers in these groups are no more than 6958, and the average DS number is only 19.1% of PS. However, the procedures to take DS selection include homogeneous pixel calculation, phase optimization, and coherence thresholding, which need much computation. Although the DS technique in GB-SAR interferometry is feasible, the DS density is rather small. A lower coherence threshold could be set to select more DSs, but larger deformation deviations are unavoidable. Therefore, the necessity of DS selection in GB-SAR applications is low, unless when the PS number is rather few for a natural slope. Figure 17 shows the cumulative deformation map during the monitoring period of 10 days. Since the monitored slope is a creep landslide, no large deformation area is found during a relatively short period, but some pixels with measurement errors of several mm are unavoidable. No other instruments with deformation measurement ability were available for the creep slope during the monitoring period. Only qualitative explanations can be made to illustrate the effectiveness of the improved method. Figure 17 shows the cumulative deformation map during the monitoring period of 10 days. Since the monitored slope is a creep landslide, no large deformation area is found during a relatively short period, but some pixels with measurement errors of several mm are unavoidable. No other instruments with deformation measurement ability were available for the creep slope during the monitoring period. Only qualitative explanations can be made to illustrate the effectiveness of the improved method.

Discussion
Pixels with high phase equality are essential and important for radar interferometry measurement. A common method is based on the ADI criterion to select PSs, as shown in Figure 1, which typically refer to deterministic point-wise scatterers, such as rocks. However, for natural scenes, the absence of point-wise PSs could severely affect interferometry measurement. Aiming at increasing the spatial density of HQP applied for natural scenes, an improved method to jointly select PS, QPS, and DS has been proposed.
With 30 GB-SAR images of a natural slope, 25,575 pixels with the ADI values A D lower than 0.25 are firstly selected as the PSs, and those pixels with < ≤ A 0.25 0.45 D are selected as QPS CS , as shown in Figure 11. Then based on these PSs, the K-means clustering and spatial interpolation are utilized to calculate the TPC γ TPC of every pixel belong-  Figure 13. Lastly, the two-sample KS test is utilized to calculate the homogeneous pixels of every pixel belonging to DS CS , and a number threshold is set to select DS. Based on (7) and (8) To validate the effectiveness of the proposed method, analyses are made from two aspects. As for the 30 images acquired during 3 h, the standard deviation of every HQP's deformation is calculated. As shown in Figure 15b, most PSs are with deformation deviations around 0.3 mm, while for QPSs and DSs, their deformation deviations are both around 0.42 mm. While slightly sacrificing phase quality, the HQP number has been significantly increased. As for the 2370 images acquired during 10 days, 79 image groups are taken HQP selection. Statistical results in Figure 16a show that the average numbers of PSs, QPSs, and DSs are 17,702, 10,052, and 3389, respectively.

Conclusions
This paper proposes an improved method to take HQP selection applied for natural scenes in GB-SAR interferometry. In order to increase the spatial density of HQP for phase measurement, three types of HQPs including PS, QPS, and DS, are selected with different criteria. The conventional ADI method which evaluates a pixel's amplitude stability is firstly utilized to take PS selection. In order to select those pixels with high phase stability

Discussion
Pixels with high phase equality are essential and important for radar interferometry measurement. A common method is based on the ADI criterion to select PSs, as shown in Figure 1, which typically refer to deterministic point-wise scatterers, such as rocks. However, for natural scenes, the absence of point-wise PSs could severely affect interferometry measurement. Aiming at increasing the spatial density of HQP applied for natural scenes, an improved method to jointly select PS, QPS, and DS has been proposed.
With 30 GB-SAR images of a natural slope, 25,575 pixels with the ADI values D A lower than 0.25 are firstly selected as the PSs, and those pixels with 0.25 < D A ≤ 0.45 are selected as CS QPS , as shown in Figure 11. Then based on these PSs, the K-means clustering and spatial interpolation are utilized to calculate the TPC γ TPC of every pixel belonging to CS QPS . A total of 20,364 pixels with γ TPC ≥ 0.91 are selected as QPSs, and those pixels with γ TPC < 0.91 are taken as CS DS , as shown in Figure 13. Lastly, the two-sample KS test is utilized to calculate the homogeneous pixels of every pixel belonging to CS DS , and a number threshold is set to select DS. Based on (7) and (8), γ DS of every DS is estimated. With a coherence threshold of γ T = 0.91, 4260 DSs are reserved.
To validate the effectiveness of the proposed method, analyses are made from two aspects. As for the 30 images acquired during 3 h, the standard deviation of every HQP's deformation is calculated. As shown in Figure 15b, most PSs are with deformation deviations around 0.3 mm, while for QPSs and DSs, their deformation deviations are both around 0.42 mm. While slightly sacrificing phase quality, the HQP number has been significantly increased. As for the 2370 images acquired during 10 days, 79 image groups are taken HQP selection. Statistical results in Figure 16a show that the average numbers of PSs, QPSs, and DSs are 17,702, 10,052, and 3389, respectively.

Conclusions
This paper proposes an improved method to take HQP selection applied for natural scenes in GB-SAR interferometry. In order to increase the spatial density of HQP for phase measurement, three types of HQPs including PS, QPS, and DS, are selected with different criteria. The conventional ADI method which evaluates a pixel's amplitude stability is firstly utilized to take PS selection. In order to select those pixels with high phase stability but moderate amplitude stability, a TPC criterion is adopted to take QPS selection. In addition, the feasibility of the DS technique in GB-SAR interferometry is explored. With a GB-SAR system, 2370 images of a natural slope were acquired during an experimental period of 10 days. As for the 30 images acquired during 3 h, statistical results of the deformation deviations prove that the HQP number can be significantly increased, while slightly sacrificing phase quality. As for the 2370 images acquired during 10 days, statistical results of the HQP numbers of 79 groups prove that, contrast with the average PS number, an average increase of 56.8% and 19.1% can be got with the QPS and DS selection.