Abstract
Hybrid wireless mesh networks are suitable for complex environment communication in coal mine. Mesh clients with application service and routing function in hybrid wireless mesh networks can form a highly robust hybrid network with mesh routers. The processes of nutrient flux transfer and path choice in Physarum network are similar to data transmission and routing selection in hybrid wireless mesh networks. In this article, we use Physarum-inspired autonomic optimization model to design a Physarum-inspired multi-parameter adaptive routing protocol to improve the service quality of coal mine hybrid wireless mesh networks. Physarum-inspired multi-parameter adaptive routing protocol has achieved distributed routing decision by drawing the hybrid wireless mesh network parameters into Poisson’s equation of Physarum-inspired autonomic optimization model to measure the quality of link and implements two adjustment strategies to make the protocol more adaptive. The resource-dependent adjustment, which considers the irreversible energy consumption and recoverable buffer occupation, makes the energy consumption problem prominent when there is a lack of energy. The position-dependent adjustment makes routing decision efficient according to the load of different positions, which is caused by many-to-one data transmission model in coal mine. Based on NS2, simulation experiments are performed to evaluate the performance of Physarum-inspired multi-parameter adaptive routing protocol, and the results are compared with those of ad hoc on-demand distance vector, HOPNET, ANT-DSR, and Physarum-inspired routing protocols. The experimental results show that the route path selected by Physarum-inspired multi-parameter adaptive routing protocol is better than those selected by the other four protocols in the performance of average end-to-end delay and delivery ratio. The balance of energy consumption and network load is achieved and the network lifetime is effectively prolonged when using Physarum-inspired multi-parameter adaptive routing protocol.
Introduction
The coal mine is a mobile and complex working environment due to coal mining and working-face driving, and wired communication networks cannot reach these special areas, which makes it easy for blind monitoring areas to form. 1 Therefore, it is urgent to pay more attention to the theory and technology of new wireless communication systems and construct a new type of broadband wireless network to adapt to the special underground working environment. Wireless mesh networks (WMNs), as a new kind of broadband wireless network, can effectively solve the above problem, which the wired network and traditional wireless network cannot solve. 2
According to the type and organization mode of the nodes, WMNs can be divided into three categories: backbone WMNs, client WMNs, and hybrid WMNs. 3 In backbone WMNs, mesh clients (MCs) cannot directly communicate with each other and must connect to mesh routers (MRs) to achieve broadband network access. Client WMNs are composed of MCs, which need to support data routing and network configuration, as well as provide end-user applications to customers. Hybrid WMNs are composed of MCs and MRs. In this kind of networks, MCs can perform the routing functionality as MRs, which results in hybrid wireless mesh networks (HWMNs) with strong self-organization and self-healing abilities, as well as high robustness. Therefore, HWMNs can be used to transmit information in special areas and on special occasions in an underground mine.
In HWMNs, MCs are powered by batteries, which are energy limited. The energy supply of MRs is sufficient because of the wired power supply. In addition, the computing power and buffer space of MRs are often stronger than those of MCs. Therefore, we should take full account of the characteristics of different types of nodes in routing decisions. The effective working times of nodes directly affect the network lifetime. The energy balance should be considered in the design of the routing protocol, as well as the load balance. Researchers have put forward some methods to solve the above problems. These methods are based on the performance parameters of nodes, such as energy and buffer, which are used as the routing decision metric by adding appropriate proportion adjustment factors. However, these factors are given fixed values, resulting in poor adaptability to the changing network state.
In recent years, researchers have begun to study the bionic algorithm, and the Physarum-inspired autonomic optimization model (PAO) has become a hot research topic. The process of nutrient delivery and the path choice in the Physarum network are very similar to data transmission and route selection in HWMNs. Thus, the characteristics of the Physarum network can inspire us to apply such a model to routing decision in HWMNs. However, the PAO is not appropriate for use in HWMNs directly because it is calculated based on the global information and has high complexity and overhead.
The algorithm proposed in this article views the tubes of the Physarum as network links in HWMNs. Parameters such as energy consumption, buffer occupation, and hop count to the gateway (node depth) are used as the routing metrics. In HWMNs, MCs are powered by limited energy and the energy consumption process is usually irreversible. The buffer occupation of nodes is a kind of “elastic” network parameter, changing with the network load. Therefore, the proportion adjustment factors of energy consumption and buffer occupation should be able to adapt to the different performance characteristics of these two parameters. In addition, for the nodes at the edge of the network, the distance to the gateway should be prioritized so the routing process will converge as quickly as possible, while more consideration should be given to load balance in the routing decision for the nodes close to the gateway with high load pressure. Based on the PAO and the analysis of the parameters’ characteristics, energy and load balance is implemented in the algorithm with low overhead and network state adaptation as well.
The remainder of this article is organized as follows. Section “Related work” reviews the related work. Section “HWMNs in coal mines” introduces HWMNs in the coal mine. The Physarum-inspired multi-parameter adaptive routing (PIMAR) protocol for coal mine HWMNs is proposed and described in detail in section “PIMAR protocol.” Section “Implementation of PIMAR” gives the implementation process of PIMAR. Section “Simulation experiment” presents the results of the simulation experiments. Finally, the conclusions of the study are presented in section “Conclusion and future studies.”
Related work
With the development of the WMN, the research on routing protocols is also advancing. There are some classical routing protocols currently suitable for WMNs, such as destination-sequenced distance vector (DSDV), ad hoc on-demand distance vector (AODV), and hybrid wireless mesh protocol (HWMP).4–6 In a WMN, MC energy, computing power, and buffer are often limited, which affects the network performance directly. Researchers have proposed many routing protocols based on these performance parameters to solve the problem of load and energy balance to prolong the network lifetime and improve the network performance.
A power-aware routing protocol, which implements energy balance in the network, was proposed in Toh. 7 A method to measure the load level of the media access control (MAC) layer was introduced, and load balance was implemented in Kim et al. 8 Because the network performance is affected by many factors, the routing metric often contains multiple parameters. A routing decision method combining node mobility and congestion was proposed, and a congestion avoidance routing protocol was implemented in Basarkod and Manvi, 9 which can adapt to frequent changes in network topology. Jabbar et al. 10 used the idle time of the node, residual power, and queue length as the metrics to measure the link quality between nodes and implemented energy and load balance. A metric based on the expected transfer time of the network and the buffer queue of the nodes was proposed to implement load balance in Ma and Denko. 11 It is necessary to introduce the parameter adjustment factors to adjust the proportions of different parameters in the multi-parameter routing metric. However, the factors often take fixed values and the adaptability of the fixed factors is poor due to the frequently changing network state.
In recent years, many fields have drawn attention to bionics and intelligent algorithms, such as the ant colony algorithm and swarm intelligence algorithm, used in research on topology and routing in WMNs. The ant colony optimization (ACO) algorithm is based on ant colony foraging to find the optimal path. 12 Aissani et al. 13 proposed ant dynamic source routing (ANT-DSR), a proactive routing protocol based on ACO. ANT-DSR uses unicast ant agents to find routes and always maintains a certain number of ants for routing optimization. This approach can effectively shorten the routing time but leads to relatively large network overhead. Wang et al. 14 proposed a hybrid ant colony algorithm called hybrid ant colony optimization routing algorithm for mobile ad hoc network (HOPNET) which contains four types of ants: intra-zone forward ants, intra-zone backward ants, inter-zone forward ants, and inter-zone backward ants. However, HOPNET is difficult to adapt to a larger network because of the lack of good scalability.
With the development of bionics and some scientific discoveries, cellular computation has become an important branch of biological heuristic algorithms and has attracted increasing attention to the calculation of the Physarum. The Physarum is an ancient protoplasmic slime with a tubular network to absorb nutrition. In the process of finding and absorbing nutrients, the Physarum can retract to an efficient tubular network. Nakagaki et al. 15 investigated the maze-solving behavior of the Physarum and observed that the Physarum can find the shortest path between a food sources at the maze exit and the entrance without any control.
Tero et al. 16 established a geographical model of Tokyo and the area around the main city using oatmeal and cultured Physarum on it. A biological network that connects all food sources had formed by the end of the experiment and the network was comparable to the real railway network of Tokyo in the aspects of efficiency, reliability, and cost. PAO laid the foundation for future research in Tero et al. 17 A method using PAO to solve the Steiner tree problem was presented in Liu et al. 18 New ideas for the construction of more efficient transportation networks by applying PAO to traffic network optimization were provided in Marwan 19 and Tsompanas et al. 20 Some methods of optimizing a logistics distribution and supply chain network based on the model were presented in Chen et al. 21 and Zhang et al. 22 A method of solving the shortest path problem was proposed in Lu. 23 Ma and coworkers applied PAO to the distribution of wireless sensor networks to optimize the detection range. 24 PAO is an adaptive algorithm that can adjust to changes in the environment to obtain a more stable network. However, the process of finding the best path with iterative calculations in PAO has high complexity and overhead. Therefore, PAO is not suitable for direct application in WMNs, which have limited resources.
In order to apply PAO to routing decision, we need to migrate the model to WMNs. Zhang et al.25,26 applied PAO to location-aware routing protocol and proposed a protocol called P-iRP, which considered the distance, energy, and location of the next hop. Zhang et al. 27 provided a novel bio-inspired hybrid routing protocol based on zone routing framework, ACO, and PAO, and the network performance has been improved significantly. PAO has also been used in energy-harvesting wireless sensor networks, and three distributed routing algorithms with low algorithm complexity have been proposed in Tang et al. 28
We propose a new HWMN routing algorithm, called PIMAR, based on PAO. PIMAR uses the information of local nodes, such as energy consumption, buffer occupation, and node depth, in the metric for route selection. The algorithm has low complexity, which ensures the feasibility in the HWMNs. In addition, the algorithm incorporates the proposed strategy of adaptive parameter adjustment based on the node resources and network topology, which increases the adaptability of routing decisions to changes in the network state. Moreover, the algorithm implements energy and load balance and has a high throughput and a small transmission delay.
HWMNs in coal mines
Network structure
There are many narrow tunnels in coal mines. With coal mining and tunnel excavation, the tunnels extend gradually and form a band structure with a length of hundreds of meters, or even a thousand meters. WMNs in coal mines have a long-chain structure (backbone WMN) adapted to the long-chain structure of the tunnels. Data packets are transmitted to MRs and are forwarded to the gateway node, as shown in Figure 1.

