Continued Monitoring and Modeling of Xingfeng Solid Waste Landfill Settlement, China, Based on Multiplatform SAR Images

Abstract

Continued settlement monitoring and modeling of landfills are critical for land redevelopment and safety assurance. This paper adopts a MTInSAR technique for time-series monitoring of the Xingfeng landfill (XFL) settlement. A major challenge is that the frequent and significant settlement in the initial stage after the closure of landfills can affect the coherence of interferograms, thus hindering the monitoring of settlement by MTInSAR. We analyzed the factors that can directly affect the temporal decorrelation of landfills and adopted a 3D phase unwrapping approach to correct the phase unwrapping errors caused by such deformation gradient. SAR images from four platforms, including 50 Sentinel-1A, 12 Radarsat-2, 4 ALOS-2, and 2 TerraSAR-X/TanDEM-X images, are collected to measure the settlement and thickness of the landfill. The settlement accuracy is evaluated by a cross-evaluation between Radarsat-2 and Sentinel-1A that have similar temporal coverages. We analyzed the spatial characteristics of settlement and the relationship between the settlement and thickness. Further, we modeled the future settlement of the XFL with a hyperbolic function model. The results showed that the coherence in the initial stage after closure of the XFL is primarily affected by temporal decorrelation caused by considerable deformation gradient compared with spatial decorrelation. Settlement occurs primarily in the forward slope of the XFL, and the maximum line-of-sight (LOS) settlement rate reached 0.808 m/year from August 2018 to May 2020. The correlation between the settlement and thickness is 0.62, indicating an obvious relationship between the two. In addition, the settlement of younger areas is usually greater than that of older areas.

1. Introduction

Landfills are a common destination for the disposal of municipal solid waste (MSW) and are a valid method for the recycling of land. For example, municipal parks and buildings can be constructed atop municipal solid waste landfills (MSWLs) when a certain stability threshold is reached [1,2,3]. However, an MSWL can experience long-term settlement as a result of mechanical compression and biodegradation after the closure [4,5,6], and uneven subsidence caused by such biodegradation may threaten the safety of the redevelopment [7,8]. A major threat is landslides [9,10]. Hence, monitoring and prediction of the long-term deformation of MSWLs are vital for land redevelopment and safety assurance.

In situ measurements obtained by instruments such as piezometers, inclinometers, and settlement plates are usually adopted for the monitoring of MSWLs. Although these methods are highly accurate, their spatial coverage is limited. In contrast, the differential interferometric synthetic aperture radar (DInSAR) provides wider coverage. It uses the phase difference between two acquisitions to measure the deformation in the line-of-sight (LOS) direction. Consequently, DInSAR has been widely used to monitor the deformation caused by land subsidence [11,12,13], landslides [14,15,16,17], volcanoes [18], and earthquakes [19,20]. However, two major factors can hinder the use of DInSAR for high-accuracy deformation monitoring: decorrelation and the atmosphere. To eliminate the influences of these two factors, multitemporal InSAR (MTInSAR) has been proposed for deformation time-series monitoring. Over the past few decades, many MTInSAR methods, such as PS-InSAR [21], SBAS-InSAR [22], SqueeSAR [23], StaMPS [18,24], TCPInSAR [25], and GMTInSAR [26], have been proposed for millimeter-level deformation monitoring in urban and mountainous areas.

