In a multiple degrees of freedom motion system with rigid-flexible coupling, the flexible internal dynamics have a significant negative impact on performance because the degrees of freedom are coupled. Inspired by this problem, a multi-input multioutput state-space model based on modal coordinates is proposed to decouple the rigid body and flexible modes. The closed-loop subspace identification method based on orthogonal is utilized to develop an unbiased standard state-space model. Based on the similarity principle, a modal analysis method is proposed to transfer the standard state-space model to the proposed modal-decomposition-dependent state-space model. Controllability and observability criteria are met to guarantee the minimum realization of the proposed state-space model. Finally, a modeling and modal analysis experiment is conducted on a developed wafer stage of lithographic tool. The results verify the effectiveness and accuracy of the proposed modeling method, as well as the controllability and observability of the proposed model.

1. Introduction

A motion system is the motion mechanism in a servo mechanism, and motion control technology is used to achieve precise positioning or motion control. It is widely used in robotics [1, 2], computer numerical control machine tools [3], and other fields [46], and it is the core subsystem in precision and ultraprecision high-end equipment [7]. The requirements for motion systems include high resolution, high precision, stability, fast response times, high bandwidth, and robustness. However, vibration is a challenge that must be overcome to meet these requirements. If only the rigid dynamics of the system are considered, then the flexible dynamics are not addressed. Therefore, the system should be treated as a rigid-flexible coupled system [1, 5, 6]. In high-end equipment, the motion system is typically a multiple-input multipleoutput (MIMO) system with multiple degrees of freedom (DOF), and the flexible internal dynamics will produce coupling effects between different motion axes [8, 9]. Therefore, it is necessary to study the modeling of rigid-flexible, coupled, multi-DOF motion systems.

The behavior of this type of system is governed by second-order ordinary differential equations, which can be obtained using the finite element method [1012]. Since the goal of this paper is to develop a model for control, the traditional way is to build an input-output (IO) frequency-domain model, which is in the form of a transfer function [13, 14]. IO model is employed to represent the system’s input-output relationships. For MIMO systems, the MIMO relationship is expressed as a generalized transfer function matrix model [7, 9, 15]. The elements of the transfer function matrix are independent transfer functions, but for a flexible MIMO system, the elements have the same denominator polynomial [16, 17]. It is helpful to reduce invalid number of parameters of model by using the same denominator polynomial. However, the above two models are frequency-domain IO models, which make it difficult to describe the internal dynamics of the system. If an IO model is utilized to control the flexible modes, the flexible mode transfer functions need to be abstracted [15]. It requires increasing the number of IOs, which is difficult to achieve or extremely costly for electromechanical systems with fixed drivers and measuring sensors [18]. Therefore, we aim to adopt the model form without increasing the number of IOs.

In this paper, the closed-loop subspace identification method [19, 20] is adopted to obtain a standard state-space model without deviation through orthogonal and oblique projection. Based on the similarity principle, the proposed modal-decomposition-dependent state-space model is obtained. Since the realization of the state-space model describing the same dynamics of IO is not unique and can be infinite, the minimum realization is needed, which means that the states can be controlled and observable. Thus, the controllability and observability criteria of the proposed state-space model are deduced, which are also a useful foundation for control. Finally, the proposed modeling and closed-loop identification methods are applied to developed lithography tools wafer stage system. The experimental results demonstrate that the proposed modeling and modal analysis method can directly obtain a modal-decomposition-dependent state-space model with rigid body modes and flexible modes decoupling. Furthermore, the obtained model can be controlled and observed, and the modeling accuracy is high, indicating this method is feasible and effective.

The main contributions of this paper have two aspects: (1) an MIMO modal-decomposition-dependent state-space model is proposed to describe the internal dynamics of a rigid-flexible, coupled, multi-DOF motion system with an inherent number of IOs; (2) the corresponding controllability and observability criteria are proposed and proved, which are important for active vibration control.

The following sections of this paper are organized as follows. In Section 2, the theoretical modeling for rigid-flexible, coupled, multi-DOF motion systems is presented, while the modal-decomposition-dependent state-space model and its identification method are proposed. Then, the modal analysis method is proposed in Section 3, and the controllability and observability criteria are introduced in Section 4. Finally, the experiment carried out is described in Section 5, and Section 6 concludes this paper.

