Statistical Assessments of InSAR Tropospheric Corrections: Applicability and Limitations of Weather Model Products and Spatiotemporal Filtering

Abstract

Tropospheric correction is a crucial step for interferometric synthetic aperture radar (InSAR) monitoring of small deformation magnitude. However, most of the corrections are implemented without a rigorous evaluation of their influences on InSAR measurements. In this paper, we present three statistical metrics to evaluate the correction performance. Firstly, we propose a time series decomposition method to estimate the tropospheric noise and mitigate the bias caused by ground displacement. On this basis, we calculate the root-mean-square values of tropospheric noise to assess the general performance of tropospheric corrections. Then, we propose the use of semi-variograms with model-fitted range and sill to investigate the reduction of distance-dependent signals, and Spearman’s rank correlation between phase and elevation to evaluate the mitigation of topography-correlated signals in hilly areas. The applicability and limitations were assessed on the weather model-derived corrections, a representative spatiotemporal filtering method, and the integration of the two mainstream methods. Furthermore, we notice that the persistent scatter InSAR processing resulted in two components, the primary and secondary images’ contribution to the tropospheric and orbit errors. To the best of our knowledge, this paper for the first time analyzes the respective roles of the two components in the InSAR tropospheric corrections.

1. Introduction

A number of time series interferometric synthetic aperture radar (InSAR) techniques have been developed to measure slow ground motion, such as persistent scatter (PS) InSAR [1,2], small baseline subset (SBAS) InSAR [3], and SqueeSAR [4]. However, in many situations, the ability of InSAR to retrieve small deformation magnitude is severely limited by the tropospheric phase delays. The tropospheric phase is a major noise source in repeat pass InSAR measurements, which results from the spatial and temporal variations of pressure, temperature, and relative humidity in the troposphere. It can be in the order of several centimeters and contaminate small deformation signals, such as land subsidence [5], inter-seismic deformation [6], volcanic deformation [7], etc. Thus, tropospheric correction is a vital step to ensure accurate InSAR monitoring of slow ground motion.

The tropospheric phase consists of turbulent and stratified delays [8]. Various methods have been implemented to correct tropospheric phase delays in InSAR measurements and can be categorized into three groups. The first group is statistical approaches assuming the random distribution of the tropospheric phase in the time dimension, such as the pair-wise logic approach [9], the stacking approach [10,11], and the commonly used spatiotemporal filtering methods [1,2]. The second group is the empirical methods aimed to estimate the stratified tropospheric phase, based on the relationships between the stratified phase and topography, such as the linear model [12,13], the power-law model [14], and the quadtree-aided joint model [15]. The third group utilizes external data to estimate both the turbulent and stratified components. These auxiliary data include global positioning system (GPS) measurements of zenith total delay (ZTD), precipitable water vapor (PWV) measured by spectrometers such as the Medium Resolution Imaging Spectrometer (MERIS) and the Moderate Resolution Imaging Spectroradiometer (MODIS), as well as weather model products such as the ERA-Interim, the ERA5 reanalysis [16,17], and the Generic Atmospheric Correction Online Service (GACOS) [18,19,20].

Despite the numerous studies on tropospheric correction methods, the evaluation of their performance remains a hot research topic. Although a small number of studies verified the accuracy of tropospheric correction using in situ data, such as the GPS measured ZTD [21] or leveling data [22], the data availability is often limited due to factors such as the lack of GPS stations or sparse coverage of ground-based measurements. A vast majority of validations employed the spatial standard deviation (StdDev) or root-mean-square (RMS) values of individual interferograms [8,23,24,25], assuming that there is no surface deformation in the study area. However, in many cases, ground deformation occurs within the region of interest and may even encompass a considerable area. Basic statistics such as StdDev or RMS, obtained under the none-deformation assumption, could potentially exhibit bias. Nonetheless, some effects of the physical characteristics of the troposphere on InSAR measurements, such as the distance-dependence and topography-dependence [26], may not be revealed by the StdDev or RMS. Some recent studies have qualitatively examined the topography-dependence using phase-elevation scatterplots [22] or have expanded the analysis using the linear correlation of phase and elevation [25]; however, they were not be able to capture the full complexity of the topography-dependence. Murray et al., (2019) calculated the semi-variance of semi-variograms to analyze the distance-dependence [24]. These studies have contributed valuable insights into the performance of tropospheric correction. However, their evaluation is limited to individual mainstream approaches.

The performance of combined approaches, i.e., integration of the abovementioned different mainstream methods, has not been systematically evaluated, despite their increasing use in practice. In particular, weather model products suffer from coarse spatial and temporal resolution, resulting in limited performance in InSAR tropospheric corrections. Statistical approaches assume the temporally uncorrelated characteristics of tropospheric signals, which does not hold for the stratified tropospheric delays. In recent years, the spatial and temporal resolutions of weather model data have been significantly improved. To what extent can these features improve the accuracy of InSAR tropospheric corrections, especially in coastal regions with large spatiotemporal variations of water vapor? What are the pros and cons of the weather model-derived corrections, the statistical methods, and the incorporation of the two mainstream techniques? These are the questions this study also intends to answer.

Furthermore, little attention has been paid to the different roles of the respective contribution from the primary (reference) image and secondary images in the InSAR tropospheric corrections. While the tropospheric phase is typically assumed to be uncorrelated in time in StaMPS (Stanford Method for Persistent Scatterers) spatiotemporal filtering, one of the most commonly used and representative statistical methods, the use of a common reference image (i.e., primary image), introduces a temporally correlated component referred to as the primary image’s contribution to the tropospheric and orbit error (TOE) [2]. The contribution from secondary images (hereinafter referred to as secondary TOE) remains uncorrelated in time. The primary and secondary TOE exist in the interferograms corrected by spatiotemporal filtering or its integration with weather model products. The contribution of the primary and secondary TOEs in InSAR tropospheric corrections is still unclear and requires further investigation.

