Tropospheric Delay Correction Based on a Three-Dimensional Joint Model for InSAR

: Tropospheric delays in spaceborne Interferometric Synthetic Aperture Radar (InSAR) can contaminate the measurement of small amplitude earth surface deformation. In this paper, a novel TXY-correlated method is proposed, where the main tropospheric delay components are jointly modeled in three dimensions, and then the long-scale and topography-correlated tropospheric delay components are corrected simultaneously. Moreover, the strategies of scale ﬁltering and alternative iteration are employed to accurately retrieve all components of the joint model. Both the TXY-correlated method and the conventional phase-based methods are tested with a total of 25 TerraSAR-X/TanDEM-X images collected over the Chaobai River site and the Renhe Town of Beijing Shunyi District, where natural scenes and man-made targets are contained. A higher correction rate of tropospheric delays and a greater reduction in spatio-temporal standard deviations of time series displacement are observed after delay correction by the TXY-correlated method in both non-urban and urban areas, which demonstrate the superior performance of the proposed method.


Introduction
Spaceborne Interferometric Synthetic Aperture Radar (InSAR) is an efficient and powerful tool for ground surface deformation mapping with high resolution, high precision, all-weather and all-time capabilities [1,2]. However, during the propagation of radar signals, atmosphere may result in phase delay or advance [3,4], which becomes an unfavourable factor for the detection of small amplitude deformations. Therefore, removing the atmospheric effects is indispensable for accurate deformation estimation in InSAR.
Atmospheric delays mainly consist of ionospheric delays and tropospheric delays. The ionosphere causes phase advance, especially for long wavelength radar signals, such as P-band or L-band SAR When it comes to C-band or X-band, ionospheric delays can be neglected [5]. The troposphere leads to phase delays due to the variations of pressure, temperature and water vapor content in space and time. The delays may cause large fluctuation of deformation, which covers up the real deformation signals [6]. The tropospheric delays are mainly composed of short-scale, long-scale and topography-correlated components [7,8]. The short-scale delay, also called the turbulent delay, is induced by turbulence in air moisture. Both the long-scale and the topography-correlated components belong to the stratified delay. The long-scale component is mainly related to lateral variation of pressure and temperature and relative humidity in space, while the topography-correlated component is mainly related to vertical variation of pressure, temperature and relative humidity in height [8]. Besides, the stratified delay can also be divided into the hydrostatic delay and the wet delay physically [9]. The hydrostatic delay is by the TXY-correlated method in both non-urban and urban areas. Thus, the proposed method can correct the tropospheric delays more accurately than the conventional one.