To continuously monitor MSWLs with SAR images, different strategies have been proposed in recent years and can be divided into three categories. The first category is based on the amplitudes of SAR images to identify and track the areas of landfills undergoing change [27]. However, these techniques cannot obtain surface deformation. The second category uses the interferometric phase to retrieve the deformation of landfills with a small temporal baseline. For example, volumetric changes were obtained by DInSAR-derived DEMs [28,29], and deformation time-series were calculated by NSBAS [30], SBAS [31,32], SqueeSAR [33], PSInSAR [34], and multiple aperture interferometry (MAI) [31] approaches. The third category adopts ground-based InSAR (GBInSAR) to generate deformation time-series with a high temporal resolution [2]. Although the methods in this category evaluate the detectability of MSWLs with MTInSAR, these studies were conducted on MSWLs that were quite far behind the closure phase, in which the deformation rate is usually less than 5 cm/year; as a result, the deformation occurring in the initial stage after closure is missed. Nevertheless, it is crucial to understand the deformation after closure to help model MSWLs. Fortunately, the free distribution of Sentinel-1A/B SAR images enables the daily deformation monitoring of MSWLs. However, the deformation monitoring of MSWLs is different from the traditional deformation monitoring in natural or man-made environments, such as bare land, buildings, and other mixed scatterers, as the scattering properties of MSWLs surfaces are more variable during their life cycle. MSWLs are affected by frequent dumping activities in the operational phase, and vegetation growth occurs in conjunction with considerable deformation in the post closure phase. As a result, temporal decorrelation is the predominant factor directly affecting the coherence of interferograms [1,35,36]. However, the factors leading to the temporal decorrelation of MSWLs have not been thoroughly analyzed thus far. Generating the features of coherence of MSWLs can give guidance on the data processing of MTInSAR and can improve the accuracy of deformation.

In addition to the continued monitoring of settlement, modeling of MSWLs is also important for prediction. The settlement of MSWLs can usually be divided into three stages: initial settlement, primary settlement, and secondary settlement [6,37,38]. The first stage is caused by an external load, while the second stage is due to the dissipation of pore water and gas and biological decomposition of solids, and the third stage is caused by the creep of the refuse skeleton and biological decay. Among these three stages, the third stage contributes the most to the settlement [39,40]. To analyze the settlement of MSWLs, four kinds of settlement models have been proposed: soil mechanics-based models, rheological models, biodegradation-induced settlement models, and empirical models [4,6,41]. The first three model categories are usually restricted by the need to acquire in situ measurements of various parameters, such as the void ratio, overburden pressure, compressive stress, and density of biodegradable solids, which cannot be measured over wide areas. In contrast, empirical models, such as the logarithmic function, power function, hyperbolic function, and multiple linear regression function, are often selected to analyze the settlement characteristics of landfills. Among these models, the hyperbolic function has been applied to the analysis of 9 landfills with different fill ages [41] and the settlement of soft soil [32,34,42,43]. Hence, the hyperbolic function is used to predict the future settlement of MSWLs in this study. This work aims to analyze the existed data processing problems in the long-term monitoring of MSWLs and improve the accuracy of time-series deformation, analyze the settlement mechanism of MSWLs, and predict the future settlement. Studying these issues can provide an important indication for uneven subsidence and safety assurance and can play a significant role in the guiding of the land redevelopment.

In this study, the continued settlement of MSWL is calculated based on a 3D phase unwrapping MTInSAR technique. The coherence in the MSWL case is analyzed to guide the selection of SAR images and interferometric pairs. The spatial characteristics of the deformation are analyzed to derive the relationships between the thickness and settlement, and a hyperbolic function is adopted to assess the future stability of the MSWL.

2. Study Area and Datasets

2.1. The Xingfeng Landfill

The Xingfeng landfill (XFL), the largest landfill in China, is located southwest of Guangzhou city, Guangdong Province. Some residential areas and the Guanghe highway that connects the cities of Guangzhou, Huizhou, and Heyuan are distributed around the landfill. The area and the design capacity of the XFL are approximately 917 thousand m² and 42.85 million m³, respectively [3,44]. The landfill started to receive MSW in August 2002 with a design life of 25 year, and the XFL was divided into seven districts from August 2002 to December 2018, as shown in Figure 1 [45]. Districts 1–5 began operating in August 2002, July 2003, September 2004, and June 2006 and closed in June 2012, while Districts 6 and 7 were opened in May 2012 and July 2015, respectively. The daily disposal capacity and the maximum height of the XFL are approximately 3000 tons and 180 m, respectively. A DEM (5 m resolution) of the XFL derived from TerraSAR/TanDEM-X (TSX/TDX) on July 2018 before the closure is shown in Figure 1b, at which time the maximum height was 174 m. The capacity of the XFL in 2014 was 11,000 t/day, nearly four times the design capacity. Hence, a vertical expansion strategy has been utilized since March 2017, and MSW has been compacted since October 2017 [8]. These expansions are mainly located in Districts 4 and 5 and part of District 3, as shown in Figure 1a. Meanwhile, a reinforced slope was built in March 2017 to provide support for the additional MSW. The composite of the fresh MSW including food (54.02%), plastic (19.89%), paper (12.26%), textile (4.62), wood (3%), glass (2.45%), cinder and dust (1.62%), metal (0.89%), and others (1.25%) [8].