Structure of a backbone WMN in a coal mine.
In the long-chain backbone WMN, MCs are overwhelmingly dependent on MRs, which will lead to the decline of network performance and even paralysis once failures of the MRs occur. Therefore, a network with this structure has low robustness and poor self-repairing ability. 29 In an HWMN, MCs can communicate with one another and data packets can be forwarded by MCs when MRs are broken down. The HWMN has higher reliability, as shown in Figure 2.

Routing repair of an HWMN in a coal mine.
Node type
MCs in a coal mine have limited energy, unlike in a WMN on the ground. There is no energy restriction for MRs because they are supplied by wired power. MCs are the main nodes in an HWMN, so the energy limitation problem will directly affect the network lifetime and is even more serious in mine emergency communication.30,31 MRs have higher performance capabilities than MCs in terms of buffer space and data processing speed. It is necessary to optimize the utilization ratio of different types of nodes to prolong network lifetime and improve network performance.
PIMAR protocol
Parameter mapping
There are many junctions in the Physarum network, such as i and j in Figure 3. We define the tube from i to j as

Physarum network.
As mentioned in Tero et al., 17 the flux of nutrient can be expressed by the Poisson equation
where
The shortest path can be chosen by PAO. Therefore, the length
The process of data flowing from one node to another in an HWMN is similar to the process in fluid mechanics. The direction of movement of the fluid is determined by the pressure difference
Resource-dependent adjustment strategy for pressure difference
The factors such as energy consumption and buffer occupation with different characteristics are indispensable when choosing the next-hop node. MCs are powered by batteries and the energy consumption is irreversible. The buffer occupation is influenced by the network load, which is an “elastic” factor. The energy consumption ratio of the next-hop node is denoted as
where
The parameters of energy consumption and buffer occupation are the measures of node state: the more energy and buffer that remain, the better the node state. The correspondence between the pressure concept of PAO and the node state in an HWMN can be expressed as follows: the better the state of the next-hop node, the greater the pressure difference between this node and the next node. If the residual energy and buffer are viewed as the driving forces of the data flow, the energy consumption and buffer occupation can be regarded as the resisting forces. We define the pressure difference between i and j as
where M is a constant.
where
If the values of
The value of
On the other hand, when the residual energy of nodes is sufficient, the buffer space should be given more weights. We set
When the energy consumption is less, the value of
Position-dependent adjustment strategy for pressure difference
The load of the node is related to the position in the coal mine network because of its many-to-one data transmission mode. The nodes close to the gateway carry higher loads than nodes that are farther away. To ensure that the route to the gateway converges as quickly as possible, we should focus on the hop count when the node is far from the gateway. A node close to the gateway will carry a higher load and more attention should be paid to its energy and buffer space.
To achieve the position-related adjustment strategy, we have introduced hop factor in our model. In the initial stage of the network, a broadcast message will be sent by the gateway to all the nodes to determine the node depth of the network. Next, the nodes can confirm the hop resistance to neighbors.
The hop resistance
where
After adding the hop factor,
where the value of
The values of