In this study, firstly, we propose a time series decomposition method to estimate the overall noise after InSAR tropospheric corrections, and mitigate the bias resulting from ground deformation. On this basis, the time series RMS of the estimated tropospheric noise is calculated as a metric to evaluate the overall performance of the InSAR tropospheric correction. This method can be easily applied to regions prone to ground displacement without relying on the none-deformation assumption. Based on the impacts of the physical properties of the troposphere on the large-area InSAR monitoring, we adapt two other metrics from recent works, i.e., the semi-variogram of interferograms and the phase-topography correlation [24,25]. The semi-variogram is employed to assess the reduction of spatially dependent signals after tropospheric corrections. The phase-topography correlation is used to assess the reduction of elevation-dependent signals, where we propose the use of Spearman’s rank correlation coefficient as a more general expression of the phase-elevation relationship strength.

This study employs advanced weather model products that are widely used in the InSAR community, including ERA5, ERA-Interim, and GACOS. The pros and cons of the sole use of weather model-derived corrections, the StaMPS spatiotemporal filtering (a representative statistical method), or the combined approaches are assessed using the proposed statistical metrics. In addition, we investigate the respective roles of primary and secondary TOE in InSAR tropospheric corrections, in terms of the decrease of overall noise, reduction of distance-dependent signals, and suppression of topography-dependent signals.

2. Materials

2.1. Study Area

The experiments were conducted on the east bank of the Pearl River Estuary (Figure 1), China, covering most areas of the city of Dongguan, Shenzhen, and the entire Hong Kong Special Administrative Region.

This area is featured by a typical subtropical ocean monsoon climate, with long flood seasons from April to September. The elevation ranges from −79 m to 1260 m. The intricately distributed hills, mountains, and plains have brought complex natural geographical conditions. The study area still has a few mountainous locations, the majority of which are over 500 m and even reach 1260 m above sea level.

In the process of rapid urbanization, there have been reported subsidence cases due to land reclamation and subway construction [27,28]. In recent decades, extreme weather events have occurred frequently in the study area. These characteristics pose a great challenge to InSAR tropospheric correction.

2.2. Data

2.2.1. Sentinel-1 SAR Imagery

The Sentinel-1 (S1) Interferometric Wide (IW) Terrain Observation by Progressive Scans (TOPS) Single Look Complex (SLC) data are employed in this research (

Table 1. Parameters of Sentinel-1 (Interferometric Wide) IW Single Look Complex (SLC) data.

Sentinel-1 IW SLC Data
Timespan	15 June 2015~7 March 2018
Revisit cycle (days)	12
Polarization	VV
Incidence angle (°)	41.69~46.11
Wavelength (cm)	5.5
Slant range spacing (m)	2.33
Azimuth spacing (m)	13.92

), consisting of 68 images acquired from June 2015 to March 2018, mostly with a 12-day interval. We selected the image acquired on 12 March 2017 as the primary (reference) image, and the rest as the secondary images. The temporal baseline varies from 12 to 372 days, and the normal baseline ranges from 3.7 to 82.3 m.

2.2.2. Digital Elevation Model (DEM)

The Advance Land Observing Satellite (ALOS) World 3D-30 m DEM, abbreviated as AW3D30, is used to remove the topographic phase from interferograms. The AW3D30 DEM was generated using the archived data of the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) onboard ALOS during its operation from 2006 to 2011.

2.2.3. Weather Model Products

In this study, three weather model products are used for tropospheric corrections, including the ERA-Interim, ERA5 reanalysis data, and GACOS.

ERA-Interim (ERA-I) is a global atmospheric reanalysis product released by ECMWF. It covers the period from 1 January 1979 to 31 August 2019, and contains 60 pressure levels to 0.1 hPa (about 64 km altitude), with a temporal resolution of 6 h and a spatial resolution of 0.7° × 0.7° (~78 km).

ERA5 is the fifth generation of atmospheric reanalysis produced by ECMWF. It is available from 1950 to the present. As the successor to ERA-I, ERA5 provides hourly estimates of the global atmosphere, land surface, and ocean waves at a spatial resolution of 0.25° × 0.25° (~28 km). It resolves tropospheric parameters for 137 vertical levels up to 0.01 hPa.

GACOS employs the ECMWF HRES data as input. It utilizes the Iterative Tropospheric Decomposition (ITD) model [19] to decouple stratified and turbulent delays, and generates ZTD maps for InSAR tropospheric correction. The ECMWF HRES is an operational climate analysis model, estimating atmospheric variables with a temporal resolution of 6 h, and a horizontal resolution of 0.125° × 0.125° (~14 km), for 137 pressure levels up to 0.01 hPa.

3. Methods

We propose the use of three statistical metrics to assess the performance of InSAR tropospheric corrections. All three metrics are applied on time series interferograms processed by the PS InSAR approach, after the steps of co-registration, removal of the topographic phase, PS candidate selection, correction of spatially uncorrelated phase, and phase unwrapping.

3.1. Evaluation Metrics for InSAR Tropospheric Correction

3.1.1. Tropospheric Noise Estimated by Time Series Decomposition

In previous studies, the root-mean-square (RMS) value or standard deviation (StdDev) on individual interferograms is a widely used metric to assess the residual noise after InSAR tropospheric corrections, assuming no surface deformation occurs [8,24]. In this research, instead of the spatial statistics of interferograms, we propose to analyze the time series noise of interferograms, which reveals the variance of the residual signals relative to the rest of the interferograms. It better captures the tropospheric variability before and after correction than a simple average of the spatial variance of individual interferograms. However, there have been reported subsidence cases in the study area. Instead of a simple calculation of the time series RMS of interferograms, we propose a time series decomposition method to estimate the tropospheric noise, and mitigate the bias caused by ground displacement.

For each PS point, we decompose its corresponding values extracted from time series interferograms into three terms:

$$f\left(t\right)=\left(a{t}^2+bt+c\right)+A\sin \left(\omega t+\phi \right)+\epsilon$$

(1)

where the first term $\left(a{t}^2+bt+c\right)$ is a quadratic polynomial component modeling possible ground deformation, $t$ represents different acquisition time, and $a$, $b,$ and $c$ are the quadratic coefficients to be solved for. For PS points on the stable ground, the first term would be negligible and yield very small coefficients. For PS points in the deforming areas, according to prior knowledge and previous studies, the ground displacements are often related to land reclamation from the ocean, groundwater exploitation, tunneling works, etc., which may exhibit nonlinear deformation tendencies in addition to the usual linear displacement. Therefore, we use a quadratic polynomial term to model the possible ground deformation. For linear deformation, the estimated coefficient $a$ would become negligible, and the quadratic polynomial term is approximately equal to a linear term.

Global weather models indicate the existence of seasonality in the stratified delay, which is characterized by periodicities in the form of inter-annual or annual sinusoidal variations [8,29,30]. Hence, the second term $A\mathit{sin}\left(\omega t+\phi \right)$, denoting a sine function, is used to simulate the residual periodic oscillation due to stratified tropospheric phase, $A$ and $2\pi /\omega$ represent the magnitude and period of the oscillation, and $\phi$ is the initial phase of the sine function. The fitted sinusoidal term is then taken into account in the estimation of tropospheric noise. In the region of interest, there are some hilly areas with an elevation ranging from 500 m up to 1260 m, where the stratified phase delays can result in significant periodic effects. In low-altitude flat terrain, the tropospheric phase is mainly dominated by the turbulent component and the stratified phase delays can be ignored. In this case, the fitted sinusoidal function term becomes negligible, which does not affect the analysis of the tropospheric noise. Based on the above considerations, we retain the second term in the time series decomposition.

The third term, $\epsilon$, denotes the unmodeled noise.

All unknown coefficients $\left(a,b,c,A,\omega ,\phi \right)$ in Equation (1) are solved by a nonlinear curve-fitting method in a least-squares sense, as follows:

$$\underset{coef}{min}{\Vert F\left(coef,tdata\right)-disp\Vert }_2^2=\underset{coef}{min}{{\displaystyle \sum }}_i{\left(f\left(coef,t\right)-dis{p}_i\right)}^2$$

(2)

where $\mathbf{c}\mathbf{o}\mathbf{e}\mathbf{f}=\left[a,b,c,A,\omega ,\phi \right]$ is a vector of coefficients to be solved, $\mathbf{t}\mathbf{d}\mathbf{a}\mathbf{t}\mathbf{a}$ denotes the vector of acquisition time, $\mathbf{d}\mathbf{i}\mathbf{s}\mathbf{p}$ represents the vector containing PS time series values extracted from interferograms, $i$ is the index of an acquisition time $t$ in the time vector $\mathbf{t}\mathbf{d}\mathbf{a}\mathbf{t}\mathbf{a}$, and $F$ is a vector-valued function of the same size as $\mathbf{t}\mathbf{d}\mathbf{a}\mathbf{t}\mathbf{a}$. We analyzed the time series of the zenith total delay (ZTD) obtained from the GACOS product at PS points with an elevation difference of more than 600 m from the reference altitude, and found that the sinusoidal variation with the longest period is less than 2 years. Therefore, we imposed constraints on the upper and lower bounds of the variable ω in the nonlinear least-squares solver. Specifically, we restrict the 2π phase change to correspond to a time period of 12 to 20 months. While this constraint reduces the impact of periodic errors within this range, longer-term periodic fluctuations cannot be captured by the model fitting of the sine function and are attributed to unmodeled noise.

After the time series decomposition, the sinusoidal component (the second term in Equation (1)) and the unmodeled noise $\epsilon$ (the third term in Equation (1)) are summed up as an estimate of the residual noise after InSAR tropospheric corrections. For each PS point, an RMS value is then computed along the timeline over $N$ acquisition dates:

$$RMS=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N\Delta d_{tropo,i}^2}$$

(3)

where $\Delta d_{tropo,i}$ represents the estimated tropospheric noise of a PS point on the ${ith}^{}$ interferogram, as derived by the time series decomposition method.

It should be noted that the quadratic term simulates the main trend of ground deformation, but certain forms of land subsidence resulting from periodic changes of precipitation or groundwater levels may not be fully accounted for. However, the residual deformation signals resulting from such variations can potentially be detected by the sinusoidal term and mistakenly attributed to tropospheric noise. Due to the unavailability of long-term measurements of groundwater table and precipitation, it is difficult to distinguish between the impacts of periodic variations in precipitation, groundwater table, and tropospheric phase delays on urban subsidence. Hence, the time series decomposition method does not entirely mitigate the bias originating from ground deformation due to the possible presence of periodic deformation signals.

It is important to mention that the time series decomposition technique is exclusively used for the interferograms corrected by spatiotemporal filtering or the combined approaches. In the original interferograms and those corrected by weather model products, high levels of time-independent noise often present, which are resulted from turbulence. As a result, the deformation signal is inevitably masked, and its effect on the time series noise statistics is negligible. Hence, for these interferograms, the conventional time series RMS method is employed to analyze the data.

3.1.2. Semi-Variograms with Model Fitted Range and Sill

As the tropospheric phase is spatially correlated, we utilize semi-variograms and compute their model fitted ‘range’ and ‘sill’ values to quantitatively analyze the reduction of spatially dependent signals after tropospheric corrections.

A semi-variogram is defined as [26]:

$$S\left(\stackrel{⃑}{\mathbf{r}}\right)=E\left[{\left(f\left(x\right)-f\left(x+\stackrel{⃑}{\mathbf{r}}\right)\right)}^2\right]$$

(4)

where $\stackrel{⃑}{\mathbf{r}}$ is a vector, $x$ is the location of a PS point, $f\left(x\right)$ denotes the signal after tropospheric correction, and $E$ indicates the expected value. Under the isotropic assumption, $\stackrel{⃑}{\mathbf{r}}$ in Equation (4) can be replaced with a scalar $r$ to only represent the distance from $x$. We adopt an empirical semi-variogram approach, which bins the randomly sampled points by distance to apply Equation (4).