2.2. Datasets

We collected 68 SAR images to obtain the deformation time-series of the XFL, including 50 C-band Sentinel-1A images from October 2018 to July 2020, 12 C-band Radarsat-2 images from April 2019 to December 2019, 4 L-band ALOS2 images from Jul. 2015 to Nov. 2017, and 2 TSX/TDX images were acquired on 24 October 2013 and 17 July 2018, as shown in Figure 2,

Table 1. Parameters of the ALOS2 and TSX/TDX pairs.

ALOS2 Image Pairs (Ascending)			TSX/TDX Pairs (Descending)
Image Pairs	Perpendicular Baseline (m)	Temporal Baseline (Days)	Image Pairs	Perpendicular Baseline (m)
20150715–20151202	−53.8	140	20131024–20131024	−169.2
20170715–20171115	−2.0	126	20180717–20180717	−252.4

, and Tables S1–S3 (See Supplementary Materials). The Sentinel-1A data are used to generate the long-term deformation time-series of the XFL, while the Radarsat-2 data are utilized to evaluate the accuracy of the Sentinel-1A results. TSX/TDX is used to obtain DEMs of the XFL with high resolution, and high accuracy from 2013 to 2018, and the ALOS2 data are used to analyze the decorrelation of the XFL.

In addition to the SAR images mentioned above, we also collected time-series optical images from Google Earth, as shown in Figure 3. The changing areas of the surface were generated by visual interpretation. Anthropogenic activities were almost complete by October 2018, which is consistent with the closure time of the XFL. Hence, we selected SAR images after 8 August 2018 to avoid the influence of anthropogenic activities on the coherence calculation.

3. Method and Data Processing

3.1. SBAS-InSAR with 3D Phase Unwrapping

SBAS-InSAR is a commonly used MTInSAR method that has been successfully applied for deformation monitoring in many fields. This technique was developed based on the traditional DInSAR technique, whose differential phase (${\phi }_{diff}$) is a combination of deformation (${\phi }_{def}$), residual topography (${\phi }_{topo}$), orbit (${\phi }_{orb}$), atmosphere (${\phi }_{atm}$), and systemic noise (${\phi }_{noi}$) components, as shown in Equation (1):

$${\phi }_{diff}={\phi }_{def}+{\phi }_{topo}+{\phi }_{orb}+{\phi }_{atm}+{\phi }_{noi}$$

(1)

Usually, SBAS-InSAR is executed based on the unwrapped phase of each differential interferogram referred to as the original phase of SBAS-InSAR. However, individual interferograms may contain phase unwrapping errors when the required point density and phase quality are not met. In this study, a 3D phase unwrapping method was adopted to obtain the original phase of SBAS-InSAR [46]. The phase unwrapping operation was executed in 2D spatial dimensions and 1D temporal dimensions. A flowchart of the technique is shown in Figure 4.

