Next Article in Journal
Cooling Effects of Urban Vegetation: The Role of Golf Courses
Previous Article in Journal
Identification of the Spring Green-Up Date Derived from Satellite-Based Vegetation Index over a Heterogeneous Ecoregion
Previous Article in Special Issue
Bayesian Estimation of Land Deformation Combining Persistent and Distributed Scatterers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

InSAR Atmospheric Delay Correction Model Integrated from Multi-Source Data Based on VCE

1
Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China
2
School of Astronomy and Space Science, University of Chinese Academy of Sciences, Beijing 100049, China
3
Institute of Earthquake Forecasting, CEA, Beijing 100036, China
4
Shanghai Key Laboratory of Space Navigation and Positioning Techniques, Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(17), 4329; https://doi.org/10.3390/rs14174329
Submission received: 4 July 2022 / Revised: 20 August 2022 / Accepted: 26 August 2022 / Published: 1 September 2022

Abstract

:
With the rapid development of interferometric synthetic aperture radar (InSAR) measurement technology, its measurement accuracy requirements are increasing. Atmospheric delay errors must be corrected, especially in the case of crustal deformation monitoring, the 20% variation of tropospheric water vapor among InSAR pairs generally produces range from 10 cm to 14 cm deformation errors. Such errors can be of the same magnitude as the annual changes in crustal deformation, or even greater, masking crustal deformation information and seriously affecting the results of crustal deformation monitoring. Therefore, in order to obtain a more accurate InSAR atmospheric delay correction model, this paper calculated and integrated atmospheric delays that were estimated by different sources, including the 37 pressure levels of the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF)) numerical weather prediction model, ECMWF Reanalysis v5 (ERA5), and Global Navigation Satellite System (GNSS) measurement data from the crustal movement observation network of China, based on the variance component estimation (VCE) weighting method. The results showed that the integrated model, based on the VCE method, is better than the generic atmospheric correction online service (GACOS) model for InSAR measuring of crustal deformation. The precision in monitoring crustal deformations was improved by approximately 5 mm, the correlation coefficient of atmospheric delay errors and crustal deformations improved from 0.287 to 0.347, and accuracy improved by approximately 25%. However, the improvement in accuracy was limited because of system error decoherence that was induced by atmospheric noise caused by abundant vegetation or snow cover. Therefore, in order to achieve more accurate results, we recommend the adoption of the multi-source integrated atmospheric delay correction model, based on the VCE method, for InSAR high-precision measuring of crustal deformation and seismic activities.

Graphical Abstract

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.

3. Methodology

3.1. Numerical Meteorological Forecasting Atmospheric Delay Calculation Model