2. Modeling for Rigid-Flexible, Coupled, Multifreedom Motion System

A rigid-flexible, coupled, multi-DOF motion system has multiple-DOF moving parts with a flexible structure. In each DOF, the dynamics are rigid-flexible and coupled. The wafer stage of the lithographic tool is a typical rigid-flexible, coupled, multi-DOF motion system, as shown in Figure 1.

2.1. Theoretical Modeling

The system can be accurately described as a distributed parameter model [21] having distributed mass, damping, and stiffness. Because motion states should be described in space-time coordinates, the system’s equation is in the form of a partial differential equation, as follows:where and are space and time coordinates, respectively. is the distributed mass parameter, while is the displacement, is the linear conjugate differential operator, and is the distributed force.

However, a partial differential equation, such as (1), has analytical solutions in only a few simple and special cases. As the actual conditions of a plant are always complicated, the analytical solution of a partial differential equation cannot be obtained directly in engineering practices. Therefore, the actual infinite-dimensional, distributed parameter model is often approximated as a finite-dimensional, multi-DOF, discrete-system model [22].

Second-order ordinary differential equations are often used for the dynamic modeling and analysis of multi-DOF systems, which are comprised of separation mass, damping, and stiffness matrices. An -DOFs discrete motion system with inputs and outputs can be represented as follows:where is the vector with -DOFs, is the input vector, is the output vector, is the matrix of input, is the matrix of displacement variable, is the matrix of velocity variable, is the positive-definite mass matrix, is the positive semidefinite damping matrix, and is the positive semidefinite stiffness matrix.

Although (1) and (2) have different mathematical expressions, they have the same dynamic properties. A continuous system can be regarded as a special case where the dimension of the discrete system reaches infinity, while a discrete system is the approximate representation of the finite dimension of a continuous system. Equation (2) realizes the discrete approximation of the continuous model as (1). The vector is based on physical coordinates, which can be utilized for multifreedom motion control on physical coordinates. However, there is a coupling relationship between the physical coordinates of each DOF in the multifreedom motion system. The modal coordinates are mutually orthogonal and perfect for the uncoupled decomposition. Thus, the paper proposes a modal-decomposition-dependent state-space modeling and modal analysis method, which is well-suited for rigid-flexible, coupled, multi-DOF motion system control.

2.2. Modal-Decomposition-Dependent State-Space Modeling

As stated in Section 2.1, modal coordinates are often used in the dynamic analysis and control of flexible modes. In a modal coordinate system, all the modal coordinates are orthogonal; thus, there are no coupling relationships between modal coordinates. Transforming the physical and modal coordinates is performed as follows:where is the vector of modal coordinates. is the transformation matrix, which includes eigenvectors as , and , is the -th normalized eigenvector, which is , where is the -th eigenvector.

From (2) and (3), an original dynamical model of modal coordinates can be generated:where and are the diagonal damping matrix and diagonal stiffness matrix, respectively, and and are the -th modal damping factor and eigenfrequency.

The diagonal matrices realize the internal decoupling between the various coordinates. Under the decoupling of the modal coordinates, the -th modal can be written aswhere denotes the modal coordinate position of -th mode and is the eigenfrequency. If the for the rigid mode, is the damping ratio. If the is small for the light damping flexible mode, is the input row vector, while is the real input vector, and is a scalar as equivalent input.

Equation (5) achieves decoupling between modal coordinates, and modal control can be realized based on this model. However, the second-order differential equations are not the standard model form for modern control theory. Thus, (5) can be rewritten as the form of state space, and the state-space model can be normalized to the Jordan standard type. It also can be regarded as the states of the rigid or flexible modes and is defined as . The corresponding -th state subspace iswhere is the output matrix. It can also be considered the controllable canonical form. It is easy to judge whether a system is controllable if , according to the controllability judgement.

From the mode superposition method, a nominal modal state-space model, which is composed of the rigid and flexible modes, is expressed aswhere , and are the state matrix, input matrix, output matrix, and feedthrough matrix, respectively, where

Equation (7) is an MIMO state-space model that achieves decoupling between rigid modes and flexible modes. It can also express flexible internal dynamics by the states of flexible modes, which is more useful in rigid-flexible, coupled, multi-DOF motion system research than the IO model of DOF.

2.3. Identification of an MIMO State-Space Model

