#### Abstract

In cold regions, the long-term stability of engineering facilities is unavoidably influenced by the negative temperature, freeze-thaw process, dry-wet process, and dynamic loading conditions induced by earthquakes and traffic loads. In order to investigate the effects of different cyclic stress paths on the evolution of dynamic mechanical properties of frozen silt clay, a series of cyclic triaxial tests with variation confining pressure (VCP) or constant confining pressure (CCP) were performed. Triaxial low-temperature apparatus (MTS-810) was taken advantage of to simulating various cyclic stress paths by changing cyclic loading conditions of axial stress and confining pressure. In this paper, the evolution features of the axial resilient modulus, damping ratios, and the shape of hysteresis loops with an increase in the number of load cycles under different dynamic stress paths are comprehensively studied. The results show that the loading angle of cyclic stress path and the phase difference between cyclic axial stress and confining pressure are the main factors that remarkably affect the development characteristics of the resilient modulus and damping ratio. With increasing of the loading angle and phase difference, the resilient modulus increases, but damping ratio increases with increasing of loading angle and with decreasing of phase difference. With the continuous increase in the number of loading cycles, the samples of frozen soil show compacting and hardening characteristics. With an increase in the number of load cycles, the shape of hysteresis loop becoming narrows, the resilient modulus decreases at the initial stage and then gradually increases, and the damping ratio stably decreases. According to contrastively analyzing the evolution of dynamic parameters and the shape features of hysteresis loops under various cyclic stress paths, it can be clearly discovered that the evolution of sample microstructure and the development of dynamic characteristics of frozen samples have obvious dependence on the cyclic stress path. Therefore, the effects of variable confining pressure (VCP) on dynamic behaviors of frozen soils are nonnegligible in practical cold region engineering.

#### 1. Introduction

With the continuous expansion of human living space, more and more engineering facilities (such as railways, highways, cannels, and pipelines) have been built in cold regions [1–4]. Meanwhile, the serviceability of these infrastructures is inevitably affected by the environmental change and geological hazard such as global-warming-related deformation, frozen-thaw, dry-wet process, and dynamic loading (e.g., earthquakes, traffic loads, or ocean wave storms) [5–9]. Especially dynamic load that is induced by earthquake, such as the 8.1-magnitude Kunlun Mountains earthquake in November 2001, caused serious cracks on some sections of the Qinghai-Tibet highway and the Qinghai-Tibet railway [10], seriously threatening the safety of the infrastructures. Long-term dynamic load may lead to vertical deformation exceeding the acceptable limit, and the infrastructures may lose serviceability [11]. Therefore, many investigations on the dynamic properties of frozen soil have been conducted [11–15]. But, those works mainly focus on dynamic properties of the frozen soil under constant confining pressure (CCP) condition. For seismic load, the S-wave and P-wave acting on the subgrade will be in random combination with different amplitudes, frequencies, and loading orders, while the CCP condition can only simulate a specific combination of normal stress and shear stress. Also, for the case of vehicle load acting on the surface structure, the lateral stress will also change passively with the change of the axial load, instead of remaining constant due to the Poisson effect of geotechnical materials. Therefore, it is very limited to simulate complex in situ dynamic stress field by axial vibration loading under CCP condition. Cai et al. [16] compared the one-way cyclic triaxial behaviors of saturated clay under several stress paths with VCP and CCP. The study found that VCP tests were more appropriate than CCP tests in the simulation of traffic loading conditions because the VCP stress paths can provide the coupling effects of the variable deviatoric stress and variable confining pressure. The study by Gu et al. [17] indicated that vibration sources applied a repeated loading on the subsurface soil layers may lead to complex stress field, in which all the stress states of the soil elements may generally vary in long-term workings. However, the case of VCP is more suitable for simulating the actual stress field [18, 19]. But, similar research methods are never used in the laboratory test for frozen soils. Due to the existence of the ice, the mechanical properties of frozen soil are more complex than those of unfrozen soils [20–26]. In order to further develop the understanding for the dynamic mechanical properties of frozen silt clay under complex stress field, the triaxial low-temperature apparatus (MTS-810) was used to simulate various cyclic stress paths. To some extent, the evolution law of resilience modulus with the number of load cycles increase under a special dynamic stress path can partly reflect the development characteristics of mechanical properties of the frozen soil samples. Therefore, the law can further reflect the stress path dependence of mechanical property evolution under dynamic stress path. The damping ratio is an important index to represent the soil’s ability of absorbing vibration energy. The hysteresis curve reflects the characteristics of elastic-plastic deformation, stiffness, and energy dissipation of soils under dynamic load. The hysteresis curve is the basis of determining and analyzing nonlinear dynamic response [27]. Meanwhile, the microscopic damage characteristics of rock and soil materials can be analyzed by the development of the macroscopic shape of the hysteretic curve. Therefore, this paper emphatically studied the evolution law of the dynamic parameters and the shape of hysteresis loop with increasing number of load cycles under different cyclic stress paths. According to the work in [11, 28], the dynamic modulus and damping ratio of the frozen soil are influenced by the many factors, such as the dynamic stress (strain) amplitude, temperature, and moisture content. Ling et al. [29] studied the relationship between the dynamic elastic modulus ratio and normalized shear strain under varies different conditions through loading step by step. But, for multistage loading, the mechanical properties of the frozen soil are inevitably influenced by the previous stress history. Zhang et al. [30] carried out a series research about dynamic resilient modulus and damping ratio with axial strain development, and the evolution law of mechanical properties of frozen soil samples during increasing of axial strain was revealed. But, few investigation works were focused on the evolution of dynamic parameters and the hysteresis shape of frozen soils with increasing number of loading cyclic under different dynamic stress paths.