The Three-Dimensionally Joint Model
When electromagnetic waves pass through the troposphere, the change in refractive index is responsible for the tropospheric delays. Thus, the line of sight (LOS) single path tropospheric delays δL LOS (z) are expressed physically as the integral of the air refractivity between the ground target elevation z 0 and an reference elevation z re f [19], which can be modeled as [32]: where θ is the incidence angle, g m is the average gravity acceleration between z 0 and z re f , R d = 287.05 J kg −1 K −1 is the dry air specific gas constant, R v = 461.495 J kg −1 K −1 represents the water vapor specific gas constant, P is the total air pressure in Pa, e is the partial pressure of water vapor in Pa, T is the temperature in K, k 1 = 0.776 K Pa −1 , k 2 = 0.716 K Pa −1 , and k 3 = 3.75e3 K 2 Pa −1 are empirical constants [33]. This formulation mainly considers the stratification component of tropospheric delays, and does not account for the turbulence component. Once the temperature, air pressure and water vapor partial pressure are acquired by meteorological observation, such as the global atmospheric models (GAM) [11] or MERIS/MODIS [15,16], the absolute tropospheric delays at the observation time can be calculated. Thus, the interferometric tropospheric delay phase ∆φ t 1 t 2 LOS (z) can be computed by combining two acquisitions at times t 1 and t 2 as where 4π/λ is the phase conversion factor, and λ is the wavelength of electromagnetic wave. However, the application of meteorological data is often limited in time and spatial resolution, and the acquisition time of meteorological data is often inconsistent with the observation time of SAR, which reduces the accuracy of tropospheric delay correction in SAR interferogram. Therefore, the tropospheric delay correction model based on the interferogram itself for InSAR is the focus of this paper. The Stanford Method for Persistent Scatterers (StaMPS) proposed by Hooper [34,35] is used to perform InSAR time series analysis. The key steps of StaMPS are shown in Figure 1 marked with black dotted boxes. Suppose a total of K interferograms are generated with respect to a selected master image, and a total of P permanent scatterer (PS) points, defined as a series of pixels with stable phase over a long interval, are selected. For the k-th interferogram, the unwrapped phase ∆φ k i,unwrap after 3D phase unwrapping, Digital Elevation Model (DEM) error correction and orbit error correction at the i-th PS point can be decomposed into [34] where ∆φ k i,de f o is the phase caused by temporal displacement of the target between the acquisitions, and ∆φ k i,noise is the contribution caused by thermal noise and preprocessing residual errors. ∆φ k i,tropo is the interferometric phase caused by tropospheric delay with the following three components [8] where ∆φ k i,long , ∆φ k i,t−c and ∆φ k i,short are the long-scale, topography-correlated and short-scale components, respectively.
The long-scale tropospheric delay has a stable characteristic and mainly exists in flat areas [8]. The topography-correlated tropospheric delay is correlated with topography and mainly exists in mountainous regions [7]. The short-scale tropospheric delay is generated in the near-ground surface troposphere, and it exists in both flat and mountainous areas. Due to significant variation of tropospheric water vapor content over a short period of time, it is impossible to describe the short-scale component with a deterministic phase model. However, the short-scale component is generally random in time and space, which can be eliminated directly by temporal low-pass filtering or smoothing in time series analysis [26,28]. Thus, in this paper we mainly model the long-scale and topography-correlated tropospheric delays.
Due to spatial correlation of the long-scale tropospheric delay, a traditional two-dimensional XY-correlated model has been proposed by Ferretti et al. [26]. The long-scale tropospheric delay phase ∆φ k i,long of the i-th PS point for the k-th interferogram can be expressed as where ξ i is the azimuth coordinate of SAR image, and η i is its slant range coordinate, A k and C k are the corresponding estimated slope factors which are global constants in an interferogram, and B k is a constant deviation to the full interferogram that can be neglected. However, this model limited in small areas [27], since the spatial correlation of the long-scale tropospheric delay decreases when the studied regions become large. To obtain more accurate estimation results over a relatively large region, the block processing strategy is introduced to calculate slope factors for each block. When values of each block are obtained, local slope factors for each PS point can be acquired by multi-weight interpolation. Therefore, the corresponding model can be modified as where A k i and C k i are slope factors for each PS point. Once the long-scale tropospheric delay phase of the PS points has been estimated, they can be interpolated to obtain a uniform image grid.
The phase variation caused by the topography-correlated tropospheric delay can be empirically described by a linear relationship with topography [29,30]. Thus, the topography-correlated tropospheric delay phase ∆φ k i,t−c of i-th PS point for the k-th interferogram can be represented as where ∆φ k 0 denotes a constant shift applied to the whole interferogram, and K k ∆φ is the transfer function between topography h i and phase ∆φ k i,t−c . The transfer function K k ∆φ is linearly fitted better by a global regression than by a local one [36].
Based on the above models, the three-dimensional joint model, representing the tropospheric delay phase, is proposed as where D k represents an overall bias of the full interferogram. The main improvement of the proposed model is that the major tropospheric delay components are all considered. As the joint model includes both global and local parameters, it is not realistic to apply linear regression directly to solve the parameters. Moreover, the tropospheric delay phase does not seem linearly correlated with topography and the spatial correlation of long-scale tropospheric delay can be reduced when multiple confounding effects exist, such as ground displacement, atmospheric disturbance and noise [26,31]. In other words, the presence of multiple signals can affect the estimation accuracy of the topography-correlated and long-scale components respectively. Thus, the strategy of simply using the XY-correlated method [26] and the T-correlated method [29] one by one to correct corresponding components may cause a less robust results. Therefore, a joint estimation method is proposed in next by considering that different tropospheric delay components have different spatial wavelength scales [25].