Values of hop resistance factor.
As shown in Figure 4, if we take a larger value of
According to the value range of each variable, we can determine that the data range of
Thus, the final Poisson equation for HWMNs is defined as follows
Routing metric of PIMAR
The flux
There are limitations in considering the performance of only the next-hop node when making routing decisions; therefore, PIMAR uses regional flow information of nodes as a routing metric.
Suppose
The method of calculating the regional flow value is defined as follows
where
Implementation of PIMAR
Node depth initialization
In the initial stage of the network, the initial packet (INIT_MSG) is broadcasted by the gateway node. The structure of INIT_MSG is shown in Table 1.
Structure of INIT_MSG.
TTL: time to live.
Packet Type represents the type flag, the value of which is INIT_MSG. Hop Count refers to the distance from the current node to the gateway node. TTL is the time to live of INIT_MSG, which is set to 1. Source IP Address is the IP address of the node that sends INIT_MSG and Destination IP Address is always set to IP_BROADCACT, which indicates that INIT_MSG is a broadcast message.
The process of determining node depth is shown in Figure 5.

Process of determining node depth.
In Figure 5, steps 1 and 2 are executed only once and step 3 is executed repeatedly by each MR and MC until the node depth is confirmed.
Maintenance of neighbor table
Neighboring nodes need to exchange information regularly to obtain the latest states of one another, so we defined a packet (NB_HELLO) as shown in Table 2.
Structure of NB_HELLO.
TTL: time to live.
Depth is the node depth, which is confirmed at the initialization state. Energy Consumption Ratio and Queue Occupancy Ratio are properties of the neighbor that sends NB_HELLO. Neighbor Average Flow represents the average flow of all j’s neighbors except i, as defined in formula (13). NB_HELLO is also a broadcast message, so Destination IP Addresses is set to IP_BROADCAST.
We set the message exchange period to 1 s in the stable network state, in which the energy consumption and node depth cannot change rapidly to minimize the network overhead caused by information maintenance. However, if the state of the node is instable, meaning that the variation range of resource potential is greater than 5%, NB_HELLO will be broadcasted to notify the neighbors immediately.
PIMAR needs parameters such as energy consumption, buffer occupation, and node depth to make routing decisions, so each node in the HWMN needs to record this information in the neighbor table. In addition, the adjustment factors,
Routing decision
Through the maintenance process of neighbors’ information, the node can obtain the parameters of its neighbors. When a data packet is received, the node will make a routing decision and select a neighbor to which to forward the packet based on the information in the neighbor table. According to the principle of PIMAR, the neighboring node with the largest
Node resource parameters vary with time in a running network, which may result in a routing loop. We let nodes record packet information such as the packet id, source IP address, and destination, which can be used to determine packet identity solely. The information is stored for a short time after the node forwards the data packet and is used to judge whether there is a routing loop. When the node receives the same packet again, there is a routing loop and the node will select another node from the neighbor table as the next-hop node to retransfer the packet.
The process of routing decision is shown in Figure 6.

