INVESTIGATION ON INSAR TIME SERIES DEFORMATION MODEL CONSIDERING RHEOLOGICAL PARAMETERS FOR SOFT CLAY SUBGRADE MONITORING

The stability control is one of the major technical difficulties in the field of highway subgrade construction engineering. Building deformation model is a crucial step for InSAR time series deformation monitoring. Most of the InSAR deformation models for deformation monitoring are pure empirical mathematical models, without considering the physical mechanism of the monitored object. In this study, we take rheology into consideration, inducing rheological parameters into traditional InSAR deformation models. To assess the feasibility and accuracy for our new model, both simulation and real deformation data over Lungui highway (a typical highway built on soft clay subgrade in Guangdong province, China) are investigated with TerraSAR-X satellite imagery. In order to solve the unknows of the non-linear rheological model, three algorithms: Gauss-Newton(GN), Levenberg-Marquarat(LM), and Genetic Algorithm(GA), are utilized and compared to estimate the unknown parameters. Considering both the calculation efficiency and accuracy, GA is chosen as the final choice for the new model in our case study. Preliminary real data experiment is conducted with use of 17 TerraSAR-X Stripmap images (with a 3-m resolution). With the new deformation model and GA aforementioned, the unknown rheological parameters over all the high coherence points are obtained and the LOS deformation (the low-pass component) sequences are generated.


INTRODUCTION
Roads and highways built on soft clay subgrade are more prone to subsidence and induced instability.The stability control is one of the major technical difficulties in the field of highway subgrade construction engineering (Xue, 2011).Therefore, monitoring long-term surface deformation of highways on soft clay subgrade is crucial for understanding the dynamics of the settlement process and preventing potential hazards (Ge, 2008).In the flow of time series InSAR algorithm, building deformation model is a crucial step, which is to model the temporal functional relationship between displacement component and deformational parameters over all the high coherence targets.Accurate and reliable deformation model can not only increase the accuracy of deformation estimation, but also control the residual phase within the range of [-π，π], consequently influence the process of deformation parameters estimation, high coherence points optimization and phase unwrapping (Zhang, 2012).It can also help to the interpretation of deformation results.However, according to the literatures published, we found most of the InSAR deformation models for deformation monitoring are pure empirical mathematical models (Yu, 2013;Dong, 2014;Kim, 2010;Hetland, 2012), without considering the physical mechanism of the monitored object.For highways built on soft clay areas particularly, the deformation characteristics are complicated and nonlinear, if we choose the single pure empirical mathematical function to model the complicated physical deformational process, the real dynamic discipline of deformation cannot be truly reflected definitely.Consequently, the accuracy of the deformation monitoring and displacement * Corresponding author prediction maybe decreased, and the corresponding long-term analysis and interpretation for displacement after the highway construction maybe misled.
According to the research in the field of Soil and Rock Mechanics, soft clay deformation is especially correlated to soil rheology.Rheology is a basic and common soil phenomenon related to time, and a temporal effect and important engineering characteristics for soft clay.Rheological parameters, including viscosity and elasticity modulus, are primary parameters indicating rheology for soft clay.Rheology exists in every kind of soft clay (Zhao 1996).If we take rheology into consideration, inducing rheological parameters into traditional InSAR deformation models, not only the accuracy of long-term deformation monitoring can be greatly increased, but also the interpretation for the estimated deformation be more reasonable.In addition, the rheological parameters estimated through InSAR technique can also be provided as a reference for the researchers in the field of Soil and Rock Mechanics, thus broadening the application of InSAR technique.
As discussed above, in order to improve the accuracy of longterm deformation monitoring after highway embankment settlement construction, and interpret the displacement result more reasonably, we establish a new InSAR time series deformation model.The rheological characteristics for soft clay are considered in traditional InSAR deformation model.
The model is established based on the mathematical relationship between strain and rheological parameters (viscosity and elasticity modulus) according to Kelvin Model, a commonly used rheological model in the field of Rheology (Liu 2006).After investigation on the relationship between SAR LOS (line of sight) deformation and strain in Kelvin model, we obtain the mathematical function for LOS surface deformation with parameters of viscosity, elasticity modulus and time, hence the time series set of SAR phase functions can be established.