The TXY-Correlated Method
Based on the model in Equation (8), the TXY-correlated method is proposed in this section to correct the main tropospheric delays. The key steps of the TXY-correlated method are shown in Figure 1, where the procedures of the conventional linear methods, including the T-correlated method and the XY-correlated method, are also compared with the proposed method. The TXY-correlated method is marked with red solid box, while the conventional methods are marked with red dotted boxes.
Before tropospheric delay correction, the DEM error and orbit error have been removed from the unwrapped interferometric phase. The topography-correlated tropospheric delay is estimated by performing band-pass filtering and linear regression, which are the strategy adopted by both the T-correlated method and the proposed method. The long-scale tropospheric delay is estimated directly by applying linear regression with the XY-correlated method. However, this strategy is only suitable for small areas. If the area size exceeds the spatial scale range of long-scale tropospheric delay, the correction accuracy will become poor. Thus, block processing and multi-weight interpolation are introduced in the proposed method to estimate the local slope factors of long-scale tropospheric delay. Moreover, to improve the robustness of the estimation results, the multi-scale filtering and alternating iteration algorithm are employed in the proposed method.
In detail, firstly, multi-scale filtering is performed to decouple the topography-correlated and long-scale components since the topography-correlated tropospheric signal is present at all wavelength scales [31], while the latter is of around 10 km scales [8]. The band-pass filter is applied to select a spatial frequency band which is relatively sensitive to the topography-correlated tropospheric signal for estimating the global transfer factor K k ∆φ , while the low-pass filter is applied to select the spatial frequency band sensitive to the long-scale one for estimating the local slope factors A k i and C k i . Then, separate strategies are adopted to estimate the topography-correlated and long-scale tropospheric delays. To obtain the global transfer factor K k ∆φ , linear regression can be adopted directly, while for the local slope factors A k i and C k i , a block processing method [20,25] is applied first. The studied area is divided into multiple blocks by moving a small window from the left bottom to the right up corner of the full interferogram. To ensure consistency between adjacent blocks, the overlap ratio between blocks is set to 50%. After the local slope factors of each block are estimated with linear regression, they are all multiple-weighted and interpolated to each PS point of interferogram. The final derived slope factors of all PS points can be written as where A and C denote the N × 1 slope factor vectors to be estimated, N indicates the number of PS on the full interferogram, A ′ and C ′ denote the M × 1 vectors of estimated slope factors over each block, M indicates the number of block, and W represents the N × M weight matrix constructed by combining the average standard errors of the estimated slope factors of all blocks and the distance between blocks and PS points, with where G and S denote the N × M matrix of the Gaussian distribution based on distance and the transformed weights of the standard error, respectively, ⊙ is the symbol for Hadamar product, and s j is the average estimated standard error of slope factors over the j-th block. The block processing guarantees spatial correlation of the long-scale tropospheric delay, and the multi-weight interpolation guarantees consistency between PS points of different blocks. Thus, the strategy for estimating the long-scale component can be applied over relatively large areas.
After applying the multi-scale filters, it is still not possible to completely separate the long-scale and topography-correlated tropospheric signals. Therefore, to improve the robustness of estimated results, an alternating iterative algorithm is proposed to optimize the joint model. The alternations are performed between the long-scale tropospheric delay ∆φ long and the topography-correlated tropospheric delay ∆φ t−c . The specific steps are listed as follows: Step 1: Estimate initial values for each component of tropospheric delays based on the corresponding strategies, ∆φ Step 2: Subtract ∆φ (0) long from the original interferometric unwrapped phase, estimate the topography-correlated tropospheric delay again, and use the estimation as the updated result, ∆φ Similarly, subtract ∆φ (0) t−c from the original interferometric unwrapped phase, estimate the long-scale tropospheric delay again, and use the estimation as the updated result, ∆φ where ∆φ ′ unwrap is the unwrapped interferometric phase after removing the DEM error and the orbit error. T{•} and XY{•} are the corresponding processing strategies.
Step 3: Repeate step 2 and compare the updated result with the last estimation until the convergence condition is satisfied: and where ε is a small constant, which can be set empirically and 0.001 is used here for a satisfactory result.
Once the convergence condition is satisfied, the iteration ends.
To avoid affecting the estimation of other tropospheric delay components, the temporal low-pass filter is employed at the last step after the long-scale and topography-correlated delays have been corrected to reduce the effects caused by the short-scale tropospheric delay.

