Research on a combinatorial control method for coaxial rotor aircraft based on sliding mode

Yi-rn Wei , Hong-bin Deng , Zhen-hu Pn , Ke-wei Li , Hn Chen

a School of Mechatronical Engineering, Beijing Institute of Technology, Beijing, China

b School of Automation, Beijing Institute of Technology, Beijing, China

c Beijing Institute of Computer and Electronics Application, Beijing, China

d Hong Kong Polytechnic University, Hong Kong, China

Keywords:Coaxial rotor aircraft Sliding mode control Position and attitude tracking Simulation Experiment

ABSTRACTAiming at the position and attitude tracking of coaxial rotor aircraft (CRA), this paper proposes a combinatorial control method of sliding mode control (SMC) coupled with proportional-integralderivative control (PIDC).Considering the complete description of flight dynamics, aerodynamics and airflow interference, the dynamical model of CRA is established.The dynamical model is simplified according to the actual flight, then the simplified dynamical model is divided into two subsystems: a fully-actuated subsystem and an under-actuated subsystem.The controller of the fully-actuated subsystem consists of a SMC controller coupled with a rate bounded PIDC controller,while the controller of the under-actuated subsystem is composed of a SMC controller.The sliding manifold is defined by combining the position and velocity tracking errors of the state variables for each subsystem.Lyapunov stability theory is used to verify the stability of the sliding mode controller, which ensures that all state trajectories of the system can reach and stay on the sliding mode surface, the uncertainty and external interference of the model are compensated.Simulation and experiment compared with the conventional PIDC are carried out, the results demonstrate the effectiveness and the robustness of the proposed control method of this paper.

1.Introduction

Unmanned aerial vehicles (UAVs) have been applied in many fields[1-8].For instance,target location and tracking[9],operation monitoring of traffic data[10],inspection of network signal line or long-distance power line[11],searching for victims and delivering relief supplies[12-14].As a kind of special small UAV,coaxial rotor aircraft has many attractive advantages,such as compact structure without tail,capability of vertical take-off and landing(VTOL),and excellent hover performance.Moreover, the coaxial-rotor exhibits better thrust and efficiency compared to single-rotor or quad-rotor configurations under the condition of the same rotor area[15-18].However, its position and attitude are easily affected by the wind,aircraft structure, variation of parameters and external disturbances [19,20].Besides, the stability and tracking control of CRA is challenging due to the strong coupling of the system, the nonlinearity of the dynamic model, under-actuated and multiple input/output variables [21-23].Therefore, maintaining the flight stability of CRA is one of the most important problems to be solved in its application.It can be foreseen that CRA will be applied in both military and civilian fields if its attitude and trajectory can be controlled satisfactorily.As a typical coaxial-rotor aircraft, the three-dimensional model of CRA used in this paper is shown in Fig.1.

In previous literature, a variety of linear/nonlinear control approaches are proposed to meet the increasing demands of aircraft.Linear control methods such as PIDC[24],linear quadratic regulator(LQR) [25,26], and dynamic output feedback control [27].It is of practical significance to use linear control method while the nonlinear degree of the actual system is not serious.However,if the nonlinear degree of the actual system is serious,linear control will be less applicable [28,29].In engineering practice, linear control may fail to meet the accuracy and real-time requirements of a complex system.In this case,a variety of nonlinear control methods are proposed [30], such as the backstepping-based control,nonlinear H∞ control and sliding mode control (SMC).The backstepping-based control has obvious advantages in the design of robust or adaptive controllers for uncertain systems, especially when interference or uncertainty does not meet the matching conditions.Unfortunately,it is difficult to reflect the actual system due to its low correlation with the controlled object [31].The H∞controller achieves the desired performance by adjusting the frequency domain characteristics of the system,and fully considering the influence brought by the uncertainties of the system, the representation of H∞control is simple, but simple H∞controller cannot satisfy the accuracy requirements, and complex H∞controller (combined with some other algorithms, e.g.neural net algorithm) will increase the complexity of the system and the amount of computations [32,33].The basic principle of SMC is to drive and constrain the trajectory of the system so that it approaches and stays near the predetermined sliding surface, it has the advantages of low sensitivity to interference and parameter variation,low requirements for precise modeling and without need for online system identification [34-36].What’s more, the SMC controller has rapid response speed and good performance even under unsteady aerodynamic [37-39].In addition to these advantages, another reason for being selected in this study is that the SMC controller can be easily implemented [40].

