Abstract
The satellite constellation-to-ground coverage problem is a basic and important problem in satellite applications. A group of judgement theorems is given, and a novel approach based on these judgement theorems for judging whether a constellation can offer complete single or multiple coverage of a ground region is proposed. From the point view of mathematics, the constellation-to-ground coverage problem can be regarded as a problem entailing the intersection of spherical regions. Four judgement theorems that can translate the coverage problem into a judgement about the state of a group of ground points are proposed, thus allowing the problem to be efficiently solved. Single- and multiple-coverage problems are simulated, and the results show that this approach is correct and effective.
1. Introduction
Currently, satellite constellations technology are widely used in many practical applications, including communication, navigation, remote sensing, and science missions [1, 2], and these applications provide very important support to many fields, such as the observation of the land, sea and air [3, 4], disaster monitoring [5, 6], and resource exploration and development [7]. To meet the requirements of these applications, some typical problems encountered when designing satellite constellations have been studied by scholars. These problems include constellation design [8, 9], constellation performance evaluation [10], satellite-ground communication [11], and satellite scheduling [12, 13]. In these typical problems, the service of the constellation to the ground region, which depends on the visibility between the constellation satellite and the ground region, is a basic and important problem. This problem is called the constellation-to-ground coverage problem.
Solutions to the constellation-to-ground coverage problem aim to calculate the coverage region and other information about the ground regions that are covered by satellite constellations during a time period or at a point in time. This problem is essentially a problem involving Boolean operations on 2D graphs [14, 15], which is commonly encountered in the fields of computer graphics and geology. However, compared with the traditional problems involving the Boolean operations on 2D graphs, the constellation-to-ground coverage problem has some differences because of the following three aspects: (1) the background manifold is a spherical or ellipsoidal surface in three dimensions rather than a 2D Euclidean space, (2) the shape of ground regions and the view field of the sensor can be arbitrary, and (3) the temporal characteristics during long simulation periods frequently arise in calculating the coverage region of sensor. These three aspects make this problem more difficult than traditional Boolean operation problems.
The solution of the coverage problem depends on the configuration of the constellation and the coverage problem itself. Different constellation configurations, including the star constellation [16], the Walker-Delta constellation [17, 18], the rosette constellation [19], and the flower constellation [20], may have some particular configuration characteristics, such as spatial symmetry [21, 22], common-track [20], regressive orbits characteristic [23], and polar orbits [24]. Proper usage of these characteristics can greatly increase computational efficiency. From the aspect of the coverage problem, owing to the different requirements of the time and space factors for different practical applications, constellation-to-ground coverage problems are different. For navigation and communication applications, the timeliness requirement is critical; however, for some remote sensing applications such as resource exploration, the timeliness requirement is minimal. Thus, the concepts of continuous coverage problems [8, 25] and discontinuous coverage problems [24, 26] have been proposed by many scholars, and thus, different constellation-to-ground coverage problems have different characteristics and different strategies for achieving a solution.
Until now, to solve constellation-to-ground coverage problems, there have been three typical methods, including the grid point method [27], the triangle mesh subdivision method [17], and the street-of-coverage method [28]. The triangle mesh subdivision method is a method that is favored by many scholars [18, 19]. This method can quickly judge whether the constellation can cover the globe at a given time. However, more detailed coverage information, such as the coverage rate, cannot be obtained. The street-of-coverage method is a method for solving the continuous coverage problem and can obtain the approximate analytic relationship between coverage and the constellation configuration so as to guide constellation design. The disadvantage of this approach is that the result is an approximate value and the method ignores phase information.
Additionally, some mathematical strategies have been applied in constellation-to-ground coverage problems, such as probability statistics [29], iterative methods [30], the differential quadrature method [31], the method based on the transformation group [21, 22], the polyhedron surrounding method [32], the method based on auxiliary space [26, 33], and the numerical fitting method [34]. These mathematical strategies can be used in some of the three aforementioned methods, which can accelerate computing efficiency. However, most of these mathematical strategies have some restrictions in that they can only work in some certain conditions and can only solve some certain problems.
In this paper, we focus on a subproblem of the constellation-to-ground coverage problem, the complete coverage problem, which judges whether a constellation can cover a ground region completely using one or multiple satellites. If the constellation cannot cover a ground region completely, we do not care about the coverage rate or other coverage properties.
In this paper, we will discuss this problem from the aspect of mathematics. From the point view of mathematics, the instantaneous coverage region of the ground region from a satellite can be seen as a spherical region. Thus, the constellation-to-ground coverage problem can be translated into a problem on the intersections of spherical regions. According to the basic characteristics of the coverage problem, a group of definitions are given, and a group of theorems as well as proofs are put forward. Then, according to the judgement theorems, highly accurate, efficient, and novel approaches for solving this problem are proposed.
The paper is organized as follows. Section 2 gives a brief review of the traditional methods. Section 3 introduces the basic concepts and appointments of this paper. Section 4 proposes a group of judgement theorems for this coverage problem. Section 5 discusses the detailed procedures of the approaches for solving this coverage problem. In Section 6, a group of experiments, including single-coverage problems and multiple-coverage problems, are conducted. Section 7 concludes the paper.
2. Background
The traditional methods for solving the complete coverage problem include the grid point method and the triangle mesh subdivision method. A brief introduction to these two methods will be given.
2.1. The Process of the Grid Point Method
The grid point method is one of the most commonly used approaches for solving the constellation-to-ground coverage problem. The main idea of this approach is discretization, which is the division of a continuous object into a group of discrete objects and then calculating the characteristics of these discrete objects to express the characteristics of the coverage problem.
For a coverage problem, two types of objects need to be discretized: the first is space, and the other is time. The discretization of space means discretizing the ground region into a group of grid points according to a regional discretization strategy, and the discretization of time means discretizing a time range into a set of time points.
For the grid point method, the process is as shown in Figure 1, and the basic steps are as follows:
Step 1. For ground region , a minimum longitude and latitude box that can completely contain is found (the longitude and latitude box corresponds to the blue rectangle in the second panel of Figure 1).
Step 2. Given a space partition precision and the chosen sampling strategy, a group of feature points are chosen and are called grid characteristic points, which are denoted by , where . corresponds to the red and greed dots in the second panel of Figure 1.
Step 3. For every grid characteristic point belonging to , whether this point belongs to the ground region is judged. Then, all grid characteristic points that belong to are gathered, and these points are written as . corresponds to the green dots in the second panel of Figure 1, while the red dots represent the grid characteristic points that do not belong to .
Step 4. For every grid characteristic point belonging to , whether this point can be covered by the constellation is judged. Then, all grid characteristic points that belong to are gathered, and these points are written as . corresponds to the black dots in the third panel of Figure 1, while the red circles stand for the coverage region of constellation .
Step 5. The number of points in and are counted, and the ratio of these two sets is regarded as the coverage rate.
This process of the grid point method does not address the temporal characterization. The temporal characterization is hidden in Step 5, the process of judging whether a grid characteristic point can be covered by the constellation. The basic approach for solving the temporal characterization also involves discretization. For the grid point method, the basic steps of discretizing the temporal elements are as follows:
Step 1. Given a time partition precision corresponding to time interval , a group of time points can be obtained from the time interval , where , , and for any .
Step 2. For any time point , the coverage state of every grid characteristic point belonging to is checked for whether each point can be covered by the constellation at the current time.
Step 3. The coverage characteristics of all time points are analyzed, and statistics are calculated for these time points. Then, the temporal characteristics of the constellation-to-ground region can be obtained.
The greatest advantage of the grid point method is that it has a wide response domain, which means that it can solve nearly all coverage problems. However, when high computational accuracy is required, this approach is highly temporally and spatially complexity. The temporal and spatial complexity grows exponentially as the required calculation accuracy increases. This feature is one of the biggest drawbacks of the grid point method and greatly restricts the practical application of this approach.
Additionally, the grid point method is good at computing the value of the coverage rate, but when solving complete coverage problems, the grid point method is not efficient as it must compute a great amount of information that will not be used when solving complete coverage problems.
2.2. The Spherical Triangle Mesh Subdivision Method
The main idea of the spherical triangle mesh subdivision method is to divide the globe into a set of spherical Delaunay triangles, where the vertexes of the set of spherical Delaunay triangles are the points beneath the satellites in the constellation, and then checking whether the constellation can completely cover a spherical Delaunay triangle given the radius of the spherical circumcircle of the spherical Delaunay triangle.
Spherical Delaunay triangulation is an important pretreatment technique in both graphical and numerical analyses. By following certain rules, the globe can be divided into adjacent and mutually disjoint triangles that have the following characteristics:(i)The minimum angle of the spherical triangles is maximized. This characteristic can make the subdivision result more uniform and regular.(ii)The triangles have empty interiors. The characteristic of empty interiors is that the interior of every spherical triangle circumcircle does not contain other vertexes.
The subdivision result is generally unique. For vertexes, the globe can be divided into spherical triangles, where each has edges.
The basic process of the spherical triangle mesh subdivision method is as follows:
Step 1. Calculate the instantaneous location of the subsatellite point of all the satellites in the constellation.
Step 2. Take the instantaneous subsatellite position of all the satellites in the constellation as the vertex set and then construct a spherical Delaunay triangle subdivision of the globe.
Step 3. For any Delaunay triangle, calculate the radius of the spherical circumcircle and obtain the maximum value, which is denoted by .
Step 4. Denote the instantaneous coverage radius on the earth of the satellite by . If , this indicates that the constellation cannot fully cover the global. If , then the constellation can completely cover the globe.
The spherical triangle mesh subdivision method is a method for solving complete coverage problems. However, the target of complete coverage problems must be the globe. However, the method does not work well for other types of targets, such as zonal targets, polar cap targets, and local area targets. Additionally, this method can only solve problems in which of every satellite is equal.
3. Basic Appointment and Definition
Appointment 1 (ground region). The ground region is appointed as a connected spherical domain and also as a closed domain.
Denote a spherical domain by , and denote the boundary of by . Because the spherical domain is appointed as a closed domain, then .
Appointment 2 (satellite coverage region). The satellite coverage region is appointed as a connected spherical domain as well as an open domain.
Let be a satellite. Denote the coverage region of by , and denote the boundary of by . According to Appointment 2, is an open domain; then , and . The union set of and is called the closure of and is denoted by .
In actual applications, we may also encounter the situation in which the ground region or satellite coverage region is not a connected domain but a nonconnected domain. Then, this nonconnected ground region or satellite coverage region can be equivalently seen as a group of connected ground regions or satellite coverage regions, and this treatment has no influence on the results of the analysis. Thus, in this paper, for the ground region or satellite coverage region, only the connected spherical domains are considered.
When is a single connected domain, consists of a connected curve. However, when is a multiple connected domain, consists of more than one single connected curve. Let be the th boundary curve of ; then , where is the number of connected boundary curves.
Because the position of satellite changes with time, is a time-varying region and has hidden time parameters , which can be written as when the time parameters are emphasized.
A constellation is a set of satellites that has a relatively stable spatial configuration. Denote a constellation by . Then, similarly, the coverage region of the constellation is written as , and the boundary of the constellation coverage region is written as . Then, the following equation holds:
However, because some points of may be covered by other satellites belonging to , which means that these points are no longer boundary points. Then, .
Generally, is not a single connected domain but a multiple connected domain. A connected domain of is called the coverage branch of . Let be the count of the connected domain of , and denote the th coverage branch of by ; then,
Then, is less than or equal to the number of satellites in constellation .
Appointment 3 (cover). A ground point denoted by is said to be covered by satellite when is an inner point of , which can be expressed by according to Appointment 2. Additionally, a ground point is said to be covered by constellation when .
According to Appointments 2 and 3, satellite cannot cover a point located on .
Generally, in actual applications, these three concepts do not need to be defined so rigorously. However, because we will solve this problem from the point view of mathematics, the rigorous definition of the concepts can make the forms of the expression and process of the proof simpler.
Then, some definitions that will be used later in this paper are given.
Definition 4 (connected satellites). Let be a time point. For two satellites and , if and belong to the same branch of , then is said to be connected to at time , which is denoted by , where is a hidden parameter of this relationship. and are called connected satellites.
For set , we can see that the connectivity of satellites can be regarded as an equivalence relationship. The equivalence class of satellite using the equivalence relationship is written as . For any satellite belonging to , belongs to the same connected branch of .
Definition 5 (connected coverage region). The connected coverage region of satellite in constellation is the union of the coverage region of satellites belonging to and is denoted by . Similarly, the boundary of is written as . Then, evidently, , and . Additionally, if , then .
Definition 6 (complete coverage). Constellation is said to completely cover the ground region at time when , such that .
Strictly speaking, complete coverage is called single-satellite complete coverage. There also exists another type of complete coverage, which is called multiple-satellite complete coverage.
Definition 7 (multiple-satellite complete coverage). Constellation is said to possess -multiple complete coverage of the ground region at time when , such that , where .
4. Coverage Problem Judgement Theorems
Theorem 8. The necessary and sufficient conditions for constellation to completely cover ground region are as follows: (1) such that .(2).
Proof. For the sufficiency, ; then can completely cover .
Because constellation can completely cover ground region , then condition (1) is clearly satisfied. For condition (2), supposing , then . Then, such that , and then cannot be covered by constellation . This condition contradicts the condition that constellation completely covers ground region , so the assumption is invalid. Thus, .
Then the theorem is proved.
Then, according to Theorem 8, if constellation can completely cover region , then constellation has one and only one coverage brunch that can contain . Additionally, if constellation can completely cover region , then, for any coverage branch , either or is held.
Theorem 8 reveals some characteristics of the complete coverage problem, and this theorem can be used to prune the branches.
Theorem 9. The necessary and sufficient conditions for constellation to completely cover ground region are as follows: (1) such that .(2).
Proof. For the sufficiency, if there exists a satellite such that , then, supposing , because , then . Then, the relationship between and may have two statuses, which are and , but . If , because is an open domain and is a closed domain, then , and then . If but , then . Then, the statuses both contradict condition (2). So, this assumption is invalid. Then, . According to Theorem 8, constellation can completely cover ground region .
For the necessity, when constellation can completely cover ground region , then condition (1) is obviously satisfied. Then, according to Theorem 8, . According to Appointment 2, it satisfies that ; then .
Generally, and are both 2-dimensional spherical domains, and the coverage problem of to is an intersecting problem of two two-dimensional spherical domains. Then, Theorem 9 translates this coverage problem into an intersecting problem of a two-dimensional domain and a one-dimensional domain, which is much easier to solve than the original problem.
However, the problem of intersecting a two-dimensional domain and a one-dimensional domain is also somewhat complex. We then introduce the boundary feature point and give another theorem, thus making the problem easier.
Before we discuss the theorem, a definition that will be used in the theorem should be defined.
Definition 10 (boundary feature points). The boundary feature points of satellite are a group of special points located on and are denoted by . The boundary feature point on the th boundary curve, which is denoted by , is divided into two conditions:
Condition 1. For , if or , then we select as the boundary feature point.
Condition 2. For , if and , we select an arbitrary point on as the boundary feature point.
Then, .
Similarly,Per the definition of a boundary feature point, we can see that consists of finite points when we do not consider the critical condition that and overlap with each other.
Theorem 11. The necessary and sufficient condition for constellation to completely cover ground region is as follows: (1) such that .(2), such that .
Proof. For the necessity, the first condition is obvious. For the second condition, per Definition 10, ; then, because constellation can completely cover ground region , the second condition holds.
For the sufficiency, from Figure 2, we can see that the boundary of the uncovered region consists of and and that the vertexes of the uncovered region are made up of the intersection of the coverage boundary of two satellites as well as the intersection of the coverage boundary of a satellite and the boundary of the region. Then, we can see that the vertexes of the uncovered region approximate the boundary feature points of the constellation. Intuitively, supposing that constellation does not completely cover ground region , if the uncovered region has vertexes, then the vertexes must be the boundary feature points, and if the uncovered region has no vertex, then boundary of the region is just Condition 2 of Definition 10, and an arbitrary point on the boundary is selected as the boundary feature point. Then, the boundary feature points are not covered by another satellite, which contradicts Condition 2. Thus, the assumption is invalid. Then, constellation can completely cover ground region .
Then, by Theorem 11, the complete coverage problem can be translated into checking the coverage state of a finite number of points, which greatly reduces the computational difficulty. More importantly, this approach has no calculation error.
Additionally, this theorem can not only be used for solving the single-coverage problem but also can solve the multiple-coverage problem. Thus, a new theorem that is used to deal with the multiple-coverage problem will be proposed.
Theorem 12. The necessary and sufficient condition for constellation to possess -multiple complete coverage of ground region is as follows: (1) such that .(2), such that , where .
The proof of this theorem is similar to that of Theorem 11. Thus, we do not cover this process here.
5. Approach for Judging the Coverage State of the Constellation-to-Ground Region
5.1. Approach for Judging the Coverage State at a Certain Time
In the actual computing process, generally, we use Theorems 11 and 12 to solve the practical problem because the computational efficiency with this approach is much higher than when using Theorems 8 and 9.
The basic steps for judging whether constellation can singly or multiply cover ground region completely at a certain time are as follows:
Step 1. Check whether there exists a satellite such that . If does not exist, then it can be asserted that the constellation cannot completely cover ground region at time .
Step 2. Calculate the boundary feature points of constellation to ground region at time , which is denoted by .
Step 3. Judge whether all points belonging to can be covered by a satellite (or satellites in -multiple-coverage problem). If the result is yes, then can completely cover at ; else, cannot completely cover at .
5.2. Approach for Judging the Coverage State for a Period of Time
Theorems 11 and 12 can only judge the coverage state at a certain time. However, in actual engineering applications, time factors are often considered, including continuous coverage [8, 25].
Definition 13 (continuous coverage). Let be a period of time. If for any , the constellation can completely cover ground region at , then we say that the constellation can continuously cover ground region at .
Because the relevant judgement theorems about time factors have not been included in this paper, when dealing with the time factors, the discretization method is adopted. The basic steps for solving continuous coverage problems are as follows:
Step 1. Given a time partition precision corresponding to time interval , can be divided into a group of time points where , , and for any .
Step 2. For any time point , if constellation cannot singly (or multiply) cover a ground region completely at time , then cannot continuously cover at . If can completely cover completely at all of the time points, then can continuously cover at .
However, different from the approach for judging the coverage state at a certain time, the approach for evaluating the continuous coverage problem may have calculation errors when is not small enough. Thus, in order to control calculation errors, is chosen to be a sufficiently small number.
5.3. Calculation of the Critical Beam Angle
Theorems 11 and 12 can only obtain a result of “yes” or “no,” which is a qualitative result. However, in actual applications, a quantitative result is more useful. So, in this subsection, we will use Theorems 11 and 12 to obtain a typical relevant quantitative result by calculating the critical beam angle of the sensor.
Definition 14 (the critical beam angle). Consider an isomorphic constellation , which means that all satellites in the constellation have the same semimajor axis, eccentricity, and orbit inclination, and all satellites carry same type sensor. Then, for ground region , the critical sensor beam angle is the minimum beam angle of sensor with which the constellation can completely cover the ground region .
For constellation , if the half-beam angle of sensor equipped on the satellite that belongs to is , then the coverage region of at is written as . Then, clearly, if , then ; that is, the results are monotonous relative to the half-beam angle of sensor. So, the critical beam angle for to can be obtained by the bisection method.
Denote the radius of the Earth by , and denote the semimajor axis of every satellite by ; then, given an elevation constraint of the sensor, the maximum half-beam angle denoted by is , and the minimum half-beam angle denoted by is 0.
Step 1. For constellation , calculate and .
Step 2. If , then the algorithm terminates and the critical half-beam angle does not exist.
Step 3. Let be equal to . Then, check the coverage state of with half-beam angle to at ; then, if , then ; else, .
Step 4. Let be the allowable error, which is given according to the requirements of the actual problem. If , then the algorithm terminates, and is the ultimate critical half-beam angle; else, jump to Step 3.
The critical half-beam angle of to at is denoted by .
Then, for a time period , the critical half-beam angle of to , which is denoted by , is defined as the maximum value of the critical half-beam angle of to at every moment in . Then,
6. Experiment
To verify the correctness and effectiveness of the theorem discussed in this paper, a group of experiments are designed.
For a certain constellation, the time for judging whether the constellation can cover a ground point is nearly a constant, so the computing time is associated with the number of ground points that need to be judged. So, to measure the effectiveness of the algorithm, we establish a performance index, which is the number of boundary feature points.
6.1. Single-Coverage Problem
Consider a Walker-Delta constellation with a configuration of 45/5/1. All satellites in the constellation have circular orbits and have the same altitude, which is 1300 km. The ground region is a latitude zonal target, the upper and lower bounds of which are 50° and 0°, respectively. Consider a simulation time period that lasts for three minutes with the start time 2017-01-01 0:0:0 and the end time 2017-01-01 0:3:0.
The variation curve of the average number of boundary feature points and the orbital inclination during the simulation time period are presented in Figure 3. Then, the variation curve in Figure 3 represents single-peak characteristics. When the orbital inclination is = 90°, the number of boundary feature points reaches the minimum. The maximum number of boundary feature points is lower than 500; that is to say, if we want to judge whether the constellation can completely cover the ground region, no more than 500 points need to be checked. Thus, the efficiency for judging whether the constellation can completely cover the ground region is very high.
The variation curve of the critical sensor beam angle with orbital inclination is presented in Figure 4. Then, from Figure 4, we can see that the curve is nearly symmetric along the axis of = 90°. The curve has three distinct troughs and four distinct peaks. The global minimum value is achieved when the orbital inclination equals 37°, at which point the corresponding critical beam angle is 54.6°.
To verify the correctness of the algorithm, we calculated this coverage problem by means of the grid point approach. We considered an orbital inclination for all satellites of 37°, and then we calculated the relationship between the continuous coverage rate during the simulation time period and the sensor beam angle of all the satellites.
Figure 5 illustrates the variation of coverage rate over orbital inclination. The -axis represents the orbital inclination and the -axis gives the corresponding coverage rate in the graph. It denotes that the constellation can continuously cover the ground region, while the orbital inclination is over 54.6°. The result shows that this approach is correct for solving the single-coverage problem.
6.2. Multiple-Coverage Problem
Consider a Walker-Delta constellation with a configuration of 18/3/0. The altitude of the satellites belonging to the constellation is 25000 km. Consider the global quadruple continuous coverage problem with the simulation time period 2017-01-01 0:0:0–0:53:0, which is the PU of this Walker constellation.
The variation of the critical beam angle with orbital inclination is presented in Figure 6. For this multiple-coverage problem, the variation of critical beam angle with orbital inclination is very similar to that in single-coverage problems.
The average number of boundary feature points during the simulation time period with the change in orbital inclination is as in Figure 7. The curve is very regular and symmetric. The average number of boundary feature points is between 270 and 290.
Figures 6 and 7 reflect that the best orbital inclination is approximately 54.2°, where the critical sensor beam angle is 11.3°.
Additionally, to verify the correctness of the algorithm for solving multiple-coverage problems, the grid point approach is also used. We considered that the orbital inclination of all satellites is 54.2°; then, we calculated the global quadruple continuous coverage problem during the simulation time period.
From Figure 8, we can see that when the sensor beam angle is greater than 11.3°, the constellation offers quadruple continuous coverage of the globe. This result shows that this approach is also correct for solving the multiple-coverage problem.
7. Conclusion
In this paper, we discussed the problem of judging whether a satellite constellation can offer complete single or multiple coverage for a ground region. To solve this problem, a group of judgement theorems and definitions are proposed.
The numerical experiments indicated that the proposed approach can solve the problem correctly and efficiently. Compared with traditional approaches to solving this problem, such as the grid point approach and the spherical triangulation method, the process of this algorithm is simple, and the amount of computation is much less, which means that this approach has a high computational efficiency. Moreover, this approach to solving coverage problems at a given time has no computational error, so computational accuracy does not need to be considered.
However, one drawback is that this approach can only judge whether the constellation can offer complete single or multiple coverage for the ground region. If the constellation cannot cover the ground region completely, the coverage rate or other coverage characteristics cannot be obtained from this approach.
Additionally, the group of judgement theorems can only judge the coverage state at a certain time. However, by means of discretization, the coverage state for a period of time can be calculated. Future studies include proposing judgement theorems for solving coverage problems for a period of time without using the discretization method. Additionally, the accumulated coverage problem as well as the discontinuous coverage problem [24, 26] should be addressed.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Zhiming Song put forward the idea and wrote the main content of this manuscript; Zhiming Song, Guangming Dai, and Maocai Wang conceived the experiments and analyzed the results. Xiangyun Hu, Guangming Dai, and Maocai Wang inspected the manuscript. All authors reviewed the manuscript.
Acknowledgments
The authors thank Lei Peng, Xin Luo, Xiaoyu Chen, and Mingcheng Zuo for assistance with finishing the experiments. This work is supported by Natural Science Foundation of China under Grant no. 41571403, Joint Funds of Equipment Pre-Research and Ministry of Education of China under Grant no. 6141A02022320, and Fundamental Research Funds for the Central Universities under Grant no. CUG160828.