On the Assessment GPS-Based WRFDA for InSAR Atmospheric Correction: A Case Study in Pearl River Delta Region of China

: The accuracy and applications of synthetic aperture radar interferometry (InSAR) are severely suppressed by tropospheric error. Numerical Weather Models (NWMs) and GPS-derived tropospheric delays have been widely used to correct the tropospheric error considering their complete spatial coverage or high accuracy. However, few studies focus on the fusion of both NWMs and GPS for the tropospheric error correction. In this study, we used the Weather Research and Forecasting (WRF) to obtain NWMs with a higher spatial-temporal resolution of 3 km and 20 s from both ERAI (79 km and 6 h) and ERA5 (0.25 ◦ and 1 h). After that, we utilized the WRF Data Assimilation (WRFDA) system to assimilate the GPS ZTD into these enhanced NWMs and generate merged NWMs products. The tropospheric correction effectiveness from different NWMs products was evaluated in a case in the Pearl River Delta region of China. The results showed that all the NWMs products could correct the stratiﬁed component in the interferogram but could not mitigate the turbulence well, even after improving the spatial-temporal resolution. As for the trend component, the merged NWMs products showed obvious superiority over other products. From the statistics perspective, the stdev of the interferogram decreased further over 20% by the merged NWMs products than other products when using both ERAI and ERA5, indicating the signiﬁcant effectiveness of GPS ZTD assimilation.


Introduction
Synthetic aperture radar interferometry (InSAR) has demonstrated a powerful measurement capability in a wide range of applications such as earthquake deformation [1], land subsidence [2], and volcanic activity monitoring [3,4]. However, similar to other spatial geodetic techniques, InSAR measurements are inevitably affected by the atmosphere, which is still an intractable limitation to high-accuracy applications and needs careful corrections [5][6][7].
There are mainly two methods to correct InSAR atmospheric delay [5], namely (1) by using the interferogram itself or (2) by using external atmospheric products (referred to EAP hereafter). The former method generally assumes atmospheric error as temporal random noise [8] or as an elevation-dependent component [9,10], which either needs a large set of interferograms to separate deformation rate and noise or relies on the assumption of correlations between the deformation and elevation. By contrast, the latter method, which is based on EAP, such as Global Positioning System (GPS) [11], Numerical Weather Models (NWMs) [12][13][14][15][16], and Moderate Resolution Imaging Spectroradiometer (MODIS) [17], can provide direct tropospheric delay correction and is more straightforward.
The temporal-spatial resolution and accuracy of EAP are vital for reliable InSAR atmospheric correction, which is, however, not satisfying for most EAP. For example, MODIS measurements are sensitive to the presence of clouds and need calibration [18]. NWMs product, such as ERA-Interim reanalysis (referred as ERAI hereafter) and operational analysis from the European Centre for Medium-range Weather Forecast (ECMWF), which has been widely used in many works of literature, can guarantee the complete spatial coverage for estimating atmospheric corrections [12,13,19,20]. Nevertheless, most current global NWMs are still insufficient in spatial resolution to capture the small-scale atmospheric water vapor turbulence affecting InSAR measurements [21]. The temporal-spatial resolution of NWMs can be further improved regionally by using tools like Weather Research and Forecasting (WRF) for InSAR atmospheric corrections [22]. However, Kinoshita et al. [23] stated that WRF only improves the spatial and temporal resolution but helps little to the accuracy of parameters that primarily relies on observations. Doin et al. [24] argued that combining mesoscale atmospheric simulations and data assimilation (DA) would be the most promising tool for correcting InSAR atmospheric delays. Yun et al. [25] assimilated meteorological data from the Global Telecommunications System (GTS) into Final Operational Global Analysis (FNL) with the WRF-3DVar data assimilation system and demonstrated the high potential and effectiveness of WRFDA in InSAR atmospheric correction.
The ground-based GPS has advantages of high accuracy, high temporal resolution, low cost, and continuous in all weather conditions for retrieving atmospheric delay or water vapor content [26]. Therefore, it is attractive to combine NWMs and GPS zenith total delay (ZTD) or water vapor content by data assimilation to generate merged NWMs. Yu et al. [13,27] applied an iterative tropospheric decomposition (ITD) method to combine NWMs with ground-based GPS and found that the combination products perform better than NWMs for InSAR corrections. However, this combination was conducted based on the interpolation of relative ZTDs but not data assimilation systems. Mateus et al. [28] carried out experiments in Lisbon city and found that WRF can reproduce the spatial atmospheric trend in the calibrated InSAR. There have also been some studies where InSAR and GPS ZTD were assimilated into NWMs by WRFDA to improve the weather forecast performance [29][30][31][32]. However, there are currently few studies assimilating GPS ZTD into NWMs for the InSAR correction. Therefore, this paper will make a preliminary study to investigate the performance of correcting InSAR atmospheric delay by combining NWMs and GPS ZTD. Three kinds of NWMs, namely, the reanalysis, reanalysis-based WRF without GPS ZTD assimilation, and reanalysis-based WRF with GPS ZTD assimilation, will be compared. The paper is organized as follows: Section 2 introduces the InSAR atmospheric delay correcting method, the used InSAR and reanalysis data as well as the WRF and WRFDA configurations. Section 3 presents the corrected results and Section 4 makes some discussions, followed by the conclusions in Section 5.

Data and Methodology
In this section, the experiment setup and the related data will be introduced first, followed by the methods for generating different correcting products.