In 2016, the European Centre for Medium-Range Weather Forecasts (ECMWF) [33] (http://www.ecmwf.int, accessed on 8 May 2021) released a revised version of the fifth generation of the global atmospheric numerical forecast reanalysis (ECMWF Reanalysis v5, ERA5); the data has been updated since 2019 [34]. This information provides 137 model levels of information on atmospheric, terrestrial, and oceanic climate changes within 80 km of height and 37 pressure levels of the relevant atmospheric data, stratified by hourly pressure values from 1 hPa to 1000 hPa on an 0.25° × 0.25° grid. The latter information is freely available. Therefore, this paper applied the pressure, temperature, specific humidity, potential volume, and relative humidity information provided from the 37 pressure levels.
There are two steps in ERA5 atmospheric zenith total delay (ZTD) calculations (subsequently referred to as ERA5 ZTDs). First, the ERA5 ZTDs are calculated for each grid corresponding to the 37 pressure levels’ atmospheric data, by integral methods, as follows [33]:
Z T D int _ s f c = 10 6 h s f c h t o p N d h
In addition,
N = k 1 P e T + k 2 e T + k 3 e T 2
e = s h × P / 0.622
k 1 = 77.604 K / Pa ;   k 2 = 64.79 K / Pa ;   k 3 = 377600.0 K / Pa
where Z T D int _ s f c denotes the atmospheric delay from the bottom to the top of the 37 pressure levels’ meteorological values; h s f c and h t o p denote the bottom and top heights of the ERA5 meteorological data, respectively; N is the total refractive index; P denotes the atmospheric pressure in hPa; e denotes the water vapor pressure in hPa; s h denotes the specific humidity; T denotes the temperature of the nth level at the corresponding grid point; the unit is K; d h = h n 1 h n ; and h n , the potential height, indicates the height of the nth level on the grid point, and is obtained by the potential divided by the local acceleration of gravity named, the unit is km; n decreases with the increase in height.
After deriving the ZTD value by the integration method from the 37 pressure levels’ meteorological data, at approximately 47 km in height, the zenith delay above the top level must be added in order to make it comparable with the ZTD of GNSS stations, based on GNSS’s high-precision data processing solution [33]. However, there is no meteorological data above the top level of approximately 47 km in height. Fortunately, the effect of the wet delay is almost negligible, due to the reduced water vapor content above the top level [35]. The ZTD value above the top level could be calculated by the Saastamoinen delay model, as in the following equation [36]:
Z T D t o p = 0.002277 × P t o p + 0.05 + 1255 T t o p + 273.15 e t o p f φ , H
In addition,
e t o p = r h × 6.11 × 10 7.5 T t o p T t o p + 273.15
f φ , H = 1 0.00266 cos 2 φ 0.00028 H
where P t o p denotes the integrated top pressure level in hPa, T t o p denotes the integrated top level Celsius temperature in °C, e t o p denotes the integrated top level water vapor pressure; r h denotes the integrated top level relative humidity (0~1); φ is the geocentric geodetic latitude of the station or grid point in rad; and H is the grid height in km above the sea level.
By adding the integrated atmospheric delay and the corresponding grid atmospheric delay obtained via the Saastamoinen model, the atmospheric zenith delay derived by the numerical weather prediction model can be expressed as follows [33]:
Z T D = Z T D int _ s f c + Z T D t o p
This Z T D value is the atmospheric delay at the grid position with the pressure of 1000 hPa. In fact, due to the existence of terrain undulations, the obtained atmospheric delay is not necessarily the surface elevation observed by InSAR measurements. Therefore, it is necessary to convert the ZTD value to the surface elevation-by-elevation correction, and the conversion formula used is the following [37,38]:
P = P 0 1 + 8.419 × 10 5 H 0 H P 0 0.190284 5.255303
T = T 0 6.5 H H 0
where P is the surface air pressure and the unit is hPa; P 0 , T 0 , H 0 are the air pressure, temperature, and elevation at the known points of the corresponding grid in hPa, K, and km, respectively; H is the elevation of the surface; T is the surface temperature; and the unit is K.

3.2. GNSS Data Estimated Atmospheric Delay Model

Using GAMIT/GLOBK software, GNSS ZTDs with high accuracy could be obtained [39]. This paper applied the continuously observed GNSS stations of the crustal movement observation network of China (i.e., the land state network, CMONC) and the IGS stations around the study area to obtain absolute ZTD values, for GNSS network stations with a maximum baseline of more than 500 km. All models, including receiver antenna files, antenna phase center models, earth orientation parameters, moon ephemeris tables, sun ephemeris tables, and ocean tide models, were updated. ITRF2014 provided the prior coordinates, with constraints of 0.1 m for CMONC stations and 0.01 m for IGS stations. The GMF mapping function was used for the tropospheric delay model and the zenith wet delays, as the unknown parameters, were solved once per hour, together with other parameters, such as station coordinates and satellite orbits.
However, ZTDs determined by the GNSS network or by ERA5 reanalysis data are often different from ZTDs that are required at the time of InSAR satellite transit in the spatio-temporal distribution. They must be interpolated and matched. This paper applied Lagrange’s interpolation to obtain the required ZTDs of the InSAR observation area at the time of SAR satellite transit.

3.3. GACOS Model

In 2018, CHEN et al. [29] proposed taking advantage of the 0.125° × 0.125° horizontal resolution global coverage and the 137 levels of vertical resolution data information provided by the 6-h ECMWF numerical meteorological model, as well as the high temporal resolution (5 min) and high accuracy measurements of discrete fixed GPS stations, to form an atmospheric delay correction model for InSAR measurements. The optimum weight between ECMWF and GPS is the lowest cross-validation RMS value of the GPS network stations within 150 km of the decorrelation distance. The test range of the ECMEF/ GPS relative weight is from 0.0 to 10 by steps of 0.1. In [29], the atmospheric delay was separated into the elevation-related components and the turbulent components by an iterative approach. Moreover, the atmospheric delay correction model for InSAR measurements was built by interpolation and exponential coefficients for the considered pixel. This GACOS model is generally applicable in both flat and upland areas. Users can download the atmospheric delay correction model for a certain region at a certain time on the GACOS website (http://www.gacos.net/, accessed on 11 May 2021) and load it into their own study area. The GACOS model is a commonly used method to correct the atmospheric delay for InSAR measurements. In this paper, it is used as the comparative atmospheric delay model.
We calculated the GNSS ZTDs by using the more than 60 GNSS stations on the Qinghai-Tibet Plateau of the Crustal Movement Observation Network of China (CMONC) and eight surrounding IGS stations. Figure 2 shows the distribution of the selected GNSS stations in the CMONC, and the location of the epicenter.

3.4. Multi-Source Data Integrated Atmospheric Delay Modeling Based on VCE

With the weather forecasting model update and the increase in the number of GNSS stations, and most importantly with the need for a more reasonable weighting method, this paper established a new multi-source data integration atmospheric delay correction model for InSAR measurements, based on the variance component estimation (VCE) weighting method [40,41,42,43,44], by GNSS high-precision real ZTD observations and high spatial and temporal resolution ERA5 ZTDs.
In order to integrate the same volume for the different sources, we converted the surface elevation of the ERA5 ZTDs to the corresponding elevation position of the GNSS station by the elevation-correction-of-the-neighborhood-grid-points method and compared them with the ZTDs at the GNSS station. The ERA5 ZTDs and the GNSS ZTDs were discussed in Section 3.1 and Section 3.2.
The integrated ZTD is a sum of the elevation-related stratified component and the turbulence component, in that the stratified component can be expressed as the exponential function. The formula is as follows [29].
Z T D = S + T + ε
In addition,
S = L e β h
where Z T D is the integrated atmospheric delay value, T is the turbulence component delay, S is the stratified component delay value, ε is the remaining residual, L is the stratified delay at sea level, β is the exponential coefficient, and h is the height.
Because the accuracy and spatio-temporal resolution of the ZTD obtained by different techniques are different, the integration from different sources needs to consider the weight assignment of each source to obtain the high-precision measurement result. In data processing, the variance component estimation (VCE) is a widely used weighting method [40,41,42,43,44] that can set an arbitrary initial weight value to provide a pre-adjustment and calculate the variance of observation; then, the new variance estimation weight is calculated to improve the previous weight value. These steps were performed repeatedly until the weight factor was convergent. It was most important to establish the error equations for the VCE method. Suppose there are two technologies; then, the VCE method is as set out below.
The functional model is as follows.
l 1 = B 1 x + Δ 1
l 2 = B 2 x + Δ 2
Then,
D l 1 = D Δ 1 = σ 01 2 P 1 1
D l 2 = D Δ 2 = σ 02 2 P 2 1
where l 1 and l 2 are the observations of two technologies, respectively; Δ 1 and Δ 2 are observation errors of two technologies, respectively; x is the unknown matrix; B 1 and B 2 are the coefficient matrices of the two technologies, respectively; D Δ 2 and D Δ 1 are the pretest variances of the two technologies, respectively; σ 01 2 and σ 02 2 are the variances in the unit weight of the two technologies, respectively; and P 1 and P 2 are the observation weights.
The established error equations are as follows.
v 1 = B 1 x ^ l 1
v 2 = B 2 x ^ l 2
Generally, the initial weight is not appropriate, while the variances in unit weight are identical for the different technologies, as follows.
σ 01 2 σ 02 2
This means that we need to establish the relationship between the residual error and the variance in unit weight. The specific formula can be found in references [40,41,42,43,44,45].
In this study, we supposed, first, that the turbulence components were zero and then we established the error equation by derivation of the exponential function. The error equation was as follows.
v = ln L β h ln S
We calculated the coefficient β and L of the stratified component by the VCE method. The iteration stopped until the β and L were convergent., During the process, the stratified and turbulence components were updated constantly by new β and L and the new stratified component was again weighted by the variance component estimation. These steps were repeated. This iteration stopped until the weight factor was convergent [45,46] and each component was also convergent.
The threshold setting of the convergence condition affects the iteration number and the accuracy of the results. the lower the stop threshold, the greater the number of iterations and the longer the procedure operation time; however, higher accuracy is not necessarily obtained. We chose different thresholds of weight factor change for our experiments to facilitate building the more accurate and appropriate atmospheric delay correction model based on the VCE. We found that the correction result was not more highly accurate when the threshold was lower, but the result was consistent with the accuracy obtained by setting the threshold at 0.001. Moreover, the coefficient convergence condition of 0.001 was sufficient to obtain the accuracy of the atmospheric delay correction values at sub-millimeter levels.
The threshold referred to was the convergence condition weight factor. In order to separate the stratification associated with the elevation and turbulent components, each component was required to converge simultaneously with the iterative process in calculating the exponential function coefficients.
We used the atmospheric delays in the study area on DOY 160 and DOY 172 in 2019 as the basis for calculating the daily atmospheric delays based on the VCE method (Figure 3c,d), and compared those results with the daily results produced by the GACOS model in the same areas at the same time (Figure 3a,b). In Figure 3, the vertical axis and horizontal axis represent the latitudinal and longitudinal ranges of the study areas, respectively, and the color-coded strips represent the calculated ZTD values. It can be seen that the two atmospheric delay models were almost the same at the same DOY times.
We found that the ERA5 weight value relative to the GNSS weight value of each considered pixel did not change significantly and was stable on DOY 160 prior to the earthquake; however, there was a significant change on DOY 172 after the earthquake, at which time the weight values varied abruptly from 0.049 to 0.053 between adjacent pixel points. The sizes of the weights on these two days were significantly different. The weight ratio between ERA5 and GNSS on DOY 172 was larger than the weight ratio before the earthquake, but it was far less than 1.0. This might have been because the ERA5 meteorological model had a certain forecast error, so that its accuracy was lower than the accuracy of the GNSS ZTD. However, the increase in the weight ratio of ERA5 relative to GNSS, as well as the weight ratio between some adjacent pixels, was unstable after the earthquake, which might be related to the co-seismic deformation of GNSS stations that led to an increase in the atmospheric delay error and a decrease in their weights.
The difference between the atmospheric delay correction model on DOY 172 and on DOY 160, obtained by multi-source data integration based on the estimated variance components and the corresponding correction by the GACOS model, is provided in Figure 4 (top). The corresponding pixel point ZTD corrected difference between the two models is provided in Figure 4 (bottom). It can be seen that the maximum difference in value between the two models was about 15 mm, which occurred in the northeast part of the study area, which had large changes in topography elevation and was far from the epicenter. The big differences were perhaps because the GNSS stations were not sufficiently dense in places with large terrain changes, leading to ZTD grid calculation errors. The other error caused by the large variations in topography may be due to the fact that the ERA5 meteorological data used in this paper included a 0.25° × 0.25° grid, while the ECMWF meteorological data used in the GACOS model included a 0.125° × 0.125° grid.

4. Results of the Atmospheric Delay Correction Model Applied to Crustal Deformation Monitoring

In order to verify that the precision of the deformation can be improved by the atmospheric delay model, based on the VCE method, the Sentinel-1A satellite data were used to obtain the raw crustal deformation information for the area of interest. The raw co-seismic deformation displacement is shown in Figure 5, which depicts an obvious regional variability; the deformation near the epicenter was drastic. It is shown as striped.
The InSAR measurement provides the line-of-sight (LOS) direction displacement of the satellite in the study area, so that the atmospheric delay correction model needs to be mapped to the LOS of the satellite by a trigonometric function on the mean elevation angle at the moment of the satellite’s transmitting signal. Then, the corrected displacement of the raw displacement field minus the multi-source data integrated atmospheric delay correction model, based on the VCE method, and the GACOS model are provided in Figure 6 and Figure 7, respectively. Generally, it is believed that the crustal deformation in the area far from the epicenter should be small or even nil. Therefore, we emphasized the eastern part of the image, which is far from the epicenter, and compared it with the GACOS model.
As shown in Figure 6 and Figure 7, the results were consistent in most regions; however, there were some differences in certain regions. On the one hand, this might be because of the uneven distribution of the GNSS stations that were used in the GACOS model, which were mainly distributed in Europe. Few GNSS stations in China could affect the accuracy of the GACOS atmospheric delay model. On the other hand, GACOS model weighting by a test step size of 0.1 could affect its accuracy in integrating the meteorological forecasting model and the GNSS data. Although this weighting method is much improved, compared with the previous equal-weight integration method, the choice of its step size was not a good method for optimal weights [47].

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 X and the displacement Y , which can be expressed as
Y = a X + b
where a is the slope and b is the intercept.
The correlation coefficient between the atmospheric delay and the raw displacement can be expressed as follows:
ρ = cov X , Y var ( X ) var ( Y ) = ( x i x ¯ ) ( y i y ¯ ) x i x ¯ 2 ( y i y ¯ ) 2
where ρ is the correlation coefficient between X and Y , cov X , Y is the covariance of X and Y , and var ( X ) and var ( Y ) represent the variance of X and Y , respectively. In this study, X and Y represent the atmospheric delay correction model and the raw displacement, respectively. x i is the atmospheric delay correction of each pixel in the study area, x ¯ is the average of the atmospheric delay correction for all pixels, y i is the raw displacement of each pixel, and y ¯ 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.

7. Summary

The atmospheric delay model for InSAR was established based on the VCE method, by integrating GNSS high-precision observation data and the high spatial—temporal resolution meteorological numerical forecasting model. ERA5. The model’s reliability and advantages were evaluated by comparing the correction results with those of the GACOS model. We found that the VCE method is an effective way to improve the accuracy of crustal deformation monitoring in the multi-source data integrated atmospheric delay correction model.

Supplementary Materials

The Sentinel-1A SAR data is available at https://scihub.copernicus.eu/dhus/#/home, accessed on 24 September 2019; The GACOS model can be found in http://www.gacos.net/, accessed on 11 May 2021; The orbit file can be downloaded in https://scihub.copernicus.eu/gnss/#/home, accessed on 27 October 2019.

Author Contributions

X.W. proposed the conceptualization for the study, supervised the progress of the study, provided advice on issues that arose, and writing—review and editing; X.W. was also responsible for project administration and funding acquisition; X.L. established the study’s methods, analyzed the results, and wrote the manuscript; X.L. was also responsible for validation, formal analysis, investigation, writing—original draft preparation; Y.C. provided suggestions for the model, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 11973073; the National Key Research and Development Program of China, grant number 2016YFB0501405; the Basic project of Ministry of Science and Technology of China, grant number 2015FY310200; and the Shanghai Key Laboratory of Space Navigation and Position Techniques, grant number 06DZ22101.

Data Availability Statement

Data used to support the findings of this study are included within the article.

Acknowledgments

We want to thank the European Space Agency (ESA) for providing free Sentinel-1A SAR data, referred to in the Supplementary Materials.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, Z.B. Research and Application of InSAR Atmospheric Delay Correction Based on Spatial Interpolation Algorithm. Ph.D. Thesis, Anhui University of Science and Technology, Anhui, China, 2018. [Google Scholar]
  2. Zebker, H.A.; Rosen, P.A.; Hensley, S. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. J. Geophys. Res. 1997, 102, 7547–7563. [Google Scholar] [CrossRef]
  3. Massonnet, D.; Feigl, K.L.; Rossi, M.; Adragna, F. Radar interferometry mapping of deformation in the year after the Landers earthquake. Nature 1994, 369, 227–230. [Google Scholar] [CrossRef]
  4. Lanari, R.; Lundgren, P.; Manzo, M.; Casu, F. Satellite Radar Interferometry Time Series Analysis of Surface Deformation for Los Angeles, California. Geophys. Res. Lett. 2004, 31, 613–620. [Google Scholar] [CrossRef]
  5. Ferretti, A.; Prati, C.; Rocca, F. Permanent Scatters in SAR Interferometry. IEEE Trans. Geosci. Remote Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
  6. Beradino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef]
  7. Li, Z.W.; Ding, X.L.; Liu, G.X. Modeling atmospheric effects on lnSAR with meteorological and continuous GPS observations, algorithms and some test results. J. Atmos. Sol.-Terr. Phys. 2004, 66, 907–917. [Google Scholar] [CrossRef]
  8. Song, X.G.; Li, D.R.; Shan, X.J.; Liao, M.S.; Cheng, L. Correction of atmospheric effect in repeat-pass InSAR measurements based on GPS and atmospheric transport model. Chin. J. Geophys. 2009, 52, 1156–1164. [Google Scholar]
  9. Williams, S.; Bock, Y.; Fang, P. Integrated satellite interferometry: Tropospheric noise, GPS estimates and implications for interferometric synthetic aperture radar products. J. Geophys. Res. 1998, 103, 27051–27067. [Google Scholar] [CrossRef]
  10. Treuhaft, R.N.; Lanyi, G.E. The effect of the dynamic wet troposphere on radio interferometric measurements. Radio Sci. 1987, 22, 251–265. [Google Scholar] [CrossRef]
  11. Li, Z.W.; Ding, X.L.; Huang, C.; Zheng, D.W. Modeling of atmospheric effects on InSAR measurements by incorporating terrain elevation information. J. Atmos. Sol.-Terr. Phys. 2006, 68, 1189–1194. [Google Scholar] [CrossRef]
  12. Li, Z.H.; Fielding, E.J.; Cross, P.; Muller, J.P. Interferometric synthetic aperture radar atmospheric correction: GPS topography-dependent turbulence model. J. Geophys. Res. 2006, 111, B02404. [Google Scholar] [CrossRef]
  13. Luo, H.B.; He, X.F. A spatial interpolation for GPS-derived atmospheric delay corrections to InSAR. J. Geod. Geodyn. 2007, 27, 35–38. [Google Scholar]
  14. Jarlemark, P.O.J.; Emardson, T.R. Strategies for spatial and temporal extrapolation and interpolation of wet delay. J. Geod. 1998, 72, 350–355. [Google Scholar] [CrossRef]
  15. Jolivet, R.; Grandin, R.; Lasserre, C.; Doin, M.P.; Peltzer, G. Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. Geophys. Res. Lett. 2011, 38, L17311. [Google Scholar] [CrossRef]
  16. Chen, Y.; Penna, N.T.; Li, Z.H. Generation of real-time mode high-resolution water vapor filed from GPS observations. J. Geophys. Res. Atmos. 2017, 122, 2008–2025. [Google Scholar]
  17. Garthwaite, M.C.; Wang, H.; Wright, T.J. Broadscale interseismic deformation and fault slip rates in the central Tibetan Plateau observed using InSAR. J. Geophys. Res. 2013, 118, 5071–5083. [Google Scholar] [CrossRef]
  18. Wang, H.; Peng, J.H. Test on InSAR atmospheric delay correction using ECMWF model. J. Guangdong Univ. Technol. 2014, 31, 74–77. [Google Scholar]
  19. Li, Z.W.; Ding, X.L.; Huang, C.; Zou, Z.R. Atmospheric effects on repeat-pass InSAR measurements over Shanghai region. J. Atmos. Sol.-Terr. Phys. 2007, 69, 1344–1356. [Google Scholar] [CrossRef]
  20. Li, Z.H.; Muller, J.P.; Cross, P.; Fielding, E.J. Interferometric synthetic aperture radar(lnSAR) atmospheric correction: GPS, Moderate Resolution Imaging Spectroradiometer (MODIS), and InSAR integration. J. Geophys. Res. Solid Earth 2005, 110, B30410. [Google Scholar]
  21. Chen, M.; Chen, J.P.; Hu, C.W. Numerical Weather Models Applied in GNSS Precise Point Positioning. In Proceedings of the China Satellite Navigation Academic Conference, Changsha, China, 18–20 May 2016. [Google Scholar]
  22. Li, Z.H. Correction of Atmospheric Water Vapour Effects on Repeat-Pass SAR Interferometry Using GPS, MODIS and MERIS Data. Ph.D. Thesis, University College London, London, UK, 2006. [Google Scholar]
  23. Li, Z.H.; Fielding, E.J.; Cross, P.; Muller, J.P. Interferometric synthetic aperture radar atmospheric correction: Medium Resolution Imaging Spectrometer and Advanced Synthetic Aperture Radar integration. Geophys. Res. Lett. 2006, 33, L06816. [Google Scholar] [CrossRef]
  24. Xu, W.B.; Li, Z.W.; Ding, X.L.; Feng, G.C. Correcting atmospheric effects in ASAR interferogram with MERIS integrated water vapor data. Chin. J. Geophys. 2010, 53, 1073–1084. [Google Scholar]
  25. Li, X.F.; Li, Y.; Zeng, Q.M.; Zhao, Y.H. Correction of atmospheric effects on repeat-pass interferometric SAR using MERIS and ASAR synchronous data. Acta Sci. Nat. Univ. Pekin. 2009, 45, 1012–1018. [Google Scholar]
  26. Hu, J.; Li, Z.W.; Zhu, J.J.; Ding, X.L.; Qian, S. Measuring three-dimensional surface displacements from combined InSAR and GPS data based on BFGS method. Chin. J. Geophys. 2013, 56, 117–126. [Google Scholar]
  27. Zhang, X.Q.; Zhou, Y.H. Research on D-InSAR Atmospheric Delay Correction with CORS and MODIS. Bull. Surv. Mapp. 2015, 8, 13–18. [Google Scholar]
  28. Zhou, W.B.; Xu, W.B.; Li, Z.W.; Wang, C.C. Elevation-dependent MERIS water vapor interpolation and its Application to Atmospheric Correction on ASAR Interferogram. Geomat. Inf. Sci. Wuhan Univ. 2012, 37, 963–977. [Google Scholar]
  29. Chen, Y.; Li, Z.H.; Penna, N.T.; Paola, C. Generic atmospheric correction model for interferometric synthetic aperture radar observation. JGR Solid Earth 2018, 123, 8295–9285. [Google Scholar]
  30. Qin, Z.P.; Liu, S.G.; Deng, B.; Li, Z.W.; Sun, W. Multiphase structural features and evolution of southeast Sichuan tectonic belt in China. J. Chengdu Univ. Technol. Sci. Technol. Ed. 2013, 40, 703–711. [Google Scholar]
  31. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Alsdorf, D. The Shuttle Radar Topography Mission. Rev. Geophys. 2007, 45, 2. [Google Scholar] [CrossRef]
  32. Rocca, F. Perspectives of Sentinel-1 for InSAR applications. In Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium (IGARSS 2012), Munich, Germany, 22–27 July 2012. [Google Scholar]
  33. Chen, Q.M.; Song, S.L.; Heise, S.; Liou, Y. Assessment of ZTD derived from ECMWF/NCEP data with GPS ZTD over China. GPS Solut. 2011, 15, 415–425. [Google Scholar] [CrossRef]
  34. Ma, Z.Q.; Chen, Q.M.; Gao, D.Z. Study on accuracy of ZTD and ZWD calculated from ERA-Interim data over China. J. Geod. Geodyn. 2012, 32, 100–104. [Google Scholar]
  35. Meng, X.G.; Guo, J.J.; Han, Y.Q. Preliminarily assessment of ERA5 reanalysis data. J. Mar. Meteorol. 2018, 38, 91–99. [Google Scholar]
  36. Saastamoinen, J. Atmospheric Correction for the Troposphere and Stratosphere in Radio Ranging Satellites. Use Artif. Satell. Geodesy. 1972, 15, 247–251. [Google Scholar]
  37. Sheng, P.X.; Mao, J.T.; Li, J.G.; Ge, Z.M.; Zhang, A.C.; Sang, J.G.; Pan, N.X.; Zhang, H.S. Atmospheric Physics; Beijing University Press: Beijing, China, 2003. [Google Scholar]
  38. Zhao, J.Y.; Song, S.L.; Zhu, Y.W. Accuracy Assessment of Applying ERA-Interim Reanalysis Data to Calculate Ground-based GPS/PWV over China. Geomat. Inf. Sci. Wuhan Univ. 2014, 39, 935–939, 1008. [Google Scholar]
  39. Jia, W.; Pan, Y.J.; He, M.L.; Qin, Y.H. Calculating Analysis of Troposphere Delay Quantity Based on GAMIT. Geospat. Inf. 2013, 11, 46–48. [Google Scholar]
  40. Cui, X.Z.; Yu, Z.C.; Tao, B.Z. Generalized Surveying Adjustment; Wuhan University Press: Wuhan, China, 2009. [Google Scholar]
  41. Nie, Y.; Shen, Y.; Pail, R.; Chen, Q. Efficient variance component estimation for large-scale least-squares problems in satellite geodesy. J. Geod. 2022, 96, 1–15. [Google Scholar] [CrossRef]
  42. Xu, P.; Shen, Y.; Fukuda, Y.; Liu, Y. Variance component estimation in linear inverse ill-posed models. J. Geod. 2006, 80, 69–81. [Google Scholar] [CrossRef]
  43. Xu, P.; Liu, Y.; Shen, Y.; Fukuda, Y. Estimability analysis of variance and covariance components. J. Geod. 2007, 81, 593–602. [Google Scholar] [CrossRef]
  44. Zhao, J.; Guo, J.F. A universal Formula of Variance Component Estimation. Geomat. Inf. Sci. Wuhan Univ. 2013, 38, 580–583, 588. [Google Scholar]
  45. Shao, F. Research on high precision SLR and ITS astronomical and Geodetic Applications. Ph.D. Thesis, Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai, China, 2019; pp. 75–79. [Google Scholar]
  46. Wu, X.Q. The Method of Weight Factor for Variance Component Estimation. Geomat. Inf. Sci. Wuhan Univ. 1989, 14, 8–15. [Google Scholar]
  47. Löfgren, J. Tropospheric correction for InSAR using interpolated ECMWF data and GPS Zenith Total Delay from the Southern California Integrated GPS Network. In Proceedings of the 2010 IEEE International Geoscience and Remote Sensing Symposium (IGARSS 2010), Honolulu, HI, USA, 25–30 July 2010. [Google Scholar]
  48. Yu, X.W.; Xue, D.J.; Wang, H.F. Deformation field extraction and simulation of Changning earthquake based on D-InSAR. Bull. Surv. Mapp. 2020, 5, 59–63. [Google Scholar]
Figure 1. The range of the area of interest and the topography from SRTM3 v4.1 DEM [31]. The red rectangular represents the area of interest, the yellow pentagram represents the earthquake location as identified by the China Earthquake Networks Center. The color-coded strip represents the elevation of the topography.
Figure 1. The range of the area of interest and the topography from SRTM3 v4.1 DEM [31]. The red rectangular represents the area of interest, the yellow pentagram represents the earthquake location as identified by the China Earthquake Networks Center. The color-coded strip represents the elevation of the topography.
Remotesensing 14 04329 g001
Figure 2. Distribution of selected GNSS stations in the crustal movement observation network of China near the area of interest. The green triangle represents the selected GNSS stations of the Crustal Movement Observation Network of China (CMONC). The yellow pentagram represents the location of the epicenter by the China earthquake networks center (CENC).
Figure 2. Distribution of selected GNSS stations in the crustal movement observation network of China near the area of interest. The green triangle represents the selected GNSS stations of the Crustal Movement Observation Network of China (CMONC). The yellow pentagram represents the location of the epicenter by the China earthquake networks center (CENC).
Remotesensing 14 04329 g002
Figure 3. Daily atmospheric delay obtained from different models: (a) the atmospheric delay of the InSAR signals on DOY 172 according to the GACOS model; (b) the atmospheric delay of the InSAR signals on DOY 160 according to the GACOS model; (c) the atmospheric delay of the InSAR signals on DOY 172 according to the variance component weighting model; and (d) the atmospheric delay of the InSAR signals on DOY 160 according to the variance component weighting model. The color−coded strips represent the calculated ZTD; the vertical axis and horizontal axis represent the north latitudinal range and the east longitudinal range of the study area, respectively.
Figure 3. Daily atmospheric delay obtained from different models: (a) the atmospheric delay of the InSAR signals on DOY 172 according to the GACOS model; (b) the atmospheric delay of the InSAR signals on DOY 160 according to the GACOS model; (c) the atmospheric delay of the InSAR signals on DOY 172 according to the variance component weighting model; and (d) the atmospheric delay of the InSAR signals on DOY 160 according to the variance component weighting model. The color−coded strips represent the calculated ZTD; the vertical axis and horizontal axis represent the north latitudinal range and the east longitudinal range of the study area, respectively.
Remotesensing 14 04329 g003
Figure 4. Comparison for atmospheric delay correction between the variance component weighting model and the GACOS model: the top shows the difference between the two models in the area of interest; the color−coded strips represent the difference between the atmospheric delay correction model, based on the VCE method, and the GACOS model; the vertical axis and the horizontal axis represent the northern latitudinal range and the eastern longitudinal range of the study area, respectively; at the bottom is the difference between the two atmospheric delay correction models in the corresponding pixels: the horizontal axis represents the pixels extracted by row, and the vertical axis represents the difference values; the unit is meter.
Figure 4. Comparison for atmospheric delay correction between the variance component weighting model and the GACOS model: the top shows the difference between the two models in the area of interest; the color−coded strips represent the difference between the atmospheric delay correction model, based on the VCE method, and the GACOS model; the vertical axis and the horizontal axis represent the northern latitudinal range and the eastern longitudinal range of the study area, respectively; at the bottom is the difference between the two atmospheric delay correction models in the corresponding pixels: the horizontal axis represents the pixels extracted by row, and the vertical axis represents the difference values; the unit is meter.
Remotesensing 14 04329 g004
Figure 5. InSAR co−seismic deformation for the area of interest, corresponding to Figure 1. The red rectangle represents the area of interest; the black pentagram represents the earthquake location according to the China Earthquake Networks Center; the color−coded strip represents the displacement of co-seismic deformation for InSAR.
Figure 5. InSAR co−seismic deformation for the area of interest, corresponding to Figure 1. The red rectangle represents the area of interest; the black pentagram represents the earthquake location according to the China Earthquake Networks Center; the color−coded strip represents the displacement of co-seismic deformation for InSAR.
Remotesensing 14 04329 g005
Figure 6. The corrected deformation displacement according to the atmospheric delay model based on the VCE method: (a) the multi−source data integrated atmospheric delay correction model based on the variance component estimation (VCE) method; (b) the corrected co−seismic deformation displacement according to the model based on the VCE. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Figure 6. The corrected deformation displacement according to the atmospheric delay model based on the VCE method: (a) the multi−source data integrated atmospheric delay correction model based on the variance component estimation (VCE) method; (b) the corrected co−seismic deformation displacement according to the model based on the VCE. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Remotesensing 14 04329 g006
Figure 7. The corrected deformation displacement according to the GACOS model: (a) the multi-source data integrated atmospheric delay correction based on the GACOS model; (b) the corrected co−seismic deformation displacement according to the GACOS model. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Figure 7. The corrected deformation displacement according to the GACOS model: (a) the multi-source data integrated atmospheric delay correction based on the GACOS model; (b) the corrected co−seismic deformation displacement according to the GACOS model. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Remotesensing 14 04329 g007
Figure 8. The corrected deformation displacement by the test model: (a) the atmospheric delay correction according to the test model; (b) the corrected co−seismic deformation displacement according to the test model. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Figure 8. The corrected deformation displacement by the test model: (a) the atmospheric delay correction according to the test model; (b) the corrected co−seismic deformation displacement according to the test model. The color−coded strip represents the atmospheric delay correction model and the displacement in (a,b), respectively.
Remotesensing 14 04329 g008
Figure 9. A comparison of the test line displacement for the before and after correction with different atmospheric delay models. Part (a) shows a raw interference displacement of the study area. We can see the position of the test line; the black pentagram is the epicenter according to the China Earthquake Networks Center; the color−coded strip represents the displacement. Part (b) shows the displacement of the black test line located in (a); the bule curve is the raw displacement, the green curve is the corrected displacement according to the GAOCS model, and the red curve is the corrected displacement according to the model based on the VCE, respectively.
Figure 9. A comparison of the test line displacement for the before and after correction with different atmospheric delay models. Part (a) shows a raw interference displacement of the study area. We can see the position of the test line; the black pentagram is the epicenter according to the China Earthquake Networks Center; the color−coded strip represents the displacement. Part (b) shows the displacement of the black test line located in (a); the bule curve is the raw displacement, the green curve is the corrected displacement according to the GAOCS model, and the red curve is the corrected displacement according to the model based on the VCE, respectively.
Remotesensing 14 04329 g009
Table 1. Performance indicators for the GACOS model and the atmospheric delay correction model based on the VCE correlation coefficients.
Table 1. Performance indicators for the GACOS model and the atmospheric delay correction model based on the VCE correlation coefficients.
ModelLinear Fit between Raw Displacement and Atmospheric DelayCorrelation Coefficient
GACOS modelY = 1.6692X − 12.80420.2870
multi-source data integrated atmospheric delay correction model based on VCE Y = 0.9495X − 6.88840.3465
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, X.; Wang, X.; Chen, Y. InSAR Atmospheric Delay Correction Model Integrated from Multi-Source Data Based on VCE. Remote Sens. 2022, 14, 4329. https://doi.org/10.3390/rs14174329

AMA Style

Li X, Wang X, Chen Y. InSAR Atmospheric Delay Correction Model Integrated from Multi-Source Data Based on VCE. Remote Sensing. 2022; 14(17):4329. https://doi.org/10.3390/rs14174329

Chicago/Turabian Style

Li, Xiaobo, Xiaoya Wang, and Yanling Chen. 2022. "InSAR Atmospheric Delay Correction Model Integrated from Multi-Source Data Based on VCE" Remote Sensing 14, no. 17: 4329. https://doi.org/10.3390/rs14174329

APA Style

Li, X., Wang, X., & Chen, Y. (2022). InSAR Atmospheric Delay Correction Model Integrated from Multi-Source Data Based on VCE. Remote Sensing, 14(17), 4329. https://doi.org/10.3390/rs14174329

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop