Abstract

The cooperative control in complex multitasks using unmanned aerial vehicle and unmanned ground vehicle (UAV/UGV) is an important and challenging issue in the multirobot cooperative field. The main goal of the task studied in this paper is to minimize the time and energy consumed by the system to complete the assigned tasks. In this paper, a complex multitask problem using the hybrid UAV/UGV system is studied, which is divided into three stages, namely, the stage of finding the optimal locations of the relay stations for the UGV; the stage of solving the path planning problem for the UGV; and the stage of the task assignment for multi-UAVs. Furthermore, an improved integrated method is proposed to deal with the cooperative control problems in these three stages. Firstly, an adaptive clustering method is proposed to determine the locations of the UGV relay stations, and then an improved cuckoo search algorithm is used to find the shortest path for the UGV. Finally, a grouping method is presented to solve the multitask allocation problem of UAVs, based on an improved dynamic programming algorithm. In addition, some simulations are carried out and the results show that the proposed method has better performance when it comes to the time and energy consumption and can effectively guide the hybrid UAV/UGV system to carry out the complex multitasks.

1. Introduction

The cooperative control of multirobots is a research hot spot in the field of robotics [14]. Recently, there have been more and more research and applications of the hybrid UGV/UAV system for the complex multitasks in some large-scale environments [57]. The heterogeneous collaboration technology can broaden the application field of UAV or UGV and improve the effectiveness of the detection, the search, and the rescue task [8, 9]. However, the heterogeneous cooperation of the multiple unmanned aerial vehicles and unmanned ground vehicles is a new and challenging field [10, 11], and the main problem is how to make full use of the existing advantages of the unmanned aerial vehicles and unmanned ground vehicles in different tasks.

There are multiple research results in the field of the hybrid UGV/UAV system. For example, Yu et al. [12] presented a cooperative path planning algorithm using a UAV/UGV system to track moving targets, where the vision occlusions due to obstacles in the environment were considered. Arbanas et al. [13] studied a symbiotic aerial vehicle-ground vehicle robotic team, where UAVs were used for aerial manipulation tasks, while UGVs aided and assisted them. Li et al. [14] proposed a hybrid algorithm for the automatic ground map building and efficient path planning in UAV/UGV cooperative systems. Cantelli et al. [15] presented architecture to allow cooperation between a ground robot and a quadrotor UAV by using the image processing algorithm, for surveying operations in humanitarian demining. In these applications of the UAV/UGV system, lots of new and difficult problems are met, such as the task assignment problem and the hybrid path planning problem.

To deal with these problems in the cooperative control of the UAV/UGV system, various methods have been proposed. For example, Manyam et al. [16] presented a branch-and-cut algorithm to solve the path planning problem for the ground and air vehicles, which was formulated as a mixed-integer linear program (MILP). Kamel et al. [17] proposed a fault-tolerant cooperative control strategy based on the Hungarian algorithm for a team of UAVs and UGVs in the presence of actuator faults. Ropero et al. [18] proposed a path planning algorithm for cooperative UAV/UGV exploration. Christie et al. [19] presented a system of radiation search using UAV and UGV by employing semantic scene segmentation. Those methods introduced above deal with some special problems in the UAV/UGV system, but few of them consider the cooperative control problem of the UGV/UAV system as a whole, which will affect the whole performance of the UAV/UGV system sometimes.

To study the cooperative control problem by the UAV/UGV system, a complex multitask scenario is presented in this paper. And the cooperative control problem of UAV/UGV in this multitask scenario is divided into three stages. The first stage is to find the proper relay stations for the UGV. In this stage, an adaptive clustering method is proposed. The second stage is to find the shortest path for the UGV, which is a traveling salesman problem (TSP). To deal with this problem, an improved cuckoo search algorithm is proposed. The last stage is to reach all the target areas and complete the tasks by UAVs, which is a task assignment problem. In this stage, an improved dynamic programming algorithm is proposed.

The main contributions of this paper are as follows: (1) A cooperative control scenario for multitasks using the hybrid UAV/UGV system is presented, where only one UGV and many UAVs are used. The proposed scenario resembles the real-world applications. (2) An improved integrated solution is proposed to deal with the problems in different stages of the cooperative control task. (3) Various simulations are conducted, and the results show that the proposed method can deal with the cooperative control problem efficiently.

