#### Abstract

In the case of dynamic motion such as jumping, an important fact in sEMG (surface Electromyogram) signal based control on exoskeletons, myoelectric prostheses, and rehabilitation gait is that multichannel sEMG signals contain mass data and vary greatly with time, which makes it difficult to generate compliant gait. Inspired by the fact that muscle synergies leading to dimensionality reduction may simplify motor control and learning, this paper proposes a new approach to generate flexible gait based on muscle synergies extracted from sEMG signal. Two questions were discussed and solved, the first one concerning whether the same set of muscle synergies can explain the different phases of hopping movement with various velocities. The second one is about how to generate self-adapted gait with muscle synergies while alleviating model sensitivity to sEMG transient changes. From the experimental results, the proposed method shows good performance both in accuracy and in robustness for producing velocity-adapted vertical jumping gait. The method discussed in this paper provides a valuable reference for the sEMG-based control of bionic robot leg to generate human-like dynamic gait.

#### 1. Introduction

Surface Electromyogram (sEMG) is the electrical manifestation of muscular contractions, which reflexes plentiful neural control information. For its strong relationship with human motion pattern, sEMG has been taken as an ideal noninvasive control signal for exoskeletons [1], myoelectric prostheses [2], and biorobot [3] and moreover for the development of rehabilitation robots [4]. An important factor presenting in the sEMG-based biomechanical leg control, during the dynamic motion such as jumping and running, is the fact that multichannel sEMG signals contain mass data and vary greatly with time. Such makes it more difficult to generate compliant motion. It is more desired to generate compliance gait with sEMG signal, for exoskeletons, myoelectric prostheses, and rehabilitation gait.

It is hypothesized that the CNS (central nervous system) coordinates groups of muscles with specific activation balances and temporal profiles, to simplify the generation of intricate movements [5]. These building modules, known as muscle synergies, can be used as a small number of coactivation patterns to imitate the performance of movement [6]. It is very attractive to point out that these synergies make it possible for the motor intentions to be rapidly translated into muscle activation and the systems can learn and plan movements so fast [7, 8]. From the computational perspective, with muscle synergies leading to dimensionality reduction that simplifies motor control and learning, such observation has recently raised the interest of many researchers to develop control strategies in robotic and biomechanical application [9]. So, we suppose that it is reasonable to apply the muscle synergies to simplify the generation of compliant dynamic gait, such as hopping, based on sEMG for biomechanical leg control.

The generality of muscle synergies across different motor tasks has been illustrated in human walking, running [10], and forward and backward pedaling movement [11]. For example, previous work on human balance control has shown that the same muscle synergies can account for balance responses under different dynamic conditions: stepping and nonstepping postural responses [12]. Five modules satisfy the human walking in sagittal plane, while sixth module, which contributes primarily to mediolateral balance control and contralateral leg swing, is needed to satisfy the additional nonsagittal plane demands of 3D walking [13]. Although the same muscle synergies are used across multiple tasks, in some instances, new synergies may be recruited to accomplish a new behavior goal [14]. Like in the motion of human finger spelling, the recruitment of muscle synergies is correlated with common hand postures [15]. Thus, the first problem that would be discussed in this paper is whether the same set of muscle synergies can explain the different phases of dynamic gait with various velocities.

To generate the gait pattern with muscle synergies is related to the problem of identifying the forward relationship between sEMG and resultant joint movement. In the past few years, several contributions were proposed to predict the forward relationship between synergies and joint angles, joint torques, or end force. Like the approaches based on the linear models, Artemiadis and Kyriakopoulos used Linear Time-Invariant (LTI) model to take PCA synergy features as inputs to relate synergy features to anthropomorphic joint movements in three-dimensional space [16]. In order to predict the wrist intended activation of natural movement, synergies of both DOF (Degree of Freedom) were extracted at once, and three synergies can achieve simultaneous control when applied separately on each DOF [17].

Nonlinear model based approaches are not as reliant on robust synergy features as linear models and therefore are able to represent rational complex relationship between synergies and desired outputs. Hahne et al. systematically compared linear and nonlinear regression techniques for an independent, simultaneous, and proportional myoelectric control of wrist movements and got the results that the kernel ridge regression outperformed the other methods, but with higher computational costs [18]. A hybrid time-delayed artificial neural network was investigated to predict clenching movements during mastication from sEMG signals. Actual jaw motions and sEMG signals from the masticatory muscles were recorded and used as output and input, respectively, [19]. Muceli and Farina also adopted multiplayer perceptron (MLP) artificial neural networks to estimate kinematics of the hand wrist from EMG signals of the contralateral limb during mirrored bilateral movements in free space [20].