A semi-variogram describes the spatial correlation of paired sample points. It usually reaches a peak and flattens at a certain distance, referred to as the ‘range’. Samples separated by distances closer than the range are spatially correlated, otherwise decorrelated. The semi-variance at the range is called the ‘sill’, the square root of which is a measure comparable to the spatial RMS of the interferogram, an indicator of the spatial noise. Here we obtain the range and sill by Gaussian model-fitting on semi-variograms [31]. It should be noted that the PS points in deforming areas were not discarded prior to the calculation, thus the sill statistics can be regarded as a biased and rough estimate of spatial noise in individual interferograms.

3.1.3. Spearman’s Rank Correlation between Phase and Elevation

The phase-elevation dependence is used as the metric to assess the reduction of topography-dependent signals after InSAR tropospheric corrections. The vertical refractivity differences between different altitudes can result in a stratified tropospheric delay. A 500 m difference in elevation can cause a phase difference of up to 10 mm [32]. There are still some hilly areas in the study area, most of which lie at an altitude of over 500 m, while some can even reach 1260 m above sea level. Thus, it is necessary to inspect the reduction of stratified delays, particularly in regions with significant variations in elevation, such as hilly areas. Because a considerable portion of the study area is situated in low-altitude flat terrain, the analysis was only performed on PS points in hilly areas (drawn in red boxes in Figure 1).

The relationship between the stratified phase and elevation is often modeled by a linear or power-law function [14,33], but a previous study also suggests it is too complicated to be fully simulated [34]. A linear correlation coefficient was used to assess the phase-elevation dependence [25]. In this study, we propose the use of Spearman’s rank correlation coefficient (SRCC) as a more general measure.

SRCC is a nonparametric measure of statistical dependence between the rankings of two variables, which aims to assess the monotonic relationship between two variables, even if the relationship is nonlinear. In addition, SRCC is more robust to outliers than the conventionally used Pearson’s linear correlation coefficient [35]. For a sample of size $N$, the $N$ raw scores $X_i$, $Y_i$ are converted to ranks $R\left(X_i\right)$, $R\left(Y_i\right)$. A rank is assigned to each value based on its position in the ordered set of values of a variable. Their Spearman correlation coefficient $r_s$ is:

$$r_s=\frac{\mathit{cov}\left(R\left(X\right),R\left(Y\right)\right)}{{\sigma }_{R\left(X\right)}{\sigma }_{R\left(Y\right)}}$$

(5)

where $\mathit{cov}\left(R\left(X\right),R\left(Y\right)\right)$ denotes the covariance of the ranks, and ${\sigma }_{R\left(X\right)}$ and ${\sigma }_{R\left(Y\right)}$ are the standard deviation of variable ranks.

The correlation coefficient $r_s$ is a value between −1 and +1, indicating a perfect negative and positive correlation, respectively. A $r_s$ value of $0$ represents no correlation. Following the guide by Ref. [36], the correlation strength can be categorized using the absolute value of $r_s$, as detailed in

Table 2. The correlation strength described by the absolute value of r_s.

|rs|	Correlation Strength
0.00~0.19	Very weak
0.20~0.39	Weak
0.40~0.69	Moderate
0.70~0.89	Strong
0.90~1.00	Very strong

.

The statistical significance of SRCC can be revealed by the critical probability ($p$) values. In this study, a p-value of 0.05 (5%) or less is considered statistically significant to reject the null hypothesis ($H_0$). In addition, it has been suggested to have at least ten samples for Spearman’s correlation analysis [36,37]. Using samples fewer than this will more likely obtain a result of chance. To ensure the reliability of the phase-elevation correlation analysis, a window is determined valid by satisfying both of the following conditions: (1) containing no less than 10 points; (2) obtaining a p-value below 0.05. Then, on the original interferograms, the valid windows yielding SRCC over 0.4 (moderate correlation as defined in [36]) are preserved for analysis.

3.2. Analysis of Primary and Secondary Images’ Contribution in Tropospheric Corrections

In this study, the following process (Figure 2) is used to analyze the respective roles of primary and secondary TOEs in InSAR tropospheric corrections.

It has been pointed out by previous studies that interferograms that share a common date share the same propagation delay [2,38]. Therefore, in PS-InSAR processing, using the same reference image (i.e., the primary image) results in a time-correlated component in the interferometric phase, which is referred to as the primary image’s contribution in the StaMPS approach [2]. As stated in [2], for a PS point, after removal of a spatially uncorrelated look angle error and phase unwrapping, the differential phase can be expressed as:

$$\widehat{\phi }={\varphi }_{def}+\Delta {\varphi }_{tropo}+\Delta {\varphi }_{orb}+\Delta {\varphi }_{\theta }^c+\Delta {\varphi }_{noi}+2k\pi$$

(6)

where ${\varphi }_{def}$ is the phase change in the satellite line-of-sight due to ground deformation. $\Delta {\varphi }_{tropo}$ is the tropospheric phase; $\Delta {\varphi }_{orb}$ is the phase due to satellite orbit inaccuracies, which is negligible after applying the precise orbits by SNAP (Sentinel Applications Platform); $\Delta {\varphi }_{\theta }^c$ is the spatially correlated look angle error; and $\Delta {\varphi }_{noi}$ represents noise related to other unmodeled sources, such as the ionospheric effects, the variability in scattering, thermal noise, co-registration errors, etc. The $2k\pi$ denotes the remaining integer ambiguity.

The three spatially correlated nuisance terms $\Delta {\varphi }_{tropo}+\Delta {\varphi }_{orb}+\Delta {\varphi }_{\theta }^c$ in Equation (6) can bias the results of InSAR deformation monitoring. These terms can be estimated utilizing the characteristics of their spatiotemporal correlation. The StaMPS approach divides the three terms into the part correlated in time and the part that is not. The former is referred to as the primary TOE, as it presents in each interferogram due to the use of the same reference image. The latter consists of the spatially correlated look angle error $\Delta {\varphi }_{\theta }^c$ and the secondary TOE. The $\Delta {\varphi }_{\theta }^c$ is proportional to the perpendicular baseline and can be removed separately. The primary TOE and secondary TOE are derived by spatiotemporal filtering, as displayed in the first block diagram in Figure 2. The combination of spatial and temporal low-pass filters obtained the estimates of the primary TOE, whilst the integration of spatial low-pass filter and temporal high-pass filter resulted in the secondary TOE.