This paper is organized as follows. Section 2 presents the problem statement of the cooperative control using the hybrid UAV/UGV system, and an integrated solution for it is proposed in Section 3. The simulations are given in Section 4. Section 5 discusses the parameter selection and the performance of the proposed approach by some comparison experiments. Finally, the conclusions are given in Section 6.

2. Problem Statement

In this paper, the cooperative control problem of a hybrid UAV/UGV system for complex multitasks is studied. The complex multitasks studied in this paper are similar to the real-world applications, such as the goods delivery and people search and rescue during wars or the rescue after disasters. In the real implementation, the cooperative controller can be embedded on the UGV. For ease of implementation, this problem is defined as follows. (1)There are a large number of target areas distributed in a large-scale environment, which is denoted by , where is the total number of the target areas(2)Each target area has some tasks, which require UAVs with different capabilities to complete. Then the -th target area can be denoted by ), where is the number of tasks in this target area. One task may need more than one capability to complete, which is denoted by , and is the number of the capability types required in the task. The distance of tasks in one target area is ignored, in view of the high speed of the UAV(3)The hybrid UAV/UGV system (denoted by ) is used to complete the tasks, which includes one UGV and many UAVs. The UAVs are labelled as , and is the number of the UAVs used in the tasks. The UAVs are given different capabilities used to complete the concrete tasks. Each UAV has just one type of capability, which is denoted by , where is the number of the capability types of UAVs. In this study, is equal to , and UAVs with the same capability in the UAV/UGV system are equal in number(4)The energy limit problem of a UAV can be represented by the longest distance (denoted by ) that the UAV carrying a sufficient amount of power can travel. For simplification without losing generality, the details of how the tasks to be solved by the UAVs are ignored and the time to fulfill a specific task by the UAV is negligible in this paper. Then, the maximum active range of a UAV can be expressed as a circle, with the (5)The UGV used in the system is a mobile platform that can carry multiple UAVs and charge them at the same time. In this paper, the UGV is assumed to have enough energy to support the whole system during the task. And the charging speed of the UAVs is assumed to be very high, which means that the charged process of all the UAVs can be finished during the transportation process(6)The cooperative control process to finish the tasks can be described as follows: at the beginning, a complex task is given to a hybrid UAV/UGV system . This task can be denoted by . The environment is large, which is denoted by . The UGV will carry all the UAVs to some proper positions (which are called the relay stations for UAVs), and the UAVs begin to perform the tasks in their range of capabilities around this position. After all target areas around this relay station are reached and all the tasks in these target areas are completed, the UGV will carry all the UAVs to the next relay station. When all the tasks in the environment are completed, the UGV needs to carry all the UAVs to the starting position then this cooperative control task is completed. An example of the cooperative control task studied in this paper is shown in Figure 1

The task in this paper is similar to those tasks studied in some related literature. For example, Mathew et al. [20] addressed a multirobot scheduling problem in which autonomous unmanned aerial vehicles (UAVs) must be recharged during a long-term mission. In their scenario, UAVs can dock with a separate team of dedicated charging robots (UGVs) in order to recharge. Yu et al. [21] presented a scenario, where the UAV can be recharged along the way either by landing on stationary recharging stations or on UGVs acting as mobile recharging stations.

However, there are some differences between the literature and the problem studied in this paper. For example, Ropero et al. [18] presented a task scenario. In their scenario, the UAV/UGV system was composed of one UGV and one UAV, and only one UAV was used to search for the targets. In the scene studied in this paper, there are two main problems to be solved. The first one is to decide the optimal positions of the relay stations for the UGV. The next one is the task assignment of targets for the UAVs. And the two tasks must be considered as an integer.

3. Proposed Approach

In order to accomplish the multitasks introduced in Section 2, the cooperative control process is divided into three stages: (1) Find the most suitable positions of the relay stations; (2) Plan the best path for the UGV; (3) Carry out the task assignment of the tasks around the relay stations for UAVs.

For each stage, a special solution is proposed to solve a separate problem. In the proposed approach, an adaptive clustering method is used to deal with the problem in the first stage, which is the basis of the proposed method. Then, an improved discrete cuckoo search algorithm is used to solve the problem in the second stage. Finally, a task assignment algorithm based on a grouping strategy is presented. The flow diagram of the cooperative control process for a hybrid UAV/UGV system based on the proposed method is shown in Figure 2. The proposed cooperative control method is presented in detail as follows.