However, during the dynamic motion such as jumping and running, the multichannel sEMG signals vary greatly with time which may result in high risk for overfitting models to training data and frequent retraining [21]. Moreover, for the application of exoskeletons, myoelectric prostheses, and rehabilitation gait, it is always desirable to generate self-adapted gait with limited experimental data. Therefore, the second problem that needed to be solved in this paper is how to generate flexible continuous dynamic gait with limited sets of experimental data based on these extracted muscle synergies, while avoiding overfitting model and alleviating model sensitivity to sEMG transient changes.

Since vertical jumping is the fundamental movement pattern of dynamic motion, such as jumping, bouncing, and running, based on the above discussion, we will focus on generating flexible gait of vertical jumping with sEMG signals for biomechanical leg. The paper is organized as follows: In Section 2, experimental protocol and vertical jumping motion are introduced. How to extract muscle synergies is explained in Section 3. To generate the coordinating variation of joint angles, a synergy-based computational framework built on the Fuzzy Wavelet Neural Networks is introduced in Section 4. Section 5 shows the experimental simulation results of the generalization gaits for vertical jumping motion. The advantages of this approach are the possibilities to generate a self-adaptive gait pattern of vertical jumping with limited number of experimental data, which is very meaningful for the sEMG-based robotic leg application.

#### 2. Experimental Protocol for Vertical Jumping

Six male volunteers with no lower-limb functional limitations or neuromuscular disorders participated in the study. The participants’ physical characteristics are the following: age, years; height, m; and body mass, kg (mean ± SD). The protocol was approved by the local ethical committee and accorded with the guidelines set out in the Declaration of Helsinki (1983).

##### 2.1. Experimental Protocol

Six volunteer subjects were trained for vertical jumping with uniform and variable speed before the experiment. The experiment contained 16 trials. For the first 8 trials, each containing 15 jumps, the subjects were asked to hop with fast speed for 4 trials and continuously jump with slow speed in the next 4 trials. In the following 8 trials, each containing 20 jumps, for the first 4 trials the subjects were asked to continuously hop with five fast and five slow hops, which were carried out alternatively. For the last 4 trials, slow speed and fast speed hops were carried alternatively in the same way. In this case, not only the muscles coactivities of jumping motion with different velocities but also the muscle activation for changing velocity could be recorded.

Three-dimensional kinematic data were collected using VICON 10-camera motion capture system at a frequency of 500 Hz. AMTI 3D force platform was used to keep track of plantar force with sampling frequency of 500 Hz. The BioVision 8 channel sEMG device with sampling frequency of 1000 Hz was used to record the activities of seven leg muscles: tibialis anterior, gastrocnemius, soleus, vastus medialis, rectus femoris, gluteus maximus, and biceps femoris. The eighth channel of sEMG equipment was used for synchronization with motion capture data. The experimental environment is shown in Figure 1.

##### 2.2. Vertical Jumping Movement

###### 2.2.1. Phase Analysis of Vertical Jumping Movement

Because of the changing of impact force imposed on the pelma, the jumping movement is usually divided into stance phase and flight (swing) phase. Figure 2 describes the transition of each phase. During the stance phase, all the joint angles firstly compress. This compression stage begins at the moment of foot touching the floor and lasts until the vertical velocity equals zero. This process ensures the leg absorbing the ground impact under high speed and being ready for the next hopping movement. Then, the leg extends to provide enough force for the requirement of flight phase. This extension stage ends until there is no ground reaction force. Based on the experiment results, it is important to point out that the maximum plantar force always appears after the moment of maximum joint compression. During the flight (swing) phase, the joint angles change to maintain the whole body in balance state and jump with certain velocity.

###### 2.2.2. Influence of Jumping Velocity on Vertical Jumping Gait

To evaluate if the velocity exerts important influence on jumping gait, we will discuss the relevance of jumping rhythm with the joint angles. Taking the motion capture data of one volunteer as an example, the available vertical jumping rhythm is between 430 ms and 660 ms after all the experimental trials finished. The joint angles were grouped according to the same rhythm for each 10 ms interval, and all the experimental data could be divided into 23 groups.

