A New Weighting Method by Considering the Physical Characteristics of Atmospheric Turbulence and Decorrelation Noise in SBAS-InSAR

: Time series of ground subsidence can not only be used to describe motion produced by various anthropocentric and natural process but also to better understand the processes and mechanisms of geohazards and to formulate e ﬀ ective protective measures. For high-accuracy measurement of small baseline subset interferometric synthetic aperture radar (SBAS-InSAR), atmospheric turbulence and decorrelation noise are regarded as random variables and cannot be accurately estimated by a deterministic model when large spatio-temporal variability presents itself. Various weighting methods have been proposed and improved continuously to reduce the e ﬀ ects of these two parts and provide uncertainty information of the estimated parameters, simultaneously. Network-based variance-covariance estimation (NVCE) and graph theory (GT) are the two main weighting methods which were developed on the basis of previous algorithms. However, the NVCE weighting method only focuses on the inﬂuence of atmospheric turbulence and neglects the decorrelation noise. The GT method weights each interferogram in a time series by using the Laplace transformation. Although simple to implement, it is not reasonable to have an equal weight for each pixel in the same interferogram. To avoid these limitations, this study presents a new weighting method by considering the physical characteristics of atmospheric turbulence and decorrelation noise in SBAS-InSAR images. The e ﬀ ectiveness of the proposed method was tested and validated by using a set of simulated experiments and a case study on a Hawaiian island. According to the GPS-derived displacements, the average RMSE of the results from the new weighting method was 1.66 cm, indicating about an 8% improvement compared with 1.79, 1.80 and 1.80 cm from the unweighted method, the NVCE method and the GT method, respectively.


