Abstract
The deformation behavior of rockfill is significant to the normal operation of concrete face rockfill dam. Considering both the nonlinear mechanical behavior and long-term rheological deformation, the E-ν model and modified Burgers model are coupled to describe the deformation behavior of the rockfill materials. The coupled E-ν and Burgers model contains numerous parameters with complex relationship, and an efficient and accurate inversion analysis is in demand. The sensitivity of the parameters in the coupled E-ν and modified Burgers is analyzed using the modified Morris method initially. Then, a new approach of parameter back analysis is proposed by combining back-propagation neutral network (BPNN) and Cuckoo Search (CS) algorithm. The numerical example shows that parameters , , and as well as are more sensitive to the deformation of the rockfill body. The inversion analysis for these four parameters and , , and as well as in modified Burgers model is performed by the CS-BPNN algorithm. The numerical results demonstrate that the parameters obtained with the proposed method are reasonable and its feasibility is validated.
1. Introduction
Concrete face rockfill dam (CFRD) is one of the most extensively constructed dams due to its strong topographic adaptability, economical engineering plan, and simple construction technology. The safe operation of concrete face slab, the most important water-retaining and antiseepage structure in CFRD, is closely related to the deformation behavior of the rockfill body. As the support structure of the face slab, the excessive and uneven deformation of rockfill may lead to interfacial hollowing or even cracking of the face plate, thus causing excessive leakage and even endangering the overall safety of the dam [1–3].
The rockfill material is a kind of granular mixtures, whose deformation behavior is the result of the interaction of instantaneous deformation and rheological deformation. Rational selection of constitutive model for rockfill mechanical behavior and determination of model parameters are the basis of rockfill deformation analysis. The instantaneous deformation of rockfill shows strong nonlinear characteristics. The E-ν [4] and E-B [5] models proposed by Duncan et al. have been used to describe this characteristic of rockfill in numerous studies and exhibit favorable consistency with the experimental results. Many scholars [6–9] have proposed models to describe the rheological behavior of rockfill, among which the Burgers model is advantaged for its simple expression and accurate description of viscoelasticity of rockfill materials. In consequence, the E-ν and modified Burgers models are coupled to describe the mechanical properties of the rockfill material in this study.
Calculated displacement from the coupled model may show higher dependence on certain parameters, i.e., subtle change of these parameters will lead to dramatic variation of calculated displacement, while other parameters present indistinctive effect on displacement, which can be determined by experimental results. Therefore, it is necessary to conduct sensitivity analysis before parametric inversion to eliminate insensitive parameter and promote computational efficiency. The modified Morris method [10–12] is adopted to identify the parameters of high sensitivity for parametric inversion of the coupled E-ν and Burgers model.
The coupled E-ν and Burgers model contains numerous parameters with complex relationship, whose inversion analysis is a typical nonlinear problem [13–18]. The traditional methods for parametric inversion are clear, but there are still shortcomings like inefficient calculation and poor accuracy. With the development of artificial intelligence, various intelligent algorithms have been introduced into the engineering reverse analysis and fruitful achievements have been presented [13–15, 19–22]. Guo et al. [15] performed deformation back analysis based on the response surface method and genetic optimization theory to sequentially calculate parameters of the creep and Duncan–Chang models for the Pubugou gravelly soil core rockfill dam. Gan et al. [20] proposed a new deformation back analysis method called MPSO-BP, which integrates a modified particle swarm optimization algorithm and neural network simulator, to reverse the creep model parameters of Jiudianxia CFRD. Zhou et al. [21] modified the genetic algorithm to solve the high-dimension multimodal and nonlinear optimal parameters inversion problem and validated this method in parametric analysis of E-B and Merchant creep model. The back-propagation neutral network (BPNN), with strong nonlinear capability, is employed to express the complicated relationship between the model parameters and simulated displacement. To overcome the disadvantages of BPNN like frequent local optimum, slow rate of convergence and possible overconvergence [23], Cuckoo Search (CS) is introduced in this study to optimize the thresholds and weights of BPNN to improve the efficiency and accuracy. Consequently, the CS-BPNN is established for parametric inversion analysis of rockfill materials.
The rest of this paper is decomposed as follows. A brief introduction of coupled E-ν and Burgers model is given in Section 2; then, the sensitivity analysis is conducted with the modified Morris method in Section 3; in Section 4, the detailed principles and procedures of CS-BPNN are illustrated and a verification example is presented for this method; parametric inversion and numerical simulation of Langyashan CFRD was conducted in the last section to validate the proposed method and it is concluded that the CS-BPNN is feasible in the inversion analysis of rockfill materials.
2. Coupled E-ν and Burgers Model
2.1. Duncan–Chang (E-ν) Model
The rockfill materials demonstrate strong nonlinearity on its stress-strain relationship. The E-ν model can be expressed by the following equations:where and are the tangent modulus and the bulk modulus of rockfill materials, respectively; , , and are the loading modulus number, damage ratio, and bulk modulus number, respectively; and are the shear modulus exponent and bulk modulus exponent, respectively; is the stress level; is the cohesion force; is the angle of internal friction; and are the nonlinear indexes of the friction angle; is the tangent Poisson ratio; is the consolidation stress; is the atmospheric pressure; and is the model parameter. In consequence, there are eleven parameters needed for E-ν model.
2.2. Modified Burgers Model
The rockfill body presents obvious rheological deformation during its construction and operation stages and the Burgers model is extensively used in its rheological calculation. The Burgers model [8] gets a Kelvin and a Maxwell model in series, as shown in Figure 1(a), and its creep equations are presented as follows:where is a constant; and are the elastic modulus of series and parallel springs, respectively; and are preset coefficients for the modified Burgers model; is the initial time; and . However, the Burgers model employs a linear function of time to describe the permanent deformation of rockfill and neglects the reduction of deformation increment caused by consolidation effect.
(a)
(b)
The modified Burgers model is employed to simulate the rheological property of rockfills in this study. Compared with Burgers model, the modified model adjusts its first damper with nonlinearity, i.e., the external damper that characterizes the viscous deformation is extended to a generalized one with its viscosity given as follows:
The schematics of modified Burgers model is shown in Figure 1(b), and its creep equations are as follows:
The modified Burgers model takes the consolidation effect of rockfill into consideration and compensate for the deficiency of conventional Burgers model, which will characterize the rheological property of rockfill with higher accuracy.
2.3. Coupling Mechanism
In the present study, the E-ν model and modified Burgers models are coupled to take both the nonlinear mechanical behavior and long-term rheological deformation into consideration. Therefore, the elastic modulus of series spring in modified Burgers model, , is given by the E-ν model; i.e., can be calculated by equations (1) and (2) for loading and unloading conditions, respectively. In this way, there are fifteen parameters needed for the coupled E-ν and Burgers model.
3. Parametric Sensitivity Analysis Based on the Modified Morris Method
The modified Morris method is a qualitative global screening method to identify the parameters with significant impact on the outputs of various models, whose feasibility has been validated in previous researches [10–12]. The key procedure of this method is to determine the input sample of parameters with rational variations in certain threshold range and calculate the sensitivity indexes of dependent variables to parametric changes. Morris method is advantaged in complex models with numerous parameters for its small computational costs and simple operation procedures.
An elementary model is defined with independent input parameters, , which varies in a -dimensional unit cube across selected levels. The elementary effect of the input factor with given is defined as follows:where is any value in and and are a random sample and its transformation within the parameter space, respectively.
Two sensitivity indexes, and , which represent the global and the higher order effect of each input parameter on output from the model, respectively, were proposed initially to identify the input parameters with appreciable significance. These two indicators can be obtained by the following equations:where is the number of sampling points in design space. Campolongo et al. [24] improved with the following equation:
Francos et al. [25] proved that the sensitivity indexes exhibit relatively higher accuracy when parameters change at a fixed step. In consequence, the modified Morris method adopts fixed variations of input parameters and improves the sensitivity indicators with the mean of multiple indexes. The sensitivity index demonstrating the sensitivity of each parameter to simulated displacements in coupled E-ν and Burgers model is defined as follows in this study:where is the sensitivity indicator of parameters to displacement of the gauging point; the positive value denotes positive correlation between the parameter and settlement and a bigger value indicates more important effect of the parameter on the simulated displacement, while the negative value means the opposite; is the calculated displacement of the simulation; is the simulated displacement with benchmark parameters; and is the number of simulations.
4. Back Analysis for Parameters Based on CS-BPNN Algorithm
Considering the numerous parameters for coupled E-ν and Burgers model and the complicated structure of practical engineering, it is difficult to obtain an explicit function for the model parameters and the simulated displacements. BPNN has demonstrated its strong expression capacity in nonlinear problems, which is therefore selected in this study to establish the nonlinear mechanical relationship between the mechanical parameters and calculated displacements. Meanwhile, aimed at the disadvantages of BPNN-like frequent local optimum, slow rate of convergence, and possible overconvergence, CS is introduced in this study to optimize the thresholds and weights of BPNN to enhance its robustness and promote the convergence.
4.1. BP Neutral Network
BPNN is a multilayered feedforward network that is capable of realizing complex nonlinear mapping, and it generally consists of an input layer, hidden layers, and an output layer with the thresholds and weights connecting neighboring layers [26]. The input layer is the calculated displacements at each gauging point, and the output layer is the mechanical parameters in the coupled E-ν and Burgers model in this study. The number of the hidden layers and their neurons are decided by trial calculation. The process of network training includes the forward propagation of input signal for actual output and the reverse transmission of error signal for weight correction, and this process is realized by continuous and recursive iterations until the accuracy of the outputs meets the requirements. The global learning rate with an added momentum is adopted in this study to speed up the learning process and avoid induced oscillations, i.e., a coefficient proportional to previous weighted variable is multiplied to every weighted regulator. The weighted regulating equation with added momentum is presented as follows:where is the iteration times, is the step length, is the total error of samples, is the connecting weight from node to node , and is the momentum term.
4.2. Cuckoo Search
CS is a latest metaheuristic optimization algorithm, which is enlightened by the obligate brood parasitic behavior of cuckoo species in combination with the Levy flight behavior of certain birds and insects and pioneered by Yang and Deb [27] in 2009. Levy flight [28, 29] essentially provides a random walk process with a power-law step-length distribution with a heavy tail for the global searching and the schematics for typical searching processes with 500 times of Levy flight and random walk are compared in Figure 2. It can be seen that the randomization with Levy flight is more efficient as the step length varies alternately; the long and short steps are responsible for searching global optimum and improve searching accuracy. Besides, CS is able to find the global optimum simultaneously if the number of nests is much higher than the number of local optimums. The basic steps of CS can be summarized as follows: Step 1: the objective function is established and input parameters of CS are determined, including the range of threshold, number of iterations, accuracy, etc. Due to different dimensions and magnitudes of parameters for back analysis, the relative error is taken as the objective function and expressed as follows: where is the number of simulations, is the number of parameters for back analysis, and and are the simulated and anticipant output parameters of the th simulation, respectively. Step 2: the initial nests, , are randomly generated. Each nest represents a set of weight and threshold of BPNN model, which is expressed as follows: where is the number of nests. The values of the objective function for each nest are calculated and the initial optimum and corresponding nest, , are decided. Step 3: all the nests are regenerated randomly by imitative Levy flight, which can be expressed as where and are the calculated values for the parameter in and iterations; is the step size and equals one in most cases; represents the entry-wise multiplication operator; is random searching vector following a Levy distribution, which can be given as where is the optimal nest in the iteration and and obey normal distribution. Equation (14) is essentially the stochastic equation for random walk. In general, a random walk is a Markov chain whose next location only depends on the current location (the first term in equation (14)) and the transition probability (the second term in equation (14)). Step 4: the egg laid by a cuckoo is discovered by the host bird with a probability . In this case, the host bird can either throw the egg away or abandon the nest and build a completely new nest. In consequence, a random number, , which is uniformly distributed in , is generated for every nest, and this process is simulated by the following equation: where and are two randomly selected and different solutions. In this way, a fraction of worse nests is abandoned and new ones are built, making sure the system will not be trapped in a local optimal. Step 5: then the values of the objective function for renewal nests are calculated and the current optimum and corresponding nest, , are decided. Then, the value of the objective function is compared with the previous one, and previous optimal nest, , will be replaced by the updated one, , if the relative error decreases; otherwise, the will be kept and will be abandoned. Step 6: steps 3 and 4 are recursively proceeded until the termination criteria or the maximum number of iterations is satisfied.
4.3. Verification Example
To verify the capability of the proposed inversion analysis method using the CS-BPNN algorithm, the vehicle flow and freight volume forecast model is selected.
Twenty groups of input and output data are chosen for the network training and the remaining 2 groups for prediction. In the training process, the topological structure of BPNN is and the Sigmoid function is employed as both the transfer and training functions. The CS algorithm is adopted to optimize BPNN, and its parameters are and , and the search range is set to .
To compare the accuracy between the CS-BPNN and the traditional BPNN, the latter is adopted with the same parameters. The convergence curves of the two algorithms are shown in Figure 3.
As illustrated in Figure 3, the mean square error (MSE) of the CS-BPNN and the traditional BPNN at the last iteration is 0.5483 and 0.5844, respectively, which indicates that the training of the samples using the CS-BPNN is more accurate than that using BPNN. Furthermore, the final iteration of the CS-BPNN and BPNN is 21 and 41, respectively, which means that the CS-BPNN has a faster convergence rate.
5. Process for Inversion Analysis
The procedures of back analysis for parameters in the coupled E-ν and Burgers model describing rockfills in CFRD based on CS-BPNN are summarized as follows, and the flowchart is illustrated in Figure 4: Step 1: the measuring deformation data are analyzed firstly and the characteristic measure points and time period are selected for inversion analysis Step 2: a finite element model is established and three-dimensional numerical simulations are carried out Step 3: parametric sensitivity analysis for the coupled E-ν and Burgers model is conducted based on the modified Morris method and the parameters with significant influence are selected for back analysis Step 4: the parameter samples for inversion are preprocessed and corresponding displacements at preselected measure points are calculated in an identical simulation condition Step 5: the simulated displacements and parameter samples obtained from Step 4 are input as the learning sets for network training to establish CS-BPNN for the back analysis of mechanical parameters Step 6: the measured displacements are imported to the established CS-BPNN, and the output is the identified parameters
6. Engineering Application
6.1. Engineering Overview
Langyashan pumped-storage power station is located in Anhui Province, the eastern part of China, which is responsible for peak regulation, frequency modulation, emergency reserve, and load following in power system. The comprehensive project consists of the reinforced concrete faced rockfill dam of the upper reservoir, channel system, underground powerhouse, exit trench, and the concrete gravity dam of downstream reservoir. Based on the storage capacity demand and topographic condition of the foundation, the axis of the reinforced concrete faced rockfill dam was designed as mansard and its general layout is illustrated in Figure 5. The crest elevation is 174.5 m with the maximum dam height of 64.5 m and the slope ratio of the upstream and downstream faces is 1 : 1.4. The dam is divided into slab area (1A), cushion area (2A), transition area (3A), upstream rockfill area (3B), downstream rockfill area (3C), and downstream protection area (3D), as shown in Figure 6. The upper reservoir is a daily regulation reservoir with the water level approximately ranging from 161.0 m to 171.80 m. Its normal water level, design flood level, check flood level, and dead water level are 171.80 m, 172.40 m, 172.60 m, and 150.00 m, respectively.
6.2. Observational Analysis of Dam Deformation
In this paper, the settlements of gauging points are used for back analysis of the mechanical parameters of the rockfill body. The settlement measurement for Langyashan CFRD is conducted by the hydraulic overflow settlement gauges, whose distribution is shown in Figure 5. The measured settlement process lines of gauging points LD1-1 to LD10-1 during the operation period are presented in Figure 7. The settlement distributions of Langyashan CFRD on selected megathermal and microthermal days in 2010 and 2016 are presented in Figure 8.
(a)
(b)
(c)
(d)
It can be seen as follows: (1) the measured settlements at gauging points increase during the operation period, which clearly reveals the rheological characteristics of rockfill body; (2) the measured settlement fluctuates periodically in annual cycle due to the influence of temperature, and the variation amplitude ranges from 6.5 mm to 8 mm; it should be mentioned that the thermal response of gauging points LD1-1 to LD10-1 is relatively more remarkable and they have smaller measured settlement because they are installed on the concrete wave wall of dam crest, which is more sensitive to temperature changes; (3) due to the construction quality and geological and topographical conditions, the dam crest close to left bank is obviously subsiding while dam crest close to the right bank has the tendency of slightly rising; (4) measured settlements of gauging points on rockfill body at the river bed are greater than that at the both bank slopes, as a result of larger dam body and stronger water pressure; and (5) the maximum settlement was measured at gauging point LD4-1 and reached 52.59 mm on December 14, 2018.
Considering the bigger disturbance on the gauging points at rockfill body close to right bank, measured settlements of LD1-1 to LD6-1 are selected for the inversion analysis of mechanical parameters, and the measured series before January 6, 2007, are abandoned to avoid the uncertain impact at construction stage.
6.3. Numerical Simulation Methods
A three-dimensional finite element mesh of Langyashan CFRD is established, which consists of 29490 elements and 49560 nodes, as shown in Figures 9 and 10. The element number of rockfill body and face slabs is 7806 and 13777, respectively. The bottom of the slope is fully fixed, and the lateral boundaries are constrained by vertical rollers. The coupled E-ν and Burgers model can be user-defined by UMAT subroutines for Abaqus/Standard solver. In addition, the construction process of Langyashan CFRD is simulated with the technique of death-birth meshing scheme built in Abaqus, which divides the construction into rockfilling, consolidation, and concreting stages. The construction time is decided according to Table 1 to restore the practical deformation filed of Langyashan CFRD. Furthermore, the general contact model is employed to describe the interfacial behavior between face slabs and cushion materials.
6.4. Sensitivity Analysis for Parameters of Rockfill Material
Considering that the rockfill body dominates the deformation of CFRD, this work only carries out back analysis on mechanical parameters of main rockfill materials. Due to the numerous parameters needed for the coupled E-ν and Burgers model, inversion analysis of all the parameters would be time-consuming, inefficient, and easily to trap into local optimum. Therefore, a sensitivity analysis based on the modified Morris method is conducted to improve the effectiveness and efficiency of parametric inversion.
Since those rockfill materials are generally granular mixtures, the value of cohesion, , is assigned as zero and denotes the atmospheric pressure and is a constant. Because the four parameters in modified Burgers model are significant to characterize the rheology property of rockfill, they are necessarily included in the inversion analysis. In consequence, the modified Morris method is adopted to analyze the sensitivity of the remaining 9 parameters in the E-ν model, and simulated settlements of gauging points LD1-1 to LD10-1 are taken as reference data.
The initial reference parameters are obtained with engineering analogy and experimental results of rockfill materials, which are shown in Table 2, and corresponding settlements at reference gauging points are calculated by numerical simulation. Then, each parameter is assigned another eight values changing at a fixed step to conduct the numerical simulation and sensitivity analysis. The sensitivity index of each parameter is calculated by equation (10), and the sensitivity indexes of parameters in E-ν model to settlements in reference gauging points are shown in Figure 11.
It can be clearly seen that , , and as well as are much more sensitive to settlements of the dam crest than the other parameters. These four parameters play dominant roles in simulated settlements at reference gauging points, of which takes the first place in sensitivity ranking. The absolute value of sensitivity index for the remaining five parameters approximately equals zero, which indicates their neglectable influence on the settlements of the rockfill body. In consequence, , , , and in E-ν model as well as , , , and in Burgers are selected for the following inversion analysis. In this way, the efficiency and accuracy of inversion analysis can be promoted dramatically.
6.5. Back Analysis for Parameters of Rockfill Material
Based on the reference parameters and corresponding variation range shown in Tables 2 and 3, eighty groups of selected parameters and corresponding simulated settlements at reference gauging points are obtained as the training sample of BPNN. In consequence, the complex nonlinear relationship between simulated displacements and material parameters can be established. The topological structure of BP neural network is as follows: the input layer has six nodes, i.e., the settlements of six reference gauging points; the hidden layer has seven neurons; the output layer has eight nodes, i.e., the parameters of rockfill materials to be inverted; and the transfer function between adjacent layers is Sigmoid function.
Then, CS is used to optimize the weights and thresholds of BPNN, and the detailed parameters of this algorithm are set as follows: the number of nests, ; the maximum number of iterations, ; and the accuracy . Eventually, the measured settlements of gauging points LD1-1 to LD6-1 shown in Table 3 are input into the established CS-BPNN, and the output is the inverted parameters, as shown in Table 4.
6.6. Analysis of Numerical Result
6.6.1. Analysis of Rockfill Deformation
Based on the inversed parameters in Table 4, displacements of gauging points LD3-1 to LD6-1 on January 19, 2016, and July 17, 2016, are calculated with numerical modeling. As shown in Figure 12, they are basically identical to the measured data. Figure 13 shows both the simulated and the measured value of gauging points LD3-1 to LD6-1 settlements and horizontal displacements during the operation stage. It reveals that simulated results closely mirror the measured results, in turn demonstrating that the parameters in coupled E-ν and Burgers model are reasonable. In addition, during the first 10 years of operation, the settlement rate of the dam body is relatively high, while the settlement rate gradually decreases and tends to be stable after 2015.
(a)
(b)
(a)
(b)
(c)
(d)
Figure 14 presents the settlement distributions of section A-A at the end of construction and first impoundment to normal water level. Figures 15 and 16 present the settlement distribution of section A-A and section B-B during the operation period, respectively. As it can be seen, the maximum settlement occurs at approximately two-thirds of the dam height rather than at the top of the dam due to the step loading of dam construction. In addition, the impounded water imposes hydraulic pressure on the face slabs and results in larger settlement of the cushion material and rockfill at upstream side, which exhibits the greater influence of instantaneous deformation. The settlement in both sections A-A and B-B increases with the running time, which embodies the rheological deformation of the rockfill and is consistent with the deformation law reflected by the process lines.
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
6.6.2. Analysis of Slab Stress
Figure 17 shows the distribution of stress along (S1), cross (S2), and perpendicular to (S3) the slope on the upper surface of the face slab on January 19, 2019, respectively, and Figure 18 shows the distribution of maximum and minimum principal stresses of face slabs on January 19, 2019, respectively. Table 5 presents the values and locations of the maximum tension and compression of aforementioned stress states. The slabs are numbered #1 to #44 from the left to right bank.
(a)
(b)
(c)
(a)
(b)
It can be seen that most area of face slabs suffers from compressive stress (−), while the partial area of the sides, top, and bottom of slabs suffers from tensile stress (+) in all the three directions. It should be mentioned that the tensile stress occurs at the connected regions of several slabs in S1 and the compressive stress tends to increase going downward along the slope at slabs of the riverbed. Stress distribution of S2 and S3 is more even except the tensile areas. The maximum tensile and compressive stress all occur at the bottoms of the concrete slabs, which can be seen from both Figure 17 and Table 5; while the distribution of principal stress is quite even at the upper areas of the slabs.
7. Conclusion
In this study, inversion analysis for coupled E-ν and Burgers model parameters is carried out by CS-BPNN based on the sensitivity analysis with the modified Morris method. The main points that can be concluded from this study as follows:(1)The E-ν model and modified Burgers model are coupled to describe rockfill deformation of CFRD to take both the nonlinear mechanical behavior and long-term rheological deformation into consideration.(2)Sensitivity analysis of parameters in coupled E-ν and Burgers model is conducted with the modified Morris method and numerical simulation in FEM. , , , , , , and as well as are selected for back analysis due to relatively stronger sensitivity.(3)The CS algorithm is used to optimize BPNN and the CS-BPNN algorithm is applied in parameter inversion of rockfill materials in Langyashan CFRD. The result demonstrates the feasibility of the proposed algorithm. Furthermore, coupled E-ν and Burgers model and its parameters are proved to be reasonable with the numerical results.
Abbreviations
| CFRD: | Concrete face rockfill dam |
| CS: | Cuckoo search |
| BPNN: | Back-propagation neutral network. |
Data Availability
The data are not available because the project is confidential.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
CS was responsible for conceptualization of the study, CS and YC were responsible for methodology, CS and XQ performed software analysis, YC, CS, and XQ were involved in validating the study, XQ was responsible for formal analysis, YC and XQ performed data curation, YC and XQ were involved in preparing the original draft, YC and XQ reviewed and edited the manuscript, and CG supervised the study.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 51739003, 51579085, 51779086, 51579086, 51379068, 51579083, and 51609074), National Key R&D Program of China (2018YFC0407104, 2018YFC1508603, 2018YFC0407101, and 2016YFC0401601), Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (YS11001), Postgraduate Research and Practice Innovation Program of Jiangsu Province (Grant no. KYZZ16_0283 and KYCX17_0424), Special Project Funded of National Key Laboratory (20165042112), and Key R&D Program of Guangxi (AB17195074).