NL-MMSE: A Hybrid Phase Optimization Method in Multimaster Interferogram Stack for DS-InSAR Applications

When the distributed scatterer interferometric synthetic aperture radar (DS-InSAR) technology is used for surface deformation monitoring, the accuracy strongly depends on the quality of phase optimization. Especially, in the low-coherence region, how to conveniently and effectively improve the quality of phase optimization has been a hot and difficult research topic in recent years. This article proposes a hybrid phase optimization method for the DS-InSAR technology, which chooses either nonlocal (NL) or minimum mean square error (MMSE) method for each pixel according to the distribution of statistically homogeneous pixel in two windows of different sizes. This hybrid method (NL-MMSE) is not limited by the number of interferogram and is not influenced by multimaster images. The NL-MMSE method was applied to the deformation monitoring in the Jinsha River basin, Tibet, China, using 28 Sentinel-1A SAR images acquired between February and December 2020. Compared with the NL and MMSE methods, the NL-MMSE method provides better phase quality of the interferogram stack and is able to extract more temporal coherence points for subsequent deformation inversion. The deformation monitoring results showed there are three obvious large-scale landslide deformation and dozens of small-scale deformation in the study area. It is demonstrated that the NL-MMSE method based on DS-InSAR technology can accurately monitor the detailed deformation characteristics of the low-coherence surface, and can be applied and promoted as an effective means of identifying and monitoring geological hazards.


