The SAR data contains one stack of C-band ENVISAT ASAR (ASAR-desc) images and two stacks of C-band Sentinel-1 (S1-desc and S1-asc) images. The former is acquired in stripmap (SM) mode from april 2007 to February 2009 in descending track 170. The latter are acquired in interferometric wide swath (IW) mode from May 2015 to July 2018 in descending track 71 and ascending track 64, respectively. The white rectangles in Figure 1 mark the spatial coverages of the three SAR data stacks. The characteristics of the three SAR datasets are elucidated in Table 1. The Shuttle Radar Topography Mission (SRTM) DEM of approximately 90-m resolution is used for the removal of topographic phases.

The Southern California Integrated GPS Network (SCIGN), one of the most densely spaced GPS networks in the world (Hudnut et al., 2001), is constructed to monitor the seismic activities in entire Southern California. Figure 1 shows the locations of 125 permanent GPS stations indicated by black rectangles. The GPS deformation results can be found from the website (http://www.scign.org), which are used to validate the deformation results of InSAR with tropospheric delay corrections.

Four meteorological reanalysis models (ERA5, GACOS, MERRA2, and NARR) and one numerical weather prediction model (WRF) are used to obtain the tropospheric propagation delays. Table 2 lists the parameters of the four meteorological reanalysis models.

As a new generation of mesoscale numerical weather forecast model, the Weather Research and Forecasting (WRF) model is dedicated to atmospheric simulation and weather forecasting research (Skamarock and Klemp, 2008). The WRF (ARW) 4.0 prediction model developed by National Center for Atmospheric Research (NCAR) is used to simulate atmospheric conditions. The Final Operational Global Analysis (FNL) data of the National Environmental Forecast Center provides initial values and boundary conditions. The experimental simulation runs for 18 h. To enhance the resolution, we set up a three-layer nesting scheme for the WRF simulation experiment, and from the initial grid resolution to the inner nested domain resolution are 27, 9, and 3 km, respectively. The tropospheric delay is estimated using data with a spatial horizontal resolution of 3 km and 37 layers in the vertical direction (1,000–100 hPa).

The Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua satellites of NASA is a passive imaging spectrometer including five near-infrared bands that are sensitive to water vapors. The spatial resolution of water vapor products is about 1 km × 1 km, and the water vapor accuracy is claimed to be 5–10%. If the daytime and clear weather are not met, the data accuracy will be worse (Gao and Kaufman, 2003). Given the wet delay estimated by MODIS is overestimated compared with the GPS zenith delay, they need to be corrected before application (Li et al., 2005).

When predicting tropospheric phases using the weather models, the interpolations in space and time generally result in biased tropospheric delay estimation. We introduce the SVS algorithm to alleviate the biases and to make it closer to the correct value of tropospheric delay for large-scale areas (Shen et al., 2019). The SVS algorithm is repeated simply here to keep the methodology intact.

We divide each interferogram into N regular non-overlapping windows. The window size is set as 50 km here. The relationship between unwrapped phase ∆ϕ _unw and tropospheric delay phase ∆ϕ _tropo predicted from weather model data in a single window is taken as linear: Δ ϕ u n w = K i Δ ϕ t r o p o + C i                 ( i ∈ N ) (2) where K _i and C _i indicate the scaling factor and constant offset at window i, respectively.

The scaling factors for each window are smoothed through a weighting strategy, which is carried out as follows:

(1) The first step is to weigh each window through a signal-to-noise ratio. Var(ϕ _tropo (i)) represents the variance of tropospheric delay phase in window i, i.e. signal. Var(ϕ _res (i)) is the variance of residual phase, i.e. noise. The residual phase equals to unwrapped phase subtracted by the tropospheric delay phase. The signal-to-noise ratio of tropospheric delay phase in window i is expressed as:

The smoothed scaling factor K (x,y) is used to replace the one in Eq. 2 to obtain the scaled tropospheric delay phase for each weather model.

The tropospheric delay phases predicted from the weather model data and spectrometer observation products are refined in Scaled tropospheric phases estimated by weather model data. Since these methods present unstable performances for different SAR data at different regions, Dong et al. (2019) proposed an optimal weighting fusion of multiple tropospheric delays using a single weight. However, in the case of wide coverages (larger than tens of kilometers), the single weight is not sufficient to describe the spatially lateral variation of tropospheric delays for the entire scene (Bekaert et al., 2015a).

We divide the large scene into N windows using the same window size and layout as described in Scaled tropospheric phases estimated by weather model data and derive weights for each window. The final tropospheric delay phase ∆ϕ _tropo (x,y) on pixel (x, y) can be obtained through weighting fusion of M tropospheric delay corrections: Δ ϕ t r o p o ( x , y ) = ∑ j = 1 M w j ( x , y ) ∗ Δ ϕ t r o p o j ( x , y )         ( j ∈ M ) (6) where ∆ϕ ^j _tropo (x,y) denotes the tropospheric delay phase estimated by correction model j with corresponding weight w ^j (x, y). The weight is calculated from the root-mean-square (RMS) value δ r m s of the unwrapped phase subtracted by the tropospheric delay phase. The standard deviation (STD) σ _rms of RMS values of all correction models is also calculated for each interferogram to remove outliers. The weight w ^i,j of correction model j in window i is derived as follows: w i , j = { 0 , δ r m s i , j > σ r m s i 1 / δ r m s i , j   p ∑ j M 1 / δ r m s i , j   p , δ r m s i , j ≤ σ r m s i ,     ( i ∈ N ) (7) where the p is set as two here. If the RMS value δ r m s i , j is larger than the standard deviation σ r m s i , the weight becomes zero. Therefore, the tropospheric delay estimated by various correction methods is weighted according to the RMS of residual phase after correction. The weighting strategy achieves the adaptive fusion of multiple delay estimations. Once the weight is determined for each window, we use the same smoothing strategy in Scaled tropospheric phases estimated by weather model data to obtain the smoothed weight w ^j (x, y) for each pixel to prevent jumps of tropospheric delays between adjacent windows. We tested Gaussian filter widths of 30 km, 50 km, and 70 km, respectively, and found the optimal width of 30 km. This implies that the closer pixels have greater contributions to the weight.

We made a statistical analysis of the RMS values of residual unwrapped phases corrected by all tropospheric delay correction methods to assess their performance. In terms of S1-desc data (Figure 5), the mean RMS values for ERA5, GACOS, WRF, MERRA2, NARR, MODIS, Linear, and Powerlaw range from 3.24 rad to 10.42 rad, and the STD of RMS values range from 1.02 rad to 5.21 rad. The AFMTD method has 320 minimum RMS values out of all 322 interferograms. The mean and STD of RMS values are 2.15 rad and 0.56 rad, respectively. Both the mean and STD declined by more than 30% compared to any single model.

The statistical histograms of RMS values on the S1-asc data in the supplementary material (Supplementary Figure S1). There are 450 out of 468 interferograms having the minimum RMS values after AFMTD correction. The mean and STD are 3.09 rad and 0.98 rad, respectively, which decreased by more than 25% compared to any single model.

The STD of the S1-asc RMS values is larger than the S1-desc, indicating stronger fluctuation of the former. Besides, the S1-asc RMS presents an approximately periodic fluctuation (Supplementary Figure S2), with the maximum and minimum values appearing in summer and winter, respectively. The acquisition time of S1-asc data is 01:50 UTC, corresponding to the local time of 17:50 at dusk. The S1-desc data is obtained at 13:52 UTC, corresponding to the local time of 05:52 a.m. at dawn. The stronger variation of water vapor at dusk than at dawn may cause the discrepancy between S1-desc and S1-asc data.

We evaluate the effects of tropospheric delay correction on the InSAR annual mean linear deformation rate. The original deformation rates of S1-desc data and the ones corrected using ERA5, GACOS, WRF, MERRA2, NARR, MODIS, Linear, Powerlaw, and AFMTD are plotted in Figure 6. There are several active geological faults in the study area, among which the most active San Andreas Fault (SAF) is the boundary between plates, driving much of the geological disasters (Tong et al., 2013). The deformation rate reflects the crustal movement around the SAF fault and the subsidence caused by non-tectonic activities, such as the subsurface fluid extraction.

The effectiveness of AFMTD method for InSAR deformation was validated by continuous GPS observations. Excluding the reference station, we collected the observations at 125 GPS stations from 2015 to 2018 from the website of http://www.scign.org. The three-dimensional GPS time series were projected onto the Sentinel-1 LOS direction (Hanssen, 2001). We averaged the InSAR deformation rate within a 300 m radius at each GPS station and compared them with the GPS observations.

We estimated the GPS deformation rate along the SAR LOS direction on the 125 GPS stations by the least square method, as shown in Figure 6. The InSAR deformation rate corrected by AFMTD method shows the best consistency with GPS in terms of both pattern and magnitude. Figure 7 shows the scatterplots of InSAR and GPS deformation rates on 125 GPS stations, as well as the RMSE and R2 provided in each sub-plot. The AFMTD method has a minimum RMSE of 3.25 mm/yr. The RMSEs of other models vary from 3.29 mm/yr to 7.92 mm/yr.

Although the NARR and Powerlaw methods have high R2 of 0.87 and 0.89, there are large discrepancies between their deformation rates and GPS observations as shown in Figures 7F,I. The ERA5 and NARR overestimated the deformation rate in the north side of SAF fault (Figures 6B,F and Figures 7B,F), while the WRF and MERRA2 underestimated the deformation rate (Figures 6D,E and Figures 7D,E). These may be caused by different assimilation algorithms and spatio-temporal resolutions of each weather model data. Different from WRF and MERRA2, the rate underestimation by the linear method is located in mountainous areas (Figure 6H and Figure 7H). The deformation rate discrepancy between GPS and MODIS correction is caused by the inconsistent acquisition time of MODIS data with SAR data (Figure 6G and Figure 7G). The deviation of deformation rate corrected by Powerlaw method from GPS may result from the inaccurate powerlaw coefficient estimation (Figure 6I and Figure 7I).

The original deformation rate and the ones corrected by every single model and AFMTD for the S1-asc data and the scatterplots between InSAR and GPS are put in the supplementary material (Supplementary Figure S3 and Supplementary Figure S4). Similar to the S1-desc result, the S1-asc deformation rate after AFMTD correction is most accurate relative to the GPS observations.

We calculated the RMSE of InSAR displacement time series with and without tropospheric delay corrections relative to the GPS observations. 123 out of 125 GPS stations are presenting reduced RMSE values for the S1-desc (Figure 8A). Figure 8B shows the improvement of RMSE corresponding to each GPS site. The positive and negative values indicate decreased and increased RMSE after AFMTD correction. The AFMTD method achieved an average RMSE reduction of 35.79%, and the biggest improvement is more than 75%.

We compared the displacement time series between InSAR and GPS observations on 12 GPS stations (black squares in Figures 9A,B). Compared with the original InSAR displacement time series, the displacement time series after AFMTD correction becomes much closer to the GPS observations, which can also be confirmed by the RMSE marked in each sub-plot (Figure 9C).

For the S1-asc data, 116 out of 125 GPS stations are showing decreased RMSE after AFMTD correction (Supplementary Figure S5). The average and maximum improvements of RMSE values are 36.28% and 73.16%, respectively. The comparison of displacement time series on 12 GPS stations for S1-asc data is shown in Supplementary Figure S6. The InSAR time series corrected by AFMTD shows good agreement with GPS observations, especially, the seasonal variation of deformation time series was effectively reduced on sites BSRY and RAMT after the AFMTD correction. The InSAR time series fluctuations of S1-asc data are greater than the S1-desc data. This may be explained by the fact that the periodic phenomenon of tropospheric delay in the S1-asc acquisition time is more obvious, which can be demonstrated using the phenomenon in Supplementary Figure S2.