InSAR time series model considering rheological parameters
Kelvin model in theories of Rheology can be written as (Yu, 2003): where ε defines the strain for material，  is a constant, which defines the invariable outside load.It can be obtained by the highway structure characteristics investigation and soil gravity experiment of the upper soil over the soft clay.E defines elasticity modulus for material，η defines viscosity，E, η are important parameters in theories of Rheology.They are the unknowns of our model, which will be illustrated thereafter.The functional relationship between the vertical subsidence S  and strain S v can be written as: where H is the mean thickness of highway, which can be obtained from the highway designed materials.t 1 , t 2 are the beginning and ending time for deformation, respectively.With (1) and (2), S  can be written as: The unwrapped interferometric phase can be written as (Zhang, 2012): where ,  define index of coherent point and interferometric pair, respectively.S LOS is the deformation along the length of sight (LOS).λ is the wavelength, B m is perpendicular baseline of the interferometric pair,  is the incident angle, and  is the distance between the radar sensor and coherent point, ∆H is the height correction, another unknown in our model.∆φ i res is the residual phase, including phase noise, atmospheric delay, and high-pass (HP) deformation components.When the horizontal component is ignored, S LOS can be written as: According to (3) -( 5), the functional relationship between unwrapped interferometric phases φ i m and rheological parameters E, η can be obtained, which can be written as: where E, η, ∆H are the unknowns.

Rheological parameters estimation
According to (6), if we have  interferometric pairs, then  functions with three unknowns can be generated.It becomes a system of nonlinear equations with three unknown parameters.
The nonlinear least squares adjustment theories can be utilized here to estimate the unknown rheological parameters.We choose three algorithms: Gauss-Newton(GN) (Panl, 2016), Levenberg-Marquarat (LM) (Yang, 2010), and Genetic Algorithm(GA) (Velez, 2011), to determine the unknown parameters.GN method is based on Newton iteration algorithm, which is using the linear function to model the nonlinear characteristics of equations, and the traditional adjustment formulas can still be used, which keeps convenience in the calculation.LM algorithm is another improved algorithm based on Newton iteration.It can adapt to the solvation of singular or ill-condition problem of Jacobi matrix.
GA is an overall optimization algorithm, which does not rely on the concrete functional relationship between observations and unknowns.If the functional relationship is too complicated, GA will show considerable advantage.In GA, the filter function should be established following the minimum residue phase principle f =||∆φ res || , and then the population size, initial population range and the condition of iteration termination should be set up.After the variation and intersection, the final individuals of population which generate the minimum residue phase will be exported as the final solutions for (6).The Simplex Algorithm is used to optimize the GA estimated solutions (Yang, 2017).We can utilize the solutions as the initial inputs for Simplex Algorithm.After the iteration, the final searching results can be used as the optimized solutions.

Algorithm flow for soft clay subgrade subsidence monitoring
As discussed above, the InSAR time series model considering rheological parameters can be applied in the subsidence monitoring for soft clay subgrade, and the estimations of viscosity and elasticity modulus can be determined through our algorithm.The integrated algorithms flow is shown in Figure 1.

EXPERIMENT
To assess the feasibility of our proposed InSAR time series deformation model and the parameter estimation algorithms, both simulated and real data experiment are designed and implemented.