To investigate the different impacts of the primary and secondary TOEs on the InSAR measurements, we implement the three metrics to compare the different roles of the two components in the reduction of the tropospheric noise, attenuation of medium to long wavelength signals, and mitigation of topography-dependent signals.

4. Experiments and Results

4.1. Experimental Settings

The Sentinel-1 TOPS SLC data were first processed by the SNAP-ESA Sentinel Application Platform (http://step.esa.int, accessed on 12 March 2022) using the PS-InSAR strategy. The co-registration was performed using the Enhanced Spectral Diversity (ESD) method [39]. The topographical phase in each interferogram was removed using the AW3D30 DEM. The differential interferograms were then processed by the PS-InSAR workflow. The spatial reference for resolving the phase time series was set to a stable zone within the SAR coverage, with longitude between 113.940320° and 113.952058°, and latitude between 22.565002° and 22.571356°. The reference phase was determined by the mean phase of PS points within the stable zone. The selected reference area is shown in Figure 3. All interferograms were unwrapped by a 3D phase unwrapping approach [40], followed by the tropospheric correction.

For InSAR tropospheric corrections, we applied three categories of methods in the experiments, including the weather model-derived corrections utilizing ERA5, ERA-I, and GACOS products, a commonly used and representative statistical method, i.e., StaMPS spatiotemporal filtering, and the integration of the two mainstream techniques. For simplicity, the term “spatiotemporal filtering” in this paper refers to the StaMPS spatiotemporal filtering method.

The ERA-I and ERA5 reanalysis were processed by the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [41] to estimate both the turbulent and stratified delays. In TRAIN, lateral and vertical spline interpolations are applied to the atmospheric parameters provided by ERA-I and ERA5. Linear interpolation is performed in the time domain to match the SAR acquisition time [41]. Then, the tropospheric delays can be computed.

As for the GACOS product, ECMWF HRES was employed as the input atmospheric data, and the ITD model [19] was utilized to separate stratified and turbulent components to generate zenith total delay maps for tropospheric correction.

In the experiments, for a fair assessment of different tropospheric correction methods, the estimated $\Delta {\varphi }_{\theta }^c$ was subtracted from interferograms before comparison. We did not remove a best-fit ramp before or after tropospheric correction. Instead, the possible orbit error, if there is any, was estimated along with the tropospheric signals by the spatiotemporal filtering.

By applying the three statistical metrics to the time series interferograms before and after tropospheric corrections, we obtained the results in Section 4.2, Section 4.3 and Section 4.4. To investigate the roles of primary and secondary TOE in InSAR tropospheric corrections, we also compare the statistical metrics obtained with and without removing the primary TOE in tropospheric corrections using spatiotemporal filtering or combined approaches (shown in Section 4.2, Section 4.3 and Section 4.4), followed by an analysis in Section 4.5. In Section 4.6, we present case studies of local subsidence that incorporate tropospheric corrections achieved by integrating spatiotemporal filtering and weather model products.

4.2. Elimination of Overall Tropospheric Noise

For each PS point, we extracted its corresponding phase values from time series interferograms and converted them to length in the unit of millimeters, so the phase is directly related to the distance that the wave has traveled. Each PS time series was then decomposed using Equation (1). The coefficients (a, b, c, A, ω, φ) in Equation (1) are solved by a nonlinear curve-fitting method in a least-squares sense, as described in Section 3.1.1.

To illustrate the results of time series decomposition, we selected four representative PS points, herein referred to as “P1”, “P2”, “P3”, and “P4”, to serve as examples. The four points are marked in Figure 4, with latitude, longitude, altitude, and topography detailed in

Table 3. Details of the four representative persistent scatter (PS) points.

	Latitude	Longitude	Altitude	Topography
P1	N 22.523119°	E 113.889458°	−3.6 m	Low altitude, flat terrain, ocean-reclaimed area
P2	N 22.313656°	E 113.917023°	−1.1 m	Low altitude, flat terrain, ocean-reclaimed area
P3	N 22.571257°	E 114.188890°	639 m	High altitude, hilly area
P4	N 22.568163°	E 114.171921°	39.4 m	Low altitude, flat terrain, metro tunneling area

. The example results are displayed in Figure 5.

As shown in Figure 4, point “P1” (a) is located in ocean reclaimed land, i.e., the Qianhai sub-area, Shenzhen. Point “P2” is located in Hong Kong International Airport (HKIA), also built on ocean reclaimed land. Point “P3” is on Wutong Mountain, Shenzhen, at an altitude of over 600 m. Point “P4” is in a residential area nearby the Liantang metro station in Shenzhen, which was affected by tunneling works of metro construction from 2015 to 2020. There has been reported subsidence occurring in sub-areas (a), (b), and (d), while the sub-area (c) in Wutong Mountain is stable.

In Figure 5, the examples of time series decomposition were derived from interferograms after tropospheric correction using ERA5 and spatiotemporal filtering. In (a), (b) and (d), we can see that “P1”, “P2” and “P4” all exhibit significant subsidence trends with weak periodic effects, as the three points are located in low-altitude flat terrain. The “P3”, at an altitude over 600 m, shows periodic oscillation, albeit with a relatively low amplitude. Upon analyzing the time series of phase differences between the new “P3” and its adjacent point, it was found that there were no occurrences where the phase differences exceeded the value of $\pm \pi$, indicating the relatively reliable nature of the calculation results for this point. Thus, we are inclined to attribute them to the influence of residual stratified phase delays. Since we imposed constraints on the upper and lower bounds of the variable ω in Equation (1) in the nonlinear least-squares solver, longer-term fluctuations such as the residual periodic error observed at P4 was not captured by the model fitting of the sinusoidal term and was attributed to unmodeled noise. The observed residual periodic signal in the fitting outcomes of P4 cannot be accounted for by the stratification, given its low altitude position. The longer-term sinusoidal variations may be caused by other factors such as variations of precipitation or groundwater table. In addition to the aforementioned phenomenon, it is observable that all four time series curves display a satisfactory fit, and the unmodeled noise generally fluctuates around zero (Figure 6) after time series decomposition. From the above results, it is evident that the time series decomposition method is able to reduce the bias caused by ground displacement to obtain a more accurate estimate of residual tropospheric noise.

Then, we estimate the tropospheric noise by the sum of the sinusoidal component and the unmodeled noise term. A time series RMS value is computed for each PS point on the tropospheric noise derived by the time series decomposition method.

Taking the Qianhai sub-area as an example, we compare the time series RMS calculated based on the time series decomposition method with those derived using the none-deformation assumption, as displayed in Figure 7. The color bar in Figure 7a was determined using the “natural breaks classification”. The color bar in Figure 7b was manually defined according to the break values used in (a) and adjusted using its own value range, so as to make the two color bars comparable. In Figure 7, it is evident that the proposed time series decomposition method yields a smaller value range of RMS, especially on PS points prone to land subsidence (highlighted by the white rectangle) identified in the previous work [27].

To evaluate the residual tropospheric noise, we calculate the percentiles of PS time series RMS derived by the time series decomposition method over the whole study area, as shown in Figure 8. Due to the big differences between value ranges of different tropospheric correction methods in the vertical axis, the percentile plots were split into two subplots.

In Figure 8a, the results suggest that the weather model-derived corrections help to remove part of the noise in the original interferograms. ERA5-derived correction achieved the most reduction of RMS among the three weather model products. The spatiotemporal filtering suppressed the noise better than the weather model data, perhaps because it was implemented with a much higher spatial resolution (tens of meters). The integration of weather model products with spatiotemporal filtering outperforms either of the two. In addition, it is noticed that the removal of primary TOE further reduces the time series RMS values, but not significantly. This indicates it is necessary to subtract the primary TOE, but it is not a dominant source of tropospheric noise over the whole study area.

4.3. Mitigation of Distance-Dependent Signals

Firstly, we calculated the empirical semi-variograms using 200 bins (with a bin tolerance of about 386 m) on interferograms before and after tropospheric corrections. Figure 9 displays the example semi-variograms plotted for the 20170312–20170815 interferogram. Due to the large disparity between the value ranges of different groups, we display them in two subplots. In Figure 9, it is intuitive that the semi-variogram values are dramatically reduced by the correction using spatiotemporal filtering or the combined approaches. It is worth noting that after removing the primary TOE, the mean semi-variogram values decrease by more than one order of magnitude.

To quantitatively assess the reduction of distance-dependent signals, all semi-variograms were fitted using a Gaussian model using the least squares estimation method [42]. The output of each includes a decorrelation range, a sill, and an R square value describing the goodness of fit [43]. Only the fits with an R square value over 0.6 were kept for analysis. Then, grouped by different tropospheric correction methods, a weighted mean is calculated using all available fits in each group using the corresponding R square values as the weights. The results are shown in

Table 4. The model-fitted mean decorrelation range and sill of semi-variograms before and after tropospheric corrections. Entries denoted by ‘*’ preserve the primary TOE in the interferograms.

		Weighted Mean Range (km)	Weighted Mean Sill (mm)
Original	No correction	53.20	761.39
Group 1	GACOS	48.19	561.66
	ERA-I	50.57	489.88
	ERA5	51.50	563.68
Group 2	Spatiotemporal filtering *	63.57	223.61
	GACOS and filtering *	76.61	71.92
	ERA-I and filtering *	61.41	63.82
	ERA5 and filtering *	63.57	134.54
Group 3	Spatiotemporal filtering	15.97	2.79
	GACOS and filtering	16.75	3.34
	ERA-I and filtering	13.72	2.86
	ERA5 and filtering	13.78	3.50

.

According to the results in Table 4, considering the sparse spatial resolution of the weather model data and the uncertainties in the model fitting of semi-variograms, we divided the different correction methods into three groups to compare the performance in the mitigation of spatially correlated signals, i.e., the weather model-derived corrections, the combined approaches with the primary TOE retained, and the combined approaches after removing the primary TOE. In Table 4, the original semi-variograms (derived from interferograms without tropospheric correction) yield the largest sill value, the square root of which is a rough measure comparable to the spatial noise of individual interferograms. The semi-variograms corrected by the three weather model products exhibit a slightly reduced mean decorrelation range and a considerable decrease in sill values. It is found that the primary TOE significantly affects the decorrelation range and sill. With the primary TOE retained, the spatiotemporal filtering and the combined approaches, although further reducing the mean sill compared with the weather model derived corrections, all yield even larger decorrelation ranges. After subtracting the primary TOE, the corresponding range and sill values decreased dramatically.

4.4. Reduction of Phase-Elevation Dependence

As described in Section 3.1.3, Spearman’s rank correlation coefficient was used to assess the phase-elevation dependence. Since a considerable part of the study area is located in low-altitude flat terrain, the analysis was conducted only on PS points in hilly areas with varying altitudes (highlighted by the red boxes in Figure 1), so as to capture the phase changes with height.

Figure 10 compares Spearman’s rank correlation and Pearson’s linear correlation coefficients derived from the same interferograms. In Figure 10, the subplots (a) and (b) each display a scatterplot of phase and elevation with two least-squares reference lines. Although Spearman’s rank correlation has no corresponding lines, to explicitly show its difference from Pearson’s linear correlation, we plot the two reference lines. The slope of each line is equal to the displayed Spearman’s or Pearson’s correlation coefficient. Outliers have been marked with black circles or ellipses. In Figure 10a, outliers can be found in both the top right and bottom left corners of the subplot. We can see that the slope of the blue line, which corresponds to Pearson’s correlation coefficient, is greater than that of the red line, which corresponds to Spearman’s correlation coefficient. Figure 10b shows that the slope of the blue line is much smaller than the slope of the red line, as the outliers are all located at the bottom of the image. The results suggest that Spearman’s rank correlation is less sensitive to outliers, while Pearson’s linear correlation is prone to larger deviations when outliers appear.

To generally analyze the reduction of topography-dependent signals, for all of the PS points situated in hilly areas (delineated by red boxes in Figure 1), we calculated the phase-elevation rank correlation using different spatial scales (window sizes) varying from 10 to 80 km. For each spatial scale, a mean correlation coefficient is computed over all interferograms, as shown in Figure 11.

To investigate the potential influence of disparate primary images on experimental outcomes, we employed an alternative primary image acquired on 4 February 2017 for cross-validation. The results area is displayed in Figure 12.

It should be noted that weather model-derived corrections do not involve spatiotemporal filtering, with both primary and secondary TOE preserved. From Figure 10b and Figure 11b, we can see that the use of a different primary image did not alter the comparison relationship of phase-elevation correlation between the original interferograms and the interferograms corrected by weather model products. After the weather model-derived corrections, the phase-elevation correlation coefficients all decreased. Among the three weather model products, GACOS achieved the most reduction in the phase elevation correlation.

Upon comparing Figure 11a and Figure 12a, we observed that the use of different primary images resulted in significant changes to the phase-elevation correlation coefficient in interferograms corrected through spatiotemporal filtering or its integration with weather model products. This outcome is expected as the interferometric phase contains primary TOE that is unique to the primary images acquired on different dates.

From Figure 11 and Figure 12, it is evident that the removal of primary TOE led to a substantial reduction in the phase-elevation dependence, regardless of the primary image used. In both Figure 11b and Figure 12b, after removing the primary TOE, the combined approach using any of the three weather model products reduced the phase-topography correlation to a similar extent, and GACOS and spatiotemporal filtering exhibits a slightly superior performance in comparison to the other combined approaches. In Figure 11b, it can be observed that the spatiotemporal filtering method still displays a slightly higher phase-elevation correlation than its combination with weather model products at the spatial scales of 10 km and 20 km. In contrast, this trend is evident across all spatial scales above 30 km in Figure 12b.

4.5. The Roles of Primary and Secondary TOE in Tropospheric Corrections

As revealed by the decorrelation range of semi-variograms (Table 4), the primary TOE significantly affects the decorrelation range of semi-variograms in the corrections using spatiotemporal filtering or combined approaches. As pointed out by previous studies, tropospheric delays are composed of a short-scale (few kilometers) turbulent component due to strong variation of water vapors, a long-scale turbulent component (tens of kilometers) attributed to the lateral variations of temperature, pressure, and humidity, and a long-scale stratified component resulted from vertical variations of temperature, pressure, and humidity [26,41]. We can reasonably surmise that the primary TOE is an estimate of long-scale signals, while the secondary TOE is more related to short-scale turbulence. Hence, removing the secondary TOE led to a considerable decrease in mean sill values due to the reduction of short-scale signals. When the primary TOE was preserved, the long-scale signals remained, so the decorrelation range did not decrease.

According to the phase-elevation correlation analysis (Figure 11 and Figure 12), the primary TOE plays a leading role in reducing the phase elevation correlation. In the combined approaches, the dependence decreased dramatically after removing the primary TOE. The reason can be inferred as follows: the primary TOE is obtained by temporal low-pass filtering from the spatially correlated nuisance terms. At low altitudes and flat terrain, the primary TOE corresponds to the temporally correlated part of the tropospheric delay introduced by the use of a common reference (primary) image in the PS InSAR approach. In middle-to-high altitude hilly areas, given that the stratified component is temporally correlated, the primary TOE actually estimates part of the stratified delays.

Regarding the statistics of PS time series RMS (Figure 8), the removal of primary TOE reduces the tropospheric noise, but not significantly, suggesting that the primary TOE is not the main source of tropospheric noise over the whole study area. However, it still substantially affects distance-dependent signals and topography-dependent signals. In the conventional view, the statistical correction methods do not take into account the stratified delay, which is distance-dependent and topography correlated. Through our investigation, we found that, in fact, the primary TOE in StaMPS spatiotemporal filtering estimates longer wavelength signals, and undertakes part of the estimates of the stratified delays in the middle to high-altitude areas.

4.6. Local Subsidence Maps Derived after Tropospheric Correction

In this section, we present examples of the local subsidence maps (Figure 13) that were obtained after the tropospheric correction using ERA5 and spatiotemporal filtering.

Figure 13a shows the reclaimed Qianhai area, which encircles Qianhai Bay, and is also known as the Qianhai Free Trade Zone. The Qianhai area underwent substantial land reclamation predominantly from 1979 to 2010, with minimal alteration to the shoreline observed after 2010. The annual displacement rate in Qianhai area reaches −13 mm/yr, with the most significant subsidence outlined in the black circle. Figure 13b illustrates the land subsidence observed in HKIA. The HKIA, as one of the largest artificial islands in the world, was reclaimed on the foundation of two original rocky islands (Chek Lap Kok drawn with white lines and marked as “A”, and Lam Chau outlined with white lines and marked as “B”) that formerly located off the northern coast of Lantau Island. Local settlements mainly present in the reclaimed area, with mean deformation rate ranges from 0 to −13 mm/yr. Figure 13c presents the subsidence monitored in a residential area during the construction of Metro Line 2 in Shenzhen. The construction of the metro tunnel spanned from 2015 to 2020, with the excavation of the shaft beginning in the latter half of 2016. The location of the shaft excavation is demarcated by the black pentagram, and the segment of the metro tunneling construction is indicated by the gray line. We have observed a significant subsidence bowl at the site of shaft excavation. In addition, more subsidence spots have been identified near the construction segment of the metro tunnel.

5. Discussion

5.1. The Applicability and Limitations of Different InSAR Tropospheric Correction Methods

Reflected by the time series RMS metric (Figure 8), the tropospheric noise was reduced by the weather model-derived corrections, but still higher than that corrected by the spatiotemporal filtering. This is likely due to the coarse spatial resolution of weather model products, which currently limits the ability of weather model-derived corrections to mitigate the short-scale turbulence with a wavelength smaller than the horizontal grid spacing of weather model data [44], leading to higher residual tropospheric noise. The combined approaches achieved the greatest reduction of tropospheric noise.

The statistics of the semi-variograms (Table 4) suggest that weather model-derived corrections mitigated part of the long-scale signals, but not completely. This may be due to the inaccuracies of the weather model corrections, as the spatial resolution of the weather model data is too coarse compared to the decorrelation range of the original semi-variograms (~53 km). In contrast, the primary TOE derived by the spatiotemporal filtering provides a more accurate estimate of long-scale signals, which should be attributed to the small window size (typically set to 800 m~1 km) used in the StaMPS spatial low-pass filter.

The analysis of phase elevation correlation (Figure 11) indicates that, in spatiotemporal filtering, the primary TOE has an important impact on phase elevation dependence. In the past, it was generally accepted that statistical methods (e.g., stacking, spatiotemporal filtering) assumed a random distribution of tropospheric phase delays along the timeline. However, in the employed StaMPS PS-InSAR approach, only the secondary TOE has this characteristic. It is not surprising that the removal of secondary TOE did not reduce much of the topography-dependent signals. The primary TOE, which is assumed correlated both in time and space, actually attributes to the estimate of topography-related signals in the middle-to-high altitudes. Despite the significant reduction of topography-dependent signals achieved by using weather model-derived corrections, our findings indicate that their integration with spatiotemporal filtering can lead to a further reduction of the phase-elevation dependence.

One direction for our future work is to use mesoscale and/or regional atmospheric reanalysis where applicable. This is expected to further improve the performance of weather model-derived corrections and the combined approaches.

5.2. Comparison of ERA-I, ERA5, and GACOS for the Tropospheric Correction in Coastal Areas

As revealed by the time series RMS metric (Figure 8), the GACOS correction exhibits the highest residual tropospheric noise among the weather model-derived corrections (Figure 8). The GACOS correction only reduces the noise of original interferograms when it is above about 40% percentile. We consider that lower percentile points, i.e., smaller time series RMS, correspond to calm tropospheric conditions. For higher percentile points, i.e., interferograms with stronger signals, GACOS correction exhibits positive results. Therefore, it can be inferred that the performance of GACOS correction is sensitive to the magnitude of the tropospheric signal. Another factor could be that the current version of GACOS uses HRES ECMWF analysis as input to generate zenith total delay maps. ERA-I and ERA5, as reanalysis, both benefit from more observations and more rigorous quality control than analysis data. In addition, we notice that the tropospheric noise after ERA-I correction is slightly lower than that corrected by ERA5 below the 50% percentile, while it is higher above the 50% percentile. These observations indicate that under strong tropospheric signals, the ERA5-derived correction is more effective in mitigating the tropospheric noise.

In the inspection of residual topography-correlated signals (Figure 11 and Figure 12), GACOS minimized the phase-elevation dependence among the three weather model products. This finding suggests that GACOS, leveraging the Iterative Tropospheric Decomposition (ITD) model, has superior performance in the decoupling of turbulent and stratified tropospheric delays, resulting in more accurate removal of the topography-dependent signals. In the future, by incorporating the ITD model with new global atmospheric reanalysis, the GACOS product is expected to achieve better performance of InSAR tropospheric correction.

6. Conclusions

In this paper, we come up with a guideline for a statistical assessment of InSAR tropospheric corrections. Firstly, a time series decomposition method was proposed to reduce the bias due to ground displacement, so as to facilitate a more accurate estimate of tropospheric noise from time series analysis of interferograms, without relying on the none-deformation assumption. On this basis, three statistical metrics were adapted to evaluate the performance of tropospheric corrections: (i) the PS time-series RMS value to assess the elimination of overall tropospheric noise, after the time series decomposition; (ii) the semi-variograms with model-fitted decorrelation range and sill to reflect the reduction of distance-dependent signals; (iii) Spearman’s rank correlation between phase and elevation to evaluate the mitigation of stratified tropospheric delays in hilly areas.

The applicability and limitations of the weather model-derived corrections, the widely used StaMPS spatiotemporal filtering method, and combined approaches were systematically assessed. Experiments were conducted on the east bank of the Pearl River Estuary, employing the advanced and globally available weather model products, including ERA5, ERA-Interim, and GACOS. It was found that the employed global weather model products cannot remove the short-scale turbulence due to their coarse grid spacing, but are effective in the mitigation of long-scale signals and stratified delays. The StaMPS spatiotemporal filtering works well in the mitigation of short-scale water vapor. The combined approaches outperform either of the two mainstream methods in the elimination of overall tropospheric noise. Comparing the three weather model products, ERA5 and ERA-I achieved better performance in the suppression of overall tropospheric noise, while GACOS outperformed the other two in the mitigation of topography-dependent signals.

Furthermore, we highlight that the primary and secondary TOEs, resulting from the StaMPS spatiotemporal filtering embed in PS InSAR processing, take on different roles in the InSAR tropospheric corrections. The primary TOE estimates long-scale tropospheric signals, while the secondary TOE is related to the estimation of short-scale turbulence. In hilly areas, the primary TOE, as derived by its time-dependent features, undertakes part of the estimates of long-scale stratified delays.

With the success of the European Space Agency’s Sentinel-1 and the subsequent support from upcoming missions, the volume of satellite data will be continuously increasing. Numerous automated/intelligent InSAR processing systems have been under development. The presented statistical metrics will enable the comprehensive and quantitative evaluation of InSAR tropospheric correction, and facilitate the intelligent system to concentrate on the interpretation of ground surface displacement. This, coupled with the availability and improvements of regional/global atmospheric reanalysis, provides great opportunities for intelligent InSAR monitoring of the earth surface from local, regional to large-area applications.