Fig.1.Three-dimensional model of CRA used in this paper.

In recent decades, SMC is recognized as an effective tool to design robust control law of aircraft in the presence of uncertain conditions[41,42].In the early stage,controllers based on SMC are designed in Refs.[36,37]to stabilize the attitude of quadrotor aircraft.Simulation results show that the control performance is accepted, but due to the chattering phenomenon, the flight performance is not as good as expected.To alleviate this phenomenon,adaptive control is added to SMC in Ref.[38], simulation results show that the adaptive SMC controller has good performance under the noise environment.With the development of SMC, more and more new control methods based on SMC are proposed.For instance, a second-order SMC with the embedment of a super twisting algorithm is proposed in Ref.[39]to address the trajectory tracking control problem,by comparing simulation results with the pioneering work of[36],the superior performance of this algorithm is demonstrated.A nonlinear control based on SMC for the position and attitude tracking control of a small coaxial-rotor UAV subjected to uncertainties and aerodynamic disturbances is proposed in Ref.[37].Simulation results show that this control method can achieve good performance, but the controller parameters are difficult to obtain only by experience, and the controller accuracy needs to be improved.With the emergence of a variety of optimization algorithms, some SMC control algorithms combined with optimization algorithms are proposed.For example, to reach the asymptotic convergence of SMC,a technique of SMC for finite-time stabilization of UAV systems with parametric uncertainties is designed in Ref.[40],simulation results are presented to exhibit the helpfulness of the offered technique compared to the previous methods.Radial Base Function Neural Network (RBFNN), Fuzzy Logic Control(FLC)and SMC are used to design a controller for a Six Degrees of Freedom (6-DOF) coaxial octorotor helicopter in Ref.[41].The simulation results show that the control method can effectively alleviate the chattering phenomenon and has good tracking performance.

As a variety of SMC algorithms are proposed, here are some issues to take into consideration:

●The accuracy of control stability cannot be guaranteed if the SMC algorithm is too simple;

●The complex control algorithm is difficult to realize in practical projects;

●Most of the algorithms are verified only by numerical simulation, not by experimental tests.

The main contribution of this paper consists in proposing a simple, efficient, and practical SMC-PIDC algorithm for stabilizing CRA, which is listed as follows:

●Aiming at the stability control of CRA,the SMC-PIDC control algorithm is proposed.This algorithm has low computational complexity and strong robustness, easy to be used in the practical engineering.It is of great significance to improve and enhance the flight stability of CRA.

●Construction of a new controller based on SMC-PIDC for the position and attitude control of CRA.Under the action of this controller,the system has good stability,fast response speed,small overshoot and high steady-state precision.The robustness of the designed controller is demonstrated, and the effectiveness of the proposed control scheme is justified.

●From the comparison with the PIDC,the SMC-PIDC is proved to have a strong ability to reduce the chattering and the tracking error,maintain the stability of flight and obtain good control performance.

This paper is organized as: the mathematical model of CRA is presented in Section 2.The controller design and stability proof are presented in Section 3.The numerical simulation results are presented in Section 4.The experimental test results are presented in Section 5.Finally, the conclusions are presented in Section 6.

2.Mathematical model of CRA

The mathematical model is essential for designing a control algorithm with satisfactory performance.CRA has two main rotors driven by brushless DC motors, symmetrically installed between the rotors.The steering mechanism adopts a semi-differential mode to realize course control by changing the lift plane of the lower rotor.From top-down, the upper rotor rotating by counterclockwise while the lower rotor rotating by clockwise, gyroscopic effects and aerodynamic torques tend to achieve a dynamic balance in hover state[42,43].The connection of the control mechanism is shown in Fig.2.

Fig.2.The control mechanism.

Under the control of input signal generated by the steering engine, the rotating part and the non-rotating part of the swashplate tilt in the expected direction together [44].There is a connecting rod between the rotating part and the variable pitch hinge of the blade, so the tilt of the swashplate will lead to a periodic change of blade pitch, making the aerodynamic asymmetry in the direction of rotation, the thrust vector of the rotor tilts in the expected direction,to achieve controlling the lateral and longitudinal flight of aircraft,while the altitude of the aircraft can be controlled by the rotating speed of the upper and lower rotors[45].

