Radar Interferometry Time Series to Investigate Deformation of Soft Clay Subgrade Settlement—A Case Study of Lungui Highway, China

Monitoring surface movement near highways over soft clay subgrades is fundamental for understanding the dynamics of the settlement process and preventing hazards. Earlier studies have demonstrated the accuracy and cost-effectiveness of using time series radar interferometry (InSAR) technique to measure the ground deformation. However, the accuracy of the advanced differential InSAR techniques, including short baseline subset (SBAS) InSAR, is limited by the temporal deformation models used. In this study, a comparison of four widely used time series deformation models in InSAR, namely Multi Velocity Model (MVM), Permanent Velocity Model (PVM), Seasonal Model (SM) and Cubic Polynomial Model (CPM), was conducted to measure the long-term ground deformation after the construction of road embankment over soft clay subgrade. SBAS-InSAR technique with TerraSAR-X satellite imagery were conducted to generate the time series deformation data over the studied highway. In the experiments, three accuracy indices were applied to show the residual phase, mean temporal coherence and the RMS of high-pass deformation, respectively. In addition, the derived time series deformation maps of the highway based on the four selected models and 17 TerraSAR-X images acquired from June 2014 to November 2015 were compared. The leveling data was also used to validate the experimental results. Our results suggested the Seasonal Model is the most suitable model for the selected study site. Consequently, we analyzed two bridges in detail and three single points distributed near the highway. Compared with the ground leveling deformation measurements and results of other models, SM showed better consistency, with the accuracy of deformation to be ±7 mm.


Introduction
The stability management of soft clay subgrade, is one of the main challenges in the field of highway engineering [1,2]. Man-made infrastructures built on soft clay subgrade are more prone to settlement and subsidence due to the geological characteristics of soft clay with high water content, high compression and weak strength. Therefore, long-term subsidence monitoring for highways built on soft clay subgrade is crucial to prevent potential hazards [3].
Differential Interferometric Synthetic Aperture Radar (D-InSAR) has been widely applied to detect and monitor ground surface movements caused by natural geological phenomena, e.g., . Each of these N interferometric pairs is processed by two-pass D-InSAR to identify the ground deformation by referring to a same stable reference point. The wrapped differential interferometric phase at the pixel with coordinates (x, y) in azimuth and range directions, in the corresponding ith interferogram, can be expressed as [11]: where λ is the radar wavelength (3.2 cm for X-band TerraSAR-X data); t A and t B represent the acquisition time for the two SAR images in corresponding ith interferometric pair; d(t B , x, y) and d(t A , x, y) are the line of sight (LOS) cumulative deformations at times t A and t B with respect to the time t 0 (which is assumed as a reference time). ∆ϕ topo i is the phase contribution related to DEM error, which is proportional to the perpendicular baseline of the interferometric pair and inversely proportional to the sine of incident angle. ∆ϕ res i (x, y) is the residual phase, including phase noise, atmospheric delay, and high-pass (HP) deformation component.
Here, the phase component ∆d = d(t B , x, y) − d(t A , x, y) is of the main interest. In MVM, the functional relationship between ∆d and the deformation parameters can be written as where r, q define the index of master image at time t A and slave image at time t B , respectively for i-th interferometric pair. v j defines the linear velocity for each temporal unit, it varies in different temporal units, so we need to estimate M number of vs as well as the unknown DEM error parameter (in total M + 1 unknown parameters) in N generated functions. To solve the singular mathematical problem, a Singular Value Decomposition (SVD) algorithm is used here to carry out and solve the unknown Equation (2) [33]. However, the drawbacks of MVM model are the deficiencies in the process of unknowns estimation, and the temporal discontinuity when the high quality interferometric temporal period cannot cover all M + 1 temporal points, which has been referred as temporal gap and analyzed in detail in [34]. The simulated functional relationship and corresponding temporal gaps are shown in Figure 1a,b.
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 23 thresholds, where satisfies the inequality ≤ ≤ ( ) . Each of these N interferometric pairs is processed by two-pass D-InSAR to identify the ground deformation by referring to a same stable reference point. The wrapped differential interferometric phase at the pixel with coordinates ( , ) in azimuth and range directions, in the corresponding ith interferogram, can be expressed as [11]: where λ is the radar wavelength (3.2 cm for X-band TerraSAR-X data); and represent the acquisition time for the two SAR images in corresponding ith interferometric pair; ( , , ) and ( , , ) are the line of sight (LOS) cumulative deformations at times and with respect to the time (which is assumed as a reference time). ∆ is the phase contribution related to DEM error, which is proportional to the perpendicular baseline of the interferometric pair and inversely proportional to the sine of incident angle. ∆ ( , ) is the residual phase, including phase noise, atmospheric delay, and high-pass (HP) deformation component.
Here, the phase component ∆ = ( , , )-( , , ) is of the main interest. In MVM, the functional relationship between ∆ and the deformation parameters can be written as where , define the index of master image at time and slave image at time , respectively for i-th interferometric pair.
defines the linear velocity for each temporal unit, it varies in different temporal units, so we need to estimate number of s as well as the unknown DEM error parameter (in total + 1 unknown parameters) in generated functions. To solve the singular mathematical problem, a Singular Value Decomposition (SVD) algorithm is used here to carry out and solve the unknown Equation (2) [33]. However, the drawbacks of MVM model are the deficiencies in the process of unknowns estimation, and the temporal discontinuity when the high quality interferometric temporal period cannot cover all M + 1 temporal points, which has been referred as temporal gap and analyzed in detail in [34]. The simulated functional relationship and corresponding temporal gaps are shown in Figure 1a,b. Another widely used linear model is PVM, which is proposed for Permanent Scatter InSAR (PSInSAR) technique. This method considers a fixed deformation rate which substitutes for all the Another widely used linear model is PVM, which is proposed for Permanent Scatter InSAR (PSInSAR) technique. This method considers a fixed deformation rate which substitutes for all the different v i in Equation (2), as the functional relationship shown in Figure 1c, thus the number of unknowns decreases in contrast to that of MVM.