In this paper, dynamic tests with three types of cycle stress paths were simulated to present the evolution characteristics of resilient modulus, damping ratios, and shape of hysteresis loops of frozen samples with increasing loading cyclic under different cycle stress paths. Moreover, the main influencing factors of mechanical properties for frozen soil are deeply analyzed.

#### 2. Experimental Program

##### 2.1. Test Apparatus, Test Materials, and Sample Preparation

Dynamic tests with variable confining pressure (VCP) or constant confining pressure (CCP) were conducted by the low-temperature triaxial apparatus (MTS-810) at the temperature of −6°C. The apparatus can control the amplitudes, frequency, and waveform of the axial stress and confining pressure, and the initial phase difference between the axial loading and confining pressure can be set up independently. The hydraulic oil was used in the confining pressure system. The outer wall of the chamber is wrapped around pipes that allow the cooling liquid to flow. So, the hydraulic oil not only provides pressure medium but also ensures if the required negative temperature environment is stable. The main parameters about the test systems are described as follows: the maximum axial load and confining pressure are 100 kN and 20 MPa, respectively, the range of the axial frequency is 0–50 Hz, the maximum displacement is 85 mm, and the range of temperature is room temperature of −30°C. The experimental data are automatically collected by the computer system.

In this study, the test material is silty clay collected from the site along the Qinghai-Tibetan Railway, China. The basic physical parameters are given in Table 1. The grain size distribution of the soil is shown in Figure 1.

Before the test, soils are placed under the sun for 3-4 days and then pulverized and sieved through a 2.0 mm screen to generate the prepared sample. Afterwards, appropriate amount of soil and deionized water were utilized to prepare the spare soil with a water content of 20.5%, which was placed in a sealed bag for more than 24 h. After that, appropriate amount of spare soil was taken to prepare samples with a height of 125.0 mm (the error is ±0.5 mm) and a diameter of 62.0 mm. The quality of the spare soil determines the dry density level of the sample (1.64 g/cm^{3}). Then, the specimens were put into a copper mold, with permeable stone on top and bottom, and fixed with an iron frame. Then, the samples with the mold were placed in a jar with an air extractor (shown in Figure 2(a)) and were held for 4 h under pressure 600 kPa. Then, the deionized water was injected into the jar through the pressure, and the samples were soaked for at least 24 h. Subsequently, the samples were taken out from the jar and were placed in a refrigerator with a temperature of −30°C for more than 24 h. Then, all the samples were carefully trimmed with a height of 125.0 mm and diameter of 61.8 mm on a lathe in a cold room under the temperature of −6°C, and then, each sample was sealed with a piece of rubber membrane as well as two resinous end caps as shown in Figure 2(b). Prior to the start of the experiment, all samples were placed in a thermostat at −6°C. The sample contrast before and after the experiment is shown in Figure 2(c).