The aircraft is illustrated in Fig.3.Its dynamical model is set up by the body-frame（） and the earth-frame（）.is the distance between the lower rotor and the center of gravity of the aircraft,is rotor radius,，v，are velocity components along aircraft axes.

2.1.Model without rotor interference or flapping motion

Take V =（，，，，v，）,η =（φ，θ，ψ）,X =（，，）,where,，，are the angular velocity components along aircraft axes,，v，are the velocity components along aircraft axes, φ，θ，ψ are Euler angles.

Based on the Euler-Poincare Equation[42],the dynamical model of CRA can be established.The general form of the equation is given by

where,M is inertia matrix, C is a transformation matrix, F is the sum of the aerodynamic in the body coordinate system,gravity,and control inputs.To simplify writing, make the following simplification:s（） =sin（）,c（） =cos（）.Partial matrices in Eq.(1)is simplified as

Fig.3.Coordinate system and structure of CRA.

The transformation between the frame coordinate system and the ground coordinate system is realized through the Euler matrix

The thrust situation at the barycenter of the aircraft is given by

where=+,=,=,= ξ（-）,is acceleration due to gravity.Among them,is the vertical height input control variable,is the roll input control variable,is the pitch input control input,is the yaw input control variable,is the distance between the thrust plane and the origin of the coordinate system and ξ is the rotary torque coefficient,which is related to the resistance of the rotor and the friction coefficient of the motor shaft.

According to Blade-Element theory [43], assume that the angle of attack is a fixed value.The thrust of the rotor is established by

where,α is blade section angle of attack,ρis air density,()is a function of blade chord varying with blade element position,is thrust coefficient, Ω is rotor angular velocity.According to Eq.(5),the dynamic model of the system can be obtained as

The K（=1～6）given above are the drag coefficients of different motion forms since these drags must be considered at high speeds[46].

From Eq.(7), it can be found that the position of the aircraft is simultaneously affected by the three input controls,and,which means strong coupling in the control model of the aircraft.Besides,there are some inevitable problems during the movement of the aircraft, such as the aerodynamic interference between rotors, dynamic balance of torque cannot be guaranteed when the aircraft is unstable.As these problems aggravate the control error caused by the coupling of the aircraft control system,the coupling problem of the control system must be addressed.

In practical application,roll angle and pitch angle are small,the deflection angle of the lower rotor plane is very small, that is δ≤5,then the corresponding projection component ofin the three directions

The k（=1 ～6） given above are the coefficients related to the parameters of CRA,while= -/,= -/,= -/,= -/I,= -/I,= -/I.

To further simplify the model, take=/,=/ I,=/I,=/I.Eq.(9) can be written as

2.2.The model with rotor interference and flapping motion

To reveal the aerodynamic interference between rotors, based on the Blade-Element theory and Pitt-Peters dynamic inflow model[43,44], the motion model of inflow disturbance and flapping is established.Compared with the analysis methods base on generalized vortex theory and CFD, this method requires less computation and the obtained result is in good agreement with the experimental data[47].

It should be noted that the established model is based on the following assumptions:(1)the blades are considered to be rigid;(2)the truncation effect at the blade root,tip loss,and extension of the flapping hinge are neglected; (3) the unsteady effect is neglected;(4)the inflow velocity is considered to be uniformly distributed in the plane of the blade disk.

The approximate inflow ratio at a blade element can be defined as

where, λ，λare inflow ratio corresponding to the induced velocity and the aircraft velocity, respectively.The subscripts u and l stand for upper rotor and lower rotor in coaxial configuration.vis the projection of the velocity of the aircraft along the-axis in the frame of reference, vis the induced velocity from infinity to a blade element.

According to the dynamic inflow model [48], the relationship between induced velocity and tension coefficient is

where, λ= （λλλ）,= （）,is the roll moment coefficient andis the pitch moment coefficient.M，V，Lare the inertia matrix, mass flow parameter matrix,and static transfer matrix of inflow dynamics, respectively.

The interaction of induced velocity is

where,，∈｛u，l｝≠.λ，λ，λare time-averaged, one order horizontal, and one order longitudinal components of the induced inflow ratio, respectively.