I. INTRODUCTION
T IME-SERIES interferometric synthetic aperture radar (In-SAR) technology has been increasingly used to obtain long-term and high-precision deformation monitoring by inverting multiple scenes SAR images [1]. In high-coherence regions that are dominated by buildings, bridges, and exposed rocks, the persistent scatterer InSAR (PSInSAR) [2] and Interferometric Point Target Analysis [3] work well. However, in low-coherence regions that are dominated by farmland, woodland, and mountainous areas, Small Baseline Subset [4] and SqueeSAR [5] methods are used. In recent years, the distributed scatterer InSAR (DS-InSAR) technique, developed on the basis of the SqueeSAR method, has been increasingly used in low-coherence regions [6], [7], [8]. In the low-coherence region, the feature backscatter is low and lacks high temporal coherence points. Therefore, the techniques based on single master images, such as PSInSAR, can hardly obtain detailed deformation information. In contrast, the DS-InSAR technology selects statistically homogeneous pixel (SHP) within a certain size window and uses it to optimize the phase of interferogram stack, to improve the temporal coherence of feature targets, and improve the point selection efficiency in low-coherence areas [9], [10].
In low-coherence areas, time-series SAR data covering identical areas are inevitably subject to interference phase noise in data processing. Interferometric phase noise is composed of a number of inherent factors that can be broadly classified as system noise, decorrelation noise, and inaccurate signal processing noise [11]. In fact, decorrelation noise, with spatial and time baselines, as well as volume decorrelation, has the greatest effect on the interference phase. The phase noise level due to temporal decoherence is one of the main variables affecting the interferogram and is an inherent characteristic present in the InSAR technique [12]. In addition, the noise in the interferogram stack increases the difficulty of phase unwrapping and affects the accuracy of the deformation calculation. Therefore, optimizing the phase of the interferogram stack is a very important step in the interferometric processing, which has received a lot of attention [13], [14], [15].
Numerous phase optimization methods based on interferogram stack have been developed [16], [17]. The simplest and original one is to apply multilook filtering to the interferograms one by one, which is a simple averaging of adjacent pixels This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in a window of a certain size [18]. Although it is simple and easy to implement, the multilook filtering method causes resolution reduction and phase edge deflection in the areas with high topographic relief and high heterogeneity. Later, the minimum mean square error (MMSE) method, nonlocal (NL) filtering, eigenvalue decomposition, and maximum likelihood estimation (MLE) were used for phase optimization [19], [20], [21], [22]. In the DS-InSAR technique, the phase optimization method is related to the number of the selected SHP [5]. The selection of SHP identifies the feature points with the same statistical behavior within a certain size window. SHP selection methods include nonparametric selection algorithms, like the two-sample Kolmogorov-Smirnov test and the Baumgartner-Wei-Schindler (BWS) test, and parametric selection algorithms like the likelihood ratio test and the fast SHP [23]. According to the SHP selected by each pixel of the interferogram, the interference phase is accurately optimized while maintaining the spatial resolution.
For DS-InSAR technology, there are many traditional and emerging phase optimization methods that exhibit great promise [24]. The NL method assigns weights by the patch similarity between the target pixel and its SHP, and executes a weighted average for each SHP to complete phase optimization. And, a rectangular window is used instead of individual pixels in the sample selection strategy. The eigenvalue decomposition method estimates the coherence matrix of N SAR images based on the SHP of the identified target pixel, and then performs eigenvalue decomposition on the coherence matrix and selects the eigenvector with the largest eigenvalue as the temporal optimization value [25]. In addition, considering the results of SHP selection, N SAR images are multiplied by conjugate or weighted conjugate to create a coherence matrix, and the phase is iteratively optimized by MLE theory to improve the interferogram coherence. Based on the above principles, several extension methods have been proposed recently, such as phase decomposition-PSInSAR [26], nonlinear optimization estimation-Stanford method for PS (NLE-StaMPS) [27], coherent scatterer-InSAR [28] and joint scatterer-InSAR [29]. However, when phase optimization is performed using the MLE principle, N SAR images have to deal with N(N-1)/2 scenes differential interference stacks, rather than simply N scenes, or any scene between N and N(N-1)/2. With the continuous accumulation of operational satellite data, the radar remote sensing technology has ushered in the era of SAR big data. The long time-series of massive SAR data increases the dimensionality of the coherence matrix and makes its phase optimization more difficult, which drastically reduces the computational efficiency of data processing. In contrast, the spatial domain phase optimization methods, such as MMSE and NL, which are not affected by the increased dimensionality of the coherence matrix, have their unique application advantages. In addition, within the interferogram domain, the NL method takes the similarity between windows as weights and applies the weighted average to estimate the optimized interferogram and denoising phase. Unlike the traditional spatial domain phase optimization method, the NL method efficiently extracts and optimizes the phase by taking the window as a unit, which in turn preserves the fringe texture information of the interferogram [30].
The mountainous regions are usually densely vegetated, resulting in low-coherence feature scatterers. When the DS-InSAR technology is applied to monitor ground deformation in such areas, the interferogram stack has obvious decoherence [31]. Furthermore, for N scenes SAR images, the decoherence between interferograms becomes more significant as the growth of the time baseline, except for that between the interferograms of one or two adjacent scenes that have good coherence. To monitor the deformation and disaster in such areas, we need efficient, effective, and flexible phase optimization methods. As the NL and MMSE methods [32] can be used in a multimaster interferogram stack with the image number between N and N(N-1)/2, here we propose a hybrid phase optimization method (referred to as the NL-MMSE method), combining the NL and MMSE methods, for monitoring low-coherence surface deformation. The surface of the Jinsha River Basin in Tibet, China, was selected to validate the proposed NL-MMSE method.
The overall framework of this article is as follows. Section II briefs some existing phase optimization methods. Section III describes the NL-MMSE method. Section IV applies the proposed phase optimization method to the deformation monitoring of the study area, using 28 Sentinel-1A SAR images. The results are compared with that of the MMSE method and the NL method. Finally, Section V concludes this article.

A. NL Method
The framework of the NL-filter is to assign weighting strategy based on window similarity [33]. The application of NL-filter to phase optimization in InSAR technology has been one of the classic and potent methods [11], [30]. It traverses the pixels in the interferogram for phase optimization and its universal formula is denoted asĨ whereĨ u is the optimized phase value of pixel u, and I t is the phase value of homogeneous pixel t in window range Ω. w(u, t) is the weighted value that takes into account the phase composition from the homogeneous pixel t. The weight w(u, t) is given by the similarity between windows C u and C t , which is called the window-based model. And w(u, t)satisfies the conditions 0 ≤ w(u, t) ≤ 1 and t w(u, t) = 1. Usually, w(u, t) is computed using an exponential function: where d(N u , N t ) represents the Euclidean metric, which is obtained by calculating the similarity between the windows; Z(u) represents the parameter obtained by normalizing the Euclidean metric; h is the variable coefficient that controls the range of the parameter. Their expressions are as follows: where a is the Gaussian kernel when computing the minimum norm; D A is the amplitude deviation index of the window. The value of d(N u , N t ) is positively correlated with the weight. h is usually given by considering the level of interferogram gradient change. In the areas with rough terrains, the phase gradient variation is generally significant, so a smaller h is assigned to preserve the nonnoise phase as much as possible; in the areas with a flat terrain, the phase gradient change is generally gentle, so a larger h is assigned to remove the noise phase as much as possible.