Non-linear SBAS Models
Both MVM and PVM treat the deformation deriving process as linear temporally. Obviously, it may not be suitable for nonlinear displacement. Seasonal Model (SM) and Cubic Polynomial Models (CPM) are designed specifically for nonlinear deformations. The SM functional relationship can be written as [32]: where T defines 365 days for a year, t defines the accumulated time with respect to the reference time t 0 , thus the unknown periodical coefficients are s 1 , s 2 , s 3 . This model is widely used for observations over regions where there are significant seasonal variations. The functional relationship between displacement and time is shown in Figure 1d. The Cubic Polynomial Model (CPM) is proposed according to the Weierstrass fitting principle. The deformation component can be written as [28]: where the unknown deformation coefficients are c 1 , c 2 , and c 3 , respectively. Earlier study [28] demonstrated that this CPM model performs well with the residual phase within the region ∆ϕ res i < π. The simulated functional relationship between deformation and time for a single high coherence point of the 4 models aforementioned are illustrated in Figure 1, where x axis defines the 10 simulated temporal points (10 ordinal SAR acquisition time with reference to the first SAR image), y axis is the simulated deformation values.

Study Area
The selected structure built on soft clay subgrade in our study is Lungui highway, in Shunde District, Fuoshan City of Guangdong Province, in southern China. Figure 2 shows the study area features at different scales. As Figure 2a shows, the red rectangle outlines the spatial coverage of the selected TerraSAR-X images, while the green rectangle shows the sub-set for generating the interferometric results. Lungui highway, with a total length of approximately 12 km. It connects the Longzhou, Nanguo and Hejiu Highways, becoming one of the three major South-North lines in Fuoshan. The subgrade was designed to be 60 m in width, with soft clay being the primary rock-soil substrate along the route (see Figure 3a). This place locates in delta fluvial plain. The soil along this road consist of mucky clay, silty clay, and mealy sand. Mucky clay, as the main component, has the characteristic of high-water content, high compressibility, low mechanical strength and high sensitivity. The earthiness for this kind of soil is not homogeneous, thus prone to the geological disasters of slide skid, uneven subsidence, creep and soft clay seismic subsidence.
The characteristic of ground layer distribution is pseudo-viscosity plain fill and quaternary system of brand-new sea-land cross stratum. The surface quaternary is comprised of mealy sand and mucky clay. From the longitudinal distribution characteristic, the soft clay of the section north to Ronguite Bridge is mainly continuously distributed with silt, and the average thickness is 12.93 m~19.50 m; the section between Ronguite Bridge and Xiti Fouth Road, with the main soil of mucky clay and silty clay and thickness of 6.46 m~9.90 m; the section south to Xiti Fouth Road, with the main soil of mucky clay and thickness of 3.57 m~7.62 m. Figure 4 shows the transversal profile of test highway at the location LL'(showed as the red line in Figure 2b). The average soft clay thickness of this location is 4.5 m. The characteristic of ground layer distribution is pseudo-viscosity plain fill and quaternary system of brand-new sea-land cross stratum. The surface quaternary is comprised of mealy sand and mucky clay. From the longitudinal distribution characteristic, the soft clay of the section north to Ronguite Bridge is mainly continuously distributed with silt, and the average thickness is 12.93 m~19.50 m; the section between Ronguite Bridge and Xiti Fouth Road, with the main soil of mucky clay and silty clay and thickness of 6.46 m~9.90 m; the section south to Xiti Fouth Road, with the main soil of mucky clay and thickness of 3.57 m~7.62 m. Figure 4 shows the transversal profile of test highway at the location LL'(showed as the red line in Figure 2b). The average soft clay thickness of this location is 4.5 m.    It should also be noted from Figure 2 that the Lungui highway is surrounded with dense pounds and rivers. It locates close to three hydrological systems: Xi River, Rongui and Shunde Branch River. Therefore, Lungui highway is surrounded by plenty of pounds and lakes and large amount of silt. Most of these pounds are local cultivation. According to the geological material we collected, this area has extensive aquifers. The underground water level is influenced greatly by seasonal climate factors, with maximum levels varying by up to 1 m. Figure 2b is the radar amplitude image of the study area, with the highway region of interest outlined in the yellow rectangle. There are two main bridges along the highway, Rongguite Bridge and Anlite Bridge, which are analyzed in detail in Section 4. Figure 3b,d are the ground pictures of the two bridges, respectively.

SAR Acquisitions and Data Processing
In total 17 repeat-pass TerraSAR X-band Stripmap descending images, provided by the German Aerospace Centre, were collected in this study. These acquisitions covered the period from 17 June 2014 to 27 November 2015. The parameters of these TerraSAR-X images are listed in Table 1. The Pixel Spacing of selected images along Range direction is 2.198 m, along azimuth direction is 1.965 m. The subset for Interferometric processing is about 18 km × 15 km (green rectangle in Figure 2a). It should also be noted from Figure 2 that the Lungui highway is surrounded with dense pounds and rivers. It locates close to three hydrological systems: Xi River, Rongui and Shunde Branch River. Therefore, Lungui highway is surrounded by plenty of pounds and lakes and large amount of silt. Most of these pounds are local cultivation. According to the geological material we collected, this area has extensive aquifers. The underground water level is influenced greatly by seasonal climate factors, with maximum levels varying by up to 1 m. Figure 2b is the radar amplitude image of the study area, with the highway region of interest outlined in the yellow rectangle. There are two main bridges along the highway, Rongguite Bridge and Anlite Bridge, which are analyzed in detail in Section 4. Figure 3b,d are the ground pictures of the two bridges, respectively.

