Utilization of WRF 3D Meteorological Data to Calculate Slant Total Delay for InSAR Atmospheric Correction

The atmosphere introduces excess delays into the synthetic aperture radar (SAR) signal trajectory, especially in the troposphere. InSAR atmospheric correction methods include the use of SAR data and external water vapor products. The latter is more effective. However, since the removal of atmospheric effects should use atmospheric delay products in the direction of the line of sight (LOS), it is necessary to convert the zenith total delay to slant delay in the LOS direction. Conventionally, the zenith delay is divided by the cosine of the average incident angles to obtain slant phase delays. But this method could cause large errors because it ignores the atmospheric horizontal gradient change and the small-scale vertical structure. These problems can be solved by using three-dimensional atmospheric data simulated by numerical models, especially in the case of intense weather changes or complex terrain. However, few scholars paid attention to the application into InSAR atmospheric correction, because of the computation complexity and low efficiency. As the requirement for higher accuracy and the introduction of large errors caused by increasing incidence angles, it is significantly imperative to make the utmost of this method. Weather Research Forecast (WRF) model can provide the precipitate water vapor (PWV) and refraction index at different levels in the three dimensions, and then the slant total delay can be obtained for removing the atmospheric effect on the InSAR process. The results demonstrate that using 3D data can obtain more accurate slant total delay and improve the accuracy of surface deformation from InSAR technology.