B. MMSE Method
For N scenes time-series SAR images, set φ z as the central pixel value of a window in one image, and φ x as the noisefree image pixel we wish to recover. v is the noise and obeys a distribution with average value of 1 and variance of σ 2 v . So we have [19]: Assume that φ x and v are two independent statistics andφ x is estimated using the prior mean of the local window. The linear combination of φ z andφ x iŝ Due to the lack of a priori mean,φ x is actually unknown. From (1), we getφ z = E[φ z ] =φ x . Therefore, in the application, the mean value of the local window is set asφ z . Parameters a and b satisfy a = 1 − b. Also, the parameter b is defined as where var(·) is the calculated variance. Combining (6), var(φ x ) is derived and given by In the above equation, the result of var(φ x ) may be negative due to sample insufficiency or inappropriate value of σ 2 v . In such case, var(φ x ) should be set as zero to ensure 0<b<1. Combining (7) and applying a = 1 − b,φ x can be rewritten aŝ whereφ x is the central pixel value after optimization. b is a weight coefficient between the original central pixel value φ z and the local window mean valueφ z . Therefore, this phase estimation algorithm is adaptive and can achieve optimal filtering results.

C. Screening Coherence Points
According to the results of phase optimization, we calculate the temporal coherence of the interferogram stacks. Using the differential interference phase values before and after optimization, the goodness-of-fit metric for a time dimension is expressed as [5] where φ nk is the original differential interferogram; ϑ n and ϑ k are the optimized and original differential interferometric phases, respectively. According to the actual number of interferograms, γ can be also written as whereφ nk is the optimized differential interferogram; M is the number of differential interferograms. For DS-InSAR, the goodness-of-fit γ can also be considered as the temporal coherence, which is used to quantify the permanence of the backscattering of the pixel in the time dimension [34]. Considering the actual situation of each study area, a reasonable temporal coherence threshold is set. The highly temporal coherent points with a suitable distribution density are selected as the monitoring points to invert the deformation.

D. Phase Quality Evaluation
To quantitatively assess the proposed method, the phase differences (PD) and the phase standard deviation (PSD) are used [35], [36]. PD is the mean value of the sum of phase differences for each pixel in the entire interferogram. First, the average phase difference (APD) between a certain pixel and the surrounding eight adjacent pixels is calculated [see (13)]; then, the PD is calculated by traversing the whole interferogram. Theoretically, PD should be approaching to zero. In practice, the interferogram often has phase noise, making the cumulative calculation result larger. The APD is calculated as follows: The PD of the whole interferogram is given as follows: PSD is a metric that reveals the smoothness of the interferogram. A sliding window is established with a pixel as the center and traverses the entire interferogram. By calculating the standard deviation of the phase in the sliding window, then the PSD of the whole interferogram can be derived as follows: where φ(i, j) is the phase of the pixel in the sliding window, φ(i, j) is the phase mean of the sliding window, and s is the number of pixels in the sliding window. Smaller values of PD and PSD indicate better quality of the interferogram stack.

