A New Method for InSAR Stratiﬁed Tropospheric Delay Correction Facilitating Reﬁnement of Coseismic Displacement Fields of Small-to-Moderate Earthquakes

: Focusing on stratiﬁed tropospheric delay correction in the small-amplitude coseismic displacement ﬁeld of small-to-moderate earthquakes (<Mw 6.5), we develop a Simple-Stratiﬁcation-Correction (SSC) approach based on the empirical phase-elevation relationship and spatial properties of the troposphere, via an equal-size window segmentation. We validate our SSC method using 23 real earthquakes that occurred from January 2016 to May 2021 with a moment magnitude (Mw) ranging from 4.5 to 6.5. We conclude that SSC performs well according to the amount of reduction in semi-variance and the root-mean-square value. This method primarily focuses on stratiﬁcation delay correction; thus, it is especially useful in regions with complex terrain, while it can mitigate partial large-scale turbulence signals. We investigate three parameters that are empirically setup in the correction working ﬂow and inspect their optimal settings, when implementing SSC for quick response after earthquake. Our method is ready to be integrated into an operational InSAR processing chain to produce a reliable atmospheric phase screen map, which can also serve as an auxiliary product to quickly and timely quantify stratiﬁcation delays in coseismic interferograms. Through improved accuracy of the coseismic displacement ﬁeld, the focal mechanism could be better constrained to facilitate the building and expansion of the geodesy-based earthquake catalogue.


Introduction
Interferometric Synthetic Aperture Radar (InSAR) observations are sensitive to the magnitude and location of earthquakes and is widely applied to recognize earthquakes with small moment magnitude (Mw < 6.5) but still having observable surface deformation [1]. Many previous studies have found that InSAR-derived coseismic displacement can be used to determine the location, depth, and coseismic slip distribution, providing more accurate results than those constrained by traditional seismic data [2][3][4][5][6][7]. It is especially useful for monitoring earthquakes that occurred in remote areas, where the density of the in situ observation network is limited. Constrained by InSAR-derived coseismic displacements, various algorithms and community tools have been developed to retrieve the kinematic and dynamic source parameters of earthquakes [8][9][10][11][12], greatly extending our knowledge of earthquake mechanisms. Additionally, the expansion of SAR satellite missions and the development of robust processing techniques in recent years led to increasing interest spatial spectrum decomposition, and then performed the phase-elevation correction in each segment. For example, windows-based algorithms may split images into multiple windows with constant size [41,43], or adjust the window size based on local elevation gradients [44]. Lin et al. [38] developed a multi-scale approach via band-pass filtering to estimate and remove the tropospheric delay components. Bekaert et al. [42] developed a power-lawbased phase-elevation model and applied a band filter to reduce the displacement signal before modeling tropospheric phases. Murray et al. [40] segmented each interferogram using a K-means clustering algorithm and constructed a tropospheric delay screen from each cluster.
To improve the accuracy of the observed coseismic displacements and better serve various geophysical applications, our study aims to develop a new tropospheric artifact mitigation method, drawing on the empirical phase-elevation relationship and spatial variant atmosphere. In the following discussion, a brief introduction to 23 small-to-moderate earthquakes in our test, InSAR datasets and the details of our proposed method will be described. The quality of our correction algorithm will be subsequently validated by applying it to multiple earthquakes and comparing the results to GACOS products. Finally, we discuss various parameter settings in our method and related them to atmospheric physical properties. The discussion will help users to better apply our method and to understand the composition of the modeled tropospheric phases.

Atmospheric Stratification Delay Estimation
The changes of tropospheric reflectivity between two SAR image acquisitions causes a phase shift that would interfere with surface deformation signals. According to the physical origin of tropospheric signal, interferometric phases include contributions from the turbulent mixing component and the vertical stratification: 1.
The turbulent delay component typically results from turbulent processes in the troposphere and leads to three-dimensional spatial heterogeneities in refractivity [20]. The stochastic model of tropospheric spatial variability shows a power-law dependence on frequency [36]. The contribution from this type of component could be reduced by averaging independent interferograms [19].