Study Area and InSAR Data Processing
A coastal area in the Pearl River Delta region, Guangdong province in the south of China, is selected as the study area, as shown in Figure 1, where atmospheric turbulences are usually stronger and harder to predict in NWMs than inland regions [33]. The SAR images of Sentinel-1A on 11/01/2017 and 23/01/2017 are selected to generate the interferogram. The acquisition time is 10:33 UTC, the perpendicular baseline is −17 m, and the flight direction is ascending. Terrain deformations are negligible considering the temporal baseline of only 12 days. The topography of the study area. The black square covers the blue footprint of Sentinel1-A, and the red triangles represent the used GPS stations.
The interferogram was produced by using the ISCE software [34]. The USGS SRTM 30m digital elevation model (DEM) generated and corrected the topographic phase component. Precise orbits from the Copernicus POD (Precise Orbit Determination) Service were applied to eliminate the orbital error. The interferogram was unwrapped using the integrated correlation and unwrapping (ICU) package [35]. Additionally, the interferogram was multi-looked by 100 × 25 looks in azimuth and range directions, respectively, and sampled to a final resolution of ~300 m.

Atmospheric Delay in InSAR
The phase observation of the interferogram, , can be written as follows: where is the phase caused by the topography, which can be reproduced by using the DEM and removed; denotes the surface deformation component; is the component caused by the satellite orbital error, which is negligible when using the precise orbit ephemeris; represents the atmospheric delay where the ionospheric delay can be ignored for C-band InSAR such as the Sentinel-1 used in our study, and therefore only the tropospheric delay needs to be considered [36]; denotes the un-modeled residual. The tropospheric delay can be divided into dry and wet delay following its physicochemical property or separated into the stratified component and the turbulence based on whether it is elevation-related or not [33]. The tropospheric delay in InSAR ( ) is the difference of delays at two SAR acquisition times and can be expressed as, where the subscript r and s represent reference and secondary interferogram. The tropospheric delay (∆ ), which is also called slant total delay (STD), can be obtained by projecting ZTD onto the oblique path using the cosine function of incidence angle in InSAR [37]. Therefore, the key to acquiring reliable atmospheric delay is to calculate ZTD for each pixel at the acquisition time accurately. The topography of the study area. The black square covers the blue footprint of Sentinel1-A, and the red triangles represent the used GPS stations.
The interferogram was produced by using the ISCE software [34]. The USGS SRTM 30m digital elevation model (DEM) generated and corrected the topographic phase component. Precise orbits from the Copernicus POD (Precise Orbit Determination) Service were applied to eliminate the orbital error. The interferogram was unwrapped using the integrated correlation and unwrapping (ICU) package [35]. Additionally, the interferogram was multi-looked by 100 × 25 looks in azimuth and range directions, respectively, and sampled to a final resolution of~300 m.

Atmospheric Delay in InSAR
The phase observation of the interferogram, φ obs , can be written as follows: where φ topo is the phase caused by the topography, which can be reproduced by using the DEM and removed; φ de f denotes the surface deformation component; φ orb is the component caused by the satellite orbital error, which is negligible when using the precise orbit ephemeris; φ atm represents the atmospheric delay where the ionospheric delay can be ignored for C-band InSAR such as the Sentinel-1 used in our study, and therefore only the tropospheric delay needs to be considered [36]; φ ε denotes the un-modeled residual. The tropospheric delay can be divided into dry and wet delay following its physicochemical property or separated into the stratified component and the turbulence based on whether it is elevation-related or not [33]. The tropospheric delay in InSAR (δL LOS ) is the difference of delays at two SAR acquisition times and can be expressed as, where the subscript r and s represent reference and secondary interferogram. The tropospheric delay (∆L LOS ), which is also called slant total delay (STD), can be obtained by projecting ZTD onto the oblique path using the cosine function of incidence angle α in InSAR [37]. Therefore, the key to acquiring reliable atmospheric delay is to calculate ZTD for each pixel at the acquisition time accurately.