Simulated Experiment
We need to simulate each phase component on the right side of (6).With the known SAR sensor parameters (TerraSAR Stripmap sensor is used here), including spatial and temporal baselines of each interferometric pair, the topographic phase component ∆ over all the pixels.Values of all the rheological parameters over 200 high coherent pixels can be picked up from the simulated field, which are used as true values in the following validation.Consequently, the corresponding phase value on the left hand side of ( 6) can be obtained as observations in the following parameter estimation process.The noise added here is 0.1rad with Gauss function (into the ∆   component).We carried out three groups of simulation experiments to test the three estimated algorithms (GN, LM and GN) aforementioned respectively.After the estimation, newly obtained parameter solutions can be compared to the original simulated real parameter value, thus the accuracy of each model can be evaluated.Table 1.Quantitative comparisons of accuracy for three algorithms in simulation In the processing of calculation, the calculation time of GN is the best, only 144.394 seconds to export the final results.However, results of "NUN" appeared on 3 pixels (with indexes of 76, 101 and 156 respectively).GA shows the second promising calculation time, with 662.816 seconds to export the final results, whereas LM methods with a total of 1223.742seconds.
After eliminating the NUN value from GN results, we compare all the estimations to true values, which is shown in Figure 2.
From Figure 2 we can see estimations over most pixels are close to the true values, indicating the three algorithms are feasible in this case.However, except the three eliminated pixels, there still exist incorrect estimation (with index of 191) in GN estimations.
Table 1 is the quantitative root mean square error comparison of accuracy for three algorithms, RMSE of GN method is not reliable due to the incorrect pixel estimation.The reason for this phenomenon is GN algorithm show great incompatibility for illconditioned equation, which means if the equation is not illconditioned, the solutions are reliable, whereas with characteristics of ill-condition or singularity, the results are incorrect.LM method can compensate for this deficiency, and the results in Figure 2 and Table 1 are satisfied.However, we tried to change the initial values of unknowns to test the calculation time for LM method, and we find that if the initial values are closer to the true ones, the calculation time can be reduced to 225.68s, whereas the initial value are significantly far away from the true ones, the iteration processing cannot even converge.This indicates that LM method shows higher reliance to initial values of unknows.

Real-data experiment
We conduct a real data case study to test the feasibility and reliability of the InSAR deformation time series model considering rheological parameters.As discussed in 3.1, we decide to use LM and GA methods to estimate the rheological parameters respectively.

Study area:
The selected structure built on soft clay subgrade in our study is called Lungui highway, in Shunde District, Fuoshan City of Guangdong Province, in southern China (see Figure 3).Lungui highway, with a total length of 11.95 km, was designed to be the main connection from Shunde west city district of Fuoshan, extending north to G321 highway, and south to Hejiu West highway.It connects the Longzhou, Nanguo and Hejiu Highways, becoming one of the three trunk South-North lines in Fuoshan.As Figure 3a shows, the red rectangle outlines the spatial coverage of the selected TerraSAR images, while the green solid rectangle is the sub-set for generating the inteferometric results.It should also be noted that the Lungui highway is located close to three hydrological systems: Xi River, Rongui and Shunde Branch River.Plenty of pounds distributed around Lungui highway as well, which contain large amount of silt.According to the geological material we collected, this area has extensive aquifers.Figure 3b is the amplitude image of the study area, the yellow rectangle demonstrates the highway region of interest.
There are two main bridges along the highway, Rongguite Bridge in the middle part and Anlite Bridge in the north-east corner.
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 June 17, 2014 to November 27, 2015.The Pixel Spacing of selected images along Range direction is 2.198m, along azimuth direction is 1.965m.The subset for Interferometric processing is about 18 km × 15 km (green rectangle in Figure 3a), which covers all the Lungui highway.

Data processing:
Two-pass DInSAR processing was adopted to generate the differential interferograms with the 17 TerraSAR images.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.Figure 4 shows the temporal and perpendicular baselines of those interferometric pairs.The number index defines the index of each obtained image, and number 7 (acquired on February 14, 2015) is the most suitable master image among all 57 pairs.To remove the topographic phases, a 3-arc-second Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM: ~90 m spacing) provided by NASA was utilized.
In addition, multi-looking operation with 1 look in the range and 5 in the azimuth direction is conducted.The interferograms were filtered with the Goldstein filter to further suppress the phase noise.Then, the Minimum Cost Flow (MCF) was applied to unwrap the interferometric phases.SARScape 5.3 and Envi 5.3 were used to implement all the interferometric processes.The coherence threshold of 0.3 is adopted through the processing in SARScape software to check the main coherence distribution.
We selected the images carefully to ensure that most coherent points are distributed over the highway among the 57 interferometric pairs.The images with bad coherence and less points on the test highway are deleted.Consequently 23 high quality pairs, with densely distributed coherent pixels over the highway region, were generated.After the unwrapped interferometric pairs were generated, we used Matlab to carry out our integrated algorithm mentioned in 2.3, including high coherence target selection, deformation modelling and rheological parameter estimation.Due to the large number of points distributed in interferometric pairs, the programming would be extremely time-consuming.We masked highway area (yellow rectangle in Figure 3b), considering both intensity, amplitude deviation and spatial coherence index in the coherence target identification process, finally generating a total of 21099 candidates.