2.
Vertical stratified delay is caused by the different vertical refractivity profiles of two acquisitions, assuming no heterogeneity within the horizontal layers [20]. In contrast to the turbulent component, the vertical stratified delay is correlated with topography and affects regions with mountainous terrain conditions [20,21]. Additionally, as revealed by the ERA5 weather model dataset, the stratified delay features a seasonal fluctuation [35]. Many methods rely on a suite of assumptions about the spatialtemporal characteristics of InSAR signals (e.g., the linear progression of deformation over time and the zero-mean Gaussian nature of atmospheric phases) that are often not valid in cases of the stratification delay. This is of great significance for many earthquake studies using InSAR. The coupling in the temporal domain further complicates the correction of stratification delays. Thus, the estimated surface deformation may be ambiguous, and inappropriate correction may lead to an erroneous interpretation of the earthquake deformation.
Mitigation of the stratification delay can be performed based on the empirical relationship between the interferometric phase and elevation [31,[38][39][40]45,46]. Cavalié et al. [39] approximated the relationship as a linear function between height and phase: where K denotes the slope relates to topography H and phase ϕ stra , and C denotes the constant parameter. In order to better approximate the phase spatial variation at a higher part of the atmosphere, Bekaert et al. [41] described phase-elevation dependence using a power-law decaying model: where A is a scale factor that relates to the elevation, α is a constant factor describing the power-law decaying, h 0 is the reference height where the tropospheric phase converges to zero (h 0 suggested to be~7-13 km [41]), and thus h 0 − h denotes a height difference relative to the reference height at each pixel.

Correction Working Flow
Focusing on the applications to earthquake studies, especially with respect to small-tomoderate earthquakes, we introduced an empirical phase-elevation-model-based tropospheric delay correction method, Simple-Stratification-Correction (termed SSC hereinafter). Two challenging factors should be considered are (1) avoiding partial or complete removal of real deformation signals, because the deformation may be correlated to elevation [38,40], and (2) preparing for significant spatial variation of the tropospheric properties at a severalkilometer scale [41,43,44]. Here, we masked the main deforming zone for each earthquake and segmented the images into small windows to deal with these two problems. The detailed steps are described below:

1.
We masked the coseismic displacement area in individual interferograms. The mask was generated based on the reported epicenter. We aimed to avoid participation of pixels dominated by coseismic displacements in estimating phase-elevation model parameters.

2.
We then segmented each interferogram into multiple small patches. We cropped the interferogram into M × 2 N by M × 2 N dimensions and the coseismic zone had to be included in this step. We aimed to roll-change the window numbers by a factor of 2 in a more automatic manner. The choice of the number of segmented windows was related to the scale of the estimated stratifications and the size of the earthquakes. We found that a patch number of 8 or 16 in both range and azimuth was appropriate for most of the earthquakes tested here. More scenarios on properly splitting windows will be discussed in Section 4.

3.
Each window contained a cluster of pixels, which were subsequently utilized to estimate the phase-elevation model parameters via the empirical linear model (Figure 1). For a window partially impacted by coseismic deformation, the percentage of deforming pixels was computed in conjunction a with previously determined mask of the coseismic zone. The empirical threshold was 60%, meaning that if >60% pixels were not masked, they were used to estimate the phase-elevation model parameters in this window. However, if the percentage was smaller than 60%, this window was recognized as a masked one and the corresponding parameters estimation were skipped here. 4.
The next step was to fill the masked windows that failed in the direct estimation in the previous step. We calculated the semi-variogram structure from the input phase over the non-deforming zone. This step aided the kriging solution in predicting the model parameters at those windows dominated by coseismic deformation zones.