Introduction
The small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) has been proven as a powerful technique for mapping precipitable water vapor (PWV) in time series [1][2][3][4][5], recovering digital elevation maps (DEM) [6] and measuring the ground surface deformation over large areas caused by crustal movement (e.g., earthquakes, volcanoes and landslides) and human activities (e.g., improve the measurement accuracy of SBAS-InSAR but also the assessment of the uncertainty of estimated parameters.
The rest of this paper is structured as follows: in Section 2, the weighting method considering atmospheric turbulence delays and decorrelation noises is proposed, after a brief introduction of SBAS-InSAR. Then, Section 3 presents simulated and real experiments wherein the proposed new weighting method was applied. Finally, some discussions and conclusions are drawn in Sections 4 and 5, respectively.

Review of SBAS-InSAR Technique
Consider that N SAR images (Figure 1a) are connected to M (where N 2 ≤ M ≤ N(N−1) 2 ) interferometric pairs (i.e., Figure 1b) and each SAR image contains T pixels. According to the principle of standard differential InSAR (D-InSAR), the differential phase (i.e., ∆ϕ h,i,j di f f in Equation (1)) of the hth pixel in an interferogram with master acquisition t i and slave acquisition t j can be written as [14]: where ∆ϕ h,i,j de f o is the surface deformation phase between the ith and the jth SAR acquisition; ∆ϕ h,i,j topo is the topographic phase, which can be compensated by an external digital elevation model (DEM). ∆ϕ h,i,j orb is the orbit error phase and can be separated from ∆ϕ h,i,j di f f by using a surface trend model [23]. ∆ϕ h,i,j atm and ∆ϕ h,i,j dec are the atmospheric and decorrelation noise phase, respectively. As stated above, ∆ϕ h,i,j atm includes the vertical and turbulence part (i.e., ∆ϕ h,i,j atm_vert and ∆ϕ h,i,j atm_turb , respectively) [14]. Therefore, after removing the vertical stratification atmosphere using a functional model, only the atmospheric turbulence and decorrelation errors are the main factors that affect the weighting matrix for each pixel and limit the measurement accuracy of geophysical parameters.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 21 improve the measurement accuracy of SBAS-InSAR but also the assessment of the uncertainty of estimated parameters. The rest of this paper is structured as follows: in Section 2, the weighting method considering atmospheric turbulence delays and decorrelation noises is proposed, after a brief introduction of SBAS-InSAR. Then, Section 3 presents simulated and real experiments wherein the proposed new weighting method was applied. Finally, some discussions and conclusions are drawn in Section 4 and Section 5, respectively.

Review of SBAS-InSAR Technique
Consider that SAR images (Figure 1a) are connected to (where 2 ≤ ≤ interferometric pairs (i.e., Figure 1b) and each SAR image contains pixels. According to the principle of standard differential InSAR (D-InSAR), the differential phase (i.e., ∆ ℎ, , in Equation (1)) of the ℎ ℎ pixel in an interferogram with master acquisition and slave acquisition can be written as [14]: where ∆ ℎ, , is the surface deformation phase between the ℎ and the ℎ SAR acquisition; ∆ ℎ, , is the topographic phase, which can be compensated by an external digital elevation model (DEM). ∆ ℎ, , is the orbit error phase and can be separated from ∆ ℎ, , by using a surface trend model [23]. ∆ ℎ, , and ∆ ℎ, , are the atmospheric and decorrelation noise phase, respectively. As stated above, ∆ ℎ, , includes the vertical and turbulence part (i.e., ∆ _ ℎ, , and ∆ _ ℎ, , , respectively) [14]. Therefore, after removing the vertical stratification atmosphere using a functional model, only the atmospheric turbulence and decorrelation errors are the main factors that affect the weighting matrix for each pixel and limit the measurement accuracy of geophysical parameters. interferograms. (c) Spatial structure of SBAS-InSAR. Each grid in (a) and (b) represents a pixel. The red pixel denotes the reference point in SBAS-InSAR. Black dots (e.g., ) and black edges (e.g., ) in (c) represent SAR images and interferograms formed by the corresponding SAR images. Bp and Bt are perpendicular and time baselines, respectively.
After removing the residual topographic phase, orbit error and vertical stratification atmosphere (i.e., ∆ ℎ, , , ∆ ℎ, , and ∆ _ ℎ, , , respectively), Equation (1) can be written in the form of matrices: where is the differential interferometric phase vector of pixel ℎ and the dimension of is × 1. is the × 1 unknown parameter vector (i.e., the deformation at each epoch); is the noise which contains the atmospheric turbulence and decorrelation noise. is an incidence matrix which has rows and columns. If the older SAR image be adopted as the master, matrix can be filled in each row with 1 for the former epoch, −1 for the latter and 0 for the rest.
where Y is the differential interferometric phase vector of pixel h and the dimension of Y is M × 1. X is the N × 1 unknown parameter vector (i.e., the deformation at each epoch); e is the noise which contains the atmospheric turbulence and decorrelation noise. G is an incidence matrix which has M rows and N columns. If the older SAR image be adopted as the master, matrix G can be filled in each row with 1 for the former epoch, −1 for the latter and 0 for the rest.

The Variance-Covariance Matrix of Atmospheric Phase in SBAS-InSAR
The vertical stratification atmosphere can be fitted and removed by a model related to the topography [14,23]; therefore, only the turbulence part is analyzed in this study. The atmospheric turbulence is mainly caused by the turbulent mixing of heat and humidity within the atmospheric boundary layer. According to Hanssen (2001) [14], the effect of the atmospheric turbulence is generally less than fifteen minutes with the horizontal scale being less than 1 km. Due to the large spatial and temporal variability, the atmospheric turbulence is hard to estimate by using a deterministic model. Fortunately, Hanssen pointed out that the atmospheric turbulence is relative to the distance between an arbitrary pixel to the reference one. Based on this physical characteristic, the atmospheric turbulence can be described by the variogram or structure function (i.e., Equation (4)) based on the intrinsic stationarity hypothesis [44].
where Var atm (r) is the atmospheric turbulence variance with a spatial distance of r, < · > is the statistical expectation operation and Z(·) denotes the observed value. Based on Equation (4), we can obtain a series of variances caused by atmospheric turbulence with different spatial distance r. To obtain the continuous result, namely, variances with any distance, the spherical model (i.e., Equation (5)) is commonly used to approximate the variogram of spatial turbulence delays [42].
where a is the largest spatial correlation distance and the so-called range. d 0 represents the nugget which mainly caused by the spatial variation of atmospheric turbulence in interferogram. d denotes the sill value of the model in Geostatistics and d 0 + d is commonly considered as the variance of interferogram.
To avoid the effects of deformation, the stacking technique is first used to estimate the subsidence, and the corresponding location of deformation in each interferogram is masked [45]. Then, based on the masked interferogram and Equation (4), the unknown parameters in Equation (5) (i.e., d 0 , d, a) can be estimated. According to Equation (5), the variance of the hth pixel with the reference point in lth interferogram can be obtained.
Based on the spatial graph (i.e., Figure 1c), the relationship of atmospheric turbulence between jth interferogram (i.e., the I j th edge in Figure 1c) and the corresponding SAR images (i.e., the s i th, s k th vertices in Figure 1c) is [43]: where TA is the atmospheric turbulence. It should be noted that the physical characteristic of atmospheric turbulence is uncorrelated in the time domain [14]. Therefore, combining Equation (6) and the covariance propagation formula, the relationship of the atmospheric variance between the M interferograms and N epochs in Figure 1c can be written as: Remote Sens. 2020, 12, 2557 Var atm (I 1 ) = Var atm (s 1 ) + Var atm (s 2 ) Var atm (I 2 ) = Var atm (s 1 ) + Var atm (s 3 ) . . .
where Var atm (I l ) is the variance of the lth interferogram; and l = 1, 2, · · · , M, Var atm s q and Var atm (s t ) are the atmospheric variances of corresponding SAR epoch which are connected the lth interferogram, q, t ∈ {1, 2, · · · , N}. Equation (7) can be written in the matrix form: where V atm (I) and V atm (s) are the atmospheric variances of the M interferograms with the dimension of M × 1 and N epochs with the dimension of N × 1 in Figure 1c, respectively. B denotes the design matrix with the dimension M × N, which is determined by the spatial structure of the graph; namely, Figure 1c and can be written as: where abs represents absolute value. In Equation (9), each row represents an interferogram. Meanwhile, B is filled with 0 except the locations of SAR epochs which are filled with 1. Combing Equations (8) and (9), the atmospheric variance of the selected pixel h in time series (i.e., V atm (s)) can be calculated by the weighted least-squares method.
where P denotes the weight of the atmospheric variances of the M interferograms (i.e., V atm (s) in Equation (10)), which is taken as the identity matrix for simplification in this study. After obtaining the atmospheric variance of the pixel h (i.e., V atm (s)), the next step is to calculate the atmospheric covariance matrix in time domain. Similar to Equation (2), the relationship of the atmospheric turbulence between the M interferograms and N epochs in Figure 1c (i.e., Equation (6)) can be written as: where TA I and TA s are the atmospheric turbulence of the M interferograms with a dimension of M × 1 and N epochs with a dimension of N × 1, respectively. Then, the atmospheric variance and covariance matrix of the pixel h in time domain can be obtained by the variance and covariance propagation law: where C atm and C atm (S) is the covariance matrix of the M interferograms with a dimension of M × M and the covariance matrix of N epochs with a dimension of N × N, respectively. G is an incidence matrix and it has been defined in Equation (3). T represents transposition. As the correlation of atmospheric turbulence in each SAR image closes to zero, C atm (S) can be described by a diagonal matrix whose diagonal elements are the corresponding variance of atmospheric turbulence.

The Variance-Covariance Matrix of Decorrelation Noise in SBAS-InSAR
The decorrelation noise in SBAS-InSAR can be divided into five parts [46]. First, the temporal decorrelation: the result of changes in backscattering characteristics with time in a resolution cell causes a correlation between two interferograms which have an overlap period of time. Second, baseline decorrelation, the result of varying incidence angles between different radar paths at the resolution cell, causes a correlation between two interferograms which have same incidence angle. The third is system noise decorrelation which depends on the signal-to-noise ratio, and the fourth is a processing-inducing noise decorrelation such as coregistration and interpolation noise in interferogram formation. Last is Doppler-centroid decorrelation, the result of varying viewing geometry in azimuth. Due to the fact that the time series of SAR images over the same area from the same track have the same imaging geometry, it can cause the correlation in decorrelation noise between any two interferograms, even if there is no overlap in temporal baseline or no shared image. Therefore, the covariance matrix caused by decorrelation noise of the M interferograms can be expressed as: where C dec is the total variance and covariance matrix caused by decorrelation noise; σ dec I i , I j is the covariance between the ith and jth interferograms calculated based on an analytical approximation using nonlinear error propagation (i.e., Equation (15)). If the ith and jth interferograms are formed by the s 1 th, s 2 th and s 3 th, s 4 th imaging time epochs, respectively, then where L is the number of independent looks, and γ s i s j is the coherence value between the s i th and s j th epoch [47].

The Weight of Each Pixel in SBAS-InSAR
Based on Equations (12) and (14), the total VCM of the selected pixel (i.e., C T ) is: The corresponding weight of the pixel (i.e., W T ) is: Then, the optimal linear unbiased estimator of the unknown parameters X and the covariance matrix C XX of Equation (3) can be obtained by: Remote Sens. 2020, 12, 2557 7 of 21 Since each pixel has the same spatial structure graph with the selected pixel (i.e., the pixel h), the weight of each pixel in SBAS-InSAR can be calculated from Equations (1)- (19) and we can obtain the deformation parameter.

Synthetic Test and Results
According to Equations (18) and (19), the weight for each pixel affects not only the accuracy of SBAS-InSAR measurement, but also the uncertainty of estimated parameters. In order to assess the performance of the new weighting method, a set of simulated experiments was conducted. At first, a total of 23 independent SAR images were simulated with a realistic spatial and temporal baseline distribution similar to the Sentinel-1 datasets over the Los Angeles Basin, South California. The spatial dimension of each simulated image was 100 × 100 pixels. As shown in Figure 2a, a subsidence funnel with 5 cm/a velocity was simulated. Subsequently, an atmospheric turbulence signal was simulated for each image based on the Kolmogorov turbulence theory [14]. The number of generated interferograms was 156 with temporal and spatial baseline thresholds of 342 days and 152 m, respectively. The spatial structure of each pixel in SBAS-InSAR is shown in Figure 2b; each vertex represents the simulated time epoch of image and each edge denotes the interferogram. In addition, five decorrelation factors (i.e., thermal decorrelation, Doppler-centroid decorrelation, coregistration induced decorrelation, geometric decorrelation and temporal decorrelation) were simulated for each interferogram. Under these conditions, the unweighted method, the GT method, the NVCE method and the new method were used to estimate the mean deformation velocity of each pixel, as shown in Figure 3.
Based on Figure 3a,d,g,j, we show that the four weighting methods have a good agreement with the simulation (i.e., Figure 2a). However, compared to Figure 3b,e,h,k, the estimation from unweighted method displays a relatively large deviation. On the contrary, due to considering the atmospheric delay phase and decorrelation noise in each interferogram, the results derived from the new method (i.e., Figure 3k) show the minimum deviation. To make it clearer, the corresponding histograms of the second column are shown in the third column of Figure 3. It can be observed that large standard deviations (std) of the unweighted and GT methods are significant.
To clearly show the superiority of the proposed method, Table 1 shows the quantitative comparison between the new and NVCE methods. It can be seen that the new method has an improvement of 4.98% in std. Meanwhile, the histograms of new method (i.e., Figure 3l) are closer to the Gaussian distribution in kurtosis and skewness than the NVCE method. These illustrate the necessity of considering the decorrelation noise in the spatial domain. the Gaussian distribution in kurtosis and skewness than the NVCE method. These illustrate the necessity of considering the decorrelation noise in the spatial domain.

Real Test Case Example: Big Island of Hawaii
For validation purposes, 24 descending Sentinel-1 SAR images over the Big Island of Hawaii were acquired from 5 January 2018 to 12 December 2018, as shown in Figure 4a. The basic parameters of the used Scentinel-1A SAR acquisitions are shown in Table 2. From the available images, 163 interferograms were generated under the conditions with a maximum temporal baseline of 145 days and a perpendicular baseline of 100 m (see Figure 4b). In addition, there are 28 Global Positioning

Real Test Case Example: Big Island of Hawaii
For validation purposes, 24 descending Sentinel-1 SAR images over the Big Island of Hawaii were acquired from 5 January 2018 to 12 December 2018, as shown in Figure 4a. The basic parameters of the used Scentinel-1A SAR acquisitions are shown in Table 2. From the available images, 163 interferograms were generated under the conditions with a maximum temporal baseline of 145 days and a perpendicular baseline of 100 m (see Figure 4b). In addition, there are 28 Global Positioning System (GPS) stations within this study area. Meanwhile, the MKEA site (i.e., the red star in Figure 4a) was used as the reference point for each interferogram, and the remaining 27 GPS stations (as shown in Appendix A Figure A1) were used for validation.  The external digital elevation model (30 m) from the Shuttle Radar Topography Mission (SRTM) and the precise orbit information were used to remove the residual topographic and orbital phase, respectively. Meanwhile, a second order polynomial related to the elevation was used to correct the vertical atmospheric phase [14,23]. Meanwhile, the corresponding 163 coherence maps could be obtained according to the standard differential InSAR (D-InSAR) technique. To ensure the highquality solution, the mean coherence can first be estimated from the 163 interferograms and pixels  The external digital elevation model (30 m) from the Shuttle Radar Topography Mission (SRTM) and the precise orbit information were used to remove the residual topographic and orbital phase, respectively. Meanwhile, a second order polynomial related to the elevation was used to correct the vertical atmospheric phase [14,23]. Meanwhile, the corresponding 163 coherence maps could be obtained according to the standard differential InSAR (D-InSAR) technique. To ensure the high-quality solution, the mean coherence can first be estimated from the 163 interferograms and pixels with the mean coherences less than 0.3 are masked out. Then, the new method is used to weight each pixel. To better understand, the GPS station of MOKP (i.e., blue dot in Figure 4a) is selected to show the VCM caused by atmospheric turbulence (i.e., Figure 5a), the VCM caused by decorrelation noise (i.e., Figure 5b), the total VCM (i.e., Figure 5c) and the corresponding weight (i.e., Figure 5d). Obviously, there is a remarkable similarity between Figure 5a,c. It indicates that the atmospheric phase of this pixel remains the main source of uncertainty for InSAR measurements.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 21 Obviously, there is a remarkable similarity between Figure 5a,c. It indicates that the atmospheric phase of this pixel remains the main source of uncertainty for InSAR measurements. Based on Equations (18) and (19), we can obtain the time series displacements of MOKP site from this new weighting method (i.e., Figure 6d). For comparison's sake, the results from the unweighted method, the GT method and the NVCE method are shown in Figure 6a-c, respectively. The error bars in Figure 6b-d represent the std of the estimated deformation derived from the corresponding variance matrix in Equation (17). Compared with the former three methods, the best agreement between the new method (i.e., Figure 6d) and the GPS result (i.e., gray dots). The RMSE of the new method is only 1.023 cm, while those of the former three methods amount to 1.477, 1.481 and 1.435 cm, respectively. The proposed new weighting method therefore represents an improvement of 29-31% over the other three weighting methods. Based on Equations (18) and (19), we can obtain the time series displacements of MOKP site from this new weighting method (i.e., Figure 6d). For comparison's sake, the results from the unweighted method, the GT method and the NVCE method are shown in Figure 6a-c, respectively. The error bars in Figure 6b-d represent the std of the estimated deformation derived from the corresponding variance matrix in Equation (17). Compared with the former three methods, the best agreement between the new method (i.e., Figure 6d) and the GPS result (i.e., gray dots). The RMSE of the new method is only 1.023 cm, while those of the former three methods amount to 1.477, 1.481 and 1.435 cm, respectively. The proposed new weighting method therefore represents an improvement of 29-31% over the other three weighting methods. By repeating the four weighting methods for each pixel in the study area, the mean deformation velocity in LOS direction can be obtained, as shown in Figure 7. It is obvious that the unweighted result (i.e., Figure 7a) shows relatively large differences compared with other three weighting methods (i.e., Figure 7b-d) and the differences are mainly located in the Kilauea caldera which erupted on 3 May 2018 and caused the serious decorrelation in this area [48]. Kilauea is one of the most active volcanoes on Earth and has been in a state of continuous eruption since 1983 at the vent of rift zone in the east of the volcano [49]. During eruption, large volumes of lava flow out from the volcano's summit. It causes significant subsidence of the ground surface in short time.   By repeating the four weighting methods for each pixel in the study area, the mean deformation velocity in LOS direction can be obtained, as shown in Figure 7. It is obvious that the unweighted result (i.e., Figure 7a) shows relatively large differences compared with other three weighting methods (i.e., Figure 7b-d) and the differences are mainly located in the Kilauea caldera which erupted on 3 May 2018 and caused the serious decorrelation in this area [48]. Kilauea is one of the most active volcanoes on Earth and has been in a state of continuous eruption since 1983 at the vent of rift zone in the east of the volcano [49]. During eruption, large volumes of lava flow out from the volcano's summit. It causes significant subsidence of the ground surface in short time. Furthermore, the magmatic activity resulted in an Mw 7.2 earthquake in the south flank of Kilauea, Hawaii, on 4 May [48]. By repeating the four weighting methods for each pixel in the study area, the mean deformation velocity in LOS direction can be obtained, as shown in Figure 7. It is obvious that the unweighted result (i.e., Figure 7a) shows relatively large differences compared with other three weighting methods (i.e., Figure 7b-d) and the differences are mainly located in the Kilauea caldera which erupted on 3 May 2018 and caused the serious decorrelation in this area [48]. Kilauea is one of the most active volcanoes on Earth and has been in a state of continuous eruption since 1983 at the vent of rift zone in the east of the volcano [49]. During eruption, large volumes of lava flow out from the volcano's summit. It causes significant subsidence of the ground surface in short time.   To better exhibit the local deformation and further validate the new algorithm clearly, the time series displacements of these remaining 26 GPS stations were calculated from the unweighted, the GT, the NVCE and the new method (as shown in Appendix A Figure A1). Similar to the figure of mean deformation velocity (i.e., Figure 7), most GPS stations kept stable in time series except the points near the Kilauea volcano (e.g., the NUPM site). The maximum subsidence occurs at the NUPM station and is approximately as high as 20 cm from May to September. The result illustrates the effect of the large eruption of Kilauea volcano during the three months of the SAR image acquisition.

The Necessity of Considering Decorrelation Noise
As shown in Figure 5, the contribution of decorrelation noise (i.e., Figure 5b) is less significant than that of atmospheric turbulence (i.e., Figure 5a) in the total VCM of MOPK (i.e., Figure 5c). However, the coefficient matrix (i.e., matrix G) is rank-deficient from Equation (3). Combining Equations (12) and (13) makes the VCM of atmospheric turbulence (i.e., matrix C atm ) also rank-deficient. If the VCM of decorrelation noise (i.e., matrix C dec ) is ignored, Equation (16) can be rewritten as: The corresponding weight of the pixel (i.e., matrix W T ) is: Although the generalized inverse method can be used to obtain the inverse of the rank-deficient matrix (i.e., matrix W T ), the generalized inverse of matrix C atm is not unique. In addition, Gui et al., has pointed the fact that the computational stability of the generalized inverse method is rather poor and a small disturbance of matrix C atm can lead to a large deviation in the estimated matrix W T [50]. This would cause some bias of the estimated deformation in SBAS-InSAR.
Fortunately, matrix C dec is of full rank from Equation (14). The matrix C T obtained from the new weighting method is the sum of matrix C atm and matrix C dec . Therefore, matrix C T is also of full rank and can be used to obtain the inverse matrix directly. It means considering the effect of decorrelation noise is not only necessary for a comprehensive analysis of noises in interferograms, but also avoids the inverse of the rank-deficient matrix. These help to improve the accuracy of the estimated deformation from SBAS-InSAR.

Average Performance
For the perspective of quantitative evaluation, we calculated the mean, std, root mean-square error (RMSE) and correlation (corr) between the simulated velocity and the estimated results, as shown in Figure 8. Compared with the unweighted, the GT and the NVCE methods, the new method had improvements of 32.21-44.35%, 4.26-28.82%, 9.52-26.52% and 0.42-5.28% in mean, std, RMSE and corr, respectively. These further demonstrate the new method has a better performance in the weighting of each pixel.
ote Sens. 2020, 12, x FOR PEER REVIEW 13 of Figure 8. Statistical quantitative comparison: the vertical movement estimations derived from the unweighted method (blue), the GT method (green), the NVCE method (yellow) and the new method (red) in mean, std, RMSE and corr, respectively.

. Validation of the Performances with GNSS Datasets
In order to intuitively analyze, regarded GPS-derived truth and RMSEs of each GPS point unde ese four weighting methods can be obtained. Figure 9 describes the comparison results from thes r weighting methods. It is obvious that the RMSE of the new method is smaller than the othe ree weighting method in most GPS stations. The average RMSEs of these four weighting metho e 1.79, 1.80, 1.80 and 1.66 cm, respectively. This also illustrates the improvement of the new ighting method. This may be because the new weighting method does not take into account th ect of covariance of atmospheric turbulence in spatial domain. As mentioned in previous study, regarded as a big challenge at present [41]. In general, the accuracy of the results without weighting was the lowest among the four method

Validation of the Performances with GNSS Datasets
In order to intuitively analyze, regarded GPS-derived truth and RMSEs of each GPS point under these four weighting methods can be obtained. Figure 9 describes the comparison results from these four weighting methods. It is obvious that the RMSE of the new method is smaller than the other three weighting method in most GPS stations. The average RMSEs of these four weighting method are 1.79, 1.80, 1.80 and 1.66 cm, respectively. This also illustrates the improvement of the new weighting method. This may be because the new weighting method does not take into account the effect of covariance of atmospheric turbulence in spatial domain. As mentioned in previous study, it is regarded as a big challenge at present [41]. igure 8. Statistical quantitative comparison: the vertical movement estimations derived from the nweighted method (blue), the GT method (green), the NVCE method (yellow) and the new method red) in mean, std, RMSE and corr, respectively.

alidation of the Performances with GNSS Datasets
n order to intuitively analyze, regarded GPS-derived truth and RMSEs of each GPS point u four weighting methods can be obtained. Figure 9 describes the comparison results from t eighting methods. It is obvious that the RMSE of the new method is smaller than the o weighting method in most GPS stations. The average RMSEs of these four weighting me .79, 1.80, 1.80 and 1.66 cm, respectively. This also illustrates the improvement of the ting method. This may be because the new weighting method does not take into accoun of covariance of atmospheric turbulence in spatial domain. As mentioned in previous stu arded as a big challenge at present [41]. n general, the accuracy of the results without weighting was the lowest among the four met . The GT method makes full use of the geometric relationship between the SAR images erferograms, and this algorithm is simple and easy to implement. However, this me ders that each pixel in the same interferogram is polluted by noise to the same degree. In general, the accuracy of the results without weighting was the lowest among the four methods tested. The GT method makes full use of the geometric relationship between the N SAR images and M interferograms, and this algorithm is simple and easy to implement. However, this method considers that each pixel in the same interferogram is polluted by noise to the same degree. This assumption makes all pixels in the same interferogram have an equal weight. Based on the GT method, the NVCE weighting method focuses on estimating the influence of atmospheric turbulence on the weight of each pixel. However, this method ignores the correlation of decorrelation noise in time domain. The proposed new weighting method avoids this limitation, and the accuracy of the results by using this new method is increased by 7.49%, 7.91% and 8.24%, respectively.
Although the performance of the new method outperforms than traditional methods and it can be applied in SBAS-InSAR without external data, the effect of atmospheric turbulence covariance in spatial domain is neglected. It is a big challenge at present as many previous studies pointed out. Therefore, we will focus on calculating the VCM of atmospheric turbulence in temporal and spatial domain by developing optimization algorithms in the future.

Conclusions
We propose a new weighting method which accounts for the atmospheric turbulence delay and inter-epoch decorrelation noise in SBAS-InSAR. Simulated and real experimental results demonstrate that the proposed method can effectively improve the accuracy of SBAS-InSAR measurement and the uncertainty information of the estimated parameters. Meanwhile, the new method relies entirely on the interferogram itself instead of external data to reduce the influence of atmospheric turbulence delay and decorrelation noises. Furthermore, the new method is not restricted by conditions such as the topography and climate of the study area. Therefore, it can be further implemented as a routine work and inserted in the current SBAS-InSAR algorithm.
In addition, the surface deformation of Hawaii caused by crustal movement (i.e., volcanoes and earthquakes) has been monitored by using the new weighting method and it is more consistent with GPS measurements than traditional methods. This illustrates the reliability of the algorithm proposed in this study; on the other hand, it can provide a certain reference value for Hawaiian volcanic activity. and InSAR-derived time series displacement (red triangles) from four different weighting methods at 26 GPS stations in Figure 8. Each row represents a GPS site and the four columns represent the results from the unweighted, the GT, the NVCE and the new method, respectively.