Introduction
Since synthetic aperture radar (SAR) was developed, its applications have covered many aspects such as earthquake (e.g., [1][2][3][4]), volcano (e.g., [5][6][7]), mine collapse (e.g., [8,9]), resource exploration (e.g., [10]), and urban ground subsidence measurement (e.g., [11][12][13][14]). All of these applications refer to surface deformation, which is difficult to be measured especially for subtle changes. In terms of the factors affecting the surface deformation accuracy, the atmospheric effect on the SAR signal is considerable, in particular, the water vapor [15]. Zebker et al. demonstrated that 20% spatiotemporal changes of relative humidity in the atmosphere could cause the 10-14 cm surface deformation [16]. And the effect of the uncertainty of 1.0 mm atmospheric precipitable water vapor (PWV) on the standard deviation of Interferometric SAR (InSAR) surface deformation in the direction of line-of-sight (LOS) is about 15 mm [17]. Besides this wet delay caused by the precipitable water vapor, the hydrostatic delay related to the dry air and hydrostatic pressure accounts for about 15% of the total delay, and so it is not negligible [18]. Therefore, eliminating the atmospheric effect on the InSAR results is an important and challenging issue for high-precision InSAR measurement and application. Therefore, spatiotemporal variations of the atmosphere especially the troposphere have a significant impact on InSAR accurate geodetic applications.
To quantify the effects and achieve InSAR atmospheric correction, many researchers introduced external water vapor datasets from MEdium Resolution Imaging Spectrometer (MERIS) ( [19][20][21]), MODerate resolution Imaging Spectroradiometer (MODIS) ( [19,21] Global Positioning System (GPS) stations ( [22][23][24]), and numerical weather models (NWMs) (e.g., [6,18,25,26]). Among them, NWMs are more feasible and potential, including the ERA-Interim from the European Centre for Medium-Range Weather Forecasts (ECMWF) [25,27], the mesoscale analysis model from the Japanese Meteorologic Agency (MANAL) [28], Fifth-Generation PSU/NCAR Mesoscale Model (MM5) ( [18,29,30]), and Weather Research& Forecast (WRF) model [6,26,[31][32][33]. However, these data are not along the SAR signal trajectory (except the advanced synthetic aperture radar (ASAR) and MERIS), and there exists a degree difference between the direction of these data and that of the SAR signal. Therefore, a corrective scheme is applied to transform vertically integrated delay into the slant direction (i.e., the SAR signal trajectory between the satellite and the ground) by dividing the cosine of the mean incidence angle. This simplified scheme could be reliable and applicable when the incidence angles of the SAR satellites are small enough (e.g., less than 32°according to Li [34]) or close to elevation angels [28].
However, atmospheric conditions in zenith and slant direction can vary significantly even over small scales, and if atmospheric noise is comparable to the geophysical signals in the spatial or temporal domains, the simplification method would not work well and can easily evoke errors. Therefore, some mapping functions, including empirical slant delay models and other improved models, are proposed to process observations from radio space geodetic techniques. Boehm et al. proposed the Global Mapping Function (GMF) model based on the ECMWF numerical weather model and demonstrated that their regional height biases and annual errors were significantly reduced with GMF instead of the Niell Mapping Function (NMF). GMF's coefficients were obtained from an expansion of the Vienna Mapping Function (VMF1) parameters into spherical harmonics on a global grid. Similarly, Lagler et al. showed that their Global Pressure and Temperature 2 (GPT2) could yield a 40% reduction of annual and semi-annual amplitude differences in station heights, overcoming the weakness of the GMF and VMF1, limited spatial, and temporal variability [35]. And Yuan et al. proposed the forecast VMF1 (VMF1-FC) model based on the ECMWF data and showed that VMF1-FC performs better than empirical models such as GPT2 [36]. However, all of these mapping functions are widely utilized in real-time GPS, Global Navigation Satellite System (GNSS), and Very-Long-Baseline Interferometry (VLBI) applications, not in the InSAR measurement due to their short wavelength. But Hobiger et al. applied their mapping function into the InSAR measurement. According to Hobiger [28,37], changes of wet refractivity or temperature in different height and horizontal direction introduce millimeter order biases when the simplified method is applied instead of their exact raytracing method. They concluded that their ray-tracing approach differed from the simple mapping strategy by up to 10 mm and had a crucial impact on improving the interpretation of ground motion obtained from InSAR.
However, in recent years, no scholars utilized this threedimensional (3D) ray-tracing approach or other mapping functions to model the atmospheric total delay and achieve the prevailing InSAR atmospheric correction. The reasons could be as follows. Firstly, it is the complex and timeconsuming computation that prevents this method from general application in InSAR technology. Secondly, considering the interdisciplinary research between meteorology and satellite science, it remains a challenge to obtain the 3D atmospheric parameters and calculate the wet delay and hydrostatic delay. Thirdly, Hobiger et al. tried to study the effects from a single profile to find whether ray-tracing for the InSAR correction is necessary. A single profile could evoke errors because of jump-like changes of wet refractivity among voxels of a SAR image [28]. Last but not least, the incidence angles of most satellites were less than 40°, not true for Sentinel, ALOS-2, and so on. The details of the prevailing incidence angles of SAR satellites are shown in the Appendix, Table 3. Therefore, due to the major impact on improving the accuracy of the reconstructed atmospheric delay, it would be imperative and significant to make further research about obtaining more precise 3D atmospheric data and applying them into general InSAR correction. And then, a faster and powerful tool to calculate the atmospheric delay for InSAR would be a better option for SAR scholars. Zhu et al. also noticed this problem and proposed a method that the troposphere is divided into three vertical layers and determined the weights of all pixels on each layer to obtain the Zenith Path Delay Difference Map (ZPDDM). They perceived that their method considered the atmospheric turbulence and heterogeneity and could improve the InSAR results compared with the conventional approach in the Beijing area. However, MERIS products are easily contaminated by the cloud. Besides, there are only three vertical layers, and actually, they projected the ZPDDM into the slant total delay (SLD) along the LOS direction by dividing by the cosine of the radar incidence angle. Despite these limitations, the great significance and necessity have attracted the attention of many scholars.
As numerical weather models are continuously improving their spatial resolution and updating the topographical and land cover datasets, it becomes more and more attainable to replace the simplification scheme with 3D atmospheric data to obtain a more accurate total delay along the SAR signal trajectory. In this study, we present that WRF simulates 3D atmospheric parameters including pressure, temperature, and water vapor mixing ratio and calculates the atmospheric delays at the acquisition time of the InSAR images along the SAR signal trajectory, which are then applied to eliminating atmospheric effects on the D-InSAR process.
This problem deserves specified and further research and description. What is more, it remains a problem whether and when one has to use the 3D meteorological tropospheric delay to avoid evoking errors considering the specific threshold or range of incidence angles. To solve this problem, different incidence angles are modeled according to the current satellites (see the Appendix), ranging from 10 to 60°, and calculated the errors caused by the difference between the local incidence angles and mean incidence angles.
The paper is constructed as follows: the second part will introduce the atmospheric anisotropy from its horizontal and vertical direction. And then, the principles of calculating the 3D total delay along the SAR signal path are described in the third part, followed by the results and discussion in the fourth part. Last, it is the conclusions.

Atmospheric Anisotropy Analysis
In the process of InSAR atmospheric correction, concerning the simplification method, the propagation trajectory is regarded as a straight line. However, due to the anisotropy of the atmosphere and refraction effect, the meteorological parameters change in the horizontal and vertical direction, especially under the situation of extreme weather. The refraction calculation should be performed on each pixel of the SAR images theoretically, which is relatively cumbersome and difficult. Hence, we propose to use short lines (see orange lines in Fig. 1) to calculate the delay, as Snell's law can be extended to the atmosphere in which the refractive index changes continuously with the height. Assuming that the atmosphere is spherically layered, the refractive index changes only with height, and the atmosphere is divided into layers of concentric spherical layers. The refractive index in each layer is constant. Propagating along a straight line in each layer and refraction at each interface (see Fig. 1), the line segment trajectory equation in this spherically layered atmosphere can be derived: n × cos θ = constent. Therefore, it is reasonable and attainable to utilize this polyline to achieve the calculation of the slant delay instead of voxels integration. Figure 2a shows the atmospheric pressure has a similar tendency to the refraction index. What is more, the difference is noticeable in horizontal direction especially near the surface with big gradients about a maximum value of 5 kPa within 5 km. Figure 2b illuminates more considerable differences in the parallel (see purple lines). The values in the green line are different from those in the red line, with a maximum value of 10 mm on the 12th level. The lateral gradients are too big to be neglected. Using the green line instead of the red line would inevitably evoke great errors in the calculation of the total delay, particularly in the bottom levels.
As for the vertical variations, as the altitude increases, the atmospheric pressure that affects the static delay of atmospheric fluid changes with height. As shown in Fig. 2a, it can be seen that when the height reaches 23 layers (corresponding to a height of 15 km), the atmospheric pressure intensity changes slowly. Similarly, the atmospheric humidity delay changes (Fig. 2b). It is getting smaller and smaller. At the height of the 16th level (about 10 km), the atmospheric wet delay is no longer change and less than 2 mm. The standard deviation (std) of zenith total delay (ZTD) in each layer is shown in Fig. 3. As we can see, if the accuracy of 2 mm is required, setting the number of layers to 20 can meet the requirement. Therefore, the model setting layer does not need too much, and the elevation setting can be reduced to the atmospheric troposphere range accordingly.

The Calculation of Atmospheric Phase Delay
Atmospheric refractivity N is computed by Eq. (1) as follows as demonstrated by Smith and Weintraub [38].
where P h and P w are the partial pressures of the hydrostatic air and water vapor in hPa, respectively, and T is the absolute temperature. k 1 = 77.604 KhPa −1 , k 2 = 70.4 KhPa −1 , k 3 = 373.900 K 2 hPa −1 according to Bevis et al. [39]. Snell's law can be used to calculate the relationship between the refractive index and angles of incidence and refraction when the SAR signal goes through the different density of the atmosphere boundary. On the boundary line between different layers, it is necessary to calculate the refraction effect by Eq. (2). We calculate n k at different levels according to formula (1) and then calculated the apparent zenith angle θ k by Snell's law (see Formulas (2) and (3)).
With respect to atmospheric parameters, the 3D fields of the temperature, atmospheric pressure, water vapor mixing ratio, geopotential, and hydrostatic pressure can be obtained from WRF, and PWV at the acquisition times of SAR images can be calculated by Eq. (4) as follows: where ρ water is the water density, k is the vertical layer, and R d = 287.0583 J/(K · kg). P k is the atmospheric pressure, T vk is the virtual temperature, Q k vapor is the water vapor mixing ratio, and Δz is the total geopotential height in the kth layer, less than 700 m. Then, the zenith wet delay (L ZWD ) can be calculated by Eq. (5). And the Π −1 is calculated by Eq. (6).
where R w = 461.495 J/(K · kg), and T s is the atmospheric temperature on the surface in K. The hydrostatic delay at each level (L k ZHD ) is computed by formula (7), which was improved by Elgered's research [40].
P sk , λ k , and H k are the total pressure in millibars, the latitude, and the height in kilometers, respectively, at the profile surface at each level.
And then, the vertical ZTD could be by adding the L ZWD and L ZHD at each level. Therefore, we can obtain the SLD by formula (8). Finally, integrate the SLD along the SAR signal trajectory and obtain more accurate values. It is noted that the SLD (k, i, j) could be SLD (k − 1, i + 1, j) if the horizontal distance is more than the spatial resolution 1 km when integrating along the SAR signal path (Fig. 4).
The WRF resolution is 1 km while the SAR spatial resolution is meter level. To address this difference, we use kriging

Interferometric Process and Atmospheric Correction
Two SAR images can be interfered and generated into an interferogram pair, and the interferogram includes five-phase components (see Eq. (9)) [26].
where Δ means the difference between the master image and the slave image. Δφ InSAR is the differential interferometric phase between two acquisition times, Δφ orbit is the phase resulting from the curved geometry of the Earth that is eliminated by using the precise orbit data, Δφ topo is the topographic phase which can be removed by subtracting a topographic phase simulated from a digital elevation model (DEM), Δφ noise is the phase noise resulting from the decorrelation of the InSAR signal because of land coverage which can be eliminated by filtering, Δφ def is the phase caused by surface deformation along the LOS direction, and Δφ atmo is the phase contribution of the atmosphere. After conventional D-InSAR processing using the software GMTSAR, only the surface deformation term Δφ def and the atmospheric contribution term Δφ atmo are left on the right side of Eq. (9).
The atmospheric phase delay is calculated by formula (10).
where λ is the wavelength of satellites. SLD m and SLD s represent the atmospheric delay at the acquisition time of the master image and slave image, respectively. As for InSAR atmospheric correction, the atmospheric phase delay is then subtracted from the unwrapped interferograms to eliminate the atmospheric effects (see Eq. (11)).
And the deformation result obtained from D-InSAR can be calculated utilizing Formula (12).
4 Results and Discussion

Data and Study Area
The study area is located around the Haiyuan Xian, in the western part of China. It is a topographically complex region dominated by multiple large, steep hills (see Fig. 5 (HY)). Atmospheric flow is strongly affected by this landscape fragmentation. Moreover, it is characterized by large variability of precipitation in summer because of its high and complex topography and its continental monsoon climate. The other study area is Taiwan island (see Fig. 5 (TW)), and it is surrounded by oceans and is considerably affected by the atmospheric water vapor. And more, its topography is complex

Estimation of the SLD Obtained by Exact Profile Transformation
According to Section 2, we computed the atmospheric slant delay at each level profile and implement the integral along the SAR signal trajectory.
For the six interference pairs mentioned above, the atmospheric slant delay is calculated and compared with those of the classic simplified method. Take the interference pairs of 201900504-20190516 in Haiyuan and 20160228-20160323 in Taiwan for example (see Figs. 6 and 7). The simulated incident angle uses a maximum of 60°, taking all of the current satellites' incidence angles into account (see the Appendix). It can be seen from the figures that the maximum difference between the two methods is more than 130 mm (see Fig. 7b-a). Theoretically, the classical simplified method ignores the complex atmospheric turbulence and horizontal gradient changes caused by the terrain changes. The method using three-dimensional data and polyline approximation is closer to the actual SAR signal propagation path and can better reflect the actual atmosphere delay effects on SAR signals.
To further quantify and evaluate the effectiveness of the proposed method, the proposed method is used to simulate the slant delay of different angles incidence, ranging 25-60°, and then was compared with those obtained by the classical simplified method and we use 35°as the value of the average angle of incidence. The atmospheric delay differences between the two methods (the value obtained by the method proposed in this paper minus the value obtained by the classic simplified method) are averaged, std, and the range of the maximum and minimum values. The results are shown in Table 1. It can be seen that, in general, the average values of the difference between the two approaches are positive, that is, the atmospheric slant Secondly, to clearly show the difference, the mean values of the discrepancies between SLDs obtained from the two methods at different local incidence angles are shown in Fig. 7 The corresponding slant delay of the 20160228-20160323 SAR interferogram in Taiwan and the difference between the two methods. Here, the maximum incidence angle is 60°, and the mean incidence angle is 35°F ig. 6 The slant delay of the 20190504-20190516 SAR interferogram in Haiyuan and the difference between the two methods. Here, the maximum incidence angle is 60°, and the mean incidence angle is 35°3  Fig. 8. It is illuminated that the difference between the two becomes larger as the incident angle becomes larger, and when the incidence angle is greater than 35°, a sharp increase phenomenon occurs. This phenomenon shows that when the incident angle is small, the influence of the simple method on the result is relatively small, but after the incident angle is greater than 35°, the error starts to increase sharply and cannot be ignored. Moreover, it is worth noting that the blue curve and black curve represent a larger difference than the others, which corresponds to the differential atmospheric phase delay of 20160228-0323 and 20190904-0916 in Taiwan, respectively. This is because it was a moderate rain on February 28, 2016, and a light rain on the September 4, 2019. For severely changing weather conditions, the results obtained by using 3D data to calculate the atmospheric slant range delay are more accurate and reliable. Thirdly, the seasonal variability in Haiyuan (HY) can be found in Fig. 8 because the other three pairs (HY20190504-0516 (red line), HY20191112-1124 (yellow line), and 20200216-0228 (green line)) are smaller than the pair HY20180801-0813 (orange line). It is related to the dry climate of the Haiyuan, only more rain in July and August. This phenomenon is not obvious in Taiwan. We perceive that Taiwan island is significantly affected by the ocean. Hence, its atmospheric condition changes quickly especially regarding the water vapor. It could be demonstrated that both the atmospheric water vapor and seasonal climate have a noteworthy impact on the SLD discrepancies. And the conclusion is the same as the above that it is necessary to utilize the 3D data and profile transformation to compute more precise SLD values.  Taiwan region for example. The deformation results from differential InSAR based on the WRF-SLD with the top incidence angle 46°are shown in the Figs. 9, 10, and 11, respectively. Because of the very short period (12 days or 24 days) and no geological activity during the intervals, their InSAR surface deformations should be 0 mm in theory. All of the results illuminate that the surface deformation maps with 3D scheme atmospheric correction (3D-AC) are closer to 0 mm than with conventional simplified scheme atmospheric correction (CM-AC).
To evaluate the effectiveness of this approach, we assessed the error of the surface deformation from the D-InSAR results with the two different methods. The D-InSAR error results of the interferograms (Ifgs) are displayed in Table 2, and it can be seen that the mean, standard deviation (std), and root mean square (RMS) errors after atmospheric correction with 3D slant total delay (3D-AC) are much smaller compared with those with CM atmospheric correction (CM-AC). This demonstrates that we simulated more accurate SLD and improve the InSAR surface deformation accuracy by using the level transformation.

Conclusions
3D data were utilized to obtain more accurate atmospheric total delay values along the SAR signal trajectory from the vertical ZTD at each WRF level. Firstly, the vertical ZTD and the refractivity index of each pixel at each level are calculated, and then at different levels, the vertical ZTD values are transformed into slant delay based on different refractivity. Ultimately, all of the slant delays at each level are integrated to obtain the atmospheric total delay. This approach makes the most of the 3D atmospheric refractivity index and solves the horizontal gradients between the two adjacent pixels of the WRF model. Meanwhile, it could reduce the errors caused by the atmospheric refraction in the horizontal and vertical directions. Compared with the conventional transformation scheme by dividing the cosine of the mean incidence angle of the satellite, considering the detailed atmospheric Moreover, we simulated slant delay at different incidence angles and found that it is necessary to use this proposed method instead of the traditional method to avoid errors and misinterpretation of the D-InSAR surface deformation when the incidence angle is larger than 35°. In terms of the evoke errors, the maximum error is about 40 mm. Therefore, it is more suitable to utilize 3D data to simulate the atmospheric slant delay and mitigate the atmospheric effects on the D-InSAR measurement.
As the development of numerical models is gradually improved, it is more imperative to use 3D atmospheric data to calculate SLD instead of the simplification scheme for the InSAR application. Integrating the delay along the line-ofsight between the ground point and the satellite is attainable and potential considerably. Applying numerical models into InSAR technology provides researchers with a reliable source of atmospheric data and an effective method in the field of geodetic technology.

0~60°4
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.