**(a)**

**(b)**

**(c)**

##### 2.2. The Dynamic Stress Path

In this paper, three types of cyclic stress path are adopted in this investigation: (1) three directions of cycle stress path are selected on the stress plane of normal stress and shear stress ( plane) to perform three groups of experiments under approximate equal stress path length, as shown in Figure 3(a); (2) as shown in Figure 3(b), the nonlinear cyclic stress path on the plane is obtained by changing the phase difference between principal stress *σ*_{1d} and *σ*_{3d}; and (3) cyclic stress paths with roughly equal deviatoric stress amplitude are simulated on the plane of mean principal stress and deviatoric stress ( plane) to observe the influence of mean principal stress on the dynamic characteristics of frozen soil samples, described in Figure 3(c). To distinctly express the characteristics of different cyclic stress paths, the parameters’ amplitude ratio and the half-length value of cyclic stress path A and loading angle for cyclic stress paths *α* on the stress plane of and amplitude ratio *ζ*^{ampl} and half-length of the cyclic stress path B on the stress plane of , as well as the principal stress *σ*_{1d} and *σ*_{3d} vibration equations, were introduced in [31]. Also, the loading angle for the cyclic stress path on the stress plane of is expressed as follows:

**(a)**

**(b)**

**(c)**

The condition of all the test samples are summarized in Tables 2 and 3.

#### 3. Test Results and Discussion

##### 3.1. Dynamic Resilient Modulus

In order to accurately describe the evolution characteristics of resilient modulus of frozen silty clay with an increase in the number of load cycles under different cyclic stress paths, the calculation method is adopted as in Figure 4. When the deviator stress is partly unloading Δ*q*, it will produce the corresponding axial strain rebound Δ*ε*_{r}. Therefore, we define the axial resilient modulus as follows:

**(a)**

**(b)**

**(c)**

Based on the measured values, the relationship between the standard logarithm of cycle number and resilient modulus under all stress paths is a fitted quadratic. The fitting parameters *a*, *b*, and *c* and *R*^{2} are shown in Table 4, and the expression is given by the following equation:

The relationship between the number of load cycles and resilient modulus with different half-length of cycle stress paths value A in the same direction and different loading directions with approximately equal half-length are shown in Figure 5. The solid lines represent the optimal fitting of the experimental data. It could be seen that the resilient modulus decreases rapidly with the increase of cyclic loading before the initial vibration stage. After that, the resilient modulus retains a slightly upward trend. With an increase in the half-length of the cyclic stress path on the same cyclic loading direction, the initial stage resilient modulus decreases. Also, the difference of resilient modulus is obvious in this stage. Subsequently, it is clear that, with an increase in the number of cyclic loading, the resilient modulus level at the same loading direction is getting closer and closer. This phenomenon precisely illustrates that, with the continuous dynamic loading, the samples are strengthening. For the same cycle stress loading direction, although the stress path length is different, the samples were induced the same anisotropic behavior structure to adapt the stress field. Obviously, with an increase in the loading angle *α* of cycle stress path, the resilient modulus increases. Compared to the half-length of the cyclic stress path, the loading angle of the cyclic stress path has a more obvious influence on the resilient modulus, which indicates that the loading angle of the cyclic stress path is the main factor affecting the resilience characteristics of frozen soil.

