1. Introduction
Interferometric synthetic aperture radar (InSAR) is an earth observation (EO) technique that combines synthetic aperture radar (SAR) technology and interferometry. After periodically obtaining and aligning images of the same area, interferometric images are generated based on the imaging geometry relationship between the SAR satellite and the ground targets; then, the monitoring of crustal deformation and seismic activities can be carried out [
1]. InSAR technology has attracted much attention, due to its advantages of large coverage and high spatial resolution; in addition, it can provide all-weather, all-day, high-precision deformation monitoring, especially in some geographical areas that are not covered by the traditional measurement methods. Therefore, the data processing accuracy of InSAR technology for atmospheric delays is improving and the technology is receiving increased attention (see the
Supplementary Materials). In 1997, Zebker et al. [
2] indicated that the deformation errors caused by atmospheric delays attributed to water vapor alone can reach 10 cm to 14 cm, a degree of error that is difficult to ignore in monitoring millimeter-scale crustal deformation. Therefore, atmospheric delay correction must be carried out for InSAR data processing (see the
Supplementary Materials).
We investigated the developmental history of atmospheric delay correction techniques via InSAR measurement, from the use of InSAR’s own image data to the use of external data to correct the atmospheric delay effect on deformation monitoring. We also examined single-technology and multi-technology fusion correction methods. In each step, we witnessed the persistence of researchers with respect to the atmospheric delay correction of InSAR measurements and the gradual improvement of accuracy.
In this paper, we built a new atmospheric delay correction model by fusing the newer numerical meteorological forecasting model and the GNSS stations of the crustal movement observation network of China, based on a more reasonable weighting approach, and verified the accuracy of the new model by comparing it with the generally used generic atmospheric correction online service (GACOS) model for InSAR (
http://www.gacos.net/, accessed on 11 May 2021), under the same experimental conditions and performing the deformation validation on GNSS stations.
In 1994, Massonnet [
3] studied the Landers earthquake and identified the presence of atmospheric signals in the interferogram of repeated tracks. After that, there were more and more studies on the correction of atmospheric delay errors in InSAR measurements. The methods that were used can be divided into two main categories: one category is the correction of errors by on the basis of a method’s own image data without the involvement of meteorological parameters; the other category is the use of external data in making corrections. The first category includes the phase accumulation method (stacking) [
4], the permanent scatterer (PS) technique [
5], and the small baseline subset (SBAS) technique [
6]. Among these three methods, the phase accumulation method may reduce the temporal resolution of InSAR measurements; the PS and SBAS techniques may waste a large amount of SAR image data, and the processor’s subjective judgment is arbitrary. The second category of methods utilizes external data from the ground-based global navigation satellite system (GNSS), satellite-based medium resolution imaging spectrometer (MERIS) and moderate-resolution imaging spectrometer (MIRIS) data, moderate-resolution imaging spectroradiometer (MODIS) observations, numerical weather prediction models, ground-based meteorological station observations [
7], and wireless sounder observations [
8]. Among these methods, the use of GPS data is the primary source of information.
In 1997, Zebker et al. [
2] pointed out that when spatial and temporal variations in atmosphere relative to humidity reaches 20%, deformation measurement errors of 10 cm to 14 cm and elevation measurement errors of 80 m to 290 m can result. In 1998, Williams et al. [
9] used GPS water vapor data obtained from the GPS monitoring network in southern California to reduce the influence of the atmosphere on InSAR interferogram maps. Their results confirmed that GNSS-derived water vapor data could provide a more thorough and rigorous modeling of the effects of water vapor on InSAR interferometric data after a certain spatial interpolation, and that the GPS zenith delay was correlated with elevation, which was consistent with the conclusion of Zebker [
2] and the statistical model of Treuhaft-Lanyi [
10,
11]. Li et al. [
12,
13] and Luo et al. [
14] proposed different interpolation models and investigated the relationship between atmospheric delay and topography. They found that the use of the topography-dependent turbulence model (GTTM) [
15,
16] could result in the interpolation correction map being more accurate.
In 2011, Jolivet et al. [
17] conducted InSAR atmospheric delay correction in the Kunlun Mountains region, using reanalysis information from the European Centre for Medium-Range Weather Forecasts (ECMWF). They calculated the ratio of phase/elevation before and after atmospheric delay correction within a certain displacement window range. They found the monitoring accuracy was improved by 2.0 rad/km. However, their study also found that the atmospheric delay values obtained from meteorological models alone were not necessarily reliable [
18].
In 2014, Wang et al. [
19] used the ECMWF meteorological forecast model to correct the effects of INSAR atmospheric delay errors and found that the ECMWF model could not accurately reflect the details of local atmospheric changes, especially in regions with large changes in relative humidity. The ECMWF model was more suitable for regions with relatively small changes in large-scale meteorological conditions [
20]. Regardless of whether GNSS data or other external data, such as data derived from numerical weather prediction models, are used, they need to match the temporal and spatial resolution of InSAR-obtained image data. However, there are no data that can achieve complete agreement with InSAR-obtained image data, in terms of spatial and temporal resolution. Therefore, it is necessary to combine multiple data to develop a better InSAR atmospheric delay correction model.
In 2007, Li et al. [
21] analyzed the Gaussian morphology, anisotropy, and energy spectrum characteristics of atmospheric noise by four interferograms in the Shanghai area, and theoretically confirmed that an external data spatial resolution of at least 0.3 km is required in order to correct 90% of atmospheric delay errors in InSAR measurements. However, GNSS stations are discrete and unevenly distributed, requiring the use of other data to achieve higher spatial resolution. Thereafter, many scholars [
22,
23,
24,
25,
26,
27,
28] studied the integration of GPS and MODIS/MERIS water vapor data for better atmospheric delay correction. However, both MODIS and MERIS are near-infrared water vapor measurements, which can only describe the amounts of column water vapor on clouds in clear daytime land areas, as well as in cloud-free land and ocean areas. Currently, MERIS water vapor data are not freely available.
Therefore, we built an atmospheric delay correction model for InSAR measurements by integrating newer numerical weather prediction model reanalysis data and ground-based GNSS observation data, based on a variance component estimation (VCE) weighting method. The model was validated and compared with the current widely used generic atmospheric correction online service (GACOS) model for InSAR measurements [
29], and then applied to deformation monitoring of the 2019 Changning M6.0 earthquake that occurred in Sichuan, China, to evaluate its monitoring accuracy.
We emphasize that atmospheric delay errors exist in a variety of measurement techniques; such errors are among the important sources of error that affect the accuracy of crustal deformation monitoring. Our atmospheric delay correction model can be extended and used in the data processing of various types of crustal deformation monitoring. Our model is not limited seismic deformation applications. Although the magnitude of the seismic deformations is generally larger than that of atmospheric deformations, atmospheric delay is significant in monitoring crustal deformation via InSAR measurements, especially in regions with a large tropospheric influence, such as those regions near the equator. Atmospheric delay is also significant in monitoring urban surface subsidence.
2. Study Area and InSAR Data
Changning is a region of complex terrain in the Sichuan province of China. It is located on a secondary fault near the Changning-Shuanghe back-slope tectonics on the edge of the Sichuan Basin, with a northwest–southeast trending stress field that is prone to triggering landslide hazards [
30]. On 17 June 2019, there was an M6.0 earthquake with an epicenter located at 28.34°N (north latitude) and 104.90°E (east longitude). It had a 16 km source depth (
http://www.ceic.ac.cn/history, accessed on 29 January 2021). The study area of interest is near the epicenter.
Figure 1 shows the range and the topography of the study area. We monitored co-seismic deformation by InSAR measurements from the Sentinel-1A satellite and calculated the atmospheric delay correction model for InSAR measurements in the study area in order to improve the accuracy of crustal deformation monitoring.
The Sentinel satellite is another epoch-making C-band SAR satellite, after ERS and ENVISAT satellites. It can provide all-weather, all-day, high spatial resolution remote sensing data with good spectral quality. Due to its short wavelength, the C-band can be ignored in considering the effect of ionospheric delay in atmospheric delay, and only the tropospheric delay in the signal propagation process needs to be considered. The Sentinel-1 series satellite opens a new era of free data accessibility [
32]. The measurement data of the Sentinel-1 satellite can be downloaded freely from European Space Agency (ESA) data (
https://scihub.copernicus.eu/dhus/#/home, accessed on 24 September 2019) by user registration and login; the orbit file, with measurement data, can be downloaded freely in the other website of the ESA (
https://scihub.copernicus.eu/gnss/#/home, accessed on 27 October 2019) via different user information. The Sentinel-1 satellite has two-track data, three-track data, and four-track data, the two-track data is generally used in different InSAR data, i.e., D-InSAR (differential interferometric synthetic aperture radar).
Due to the short wavelength and high frequency of the C-band, the effect of ionospheric delay is very small and the effect of atmospheric delay is mainly that of the tropospheric delay in the signal propagation process.
Usually, InSAR data processing is used to extract DEM with high accuracy and D-InSAR data processing is used to obtain the surface deformation displacement caused by an earthquake, a volcano, or city surface subsidence, and to carry out interference stake time series analyses, such as persistent scatters (PS) and Small baseline subsets (SBAS). For simplicity, in this paper, we refer to D-InSAR as InSAR.
In processing the InSAR measurements for crustal deformation displacement, co-polarization is usually used, because its penetrability is better than that of cross-polarization. Sngle look complex (SLC) radar image data, before and after the interested event, is taken as input data (e.g., Sentinel-1A VV polarization image data). One image is used as the master image and another is used as the slave image. The quality of the interference image pairs is evaluated by a baseline estimation and multi-look processing by azimuth resolution and range resolution to suppress the speckle noise of the SAR images. Then, an interferogram could be obtained, based on STRM DEM [
31], in order to reduce the influence of terrain error. This paper adopted a DEM of 90 m × 90 m resolution. The phase unwrapping is carried out by the Goldstein filtering method, using the Delaunay minimum cost flow theory to solve the problem of phase ambiguity. The corresponding precision orbit is used to adjust the orbit and re-flatten it by selected ground control points (GCP). Finally, the line-of-sight (LOS) co-seismic deformation displacement in the WGS84 coordinate system is obtained via geocoding, and the initial or raw displacement field before atmospheric delay correction is obtained.
We chose the repeat track data of the Sentinel-1A satellite for a 12-day interval period before and after the Changning M6.0 earthquake in the interest area, i.e., SLC image data for 9 June 2019 (day of year, DOY 160) and 21 June 2019 (DOY 172) were downloaded, respectively; we obtained the co-seismic displacement field of the area of interest by D-InSAR measurements.
5. Discussion
In order to determine the main factor affecting the precision of the atmospheric delay correction model for the uneven GNSS stations and the weighting method, we built a new atmospheric delay correction model by using the same sources of the GNSS stations’ ZTD and the ERA5 ZTD as those of the model based on the VCE method, and the same weighting method as that of the GACOS model. This was the test model. The results were as follows.
As shown in
Figure 8, the atmospheric delay correction model was no better than the GACOS model, especially in areas where the elevation change in the northeast corner was dramatic and far from the epicenter. There was reason to believe that the weights method accounted for the more dominant factor influencing the accuracy of the atmospheric delay correction model.
The RMS could be used as an indicator to evaluate the accuracy of the atmospheric delay correction models. Therefore, the RMS values of the overall displacement fields in the study area were calculated on the raw displacement and the corrected displacement by different atmospheric delay correction models that were based on the VCE method and the GACOS model, respectively. The results showed that the RMS value was 19.8 mm on the raw displacement field, without atmospheric delay correction, in the study area. The RMS value was 22.7 mm on the corrected displacement field according to the GACOS model, with deterioration of accuracy. The RMS value was 21.4 mm on that corrected displacement according to the test model, which was almost identical to that of the GACOS model. The RMS value of the corrected displacement according to the multi-source data integrated atmospheric delay model based on the VCE method was 17.4 mm, which was a 5.3 mm improvement on that of the corrected displacement by the GACOS model, with accuracy improvement of approximately 25%, and a 2.4 mm improvement on that of the raw displacement, with accuracy improvement of approximately 25%. These were important improvements for highly accurate deformation monitoring.
Generally, the crustal deformation was small in most the regions, except in circumstances of sudden change. For example, near the earthquake center there was large displacement. In order to reduce the effect of the larger deformation by the earthquake, we tested the accuracy of the regions far from the epicenter. The results showed that the RMS was 21.4 mm, 24.4 mm, and 15.1 mm for the raw displacement, the displacement that was corrected via the GACOS model, and the displacement according to the new model, respectively. The results showed that the accuracy of the model based on the VCE method was better by providing 6 mm improvement compared with that of the raw displacement and provided 9 mm improvement in comparison with the displacement corrected by the GACOS model.
Because of the rich vegetation and complex topographic features of the area of interest, the raw phase measurement has low coherence, and other noise masked some of the accuracy of the atmospheric delay correction model. Accordingly, while the precision of the corrected result according to the integrated atmospheric delay correction model based on the VCE method was improved, the improvement was limited.
Another factor to verify the quality of the atmospheric delay correction model is the correlation between the raw displacement and the atmospheric delay correction [
29]. The InSAR phase measurement is highly correlated with the atmospheric delay correction. This means that the atmospheric delay model is good when the atmospheric effect can be most reduced from the raw displacement. Therefore, the higher the correlation coefficient between the raw displacement and the atmospheric delay correction model, the better the accuracy of the atmospheric delay correction model. Generally, the calculation of the correlation coefficient is based on the degree of linear correlation between the atmospheric delay
and the displacement
, which can be expressed as
where
is the slope and
is the intercept.
The correlation coefficient between the atmospheric delay and the raw displacement can be expressed as follows:
where
is the correlation coefficient between
and
,
is the covariance of
and
, and
and
represent the variance of
and
, respectively. In this study,
and
represent the atmospheric delay correction model and the raw displacement, respectively.
is the atmospheric delay correction of each pixel in the study area,
is the average of the atmospheric delay correction for all pixels,
is the raw displacement of each pixel, and
is the average of the raw displacement for all pixels.
The correlation coefficients between the raw displacement and the atmospheric delay correction model based on VCE, as well as the coefficient between the raw displacement and the GACOS model, respectively, are shown in
Table 1. It can be seen that the correlation coefficient between the displacement and the model based on the VCE method was higher than that between the displacement and the GACOS model. The former improved the correlation coefficient by approximately 25%, compared with the latter.
In order to intuitively express the correction effect, the displacement of the test line was analyzed for raw displacement and corrected displacement by different atmospheric delay correction models. The results are shown in
Figure 9.
Figure 9a shows the raw deformation that was not corrected by the atmospheric delay correction model. Based on the test line in
Figure 9a, we extracted the raw displacement and corrected it in different models, as expressed in
Figure 9b.
Figure 9b showed the displacement of the test line in
Figure 9a. In
Figure 9b, the blue area represents raw displacement, the green area represents the corrected displacement by the GACOS model, and the red area represents the corrected displacement according to the model based on the VCE method. Based on the principle that the smaller the displacement in a region far from the epicenter, the better the correction effect, we concluded that the red curve, representing the correction model, was the best, and that its curve contour shape but not its value was almost the same as the blue one, i.e., the curve trend of the corrected test line displacement based on the VCE method was similar to the raw displacement, and the smaller displacement than the raw one indicated that there was a systematic constant difference between the displacement corrected by the atmospheric delay correction model based on the VCE method and the raw displacement. Moreover, the corrected displacement according to the model based on the VCE method was lower. The existence of systematic error between the corrected displacement of the VCE model and the raw displacement could be the noise factors, such as rich vegetation in the region.
Given the high precision deformation monitoring capability of the GNSS stations, we considered using a GNSS station to verify the external accuracy of the results achieved. Because the GNSS station was sparse in the area of interest, we chose and performed the deformation of the SCJU station that was close to the study area by the different models; then, we compared it with the solution of the GAMIT/GLOBK. The derived deformations against the SCJU station from the GAMIT/GLOBK solution, the correction model base on the VCE, and the GACOS model, were 0.016, 0.0377, and 0.0695, respectively. This indicated that the result of the model based on the VCE method was better than that of the GACOS model, due to the smaller deformation difference to the GAMIT/GLOBK. As there was still a distance of about 10 km between the considered pixel and the SCJU station, this external accuracy was used as the reference.
6. Conclusions
This paper’s proposed atmospheric delay correction model is better than the GACOS model in separating the elevation-related stratified components from the water vapor components, based on the VCE method. It can be applied to various geographic situations for crustal deformation monitoring. For example, the selected Changning region has large terrain fluctuations and large changes in temperature and relative humidity. Therefore, its atmospheric delays, and especially its water vapor variations, are very complex. The experimental results showed that the model based on the VCE method is better than the GACOS model and that it is also suitable for flat terrain conditions.
In this paper, the denser GNSS data for the study area and the newer meteorological weather forecasting model, ERA5, were adopted. Even if the degree of granularity of the ERA5 37-layer pressure level is poorer than that of the GACOS model’s ECMWF 137-layer spatial resolution, the accuracy of the new model, based on the VCE method, is much higher than that of the GACOS model. the model’s corrected seismic deformation displacement RMS, based on the VCE method, decreased by 5 mm and 2.4 mm, compared with the GACOS model and the raw displacement. This is important for high-precision crustal deformation monitoring. Given the large deformations due to earthquakes, we tested the corresponding RMS in an area far from the epicenter, the precision improvement of the corrected deformation by the new model was 9.3 mm and 6.3 mm, compared with that of the GACOS model and the raw displacement.
The results indicate that with variance component estimation weighting it is feasible to build the multi-source data integrated atmospheric delay model, and the correction effect of the model based on the VCE method is better than that of the GACOS model. Given the more refined meteorological numerical information that can be obtaine, with an auxiliary to denser GNSS stations and other source data, it is possible to obtain a more accurate atmospheric delay correction model for InSAR by the variance component estimation (VCE) method.
In addition, the vegetation cover in the Changning area was very rich in the study period. For the C-band Sentinel satellite, the penetration of vegetation was not as good as the longer wavelength, such as L-band; thus, there was noise in the interferogram. However, the resulting analysis showed that the multi-source data integration atmospheric delay correction model based on the VCE method had higher precision than the GACOS model. The fault was consistent with the results of the analysis by Yu et al. [
48] of the Changning M6.0 earthquake. The strike of the fault was northwest to southeast and the regional characteristics of the two sides of the fault were obviously different. The uplift of the southeastern side was up to 6 cm and the subsidence of the northwestern side was up to 8 cm. Compared with the GACOS model, the major improvement of the multi-source data integrated atmospheric delay correction model based on VCE was in the northeastern area of the study area, which had the largest elevation change. This indicates that accurate atmospheric delay correction is good for the correction effect of the large topographic elevation change.
We recommend the integration of the different multi-sources by the VCE method to obtain a more accurate atmospheric delay correction model to monitor crustal deformation. With the establishment of the China Seismic Science Experimental Site, the GNSS network in Sichuan and Yunnan will become more dense. This paper provides ideas and methods for the future integration of multi-source data for InSAR atmospheric delay correction. Such ideas and methods could be used as a reference for conducting atmospheric delay correction scientific research in different regions and different terrain conditions.