Selection of Spatial Filter Bandwidth
The spatial filter applied in the TXY-correlated method is conducted by a two-dimensional Gaussian function convoluted with the signal. It is effectively used for spatial phase filtering in StaMPS [35]. The width of the Gaussian convolution kernel, also known as the bandwidth of the spatial filter, represents the spatial scale. The low-pass images are generated by applying the Gaussian filter with different spatial scales. The band-pass images are generated by applying the Gaussian filter and taking the difference between two different scales.
The interferometric phases are comprised of multiple components, such as deformation, residual orbit error and noise. Different components have different multi-scale dependent spectrums. For example, ground deformation such as tidal loading [37], tectonic slow slip [6], is considered to be large-scale, and the land subsidence is often regarded as small-scale. Thus, an appropriate bandwidth for the employed spatial filters is not easy to obtain. Moreover, as different interferograms include different deformation features, atmospheric conditions, residual DEM errors and orbit errors, the spatial bandwidths where interested signals are present may be different. Therefore, different bandwidths should be selected for different interferograms to perform the spatial scale filtering. In this paper, a statistical analysis method is introduced to select the bandwidths based on the estimated standard error of slope factors after linear regression. The key procedures is presented in Figure 2.  Firstly, a series of spatial bandwidths are set up with their range from the resolution to the spatial size of the area, where the lower limit of bandwidth is restricted by the pixel resolution and the upper limit is restricted by the spatial extent of the area. For the topography-correlated tropospheric component, band-pass filtering with different bandwidths is performed on the interferograms and topography in space, and then the interferometric phase and topography are fitted linearly. For the long-scale tropospheric component, low-pass filtering with different bandwidths is performed on the interferograms in space, and then the interferometric phase and azimuth, slant coordinates are fitted linearly. Next, the estimated standard error of slope factors is calculated after linear regression, and the relationship map between bandwidth and estimated standard error for each interferogram is then obtained. When the estimated standard error on slope factors of the full interferogram reaches the minimum value, the corresponding bandwidth is selected as the optimum one.

Study Area and Dataset Used
Beijing is located in the north of China (39.4 • N-41.6 • N and 115.7 • E-117.4 • E) with a total area of 16,410.54 km 2 . The mountainous area of Beijing covers 10,200 km 2 , accounting for 62% of the total area, surrounding the Beijing plain from the southwest to the northeast, and the plain area is 6200 km 2 , accounting for 38% of the total area. The elevation of the Beijing plain is 20 m to 60 m, and the mountain ranges from 1000 m to 1500 m high. Beijing has a monsoon-influenced semi-arid and semi-humid continental climate, where the distribution of precipitation season is very uneven, and 80% of the annual precipitation is concentrated in the summer months of June, July and August. There are two alluvial-pluvial fans of the Yongding River and the Chaobai River on the Beijing plain [38].
The studied areas are located in Shunyi District at the northeast of Beijing, where targets with different scattering properties are distributed, such as the man-made and natural targets, which are marked by yellow boxes in Figure 3. One studied area is the non-urban area containing natural targets, such as river, lawn and bare rock, with a ground extension of 4.6 km × 3.6 km and the highest topography of 56 m. It is crossed by the Chaobai River and contained in one of the alluvial-pluvial fans, where the land subsidence issues are obvious [39,40]. The other one is the urban area belonging to the Renhe Town, containing man-made targets, such as complex buildings, with a ground area of 4.0 km × 3.2 km and the highest topography of 62 m. A stack of 25 TerraSAR-X/TanDEM-X images acquired from the ascending tracks are provided for the 2012 to 2015 period by the German Aerospace Center, which has a spatial resolution of 3 m on the ground [41]. Before tropospheric delay correction, the StaMPS processing [42] is used to generate the unwrapped interferograms. The detailed flow chart has been provided in Figure 1 marked by black dotted boxes. The external digital elevation model (DEM) applied for topography contribution removal is obtained from the Shuttle Radar Topography Mission (SRTM) DEM with a spatial resolution of 30 m. The 14th SAR image acquired on 10 October 2013 is selected as the master image and twenty-four interferograms are obtained with a maximum temporal baseline of 682 days, and perpendicular baseline ranges from −177 m to 315 m, shown as Figure 4. To select PS candidate pixels, the dispersion threshold in the amplitude dispersion index method [34] is set to 0.4. The PS points are finally selected based on the PS probabilistic model proposed by Hooper [34,35]. Once the PS points have been selected, the 3D phase unwrapping, DEM and orbit error correction can then be conducted and the unwrapped interferograms are corrected by the proposed method in this paper.