3.1. The Clustering Algorithm in the First Stage

Because the energy of the UAV is limited, it is very important to determine the position of the relay stations for the UGV, which can reduce the total time and energy consumption in the task. In this paper, an adaptive clustering method based on the -means algorithm is proposed that is able to adaptively determine the center of the cluster, namely, the location of the relay station.

The basic idea of the -means algorithm is to divide the data into different clusters by iteration, so as to minimize the objective function. The -means algorithm can generate the clusters as compact and independent as possible. However, the general -means algorithm has some defects [22]. For example, it depends on the selection of initial values. To deal with this problem, the -means++ algorithm is used in this paper [23].

The -means++ algorithm proceeds with 2 steps, namely, the allocation step and the update step. In the allocation step, the target areas satisfying the constraint condition are assigned to different clusters. The equilibrium assignment task here is modeled as an integer linear programming (ILP). In the update step, the updated center points are assigned to different clusters as their centers.

In the cooperative control task studied in this paper, it is difficult to give out the value of directly, because the tasks in the target areas are different. This is different from the clustering problem studied in other literature. In this paper, an empirical equation is obtained based on the analysis of the relationship between the tasks and , which is as follows: where is a parameter that is proportional to the size of the environment and the complexity of the target tasks (namely, the number of the tasks in each target area). The larger the environment size or the more the tasks in one target area, the larger the value of will be. So is often set as 0.1 to 0.5.

To homogenize the efficiency of UAV in the tasks of the next stage, the size of these clusters should be constrained. Therefore, in the allocation step, it is necessary to judge the distance between the target area and the clustering center. The distance from all the target areas in the cluster to the relay station (denoted by ) should satisfy .

In this paper, the minimum mean square error (MSE) is used as the optimal index, which is one of the most common criteria in the clustering algorithm. The details of the -means++ algorithm in this paper are introduced as follows:

Given target areas, the goal is to divide them into clusters, each containing to target areas, where and represent the number of those target areas that are far from the centroid and cannot be assigned to this cluster. The problem can be expressed as follows: where denotes the location of the -th target area; represents the center of the -th cluster; denotes the number of target areas in the -th cluster; represents the distance between and .

Let denote the partition matrix. indicates that the -th target area belongs to the -th cluster, and indicates that -th target area does not belong to the -th cluster. Here, . Then, the problem in (2) can be described as follows: where is the evaluation index of the cluster algorithm. In the allocation step, the goal is to minimize by changing while is fixed. The concrete strategy is where , , and are slack variables used to eliminate inequalities, which are integer numbers, and . An example of the process in the allocation step is shown in Figure 3, where five target areas need to be assigned to two clusters and the distances between all the target areas and the center of the two clusters are all less than . The weight of each edge is denoted by , indicating the distance from the -th target area to the center of the -th cluster.

In the update step, when all the target areas are assigned, the new centers have to be updated to minimize by changing , while keeping fixed. The concrete strategy is as follows:

The proposed clustering method for the position selection of the relay stations can be described in detail as Algorithm 1 (see Figure 4).

3.2. The Cuckoo Search Algorithm in the Second Stage

The goal of the second stage is to plan a path for the UGV. After the first stage, the locations of the relay stations for the UGV are obtained. Then, the UGV needs to arrive at these relay stations in a certain order, which can be seen as a traveling salesman problem (TSP). Here, the distance between two stations in our model is symmetric, namely, it is a symmetric TSP. To deal with this problem, an improved discrete cuckoo search algorithm is proposed.

The main reason of using the cuckoo search (CS) algorithm is that it has the advantages of simple structure, less parameters, and strong ability to jump out of local optimization [24, 25]. The basic ideas of it are as follows: (1) Each cuckoo produces only one egg at a time and randomly selects a parasitic nest to hatch it. (2) In the randomly selected parasitic nest, the best nest will be retained to the next generation. (3) The number of available parasitic nests is fixed, and the probability that the host of the parasitic nest will find outgoing eggs is denoted by . The cuckoo finds the bird’s nest path and updates its location by: where means point-to-point multiplication; and represents step size and obeys Levy distribution [26]; means the Levy flight (where is a parameter between 1 and 3), which essentially provides a random walk for the UGV and is used to randomly change the position of the UGV.