SAR Acquisitions and Data Processing
In total 17 repeat-pass TerraSAR X-band Stripmap descending images, provided by the German Aerospace Centre, were collected in this study. These acquisitions covered the period from 17 June 2014 to 27 November 2015. The parameters of these TerraSAR-X images are listed in Table 1. The Pixel Spacing of selected images along Range direction is 2.198 m, along azimuth direction is 1.965 m. The subset for Interferometric processing is about 18 km × 15 km (green rectangle in Figure 2a). Two-pass DInSAR processing was adopted to generate the differential interferograms with the 17 TerraSAR images. The resolution along azimuth and range directions are 3.29 m and 2.64 m respectively, and the average incident angle is 26.4 • . SBAS-InSAR with four different deformation models has been tested to estimate the surface deformation of the study area. With the selection criteria of the perpendicular baseline smaller than 130 m and the temporal baseline shorter than 300 days, 57 small baseline interferograms were generated [18]. To remove the topographic phases, a 1-arc-second Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM:~30 m spacing) provided by NASA was utilized. In addition, to constrain the resolution for our narrow ribbon highway object, multi-looking ratio is set to be 1:1, to constrain the original resolution. The interferograms were filtered with the Goldstein filter to further reduce the phase noise. Then, the Minimum Cost Flow (MCF) was applied to unwrap the wrapped interferometric deformation phases.
The interferometric process is carried out by SARScape 5.2. Coherence candidates are selected based on a coherence threshold of 0.3. We selected the images carefully to ensure that most coherent points are distributed over the highway among the 57 total interferometric pairs, and deleted the ones with bad coherence and less points on the tested highway. Consequently only 33 high quality pairs, with densely distributed coherence pixels over the highway region, were selected. Unfortunately, temporal gaps occurred due to limited suitable interferometric pairs. Figure 5 shows the temporal and perpendicular baselines of those interferometric pairs. The number indicates the image number as shown in Table 2, and number 7 (acquired on 14 February 2015) is the most suitable master image among all pairs. Due to the low coherence, especially among images number 0-7, only 6 satisfactory pairs were selected, thus 2 temporal gaps could not be covered. Figure 6 shows the mean coherence map with clear bright coherent line along the highway. The following SBAS steps, including high coherence target selection, deformation modeling and deformation parameter estimation, were all carried out through Matlab after the unwrapped interferometric pairs were generated. Due to the large number of points distributed in interferometric pairs, the programming would be extremely time-consuming. Therefore, a subset around the highway area (yellow rectangle in Figure 2b) was extracted and tested. We consider both intensity, amplitude deviation and spatial coherence in the coherence target identification process, finally generating 201,703 candidates in total.

Deformation Parameter Estimation
The results are shown in Figures 7-9 (over coherent points only and with mean amplitude image as background, in slant-range projection). Figure 7 shows the mean linear velocities for MVM and PVM, respectively, where Figure 7a illustrates the mean value of all the temporal deformation velocities between every two time-adjacent SAR acquisitions, while Figure 8b is the only deformation parameter in PVM. From the color and scale comparisons it is obvious that MVM results shows

Deformation Parameter Estimation
The results are shown in   show the estimated coefficient results for SM and CPM respectively. As we can see, although these two models have the same number of unknown deformation parameters, the estimated results are extremely different. Evidently, dark color scale appears in CPM (coefficients range from −60 to 20 mm/year), indicating the parameters of all the points have higher magnitude than SM. Besides, total blue color region can be seen in Figure 9b, which implies uplift displacement coefficients c2 for CPM. Height correction for 4 models are illustrated in Figure 10. Apparently, height corrections obtained by SM, PVM, and CPM show similar distribution, whereas MVM performs with a darker color scale and larger magnitude.  Figures 8 and 9 show the estimated coefficient results for SM and CPM respectively. As we can see, although these two models have the same number of unknown deformation parameters, the estimated results are extremely different. Evidently, dark color scale appears in CPM (coefficients range from −60 to 20 mm/year), indicating the parameters of all the points have higher magnitude than SM. Besides, total blue color region can be seen in Figure 9b, which implies uplift displacement coefficients c2 for CPM. Height correction for 4 models are illustrated in Figure 10. Apparently, height corrections obtained by SM, PVM, and CPM show similar distribution, whereas MVM performs with a darker color scale and larger magnitude. Figures 8 and 9 show the estimated coefficient results for SM and CPM respectively. As we can see, although these two models have the same number of unknown deformation parameters, the estimated results are extremely different. Evidently, dark color scale appears in CPM (coefficients range from −60 to 20 mm/year), indicating the parameters of all the points have higher magnitude than SM. Besides, total blue color region can be seen in Figure 9b, which implies uplift displacement coefficients c2 for CPM. Height correction for 4 models are illustrated in Figure 10. Apparently, height corrections obtained by SM, PVM, and CPM show similar distribution, whereas MVM performs with a darker color scale and larger magnitude.

Accuracy Assessment of Different SBAS Models
Having obtained all the deformation coefficients for the four models over coherence points, we estimated the model phase in Equation (1)