The box-whisker plot is chosen to estimate the maximum joint compression angles along with the changing jumping rhythm. In Figure 3, red line indicates the median of the max compression angle, blue box expresses the interquartile range, and the red + is the extreme outliers. The advantage of using box-whisker plot is that it can minimize the effect of extreme experimental data on the statistical analysis, express data dispersion degree, and be used for comparing examples. From the results it is obvious that the maximum compression angles of all the three joints decrease with jumping cadence increasing, which indicates that the changing velocities have greatly influenced the hopping gait. Therefore, we focus on how to generate the self-adapted jumping gait with variable velocity for biomechanical leg control.

**(a) Maximum hip compression angles**

**(b) Maximum knee compression angles**

**(c) Maximum ankle compression angles**

#### 3. Muscle Synergies Extraction and Analysis

A prominent hypothesis suggests that the biological system underlies muscle contractions during movement execution in a modular fashion by the CNS. These modular have been observed in forms of muscle synergies. Evidences observed in many cases and species show that regularities of the muscle synergies appear to be very similar across subjects and motor tasks [22, 23]. This section will discuss how to extract muscle synergies and investigate if the same group of synergies can explain the difference phases of vertical jumping motion.

##### 3.1. sEMG Signal Preprocessing and Analysis

Many noise signals will contaminate the original sEMG signals during the experiment, such as the inherent noise of equipment like DC bias, motion artifact, or firing rate of the motor units. To remove these unwanted noises is very important to obtain reliable motion intention of human leg. Specifically in this paper, the sEMG preprocessing experienced the following digital operations:(i)After applying low-pass filter with 35 Hz to the raw sEMG signals, DC offset was removed for each channel.(ii)The natural sEMG signal time series vibrates very frequently at the zero point. Through full-wave rectification, the amplitude of signals was more clearly presented.(iii)In this experiment, the sampling frequency of sEMG signal is 1000 Hz, while the motion capture system samples are at a frequency of 500 Hz; thus, the sEMG signals were subsampled to be consistent with motion signals.(iv)A high-pass filter was constructed for a cut-off frequency of 10 Hz.

After the above procedures, the noises of the raw sEMG signals are removed and their envelopes are smoothed. The reprocessing results of one channel sEMG after each digital operation are shown in Figure 4. The envelopes of sEMG signals after preprocessing for seven channels are listed in Figure 5, which can be used for joint motion prediction.

##### 3.2. Extracting Muscle Synergies from Different Phases of Jumping

###### 3.2.1. Muscle Synergies Extraction by NMF

Several feature projection techniques, such as Nonnegative Matrix Factorization (NMF) [24], Principal Component Analysis (PCA) [25], Independent Component Analysis (ICA) [26], Linear Discriminant Analysis (LDA) [27], and nonlinear projections [28], can be used to extract muscle coordination pattern. NMF is commonly used as descriptive measure of specific time-invariant muscle synergies because of relaxed constraints on orthogonality and statistical independence between each component and relative robustness to noisy data [24]. NMF is applied to extract synergies of multivariate sEMG data in the following manner. If the structure of generators, which are combined to generate command sEMG signal across tasks, is chosen as the form of spatial dimensionality [29], for generators, where are the set of signals () for task condition , is condition dependent, time-varying combination coefficient for the th generator, and is the condition-independent, time-invariant th spatial generator. Fundamentally, a muscle synergy consists of a time-invariant weighing coefficient and a time-varying activation coefficient . The weighing coefficients within a synergy determine the number of muscles along with the extent of their activation, while the activation coefficient captures when the muscles are active during a task.

###### 3.2.2. Muscle Synergies for Stance and Flight Phase