Results
The TXY-correlated method is applied to correct the major tropospheric delays of 24 TerraSAR-X/TanDEM-X unwrapping interferograms. The existing phase-based methods have been validated in previous research [26][27][28][29][30]. Thus, the TXY-correlated method modified according to the existing phase-based methods can be verified in comparison with the conventional method which combines the XY-correlated and the T-correlated methods. For fair comparison, all preprocessing steps are the same, and the conventional method also considers the correction of three tropospheric delay components. In the conventional method, the topography-correlated component is corrected firstly by using the T-correlated method, and then the long-scale component is corrected by applying the XY-correlated method, while for the proposed method, the topography-correlated and long-scale components are corrected simultaneously. Similarly, both the conventional and proposed methods correct the short-scale component by applying the temporal low-pass filter. The accuracy of the proposed method is evaluated by the residual effects of tropospheric delays and spatio-temporal standard deviations of time-series deformation after correcting the tropospheric delays, compared with the conventional method. Moreover, the universality of the proposed method is validated by the representative study areas, the non-urban area and the urban area.

Estimated Tropospheric Delays
The estimated tropospheric delays for each interferogram in the Chaobai River site and Renhe Town with the TXY-correlated method are shown in Figures 5 and 6. It is obvious that the PS points in urban area are denser than in non-urban area.    Figures 5a and 6a show the spatial distribution of estimated topography-correlated tropospheric delay in the Chaobai River site and the Renhe Town, respectively. The estimated topography-correlated delay varies from −0.75 rad to 0.75 rad in the Chaobai River site, and −0.63 rad to 0.79 rad in the Renhe Town, accounting for a relatively small portion of the tropospheric delays. Figures 5b and 6b provide the spatial distribution of estimated long-scale tropospheric delay, which varies from −3.10 rad to 3.34 rad in the Chaobai River site, and −10.02 rad to 9.90 rad in the Renhe Town. It is obvious that the long-sclae tropospheric delay is dominant in the tropospheric delays. The topography-correlated component is relatively small relative to the long-scale component. The reason is that both the Chaibai River area and Renhe Town with an average elevation of 30 m belong to the plain region. The topography-correlated tropospheric delay is often more obvious in regions with significant topography, such as mountainous areas [43], while the long-scale tropospheric delay mainly exists in flat areas, such as the plains [6].
When the absolute difference value between the updated result and the last result is less than 0.001, the iteration terminates and the number of iterations no longer increases. The convergence analysis is shown in Figure 7. Through the comparison of results in the study of both areas, the convergence speed in the Renhe Town is faster than in the Chaobai River site. The number of iterations reaches 6 in the Chaobai River site, while in the Renhe Town, the number of iterations is only 3. The reason is that the PS points in urban area are denser than in non-urban area. The increase of PS density will enhance the reliability of results and accelerate the convergence of iterative algorithm. [rad] Trend of convergence [rad] Trend of convergence

Residual Tropospheric Delays
As mentioned by the three-dimensional joint model, the magnitude of topography-correlated tropospheric delay can be expressed by the absolute value of corresponding slope factor K ∆φ . Similarly, the magnitude of long-scale tropospheric delay can also be expressed by the absolute value of corresponding slope factors A and C. As a result, the residual tropospheric delay effects could be represented by the average absolute magnitude of the residual slope factors. The residual amount of tropospheric delay phase is related to the accuracy of the tropospheric delay correction method. The smaller the residual tropospheric delay, the better the correction accuracy.

The Chaobai River Site
The average of residual tropospheric delay effects in the Chaobai River site are listed in Table 1. The residual topography-correlated and long-scale tropospheric delay effects are shown in Figures 8 and 9, respectively. In general, compared with the conventional method, the residual effects of tropospheric delays by the TXY-correlated method are smaller. For the topography-correlated component, it can be observed in Figure 8 that the TXY-correlated method performs better. Taking the 23 August 2015 interferogram as an example, the TXY-correlated method shows the reduction in slope factor K ∆φ of 0.09 rad/m, while the conventional method shows the reduction of 0.06 rad/m. However, on average, the residual slope factor K ∆φ after the conventional and TXY-correlated methods is 0.12 (±0.0028) rad/m and 0.11 (±0.0027) rad/m, as shown in Table 1. The proposed method improves the correction rate for slope factor K ∆φ with only 6.25% compared to the conventional method, which indicates the proposed method does not significantly improve the estimation of topography-correlated tropospheric delay in the Chaobai River site. For the long-scale component, since the original global model parameters are modified as local parameters, the TXY-correlated method shows the higher correction capability than the conventional method, as shown in Figure 9. Especially for the slope factor |C|, the proposed method shows the maximum reduction of 1.38 rad/km in the 22 January 2012 interferogram, while the conventional method only shows a maximum reduction of 0.35 rad/km in the 21 July 2015 interferogram. The average of |C| is reduced to 0.18 (±0.0991) rad/km, with a reduction rate of 75.8% after the TXY-correlated method, while using the conventional method, the average of |C| becomes 0.58 (±0.0990) rad/km, with a reduction rate of only 21.1%. The proposed method improves the correction rate for slope factor C with 54.7% compared to the conventional method. Different from the slope factor |C|, the performance of the proposed method for the slope factor |A| is slightly better than that of the conventional one. On average, the residual slope factor |A| after the conventional and TXY-correlated methods are 0.33 (±0.1294) rad/km and 0.31 (±0.1296) rad/km. The proposed method improves the correction rate for slope factor A with only 5.1% compared to the conventional method. However, for spaceborne InSAR, the slant range is usually relatively large, resulting in a larger proportion of tropospheric delay phase corresponding to the slope factor |C|, which indicates most of the long-scale component can be corrected. Although both the TXY-correlated and conventional method have a few overcorrected phenomena sometimes, such as the slope factor K ∆φ in the 1 July 2014 interferogram, |A| in the 18 May 2014 and 23 August 2015 interferograms, the proposed method has less error caused by overcorrection than the conventional method.