Accuracy Assessment of Different SBAS Models
Having obtained all the deformation coefficients for the four models over coherence points, we estimated the model phase in Equation (1) To estimate the reliability of the different models, three different accuracy indexes were used. These were [35]: (1), which can be calculated by subtracting the model phase component from the original unwrapped phase observations. This index implies the fitting performance for each model [36]. The smaller the residual phase is, the higher the accuracy for the model selected. RP estimated results are shown in Figure 11. The color bar was set to (−0.5 to 0.5) radians. From the color characteristic in the Figure, we can see that the color for SM distributes around 0.1(mainly light green), for MVM distributes around −0.1 (most light blue), for CPM and PVM both higher than −0.5 (deep blue). Statistically the range of temporal average RP value for the four different models in radians, varied from−0.183 to 0.183, −0.920 to 0.664, −0.157 to 0.188, and −0.811 to 0.660 for LVP, PVM, SM and CPM, respectively. In our case MVM and SM models have better performance.
(2) Temporal Coherence (TC): according to literature [37], was calculated as the temporal mean TC value for different models of all coherent points. A higher TC value defines higher accuracy of the deformation model. From Figure 12, we can find obviously MVM shows the highest TC value (heavy red color all over the total image). PVM has the lightest TC color distribution. The average of TC value for each model is shown in the first row in Table 2, which quantitatively illustrates the performance for the 4 models, consistent with Figure 12. (3) HP Deformation: according to literature [29], uses the High-Pass (HP) deformation contribution.
In Equation (1), the part of 4π is called the Low-Pass (LP) deformation component, which is modelled by the function of deformation coefficients. As all unknowns in the SM model have been estimated, the LP deformation component can be obtained. The atmospheric phase component is highly correlated in space but poorly in time, while the High-Pass deformation signal is highly correlated not only in space but also in time, so the HP-deformation component can be estimated through temporal low-pass filtering (we use triangle filter method here) followed by spatial low-pass filtering (mean filter here) [37]. Figure 13 shows HP deformation between four models on each interferogram. We can see from it that all the HP deformation is lower than 10mm for each interferogram and each model, where the deformation of 10, 12, 13 and 16th interferogram are particular higher. It indicates these 4 interferograms are with larger residual phases. We calculated the Root Mean Square (RMS) value of temporal mean HP deformation for all the high coherent points, which are listed in Table 2. accuracy for the model selected. RP estimated results are shown in Figure 11. The color bar was set to (−0.5 to 0.5) radians. From the color characteristic in the Figure, we can see that the color for SM distributes around 0.1(mainly light green), for MVM distributes around −0.1 (most light blue), for CPM and PVM both higher than −0.5 (deep blue). Statistically the range of temporal average RP value for the four different models in radians, varied from−0.183 to 0.183, −0.920 to 0.664, −0.157 to 0.188, and −0.811 to 0.660 for LVP, PVM, SM and CPM, respectively. In our case MVM and SM models have better performance. (2) Temporal Coherence (TC): according to literature [37], was calculated as the temporal mean TC value for different models of all coherent points. A higher TC value defines higher accuracy of the deformation model. From Figure 12, we can find obviously MVM shows the highest TC value (heavy red color all over the total image). PVM has the lightest TC color distribution. The average of TC value for each model is shown in the first row in Table 2, which quantitatively illustrates the performance for the 4 models, consistent with Figure 12. (3) HP Deformation: according to literature [29], uses the High-Pass (HP) deformation contribution. In Equation (1), the part of ( , , ) − ( , , ) is called the Low-Pass (LP) deformation component, which is modelled by the function of deformation coefficients. As all unknowns in the SM model have been estimated, the LP deformation component can be obtained. The atmospheric phase component is highly correlated in space but poorly in time, while the High-Pass deformation signal is highly correlated not only in space but also in time, so the HPdeformation component can be estimated through temporal low-pass filtering (we use triangle filter method here) followed by spatial low-pass filtering(mean filter here) [37]. Figure 13 shows HP deformation between four models on each interferogram. We can see from it that all the HP deformation is lower than 10mm for each interferogram and each model, where the deformation of 10, 12, 13 and 16th interferogram are particular higher. It indicates these 4 interferograms are with larger residual phases. We calculated the Root Mean Square (RMS) value of temporal mean HP deformation for all the high coherent points, which are listed in Table 2. (3) HP Deformation: according to literature [29], uses the High-Pass (HP) deformation contribution. In Equation (1), the part of ( , , ) − ( , , ) is called the Low-Pass (LP) deformation component, which is modelled by the function of deformation coefficients. As all unknowns in the SM model have been estimated, the LP deformation component can be obtained. The atmospheric phase component is highly correlated in space but poorly in time, while the High-Pass deformation signal is highly correlated not only in space but also in time, so the HPdeformation component can be estimated through temporal low-pass filtering (we use triangle filter method here) followed by spatial low-pass filtering(mean filter here) [37]. Figure 13 shows HP deformation between four models on each interferogram. We can see from it that all the HP deformation is lower than 10mm for each interferogram and each model, where the deformation of 10, 12, 13 and 16th interferogram are particular higher. It indicates these 4 interferograms are with larger residual phases. We calculated the Root Mean Square (RMS) value of temporal mean HP deformation for all the high coherent points, which are listed in Table 2.

Discussion
In the scenario of using TerraSAR-X data, the four different deformation models illustrated inconsistent results, both in amplitude and geographic distributions. This demonstrates the importance of model selection. Through the accuracy analysis, based on three assessment indices, PVM had the least accuracy. The reason is that PVM does not incorporate the nonlinearity of the deformation in the process of modeling, which may cause greater modelling error if nonlinearity exists; Although CPM improves the accuracy by suppressing the fitting error, however, if the real deformation sequences do not follow the temporal 3rd order polynomial function, the estimated time series of deformation may be incorrect. The accuracy of real data deformation inversion will turn out to be low due to an unreasonable fitting functional model. MVM shows better fitting accuracy in our case, this is because of its sectional fitting of the total deformation processing. However, as introduced in Section 4.2, we generated only 33 high quality interferometric pairs, and the corresponding temporal gaps in MVM modeling refrains us to inverse the deformation all over the acquisition dates with MVM, this is a great deficiency for understanding of the dynamical deformation process. Besides, considering the geological characteristics of high-water content and high compression of soft clay, we assume the subsidence of the highway correlated to seasonal variations during the total time series subsiding process. As a result, we decide to carry out the time series deformation inversion based on SM, with comparison to other models, which will be discussed in the following section.