CONCLUSIONS
In this study, a new time series InSAR deformation model considering rheological parameters (including elasticity modulus and viscosity) is proposed.To assess the feasibility and accuracy for our new model, both simulation and real deformation data over Lungui highway (a typical highway built on soft clay subgrade in Guangdong province, China) have been investigated with TerraSAR-X satellite imagery.Three algorithms, Gauss-Newton(GN), Levenberg-Marquarat(LM), and Genetic Algorithm (GA), have been utilized and compared to estimate the rheological parameters.According to the results of simulated experiment, we find GN shows the best calculation efficiency whereas LM the worst performance.However, due to its great incompatibility for ill-conditioned equation, RMSE of GN method is not reliable.Although LM method can compensate for this deficiency, due to its high reliance to the initial values of unknowns, the iteration processing for sampled pixels may probably not even converge.We also conducted preliminary real data experiment with use of 17 TerraSAR-X Stripmap images (with a 3-m resolution) covering part of Lungui highway, with a total length of 11.95 km in Shunde District, Fuoshan City of Guangdong Province, in southern China, as the test highway.The interferometric processing is carried out using SARScape 5.2.We finally picked up 23 unwrapped high quality interferograms, generating 21099 high coherence points distributed along the highway.With the new deformation model and three algorithms aforementioned, we acquired the unknown rheological parameters over all the high coherence points.The corresponding Low-pass time series deformation is obtained consequently.Future work will be concentrated on the high-pass deformation component and the total subsidence sequences generation.In order to verify the accuracy of our new model, we plan to use the RMSE of residue phase over all the targets as an evaluation index for the three algorithms respectively.And soil characteristic material and the field levelling results collected and will also be used to validate the rheological parameters and subsidence results.

Figure 1 .
Figure 1.Algorithm flow can be determined.200 random locations of high coherence points and 10 interferometric pairs are simulated in our simulation.Surface function is used to simulate the field of rheological parameters, and Gauss function is used to simulate

Figure
Figure 2. Simulated parameters comparison

Figure 3 .
Figure 3. Study area (a) Optical image; (b) amplitude image of the area with the highway region of interest outlined in the yellow rectangle and (c) location in China

Figure 5 .
Figure 5.Time series low-pass deformation (LOS) over test highway area3.2.3 Obtained rheological parameters and time series LP deformation:In our case study, we planned to use LM and GA methods to estimate the rheological parameters respectively.When using LM algorithm, although we changed couple of initial values for the unknown parameters, it could not show convergence on most coherence points.The failure of obtaining the estimations with LM indicates LM method may not feasible in our case study.More initial values and programming test still need to be carried out in our future work.Great reliance to initial values of unknowns may be the key limitation of LM algorithm in case of unavailability of priori deformation information.With use of GA algorithm ( with SA optimization), we implemented the estimation process, and the preliminary estimated parameters are shown in Figure4.The total calculation time is 120866 seconds.From Figure6we can see, the total range of elasticity modulus is around (0, 20000), actually the about 75% pixels are with values around(10000, 20000).Total range of viscosity is around (2000, 5000)pa*s, whereas height correction is around (-10,10) m.We use these preliminary parameters to obtain the time series of Low-Pass deformation along the direction of LOS, from 22 August 2014 to 27 November 2015, shown in Figure5.The total deformation range is around[-40 20]  mm, the 16 deformation images are all with reference to the time of June, 17 th 2014.The overall subsidence trend is nonlinear, from August, 22 nd 2014 to May, 13st 2015 subsiding is slightly increasing, from June, 26 st 2015, the subsiding is becoming serious dramatically.