Based on the general cuckoo search (CS) algorithm, there are some shortcomings, such as precociousness and so on. To make the path planning method based on the CS algorithm more efficient, some improvements are presented, which are introduced in details as follows: (1)Initialization Route. To speed up the convergence of the algorithm and prevent the algorithm from falling into precociousness, the methods of the Greedy strategy and Tabu Search algorithm are used to initialize the path. First of all, the initial position is given. Then, the taboo table is established, and the next relay station is selected according to the probability , which is calculated bywhere represents the location of the relay station visited at the -th order; represents the distance between and .

The relay station that has been visited is added to the taboo table and will be never accessed again. This process will continue till the taboo table contains all the relay stations. (2)The Fitness Function. In this paper, the fitness function is(3)The Adaptive Adjustment Operator. In this paper, an adaptive local adjustment operator combined with the 2-opt operator is proposed as the bird’s nest movement strategy, which is used to replace the Levy flight in (6) of the general cuckoo search algorithm. The process of this adaptive adjustment operator is as follows:

Step 1. Suppose is a path for relay stations, and the path is divided into part1 () and part2 (), where is the function to get the remainder value; is an integer number, and .

Step 2. If , two relay stations are randomly taken out from different segments in part1. If , two relay stations are randomly taken out from different segments in part1 and one relay station is selected in part2. If , two relay stations are randomly taken out from different segments in part1 and two relay stations are selected in part2.

Step 3. Determine whether to exchange the relay stations selected in Step 2 according to the value of . If , exchange the two relay stations, otherwise, do not change. Here, is a random number between (0,1), and is a parameter increasing with the number of iterations , namely, where is the maximum number of iterations.

An example of the process of the adaptive local adjustment operator is shown in Figure 5.

3.3. The Improved Dynamic Programming Algorithm in the Third Stage

The main goal of the third stage is to find an optimization scheme that can assign complex tasks to the right UAV, with the aim of minimizing energy consumption and balancing the workload while completing the task. To complete this task efficiently, an improved task assignment method is proposed, which includes two steps. One step is the task allocation between UAV groups, and the second one is the task allocation within the UAV group. (1)The first step

The first step is the task allocation between UAV groups, the purpose of which is to assign complex tasks to suitable UAV groups with minimum energy consumption under the premise of meeting the needs of capability. In this study, a strategy to calculate the number of UAV groups is proposed, namely, where is the number of capability types of the UAVs, and is an adjustment parameter to make sure of the success of the tasks. Then, the UAVs are divided into groups; the detailed processes are as follows: (1) The UAVs with all kinds of capabilities are distributed into each group evenly. (2) If there are still some redundant UAVs, they will be distributed randomly among groups.

After the grouping process of UAVs is finished, an improved dynamic programming algorithm (DP) is used to assign tasks to these UAV groups. The core idea of the task assignment algorithm based on DP is to decompose the original problem into subproblems and find the best scheme according to the value of the energy consumption. The processes are carried out as follows.

Firstly, the sum of the capabilities of each group is calculated by where is the capability of the -th UAV in the group . If the sum of the capabilities of the group is greater than that required by the task (denoted by ), this means that the group can be arranged the task . Then, all the tasks are randomly assigned to those groups that can be arrange them. On this basis, the tasks are assigned between groups based on the DP algorithm. The detailed processes are as follows: (1) Set a parameter for each task and the initial value of is set as 0. (2) Calculate the value of by the iterative method till the minimum is found for each task. In this study, the value of consists of two parts, namely, where represents the energy consumption of executing the current task by the group ; denotes the energy consumption of performing the subsequent tasks.

An example of this process is shown in Figure 6, where 6 UAVs are divided into 2 groups, and 7 tasks are required to be completed. In this example, two different capabilities are required by the task (denoted by different colors). First, an initial random assignment is made between the two UAV groups to meet the task capability requirements, namely, and . Based on the optimal strategy, the final assignment is and . (2)The second step

In the task studied in this paper, the target areas have multiple tasks, each of which requires specific capabilities. After the task assignment between UAV groups, the complex tasks are assigned to multiple UAV groups. Then, the members of the UAV groups should be assigned to their own tasks based on their capabilities. In this step, not only the energy consumption should be considered but also the load balance should be taken into account.

