Abstract
The cis-lunar space has been more and more attractive for human beings, and different kinds of missions have been proposed. For cis-lunar missions with long durations, the stationing-keeping is a pivotal problem. In this paper, the station-keeping problem with continuous thrust for different cis-lunar orbits, including distant retrograde orbits (DROs), near rectilinear halo orbits (NRHOs), and halo orbits, are investigated in the ephemeris model. The optimal sliding mode control (OSMC) based on the linear quadrant regulator (LQR) control is designed for the station-keeping problem. Simulations only considering the initial insertion error are conducted first to show performances of the OSMC controller, and the Jupiter gravity and solar radiation pressure (SRP) are then included as unknown perturbations to test the controller’s robustness. Then, with considerations of more practical constraints caused by the navigation and propulsion systems, Monte-Carlo simulations are carried out to provide more realistic results, and station-keeping performances are compared and analyzed for different nominal orbits. The results can provide useful references for the selection of station-keeping strategy in future long-term lunar missions.
1. Introduction
The enduring scientific exploration of the cis-lunar space is the next essential step in the development of space economy. An important measure that can greatly enhance the long-term capacity of scientific exploration is the construction of cis-lunar space stations. The National Aeronautics and Space Administration (NASA) has proposed the Artemis program [1], including the Lunar Orbital Platform-Gateway (LOP-G). The LOP-G is designed to provide valuable support for robotic and manned missions, and to serve as a stepping stone for future exploration of further targets, such as the Mars and asteroids.
Some orbits in the vicinity of the Moon with favorable characteristics, such as near rectilinear halo orbits (NRHOs) and distant retrograde orbits (DROs), are potential candidates for nominal parking orbits of long-term space stations. NRHOs, the candidate of the LOP-G, are members of the halo orbit families associated with the libration points L1 and L2 in the circular restricted three body problem (CR3BP) [2–7]. DROs are periodic orbits in retrograde motion around the Moon and are considered as another option for parking orbits of cis-lunar space stations [5, 8–13]. Although NRHOs and DROs have shown favorable stability characteristics, reliable station-keeping strategies are still required due to the long-term activity requirements and orbital perturbations. Different kinds of station-keeping techniques have been successfully developed for past missions in libration point orbits, including the Acceleration, Reconnection, Turbulence and Electrodynamics of the Moon’s Interaction with the Sun (ARTEMIS) [14, 15]. The recent program about long-term space stations is drawing much attention on the station-keeping strategies for NRHOs and DROs [16–20]. Besides the NRHOs and DROs, Condoleo et al. [21] have investigated the prospect of lunar missions for landing assistance by considering two quasi-circular orbits.
Most research about station-keeping strategies is focused on the libration point orbits, which can exist not only in the CR3BP but also in the elliptic restricted three-body problem (ER3BP) [22, 23]. Shirobokov et al. [24] have reviewed the station-keeping strategies for libration point orbits in a comprehensive manner. In their survey, the station-keeping techniques with impulsive or continuous thrust are distinguished as two categories: the methods that exploit dynamical insights of the three-body dynamics and the methods based on control theories. The Floquet mode approach, a classic impulsive station-keeping strategy in the first category, determines the unstable directions in phase space, and then finds an impulse that eliminates instability of motion in linear approximation [25–27]. Another station-keeping strategy in the first category, called Cauchy-Green Tensor (CGT) targeting approach, was developed by Guzzetti et al. [18] for NRHOs. The CGT targeting approach utilizes the state transition matrix to construct the CGT matrix, and then the accessible region under linear approximation is obtained by using eigenvectors of the CGT matrix. As a station-keeping strategy with continuous thrust in the first category, the pole-placement technique, proposed by Gurfil and Meltzer [28], exploits the Floquet theory to construct an optimal continuous control.
As for station-keeping strategies in the second category, different control theories, such as linear quadratic regulator (LQR) control and backstepping technique, have been applied [29, 30]. The LQR control was applied to the L2 halo orbits in the Earth-Moon system by Breakwell et al. [30]. To avoid solving the Riccati equation repetitively in the classic LQR control, Nazari et al. [31] combined the backstepping technique with the LQR control in the station-keeping of halo orbits. Qi and Ruiter [32, 33] extended this technique to halo orbits, vertical Lyapunov orbits, Lissajous orbits, and transfer orbits in the ephemeris model, and more practical constraints are included in their investigations. As a representative station-keeping strategy with impulsive thrust in the second category, the target point approach was developed by Howell and Pernicka [34]. By optimizing the performance function, an impulse that achieves the balance between the orbital error and station-keeping cost can be determined.
As one of the most successful modern control theories, the sliding mode control (SMC) has been widely applied in various fields, including orbital dynamics and control [35]. The typical SMC is conservatively designed for the worst scenario with the maximum system uncertainties and disturbances considered. But for space missions with limited fuel, the robustness is not the only concern, but minimizing the control cost is also essential in the mission consideration. The optimal sliding mode control (OSMC), which combines the optimal control and the SMC control, can serve as an ideal candidate for such scenarios. The optimal control can provide the balance between the orbital error and control cost, while the SMC control can provide robustness with respect to uncertainties and disturbances. The OSMC method has been successfully applied in different space scenarios, such as rendezvous [36], formation flying [37], and hovering [38], and it is also promising in the station-keeping scenario.
In this study, the station-keeping strategy for cis-lunar orbits using continuous thrust and the OSMC is investigated in the ephemeris model. Three types of cis-lunar orbits, including DROs, halo orbits, and NRHOs, are constructed in the ephemeris model. The OSMC controller based on the LQR control is designed, and the LQR control is also considered for comparisons. These two station-keeping strategies are first examined in simulations only considering the initial insertion error, and some unknown perturbations are included later to demonstrate the robustness of the OSMC controller. Then, some practical constraints caused by the navigation and propulsion systems are further included. Station-keeping performances for different nominal orbits are compared and analyzed with Monte-Carlo simulations. Finally, the practical constraints and unknown perturbations are both included to further explore the performances of the OSMC controller.
The rest of the paper is organized as follows: The ephemeris model and the nominal orbits are introduced in Section 2 as the background; in Section 3, the OSMC controller is designed based on the LQR method; in Section 4, simulations only considering the initial insertion error are carried out to present the performances of the proposed station-keeping controller; in Section 5, more practical constraints are included to verify the controller in realistic scenarios; finally, the paper is concluded in Section 6.
2. Background
In this section, the high-fidelity ephemeris model is introduced, and the nominal cis-lunar orbits of the station-keeping problem are presented.
2.1. Ephemeris Model
The Moon-centered J2000 inertial coordinates are used here to study the motion of spacecraft, as shown in Figure 1. The Earth-centered J2000 inertial coordinates are defined with the Earth’s mean equator and mean equinox at 12 : 00 terrestrial time on 1 January 2000. The Moon-centered J2000 inertial coordinates have the same coordinate axes with the Earth-centered J2000 inertial coordinates, but its origin is centered at the Moon. The spacecraft is denoted as , the central body, i.e., the Moon, is the denoted as , and other perturbing bodies, such as the Earth, the Sun, and the Jupiter, are denoted as . The relative position of each perturbing body with respect to the Moon is instantaneously obtained by employing the NAIF SPICE software and the DE 438 file [39].
The equation of motion of the spacecraft in the Moon-centered J2000 inertial frame is expressed as where and are position vectors from the central body, i.e., the Moon, to the spacecraft and the th body, respectively, and is the position vectors from the spacecraft to the th body. and are the nonspherical part of the lunar gravity and the SRP acceleration, respectively. The lunar gravity field is modeled by using the spherical harmonics truncated to the eighth degree and order, and the coefficients in the GRAIL (GRGM660PRIM) model [40] are exploited. The SRP acceleration can be obtained by [41] where and are the surface area and the mass of the spacecraft, is the position unit vector from the Sun to the spacecraft, is a characteristic constant, is the speed of light, is the distance between the Sun and the Earth, are the reflection and scattering rate of the spacecraft, and is normally considered to be zero. In this study, the mass, the surface area, and the reflectivity of the spacecraft are 1850 kg, 34.675 m2, and 0.9, respectively, which are same with the SOHO spacecraft.
2.2. The Cis-Lunar Orbits
The cis-lunar orbits are considered as the nominal orbits in this investigation, including five DROs denoted as the DRO (1)-(5), five halo orbits denoted as the halo (1)-(5) orbit, and five NRHOs denoted as the NRHO (1)-(5). The DROs and NRHOs are stable or weakly unstable, and the unstable halo orbits are considered for comparison. The 2 : 1 DRO, the halo (3) orbit, and the 9 : 2 NRHO are used as main nominal orbits. The 2 : 1 DRO is a planar orbit in retrograde motion around the Moon, and the 9 : 2 NRHO is a spatial orbit with the perilune close to the Moon, which is also the nominal orbit for the Lunar Orbital Platform-Gateway (LOP-G). It is worth mentioning that the 2 : 1 DRO and the 9 : 2 NRHO are resonant with the lunar sidereal period and synodic period, respectively. All nominal orbits in the CR3BP are shown in Figure 2 with nondimensional (nd) units, and their characteristics are listed in Table 1. The main nominal orbits are represented by the black lines, while other nominal orbits are represented by the blue lines.
The nominal orbits in a low-fidelity model, i.e., the CR3BP, cannot be used in real mission scenarios. The transition from the CR3BP to the ephemeris model is pivotal in the construction of nominal orbits. Several numerical methods are developed to achieve the transition, such as differential correction [18], two-level corrections [42], forward-backward shooting method, and receding horizon method [43]. All nominal orbits with a two-year duration are obtained by using differential correction in this study. The nominal orbits refined in the ephemeris model are shown in Figure 3, and their characteristics are listed in Table 2.
3. OSMC Controller for Station-Keeping
In this section, the OSMC controller based on the LQR method is designed for station-keeping, with the LQR controller as a comparison.
The equation of motion of the spacecraft with continuous control input can be expressed as where is the state vector of the spacecraft, is all forces considered in the dynamical model, consists of a zero matrix and an identity matrix, is the control variable, and is the perturbation vector including the external disturbances. The state deviation with respect to the nominal orbit is denoted as , and the dynamical equation with linearization can be expressed as where is the state of the nominal orbit, and since the nominal orbit is determined, we denote as for concise.
The following quadratic performance index for the classical linear quadratic regulator (LQR) control technique is considered for the dynamical system (Equation (5)) where is a positive semi-definite matrix and is a positive definite matrix. Such a performance index includes two major factors, the state deviation and the station-keeping cost. By adjusting the elements in two weighted matrixes, it is convenient to reduce the station-keeping cost or improve the tracking performance. Since the dynamical model has been nondimensionalized, the position error, velocity error, and control inputs have similar orders of magnitude. Thus, to achieve favorable tracking performance, the weighted matrixes are set as and . Then, the optimal feedback control law can be derived as follows where is determined from the following Riccati equation
With the presence of strong perturbations, the optimal control obtained from the LQR technique cannot stabilize the dynamical system, and then, the OSMC controller is proposed. The sliding surface is designed as the integral form where is the initial state, , and is nonsingular. The matrix is set to be .
The equivalent control obtained from is which is same as Equation (7). Thus, when the desired sliding surface is satisfied, the control input will be optimal for the predefined quadratic performance index. In order to improve the robustness of the controller, the control input is complemented as follows where , and its elements can be determined by considering the Lyapunov function . The derivate of the Lyapunov function is
By assuming that , one can get
Thus, the asymptotical stability can be achieved by properly choosing the values of . When the Jupiter gravity is considered as the perturbation, the magnitude of which is about in nondimensionalized unit, and the matrix is set to be . When the SRP is considered as the perturbation, the magnitude of which is about in nondimensionalized unit, the matrix is set to be .
4. Simulations with Initial Insertion Error
In this section, the performances of station-keeping strategies are presented with the consideration of the initial orbital insertion error.
4.1. Simulations without Perturbation
In this subsection, simulations without perturbation are conducted to present performances of the controllers. The position and velocity errors of initial insertion are assumed to have normal distributions in the magnitude, and random distributions in the direction. The magnitudes of position and velocity errors of initial insertion are considered as 100 km and 1 cm/s (), respectively. The main nominal orbits with a one-year duration are used here.
The simulation results of all nominal orbits are displayed in Figures 4–6, respectively. The position error vector with respect to the nominal orbit is denoted as , and the position deviation is expressed as . The velocity error vector with respect to the nominal orbit is denoted as , and the velocity deviation is . The control input vector is , and the acceleration provided by thrusters is .
As shown in Figures 4–6, there is no clear difference between the station-keeping performances by using the OSMC and LQR controllers. The position deviations and the control inputs of the 2 : 1 DRO, the halo (3) orbit, and the 9 : 2 NRHO decrease to near zero in 20, 20, and 35 days, respectively. The maximum magnitudes of the control inputs of the 2 : 1 DRO and the halo (3) orbit are about , while the maximum magnitudes of the control inputs of the 9 : 2 NRHO is about . Such significant difference is caused by abrupt peaks of control inputs near the perilune of the 9 : 2 NRHO. As shown in Figure 6, it is noticeable that the magnitudes of the velocity deviation and the control input are strongly relevant. The velocity varies more rapidly near the perilune than in other regions of the 9 : 2 NRHO. The velocity deviation near the perilune could be near 1 m/s with a nonsignificant position deviation, and then, larger thrusts are required to maintain the spacecraft near the perilune. This could explain the abrupt peaks of the velocity deviation and the control input near the perilune. In conclusion, with only the initial insertion error included, all nominal orbits can be maintained with the OSMC and LQR controllers, while larger magnitudes of the velocity deviation and the control input have been found in the case of the 9 : 2 NRHO.
4.2. Simulations with Perturbations
In this subsection, some perturbations are considered to test the robustness of the OSMC controller. Since unknown perturbations in real mission scenarios are difficult to simulate, the Jupiter gravity and the SRP are taken as the weak and strong perturbations, respectively. The dynamical models without the Jupiter gravity or the SRP are constructed, and the corresponding nominal orbits are obtained.
The simulation results with nominal orbits constructed in ephemeris models without the Jupiter gravity or the SRP are displayed in Figures 7 and 8, respectively. As shown in Figure 7, a small perturbation with a similar magnitude of the Jupiter gravity does not affect the convergence of the position deviation. Both station-keeping controllers can maintain the nominal orbits. However, as shown in Figure 8, when a perturbation with a similar magnitude of the SRP is considered, the position deviation cannot decrease to near zero with the LQR controller, whereas the OSMC controller has shown better performance and robustness. The position deviations can decrease to near zero in a short time for all the nominal orbits by using the OSMC controller. That is to say, the OSMC controller is more robust against strong perturbations than the LQR controller.
(a) 2 : 1 DRO with LQR controller
(b) 2 : 1 DRO with OSMC controller
(c) 9 : 2 NRHO with LQR controller
(d) 9 : 2 NRHO with OSMC controller
(e) Halo (3) with LQR controller
(f) Halo (3) with OSMC controller
(a) 2 : 1 DRO with LQR controller
(b) 2 : 1 DRO with OSMC controller
(c) 9 : 2 NRHO with LQR controller
(d) 9 : 2 NRHO with OSMC controller
(e) Halo (3) with LQR controller
(f) Halo (3) with OSMC controller
5. Station-Keeping under more Practical Constraints
In this section, the practical constraints caused by the navigation and propulsion systems are introduced, and the station-keeping simulations under these practical constraints are carried out with and without unknown perturbations.
5.1. Practical Constraints
In real mission scenarios, the capability of the navigation and propulsion systems are limited, and then, some practical constraints are necessary to consider in the verification of the station-keeping. The practical constraints and their parameters used in the Monte-Carlo simulations are introduced as follows.
The first constraint is that the initial orbit cannot be precisely inserted. The magnitude of initial insertion errors on position and velocity is considered to be 100 km and 1 cm/s () in the Monte-Carlo simulations. Another constraint is that the navigation system cannot provide real-time information about position and velocity, and the navigation interval is considered to be 2 days. The navigation errors on position and velocity are also nonnegligible limitations, and their magnitudes are considered to be 1 km and 1 cm/s (), respectively. The above parameters in constraints are selected by referring to previous investigations of station-keeping problems with both the impulsive and continuous thrust.
The practical constraints caused by the propulsion system are also worth consideration. The first constraint for the propulsion system is its engine limitation, and the magnitude of the thrust acceleration is limited within the range of . Thus, the real control input can be rewritten as
Since the thruster cannot precisely provide the control input, the propulsion error is also considered in the station-keeping simulations, and its magnitude is assumed to be 2%.
In the following station-keeping simulations, the initial insertion error and navigation error are assumed to have normal distributions in the magnitude and random distributions in the direction. All constraints discussed above and other parameters are listed in Table 3, which are same with those used in Qi and Ruiter (2019) [32].
5.2. Station-Keeping without Perturbation
The simulation results of station-keeping for all nominal orbits under practical constraints in the absence of perturbation are shown in Figures 9–11, respectively. Different with the results in Section 4, the position deviations cannot decrease to near zero due to the practical constraints. There is no obvious difference between station-keeping performances of the LQR and OSMC controllers. As shown in Figure 9, the steady position deviations of the 2 : 1 DRO are smaller than those of the 9 : 2 NRHO. It is observable that the peaks of control inputs for the 9 : 2 NRHO also exist.
(a) 2 : 1 DRO with LQR controller
(b) 2 : 1 DRO with OSMC controller
(a) Halo (3) with LQR controller
(b) Halo (3) with OSMC controller
(a) 9:2 NRHO with LQR controller
(b) 9:2 NRHO with OMSC controller
The Monte-Carlo simulations with parameters in Table 3 are used to obtain more reliable results. In the Monte-Carlo simulations, 200 sample points with random initial insertion errors and practical constraint parameters are chosen, and then, the results are analyzed by using statistic methods. The performance indexes for different nominal orbits are listed in Table 4. The annual cost is the average of the annual station-keeping cost of each sample, and the mean deviation and the max deviation are the average values of the mean and maximum position deviation of each sample, respectively. According to Table 4, the 2 : 1 DRO and the halo (3) orbit can be maintained by both station-keeping controllers, while some failed cases occurred for the 9 : 2 NRHO. The station-keeping performances of these two strategies do not show clear differences, which indicates that the OSMC controller possesses optimality in the quadratic performance index of the LQR controller. The station-keeping cost of the 9 : 2 NRHO is much higher than other nominal orbits, which could be caused by the thrust peaks near the perilunes.
As shown in Section 4, the position deviation of the 9 : 2 NRHO can decrease to near zero in the station-keeping simulations only considering the initial insertion error. The practical constraints have resulted in failed cases for the 9 : 2 NRHO, and then, the effects of different practical constraints are investigated. Four limitations of propulsion and navigation systems could cause the failed cases, including the engine limitation, propulsion errors, navigation interval, and navigation errors. Seven combinations of these four limitations are examined to find the main reason of the failed cases.
The case (a) has no engine limitation, which means that there are no upper or lower bounds for the control inputs. No propulsion error is considered in the case (b), which means that control inputs within the upper and lower bounds can be given precisely. The case (c) is the combination of the cases (a) and (b), which means that the propulsion system can provide any control input precisely. The case (d) has reduced the navigation interval from 2 days to 6 hours, which can reduce the state errors caused by the linearization of the dynamic model. The case (e) has no navigation error, which means that the navigation system can provide the accurate state information with a certain interval. No navigation limitation is considered in the case (f), and then, the spacecraft can have the real-time accurate state information. The case (g) has ignored all limitations caused by the propulsion and navigation systems, which is the same as the case only considering the initial insertion error. The Monte-Carlo simulation results of these seven cases are listed in Table 5.
As shown in Table 5, the results of the cases (a), (b), and (c) indicate that the improvement of the propulsion system cannot achieve a 100% success rate. There is no failed case in the case (d), indicating that the linearization of the dynamical model has caused unneglectable effect for the 9 : 2 NRHO. In the case (e), failed cases still exist, but both the station-keeping cost and position deviation have significant reductions, indicating that station-keeping performances are closely relevant to navigation errors. By ignoring both constraints caused by the navigation system, there is no failed case in the case (f), and the station-keeping performance is also desirable. There is no practical constraint in the case (g), and the station-keeping cost is much lower than other cases, indicating that eliminating the initial insertion error only requires a quite low cost. That is to say, the high station-keeping cost of the 9 : 2 NRHO is mainly used to overcome errors caused by the navigation and propulsion systems. Thus, the failures in the station-keeping of the 9 : 2 NRHO are caused by the distinct dynamical characteristic and the limited capabilities of the navigation and propulsion systems, and then, it is pivotal to improve the performances of the navigation and propulsion systems.
Then, the OSMC controller applied to all nominal orbits and their station-keeping performances are listed in Table 6. With the increase of the period of the DROs, the amplitude of the DROs increases, the station-keeping cost decreases significantly, and the average deviation and maximum deviation have no obvious changes. As the period of halo orbits decreases, their stability indexes also decrease, but their station-keeping costs and mean deviations increase gradually. The NRHOs belong to halo orbits, and their periods are shorter than the selected halo orbits. With further decrease of the period of the NRHOs, their station-keeping costs increase more significantly. It is worth mentioning that with the decrease of the period, the halo orbits gradually get closer to the Moon, and the gravity of the Moon is more significant, and then, higher thrust is required to adjust the orbital motion of the spacecraft. Thus, the variation of the station-keeping costs of the halo orbits and NRHOs indicates that the dynamical environment of the nominal orbit has more effect on the station-keeping cost, rather than their stability indexes.
5.3. Station-Keeping with Perturbations
In this subsection, some perturbations are included in the station-keeping simulations under practical constraints. The same as Section 4.1, the Jupiter gravity and the SRP are considered to represent the weak and strong perturbations, respectively. Since failed cases occurred in the station-keeping of the 9 : 2 NRHO with both station-keeping strategies, only the 2 : 1 DRO and the halo (3) orbit are used as nominal orbits here. The results of station-keeping simulations with the Jupiter gravity and the SRP as the perturbation are shown in Figures 12 and 13, respectively, and their performance indexes are listed in Table 7. As shown by Table 7, the perturbation of the Jupiter gravity has no significant effect on the position deviation and annual cost, whereas the perturbation of the SRP has significant effects.
(a) 2:1 DRO with LQR controller
(b) 2:1 DRO with OSMC controller
(c) Halo (3) with LQR controller
(d) Halo (3) with OSMC controller
(a) 2:1 DRO with LQR controller
(b) 2:1 DRO with OSMC controller
(c) Halo (3) with LQR controller
(d) Halo (3) with OSMC controller
For both the 2 : 1 DRO and the halo (3) orbit, the position deviations of the OSMC controller are lower than those of the LQR controller. The advantage in terms of the position deviation of the OSMC controller has also been verified in the Monte-Carlo simulations with parameters in Table 3. Due to the practical constraints, the advantage of the OSMC controller on the position deviation is not as significant as the simulations only with the initial insertion error. However, the station-keeping cost of the OSMC controller is higher than the LQR controller. That is to say, when strong perturbations exist, the station-keeping with the OSMC controller has a lower position deviation but a higher cost than the LQR controller.
6. Conclusion
In this paper, the station-keeping strategies with continuous thrust have been applied to the cis-lunar orbits in the ephemeris model. The cis-lunar orbits including DROs, halo orbits, and NRHOs are used as nominal orbits. The OSMC controller based on the LQR method is designed, and the LQR controller is also considered for comparison.
The simulations only considering the initial insertion errors have shown that all nominal orbits can be maintained by the OSMC and LQR controllers in the cases of no perturbation or a weak perturbation with a similar magnitude of the Jupiter gravity. When a strong perturbation with a similar magnitude of the SRP is included, the OSMC controller is more robust than the LQR controller.
The simulations with more practical constraints caused by the navigation and propulsion systems have shown that the 2 : 1 DRO and the halo (3) orbit can be maintained by both strategies, while some failed cases have occurred for the 9 : 2 NRHO. It has been found that the failure in the station-keeping of the 9 : 2 NRHO is caused by the distinct dynamical characteristics and the limited capabilities of the navigation and propulsion systems. Also, the improvement of navigation performances can greatly enhance the station-keeping performances.
The station-keeping results of all nominal orbits indicate that the station-keeping costs of DROs decrease with the increase of orbital period, while the station-keeping costs of halo orbits and NRHOs increase with the decrease of orbital period. It is also found that the dynamical environment of nominal orbits has more significant effect on station-keeping performances than their stability indexes. The station-keeping cost of the 9 : 2 NRHO is much higher than those of other nominal orbits, which is caused by abrupt thrust peaks near the perilunes. As for the robustness in the cases of perturbation, when a perturbation with a similar magnitude of the SRP is considered, the OSMC controller has a lower position deviation but a higher station-keeping cost than the LQR controller.
Data Availability
Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (11872007).