III. NL-MMSE METHOD
In order to further improve the phase quality of the interferogram stack, we fully consider the respective technical characteristics and application advantages of the NL method and the MMSE method, and propose the NL-MMSE method. For N scenes time-series SAR images, M (N − 1 ≤ M ≤ N (N − 1)/2) scene interferograms are formed through registration, multilook operation, and differential interferometric data processing. Before phase optimization, we should select SHP according to the intensity of the SAR image. The selection of SHP as a key step before phase estimation has numerous selection algorithms. Considering that multilook of SAR images has been performed in this article, the nonparametric statistical test is more preferable. Therefore, we use the BWS test, which has relatively higher test efficacy, to select SHP.
The local window size for SHP selection in this article is set as 15 × 15. Fig. 1 abstractly represents the two results of SHP selection using the BWS test. The number of SHPs (SHP 15×15 ) selected in the 15 × 15 window ranges from 0 to 255. In general, the pixels with more than 20 SHPs are selected as DS candidates in order to preserve PS information [5], [28]. So we set a 5 × 5 window closer to SHP = 20 [6]. The position of a small 5 × 5 window in a 15 × 15 window is shown in Fig. 1. According to the similarity of adjacent objects, the closer the spatial distance to the central pixel, the greater the probability of becoming a homogeneous pixel. This article introduces the idea of weight where r and c are the numbers of rows and columns, respectively, and k is the window step size. The above equation satisfies w 0 r,c = 0 and w k r,c = 1. We can take the weight distribution result as the distribution ratio of SHP for different window sizes, as shown in Fig. 2.
As Fig. 2 shows, when the window size is 5 × 5, the distribution ratio of SHP is close to 0.5. Obviously, this is only a theoretical calculation value. In practice, to satisfy a distribution ratio of 0.5 for a 5 × 5 window, the total number of SHP should be not higher than 50. So we set the following.
1) When SHP 15×15 > 50, the central pixel is considered to be located in an absolutely homogeneous region with distributed characteristics [see Fig. 1(a)]. Therefore, on the basis of the selected SHP, the central pixel can be optimized by (1) using the NL method. 2) When 1 < SHP 15×15 ≤ 50, the general consensus is that the distributed characteristics of the central pixel are not robust enough. We calculate the ratio of the selected SHP in the small window to that in the large window, and denote it as K, and use the distribution ratio 0.5 as a measure of K. Then, the distributed characteristics of the central pixel are objectively evaluated as follows: a) If K < 0.5; as shown in Fig. 1(a), the distribution of SHP for the central pixel is scattered, so the NL method is used for optimization. b) If K ≥ 0.5, as shown in Fig. 1(b), the distribution of SHP for the central pixel is concentrated. According to the principle that the features of adjacent pixels are closer, we use the MMSE method to optimize the central pixel within the small window. The filtering strategy of the MMSE estimation method is as follows: whereφ center and φ center are the optimized and original phase values of the central pixel, respectively. φ mean is the mean of all pixels in the 5 × 5 window. b is a weighting coefficient ranging from 0 to 1, which is calculated by minimizing the mean square error of the central pixel and its neighborhood pixels (all pixels inside the 5 × 5 window). The detailed calculation principle of the weighting factor b is shown in (8) and (9). 3) When SHP 15×15 = 1, the central pixel is considered not to have distributed characteristics and it is regarded as a PS candidate. We do not perform phase optimization for it, but will decide whether it could be a temporal high coherence point for subsequent surface deformation inversion according to the amplitude departure index information. Fig. 3 shows the specific execution flow of our proposed NL-MMSE method. The processing flow of the NL-MMSE method shows that the distribution characteristics of the central pixel is evaluated depending on the number of SHP, which in turn leads to the selection of a more suitable phase optimization method. In the interferogram, only a small number of pixels are phase optimized using the MMSE method, which satisfy the two requirements of a relatively concentrated SHP distribution within the 5 × 5 window and 1 < SHP 15×15 ≤ 50. A large number of pixels are phase optimized using the NL method based on the SHP results. Choosing a better method (either NL or MMSE) based on the characteristics of all pixels is more comprehensive and scientific than using only one of them. Therefore, using the NL-MMSE method can further improve the phase quality of interferogram stack, theoretically.

IV. EXPERIMENT RESULTS AND DISCUSSION
To verify the improved effectiveness of the NL-MMSE method, a vegetated area with low coherence was chosen as an example for the experiment. The experiment consists of two parts. One compares the phase optimization quality of the NL-MMSE method with that of the MMSE method and the NL method. Another is to extract DS points according to the temporal coherence after the phase optimization is completed with the NL-MMSE method, and perform deformation detection and further discussion on the study area.