Having described the theoretical modeling approach in the previous subsection, the focus changes to solving the problem of obtaining an MIMO state-space model using an experimental approach. Experimental modeling can be divided into frequency- and time-domain methods, according to different definition domains.

2.3.1. Frequency-Domain Identification

Frequency-domain identification is a nonparameter identification method, which usually refers to the spectral estimation method of frequency-domain functions [13]. The difficulty of algorithm implementation when using the frequency-domain identification method is not affected by the complexity of the model, so it is often used in dynamic models of motion systems. However, for multi-DOF MIMO systems, frequency-domain identification results in a single-input multioutput (SIMO) model, which needs multiple identifications prior to integrating it into an MIMO model [5].

A frequency-domain identification method is presented to identify the MIMO motion system with a flexible structure in our previous study [14]. In this method, a nonparametric model is obtained via a frequency-domain identification experiment, and an MIMO transfer function matrix with the same denominator is obtained through general parameters orthogonal polynomial curve fitting [16]. Finally, the state-space model of the MIMO system is obtained through the modal superposition and singular value decomposition (SVD) principle.

2.3.2. Time-Domain Identification

Subspace identification is one of the time-domain identification methods [23]. The advantages of subspace identification are as follows: no prior knowledge of the model is required; no initial conditions and initial states of the system are required; no explicit model parameters are required for the identified model. It uses an orthogonal projection or oblique projection algorithm and has better computational efficiency and numerical robustness. The order of the model is determined through SVD. In the least-squares method, the identification process does not require iterative calculations, and the state-space model can be identified directly.

Subspace identification is applicable to the identification of MIMO systems. The open-loop subspace identification process involves, first, obtaining input and output data through identification experiments. Then, Hankel matrices are constructed. Through orthogonal or oblique projection, the Kalman state sequence of the model and extended observable matrix are obtained. Finally, the least-squares algorithm is used to calculate the system matrices.

There are stability and safety problems in the open-loop identification of many precision motion systems, such as lithography machine tools. Therefore, the identification must be conducted in the closed-loop condition. In terms of closed-loop identification methods, the two-step method and projection method [19, 20] are adopted. The specific calculations process is shown in Appendix A. Through this method, a nominal state-space model is obtained, and it can ensure that the identification result is unbiased.

Compared with the frequency-domain identification method, which first identifies the SIMO model and then integrates the MIMO model, subspace identification can improve identification efficiency by simultaneously stimulating multiple inputs to obtain the system nominal state-space model ast

However, the identified result is the standard state-space model (9), in which states have no physical significance.

In this subsection, a modal analysis method using the similarity principle is presented to obtain the proposed modal state-space continuous domain model (7) with state variables in the modal coordinate system.

The subspace-identified model is a discrete domain model. Thus, the corresponding continuous state-space model is obtained by digital-to-analog conversion:where , , , and are the system matrices of the continuous domain state-space model. In this paper, a modal analysis method based on the similarity principle is proposed to obtain the standard state-space model in the form of (7). The model is diagonalized based on the similarity principle:where is the state transition matrix composed of the system eigenvectors, which is a diagonal matrix. Its diagonal elements are the eigenvalues of the system. Corresponding to the conjugate eigenvalues, the subspace is obtained aswhere and are a pair of conjugate poles of the corresponding system transfer function. The corresponding differential equation is in the form of

Equations (5) and (13) in the modal coordinate system are identical in form. However, (6) and (12) have different definitions of the state variables. Because the definition of system states is nonuniqueness, (12) is a diagonal characteristic value standard form. It is defined that the system states of the modal coordinate subspace system as require further standardization to a controllability type specification:

Through the superposition of modal subspace models, the whole system state-space model can be obtained as

Equation (15) is the final whole system state-space model, which is identical with (7) of the modal state-space model. Through the mode analysis method proposed in this section, the canonical modal state-space model is obtained by a closed-loop subspace-identified state-space model.

4. Controllability and Observability Criteria

A state-space model is proposed as in (7) to describe the internal dynamics of the system under the condition of the inherent quantity of inputs and outputs. Since the realization of the state-space model is not unique and can be infinite, describing the same dynamics of system in the same inputs and outputs condition, the minimum realization is needed to be introduced, which describes the dynamic characteristics of system in same inputs and outputs by using the minimum number of states. The system is the minimum realization if and only if the system states are controllable and observable.