As shown in Figure 6, the development tendency of is obtained under the different phase difference *φ* of principal stress *σ*_{1d} and *σ*_{3d}. The whole variation trends are identical with the former cases; decreases in the initial some hundred cycles and then enters a stable increase stage. The comparison results illustrate that increases as the phase difference increases under the same number of load cycles. According to analyzing those nonlinear stress paths in a normal-shear stress plane, the path loci for *φ* = 45° and *φ* = 270°or *φ* = 90° and *φ* = 225° are similar but rotate in opposite directions. Obviously, when the dynamic stress path rotates in an anticlockwise direction, is larger than that in the clockwise condition. For nonlinear cyclic stress paths, the loading orders of normal and shear stress are different from those for the linear condition of synchronous and asynchronous loading. As is shown in Figure 3(b), for the clockwise stress path rotation, when the shear stress is loading from the minimum value to the maximum value, the normal stress is unloading from the middle value to the minimum and then loading to the middle value. When the shear stress is unloading from the maximum value to the minimum value, the normal stress is loading from the middle value to the maximum value and then unloading to the middle value. For the anticlockwise condition, the loading order of normal stress is completely opposite to that in the former case in the process of shear stress loading or unloading. Therefore, in the process of shear stress unloading, the normal stress level during cyclic stress path clockwise rotation is higher than that during counterclockwise rotation. Also, the normal stress first increased and then decreases when the clockwise condition, but first decreases and then increases for the case of counterclockwise condition in the shear stress unloading process. This means that, in the process of shear stress unloading, along the clockwise stress path, the produced rebound elastic strain is larger than that in the anticlockwise condition.

In order to understand how the amplitude of the mean principal stress *p*^{ampl} affected the dynamic proprieties of frozen silty clay, a new series of tests are conducted under the approximately equal amplitude of deviatoric stress *q*^{ampl}. The comparison results of permanent strain accumulation with loading cyclic number are shown in Figure 7 for further confirming the effects of amplitude ratio *ζ*^{ampl} on the accumulation deformation behaviors. The absolute value of *ζ*^{ampl} decreases, namely, the length of dynamic stress path increases and the axial strain is faster accumulating under approximately equal amplitude of deviatoric stress. Through the result, it is discovered that the increase in amplitude of mean principal stress can obviously accelerate the permanent axial strain accumulation, but the loading order of deviatoric stress and mean principal stress is not obvious for the strain accumulation of this frozen soil.

The relationship between resilient modulus and is plotted in Figure 8. Obviously, as the cyclic stress loading angle increases, the resilient modulus increases. It is remarkable fact that the elastic deformation characteristic has an obvious dependence on the cycle stress path, and the loading sequence between deviatoric stress and mean principal stress is markedly affected by resilient modulus. Comparing to the influence of mean principal stress amplitude on resilient modulus, the loading sequence are more important for the evolution law of the resilient modulus. The primary reason lies in the fact that both the mean principal stress and deviatoric stress have the property of producing elastic strain. For the condition of synchronous loading and unloading, the elastic strain is produced by deviatoric stress and mean principal stress storing and releasing simultaneously, but for case of asynchronous loading and unloading, the released elastic strain because of the unloading of deviatoric stress is partially offset by the elastic strain which is produced by the mean principal stress loading process.

##### 3.2. Damping Ratio

The damping ratios of materials are an important parameter to measure the capacity of the soil to absorb the vibration energy during each cyclic loading process. For unfrozen soil, the mechanisms of material damping are contributed to friction between soil particles, strain rate effect, and nonlinear soil behavior [32]. From dynamic theory for soils mechanics, the hysteretic damping ratio of unfrozen soil is calculated as described in Figure 9, and its formulation is expressed as follows:where stands for the area of hysteresis loop and represents the area of shadow triangle. For frozen soil, damping primarily reflected the viscosity properties of ice and particles that was generally estimated by the damping ratio [8]. The calculation method damping ratio for frozen soil is shown in Figure 10, and the formulation is given bywhere stands for the area of the hysteresis loop ABDA and stands for the area of the shadow triangle ΔCEF. The damping ratio reflected the absorption ratio of energy during a single vibration. All the relationships between the standard logarithm of the number of load cycle and the damping ratios were fitted by the following equation:

The fitting parameters *a*, *b*, c, and *d* and *R*-square are given in Table 5.

Typical test results for the evolution law of damping ratio and under different dynamic stress paths with approximately equal half-length of cycle stress path or for different half-lengths with the same cycle load direction are shown in Figure 11. In the figure, the solid lines represent the optimal fitting of the experimental data. It can be discovered that the length of the cycle stress path and cyclic stress loading angle are the key elements which influence the shape of the curve, and when the cyclic stress loading angle increases and the length of cyclic stress path decreases, the damping ratio increases. Comparing with the length of dynamic stress paths, the loading angle of cyclic stress is the most important factor for affecting the damping ratio, and when the cyclic stress loading angle *α* > 90°, the dynamic damping ratio is obviously large than that when *α* < 90°. The curves of damping ratio get closer with equal amount increases in the length of the dynamic stress path under the same cyclic loading direction. With an increase in the number of load cycles, the damping ratio decreases and gets closer and closer under the same cyclic stress path direction. The results illustrated that the similar microstructures are induced by the loading in the same cyclic load direction.