A. Study Area and Data Description
The study area is located in the upper basin of the Jinsha River, bordering Jomda County and Dege County, China (see Fig. 4). The area has a typical alpine valley landscape, featured by a narrow section of the Jinsha River passing through deep valleys [7]. The geological structure is dominated by schist, slate, and shale of the Upper Triassic Lanashan Formation and crystalline tuff and feldspathic quartz sandstone of the Upper Triassic Tumugou Formation. Influenced by long-term plate movement, high winds, and extreme rainfall, the rock fissures are obvious and the rock structure is comparatively fragmented. This region grows large areas of forest and distributes a small amount of farm land and sparse buildings. The region has the highland cold temperate semihumid climate, featured by large fluctuations precipitation during the year and distinct wet and dry rainy seasons. There are occasional heavy rainfalls in some parts of the region. All these factors contribute to frequent and widely distribute geological hazards, such as landslides, debris flows, and fissures.
As shown in Fig. 4(c), in the study area, most of the normalized difference vegetation index (NDVI) values are larger than 0.4, indicating a predominant vegetation cover. This low-coherent vegetation cover of the study area has obvious distributed characteristics, suitable for testing DS-InSAR methods in comprehensively monitoring and identifying surface deformation and hazards. Moreover, the topography of the study area usually brings noise. Effectively optimizing the differential interference phase can improve the accuracy for surface deformation monitoring and geological disasters identification. Taken together, it is relatively representative and challenging to experiment low-coherence vegetation surface as a scheme to validate the NL-MMSE method in this article.
The Sentinel-1A satellite used C-band, terrain observation progressive scans mode to acquire images from track 99. The acquired images have the pixel space of approximately 2 and 14 m in range and azimuth, respectively. We selected 28 scenes of ascending orbit Sentinel-1A images covering the study area for the experiment, spanning the period from February 7, 2020, to December 27, 2020. The detailed image time information is shown in Table Ⅰ. In this article, the image on February 7, 2020, is used as the master image to register the remaining images, and are processed according to the multilook ratio of 4:1. The image stack range and azimuth size after multilook is 1100 × 1000 pixels, corresponding to the actual area of about 10 × 14 km. In addition, the 30-m resolution SRTM DEM from National Aeronautics and Space Administration (NASA) was used to remove terrain phase and geocode.