Process of routing decision.
Algorithm complexity analysis
PIMAR is a kind of context routing protocol, according to the distributed process of routing decision which depends on the neighbors’ current state. We suppose that the time of nodes to calculate
When the nodes calculate
Node depth is confirmed in the initial stage of the network. INIT_MSG is sent from gateway, and the message will be discarded or forwarded according to the initialization process (Figure 5) by the other nodes. Each INIT_MSG has the same message length and nodes cost the same time to handle it. The arrival time of INIT_MSG is determined by hop counts. We can infer that the later arrival message has passed more intermediate nodes and the Hop Count obtained in the first arrived INIT_MSG is the minimum node depth. The number of messages sent in this process is
However, the wireless transmission may be delayed in practical applications. The node may receive multiple initialization messages with smaller network depth, resulting in redundant initialization information in the network. In order to reduce redundant information as much as possible, nodes will delay a little time to send its node depth to its neighbors after receiving INIT_MSG. We suppose that
Example of PIMAR
For a more detailed introduction of PIMAR, we assume the network topology shown in Figure 7. Node g is a gateway node, b, e, and f are MRs, and the other nodes are MCs. The node depth is marked at the side of each node in the figure. In the initial state of the network, there is no energy consumption and the buffer occupation of each node is zero.

Network topology.
In this example, node m is assumed to be the data source. The process of routing decision is shown in Table 3.
Routing decision process.
DP: node depth.
Routing result of each step.
The routing result of each step is shown in bold text with “*.” From Table 3, the first hop chosen by m is f and the next one chosen by f is b. From the topology, we know that b is connected with the gateway node g directly, so the data packet will be delivered to the gateway node without calculation. The final path chosen by PIMAR is m → f → b → g, as shown in Figure 8. According to the result of path selection, PIMAR can choose the link with more routing nodes, which can ensure the stability of the link.