All the experimental data of each subject were grouped according to the jumping rhythm. The data were divided into five sections with different rhythm . If the number of muscle coordination patterns is chosen as four, the extracted muscle synergies of stance phase for different velocities are given in Figure 6 (taking jumping rhythm of 430 ms, 480 ms, 530 ms, and 580 as examples). The synergies weighting coefficients indicate that the four coactivation patterns are vastus medialis, rectus femoris, and gluteus maximus; gastrocnemius, soleus; tibialis anterior; and biceps femoris, respectively. To evaluate the quality of the extracted synergies, the variance accounted for (VAF) is usually used to calculate the percentage of variability in the sEMG dataset that is accounted for by the extracted synergies. The average VAF value of the six subjects who took part in the experiment is calculated. VAF for leg stance phase with different jumping cadence are listed in Table 1, in which all VAF values are bigger than 90%. This indicates that synergies from one velocity can well explain other velocities and the recorded sEMGs with different jumping velocities are well reconstructed by the same extracted synergies.

**(a) Jumping rhythm of 430 ms**

**(b) Jumping rhythm of 480 ms**

**(c) Jumping rhythm of 530 ms**

**(d) Jumping rhythm of 530 ms**

Next, we will test if the extracted synergies of stance phase can explain the muscle coordination pattern of flight phase. Supposing that the number of synergies is four, the VAF values for different hopping velocities are listed in Table 2. The low VAF values demonstrate that four muscle synergies are not sufficient to explain the flight phase with different velocities and it may employ more muscle coordination patterns. This conclusion is also approved by the calculation results in Figure 7, in which the synergies weighting coefficients are not the same for different velocities. On one hand, during the flight (swing) phase, the number of synergies is more than four and the muscles are more freely organized based on the above results. One the other hand, the effect of changing joint angles is to maintain the whole body in balance state since there is no contact force. Therefore, it is inappropriate to generate the hopping gait of flight phase with the same synergies of the stance phase. In the following section, the muscle synergies are employed to generate hopping gait only for stance phase, while the joint angles trajectory for swing phase is generated by interpolation of the jump-initiation state and touching-the-floor state.

**(a) Jumping rhythm of 430 ms**

**(b) Jumping rhythm of 480 ms**

**(c) Jumping rhythm of 530 ms**

**(d) Jumping rhythm of 580 ms**

#### 4. Generate Velocity-Adapted Flexible Jumping Gait

##### 4.1. Estimating the Reference Gait Pattern for Stance Phase

The wavelet neural network (WNN) is adopted to identify the relationship between the muscles coactivation patterns and coordinated variation of joint angles. The wavelet and the neural network processing can be performed separately. Firstly, the input signals are decomposed using some wavelet basis stored in the hidden layer. The hidden layer consists of neurons, which are usually referred to as wavelons, and their activation functions are drawn from a wavelet basis. Secondly, the wavelet coefficients are output to some summers, whose input weights are updated in accordance with certain learning algorithm.

For a multi-input multioutput nonlinear identification system, the output can be defined by wavelons aswhere and are the dilation and translation parameters, respectively. The parameters , , , and can be grouped into a parameter vector . The objective function to be minimized is defined asThe minimization of the above function is performed using stochastic gradient algorithm. This recursively modifies . After each sample pair , the objective function is written asThe parameters and can be fixed at initialization of the network. , which is the only parameter that needed to be adjusted, is modified in the opposite direction of , and the gradient can be computed by the partial derivatives as follows:in which

As it is analyzed in Section 3, the number of muscle synergies which can explain the muscle coactivation patterns of stance phase is four. A 4-input and 3-output wavelet neural network is constructed to predict the coordinating variation of three joint angles with muscle synergies for stance phase, taking the muscle activation coefficient as input signal.

##### 4.2. Estimating the Reference Gait Pattern for Flight Phase

Considering the flight phase of jumping motion, the number of muscle coordination patterns is more than four and the muscles are more freely organized since there is no contact force with the ground (as analyzed in Section 2.2). To estimate the covariation of joint angles with extracted synergies is a rough task and may cause computational burden. Therefore, we use cubic polynomial interpolation to generate the reference joint angles of flight phase, taking the conditions of the adjacent two jumps. Assume the time of each jump accounting for the moment of the feet touching the ground until the beginning of the next stance phase, which includes stance time and flight time . The joint angle trajectories of flight phase for the jump during continuous jumping can be expressed asin which

All the hip, knee, and ankle joint angles are uniformly expressed in vector , and is the coefficient of cubic polynomial expression. Take the end state of stance phase for the jump as the initial condition of the flight phase and the start state of stance phase for the jump as the terminating condition :in which and are the initial joint angle and angle velocity of flight phase for the jump, while and are the initial joint angle and angle velocity, respectively. Substituting (10) into (7), coefficient of cubic polynomial expression can be solved.