However, the energy consumption and the load balance are to some extent a pair of contradictory parameters. Therefore, in order to achieve global optimization, some trade-offs must be made between the two indexes. To deal with this problem, a probability-based method for the UAV selection is proposed. The selection probability is determined by the energy consumption and the load, which is defined as follows: where is the total energy consumption of the UAV for the task requiring the capability ; is the remaining energy of the ; means the UAVs that provide the capability in the group ; is the parameter used to balance the energy consumption and the load balance. A larger value of indicates more importance of the load balance.

In summary, the allocation algorithm in this paper is shown in Algorithm 2 (see Figure 7).

4. Experiments and Analysis

To verify the feasibility and effectiveness of the proposed method, some simulation experiments were carried out by a computer with 8 GB RAM and i5-3230M 2.5 GHz CPU using the MATLAB as a platform. To simplify the realization, the assumptions in this study are as follows: (1) The UAVs, UGV, and targets are assumed as points without any shapes; (2) There are not any obstacles in the environment. (3) The speed of the UAV is 120 km/h, and the speed of the UGV is 25 km/h. (4) The energy consumption of the UAV is , while the UGV’s is . (5) The energy consumption for the UAV to complete a task is . Some assumptions of the UAV are near to the real commercial drones, such as the UAV products of Dajiang Innovations Technology Co., Ltd. [27, 28]. All the parameters in these simulations are the same, which are listed in Table 1.

In this study, the following indexes are used to evaluate the performance of the proposed method in completing the task. (1)The time used in the task (). Because the running time of the algorithm is negligible compared with the time needed to conduct the task by the UAV/UGV hybrid system, can be calculated as follows:where is the time taken on the path of the UGV; and is the time spent in conducting the tasks of the UAVs. (2)The energy consumed in the task (). The total energy consumed in the task includes the energy consumption of the UAVs in conducting the tasks and the energy consumption of the UGV on the road, namely,(3)To evaluate the comprehensiveness of the proposed method in different tasks, the average energy consumption and time used in each relay station are calculated and analyzed

4.1. Simulation of Simple Cooperative Task

To test the performance of the proposed approach, this experiment is carried out, where the size of the environment is . In this simulation experiment, there are 400 target areas and 9 UAVs. The initial positions of these target areas are shown in Figure 8(a). In this experiment, the maximum number of tasks for one target area is set as 1, and the types of the capabilities required of the task are 3. The UAVs are divided into 3 groups based on the proposed division strategy, and the parameter is set as 0.22 (namely, ). The clustering and path planning results are shown in Figure 8(b). The time and the energy consumption in this cooperative control task are listed in Table 2. To evaluate the performance of UAVs based on the proposed method in the process of cooperative control task, 10 relay stations are selected randomly and the time and energy consumption of UAVs in these 10 stations are given out in Table 3.

Remark 1. In Tables 2 and 4, the values of Average1 and Average2 are calculated by:

The results in Figure 8(b) show that the proposed algorithm can find the optimal locations of the relay stations and complete the path planning task for the UGV efficiently. The results in Table 2 show that most of the time and energy are spent on the process of performing the task by UAVs. In this experiment, the maximum energy consumption of UAV is close to the average consumption per UAV (see Table 3), which means the load balance among UAVs in the system is good. Based on the proposed method for the UAV/UGV hybrid system, the entire task can be completed quickly and efficiently.

4.2. Simulation of Complex Cooperative Task

To further verify the effectiveness of the proposed method, the second simulation experiment is conducted, where the size of the environment is and the number of target areas is 4000. In this task, there are 30 UAVs, which have 5 types of capabilities to complete the tasks. In this complex cooperative control task, the maximum number of tasks of one target area is set as 10, and the types of the capabilities required by the task are set as 5. Based on the proposed method, the UAVs are divided into 5 groups, and the value of the parameter is set as 0.3 (namely, ). The initial positions of these target areas are shown in Figure 9(a). The clustering and path planning results are shown in Figure 9(b). The time and the energy consumption in this cooperative control task are listed in Table 4. The 10 relay stations of the 200 relay stations are also selected as representatives to observe the time and energy consumption of UAVs in the detailed task execution process (see Table 5).