B. Phase Optimization Results and Evaluation
Based on the multimaster principle, a total of 53 interferogram stacks were generated from the preprocessed SAR images, using the time baseline of 24 days and spatial baseline of 122 m. The distribution of interferogram stacks is shown in Fig. 5. The BWS test method was used to select SHP. MMSE, NL, and NL-MMSE were applied for phase optimization of the interferogram stack according to the SHP.
Two interference pairs were selected to compare the effects of the three phase optimization methods, which are 2020/02/07-2020/02/19 in winter and 2020/07/12-2020/07/24 in summer. The original and optimized interferograms and the corresponding enlargements of the local areas in the red boxes are shown in Fig. 6. In summer, the interferograms have a lot of noise, especially in the forest or steep terrain areas with weak backscatter, and the coherence is poor [see Fig. 6(g)]; in winter, the interferograms have less noise and the coherence is better [see Fig. 6(c)], because the grass and crops withered, and the vegetation is dominated by tree trunks. Compared with the original interferograms, the coherence of the interferograms after phase optimization by all the three methods is improved in both winter and summer, and the noise is significantly removed. Among the three methods, the MMSE method has relatively poor performance [see Fig. 6(d) and (h)], as a large area of speckle noise remain in some decoherent regions. The performance of the NL and NL-MMSE methods is similar and requires further quantitative analysis. But the comparison shows that optimizing the phase of the interferogram stack is necessary.
We used PD and PSD to quantitatively analyze the phase optimization effect of the interferogram stack. The results of PD and PSD of 53 interferograms were calculated (see Fig. 7). Fig. 7(a) and (b) shows the distribution of the PD and PSD values calculated from the MMSE, NL, and NL-MMSE optimized interferogram stacks, respectively. After the phase optimization, both PD and PSD values are reduced. The mean PD values of the 53 interferograms optimized by the MMSE, NL, and NL-MMSE methods are 0.78, 0.74, and 0.70, respectively, which are 55.4%, 57.7%, and 60.0% lower than the original value 1.75; the mean PSD values of the 53 interferograms optimized by the MMSE, NL, and NL-MMSE methods are 1.12, 0.92, and 0.87, respectively, which are 32.5%, 44.6%, and 47.6%, lower than the original value 1.66. The calculations were conducted using MATLAB R2019b software and an Intel(R) 2.20 GHz CPU. The MMSE method takes only 12 min as it does not consider the SHP when the interferogram stack is optimized. NL and NL-MMSE methods involve window-by-window sliding optimization of the center pixel according to the SHP results, so they take 389 and 382 min, respectively. Therefore, in order to pursue higher optimization quality, computational efficiency is inevitably sacrificed, but the advantage is that there is no need to worry about computer memory. Taken together, it can be concluded that the NL-MMSE method has better phase optimization capability and can provide differential interferometric phases of higher quality for subsequent data processing.
After the phase optimization is completed, the temporal coherence of the interferogram stack is calculated by (12). Considering the actual situation that the overall coherence distribution in this study area is relatively low, we set the temporal coherence threshold to 0.5 and perform the subsequent deformation inversion. The results are as follows. The MMSE method selects 565 800 points; the NL method selects 377 545 points; the NL-MMSE method selects 399 566 points, an increase of 22 021 points compared to the NL method, which improves the distribution density of monitoring points. Fig. 8 shows the characteristics of the temporal coherence distribution, the histogram statistics, the distribution of the selected monitoring points, and the local amplification results of the monitoring points for the MMSE, NL, and NL-MMSE methods respectively. The MMSE results [see Fig. 8(a)] exhibit fake high temporal coherence and increased estimation bias, as the optimization does not consider the backscattering, leading to contamination of surrounding targets by points with high coherence [37]. The MMSE method even selects monitoring points from some areas of severe geometric distortion, indicating that the selection results are not valuable for research in low-coherence areas. The NL results based on SHP optimization [see Fig. 8(b)] has better spatial distribution characteristics. The NL-MMSE results [see Fig. 8(c)] increases the temporal coherence slightly and keeps the spatial distribution characteristics. The monitoring points selected by the NL-MMSE method are more scientific and rational than the MMSE results. The density of selected monitoring points is higher compared to the NL results.