Route selection of initial conditions.
While the network is running, the energy and buffer of the nodes are consumed. The process of routing decision is shown in Table 4 after a period of time.
Routing decision process.
DP: node depth.
Routing result of each step.
In Table 4, nodes b, e, f, and j are MRs, which have wired power supplies, so the energy consumptions are all 0. We suppose that b and f, which are in m → f → b → g, have high buffer occupations. The final routing path is m → j → e → a → g, as shown in Figure 9.

Route selection after node resource change.
As shown in Figure 9, PIMAR avoids the path with a heavy load and implements the load and energy balance.
Simulation experiment
Simulation environment and experimental parameters
Based on the NS2 simulation platform, we test the performance of PIMAR and the results such as average end-to-end delay, delivery ratio, the usage of the nodes, and network lifetime are compared to those of AODV, ANT-DSR, 13 HOPNET, 14 and P-iRP. 26 The parameters of the simulation network are listed in Table 5.
Simulation parameters.
CBR: constant bitrate; MC: mesh client.
In the simulation, we set the network in a long and narrow area of 1500 × 8 m2 to simulate the tunnel in an underground mine. We have deployed 15 MRs, 135 MCs, and 1 gateway in the tunnel. We have simulated the network performance using the parameters, such as average end-to-end delay, delivery ratio, node utilization, and network lifetime, under different load pressures by increasing the number of constant bitrate (CBR) data flows.
Performance analysis
Average end-to-end delay
The trends of average end-to-end delay for four different protocols are shown in Figure 10. AODV has the highest average end-to-end delay as well as the fastest growth rate when the number of CBR flows is increasing. The average end-to-end delays of HOPNET and ANT-DSR are smaller than that of AODV. The delay values for P-iRP and PIMAR are smaller than those for HOPNET and ANT-DSR, and PIMAR has the smallest one. This is because AODV is an on-demand routing protocol that requires routing requests before sending data, which may cause time delay. In addition, AODV always selects the path with the minimum number of hops to send data. When the data flow increases, network congestion will occur and the average end-to-end delay will increase. HOPNET and ANT-DSR are proactive routing protocols that always maintain a certain number of ants for routing optimization, so the route selection will be faster than with AODV. However, when the network load is heavy, the overheads of HOPNET and ANT-DSR will become large due to routing maintenance, which requires ants for each data flow. The large overhead will also cause high end-to-end delay. P-iRP does not need much maintenance packets compared with HOPNET and ANT-DSR, because it is location aware. So the delay can be less than HOPNET and ANT-DSR when the data flows become larger. PIMAR can avoid the congested areas effectively and balance the network load using two parameter adjustments. The adjustments consider the buffer occupation of each node as well as its characteristic changing with the node depth to calculate regional flow, which can ensure that nodes with heavy loads are less likely to be chosen. In addition, PIMAR uses a lightweight neighbor information maintenance method, so that its routing maintenance overhead can be maintained at a low level. Based on these two points, PIMAR can achieve low end-to-end delay. Because P-iRP has less consideration to deal with network congestion, the delay of PIMAR is lower than that of P-iRP.