The Renhe Town
The average residual tropospheric delay effects in the Renhe Town are listed in Table 2. The residual topography-correlated and long-scale tropospheric delay effects are shown in Figures 10 and 11, respectively. Compared with the conventional method, the TXY-correlated method performs better and has smaller overcorrection rate. For the topography-correlated component, there are three interferograms overcorrected after performing the conventional method, including the 11 May 2012, 1 July 2014 and 23 August 2015 interferograms, while only one interferogram on 1 July 2014 is overcorrected after performing the TXY-correlated method. For the long-scale component, there are two interferograms overcorrected with the conventional method, i.e., the 20 May 2013 interferogram for slope factor |A| and the 18 May 2014 interferogram for slope factor |C|, while none of interferograms is overcorrected with the TXY-correlated method. Similar to the results in the study of Chaobai River, the tropospheric delays associated with slope factor |C| can be removed almost completely by the TXY-correlated method in the Renhe Town. The average of |C| is reduced to 0.07 (±0.0206) rad/km, with a reduction rate of 89.0% after the TXY-correlated method, while using the conventional method, the average of |C| becomes 0.23 (±0.0207) rad/km, with a reduction rate of only 64.1%. The proposed method improves the correction rate for slope factor C with 24.9% compared to the conventional method. Different from the results in the study of Chaobai River, the correction ability for tropospheric delays related to slope factor K ∆φ is improved by the TXY-correlated method in the Renhe Town. The average of K ∆φ is reduced to 0.021 (±0.00044) rad/m, with a reduction rate of 61.8% after the TXY-correlated method, while using the conventional method, the average of K ∆φ is 0.031 (±0.00045) rad/m, with a reduction rate of only 43.6%. The proposed method improves the correction rate for slope factor K ∆φ with 18.2% compared to the conventional method. However, the proposed method improves the correction rate for slope factor A with only 2.9% compared to the conventional method, which indicates the performance of correcting tropospheric delays related to slope factor |A| by the TXY-correlated method is slightly improved compared with the results obtained by conventional method. Moreover, the uncertainties of the estimated results in area of Renhe Town is much smaller than in area of Chaobai River site, which may be associated with the density of PS points in SAR image. The denser samples, the smaller estimation uncertainty.   The green squares represent the correlation A and C between the unwrapped interferometric phase and SAR coordinate before tropospheric delay correction. The blue circles and the red triangles represent the residual values after correction by the TXY-correlated method and the conventional method respectively. The vertical segment on each pattern represents the uncertainty of the estimation results.

Deformation Estimation Accuracy
The aim of correcting tropospheric delays for InSAR is to retrieve accurate deformation information. After removing the tropospheric delays, the remaining phase in interferograms is mainly ∆φ de f o which is the contribution caused by temporal displacement of targets. The correction accuracy of tropospheric delays will affect the estimation accuracy of ground deformation. Thus, the TXY-correlated method and the conventional method are compared in terms of the monitoring accuracy of ground deformation to demonstrate the performance of the proposed method. Lithology provides the geological background for land subsidence, so the spatial variations are similar in small regions over short periods of time, which means the ground deformation has high correlation in space and time. Thus, both the spatial standard deviations of deformation in small regions and the variation of spatial deformation standard deviations in time are small. To obtain the statistical results of time series LOS displacement, a small region is selected on interferograms at first. Then, a statistical analysis for ground LOS displacement of each interferogram after correcting tropospheric delays is made to calculate the corresponding mean values and standard deviations in the small region, so that the average ground LOS displacement curve varying with time can be obtained. Finally, the statistical results of time series LOS displacement over small regions are compared between the TXY-correlated and conventional methods.