Assuming that the induced velocity is uniformly distributed in the paddle plane.The thrust coefficient and torque coefficient can be obtained

where, σ=/（π） is rotor solidity (total solidity for coaxial rotor),θis the pitch of blade element,is the number of blades in rotor,is drag coefficient of airfoil.

Therefore,the thrust and the torque of a single rotor are

The subscriptrepresents the upper rotor()or the lower rotor(),is disk area.

Considering the rigidity of the blades,the orthogonal torques of the flapping motion in pitch channel and roll channel are

where,is blade stiffness coefficient,γ is blade lock number,is the moment of inertia of wave motion, αis the vertical angle of wave motion and αis the horizontal angle of wave motion.

After taking the rotor interference and the flapping motion into consideration, some equations in the dynamical model can be modified as follows

3.Controller design

This section mainly introduces the flight controller designed with SMC-PIDC, as shown in Fig.4.The control system is divided into a fully-actuated subsystem and an under-actuated subsystem[49].The fully-actuated subsystem is controlled by a rate bounded PID controller and an SMC controller.The under-actuated subsystem is stabilized by the algorithm of SMC[50].The whole controller is designed to maintain the required yaw angle and move to the target point while keeping the pitch angle and roll angle close to zero.

The fully-actuated subsystem is

The control objectives are:→,→,→, φ→0, θ→0,ψ→ψ.With four control inputs,the tracking ofand ψ in the fullyactuated subsystem can be realized,as well as the tracking ofand, the stability of φ and θ of the under-actuated subsystem.

3.1.The control law for fully-actuated subsystem

For the fully-actuated subsystem(18),the controller is required to driveand ψ to their desired valuesand ψ.Considering the design of SMC law for the under-actuated subsystem, to design abounded, it is necessary to ensure thatis non-singular [30].For subsystem ¨=cφcθ-+˙, the control input is given by

Taking=-, the PID controller is adopted in Eq.(20) to obtain→,and→z.At the same time,the signalneeds to be stable,and the parameters,,,K,Kand Kneed to be adjusted several times.

3.2.The control law for the under-actuated subsystem

By setting a positive value η(η>0 and η≠), one obtains

The transformation of the under-actuated system is

Fig.4.Control structure.

where x∈ℝ（=1，2，3，4） are system states; u∈ℝis the control input; fand f:ℝ→ℝb:ℝ→ℝ, it’s invertible and composed of the nonlinear smooth vector functions;d∈ℝ（=1，2，3） is disturbance.The block diagram is depicted in Fig.5.

3.3.Error and stability analysis for under-actuated subsystem

Considering the error equation

Fig.5.The block diagram of the under-actuated system.

Take E= ［e，e，e］，E= ［e，e］.It is more convenient to exclude eand/or efrom Eand Ein later definitions and stability proof.The following assumptions are also needed to satisfy the matched disturbance condition of the SMC design, as the disturbance dis bounded, ∂f/∂x=0 and ∂f/∂xare invertible[51].

Then the switching surface is defined as s =e+e+e+e, where c>0,= 1，2，3.

The SMC controller consists of two parts:the equivalent control part uand the switching control part u.The equivalent control on the manifold s=can be calculated from ˙s =, one obtains

The switching control part is designed to derive s to the manifold s =, one obtains

judgment can be obtained

where,λ（Q）is the minimum eigenvalue of the positive definite matrix Q , Q is the generalized coordinate vector of the center of mass relative to the earth coordinate system.

Substitute the control law into the above equations, then Eq.(36) can be expressed as

It can be obtained from Eq.(38)that the Lyapunov function will globally asymptotic convergent to zero.Therefore, under the control law in Eq.(33),the actuated states will be driven to the sliding mode surface,the system will reach and thereafter stay on s=0 in a finite amount of time, the effect of the disturbance will also be eliminated.

4.Numerical simulation

In the simulation, the initial state of the system is defined as=［，，，，v，，φ，θ，ψ］,the initial values of the system are taken as= ［0，0，0，0，0，0，0，0，0］.To demonstrate the performance of the SMC-PIDC, the PIDC is selected as the control group.Under the control of these two algorithms, the experiments of fixed-point flight are designed [52].The initial point is (0,0,0) and the target point is(3,3,3),the attitude instruction of the target attitude is(0,0,π/3).These two groups of experiments are carried out under the same experimental parameters [53].