The results in Figure 9(b) show that the proposed algorithm can find the optimal locations for all the target areas and avoid the occurrence of local optimization when the number of target areas is large. The results in Table 4 show that the energy consumption increases obviously in one relay station because the number of tasks of one target area increases. However, the hybrid UAV/UGV system can complete the complex cooperative task efficiently, the average time needed for one relay station is 0.164 (), which is almost equal to the time needed in the first simulation (see Tables 3 and 5). The main reason is that the proposed method considers the balance between the energy consumption and the load of UAVs, which is very important in the complex cooperative control task. The standard deviations of the UAVs’ energy consumption and time consumption in different relay stations are small compared with the average values (see Table 5), which means the task assignment of the proposed method is efficient, and the stability of the proposed method is good.

5. Discussions

The simulation results in Section 4 show the efficiency of the proposed approach in various cooperative search tasks under various situations. In this section, some key parameters and the performance of the proposed method are discussed.

The first important parameter of the proposed method is of the clustering algorithm in the first stage. To discuss the effect of the value of , different experiments are conducted on the two tasks tested in Section 4, with different . The experimental results are listed in Table 6. The results in Table 6 show that the value of total energy consumption will decrease with the increase of the value of for the same tasks. The main reason is that the energy utilization of UGV is higher than that of the UAV. However, the total time to complete the task will decrease firstly and then will increase with , and there is an inflection point (see the last column in Table 6). So to balance the cost and time, it is better to set the value of at about 30 for the task1 and 200 for the task2. Based on the proposed empirical equation, it is easy to obtain a proper value of , which is near to the needed optimal value.

Another important parameter is of the task assignment algorithm in the third stage. To discuss the effects of the parameter , some simulations are conducted with different values of . In these simulations, the tasks in Section 4 are used. The results are shown in Figure 10. Figures 10(a) and 10(b) show how the energy consumption and its standard deviation change with the parameter in the task1 and task2, respectively. The results in Figure 10 show that the energy consumption in the task gradually increases with the gradual increase of value, while the standard deviation of energy consumption of all UAVs decreases with the increase of (which means the load balancing becomes better when the value of increases). So to consider the load balancing and the energy consumption comprehensively, the value of in this paper is set as 0.4.

To show the efficiency of the proposed method in the task assignment of the third stage, some comparison simulations are conducted among the proposed method, the particle swarm optimization (PSO) method, and the greedy method (Greedy). To be comparable, the methods used in the previous two stages are the same as the proposed methods, except the method used in the third stage. The task in this experiment is the same one in Section 4.2. The data from 10 randomly selected relay stations are analyzed for comparison. The results are shown in Table 7. The results in Table 7 show that the proposed method can maintain the lowest energy consumption and complete the task faster than another two methods. And the results show that the greedy method has the worst performance. This is due to the fact that the core of the greedy-based algorithm is to find the local optimal solution without considering the global situation, which will lead to the deterioration of the overall performance. The PSO method can get a better performance than the greedy method, but it does not take into account the requirements of load balancing between UAVs, which will lead to the overuse of some UAVs, thus increasing the overall task completion time (see Table 7).

6. Conclusion

The problem of cooperative control for the complex multitasks in a large environment using unmanned aerial vehicle and unmanned ground vehicle (UAV/UGV) system is studied in this paper, and a novel integrated method is proposed. The goal is to minimize the time and energy consumption taken by the system to complete the entire task. In the proposed method, the cooperative control problem for the complex multitasks is divided into three stages. Firstly, an adaptive clustering method is proposed to determine the optimal positions of the relay stations for the UGV. Then, the improved cuckoo search algorithm is used to solve the shortest path problem of the UGV. Finally, a grouping method is proposed to solve the multitask assignment problem of UAVs. The experimental results show that the proposed method can finish the cooperative task efficiently. Some discussions are conducted on the key parameter selection and the performance of the proposed method compared with other classical task assignment methods. In the future work, the practical experiment of the hybrid UAV/UGV system cooperation will be carried out to verify the actual performance of this method. In addition, some more efficient bioinspired methods will be studied to satisfy the requirements of some more complex UAV/UGV cooperative tasks, such as the cooperation task in the environment with obstacles.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declared that they have no conflicts of interest to this work.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (61873086), and the Fundamental Research Funds for the Central Universities (2018B23214).