The Chaobai River Site
The results of reconstructed mean deformation rate of PS points in the Chaobai River site are shown in Figure 12. The statistical results of time series LOS displacement over the small region R1 marked with red rectangle in Figure 3, which belongs to the Chaobai River bridge, are compared. The date of 22 January 2012 is set as the starting point when the LOS displacement began to change. According to Figure 13, the TXY-correlated method has achieved a better result. The average of the spatial standard deviations of time series LOS displacement are 0.34 mm/year and 0.21 mm/year after correction with the conventional method and the TXY-correlated method, respectively, which means the ground deformation estimated after applying the TXY-correlated method has higher spatial correlation. The standard deviations of the first observation on 22 January 2012 and the last observations on 23 August 2015 are 0.009 mm/year and 0.67 mm/year for the conventional method, and 0.006 mm/year and 0.39 mm/year for the TXY-correlated method, which indicates that the temporal variation of deformation standard deviation is smaller after using the TXY-correlated method. Thus, the ground deformation estimated after the TXY-correlated method has higher temporal correlation. Moreover, since the studied region belongs to a subsidence area [39,40], the deformation results obtained after tropospheric delay correction using the TXY-correlated method are more consistent with the actual land subsidence.

The Renhe Town
The reconstructed mean deformation rate maps of PS points with the TXY-correlated and conventional method are shown in Figure 14, and the statistical results of time series LOS displacement over the small region R2, which belongs to the Jihui street, are shown in Figure 15.