The dynamic parameters of the aircraft are shown in Table 1.

Table 1Aircraft dynamics parameters.

The parameters used in the SMC-PIDC algorithm are shown in Table 2.

Table 2Controller parameter list.

The position information is shown in the form of a threedimensional trajectory in Fig.6.In terms of the smoothness of the trajectory and the degree of path optimization, SMC-PIDC is superior to PIDC.To clearly express the motion in three directions,the details of position information are shown in Fig.7.Under the control of PID (as shown in Fig.7(a)), the transient time ofdirection are 7 s, 8 s, 5 s, respectively, the overshoots ofdirection are 2.2 m,1 m,2.8 m,respectively.And the performance of position holding is not good, a constant position fluctuation could be found in the x-direction.While under the control of SMC-PIDC(as shown in Fig.7(b)), the transient time ofdirection are 3 s,3 s,6 s,respectively,the overshoots ofdirection are 0 m,0 m, 0.9 m, respectively.And there is no jitter in the three directions.It can be concluded from the above, under the control of SMC-PIDC, the flight control system can follow the command instructions quickly and keep high precision.SMC-PIDC has the advantages of a stable flight path, good control performance, strong stability, and so on.

The attitude information of the simulations can be seen in Fig.8.Under the control of PID(as shown in Fig.8(a)),jitter exists in pitchand roll direction, meaning that the aircraft vibrates while hovering, especially in the x-direction.While under the control of SMCPIDC (as shown in Fig.8(b)), the aircraft hovers without jitter,what’s more,the time needed to reach a stable state is shorter than that of PIDC.

5.Prototype experiment

The structure of CRA prototype used in this experiment is shown in Fig.9,which is made according to the 3D model in Fig.1.All the components are fixed on the hollow carbon fiber tube from top to bottom,power cord and signal cord are arranged inside the carbon fiber tube.To ensure the accuracy and performance of the prototype, the power transmission part is made of aluminum alloy, the connecting rod part is made of nylon, the motor pinion is made of copper, which is coordinated with the large aluminum alloy gear.The body is composed of carbon fiber board and aluminum alloy hollow tube, and equipped with a foot frame.The total weight of the machine is about 2 kg.

Although the prototype is made according to the 3D model in Fig.1., the parameters of the prototype are slightly different from those of the 3D model due to the processing and selection of parts.Prototype’s parameters are shown in Table 3.

Table 3Parameters of prototype.

Several flight states during the flight experiment are recorded in Fig.10.The left photo shows the take-off state.Four brackets are installed at the bottom of the aircraft to ensure a smooth takeoff.After the zero-seeking of the inertial sensors is completed,the takeoff command can be issued[54].Although there is unstable air flow in the experimental environment, there is no significant influence on the stability of the aircraft from the experimental data.The right photo shows the hovering state of CRA, from which it can be seen that the attitude of the aircraft is self-stable.This shows that the attitude control algorithm of the aircraft can adjust the attitude of the aircraft under the disturbance of small wind to achieve a satisfactory attitude.(The recorded wind speed is about 1 m/s~1.5 m/s at a height of 3 m above the ground, in the southerly direction.However, the wind speed at a height of 30 m may be slightly higher, but the tests are in the same place, and the wind conditions is basically similar.) During the test, it can be observed that SMC-PIDC has a better stability than PIDC, which can be explained by the test data below.

Fig.6.Position information of PIDC and SMC-PIDC.(a)PIDC; (b) SMC-PIDC.

Fig.7.Position information of PIDC and SMC-PIDC.(a)PIDC; (b) SMC-PIDC.

Fig.8.Attitude information of PIDC and SMC-PIDC.(a)PIDC; (b) SMC-PIDC.

Fig.9.Prototype of CRA in this experiment.