##### 4.3. Gait Generalization with Takagi-Sugeno Fuzzy Inference System

The proposed fuzzy inference system is based on the well know Takagi-Sugeno fuzzy inference system (TS-FIS). The TS-FIS is described by a set of fuzzy rules presented in (11). are the inputs of the FIS with dimension input space, and are linguistic terms, which are numerically defined by membership functions distributed in the universe of discourse for each input . Each output rule is a linear combination of input variables.

In our problem, each reference gait has been memorized by one wavelet neural network system, which records the relationship between muscle coactivation coefficient and joint coordinated variation. Based on the available hopping velocity from experiment data, totally six reference gait patterns are chosen and memorized.

The membership functions are shown in Figure 8, in which the desired jumping rhythm is modeled by six fuzzy sets (“VeryFast,” “Fast,” “MediumFast,” “MediumSlow,” “Slow,” and “VerySlow”). The final desired trajectory of hip, knee, and ankle joint angles is generalized by the TS-FIS and WNN architecture, on the basis of predefined membership functions: If is VeryFast then . If is Fast then . If is MediumFast then . If is MediumSlow then . If is Slow then . If is VerySlow then . The output is carried out in two stages [30].

Firstly, take the extracted coactivation coefficient as the inputs of the wavelets neural network. The coordinating variation of hip, knee, and ankle joint angles can be identified according to (2), indicating each reference gait pattern with corresponding muscle coactivation coefficient.

Secondly, the generalized joint angle is estimated using the weighted average of all wavelets neural networks, which is calculated usingwith given byand is computed with the membership function according to

By using the fuzzification of several outputs of wavelets neural network, it is possible to obtain a global generalization, which allows decreasing the memory size and computing cost using only a small set of identification system. Moreover, it makes the generalization of gait pattern for a variety of velocities possible, with only limited sets of jumping experiment results.

#### 5. Discussion

##### 5.1. Muscle Synergies for Stance Phase of Vertical Jumping

If the number of synergies is chosen as four, the simulation results in Table 1 (Section 3.2) show that the VAF of all the trials is above 90%, which indicates that the recorded sEMGs are well reconstructed by the extracted synergies. The number of synergies that can explain the muscle coactivation pattern efficiently was discussed in our previous work in [31]. Take two experimental trials as examples (uniform velocity jumping and variable velocity jumping). The sEMG signals after preprocessing are shown with the same time sequence in Figure 9. The four synergies of stance phase are gastrocnemius, soleus; vastus medialis, rectus femoris, and gluteus maximus; tibialis anterior; and biceps femoris, respectively. The muscles were grouped into one synergy activated with very similar time sequence, while the time series of muscle explosive for each synergy are different from one another. It is noticed that the biceps femoris collected by the first channel sEMG is always activated in front of other muscles and even before the moment of maximum joint compression. This observation coordinates with the conclusion in [32], in which, before the stance phase of jumping, some muscles have already been preactivated in order to decrease the leg impact force with ground.

**(a) Continuous jumping with uniform velocity**

**(b) Continuous jumping with variable velocity**

##### 5.2. Velocity-Adapted Gait Generalization

###### 5.2.1. Flexible Vertical Jumping Gait for Stance Phase

The hip, knee, and ankle joint angle defined in bionic leg model are shown in Figure 10. As designed in Section 4.1, a 4-input and 3-output WNN was constructed to predict the coordinating variation of joint angles with muscle synergies for stance phase. Taking the muscle activation coefficients as input signals, the output joint angles of neural network for reference “Gait3” defined in Section 4.3 (jumping with rhythm of 517 ms) are shown in Figure 11, in which blue line expresses the joint angles from motion capture data and red line is the estimated outputs. In order to present the performance of the proposed method more clearly, all the joints were undergoing normalization by the maximum value.

**(a) Hip angle**

**(b) Knee angle**

**(c) Ankle angle**

Root mean square error (RMSE) is computed to measure the performance of the wavelet neural network model. The RMSE is defined as follows:in which is the estimated joint angle, is the real joint angel, and is the length of the sample data. Table 3 gives the RMSEs between the real joint angles from the experimental data and the estimated ones for one subject, which does not have much difference from other five subjects. The results indicate that the proposed WNN model with muscle synergies shows good performance both in accuracy and in robustness for different hopping velocities. Six reference gait patterns of stance phase are memorized in these neural networks.