Overall Characteristics of the Surface Deformation Comparison
SBAS InSAR technique was attempted based on the aforementioned 4 models to derive and investigate the deformation process over the test highway on soft clay. As explained in Section 4.2, the HP-deformation component extracted from ∆ϕ res i (x, r) in Equation (1), can be added to LP-deformation. Consequently, the total deformation sequences can be generated. We also generate time series deformation with PVM, and CPM. Due to temporal gaps of MVM, we only obtained discontinuous deformation sequences. The deformation results between 2014/06/17-2015/11/27 are selected and compared, as shown in Figure 14. We can see similar deformation distribution characteristics in MVM, SM and PVM, including areas A (mainly subsidence) and B (slightly uplift). However, great uplift areas C and D are detected in MVM result (up to 32 mm), where are indicated as subsidence in SM and PVM results. Besides, the detected subsidence magnitude of MVM is larger than SM and PVM. The reason may be although MVM can perform high fitting accuracy to the different real deformation phenomena, when temporal gaps occurred, the generated deformation accuracy maybe decreased due to the increased singular problem in the parameter estimation [34]. From Figure 14d, we can find the deformation derived from CPM is quite different from the other three, including the opposite deformation distribution in A and the total large deformation magnitude. We can find large serious subsidence even up to 200 mm compared to the other three models. As discussed in Sections 3.2 and 3.5, CPM showed high sensitivity to noise and unsatisfied fitting accuracy. By considering the results presented using simulated data and the real TerraSAR-X time series deformation, we finally choose seasonal model in this case study. Figure 15 shows the time series of subsidence from 17 June 2014 to 27 November 2015. The overall accumulated deformations are generally following a subsiding trend (color distributed mainly ranging from orange to light blue). Spatially, more subsiding points exist in the top part to the Ronggui River. There are three main uplift zones in the Figure (area A, B, C), although the main deformation trend is subsiding. Darker red color appears where the ponds are distributed densely (where the small black rectangles are). Two bridges appear to be subsiding according to the darker color. Temporally, from 17 June 2014 to 01 January 2015, the subsidence is decreasing, with the uplifting points increasing spatially (warm color fades and light blue points dominate at 01 January 2015). As our measurements show, the uplift magnitude accumulated to a maximum of 14 mm at 01 January 2015, winter. In contrast, the deformation trend changes to subsidence from 14 February 2015, with the main color trend becoming warm. The subsidence reaches up to −76 mm by 31 August 2015, with the total color being dark red. From 22 September 2015, another period of uplifting started. It can be inferred that periodical deformation was the main deformation trend around the Lungui road over the 15 months-period. According to our qualitative analysis, both the magnitude and trend of subsidence can be related to season. The soft clay contains a large amount of water, so with temperatures rising gradually in summer, heat absorption contributes to the evaporation of a large percentage of water and results in surface subsidence. This is the key reason why the subsidence reaches the maximum in warm seasons. In contrast, in cold seasons, the temperatures drop gradually, and consolidation apparently surpasses the evaporation effect. This suppresses the ground subsidence significantly, even performing to be uplift in some areas. According to the records of Fuoshan weather bureau, from August 2015, the temperature performs to drop gradually, with an accumulated decrease of 16 • C until December. Another reason of this seasonal deformation fluctuation may be related to the increasingly serious rainfall during August to December. According to the rainfall records of Fuoshan, affected by the combined upper trough and bottom vortex, continually rainfall appeared since 27 August 2015. The average rainfall amount is 127.3 mm in Fuoshan city, with 100~250 mm at 65% automatic stations recorded, and 250 mm at Shunde automatic station. Due to increase of rainfall amount, both water content and discharge in water system increased correspondingly. With accelerated flow speed, perched ground water in the upper layer of subgrade increased due to the impact of rainfall and supply from surrounded water system. To further investigate the deformation processing characteristics over this area, we conduct detailed analysis over some ground features based on SM, which will be illustrated in the following section.

Deformation over Ground Features-investigation based on Seasonal Model
We conducted more detailed deformational investigations over two primary bridges (Rongguite Bridge and Anlite Bridge) and three feature points (see Figure 13, P1 to P3, one is located in Southern Medical University, the other two are surrounded by ponds near this highway).

Displacement Characteristics over Two Bridges
Ronguite Bridge, with a total length of 1.8 km, lying in the Magang river segment, is the most important junction for the Lungui Highway, connecting the Southern Medical University in Magang district and the Rongguidaoxiang industrial area. Total construction cost of this bridge was about 420 million yuan and it opened to traffic in 2015 (During June 2014 to December 2014, a small section of this road is still under surfacing, which is not included in the tested area in Figure 2b). The surface of the bridge is made of 6000 tons of pure iron (see Figure 3b,c). Due to the important location of this bridge, it's necessary to conduct long-term deformation monitoring. The permissible vertical postconstruction subsidence is 30 cm/yr for regular road segment, 20 cm/yr for culvert and 10 cm/yr for bridges' connections, according to the design criterion. The particular reinforcement treatment for the connection between bridge and subgrade is CFG (Cement Fly-Ash Gravel) pile and bag sandwell. For both the two bridges, the radius of CFG pile is 40 cm, the interval range is 1.3 m~1.6 m, with the material of mixed breakstone, stone chips and pulverized fuel ash. The bag sand-well is made of medium coarse sand, with muddy component lower than 3% percent and permeability coefficient higher than 5 × 10 cm/s. The length of CFG pile is designed according to in-situ calculation, should be deeper than the thickness and longer than 15 m in common.
We extracted all the high coherence points on Ronguite Bridge. After deleting the points that were wrongly detected distributed in water, 3012 points were extracted in total. Figure 16 shows the overall deformation sequences of Ronguite Bridge. The maximum subsidence over Ronguite Bridge is -54mm, which occurred on 31 August 2015. We plotted the detailed portrait deformation profiles to analyze the subsidence characteristics more clearly. The longitudinal direction profile along D1D2 is shown in Figure 16. After interpreting the sparse deformation over all the points, we extracted