5.
Based on the obtained sparse C (constant parameter) and K (scale parameter) grid, we applied the kriging interpolation to calculate the values at all pixels. For each pixel, we obtained a pair of K and C. Given the known height, we applied these resolved values back to the empirical phase-elevation model to calculate φ stra . Note that the estimated parameters in the last step were assigned to be at the center of each window. Therefore, after the spatial interpolation, it left an empty area that was not computable. This appears in Figure 1b as the region outside the dashed box. The width of this zone is the half of a segmented window. 6.
Finally, we subtracted the modeled tropospheric delay from each unwrapped interferogram to obtain the corrected coseismic displacement map. To better understand the overall procedure of SSC, Figure 1 takes the M 6.0 earthquake that occurred on 25 August 2018, near Javanrud, Iran, as an example and provides a visual demonstration. The unwrapped interferogram are segmented into 8 × 8 windows. Six phase-elevation model examples are shown in Figure 1d-i. For each plot, the color of the marks corresponds to the box with the same color in Figure 1a. They have different terrain conditions. The one in Figure 1f is partially masked, due to the nearby coseismic zone.
We calculated the standard deviation (SD) of the regional height inside each window and the coefficient of determination (R 2 ) for each linear regression. Both calculations were annotated in the corresponding plots. Of the unmasked windows, 60% had height SDs larger than 200 m, and 50% achieved R 2 values larger than 0.5. This indicated that this coseismic example had certain topography relief and the regression reached an ordinary quality. Figure 1h,i shows two windows as examples over a flat area. According to the distribution of these phase-elevation samples, the stratification relationship was weak here and the main atmospheric impact of this local zone was a likely turbulence component. Thus, the absence of terrain relief was a main factor that reduced the quality in the phaseelevation modeling.
In summary, the key objectives of our SSC solution are: (1) estimating atmospheric phase from the non-deforming pixels, implemented by introducing a proper mask for each earthquake; and (2) combining the phase-elevation dependence model with spatial property of the troposphere.
Note that we have also tested with the power-low model and integrated that test into SSC scheme. However, after the interpolation step (step 4), it failed in producing a physically reasonable atmospheric delay map. We suspect that the parameters of the power-law model are not likely correlated at the spatial scale we used in SSC. Therefore, our subsequent real earthquake test was performed using the linear phase-elevation model.