C. Deformation Detection and Analysis
Based on the NL-MMSE phase optimization method and coherence point selection method proposed above, the line-ofsight (LOS) surface deformation monitoring results of the study area were obtained by inversion after the atmospheric phase separation and phase unwrapping of the coherence points, as shown in Fig. 9(a). The coherence points are densely distributed in the study area, except for some areas with geometric distortions caused by large surface undulations. Such distribution can reflect the time-series deformation pattern of the surface. During the study period, dozens of small-scale deformation sporadically distributed and three large-scale landslides with the deformation rates of around -100 mm/yr were observed in the study area, as shown in A, B, and C in the Fig. 9(a). As no external measurement data, such as leveling and GNSS, are available to verify the accuracy, we use the standard deviation of the residual deformation due to atmospheric delay and baseline error for internal verification. Fig. 9(b) shows the histogram of the standard deviation of the residual deformation, with the maximum value of 9.4 mm and the mean value of 1.7 mm, indicating that the monitoring results have a high reliability.
In this study area, as the small-scale deformation surfaces are large in number and sporadic in distribution, we do not provide detailed analysis and discussion. But these deformation surfaces can be used as identified geological hazard potential points to provide reference for on-site inspection and emergency monitoring by local geological hazard prevention and control departments. Three larger landslide deformations A, B, and C, respectively, are shown in detail in Fig. 10, and the legend for the LOS deformation monitoring results remains consistent with Fig. 9(a). A and B are located in the villages of Woda and Gaina in Jomda County, respectively, adjacent to the Jinsha River, where the terrain is steep and the slope has slow sliding for many years; C is located in the village of Aiguo in Jomda County, where the entire slope of the mountain has undergone significant sliding. As shown in Fig. 10(b) and (d), the deformation in A and B mainly distributes on the bare surface of the villages, and is pronounced in the middle of the mountain. Once the deformation accelerates, it is very likely to develop into a catastrophic disaster event that may cause damage to local people and block the Jinsha River. As shown in Fig. 10(f), the deformation in C is most extensive and the whole slope has different sliding speeds in different parts. The sliding area covers the village as well as the surrounding woodland and grassland.
In order to understand the variation features of the surface deformation at A, B, and C, we selected three feature points [see Fig. 10(b), (d), and (f) for specific geographical locations] from each of them located at the upper, middle, and lower parts. The deformation trends of the three landslides are analyzed by showing the time-series deformation of the feature points. Fig. 11(a)-(c) shows that all the feature points have time-series deformation of different magnitude, and they all show surface subsidence. Fig. 11(a) shows that A1 on the upper slope and A3 on the lower slope have close temporal deformation trends, with cumulative deformation values of 71 and 68 mm respectively, and A2 on the middle slope has the most noticeable temporal subsidence trend, with the cumulative deformation of 90 mm. Fig. 11(b) shows that B2 on the middle slope has the largest magnitude of temporal deformation with a cumulative deformation of 73 mm, and B1 on the upper slope and B3 on the lower slope have cumulative deformation of 47 and 60 mm, respectively. Fig. 11(c) shows that C3 on the lower slope has the least obvious temporal deformation trend, with a cumulative deformation of only 31 mm, and C1 on the upper slope has slow temporal deformation until June 2020, which may be related to seasonal permafrost, and after June it starts to deform faster, but still less than C2 in the middle of the slope. The cumulative deformation of C1 and C2 are 52 and 82 mm, respectively, which may be due to its denser vegetation generating greater slip resistance. For the three landslides, we conclude that A2, B2, and C2, located in the middle of the slope, have the largest time-series deformation. Combined with the existing research results of A and B [7], it is inferred that the deformation is most obvious in the central area of each landslide. In addition, we related the average deformation rate of A1-C3 to the local weather data and the results are shown in Fig. 11(d). From February to April, the weather is dominated by snowfall, and the average deformation rate does not change much, indicating that snowfall has little effect on the surface deformation. The average deformation rate fluctuates and increases from April to October as the rainfall increases. Between August and October, which is the rainy season, the average deformation rate reaches the peak. From October to December, the weather is again dominated by snowfall, and the average deformation rate decreased accordingly. Combining the results of existing analyses of the relationship between rainfall and surface deformation in the region [7], we conclude that rainfall is one of the main influencing factors to accelerate landslide deformation.

V. CONCLUSION
In this article, we propose the NL-MMSE phase optimization method where a better method (either NL or MMSE) is chosen depending on the characteristics of each pixel. Using 28 Sentinel-1A images covering the Jinsha River basin, we generated 53 interferograms to verify and analyze the superiority of the NL-MMSE method. Some conclusions can be drawn as follows.
1) By analyzing the phase optimization results of the 53 interferograms, the mean PD value of the NL-MMSE method is 0.08 and 0.04 smaller than that of MMSE and NL, respectively, and the mean PSD value is 0.25 and 0.05 smaller than that of MMSE and NL, respectively. Using 0.5 as the temporal coherence threshold, the NL-MMSE method selects 22 021 points more than the NL method. This demonstrates that the NL-MMSE method has a better capability of phase optimization and provides a higher density of monitoring points in the low-coherence areas covered by vegetation.
2) The monitoring results show that there are three obvious large-scale landslide deformations in the study area. Detailed deformation analysis of the three extensive landslide areas was carried out by considering the topographic and weather data. Subsequently, more deformation and deformation characteristics of low-coherence areas dominated by vegetation cover will be revealed in detail by the method of this article. His research interests include InSAR, surface deformation monitoring, and geological hazard susceptibility mapping.