Deformation over Ground Features-Investigation Based on Seasonal Model
We conducted more detailed deformational investigations over two primary bridges (Rongguite Bridge and Anlite Bridge) and three feature points (see Figure 13, P1 to P3, one is located in Southern Medical University, the other two are surrounded by ponds near this highway).

Displacement Characteristics over Two Bridges
Ronguite Bridge, with a total length of 1.8 km, lying in the Magang river segment, is the most important junction for the Lungui Highway, connecting the Southern Medical University in Magang district and the Rongguidaoxiang industrial area. Total construction cost of this bridge was about 420 million yuan and it opened to traffic in 2015 (During June 2014 to December 2014, a small section of this road is still under surfacing, which is not included in the tested area in Figure 2b). The surface of the bridge is made of 6000 tons of pure iron (see Figure 3b,c). Due to the important location of this bridge, it's necessary to conduct long-term deformation monitoring. The permissible vertical post-construction subsidence is 30 cm/yr for regular road segment, 20 cm/yr for culvert and 10 cm/yr for bridges' connections, according to the design criterion. The particular reinforcement treatment for the connection between bridge and subgrade is CFG (Cement Fly-Ash Gravel) pile and bag sand-well. For both the two bridges, the radius of CFG pile is 40 cm, the interval range is 1.3 m~1.6 m, with the material of mixed breakstone, stone chips and pulverized fuel ash. The bag sand-well is made of medium coarse sand, with muddy component lower than 3% percent and permeability coefficient higher than 5 × 10 −3 cm/s. The length of CFG pile is designed according to in-situ calculation, should be deeper than the thickness and longer than 15 m in common.
We extracted all the high coherence points on Ronguite Bridge. After deleting the points that were wrongly detected distributed in water, 3012 points were extracted in total. Figure 16 shows the overall deformation sequences of Ronguite Bridge. The maximum subsidence over Ronguite Bridge is −54 mm, which occurred on 31 August 2015. We plotted the detailed portrait deformation profiles to analyze the subsidence characteristics more clearly. The longitudinal direction profile along D1D2 is shown in Figure 16. After interpreting the sparse deformation over all the points, we extracted deformation during the period from 13 May 2015 to 27 August 2015. From Figure 17 we can see, both riversides (D1 and D2 in Figure 14) show stability. The reason for this is probably the iron reinforcing tower set up on the north riverside (D2 in Figure 16), the in-situ pictures of this tower are shown in Figure 3c. The iron tower, with its height of 124 m, width of 4 m and length of 12 m, bracing both sides of the 30 m-high bridge surface, can suppress the subsidence damage for the bridge. Correspondingly the deformation on area D2 shown in Figure 16 is close to 0, more stable than most points distributed along D1D2. Without an iron reinforcement instrument in another riverside area (D1), the deformation is slightly more serious, as shown in Figure 16. As we obtained, the maximum subsidence over this  Figure 17 we can see, both riversides (D1 and D2 in Figure 14) show stability. The reason for this is probably the iron reinforcing tower set up on the north riverside (D2 in Figure 16), the in-situ pictures of this tower are shown in Figure 3c. The iron tower, with its height of 124 m, width of 4 m and length of 12 m, bracing both sides of the 30m-high bridge surface, can suppress the subsidence damage for the bridge. Correspondingly the deformation on area D2 shown in Figure 16 is close to 0, more stable than most points distributed along D1D2. Without an iron reinforcement instrument in another riverside area (D1), the deformation is slightly more serious, as shown in Figure 16.     Figure 17 we can see, both riversides (D1 and D2 in Figure 14) show stability. The reason for this is probably the iron reinforcing tower set up on the north riverside (D2 in Figure 16), the in-situ pictures of this tower are shown in Figure 3c. The iron tower, with its height of 124 m, width of 4 m and length of 12 m, bracing both sides of the 30m-high bridge surface, can suppress the subsidence damage for the bridge. Correspondingly the deformation on area D2 shown in Figure 16 is close to 0, more stable than most points distributed along D1D2. Without an iron reinforcement instrument in another riverside area (D1), the deformation is slightly more serious, as shown in Figure 16.    Anlite Bridge, with a total length of 1.0km, lies across the Shunde Brunch River. Its upper surface is made of concrete. Figure 18 shows the overall deformation time series from 17 June 2014 to 27 November 2015. We detected a total of 1547 points all over the bridge, and the deformation temporally developed following the periodical changing, with a maximum subsidence of −39 mm, occurred on 31 August 2015, similar to Rongguite Bridge. However, in contrast to Rongguite Bridge, an obvious subsidence bow can be found in Figure 18, in the middle segment of the bridge, and temporally increasing. We plotted the subsidence profile along E1E2. From our measurements, the deformation accumulated up to −39 mm until 31 August 2015. Similar to Ronguite Bridge, both the riversides show stable performance (see the purple oval E1 and E2 area in Figure 19), due to the solid bridge pier on both sides holding up the whole bridge surface. As we investigated, there is no visible crack found in situ on the bridge in 2015 when it opened to traffic, however, unlike Ronguite Bridge, without iron strengthening, and with the obvious subsidence bowl in the middle bridge part, it is necessary to consider the risk of cracking into account in the long term. Anlite Bridge, with a total length of 1.0km, lies across the Shunde Brunch River. Its upper surface is made of concrete. Figure 18 shows the overall deformation time series from 17 June 2014 to 27 November 2015. We detected a total of 1547 points all over the bridge, and the deformation temporally developed following the periodical changing, with a maximum subsidence of −39 mm, occurred on 31 August 2015, similar to Rongguite Bridge. However, in contrast to Rongguite Bridge, an obvious subsidence bow can be found in Figure 18, in the middle segment of the bridge, and temporally increasing. We plotted the subsidence profile along E1E2. From our measurements, the deformation accumulated up to −39 mm until 31 August 2015. Similar to Ronguite Bridge, both the riversides show stable performance (see the purple oval E1 and E2 area in Figure 19), due to the solid bridge pier on both sides holding up the whole bridge surface. As we investigated, there is no visible crack found in situ on the bridge in 2015 when it opened to traffic, however, unlike Ronguite Bridge, without iron strengthening, and with the obvious subsidence bowl in the middle bridge part, it is necessary to consider the risk of cracking into account in the long term.