Earthquake Catalog and Coseismic InSAR Observations
We tested our SSC method with a total of 23 earthquakes that occurred from January 2016 to May 2021 ( Table 1). The moment magnitude ranged from Mw 4.8 to 6.5. We excluded earthquakes that occurred on flat areas as we are focusing on the stratification correction. We also discarded those earthquakes that occurred at the seashore, because there were not enough coherent pixels surrounding the coseismic zones, thereby causing difficulties in determining the phase-elevation dependence. All of these earthquakes were observed by Sentinel-1 SAR images, and three of them had both descending and ascending interferograms. The majority of the coseismic interferometric pairs were downloaded from SARVIEW projects (https://sarviews-hazards.alaska.edu/, accessed on 4 March 2022). SARVIEW is a fully automatic D-InSAR processing system for Sentinel-1SAR acquisitions; it includes interferometric pair selection, phase filtering, unwrapping, geocoding, and other post-processing steps [47]. The focal mechanism parameters of the 22 May 2016, Dingjie earthquake and the 8 August 2017, Jinghe earthquake were obtained from Hou et al. [33] and Gong et al. [48], respectively, while the others were obtained from Zhu et al. [13]. We limited the temporal baseline of each interferogram to the shortest revisit period, generally 6-12 days, to maximize the coherence. It was also important to make sure that there were enough pixels surrounding the coseismic zones to ensure the sufficiency of non-deforming pixels in the atmospheric delay modeling.

Results
Our proposed atmospheric delay correction approach, SSC, was evaluated using 23 earthquakes in nature. Three of them had coseismic observations on both the ascending and descending tracks. Therefore, a total of 26 coseismic interferograms were used in the real data validation. We calculated the statistics before and after the corrections, in order to better quantify the effectiveness of such corrections. In order to compare the correction quality, we also applied correction with the Generic Atmospheric Correction Online Service (GACOS; http://www.gacos.net/, accessed on 4 March 2022). The GACOS product has global availability and provides high spatial resolution and easy-to-implement zenith total delay maps [31].
We took the epicenter shown in Table 1 as a priori information to mask the coseismic deformation zone. We solved the phase-elevation dependence, as described previously. We then reconstructed the atmospheric stratification delay for each coseismic interferogram and removed that factor from the original unwrapped phase. The resulting corrected coseismic displacement maps could be further used in slip distribution inversion in geophysical studies.
As shown in Figures 2 and 3, we took two earthquakes, the Mw 4.8 Mexico earthquake and the Mw 6.5 Iran earthquake, as examples for our demonstration. These two earthquakes had the minimum and maximum magnitudes in all our tested earthquakes, respectively. These two figures demonstrated interferograms before and after the atmospheric delay correction using our SSC method and the GACOS products. For these two examples, both the GACOS and SSC reduced the noises, but at different levels. The residual signals after SSC in the corrected interferograms (Figures 2c and 3c) were likely the localized phase noises (e.g., the small scale turbulence or localized stratification). The plots (d,e) of Figures 2 and 3 are the scale and offset parameter maps after kriging interpolation. The unit of digital elevation applied in the modeling was meter; therefore, Figure 2d shows a relatively small value range. Both Figures 2f and 3f show the displacement differences along a profile across the coseismic zone after the correction. Figures 2g and 3g are the histograms calculated from non-deforming pixels that were identified using the same masks in SSC implementation. They show an apparent change in the shape of displacement distribution, where the SSC largely reduced the tropospheric delays. The correction results for all other earthquakes are demonstrated in the Supplementary Materials (Figures S1-S21).
Additionally, we calculated the Root-Mean-Square (RMS) value before and after the atmospheric correction over the non-deforming zone as a statistic factor to evaluate the overall correction quality. This value was used to indicate an overall level of phase noise in each interferogram. The result is summarized in Figure 4. After implementing our SSC, most interferograms had a reduced RMS by 45%.   We also calculated the semivarigoram to quantify the reduction in the noise magnitude at various spatial scales. The semivariogram depicts the spatial autocorrelation of the sampled measurements and was used to quantify the atmospheric impacts in many previous studies (e.g., [20,40,49]). Note that we excluded pixels at deforming zones in the calculation.
As shown in Figure 5, SSC performs well for all tested coseismic interferograms, with a significant reduction in semi-variance. In this study, the SSC-corrected interferograms showed a reduction of 5 to 20 km at an overall spatial scale. This was likely related to (1) the spatial size of the coseismic zone; (2) the parameters used in the SSC correction, especially the number of windows; and (3) the spatial scale of the regional terrain. Here, we took the Mw 6.5 Sarpol-e Zahab, Iran, earthquake ( Figure 3) as an example. After cropping the interferogram into a 2 N by 2 N dimension, the subset was then split into 8 × 8 windows and each window corresponded to a spatial size of~17 km × 17 km. Note that we did not expect to accurately reconstruct the real phase-elevation parameters below this spatial scale (e.g., smaller than a single window size). Hence, the atmospheric signals related to local terrain or turbulence signal with small spatial scale remained as residuals in the corrected interferograms. Overall, the SSC method provided a satisfactory performance, validated by both semivariance and RMS. Another advantage of the SSC method was that it only used a single interferogram. Compared to the multi-temporal-analysis-based tropospheric correction methods (e.g., [32][33][34]), our SSC could perform as long as a single interferogram was available. In this study, the GACOS products could also reduce partial atmospheric delays. The inconsistent correction quality of GACOS may be caused by the input atmospheric reanalysis products or by the GNSS meteorological products in each single case.

Discussion
In this section, we inspect the optimal parameter settings in SSC and its performance dependence to regional terrain. We also discuss the link of the parameters of the empirical phase-elevation model to the tropospheric contributions. This analysis help us in understanding which part of tropospheric artifacts could be modeled.

Parameter Setting
To implement the SSC, there were three parameters that had to be empirically set up. The first one was the size of image subset covering the coseismic zone, which was cropped from the uncorrected unwrapped interferogram. This was not a critical step, and the SSC could be implemented into an entire interferogram. This step was conducted in order to improve computation efficiency and to exclude large water bodies (e.g., oceans). To generate this subset, we required (1) enough land area surrounding a coseismic zone, ensuring sufficient pixels for phase-elevation parameter estimation, and (2) an image subset having a dimension of M × 2 N in width and length, ensuring that the number of windows in both directions was ×2 i ( i < 16).
The second empirically determined input was the mask of the coseismic zone. The mask, or an equivalent method of reducing the estimation bias caused by coseismic deformation, could be realized via a spatiotemporal filter. For example, Murray et al. [40] built a temporal low-pass filter to identify and mask deforming pixels. Bekaert et al. [42] applied a band filter, which was constructed based on an initial analysis with InSAR and GNSS, to avoid contamination of the tectonic slow slip deformation signal. Given that the rupture length of small-to-moderate earthquakes with Mw ≤ 6.5 may be limited to 10-20 km in length or less (e.g., an Mw 6.5 event may produce a rupture of about 20 km [50]), and concentrating on a known tectonic area from a seismically determined epicenter, we built a rectangle mask for each earthquake. The spatial extent of the mask was generally controlled by the earthquake's magnitude and focal depth. In this study, the mask zone was manually set up mainly based on USGS reported epicenters. The overall length of the masks for all of the events here ranged from 3.3 to 33 km. Here, we also proposed a potential treatment to generate a proper mask in an automatic way, to allow for possible use in future applications. First, one can predict a spatial size of the coseismic displacement zone with a facility's (e.g., via U.S. Geological Survey Earthquake Hazards Program) reported focal mechanism parameters. In other words, the size of the mask could be considered as a known. Thereafter, one could search for an optimal masking location by minimizing the standard deviation of un-masked pixel phases. As suggested by a previous study [51], the average distance between an InSAR-reported epicenter and one from USGS was 3.7 km for Iranian earthquakes, and 2.7 km for Japanese earthquakes. Therefore, the search could be started from USGS-reported epicenters and bounded in a limited zone within a several-kilometers radius.
The third consideration is the number of segmented windows. As shown in Table 1, we preferentially segmented the interferogram into 8 × 8 window subsets for most of the events. For a few earthquakes with smaller magnitudes or deeper depths, we applied finer window segmentation (16 × 16). For earthquakes with similar moment magnitudes (Mw 4.5-6.5), an 8 × 8 window segmentation was sufficiently effective. To guarantee the correction quality, we suggest not splitting an image into very fine equal-size windows. Because the number of pixels of a smaller window would be reduced, the chance of them capturing stratification delays presented in areas with a sufficient signal-to-noise ratio became smaller. This would reduce the robustness of the phase-elevation model's parameter regression. In addition, it would leave more "blind" windows sitting inside the coseismic displacement zone. As mentioned above, for those windows containing insufficient number of pixels to capture the atmospheric delay, calculation could only be done by kriging interpolation instead of by direct phase-elevation regression. However, the derived windows could not be too coarse; otherwise, a single window could be controlled by a multiple trend of phase-elevation dependence.

Relationship with Regional Terrain
We compared the percentage of the RMS reduction value with the regional terrain conditions in order to quantify the terrain's dependence of correction quality. This helped us to better understand its future application scenario. We took the RMS in Figure 4 and calculated the corresponding reduction percentage in each case. Figure 6 visually demonstrates the comparison of the mean and standard deviation (SD) of the regional elevation. It shows that a reduction was >60% for most of the tested cases. For those with an average elevation <2.5 km, the reduction ranged from 30% to 90%. For the other interferograms with a higher average elevation, the correction percentage was >50%. However, we did not see a correlation between the quality of correction and the regional terrain complexity according to the SD of the height.

Regression Models and Parameters
As presented above, the linear approximation [39] between atmospheric delays and local topography can be described by a scale factor and a constant. If we apply the phaseelevation model to an entire interferogram, the constant parameter is only a shift of the interferogram [45]. The estimation of phase-elevation model parameters is implemented for each individual window in the SSC method. Thus, the estimated constant parameter is an overall phase shift of each window. This shift describes the regional long wavelength signals and should be smooth in space. Therefore, after the spatial kriging interpolation, the whole map of the constant parameter approximates the long wavelength signal, which is likely composed of both turbulence and stratification components, while the estimation of the scale parameter (K in Equation (1)) describes the relationship between the phase and the regional elevation. Both the de-correlated pixels and those pixels within the coseismic zones should be pre-masked before the estimation.
After spatial interpolation of the K field, we obtained a pixel-by-pixel map that could be used to model stratification within the interferogram subset. Therefore, according to meaning of the parameter pair in the phase-elevation model and also the choice of spatial interpolation in the SSC procedure, this method estimated not only stratification but also some amount of large-scale troposphere phases.
As suggested by Bekaert et al. [41], the linearization assumption is only valid for the lower part of the atmosphere. Bekaert et al. suggested a power-law model that predicted that the phase delay would converge to zero at larger elevations. To integrate the powerlaw model into our SSC method, we modified the power-law model by adding an offset parameter ε (see Equation (3)).
This is because the offset term can be used to represent the overall shift within each window that is caused by the atmospheric delays with a large spatial scale, as discussed above. Figure 7a-c shows the SSC resolved three parameters (scale parameter A, power parameter α, and offset parameter ε) at each window. Figure 7d-f shows the results after applying kriging interpolation. Figure 8 demonstrates the resulting stratification delay maps based on the window-wise parameters and interpolated parameter maps. Figure 8a indicates that the power-law model works well when each window is only associated with a single set of estimated parameters. However, we could not produce a meaningful phase-delay result (Figure 8b) when using interpolated parameter maps. Even the overall spatial pattern after kriging remained similar to the original window-wise parameters (see Figure 7). The failure here was likely because these variables in the power-law model remain invariant across a certain area; thus, they do not satisfy the spatially smooth assumption in the SSC strategy. Figure 7. Demonstration of power-law-based SSC method: (a-c) show the estimated windowwise scale factor, the power factor, and the offset factor in Equation (3), respectively, and each window is described by a single set of variables; (d-f) are parameter maps after spatial interpolation, corresponding to the scale factor, the power factor and the offset factor, from left to right. The dashed box in Figure 8 denotes the masked coseismic zone, while the solid box denotes windows where the nonlinear inversion failed in finding optimal values. The failure was mainly because the phase-elevation relationship was weak at these windows (e.g., a flat zone). For example, the dark blue window in Figure 8 bounds a region identical to that outline by the box with the same color in Figure 1a. Both the linear model and the power-law model show unsatisfied performance in this area. The region in the red solid box has some terrain relief (with a height SD of 180.0 m), while also failed in modeling power-law model parameters. The difficulty was caused by the mixing of multiple phaseelevation dependence in this window. Overall, the linear-model-based SSC method is preferred and performed well on the tested coseismic cases.

Conclusions
InSAR observations provide rich products for building a geodetic earthquake catalogue. Nevertheless, the accuracy and the capability of earthquake deformation observations are largely dependent on the atmospheric noisy level of the derived displacement field. Focusing on the coseismic displacement field reconstruction for small-to-moderate earthquakes (<Mw 6.5), we developed a Simple-Stratification-Correction (SSC) method based on the phase-elevation relationship and spatial smooth property of the troposphere. SSC employs a windows-based segmentation strategy, given its simplicity and efficiency. By applying the spatial-correlation property of atmospheric delays, SSC allows for the phase-elevation dependence variant in space. The developed SSC method is especially useful for zones with complex terrain, where the stratification delay is most severe. Our method is also helpful in mitigating some large-scale turbulence signals, according to the property of the offset term in the phase-elevation model and its usage in the SSC approach.
In SSC method, a proper mask is necessary to exclude impacts from the coseismic signal in the tropospheric phase modeling, which can reduce the bias from coseismic displacement signals. We also determined the optimal parameter setting to implement SSC for earthquake events with moment magnitude (Mw) 4.5 to 6.5. SSC only needs a single interferometic pair, rather than an interferogram stack, and thus can provide a quicker response after an earthquake. A single interferogram can maintain better coherence with the shortest available temporal baseline. This method can also be integrated in time series analysis that would have a broad implication for other tectonic studies [52,53], e.g., assessing interseismic displacement [54,55]. The key to initiate such applications with SSC is to provide a proper mask or other analysis schemes (e.g., integrating it in a small baseline approach [56]) that can reduce the bias in atmospheric signal modeling.
Overall, by improving the accuracy of the co-seismic displacement field, one could retrieve the seismogenic parameters with better quality that facilitate the building and expansion of the geodetic earthquake catalogue. Thus, it allows us to extend the capacity of InSAR in detecting and monitoring small earthquakes, increasing the number of the geodetic-based earthquake studies, and improving the quality of geodetic earthquake catalogues. Our method is ready to be integrated into an operational InSAR processing chain, e.g., one can produce an additional atmospheric phase screen map via SSC as an auxiliary product to help users understand the amount of stratification delays in each coseismic interferogram.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14061425/s1. This file contains Figure S1-S21: Interferograms before and after tropospheric corrections of the rest 21 earthquake events tested in this study are listed here. Figure Figure S8: Correction result for the Mw 6.14 earthquake in Torbat-e Jam, Iran, on 5 April 2017. Figure S9