First, the interferometric pairs were selected by utilizing spatial and temporal baseline thresholds. According to the coherence analysis in [1], the temporal baseline threshold is a vital factor that can directly affect the coherence compared with the perpendicular baseline threshold. In this paper, the thresholds for Sentinel-1A and Radarsat-2 were 175 m/24 days and 300 m/48 days, respectively, to account for the considerable deformation that occurred during the initial stage after the closure of the XFL. The spatial and temporal baselines of the selected interferograms are shown in Figure 5 and Table 1. In addition, multi-look operations for Sentinel-1A (4 looks in the range direction and 1 looks in the azimuth direction), Radarsat-2 (2 looks in the range direction and 3 looks in the azimuth direction), and ALOS2 (3 looks in the range direction and 8 looks in the azimuth direction) were adopted to achieve a compromise between resolution and denoising. It is worth noting that a multi-look operation with large looks may not be suitable for landfills undergoing large deformation, typically associated with phase aliasing in the differential interferograms (for further details, see Section 4.1). Traditional differential interferograms of each selected pair were generated after removing the topography and flattened phases. External DEMs were derived from an iterative TSX/TDX DEM generation method instead of the commonly used SRTM data due to the notable topographic variation. A schematic flowchart of the proposed method is illustrated in Figure 4. A small multi-look operation (2 × 2) follows a large multi-look operation (4 × 4) to help phase unwrapping in areas with an obvious topographic variation. The final DEMs were obtained after reaching the resolution and accuracy requirements, following a procedure by Du et al. [47]. High-quality points were selected based on the coherence calculated by an unbiased estimation method [48], and the threshold used in this case was 0.25. Next, the wrapped phases were unwrapped through a 3D phase unwrapping method, and a patch-based polynomial method was adopted to remove the orbit phase of each unwrapped phase [49]. Then, the residual topography was calculated by an iterative phase analysis method and removed from the original phase of SBAS-InSAR [24]. Subsequently, spatial and temporal filtering was performed to remove the atmospheric effects from each image, and a final map of the mean settlement rate was obtained by least-squares fitting of the deformation time-series. Finally, a settlement model was employed in conjunction with high-resolution optical images to analyze the spatial characteristics of the landfill.

3.2. Settlement Model

As mentioned above, we selected the hyperbolic function to model the future settlement of the XFL, as shown in Equations (2) and (3):

$$\Delta D=\frac{\Delta t}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$v_0$}\right.+\raisebox{1ex}{$ \Delta t$}\!\left/ \!\raisebox{-1ex}{$D_{ult}$}\right.}$$

(2)

$$\Delta t=t_i-t_0, \Delta D=D_i-D_0$$

(3)

where $\Delta D$ and $\Delta t$ are the change in deformation and the time interval between an arbitrary time ($t_i$) and the beginning time ($t_0$), respectively; $v_0$ is the initial settlement rate at the time $t_0$; and $D_{ult}$ is the ultimate settlement at the time $t_{\infty }$. To calculate $v_0$ and $D_{ult}$, Equation (2) can be rewritten as a relationship between $\frac{\Delta t}{\Delta D}$ and $\Delta t$, as shown in Equation (4):

$$\frac{\Delta t}{\Delta D}=\frac{1}{v_0}+\frac{\Delta t}{D_{ult}}$$

(4)

The parameters $v_0$ and $D_{ult}$ can be obtained after calculating $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$v_0$}\right.$ and $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D_{ult}$}\right.$ with a linear least-squares method. Additionally, the time corresponding to the ultimate settlement ($t_{\infty }$) can be calculated if $D_0$ and $t_0$ are given, as shown in Equation (5):

$$t_{\infty }=t_0+\frac{D_{ult}·\left(D_{ult}-D_0\right)}{v_0·D_0}$$

(5)

4. Results

4.1. Coherence Analysis of the XFL

The coherence of interferometric pairs for a single platform was affected mainly by temporal decorrelation, spatial decorrelation, and volumetric decorrelation [50]. In the case of an MSWL, spatial decorrelation was not sensitive to coherence compared with temporal decorrelation [1]. Moreover, volumetric decorrelation was not obvious in the initial stage due to the sparse coverage of vegetation. Thus, we analyzed the relationships between coherence and the temporal baseline, platform type, and processing procedure.