Discussion
The TXY-correlated method is proposed to improve the correction accuracy of tropospheric delays for InSAR by modifying the existing phase-based methods. The non-urban and urban areas from the same InSAR track are analyzed by the TXY-correlated and conventional methods. The proposed method has more accurate and robust correction capability compared with the conventional one. The improvements by proposed method are shown in three aspects. Firstly, the estimated slope factors A and C associated with the long-scale component in the conventional method are global constants. It does not consider the spatial variation characteristics of tropospheric delays, and the spatial correlation characteristic of the long-scale tropospheric delay will decrease when the studied areas become larger. Thus, the two slope factors are estimated as local variables in the proposed method which are related to the PS points. Secondly, the estimation of long-scale and topography-correlated tropospheric delays in the conventional method will affect each other, while multi-scale filtering is performed to decouple the two tropospheric delay components in the proposed method. Thirdly, the long-scale and topography-correlated components are estimated independently in the conventional method, while in the proposed method alternating iteration is conducted to jointly estimate the two components, which can improve the robustness of estimated results.
The selection of spatial filter bandwidth is associated with the deformation features and atmospheric conditions in interferogram. Different interferograms have different deformation features and atmospheric conditions. Thus, the spatial filter bandwidths on interferograms observed at different time are inconsistent. For the spatial band-pass filter, the corresponding effects caused by different bandwidths are shown in Figure 16. It can be observed that the optimum bandwidths are different for the interferograms acquired at different times and places. Whether it is in the Chaobai River site or in the Renhe town, when the lower limit of bandwidth is less than 400 m, the results fluctuate greatly with the increase of upper limit. The reason is that if the lower limit is set too low, the probability of other signals, not the topography-correlated tropospheric delay signal, included in the bandwidth will increase. When the lower limit of bandwidth is greater than 400 m, the spatial standard deviations become stable with the increase of upper limit, and the greater the lower limit, the more stable the curve. Especially, when the lower limit of the bandwidth is greater than 2 km, the results for most interferograms are all stable, which can explain why the lower limit of spatial band-pass filter bandwidth is often set to 2 km in many studies of tropospheric delays [20,25,31]. As for the spatial low-pass filter, the effects caused by different bandwidths are shown in Figure 17. The optimum bandwidth is about 50 m to 100 m in both the Chaobai River site and the Renhe Town, which is close to the empirical value presented in StaMPS [35]. Furthermore, the optimum bandwidth of low-pass filters for different interferograms is basically identical. It is probably because the smaller scale components, such as orbit error, which affect the estimation of long-scale tropospheric delay phase, have been removed before.
The slope factors A and C of each PS point are obtained by interpolation of the slope factors of all blocks. Thus, the estimation accuracy of the slope factors for each block directly affects the final results. Since the studied area is divided into several blocks by moving a small window from the left bottom to the right up corner of the full interferogram, to ensure consistency between adjacent blocks, the moving step size of the window is set to half of its size. If the window size is too small, it may contain fewer PS points, resulting in an inaccurate estimation of parameters. If the window size is too large, the correlation characteristic of long-scale tropospheric delay within the region contained in the window could be reduced, which also leads to serious estimation error. Therefore, the selection of the window size depends on the property of the studied areas. In order to avoid the slope factors of any blocks become null, each block should contain at least 2 PS points. The estimated standard errors of slope factors A and C after linear regression in each block with different window sizes are calculated to select the optimum window size. The optimum window size for each interferogram is given by what minimizes the estimated standard errors of slope factors after linear regression. Effects of different window size to the estimated standard errors of slope factors A and C are shown in Figure 17. It can be observed that the size of 4.1 km may be the best choice as window width in the Chaobai River site, and the size of 2.6 km may be the optimum in the Renhe Town.
According to the comparison results of residual tropospheric delays, the TXY-correlated method has higher correction accuracy than the conventional method, especially for the long-scale tropospheric component. The proposed method has different correction abilities on different parameters. It has the greatest correction rate for the slope factor C associated with the long-scale tropospheric delay. The reason is that the proposed method modified the XY-correlated model with respect to long-scale tropospheric delay, where the original global model parameters are modified as local parameters. Besides, for spaceborne InSAR, the slant range is usually relatively large, and as a result, the spatial correlation characteristic of the long-scale component in the slant range direction may decrease. Thus, the correction ability of the proposed method for slope factor C is stronger than for slope factor A. Both studies in the Chaobai River site and the Renhe Town have the same observation that the improvement of the correction ability for the slope factor A is not obvious. The reason may be that the spatial correlation scale of tropospheric delays is relatively long in azimuth direction, and it is shorter in slant range over the studied areas. Thus, the global slope factor A calculated by the conventional method is similar to the local slope factor A acquired by the proposed method. The correction improvement for topography-correlated component is lower than that for long-scale component. The reason may be that the same linear model about topography-correlated component is used in the TXY-correlated and conventional methods, so that the correction ability can be improved only by the alternating iterative algorithm in the proposed method. However, it can be observed in Figures 8 and 10 that the correction ability of K ∆φ is obviously improved with the TXY-correlated method in the Renhe Town. Denser samples will benefit the parameters estimation by linear regression. Thus, the performance for correcting the slope factor K ∆φ associated with the topography-correlated component can be improved with the increase of PS density in the studied area.
Overall, although accurate tropospheric delays can be obtained by data-based tropospheric delay correction methods, such as weather models, GPS or multi-spectral observations, the TXY-correlated method has been demonstrated to be an effective solution when external weather data is unavailable. However, the linear model about topography-correlated component has a major problem that the estimation error increases with the decrease of topography [20]. Thus, the performance of the TXY-correlated method can be improved by modifying the model associated with topography-correlated tropospheric delay in the future.

Conclusions
In this study, a novel tropospheric delay correction method has been proposed for the correction of tropospheric delays by jointly modeling the three tropospheric delay components. Moreover, the long-scale tropospheric delay and the topography-correlated tropospheric delay are estimated simultaneously by modified strategies. This method has two main advantages: the correction accuracy is improved by using the scale filtering and alternating iteration algorithm rather than a simple superposition of conventional methods, and the size limitation of the studied area with the conventional phase-based methods is relaxed. It can be applied to larger regions with deformation and suitable for a variety of scenarios, including natural areas and urban areas where the atmospheric condition is affected by human activities.
The TXY-correlated method is tested on 25 TerraSAR-X/TanDEM-X images in the Chaobai River and the Renhe Town of Beijing Shunyi District, where natural scences and man-made targets are both present. For different study areas, the selection of spatial filter bandwidth and window size are different, and the number of iterations of the alterative algorithm is related to the number of PS points selected in studied areas. The tropospheric delay phase estimated by the TXY-correlated method in the Chaobai River site is between −3.85 rad and 4.09 rad, and it is between −10.65 rad and 10.69 rad in the Renhe Town. Compared with the conventional method combined by the XY-correlated method and the T-correlated method, the proposed one has achieved a better delay correction performance. Moreover, the monitoring accuracy of deformation is improved due to greater reduction of spatio-temporal standard deviations for time-series LOS displacement after performing the proposed method.