GPS Data Processing
There are 59 GPS stations in the study region, as shown in Figure 1b, where three stations experienced outages during the study period, resulting in 56 stations available in this work. The mean station spacing is about 35 km. The PANDA software [38] was used to process GPS observations in 30 s intervals with the precise point positioning (PPP) method in post-processing mode. Final satellite orbit and clock products from the International GNSS Service (IGS) (https://cddis.nasa.gov/archive/gnss/products/) (access on 8 June 2021) were used. The VMF3 model [39] was applied to estimate the a priori ZHD and the a priori ZWD and provide mapping functions. Residuals for the ZWD were estimated as piece-wise constants every 10 min with a power density of 15 mm h −1/2 . The final GPS ZTDs were retrieved by summing the a priori ZHD, the a priori ZWD, and the estimated ZWD residuals. The tropospheric gradients, including the north-south and the east-west components, were estimated every two hours to account for the azimuthal asymmetry of the atmosphere. Receiver clocks were considered uncorrelated between successive epochs and were estimated every epoch. Station coordinates were estimated as daily constants. The cutoff elevation angle of satellites was set to 7 • , and an elevation-dependent weighting strategy was applied [40]. The precision index of the estimated ZTDs output by the software was used for further designing the WRFDA assimilation weighting strategies.

Reanalysis and Data WRF Simulation
Two kinds of reanalysis data: ERAI and ERA5, both from the ECMWF, will be utilized in this paper. The ERAI dataset has a global coverage (except the poles) with a horizontal resolution of about 79 km and a temporal resolution of 6 h [41]. The ERA5 dataset is the latest reanalysis released by the ECMWF. Compared with its predecessor ERAI, ERA5 improves the temporal and spatial resolution to 1 h and 0.25 • . In addition, ERA5 assimilates recent observations and newly reprocessed datasets, marking a potential improvement in overall quality and detail level [42]. The pressure-level reanalysis products, including the air temperature (T), air pressure (P), geopotential height (H), and specific humidity (Q), nearest the SAR acquisition time were used to calculate ZTD. In addition, the u/v components of wind in each pressure level, 10-m u/v wind components, skin temperature, surface pressure, and other surface parameters, were extracted to serve as the WRF initial and lateral boundary conditions. The Weather Research and Forecasting Model is a mesoscale numerical weather prediction system designed for atmospheric research and operational forecasting applications [43]. The WRF can generate three-dimensional fields for meteorological parameters at a given time, taking NWMs as initial and boundary conditions and allow nests to achieve higher resolution. The WRF in version 4.2 was used in this work. Three two-way nested domains were set with a horizontal resolution of 27, 9, and 3 km, respectively, as shown in Figure 1. The air pressure at the top was set to be 10 hPa, and a total of 49 vertical levels were used in eta vertical coordinate, with the thickness of the top layer and lowest layer of 1000 and 50 m, respectively. The Goddard 4-ice microphysics scheme [44] was used for microphysics. The Kain-Fritsch scheme [45] was applied for the cumulus parameterization on domain d01 and d02, and domain d03 was regarded as convection-permitting [25]. The Dudhia scheme [46] was used for the shortwave radiation, and the Rapid Radiative Transfer Model (RRTM) [47] was used for the longwave radiation. The Yonsei University (YSU) [48] scheme was used for the planetary boundary layer. ERA5 and ERAI were used separately to provide the atmosphere and land surface initial and lateral boundary conditions. The cold start initialization was implemented and reanalysis data from 6:00 to 12:00 were used to cover the SAR acquisition time (10:33 UTC). Time steps were set as 180 s, 60 s, and 20 s for d01, d02, and d03, respectively.
The air temperature (T), air pressure (P), height (H), and specific humidity (Q) are needed to calculate ZTD where T, P, and H are WRF output variables, and Q can be converted from the water vapor mixing ratio (QVAPOR) by, The NCAR command language (NCL) [49] is used to extract parameters from WRF outputs and to make further interpolation and extrapolation to the desired layers.

WRFDA and Configuration
The Weather Research and Forecasting model data assimilation system provides an improved estimate (the analysis) of the atmospheric or oceanic state by combining observations with NWMs. The estimation is based on the iterative minimization of a prescribed cost function as follow, where x a and x b are the analysis and the background (previous forecast), respectively. y 0 denotes the assimilated observations and H is the observation operator. B and R are the background and observational error covariance matrices, respectively. Thus, the analysis is the compromise between the observations and the background based on error statistics. The three-dimensional variational (3DVAR) and four-dimensional variational (4DVAR) data assimilation are available in WRF. The 4DVAR uses the numerical forecast model to mitigate the impact from the time difference between observations and the background. Compared with 4DVAR, however, 3DVAR requires fewer computing resources and has higher speed which has been widely used in many studies and will be utilized in this work. For a more detailed description of the 3DVAR, see Barker et al. [50]. Since the observations are assimilated regardless of the time difference with the background in an assimilation time window for 3DVAR, only GPS ZTDs at 10:30 UTC (nearest the SAR acquisition time 10:33 UTC) were assimilated. The assimilation performed in the inner nest. The generic CV3 background error statistics file was used in this paper [51]. The GPS ZTD uncertainty was derived from the inner precision index output from the GPS data processing. The inner precision index is generally too ideal with a value smaller than 2 mm. However, the nominal ZTD uncertainty estimated in post-processing mode is claimed to be about 4 mm globally by IGS (https://www.igs.org/products/#about) (access on 10 June 2021). In order to set a reasonable ZTD uncertainty in WRFDA, we first estimated a scale factor to convert the inner precision index to ZTD uncertainty as follow, where ξ is the scale factor. σ s is the empirical precision, which is set to 5 mm considering the location of the study area (coastal and low-latitude). σ i is the inner precision index with outliers (larger than 2 mm) removed and n denotes the epoch number. Table 1 presents the statistics of the scale factors estimated on two acquisition days, and we finally set the scale factor to 4 for simplicity in the later sections.

NWMs-Based ZTD Estimation
NWMs have an upper boundary with high altitude, and tropospheric delay exhibits different characters above and below that boundary. The tropospheric delay is composed of both dry and wet delays. However, almost only the dry component is left in the tropospheric delay above that boundary. Therefore, we separate ZTD into two parts, where ZTD nwm and ZTD top are the part under and above the NWMs boundary, respectively. According to Haase et al. [52], the ZTD nwm can be expressed as, where P sur f ace and P top denote the air pressure in Pa at the ground point and NWMs top layer, respectively. R d and R v are the gas constant of dry air and water vapor with values of 287.058 J K −1 kg −1 and 461.5 J K −1 kg −1 , respectively. ε is the ratio of the gas constants with the value of 0.622. g denotes the gravity acceleration in m s −2 , T denotes the air temperature in K, and Q is the specific humidity in kg kg −1 . The remaining terms are constants where k 1 = 0.7760 K Pa −1 , k 2 = 0.704 K Pa −1 , and k 3 = 0.03739 × 10 5 K 2 Pa −1 refer to Haase et al. [52]. The Saastamoinen model [53] is used to calculate ZTD top as, where ϕ is the latitude of the point and h top is the height of the top layer of the NWMs in kilometers.
The NWMs grids are generally not collocated with the interferogram pixels. Here, the NWMs parameters are horizontally interpolated to each pixel in the interferogram. The vertical difference of NWMs-based parameters between the ground point and NWMs layers is compensated by the method described in Zhang et al. [54].

Results and Comparisons
In this section, atmospheric corrections for the study case from three NWMs products, namely, reanalysis (referred to as ERA), reanalysis-based WRF without GPS ZTD assimilation (referred to as WRF), and reanalysis-based WRF with GPS ZTD assimilation (referred to as WRFDA), will be compared at both GPS station locations in Section 3.1 and interferogram pixels in Sections 3.2 and 3.3.

Comparisons at GPS Stations
Differences in the atmospheric delay between two acquisition epochs estimated from GPS ZTD at GPS stations are overlapped on the interferogram, as shown in Figure 2. The datum discrepancy between GPS and interferogram is eliminated by setting the atmospheric delay difference at station DLSH to zero. As mentioned in the previous section, after removing impacts from topography and orbit errors and neglecting deformation, the residual phases in the interferogram are mainly caused by atmospheric delay differences between two SAR acquisitions. From Figure 2, we can find a pretty good agreement between GPS-derived tropospheric delay differences and residual interferogram phases, indicating that GPS can well match the tropospheric delays in InSAR in this case. In addition, the atmosphere generally shows more substantial variations over regions near the coastline than inland regions, resulting in larger delay residuals in Figure 2 Comparison of atmospheric delay differences at 18 GPS stations between three NWMs products and GPS are illustrated in Figure 3. The linear fitting slope value can be used as an indicator to evaluate the performance of each NWMs product, where the closer the slope value to one, the better the NWMs product is consistent with GPS measurement. ERA5 shows better performance than ERAI when using the ERA method, while the situation reverses when using the WRF method. After the assimilation of GPS measurements, the WRFDA product has much better consistency with the GPS product, with slope values of 0.62 and 0.65 when using ERAI and ERA5 as input in WRFDA. Again, ERA5 is superior to ERAI in the WRFDA method in this case. The GPS and three NWMs products are also applied to correct the interferogram. Theoretically, if the EAP products are perfect, only noises remain in the corrected interferogram, with negligible differences of residual phases among different pixels. Comparisons of residual phases at GPS station locations by using different products are presented in Figure 4. The amplitude (difference between maximum and minimum value) and the Comparison of atmospheric delay differences at 18 GPS stations between three NWMs products and GPS are illustrated in Figure 3. The linear fitting slope value can be used as an indicator to evaluate the performance of each NWMs product, where the closer the slope value to one, the better the NWMs product is consistent with GPS measurement. ERA5 shows better performance than ERAI when using the ERA method, while the situation reverses when using the WRF method. After the assimilation of GPS measurements, the WRFDA product has much better consistency with the GPS product, with slope values of 0.62 and 0.65 when using ERAI and ERA5 as input in WRFDA. Again, ERA5 is superior to ERAI in the WRFDA method in this case. Comparison of atmospheric delay differences at 18 GPS stations between three NWMs products and GPS are illustrated in Figure 3. The linear fitting slope value can be used as an indicator to evaluate the performance of each NWMs product, where the closer the slope value to one, the better the NWMs product is consistent with GPS measurement. ERA5 shows better performance than ERAI when using the ERA method, while the situation reverses when using the WRF method. After the assimilation of GPS measurements, the WRFDA product has much better consistency with the GPS product, with slope values of 0.62 and 0.65 when using ERAI and ERA5 as input in WRFDA. Again, ERA5 is superior to ERAI in the WRFDA method in this case. The GPS and three NWMs products are also applied to correct the interferogram. Theoretically, if the EAP products are perfect, only noises remain in the corrected interferogram, with negligible differences of residual phases among different pixels. Comparisons of residual phases at GPS station locations by using different products are presented in Figure 4. The amplitude (difference between maximum and minimum value) and the The GPS and three NWMs products are also applied to correct the interferogram. Theoretically, if the EAP products are perfect, only noises remain in the corrected interferogram, with negligible differences of residual phases among different pixels. Comparisons of residual phases at GPS station locations by using different products are presented in Figure 4. The amplitude (difference between maximum and minimum value) and the standard deviation (stdev) for each method are summarized in Table 2, where smaller amplitude and stdev indicate better correcting performance. From Figure 4 and Table 2, we can easily find that GPS has the best performance among different methods. The amplitude and stdev after correction by the WRF products are larger than those by the ERA products for both ERA-Interim and ERA5, indicating that the WRF products show no improvement than ERA products. After the correction by the WRFDA products, the amplitude and stdev are significantly reduced, illustrating the considerable contributions from GPS. standard deviation (stdev) for each method are summarized in Table 2, where smaller amplitude and stdev indicate better correcting performance. From Figure 4 and Table 2, we can easily find that GPS has the best performance among different methods. The amplitude and stdev after correction by the WRF products are larger than those by the ERA products for both ERA-Interim and ERA5, indicating that the WRF products show no improvement than ERA products. After the correction by the WRFDA products, the amplitude and stdev are significantly reduced, illustrating the considerable contributions from GPS.

Comparisons at Interferogram Pixels for ERAI
The residual phase in the interferogram and atmospheric correction products generated by ERA, WRF, and WRFDA based on ERAI are presented in Figure 5. Atmospheric delays mainly contain two parts. One part is caused by the stratified atmosphere, which is relatively more stable and generally related to the topography, for example, areas indicated by circles in Figure 5a. The other part is caused by the atmosphere transport, including geostrophic flow, general circulation, and turbulence, leading to global and local water vapor variations [55]. For the second part, geostrophic flow and the general circulation determine large-scale air and water vapor transport. At the same time, turbulence usually exhibits considerable irregularity in the vertical and horizontal flow fields and shows small-scale distributions. In this case, the large-scale component induces a noticeable trend in the south-north direction, as shown in Figure 5a, resulting in a variation of nearly

Comparisons at Interferogram Pixels for ERAI
The residual phase in the interferogram and atmospheric correction products generated by ERA, WRF, and WRFDA based on ERAI are presented in Figure 5. Atmospheric delays mainly contain two parts. One part is caused by the stratified atmosphere, which is relatively more stable and generally related to the topography, for example, areas indicated by circles in Figure 5a. The other part is caused by the atmosphere transport, including geostrophic flow, general circulation, and turbulence, leading to global and local water vapor variations [55]. For the second part, geostrophic flow and the general circulation determine large-scale air and water vapor transport. At the same time, turbulence usually exhibits considerable irregularity in the vertical and horizontal flow fields and shows small-scale distributions. In this case, the large-scale component induces a noticeable trend Remote Sens. 2021, 13, 3280 9 of 17 in the south-north direction, as shown in Figure 5a, resulting in a variation of nearly 60 mm in the residual phase differences. The small-scale component, caused by the turbulence, can also be found in the region near the coastline where water vapor variations are usually substantial [24], as shown in Figure 5. For clarity and simplicity of discussion, the small-scale and large-scale components are referred to as turbulence and trend hereafter.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 17 60 mm in the residual phase differences. The small-scale component, caused by the turbulence, can also be found in the region near the coastline where water vapor variations are usually substantial [24], as shown in Figure 5. For clarity and simplicity of discussion, the small-scale and large-scale components are referred to as turbulence and trend hereafter. Three NWMs products can generally provide reliable corrections to the stratified component, which can be found in regions indicated by circles in Figure 5. However, due to the relatively flat terrain and considerable water vapor variations, in this case, the In-SAR atmospheric delay is dominated by turbulence and trend instead of the stratified component. Therefore, the reliability of NWMs in generating turbulence and trend determines the correction performance.
As for the turbulence component, by comparing ERA and WRF products in Figure  5d,e, we can find that WRF method can produce more small-scale details than ERA method due to the improvement of spatial resolution, which can also be found from the differences in Figure 5b. However, these details only make marginal improvements in the correction results, as shown in Figure 5g,h, indicating that the produced details are not consistent with the actual turbulence in the atmosphere. Most turbulences remain in the interferogram after corrections for both ERA and WRF methods. On the other hand, the assimilation of GPS ZTD does not make any difference in the turbulence, as shown in Figure 5c, mainly due to the spatial resolution of GPS stations. Therefore, in this case, all three methods based on ERAI cannot provide reliable turbulence corrections.
Considering the trend component, as mentioned before, the interferogram before corrections shows a trend in the south-north direction. However, the ERA and WRF products show similar spatial distribution with an incorrect trend in the southwest-northeast di- Correcting results with ERAI. (a) is the interferogram before correction, (b) denotes the difference between WRF and ERA products, and (c) denotes the difference between WRF and WRFDA products, (d-f) are atmospheric correction products generated by ERA, WRF, and WRFDA methods, respectively, and (g-i) is the interferogram after atmospheric corrections by (d), (e) and (f), respectively. The cross symbols denote the grid points of ERAI.
Three NWMs products can generally provide reliable corrections to the stratified component, which can be found in regions indicated by circles in Figure 5. However, due to the relatively flat terrain and considerable water vapor variations, in this case, the InSAR atmospheric delay is dominated by turbulence and trend instead of the stratified component. Therefore, the reliability of NWMs in generating turbulence and trend determines the correction performance.
As for the turbulence component, by comparing ERA and WRF products in Figure  5d,e, we can find that WRF method can produce more small-scale details than ERA method due to the improvement of spatial resolution, which can also be found from the differences in Figure 5b. However, these details only make marginal improvements in the correction results, as shown in Figure 5g,h, indicating that the produced details are not consistent with the actual turbulence in the atmosphere. Most turbulences remain in the interferogram after corrections for both ERA and WRF methods. On the other hand, the assimilation of GPS ZTD does not make any difference in the turbulence, as shown in Figure 5c, mainly due to the spatial resolution of GPS stations. Therefore, in this case, all three methods based on ERAI cannot provide reliable turbulence corrections.
Considering the trend component, as mentioned before, the interferogram before corrections shows a trend in the south-north direction. However, the ERA and WRF products show similar spatial distribution with an incorrect trend in the southwest-northeast direction, resulting in a new northwest-southeast trend in the interferograms after corrections. In addition, the WRF product has a larger amplitude for the trend than the ERA product, as can be seen from Figure 5d,e. The regions along the coastline are better corrected with WRF product than with ERA product, while other regions such as the northeast part have more extensive residual phases, illustrating that the larger amplitude in WRF product does not fully agree with the authentic atmosphere. The GPS ZTD can adjust the unrealistic trend in the WRF product, leading to a much better correction by WRFDA product, as shown in Figure 5i, which demonstrates the beneficial contributions from GPS for trend correction. Table 3 summarizes the statistical information, including the amplitude of the atmospheric correction products and the amplitude and stdev of the interferogram after the correction, where smaller amplitude and stdev mean better performance. The pre-corrected interferogram has an amplitude of 59.2 mm and a stdev of 16.25 mm. As shown from Table 3, although the WRF correction product has the largest amplitude of 59.2 mm, which is closest to the amplitude of the pre-corrected interferogram, both the amplitude and stdev after applying the WRF correction product show slight improvement compared to ERA product. Among the three types of products, the WRFDA product provides the best corrections to the interferogram, with amplitude and stdev reduced by 31% and 43%, respectively, mainly due to the assimilation of GPS ZTD. Table 3. Statistical information of correcting results from ERA, WRF, and WRFDA products based on ERAI, including the amplitude of the products and amplitude and stdev of the interferogram after correction. The percentage of reduction compared to the pre-corrected interferogram is given in the brackets.  Figure 6 displays the quantile-quantile (Q-Q) plot of the interferogram with respect to products derived from three methods based on ERAI. If interferogram and correction products have the same distribution, the Q-Q plot should be a straight line. It can be found in Figure 6 that the WRFDA product has the most similar distribution with the interferogram, meaning that WRFDA method can reproduce tropospheric features of the interferogram best. rection, resulting in a new northwest-southeast trend in the interferograms after corrections. In addition, the WRF product has a larger amplitude for the trend than the ERA product, as can be seen from Figure 5d,e. The regions along the coastline are better corrected with WRF product than with ERA product, while other regions such as the northeast part have more extensive residual phases, illustrating that the larger amplitude in WRF product does not fully agree with the authentic atmosphere. The GPS ZTD can adjust the unrealistic trend in the WRF product, leading to a much better correction by WRFDA product, as shown in Figure 5i, which demonstrates the beneficial contributions from GPS for trend correction. Table 3 summarizes the statistical information, including the amplitude of the atmospheric correction products and the amplitude and stdev of the interferogram after the correction, where smaller amplitude and stdev mean better performance. The pre-corrected interferogram has an amplitude of 59.2 mm and a stdev of 16.25 mm. As shown from Table  3, although the WRF correction product has the largest amplitude of 59.2 mm, which is closest to the amplitude of the pre-corrected interferogram, both the amplitude and stdev after applying the WRF correction product show slight improvement compared to ERA product. Among the three types of products, the WRFDA product provides the best corrections to the interferogram, with amplitude and stdev reduced by 31% and 43%, respectively, mainly due to the assimilation of GPS ZTD. Table 3. Statistical information of correcting results from ERA, WRF, and WRFDA products based on ERAI, including the amplitude of the products and amplitude and stdev of the interferogram after correction. The percentage of reduction compared to the pre-corrected interferogram is given in the brackets.  Figure 6 displays the quantile-quantile (Q-Q) plot of the interferogram with respect to products derived from three methods based on ERAI. If interferogram and correction products have the same distribution, the Q-Q plot should be a straight line. It can be found in Figure 6 that the WRFDA product has the most similar distribution with the interferogram, meaning that WRFDA method can reproduce tropospheric features of the interferogram best.  Figure 7 is the Q-Q plots of the Gaussian distribution w.r.t. the interferogram after atmospheric corrections with three products based on ERAI. If the interferogram after corrections has a Gaussian distribution, it is more likely to contain only random noise, and the plot will be a straight line. From Figure 7, we can also find that the Q-Q plot for  Figure 7 is the Q-Q plots of the Gaussian distribution w.r.t. the interferogram after atmospheric corrections with three products based on ERAI. If the interferogram after corrections has a Gaussian distribution, it is more likely to contain only random noise, and the plot will be a straight line. From Figure 7, we can also find that the Q-Q plot for WRFDA method is closest to a straight line, indicating that most systematic errors such as tropospheric error have been mitigated. The statistic results suggest that the WRFDA product has the most satisfactory performance while the WRF product shows no visible improvement compared to the ERA product.

Products
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 17 tropospheric error have been mitigated. The statistic results suggest that the WRFDA product has the most satisfactory performance while the WRF product shows no visible improvement compared to the ERA product.  Figure 8 presents the atmospheric correction products and correcting results based on ERA5. As for the stratified component, ERA, WRF, and WRFDA products all well capture them and make suitable corrections, for example, regions as indicated in the dashed circles in Figure 8. As for the turbulence component, the ERA product does not reproduce reliable information. It is probably because ERA5 is still sparse compared with the scale of the turbulence though it has an improved spatial resolution of 0.25°. Figure 8b presents the difference between ERA and WRF products, where the distribution shows some smallscale differences. However, these small-scale differences from the WRF product contribute little to the turbulence corrections, as seen from the post-corrected interferogram in Figure 8h. The assimilation of GPS mainly adjusts the trend component instead of the turbulence component due to the density of GPS stations, as shown in Figure 8c for the difference between WRF and WRFDA products. As a result, most turbulences remain in the corrected interferogram in Figure 8i.  Figure 8 presents the atmospheric correction products and correcting results based on ERA5. As for the stratified component, ERA, WRF, and WRFDA products all well capture them and make suitable corrections, for example, regions as indicated in the dashed circles in Figure 8. As for the turbulence component, the ERA product does not reproduce reliable information. It is probably because ERA5 is still sparse compared with the scale of the turbulence though it has an improved spatial resolution of 0.25 • . Figure 8b presents the difference between ERA and WRF products, where the distribution shows some small-scale differences. However, these small-scale differences from the WRF product contribute little to the turbulence corrections, as seen from the post-corrected interferogram in Figure 8h. The assimilation of GPS mainly adjusts the trend component instead of the turbulence component due to the density of GPS stations, as shown in Figure 8c for the difference between WRF and WRFDA products. As a result, most turbulences remain in the corrected interferogram in Figure 8i.

Comparisons at Interferogram Pixels for ERA5
As for the trend components, unlike the ERAI product in the last section, ERA5 reproduces a more consistent trend in the south-north direction with the interferogram. Although the ERA5 trend has a smaller amplitude than that in the interferogram, it still leads to a much better-corrected result than ERAI by comparing Figures 5g and 8g. Surprisingly, the WRF product presents an unrealistic east-west trend, resulting in a worse corrected performance. However, the assimilation of GPS ZTD significantly corrects the trend errors in the WRF product. The WRFDA method corrected interferogram in Figure 8i shows the best results among the three methods. As for the trend components, unlike the ERAI product in the last section, ERA5 reproduces a more consistent trend in the south-north direction with the interferogram. Although the ERA5 trend has a smaller amplitude than that in the interferogram, it still leads to a much better-corrected result than ERAI by comparing Figures 5g and 8g. Surprisingly, the WRF product presents an unrealistic east-west trend, resulting in a worse corrected performance. However, the assimilation of GPS ZTD significantly corrects the trend errors in the WRF product. The WRFDA method corrected interferogram in Figure  8i shows the best results among the three methods. Table 4 shows the statistical information of the corrections based on ERA5. Similar to Section 3.2, WRF and WRFDA products generate larger amplitudes than the ERA product. However, the WRF product underperforms the ERA product, with amplitude reduced by 9% compared with 11% from the ERA product. The WRFDA product holds the best correcting effect with a reduction of amplitude up to 32%. Situations are similar from the perspective of stdev, with decreases of 28%, 22%, and 48% for ERA, WRF, and WRFDA products, respectively. In addition, compared with the statistical results in Section 3.2 for ERAI, products based on ERA5 make more considerable reductions of amplitude and stdev than those based on ERAI except for using the WRF method. Table 4. Statistical information of correcting results from ERA, WRF, and WRFDA products based on ERA5, including the amplitude of the products and amplitude and stdev of the interferogram after correction. The percentage of reduction compared to the pre-corrected interferogram is given in the brackets.   Table 4 shows the statistical information of the corrections based on ERA5. Similar to Section 3.2, WRF and WRFDA products generate larger amplitudes than the ERA product. However, the WRF product underperforms the ERA product, with amplitude reduced by 9% compared with 11% from the ERA product. The WRFDA product holds the best correcting effect with a reduction of amplitude up to 32%. Situations are similar from the perspective of stdev, with decreases of 28%, 22%, and 48% for ERA, WRF, and WRFDA products, respectively. In addition, compared with the statistical results in Section 3.2 for ERAI, products based on ERA5 make more considerable reductions of amplitude and stdev than those based on ERAI except for using the WRF method. Table 4. Statistical information of correcting results from ERA, WRF, and WRFDA products based on ERA5, including the amplitude of the products and amplitude and stdev of the interferogram after correction. The percentage of reduction compared to the pre-corrected interferogram is given in the brackets. Q-Q plots of the interferogram w.r.t products derived from three methods based on ERA5 are presented in Figure 9. Compared with the ERA product, the WRFDA product shares a more similar distribution with the interferogram, while the WRF product presents larger biases. Q-Q plots of the Gaussian distribution w.r.t the interferogram after correction in Figure 10 show that the ERA product is off to the straight line at both ends but still closer than the WRF product. On the other hand, the distribution of residuals after correction with WRFDA product is closest to a Gaussian distribution. In addition, compared with Q-Q plots in Section 3.2 for ERAI, products based on ERA5 has better consistency with the interferogram. Q-Q plots of the interferogram w.r.t products derived from three methods based on ERA5 are presented in Figure 9. Compared with the ERA product, the WRFDA product shares a more similar distribution with the interferogram, while the WRF product presents larger biases. Q-Q plots of the Gaussian distribution w.r.t the interferogram after correction in Figure 10 show that the ERA product is off to the straight line at both ends but still closer than the WRF product. On the other hand, the distribution of residuals after correction with WRFDA product is closest to a Gaussian distribution. In addition, compared with Q-Q plots in Section 3.2 for ERAI, products based on ERA5 has better consistency with the interferogram.

Discussions
Results in Section 3 have demonstrated a significant improvement of the merged NWMs products, indicating the benefits from the assimilation of GPS ZTDs. Comparison between WRF and WRFDA products revealed that the assimilation of GPS ZTD mainly adjusted the linear trend probably due to the relatively low density of GPS stations (with mean station spacing of 35 km) in this case. If denser GPS stations are available, the WRFDA products may potentially capture more small-scale atmospheric variation signals. We also found the limited contribution of the WRF process in this study. The possible reason may be that WRF cannot predict the atmospheric conditions on this day, which can also be inferred from Section 3.1 since WRF estimates have poor correlation with GPS than ERA estimates. Mateus et al. [28] also found that WRF might fail to model the atmospheric delay under some atmospheric conditions accurately. In addition, we can find that the stdev of the interferogram decreased further over 10% and 5% by the ERA and WRFDA products based on ERA5 compared with those based on ERA-Interim. Thus, ERA5 outperforms its predecessor ERAI although it is not that significant. This comparison can give an informative reference to choosing the more proper reanalysis to make corrections depending on requirements.  Q-Q plots of the interferogram w.r.t products derived from three methods based on ERA5 are presented in Figure 9. Compared with the ERA product, the WRFDA product shares a more similar distribution with the interferogram, while the WRF product presents larger biases. Q-Q plots of the Gaussian distribution w.r.t the interferogram after correction in Figure 10 show that the ERA product is off to the straight line at both ends but still closer than the WRF product. On the other hand, the distribution of residuals after correction with WRFDA product is closest to a Gaussian distribution. In addition, compared with Q-Q plots in Section 3.2 for ERAI, products based on ERA5 has better consistency with the interferogram.

Discussions
Results in Section 3 have demonstrated a significant improvement of the merged NWMs products, indicating the benefits from the assimilation of GPS ZTDs. Comparison between WRF and WRFDA products revealed that the assimilation of GPS ZTD mainly adjusted the linear trend probably due to the relatively low density of GPS stations (with mean station spacing of 35 km) in this case. If denser GPS stations are available, the WRFDA products may potentially capture more small-scale atmospheric variation signals. We also found the limited contribution of the WRF process in this study. The possible reason may be that WRF cannot predict the atmospheric conditions on this day, which can also be inferred from Section 3.1 since WRF estimates have poor correlation with GPS than ERA estimates. Mateus et al. [28] also found that WRF might fail to model the atmospheric delay under some atmospheric conditions accurately. In addition, we can find that the stdev of the interferogram decreased further over 10% and 5% by the ERA and WRFDA products based on ERA5 compared with those based on ERA-Interim. Thus, ERA5 outperforms its predecessor ERAI although it is not that significant. This comparison can give an informative reference to choosing the more proper reanalysis to make corrections depending on requirements.

Discussions
Results in Section 3 have demonstrated a significant improvement of the merged NWMs products, indicating the benefits from the assimilation of GPS ZTDs. Comparison between WRF and WRFDA products revealed that the assimilation of GPS ZTD mainly adjusted the linear trend probably due to the relatively low density of GPS stations (with mean station spacing of 35 km) in this case. If denser GPS stations are available, the WRFDA products may potentially capture more small-scale atmospheric variation signals. We also found the limited contribution of the WRF process in this study. The possible reason may be that WRF cannot predict the atmospheric conditions on this day, which can also be inferred from Section 3.1 since WRF estimates have poor correlation with GPS than ERA estimates. Mateus et al. [28] also found that WRF might fail to model the atmospheric delay under some atmospheric conditions accurately. In addition, we can find that the stdev of the interferogram decreased further over 10% and 5% by the ERA and WRFDA products based on ERA5 compared with those based on ERA-Interim. Thus, ERA5 outperforms its predecessor ERAI although it is not that significant. This comparison can give an informative reference to choosing the more proper reanalysis to make corrections depending on requirements.
This study is significant because it reveals that introducing extra high-quality atmospheric information can be equally, if not more, important as improving the temporalspatial resolution when generating atmospheric products. Based on this concept, more promising methods can be explored in the future. For example, assimilating GPS ZTDs several hours before the SAR acquisition time may improve the initial background of the atmosphere and probably obtain better correction products; utilizing domain-dependent background error covariance such as CV5 when running WRFDA may benefit the data assimilation procedure. However, these methods and assumptions need further investigations in the future.

Conclusions
High-accuracy InSAR applications significantly suffer from errors induced by the atmospheric delay, especially the tropospheric delay. This paper compared InSAR atmospheric corrections derived from three methods, i.e., reanalysis (ERAI and ERA5), reanalysis-based WRF, and WRFDA by assimilating GPS ZTD. A case study over a region in Pearl River Delta, China, was carried out to evaluate the performance of different products on correcting the interferogram atmospheric effect. Results revealed that all three methods are capable of reproducing the stratified components well, which is consistent with conclusions from previous studies, e.g., Doin et al. [24], Jolivet et al. [33], Kinoshita et al. [23], Bekaert et al. [36]. However, as for the turbulence, none of the products can capture well, even after improving the spatial resolution up to 3 km in WRF products, and even after assimilating the GPS measurements in WRFDA, which should be mainly due to the limited number of stations. Regarding the trend component, similar to the turbulence, the WRF product shows no superiority to the ERA product, indicating that the higher spatial and temporal resolution in WRF does not help much. It is probably because the WRF method does not introduce any additional observations but only improves the resolution by interpolation. However, after the assimilation of GPS ZTDs, the WRFDA product shows considerable improvements compared with WRF and ERA products for using both reanalyses. After the correction, the stdev of the phase differences in the interferogram decreases up to 48%, and the residuals after the correction show Gaussian distribution when using the ERA5-derived WRFDA product. In addition, if assimilating denser GPS observations, the WRFDA products may be promising to correct small-scale atmospheric delay which needs further investigations in the future.