We calculated the mean coherence of the XFL for the one-revisit and two-revisit time of each SAR platform (for example, 12 days and 24 days for Sentinel-1A and 24 days and 48 days for Radarsat-2). The results are plotted in Figure 6. The mean coherence of Sentinel-1A for pairs with one-revisit time was usually larger than 0.3 and increased with time. The coherence of the Sentinel-1A interferometric pairs with a 24-day interval also increased with time, but the overall coherence was lower than that of the pairs with a 12-day interval. The smallest value of the mean coherence for Sentinel-1A was 0.28. Similarly, the coherence of Radarsat-2 increased with time for the pairs with 24-day and 48-day intervals. However, the coherence was lower than that of Sentinel-1A. For example, the maximum coherence for the 24-day and 48-day pairs of Radarsat-2 were 0.32 and 0.22, respectively. We also analyzed the coherence of the ALOS2 data; however, the smallest temporal baseline was 126 days due to the discontinuity of the archived datasets. Thus, we selected two pairs with 126-day and 140-day intervals and plotted the coherence; the results are shown in Figure 7a,b. In contrast to our expectations, the mean coherence for both pairs was 0.11, resulting in serious decorrelation despite the better penetration capability of the long-wavelength ALOS2. Therefore, the ALOS2 images were not adopted for subsequent processing. In addition to the temporal baseline, we also plotted the perpendicular baseline of each selected image pair, as shown in Figure 6, revealing that the mean coherence of both Sentinel-1A and Radarsat-2 was not sensitive to the perpendicular baseline, as was demonstrated by [1].

In addition to the coherence, we plotted the wrapped phases of the differential interferograms. As shown in Figure 7c, two obvious fringes are observed in a 12-day differential interferogram, representing 5.6 cm of settlement during a one-revisit time of Sentinel-1A. Similarly, a 24-day differential interferogram from Radarsat-2 with five obvious fringes is shown in Figure 7f, which further confirms that the XFL exhibited considerable settlement. Thus, in addition to the temporal decorrelation caused by the changes in surface features, a large deformation gradient may have been a major cause of the temporal decorrelation of the XFL. Moreover, we analyzed the multi-look operations and filtering window size for the differential interferograms. We found that larger multi-look operations, namely, 8 looks and 2 looks in the range and azimuth directions, respectively, can directly destroy fringes when considerable deformation is occurring, as shown in Figure 7d. In contrast, small multi-look operations can not only maintain the features of the fringes but also reduce phase noise. Similarly, filtering with large window size (64 × 64 pixels) can affect the completeness of fringes in the differential interferograms compared with a 16 × 16 pixels window, as shown in Figure 7e.

4.2. Deformation Time-Series

The mean settlement rate in the LOS direction of Sentinel-1A from August 2018 to May 2020 was generated with the processing procedures described in Section 3.1, and the results are shown in Figure 8. Two obvious subsidence basins were observed, marked as M and N in Figure 8. Area M spanning 2.27 km² was the area of the XFL, a settlement rate exceeding 0.25 m/year. The maximum settlement rate within area M was 0.808 m/year. The inset shown in Figure 8b is an enlarged view of area M, showing inhomogeneous settlement in the landfill. In addition, another large subsidence area, area N, was also observed. Similar to area M, area N was located on reclaimed land and was active from August 2016 to June 2019 according to the Google Earth optical images presented in Figure 9. The maximum subsidence rate of area N was 0.258 m/year, and the subsidence area was approximately 1.65 km², with a settlement rate exceeding 0.1 m/year. However, the details regarding the subsidence in area N were beyond the scope of this paper.

To analyze the spatial characteristics of settlement occurring within the XFL, we extracted the mean settlement rate of the XFL on a 3D DEM background, as shown in Figure 10, revealing that the majority of subsidence was located on the forward slope of the XFL. A total of three profiles were selected to analyze the spatial characteristics of subsidence, as shown in Figure 11. Profile S1 runs southwest to northeast and passes through multiple subsidence basins, and the settlement rate profile is plotted in Figure 11b. Several basins with different magnitudes of the settlement were observed, indicating that subsidence basins formed individually across the slope. Profile S2 runs northwest to southeast and also passes through a subsidence basin, as shown in Figure 11c. Only one obvious basin with a width of 0.78 km was observed, indicating that the subsidence in the backward slope was gentle compared with that in the forward slope. We also plotted the settlement rate along the crest of the XFL (profile S3), as shown in Figure 11d. The settlement rate on the east side of S3 was more stable than that on the west side, indicating that the settlement in District 3 was stable.