Temporal deformation over feature points
In order to show more dynamic subsidence evolution details for this soft clay area, we chose three feature points (shown in Figure 20a, red points P1, P2 and P3) and analyze their time series deformation characteristics. Figures 21-23 show the corresponding position for the three points Anlite Bridge, with a total length of 1.0km, lies across the Shunde Brunch River. Its upper surface is made of concrete. Figure 18 shows the overall deformation time series from 17 June 2014 to 27 November 2015. We detected a total of 1547 points all over the bridge, and the deformation temporally developed following the periodical changing, with a maximum subsidence of −39 mm, occurred on 31 August 2015, similar to Rongguite Bridge. However, in contrast to Rongguite Bridge, an obvious subsidence bow can be found in Figure 18, in the middle segment of the bridge, and temporally increasing. We plotted the subsidence profile along E1E2. From our measurements, the deformation accumulated up to −39 mm until 31 August 2015. Similar to Ronguite Bridge, both the riversides show stable performance (see the purple oval E1 and E2 area in Figure 19), due to the solid bridge pier on both sides holding up the whole bridge surface. As we investigated, there is no visible crack found in situ on the bridge in 2015 when it opened to traffic, however, unlike Ronguite Bridge, without iron strengthening, and with the obvious subsidence bowl in the middle bridge part, it is necessary to consider the risk of cracking into account in the long term.

Temporal deformation over feature points
In order to show more dynamic subsidence evolution details for this soft clay area, we chose three feature points (shown in Figure 20a, red points P1, P2 and P3) and analyze their time series deformation characteristics. Figures 21-23 show the corresponding position for the three points