The relationship between and with different phase differences of principal stress is shown in Figure 12. It is clearly seen that the damping ratio increases with a decrease of the phase difference of principle stress under the same cycle load numbers. As is shown in Figure 3(b), with increasing phase difference of principle stress, the rotation direction of the nonlinear cycle stress path is clockwise when *φ* < 180° and then the path direction turns counterclockwise when *φ* > 180°. When *φ* = 45°and *φ* = 270° or *φ* = 90°and *φ* = 225°, the stress path loci are similar but rotate in opposite directions, so the damping ratio has a great difference. Therefore, the loading order of normal and shear stress distinctly affects the damping ratio. The vibration energy absorption rate of frozen soil is clearly affected by the coupled effect of normal and shear stress.

Figure 13 describes the relationship between and under approximately equal amplitude of deviator stress. With an increase in the number of load cycles, the damping ratio decreases. The result indicates that the effect of the amplitude ratio on the damping ratio is not obvious when *ζ*^{ampl} > 0. But, when *ζ*^{ampl} < 0, the damping ratio is markedly lager than that in the former condition, and with an increase in cyclic loading, the angle value increases. According to the contrast the loading methods, it is found that, for synchronous loading of the mean principal stress and deviatoric stress (*ζ*^{ampl} > 0), the loading energy produced by the two loads is stored and released simultaneously. Similar as previously described about resilient modulus lager when *ζ*^{ampl} < 0, the axial deformation is mainly plastic deformation during the deviatoric stress loading process, and the recoverable elastic deformation is very limited to the deviatoric stress in the unloading process. Therefore, the comprehensive effect shows that the loading energy is mainly absorbed by the samples in the process of deviatoric stress loading and unloading process. From the findings, it can be concluded that the different vibration combination of deviatoric stress and mean principal stress is obviously affected the energy absorption rate, and this point is always ignored by investigators.

##### 3.3. Hysteresis Loop

This paper defined the hysteretic curve as the deviatoric stress-axial strain relationship under one cyclic loading. In this part, the hysteresis loops of each cyclic loading stress path when *N* = 10–11, 100–101, 1000–1001, and 10000–10001 are given in Figures 14–16. It can be seen that the area of hysteresis loops gets narrower and denser and the bottom opening gets smaller as the cycle number increases. It is evident that, for a single cycle, the ability of frozen soil to absorb vibration energy increases, the microstructure damage of the sample and unrecoverable strain decrease, and the samples are strengthened, and the deformation is the main elastic deformation with the increase of the number of load cycles. The testing result indicates that the shape of hysteresis loop varies greatly under different stress paths.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

**(m)**

**(n)**

**(o)**

**(p)**

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

**(k)**

**(l)**

**(m)**

**(n)**

**(o)**

**(p)**

The sharp of hysteresis loops of different cycle stress paths under the same half-length (*A* = 0.27 MPa) is shown in Figure 14. It is obvious that the sharps of hysteresis loops with = 1 and = −1 are obviously different. The cycle stress path of different amplitude ratios with the same length of cyclic loading paths is shown in Figure 3(a). The major difference between cyclic stress paths lies in the different combinations of normal and shear stress amplitude and order of the two loading and unloading. Shear and normal stress loading or unloading take place simultaneously when the amplitude ratio = 1, but when = −1, the two stresses’ vibrations are asynchronous. Comparing with the two kinds of hysteresis loop, it is found that the hysteresis loop is fatter when = −1 than when = 1. The curve of hysteresis loops bulge outward in the process of loading and unloading for the condition of = −1. This means that, with shear stress loading and normal stress unloading, the axial strain is increasing faster and faster. Then, with the shear stress unloading and normal stress loading, the corresponding axial strain changes to rebound. Compared with the case = 1, the slope of the hysteresis curves is relatively stable. Therefore, it may be concluded that shear stress loading and normal stress unloading process is more likely to cause microstructure damage, and the loading of normal stress in the shear stress unloading process can partly heal the damage.