Similarly, to analyze the features of the deformation time-series, three points located in areas with different magnitudes of the settlement were selected, and the corresponding deformation time-series are shown in Figure 11a. The deformation time-series along the margin of the landfill tend to be stable with a maximum deformation of 35.1 cm (P1). In contrast, the deformation time-series in the subsidence basins undergo continuous deformation with maximum subsidence values of 130.2 cm and 73.7 cm at P2 and P3, respectively.

4.3. Accuracy Evaluation

The relative accuracy was adopted in this study due to the lack of in situ measurements. We used SAR datasets from the two platforms with similar time spans to calculate the consistency. The time-span of Radarsat-2 ranged from 4 April 2019 to 24 December 2019, and that of Sentinel-1A ranged from 7 April 2019 to 27 December 2019. The geometry of the two platforms was different: ascending datasets with an incidence angle of 41.5° for Sentinel-1A and descending datasets with an incidence angle of 27.4° for Radarsat2. Hence, the LOS deformation can be decomposed to east-west (horizontal) and up-down (vertical) direction based on the assumption that there is no obvious deformation in the north-south direction [32,34,51,52,53,54], as shown in Equation (6):

$$\left(\begin{array}{c}{d}_{asc}\\ d_{dsc}\end{array}\right)=\left(\begin{array}{ll}cos{\theta }_{asc}& -cos{\alpha }_{asc}·sin{\theta }_{asc}\\ cos{\theta }_{dsc}& -cos{\alpha }_{dsc}·sin{\theta }_{dsc}\end{array}\right)·\left(\begin{array}{c}{d}_{U-D}\\ d_{E-W}\end{array}\right)$$

(6)

where ${\theta }_{asc}$ and ${\theta }_{dsc}$ are the incidence angles and ${\alpha }_{asc}$ and ${\alpha }_{dsc}$ are the satellite heading angles of ascending and descending platforms. $d_{U-D}$ and $d_{E-W}$ are the vertical and horizontal deformation.

Because of the resolution difference of the two platforms, the high-resolution results of Radarsat-2 were down-sampled to the medium-high-resolution results of Sentinel-1A. Then, valid pixels were selected for both platforms for the 2D deformation calculation. The settlement rates in the horizontal and vertical directions are shown in Figure 12e,d, respectively. The deformation in the vertical direction is more dominant than the horizontal displacement, except for some local deformation in the horizontal direction. Hence, horizontal displacement can be neglected, and the vertical deformation for two platforms are calculated with Equation (7):

$$V_{vertical}=\frac{V_{LOS}}{\cos \left(\theta \right)}$$

(7)

The vertical deformation of Radarsat2 and Sentinel-1A are shown in Figure 11a,b, respectively. The subsidence patterns of the two results are comparable with a maximum subsidence rate in District 6. A scatterplot in Figure 12c reveals the relationship between the two. A correlation coefficient of 0.84 indicates the high consistency between the Sentinel-1A and Radarsat-2 datasets. Moreover, we also plotted a scatterplot between the decomposed vertical deformation calculated by Equation (6) and vertical deformation calculated by Equation (7) of two platforms, as shown in Figure 12f. The correlation coefficients for Radarsat2 and Sentinel-1A are 0.97 and 0.94, which proves the feasibility of Equation (7) in this study.

4.4. Modeling of Settlement

Using the hyperbolic function introduced in Section 3.2, we can model the future vertical settlement of the XFL. In this study, the beginning time ($t_0$) was selected as 3 August 2018, the $D_0$ was selected as 0, and the final ultimate settlement $(D_{ult})$ of each pixel was calculated based on the deformation time-series and time interval, as shown in Figure 13a. It is noted that the time corresponding to $D_{ult}$ can also be calculated with a given $D_0 \left(D_0\ne 0\right)$ and $t_0$ through Equation (5).