6The state-space model must ensure controllability and observability. If it does not, the states are sufficient to fully describe the dynamics of the system. The existence of states that cannot be controlled and observed results in the following problems: (a) The unknown parameters of the system cannot be identified through input and output data. (b) In the state feedback control, the states cannot be observed by a state observer, so it is difficult to form the state feedback control loop. (c) The noncontrollable state cannot be controlled. (d) If this state is unstable mode variable, the system cannot achieve system stabilization through feedback control. Therefore, it is necessary to discuss the controllability and observability of the model obtained by theoretical and experimental modeling and describe the controllability and observability criteria of the model.

4.1. Theoretical Modeling, Controllability, and Observability Analysis

The proposed modal-decomposition-dependent model (7) can be written as follows:where is the states of state-space model, where are the six rigid mode DOFs’ states and is the -th flexible mode states.

4.1.1. Criterion of Controllability

Lemma 1. Criterion of controllability: the sufficient and necessary condition of controllability of the proposed modal-decomposition-dependent state-space model is that are nonzero vectors and linearly independent of each other.

Proof. Proof of Lemma 1 is detailed in Appendix B.

4.1.2. Criterion of Observability

Lemma 2. Criterion of observability: the sufficient and necessary condition of observability of proposed modal-decomposition-dependent state-space model is that the columns ofare not all zero vectors.

Proof. Proof of Lemma 2 is detailed in Appendix C.

4.2. Experimental Modeling, Controllability, and Observability Analysis

In this paper, experimental modeling and closed-loop state-subspace identification, combined with modal analysis, are used to obtain the state-space model. For subspace identification, the identification result must be both controllable and observable, so it must be the minimum realization [23].

In the process of modal analysis, the similarity principle is adopted, and the transformation matrix is a matrix composed of eigenvectors. The conjugate subspace is also transformed into a controllable canonical type. Since it is the controllable canonical type transformation of MIMO, the canonical type is not unique, and different change matrices can obtain the different canonical forms, such as the frequently used Wonham and Luenberger types. In this paper, the controllable specification type is directly calculated by (12)–(14). All of these processes involve nonsingular transformation, which means the transformation matrix is nonsingular. According to the properties of nonsingular transformation, which does not change the controllability and observability of the system, the process of modal analysis also does not change the controllability and observability of the model.

Based on the experimental modeling of closed-loop subspace identification and modal analysis methods described in this paper, the obtained model must be controllable and observable. Furthermore, the MIMO modal state-space model is the minimum realization, meaning the state-space model with the least number of parameters. The model is suitable for designing a state feedback controller, state observer, determination of system stability, and so on.

Based on the above comprehensive analysis of the theoretical state-space model, the controllability and controllability criteria of the model are proposed. It is also explained that the experimental modeling method adopted in this paper does not change the controllability and observability of the model.

5. Experiment

The wafer stage of a lithographic tool realizes multidegree ultraprecision and high speed/acceleration motion simultaneously. Due to the ultraprecision and high response bandwidth requirements, the flexible dynamics of the structure of the wafer stage cannot be ignored. Therefore, the motion system of a wafer stage is a typical ultraprecision, rigid-flexible, coupled, multi-DOF motion system, which can be used to verify the effectiveness and feasibility of the proposed method.

5.1. Experiment Setup

The proposed modeling and modal analysis method in this paper is applied to the experimental wafer stage system, as shown in Figure 2. The measuring frame is made of marble, which provides high stiffness. Meanwhile, the high-performance vibration isolator can isolate external low-frequency vibration. Their common function is to prevent the vibration in the measuring system and, thus, meet the requirement of nanoscale precision measurement. The measuring system is a nine-axis laser interferometer, which is used to measure the displacement of the wafer stage with 6 DOFs relative to the frame, and the measurement resolution can reach 0.15 nm. The motion system consists of two parts: a long-stroke (LS) moving stage and a short-stroke (SS) moving stage, which realize long- and short-stroke motion, respectively. The SS stage has 6 DOFs and is driven by voice coil motors. All signal processing and system identification are realized based on a discrete domain time system with a sampling frequency of 5 kHz. The experimental wafer stage system is detailed in [24, 25].

5.2. Implementation of Modeling and Identification