Temporal Deformation over Feature Points
In order to show more dynamic subsidence evolution details for this soft clay area, we chose three feature points (shown in Figure 20a         P1 lies on the famous Lingyun Tower located in the South Medical University. Figure 21a illustrates the accurate location on Google Earth map, Figure 21b is the in-situ picture of P1, Figure 21c is the time series deformation of P1 over the period 17 June 2014 to 31 August 2015. We can see that the total deformation evolution follows a temporal seasonal pattern. The maximum amplitude of deformation variance is 24 mm. According to our measurements, the most serious subsidence on P1 occurred on 26 June 2015, with the magnitude of −41 mm, whereas minimal subsidence occurred on 01 January 2014 in winter, with a magnitude of −17 mm. As discussed above, P1 is located in the center of the half island, with soft clay as its the main soil in this area, seasonal deformation dominates the deformation performance. In summer, as temperatures increase, large amount of water contained in soft clay evaporates, thus the surface deformation is reflected as subsiding. While in winter, as temperatures decrease, water consolidation effects dominate, suppressing the subsidence process during this period. In contrast, the temporal evolution at P2 and P3 as shown in Figures 22  and 23, perform a totally subsiding trend with seasonal fluctuating. Both P2 and P3 are located on border sides between ponds close to the road. From Figures 22c and 23c, we can see the difference of those two points is that P2 shows a totally subsiding trend during the total period, whereas P3 performs more serious uplift during the period of 13 September 2014 to 14 February 2015 in winter. P2 and P3 perform more serious seasonal fluctuations than P1, with their deformation variances of 42 mm and 37 mm, respectively. The reason for these large variances may be their closer location to ponds and the surface material is silt, which contains more water and with even higher compression. Taking P3 for instance, during winter from 13 September 2014 to 14 February 2015, as the low freezing temperature increases the consolidation of the silt, the uplift displacement accumulates up to 10 mm. In summer (16 March 2015 to 31 August 2015), as climate becomes warmer, the evaporation process occurred instead, thus the surface deformation is reflected as subsiding. The maximum subsidence for P3 is −23 mm, which occurred on 13 May 2015.

Comparison with Leveling Measurements
In order to quantitatively evaluate the accuracy of the time series deformation obtained, we collected the leveling field measurements for direct comparison. The location of the benchmarks is shown in Figure 2b (see the red circle points K3, K4, K5), distributed close to Anlite Bridge. The date of leveling measurement is from June 2014 to February 2015, in order to carry out a precise comparison. We transferred the generated LOS deformation into vertical displacement according to S LOS = S v cos θ and extracted the 8 measurements dates data coincided temporally with our SAR acquisition dates. The total reference point of leveling and all the SBAS processing methods are the same pixel, which was selected according to in-situ investigation and registered deformation material collection. The comparison results are shown in Figure 24, blue solid squares define deformation results obtained by SM, while red solid circles illustrate leveling measurement. Figure 24a-c shows the comparison on each benchmark respectively. Obviously, SM shows better consistency with leveling results. We calculate the mean RMSE over all acquisition dates on 3 benchmarks with SM method, to be ±6.6 mm. P1 lies on the famous Lingyun Tower located in the South Medical University. Figure 21a illustrates the accurate location on Google Earth map, Figure 21b is the in-situ picture of P1, Figure  21c is the time series deformation of P1 over the period 17 June 2014 to 31 August 2015. We can see that the total deformation evolution follows a temporal seasonal pattern. The maximum amplitude of deformation variance is 24 mm. According to our measurements, the most serious subsidence on P1 occurred on 26 June 2015, with the magnitude of −41 mm, whereas minimal subsidence occurred on 01 January 2014 in winter, with a magnitude of −17 mm. As discussed above, P1 is located in the center of the half island, with soft clay as its the main soil in this area, seasonal deformation dominates the deformation performance. In summer, as temperatures increase, large amount of water contained in soft clay evaporates, thus the surface deformation is reflected as subsiding. While in winter, as temperatures decrease, water consolidation effects dominate, suppressing the subsidence process during this period. In contrast, the temporal evolution at P2 and P3 as shown in Figures 22 and 23, perform a totally subsiding trend with seasonal fluctuating. Both P2 and P3 are located on border sides between ponds close to the road. From Figures 22c and 23c, we can see the difference of those two points is that P2 shows a totally subsiding trend during the total period, whereas P3 performs more serious uplift during the period of 13 September 2014 to 14 February 2015 in winter. P2 and P3 perform more serious seasonal fluctuations than P1, with their deformation variances of 42 mm and 37 mm, respectively. The reason for these large variances may be their closer location to ponds and the surface material is silt, which contains more water and with even higher compression. Taking P3 for instance, during winter from 13 September 2014 to 14 February 2015, as the low freezing temperature increases the consolidation of the silt, the uplift displacement accumulates up to 10mm. In summer (16 March 2015 to 31 August 2015), as climate becomes warmer, the evaporation process occurred instead, thus the surface deformation is reflected as subsiding. The maximum subsidence for P3 is −23 mm, which occurred on 13 May 2015.

Comparison with Leveling Measurements
In order to quantitatively evaluate the accuracy of the time series deformation obtained, we collected the leveling field measurements for direct comparison. The location of the benchmarks is shown in Figure 2b (see the red circle points K3, K4, K5), distributed close to Anlite Bridge. The date of leveling measurement is from June 2014 to February 2015, in order to carry out a precise comparison. We transferred the generated LOS deformation into vertical displacement according to S LOS =S cosθ and extracted the 8 measurements dates data coincided temporally with our SAR acquisition dates. The total reference point of leveling and all the SBAS processing methods are the same pixel, which was selected according to in-situ investigation and registered deformation material collection. The comparison results are shown in Figure 24, blue solid squares define deformation results obtained by SM, while red solid circles illustrate leveling measurement. Figure 24a-c shows the comparison on each benchmark respectively. Obviously, SM shows better consistency with leveling results. We calculate the mean RMSE over all acquisition dates on 3 benchmarks with SM method, to be ± 6.6 mm.

Conclusions
A case study was presented assessing long-term deformation for a stretch of highway, built on soft clay subgrade, after road embankment settlement construction over a 15 months-period from June 2014 to November 2015, based on comparative analyses of four widely used InSAR time series deformation models. The deformation coefficients for the different models were computed. Three accuracy indices were used to assess the residual phase. The comparison showed differences in the measured deformation in magnitude and distribution, implying different models can significantly impact the deformation inversion result over the study area. Finally, SM was recommended here for modelling long-term dynamic subsidence for this highway region over soft clay subgrade. After our comparative model analysis, we carried out seasonal time series deformation investigations based on SM over the test highway area. According to our measurements, the most obvious subsidence occurred in August, up to −76 mm, whereas the uplift effect in January, with magnitude accumulated to 14 mm. The reasons suggested for this seasonal effect are that during winter as the low temperature influence the consolidation of the mucky clay, and the increased rainfall into soil expanded the volume of soft layer, thus the subsidence was suppressed; while in summer as climate becomes warmer, evaporation effect occurred instead which make the void water discharged, thus the surface deformation appears to be subsiding. We also conducted concrete deformation monitoring and analysis over ground features, including the important Rongguite Bridge and Anlite Bridge, and found that due to an iron strengthening tower installed close to one riverside of Rongguite Bridge, it performs to be stable; unlike Ronguite Bridge, without iron strengthening, we found a subsidence bowl in the middle bridge part on Anlite bridge, implying an urgent necessity to consider the risk of cracking in the long term. Compared to field leveling deformation measurements, we found seasonal model showed good consistency, and the final accuracy of derived deformation is estimated to be ±7 mm.
Future work will concentrate on geological processes in soft clay which may influence the temporal surface displacement. A new physical temporal model considering geological characteristics of soft clay will also be explored.
Author Contributions: X.X. designed the experiments and produced the results; X.X. and H.-C.C. analyzed the precipitation data; X.X., J.Z. and Z.S. helped to collect and analyze the leveling measurement in the real data experiment; Z.Y. and L.C. contributed to the discussion of the results. X.X. drafted the manuscript. All authors contributed to the study, reviewed and approved the manuscript.