The small ultimate settlement in Districts 1–3 indicates the stability of Districts 1–3 after a long-term settlement. In contrast, settlement exceeding 3 m was observed mainly in Districts 4–6, and part of District 7 in the forward slope. The maximum settlement in District 6 reached 29.45 m, which was 28.1% of the total thickness (nearly 105 m [45]), within the range of 20–60% [41]. In addition, we plotted the hyperbolic function model of the three selected points described in Section 4.2, and the results are shown in Figure 13c–e. Points P1, P2, and P3 were located in Districts 7, 6, and 4, respectively. The linearity of the fitted line for P2 indicates continuous settlement in District 6 with an ultimate settlement of 16.67 m, while the fitted lines for P1 and P3 represent a slowing trend with ultimate settlements of 0.93 m and 1.84 m, respectively. The correlation coefficient of P1–P3 was over 0.99, indicating the reasonability of the hyperbolic function model. Moreover, we plotted the slope map of the XFL on 17 July 2018 in Figure 13b, demonstrating that the relationship between the ultimate settlement and slope was not obvious.

To evaluate the accuracy of the settlement model, we calculated the accumulated deformation of MTInSAR and the settlement model between Aug. 2018 and May 2020 with a temporal coherence larger than a given threshold (0.6 in this case) [54], as shown in Figure 14a,b. The deformation pattern of the two methods was approximate. The difference between these two results is shown in Figure 14d, with an RMSE of 6.1 cm. The difference with larger values was mainly located in Districts 6 and 7. In addition to the accumulated deformation, we also calculated the RMSE of the time-series deformation between these two results, as shown in Figure 14c. The mean RMSE was about 3.9 cm, and the percentage of pixels whose RMSE was smaller than 4 cm is 82%. The RMSE of the three selected points were 0.92 cm, 4.78 cm, and 2.16 cm for P1, P2, and P3, respectively, as shown in Figure 13c–e. Pixels with RMSE larger than 4 cm were mainly located in areas with approximate linear deformations (such as P2 in District 6, as shown in Figure 13d) and some local horizontal deformations (Figure 12e) in Districts 6 and 7. The large RMSEs were probably caused by the complex deformation in the fresh part of MSWL that cannot only be modeled by the hyperbolic function, which was verified by Park and Tan [32,42]. Moreover, the hyperbolic function model needs enough observations to calculate its parameters [32]. Hence, the usage of the hyperbolic function model may be limited by the complex deformation and limited SAR images at the onset of MSWL settlement, and the continued datasets in the future can be used to help improve the accuracy of the model.

5. Discussion

5.1. Analysis of Coherence and Processing Procedures in the Initial Stage of MSWL after Closure

The coherence and processing procedures should be taken into consideration for the continued deformation monitoring in the initial stage of MSWL after closure. The coherence increases with time for both SAR images, i.e., Sentinel-1A and Radarsat2, in the initial stage of XFL, which is different from the feature that coherence has an obvious seasonality [1]. This may be caused by the decreasing water content of the MSWL over time [40] and the lacking vegetation in the initial stage. Additionally, the lower coherence of Radarsat-2, compared with Sentinel-1A, is possibly caused by its high resolution. That is, 2 m and 3 m in the azimuth and range directions, respectively. The higher the resolution is, the greater the speckle noise, which can directly affect the coherence of interferograms. Moreover, the completeness of deformation fringes is also a vital factor that can affect the phase quality. Multi-look operations with large looks and filtering operations with large windows can destroy fringes, leading to phase aliasing and wrong interpretation of deformation. Hence, a compromise should be given between resolution and denoising when high-resolution SAR images such as TerraSAR-X or Radsarsat-2 are adopted.

5.2. Relationship between Settlement and Fill Age

The fill age of a landfill is defined as the time span between the beginning of the measurement and the midpoint between the beginning and completion of filling [41]. The measurement times and fill age for each district are described in

Table 2. Fill age of each district in the Xingfeng landfill (XFL).