Average end-to-end delay.
Delivery ratio
Figure 11 shows that PIMAR has the highest data delivery ratio, compared with P-iRP, HOPNET, and ANT-DSR, and AODV has the lowest one. Similar to the end-to-end delay, PIMAR has a congestion mechanism, which can avoid the heavy-load region to reduce the packet loss ratio. The route selection method of AODV only considers the hop count from node to gateway and the network congestion is not included. When the number of CBR flows is large, AODV will lose more packets. HOPNET and ANT-DSR adopt the route optimization mechanism of maintaining a certain number of ants all the time to obtain better routing paths, so their delivery ratios are better than that of AODV. P-iRP has lower routing maintenance cost than HOPNET or ANT-DSR, as a result of which its delivery ratio shows better performance.

Delivery ratio.
Network lifetime
First, we define the time at which the first node dies as the network lifetime. It can be seen from Figure 12 that the network lifetimes of AODV, HOPNET, ANT-DSR, and P-iRP are shorter than that of PIMAR. This is because PIMAR adopts the adaptive parameter adjustment strategy that considers the irreversible energy consumption characteristic of MCs, which are different from those of MRs. In the process of running the network, energy problem is prominent for MCs; therefore, PIMAR adjusts the weight of the energy factor with the change in energy consumption: the more energy the node has consumed, the bigger the value of

Network lifetime.
Residual energy of MCs
The residual energy of each node is collected at the time at which the first node dies in the network. We adopt the variance of residual energy to measure use balance of the nodes, and analyze the residual energy balance of the nodes under light (CBR = 16), medium (CBR = 40), and heavy (CBR = 72) loads. If the value of the variance is smaller, the use of the nodes is more balanced. It can be seen from Figure 13 that the nodes of AODV are the most unbalanced, followed by those of the ANT-DSR, HOPNET, and P-iRP protocols, and the nodes of the PIMAR protocol are the most balanced under different network loads. From this result, we conclude that the adjustment strategies based on node resources and node depth can effectively implement the energy and load balance.

Variance of residual energy percentage.
Conclusion and future studies
The HWMN is suitable for complex environments such as coal mines because of its high robustness. Many types of routing protocols have been proposed in recent years, such as bionic algorithms. Cellular computation, an important branch of biological heuristic algorithms, has attracted increasing attention as a hot research topic. There are similarities between Physarum network and HWMNs. This article proposes a PIMAR protocol to implement energy and load balance in HWMNs in coal mines. Based on PAO, we use the Poisson equation combined with the parameters of HWMNs, such as energy consumption, buffer occupation, and node depth, to calculate the network flow, which can be used to measure the quality of a link. We analyze different characteristics of energy and buffer as well as changes in the network load with node depth and propose two adaptive parameter adjustment strategies. The energy consumption of MCs is irreversible, while the buffer occupation is “elastic.” The network loads of nodes near the gateway are heavier than those of the nodes that are farther away because of the many-to-one data transmission mode in coal mines. Based on the analysis, we adjust the weights of those parameters dynamically to prioritize the energy problem and avoid regions with nodes in poor states in the network running process. In addition, PIMAR uses local information to calculate regional flow to make routing decisions, which greatly reduces the complexity of the algorithm, lowers the network cost of the protocol, and improves the scalability. On the NS2 simulation platform, we test the performance of PIMAR and the results show that PIMAR not only has better end-to-end delay and delivery ratio compared to AODV, ANT-DSR, HOPNET, and P-iRP but also prolongs the life of the network and implements energy and load balance.
Most of the nodes in HWMNs in coal mines are static, but there are also mobile nodes, such as the terminals for the personnel positioning system or the wireless sensors on the transport equipment. PIMAR does not yet consider the mobility of the nodes, which can cause frequent network topology changes and increase the difficulty of solving problems. The overhead of network topology information maintenance and ensuring stable links are the problems in the mobile network environment. In the future, we will pay more attention to node mobility and design our routing protocol to be more suitable for a network with unstable topology.
Footnotes
Handling Editor: Wei Yu
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The financial support for this work provided by the Fundamental Research Funds for the Central Universities (No. 2017XKQY077) is gratefully acknowledged.