In the implementation of modeling and identification, the closed-loop identification approach is employed, as shown in Figure 3. Particularly, the input excitation signals are added to the input signals . The input signals are 6 DOFs forces, which will be transformed into 8 motor forces in the wafer stage, and the output measurement signals are the 6 DOFs displacements decoupled by nine-axis laser interferometer signals, and the sampling data of are 10,000 in this identification testing.

A 16-order state-space model was obtained through the identification method shown in Appendix A. Combined with modal analysis proposed in Section 3, the corresponding model’s 6 DOFs input/output frequency-domain responses are shown in Figure 4, including 6 rigid body modes and 2 flexible modes as (7) form, where

It is worthy of note that an infinite number modes model is neither possible nor necessary in engineering practice, and the model is realized as the model reduction by singular value decomposition (SVD) in (A.6). The singular value determines the influence of the modes on the dynamic system. The flexible dynamics can be approximated by selecting the first few orders according to the actual situation.

5.3. Modal Analysis

Modal decomposition is realized by the above model. Using the Y direction (scanning direction for wafer stage) DOF as an example, the IO dynamic relationship can be decomposed into one rigid mode and two flexible modes, which, in turn, are the mode superposition used to obtain the whole IO model of the DOF, as shown in Figure 5.

The rigid body mode in the y DOF is a double integral form. The flexible modes have resonance frequencies of 913.8 Hz and 1032 Hz. According to the theoretical state-space model (7), the IO model of physical coordinates, conventionally expressed as the transfer function form, is decomposed into the superposition form of three modes, which can realize mode decoupling.

The state equation of the y modal coordinate is the y DOF rigid mode under the modal coordinate:

The state equations of two flexible modes under the modal coordinate are

The y physical coordinate is obtained by modes superposition:

Thus, the orthogonal mode states under the modal coordinate can be obtained from state equations of the proposed state-space model, and the 6 DOFs physical coordinate can be obtained by the modes’ superposition.

5.4. Controllability and Observability Analysis

According to the 4 section’s analysis, the final state-space model, obtained from the proposed experimental modeling method, must be the minimum realization (i.e., controllable and observable).

The final state-space model is shown in Section 5.2, in which controllability and observability are determined by the criteria presented in Sections 4.1.1 and 4.1.2:

Thus, are nonzero vectors and linearly independent from each other. This satisfies the criterion of controllability presented in Section 4.1.1:

Thus, the columns of are not all zero vectors, and are not all zero vectors. This result satisfies the criterion of observability presented in Section 4.1.2.

Therefore, the model is completely controllable and observable, as it meets the controllability and observability criteria presented in Section 4.

5.5. Model Accuracy Analysis

From the above statement, it can be seen that the proposed MIMO modal state-space model is obtained effectively from the closed-loop subspace identification and modal analysis method. The model is completely controllable and observable according to the controllability and observability criteria.

According to the modal model parameters, the natural frequency and damping of the flexible modes are shown in Table 1.

It is an effective way to improve the response speed and reduce steady-state errors by increasing the closed-loop bandwidth, but the existence of low-order resonance limits making improvements in the closed-loop bandwidth. By using this characteristic, the feedback controller is adjusted to increase the closed-loop bandwidth so that the wafer stage resonates, and the Z-direction static error signal is obtained, as shown in Figure 6. The resonance frequencies of the wafer stage motion system are 913.8 Hz and 1032 Hz through the unilateral amplitude power spectrum analysis of the z static error signal.

Compared with the parameters in Table 1, the identified model parameters, natural and resonant frequencies of the flexible modes, and extended control bandwidth result in system resonance frequency errors of 4.6 Hz and 0.5 Hz, respectively, indicating that the obtained method has high accuracy.

5.6. Experiment Result Discussion

The experiment is conducted on a developed ultraprecision wafer stage using the proposed modeling and modal analysis method, and the results validate the practicability and effectiveness of the proposed modeling and modal analysis method. The controllability and observability of the proposed model can be easily determined by the proposed criteria of controllability and observability. Additionally, the accuracy of the proposed modeling method is verified by the spectral analysis of vibration signals.

6. Conclusion