District Name	Beginning of Measurement	Beginning Time/Vertical Expansion	Completion Time/Vertical Expansion	Fill Age (year)	LOS Maximum Settlement Rate (cm/year)
1	August 2018	August 2002	June 2012	11.08	−31.38
2	August 2018	July 2003	June 2012	10.62	−28.65
3	August 2018	September 2004	June 2012	10.04	−22.55
4	August 2018	June 2006/November 2017	June 2012/October 2018	9.17/0.42	−55.65
5	August 2018	June 2006/November 2017	June 2012/October 2018	9.17/0.42	−62.30
6	August 2018	June 2012	October 2018	3.08	−80.80
7	August 2018	June 2015	October 2018	1.58	−56.40

. The landfill groups are divided by the fill age: in the young group, the fill age is <5 year, while in the old group, the fill age ranges from 5 year to 20 year [41]. In this study, Districts 1–5 are old, while Districts 6 and 7 are young. As we expected, obvious settlement is located primarily within the young districts (Districts 6 and 7) with fill ages of 3.08 year and 1.58 year, respectively. However, the settlement rates in Districts 4 and 5 are also obvious even though the fill age reaches 9.17 year; this is due mainly to the vertical expansion that occurred in November 2017 after the closure of these districts [8], resulting in a fill age of 0.42 year for parts of Districts 4 and 5. The settlement rates in Districts 1–3 are gentler than those in Districts 4–6, indicating that the rates are slowing down after 10 year of settlement. Therefore, the settlement of XFL is occurring predominantly in areas with young fill ages.

5.3. Relationship between Settlement and Landfill Thickness

To analyze the relationship between the settlement and thickness of the XFL, we generated two TSX/TDX DEMs on 24 October 2013 and 17 July 2018 with an iterative method [47]. The resolution of these two DEMs was resampled to the resolution of Sentinel-1A. A map of the thickness of the landfill between the two TSX/TDX DEMs is shown in Figure 15a. The thickness of Districts 1–3 was gentle due to their closure in June 2012, which was not within the time-span of two TSX/TDX DEMs. The landfill was more than 20 m thick, mainly in Districts 4–7, which was in accordance with the mean LOS settlement rates shown in Figure 15b. To clarify this relationship, we plotted the thickness and settlement rate in Figure 15c. The thickness of areas with a high settlement rate (<−0.4 m/year) was mainly in the range of 10 m to 20 m. A linear fit was also obtained by implementing a least-squares method, yielding $\mathrm{y}=-0.01181-0.10484\ast \mathrm{x}$. The correlation between these two factors was 0.62, indicating an obvious relationship between the thickness and settlement of the landfill, as was reported by [39]. The greater the thickness of the landfill is, the greater the settlement rate; however, the settlement rate in part of District 7 was gentle even though this part of the landfill was quite thick. This is probably due to the different types of MSW with different moisture contents, resulting in different biodegradation rates [8]. The settlement rates with the approximate thickness within the red dotted rectangle in Figure 15c are mainly located in Districts 4–7. The reason for this wide range distribution of settlement rate may be caused by the different dumping sequence (fill age) in the four districts during 24 October 2013 and 17 July 2018.

6. Conclusions

This paper employs SBAS-InSAR with 3D phase unwrapping to generate a deformation time-series of the XFL based on multiplatform SAR images. Data processing in the initial stage after the closure of XFL is analyzed to guide the correct selection of interferometric pairs in future work. The coherence of the landfill during the initial stage of closure is affected mainly by temporal decorrelation caused by considerable gradient deformation. The coherence of XFL increases with time for Sentinel-1A and Radarsat2 datasets. In addition, the size of the multi-look operation can affect the completeness of fringes, and the results indicate that a compromise should be found to balance the resolution and denoising. The settlement is located mainly on the forward slope of the XFL, with a maximum LOS settlement rate of 0.808 m/year. Moreover, the settlement mechanism of XFL is analyzed, and a hyperbolic function model is adopted to model the future settlement. An obvious relationship between the settlement and thickness of the XFL is found, and the correlation coefficient is 0.62. Additionally, prominent magnitudes of the settlement are observed in the young districts with relatively recent fill ages. The ultimate settlement of the XFL can reach 29.45 m in District 6.