To ensure that the experimental conditions of each tests are consistent.CRA takes off from(0,0,0)and hovers at(30,10,10),the coordinates are in meters.In order to verify the hovering stability of CRA with different algorithms,the altitude data and attitude data of CRA during the hovering state are sorted out.The two figures in Fig.11 are the attitude measured by the above hovering flight experiment in the presence of PIDC and SMC-PIDC,respectively.On the whole,the jitter of roll angle and pitch angle was severe,while the jitter of yaw angle is slightly, mainly due to the impact of environmental factors (wind and magnetic field, etc.) and aircraft structure (connecting mechanism, rotor), also related to the dynamic system (motor, transmission mechanism) which will produce high-frequency vibration in the process of generating and transferring power.Under the control of PIDC, the amplitudes of roll angle and pitch angle were controlled within±5,which were decreased to ±3and the jitter frequency was significantly optimized under the control of SMC-PIDC.The yaw angle changes more gently than the others, but its amplitude is larger.The reason for this phenomenon is that the aircraft can rotate at any angle in the yaw direction while flying.But the yaw angle needs to be controlled steadily during the mission.Under an unbalanced moment,the yaw angle is easy to vary greatly, but the torque on the yaw direction will not change quickly, resulting in the jitter frequency of yaw angle significantly smaller than the jitter frequency of roll angle andpitch angle.Under the control of PIDC (as shown in Fig.11(a)), the amplitude variation range of yaw angle is （-200，100）.Once a larger disturbance appears, it will be difficult to adjust the yaw angle back to steady-state,which will bring great disadvantage for the mission to be performed,therefore,there is an urgent need for algorithms that can stabilize the attitude of aircraft.It can be seen that under the control of SMC-PIDC (as shown in Fig.11(b)), the amplitude variation range of yaw angle during the whole flight is（-50，-28）,it can be obtained that the range of variation is 7.3%of that under the control of PIDC,what’s more,the jitter frequency of the yaw angle is reduced,indicating that the SMC-PIDC has a good stability control on yaw angle.The comparison between Fig.11(a)and (b) shows that the control performance of SMC-PIDC is superior, which can greatly improve the attitude stability of the CRA.

Fig.10.Prototype flight experiment.

Fig.11.Attitude information of aircraft in flight experiment.(a) PIDC; (b) SMC-PIDC.

However, vibration is inevitable during the flight, mainly because the aircraft’s power system(motor,rotor)generates highfrequency vibration during the process of generating and transmitting power,and the surrounding environment(wind,magnetic field, etc.) will affect the aircraft.It is acceptable to control this vibration within a certain range.

The position of flight test under the control of SMC-PIDC and PIDC is shown in Fig.12.The position curve of SMC-PIDC is significantly smoother than PIDC.Under the control of SMC-PIDC,all of the directions of X, Y and Z, the position fluctuation of the aircraft is less than ±0.3 m which is significantly less than that of PIDC,indicating that Under the control of SMC-PIDC,CRA is able to achieve a stable hover[55].

6.Conclusions

In this paper,based on a new type of CRA,we analyzed the flight principle and established the mathematical model of CRA, proposed the SMC-PIDC control method to improve the flight stability of CRA.To verify the performance and practical feasibility of this control algorithm, simulations and flight experiments are performed.From the comparison with the PIDC, the SMC-PIDC is proved to have a strong ability to reduce jitter and tracking error,maintain flight stability and obtain good control performance.

A simulation model consistent with the control model used in the actual flight test is established in MATLAB/Simulink, then the simulation results are compared and verified with the flight test results.In the analysis of simulation results,it was found that under the control of SMC-PIDC,the jitter in the hover state is reduced,the time it takes to reach a stable attitude is significantly shorter than that of PIDC, In the flight test, a prototype of CRA was made according to the designed parameters.Under the same experimental conditions, these two control algorithms are respectively written into the flight control chip for the fixed-point hover flight test.By comparing and analyzing the experimental data, the good control performance of the SMC-PIDC algorithm is further illustrated.

Fig.12.Position information of the aircraft.(a) PIDC; (b) SMC-PIDC.

Simulation and experiment show that SMC-PIDC is a simple and practical flight control algorithm.Under the control of SMC-PIDC,the speed of command following is accelerated, the maneuvering performance of CRA is improved, the attitude and hovering performance of CRA become more stable.The algorithm can achieve good control effect for nonlinear strong coupling systems such as CRA.In addition,the requirement of calculation is not high,it can be recognized that SMC-PIDC has great engineering application value.Future work will be carried out to optimize the algorithm based on the structure of the aircraft,further increase the attitude stability of the aircraft during flight state and hover state, strengthen the ability of the control system to resist external interference, and accelerate the application of CRA in daily life.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This work was supported by National Natural Science Foundation of china(Grant No.51774042).