The fuzzy inference system proposed in Section 4.3 is used to generalize more flexible gait with limited experimental gait pattern. Six triangle membership functions were set according to the different hopping velocity grouped based on the experimental data. Figure 12 presents the generalized gait for stance phase with jumping rhythm of 580 ms after FIS, which refers to the reference gait pattern with cadence of 560 ms and 601 ms, and fuzzy sets of “MediumSlow” and “Slow.” The hip angle of the reference gait pattern is especially shown in Figure 12(a). The difference between the real hip angle and the generalized one after FIS is bigger than the other joint angles. Furthermore, we analyze the stance phase gaits of the same hopping velocity obtained from motion capture data. It can be seen from Figure 13 that the hip angles vary more than other joint angles even for the same subject. This may explain why the generalized hip gait is slightly different from the real experimental data.

**(a) Hip angle**

**(b) Knee angle**

**(c) Ankle angle**

**(a) Hip angle**

**(b) Knee angle**

**(c) Ankle angle**

###### 5.2.2. Velocity-Adapted Continuous Vertical Jumping Gait for Bionic Leg

Next, we evaluate the proposed approach to generate self-adapted continuous vertical jumping for bionic leg with variable velocity. Figure 14 represents the motion sequences and joint angle trajectories with the jumping velocity changing from slow to fast. While the adjacent two jumps from fast to slow jumps are shown in Figure 15. From the results, it should be noticed that, in the case of slow jump of bionic leg, all of the hip, knee, and ankle angles are more compressed than fast jump. This also follows the regularities observed from the subject jump experiment. The generalized gait is compliant in the condition of changing velocity, either from slow to fast hops or from fast to slow hops.

**(a) Motion sequences**

**(b) Joint angle trajectories**

**(a) Motion sequences**

**(b) Joint angle trajectories**

Figure 16 shows generalized joint angle trajectories of bionic leg for velocity-adapted continuous jump. The results indicate that the proposed approach can generate flexible hopping gait very similar to human motion with limited experimental data. Furthermore, compared with stance phase, the joint angles in the flight phase show smaller fluctuation. That is because the jumping gait of flight phase is generalized with polynomial interpolation simply. With this proposed approach, the self-development of new flexible vertical jumping gait can be achieved only referring to limited reference gait patterns, which extremely shortens the training time of estimation model for sEMG-based biomechanical leg control.

**(a) Hip angle trajectory for continuous jumps**

**(b) Knee angle trajectory for continuous jumps**

**(c) Ankle angle trajectory for continuous jumps**

#### 6. Conclusion

Inspired by the hypothesis that CNS modulates muscle synergies to simplify the motor control and learning of coordinating variation of redundant joints, this paper proposed a novel approach for flexible gait generation of hopping motion with sEMG signals. Two questions were analyzed and discussed in the paper, the first one concerning whether the same set of muscle synergies can explain the different phases of jumping movement with various velocities. Based on the analysis of synergies weighing coefficients and variance accounted for (VAF) value, the muscle synergies were extracted for stance phase separately. The second one is about building a model for generating velocity-adapted jumping gait with muscle synergies, in which wavelets neural network is proposed to predict the reference gait pattern, while fuzzy inference system is adopted to merge these reference gaits in order to create more generalized gaits with different jumping rhythm. From the experimental results, the proposed method shows good performance both in accuracy and in robustness for producing continuous flexible jumping gait with different velocities.

The proposed method can be adopted as the decoder in sEMG-based controls [21] for bionic leg, which decodes muscle activities into intuitive control outputs by training a model on sEMG-related inputs with desired motion gait. Once the decoder was trained, it is used in real time to estimate the multijoint angles and map them to output of myoelectric interface. Moreover, linear combinations of synergies are capable of describing complex force and motion patterns in reduced dimensions [33]. Therefore the robust representations of synergies within the control scheme can generate flexible gaits for other complex motions, such as hopping and running with user’s intent.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research has been made possible by “111 Project” (Grant no. B13044), Natural Science Basic Research Plan of Shaanxi Province (Grant no. 2016JQ6009), and National Science Foundation of China (Grant no. 61603302).