It is shown in Figure 15 that the sharp of hysteresis loops which is caused by nonlinear cycle stress paths are significantly different with the linear condition. For different phase differences of principle stress, the sharps of hysteresis loops have a distinct difference, especially when the rotation direction of the dynamic stress path is opposite (such as *φ* = 45° and *φ* = 270°). With an increase in the phase difference, the area of hysteresis loop decreases under the same number of loading cycles. When *φ* = 270°, the deviatoric stress loading curve is below the unloading curve that leads to the calculated damping ratio being negative. Also, the negative increases with an increase in the number of loading cycles. Therefore, the normal stress obviously affected the shape of hysteresis loops and the state absorption or release of vibration energy. It is clear that the shape of the hysteresis loop obviously depends on the cyclic loading stress path and affected by the normal and shear stress coupling effect.

As is shown in Figure 16, which reflected the shape of hysteresis loops with different cyclic loading stress paths and cycle numbers under approximately equal deviatoric stress amplitude, the difference of the hysteresis loop is not obvious when *ζ*^{ampl} > 0. With an increase in the number of load cycles, the hysteresis loop becomes narrow and maintains a steady slope. This result illustrated that the major deformation of the sample in the last stage is elastic deformation. Similar with the former results, the shape of hysteresis loops is greatly different when *ζ*^{ampl} < 0. It is obviously fatter than the case of *ζ*^{ampl} > 0. So, it is illustrated that ability of energy dissipation is stronger during the asynchronous loading of deviatoric stress and mean principal stress than during synchronous loading.

In summary, the loading order of normal stress and shear stress or mean principal stress and deviatoric stress is a significant factor that affected the dynamic resilient modulus, damping ratio, and the hysteresis shape of frozen samples. The result clearly reflected that the initial isotropy microstructures of remold frozen soil samples induced different anisotropy structures by different cyclic loading stress paths. Therefore, the evolution of frozen soil microstructure with an increase in the number of load cycles has an obvious cyclic loading stress path dependence.

#### 4. Conclusions

In this paper, a series of cycle triaxial tests under variable confining pressure or constant confining pressure were conducted on frozen silt clay at the temperature of −6°C and frequency of 0.5 Hz to simulate various near-linear or nonlinear dynamic stress paths. The evolution of dynamic parameters and characteristics of the hysteresis loop with an increase in the number of cyclic loading under different dynamic stress paths are investigated. The main findings are summarized as follows:(1)With the continuous dynamic loading, the parameters resilient modulus and damping ratio rapidly decrease in the initial loading stage first and then go into the stability development stage. The shape of hysteresis loops is getting narrower and denser. As the cyclic loading angle of the stress path increases, the resilient modulus and damping ratio increase. Compared to the length of the dynamic stress path, the inclination angle of the cyclic stress path is the main factor to affect the resilient modulus and damping ratio .(2)The curves of and under the same cyclic stress path loading direction with different half-lengths of the cycle stress path are getting closer and closer. This result indicates that the same cyclic loading direction can induce the frozen samples to produce a similar anisotropy microstructure.(3)When the slope of the near-linear cycle stress path is opposite, or the rotation direction of the nonlinear cycle stress path is opposite, the dynamic parameters and the shape of the hysteresis loop are significantly different. This result reflects that the evolution of the dynamic characteristic of frozen soil is influenced by the coupling effect of normal and shear stress or deviatoric and mean principle stress. The mechanical property development of frozen soil has an obvious cyclic loading stress path dependence.

#### Data Availability

The data used in this study are provided in this manuscript.

#### Additional Points

1. A series of dynamic triaxial tests were conducted on frozen silt clay by using different cyclic stress paths2. Effect of cyclic stress paths on dynamic parameters and hysteresis loop characteristics of frozen silt clay were investigated in detail3. The influence of dynamic amplitude, cyclic stress path direction, phase difference, and rotation direction on dynamic parameters and hysteresis loop characteristics were systematically analyzed4. The evolution law of elastic-plastic deformation caused by normal and shear stress coupling is studied experimentally

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (grant no. 41571064) and the National Key Research and Development Program of China (no. 2017YFC0405101).