This paper studied modeling and identification, as well as modal analysis and its controllability and observability, for a rigid-flexible coupled multifreedom motion system. Some conclusions can be drawn, as follows:(1)The modal-decomposition-dependent state-space modeling and modal analysis method proposed in this paper achieves rigid-flexible mode decoupling between modal coordinates, which is adaptive for modern control theory and methods used to solve the control of rigid-flexible coupled multi-DOF motion systems.(2)The proposed criteria for judging the controllability and observability of the proposed modal-decomposition-dependent state-space model are simple and convenient to operate.(3)The experimental results show that the modal-decomposition-dependent state-space model can be obtained using the proposed modeling and modal analysis method. It has high accuracy and is controllable and observable based on the proposed criteria for controllability and observability. The feasibility and effectiveness of the proposed method are also verified by the experimental results.

The goal of this paper is to propose a modal-decomposition-dependent state-space modeling and modal analysis method for developing a rigid-flexible, coupled, multi-DOF motion system for modal control that can directly eliminate arbitrary vibration modes and improve the control performance of the system. Based on this paper’s model, the work of [25, 26] is studied to achieve the above goals.


A. Closed-Loop Subspace Identification

The calculation process of closed-loop subspace identification is stated as follows.

First, we provide some basic definitions. The exogenous inputs and input-output data are given by the identification experiment. The exogenous input block Hankel matrix is defined asand the input and output data block Hankel matrices and can be defined in the same way.

The closed-loop subspace identification uses a two-stage method and a projection method [19, 20], which could realize an unbiased identification result since there is no correlation between the unmeasured noise and exogenous inputs . The specific procedure is realized as follows:Step 1. Construct the matrices and aswhere is shorthand for the matrix projecting into the matrix and denotes the Moore-Penrose pseudoinverse of the matrix .Step 2. Identify the open-loop system using the model:where and are the new input and output vectors from and and is the new state variable vector using the same process as and .

Through (A.3) and (A.4), an open-loop deterministic identification problem is formulated. The next algorithm is identical to the deterministic open-loop system subspace identification problem. The solution to this problem utilizes a projection algorithm [19, 20]:

Calculate the singular value decomposition (SVD) of the weighted projectionwhere and are the weight matrices, which are chosen in Table 2 by different subspace algorithms, such as numerical algorithms for subspace state-space system identification (N4SID), MOESPs, and canonical variate algorithms (CVAs) [23]. Determine the order by inspecting the singular values in and partition the SVD accordingly to obtain and .

Determine the extended observability matrix as

Determine the extended observability matrix as

Determination of and is as follows:

can be calculated from as , where denotes without the last rows and denotes the without the first rows. Next, can be calculated as the first row of .

Determination of and is as follows:

From the input-output equation in [23], we can obtainwhere the lower block triangular Toeplitz matrix is defined as

Then, a least-squares method to solve and can be derived:

B. Proof of Lemma 1

Proof. Controllability eigenvalue criterion: the sufficient and necessary condition for an -dimensional linear time-invariant system to be fully controllable is that all the eigenvalues of , , are satisfied:According to the controllability eigenvalue criterion, for the rigid body modes, the eigenvalue of the rigid body mode is :The controllability of rigid body modes satisfying the above equation is nonzero vector and linearly independent.
The controllability condition of the flexible mode can be determined by the eigenvalue criterion:It is a pair of conjugate poles :By elementary transformation,The controllability of the flexible mode satisfying the above equation is a nonzero vector and linearly independent.
In conclusion, the sufficient and necessary conditions of controllability of the proposed modal state-space model are that is nonzero vector and linearly independent and is not nonzero vector and linearly independent. □

C. Proof of Lemma 2

Proof. Observable eigenvalue criterion: an -dimensional linear time-invariant system . The sufficient and necessary condition for complete observability is that all the eigenvalues of , , , satisfyAccording to the controllability eigenvalue criterion, the observability condition of the model can be denoted as follows:Using the conjugate pole relation, such as (B.4), after the elementary transformation,Thus, the sufficient and necessary condition of observability of the proposed modal-decomposition-dependent state-space model is that the columns ofare not all zero vectors. □

Data Availability

The data reported in this paper are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.


This work was supported by the National Natural Science Foundation of China under Grant 51805053, the Fundamental Research Funds for the Central Universities under Grant 2020CDJGFJX010, the National Major Science and Technology Project of China under Grant 2017ZX02102004, and the Natural Science Foundation of Hubei Province of China under Grant 2019CFB206. The authors express their gratitude for the support.