Minimum spanning tree-based digital terrain model detection from light detection and ranging points

Abstract

LiDAR is widely used to collect point data-set representing the scanning scenes. Several products can be extracted from raw LiDAR data. DTM is among the most important. For traditional DTM detection approaches such as TIN-based and parametric methods, a negative outlier detector is employed beforehand since negative outliers will drive these methods to converge into an erroneous terrain surface. However, automatic outlier removal remains challenging in related fields. To overcome the obstacle, this article proposes minimum spanning tree-based detection method. An undirected graph is constructed for LiDAR points on basis of delauney rule. Edges are weighted by absolute slope value. Then, greedy algorithm is taken to generate MSF for the undirected graph. Edges with steep slopes are omitted in this process, which is how more than one MST arises. The DTM is generally the MST with largest area. Cases validate the effectiveness.

Keywords:

• LiDAR
• DTM detection
• minimum spanning tree
• outliers

1 Introduction and motivation

LiDAR (Light Detection and Ranging) can directly acquire massive 3D point data-set representing the scanning scenes. Compared with photography measurement, LiDAR provides a flexible and efficient approach. DTM (Digital Terrain Model) that represents the bare earth is one of the most important products detected from the LiDAR raw data-set.[1, 2] DTM detection, also called filtering, is to separate the terrain points from non-terrain points. The related methods can be classified as morphological, statistical, parametric and TIN (Triangular Irregular Network)-based method.

Borrowing a series of morphological opening and closing operations from image-processing community, the morphological approach on 3D point data-set is firstly proposed in [3]. Slope-based filtering method in [4] is closely related to the erosion operation in morphology. A dual-rank filter has also been used to manage DTM generation, with an evolution of the structural element size and adaptive threshold.[5]

Skewness balancing in its basic form exploits the global statistics of a LiDAR tile in order to separate object and terrains.[6] This algorithm assumes that object points are located above terrain, which is not true in mountainous areas. Bao improved skewness balancing by incorporating the measure of kurtosis, where the author exploited the different changes of both statistical measures with respect to vegetation and ground.[7] A framework based on skewness balancing has also been built with a prediction model, in which the unknown LiDAR tiles can be categorized as “hilly” or “moderate” terrains.[8] Höhle et al. proposed a robust statistic method by computing the normalized median absolute deviation with confidence intervals in [9].

Parametric method usually starts by approximating the surface using equal weights for all points. Points with negative residuals are more likely to be on terrain. Therefore, higher weights are assigned to points with negative residuals and lower weights are assigned to points with larger positive residuals. Points with lower weights will be eliminated or not considered for further iterations. Thus, the result is that the estimated surface converges toward the points with negative residuals. Linear prediction to estimate DTM in wooded areas was introduced in [10]. An active contour-based method introduced an energy function to control the shape of the estimated surface.[11] Shan et al. provided a mathematical definition of the terrain points as a function of the slope and elevation difference between two neighbour points.[12] Based on this definition, two steps proposed are labelling and elevation adjustment. Bretar et al. generated an adaptive local geometry-based algorithm called a predictive Bayesian framework, which propagates through the point set following specific rules in order to optimize the probability of computing areas containing terrain points.[8] This method is mainly used in natural landscapes that are composed of vegetated and bare Earth landscapes.

The iterative TIN-based approach may be the most widely used DTM detection approach from LiDAR data.[13–15] It starts from a coarse triangulated surface composed by a set of seed points (usually taking the lowest points in the partial regions as seed points), densifies the TIN by adding new LiDAR points under some constraints. Then, the iteration is repeated until there is no changes in the TIN or after a specified number of iterations. This strategy has also been implemented in Terrasolid, which is a well-known software. Yet, this algorithm depends on the size of the largest non-terrain feature and is susceptible to negative outliers. Belkhouche et al. proposed an iterative TIN-based automatic interpolation filtering method, which iteratively eliminates the edges connecting terrain to object.[15] The height of object point is adjusted using an interpolation algorithm according to the estimated local slope. At the end, all the adjusted points are classified as objects and the others as terrain.

Even if these approaches already are proved to be successful in some scenarios, they suffer from some drawbacks. The morphological one remains unsatisfactory over various landscapes. Depending on the window size, large objects may remain in its final results. For algorithms of statistic, parametric and TIN-based, the obligatory operation is to remove the negative outliers since the existence of outliers below terrain surface leads to the creation of steep edges. Negative outliers will therefore be mistaken as terrain.[8, 15, 16] However, to the best of our knowledge, automatic removal of outliers from LiDAR raw data-set remains an open issue.[17, 18] Then, these implementations usually need the manual refinement of seed points or the final terrain surface in practical applications. So, the main motivation of this article is to generate a robust DTM detection method. In order to attach this aim, we resort to minimum spanning tree (MST). The main idea of MST-based approach is to generate a minimum spanning forest (MSF) from an undirected graph composed by all LiDAR raw points. The MSF is the union of MSTs for the connected components of the undirected graph under the threshold of slope. Terrain tends to be the one with the largest area.

2 Methodology

2.1 Undirected graph

A graph is a set of vertices and a collection of edges that each connects a pair of vertices. An undirected graph is a graph in which the vertices are connected by undirected arcs. Both ends of an undirected arc are equivalent—there is no head or tail, i.e. the cost from one vertex to another is same as that in reverse direction.

For LiDAR data-set, two points in neighbours are taken as a point pair and connected by an edge. The edge’s weight is approximated by the absolute slope between two points. An undirected graph of a LiDAR point data-set is then composed by regarding a point, an edge connecting a point pair and the related slope as a vertex, an edge and a weight, respectively. Therefore, the undirected graph of a LiDAR data-set containing $N$ points is $G=(V,E,W)$ with $V=\{V_i|i=1,2,...,N\}$, $E=\{e_{ij}|i\ne j\}$ and $W=\{w(e_{ij})\}$. $V$, $E$ and $W$ represent the set of vertices, edges and weights of the undirected graph $G$. $v_i$ is the $i-th$ point, $e_{ij}$ is the edge connecting $v_i$ and $v_j$, $w(e_{ij})$ is its weight.

To construct the undirected graph, two strategies are usually adopted. One is to connect the points in nearest neighbours and another is to construct TIN (Triangle Irregular Network) meshes. There are two ways to set the neighbours, one is k-nearest and another is $\mathit{\epsilon }$-ball neighbours. K-nearest strategy takes k nearest points as current point’s neighbours. For $\mathit{\epsilon }$-neighbours, distances from each point to its neighbours should be less than a threshold. According to their principles, the constructed undirected graphs are not necessarily connected especially when the LiDAR data are not uniformly distributed. For TIN-based strategy, all spatial points are firstly projected to 2D plane. The triangle meshes of those 2D points are then constructed on basis of Delaunay rule http://en.wikipedia.org/wiki/Delaunay_triangulation. The triangle meshes are formed by some edges connecting the neighbour points. After projecting these 2D meshes to the original points, the undirected graph of LiDAR points are actually obtained. In this situation, vertices in the graph are all connected if no any restriction is imposed. In the experiment of this article, TIN is adopted to construct undirected graph.

Another problem in constructing undirected graph is to determine the weight or cost of the edge. In this article, the weight of edge $e_{ij}$ is defined as the absolute slope in formula (1), where, $(x_i,y_i,z_i)$, $(x_j,y_j,z_j)$ are the coordinates of vertices $v_i$ and $v_j$.1 $$\begin{array}{c}\hfill w_{ij}=\frac{|z_i-z_j|}{\sqrt{{(x_i-x_j)}^2+{(y_i-y_j)}^2}}\end{array}$$1

2.2 Minimum spanning tree

In graph theory, a tree is a graph in which any two vertices are connected by exactly one path. Then, a single graph can have many different spanning trees. MST is the spanning tree with weight less than or equal to each of others. There are now two algorithms commonly used for MST finding, Prim and Kruskal algorithm, which are all based on greedy theory. By comparison, Kruskal algorithm is more flexible and has certain practical advantages. So, it is used in our algorithm, which consists of four steps as follows.

(1)	Sort $w$ of $G=(V,E,W)$ in numerical order.
(2)	Choose $e_{min}$ from $E$ with minimal weight $w_{min}$, remove the edge and its weight from $E$ and $W$, respectively. Add $e_{min}$ and $w_{min}$ to $E_{MST}$ and $W_{MST}$.
(3)	Remove $e_{min}$ and $w_{min}$ from $E_{MST}$ and $W_{MST}$ if a cycle in $G_{MST}=(V,E_{MST},W_{MST})$ is formed after adding $e_{min}$.
(4)	Repeat step (2) and (3) until $E$ and $W$ become empty. $G_{MST}$ is the MST of $G$. If there are edges with same weight, the MST may be not uniquely obtained. In this case, the edge firstly found is chosen as an edge to construct the MST.

Accordingly, minimum spanning forest(MSF) rather than MST is obtained when the undirected graph are not all connected. Each MST corresponds to a subgraph which is the connected part of undirected graph. MSF has the minimum cost when all MSTs are the lowest cost.

2.3 Minimum spanning tree-based filtering method

Suppose a LiDAR data-set is expressed as an undirected graph taking the absolute slope as edge’s weight and its MST is obtained through Kruskal algorithm computation. Then, object points and outliers including positive and negative outliers usually tend to be the end-vertices in the MST. However, the end-vertices may not be the outliers or objects. For instance, if there are no any man-made structures or outliers in the graph, the end-vertices in the MST still belong to terrain. To cover general situations, we should further judge the object probabilities of end-vertex by referencing the cost of the edge in which the end-vertex is located. The simple principle is illustrated as follows. Figure 1(a) is the undirected graph (triangle meshes) of some LiDAR points. Black line connects two neighbour points. Red number denotes the weight (slope) between two points. The MST of this undirected graph is then generated as the blue lines in Figure 1(b). If the threshold of absolute slope is 10, the edges whose weights are more than the threshold should connect the objects and terrains as the green lines in Figure 1(c). That is, the weights of edges connecting outliers, objects and terrains are relatively lager. This leads to the fact that outliers and objects tend to be the end-vertices in the MST. At the same time, the generated MST will be divided into three parts shown in Figure 1(d). The middle MST with largest area will be taken as terrain and others are considered as object.

Figure 1. Minimum spanning tree generation of an undirected connected graph. (a) shows the undirected graph for original points (circles) according to the delaunay rule. (b) is the generated MST which means each point pair is connected by only one path and the weight of this tree (sum of the edges’ weights) is minimum. In the MST, some edges’ weights are more than a threshold represented as dotted edges in (c). Those edges are removed from the original MST as a result of sub-MSTs (see (d)). Generally, the largest sub-MST will be taken as terrain and others are objects and outliers.

Figure 2. The Flowchart of MST-based DTM Detection.

It is obvious that the multi sub-MSTs can be obtained when some edges with larger weights are omitted in MST generation. In another word, those sub-MSTs can still be obtained by performing MST generation from the undirected graph in conditions of the weight of newly added edge is less than a threshold. Generally, the sub-MST with largest area is taken as the terrain. On basis of the object probabilities of end-vertices in MST and an estimated absolute slope threshold, our method is as follows. Firstly, we statistically compute a threshold which is the maximum absolute slope of terrains. Then, an undirected graph is constructed for LiDAR raw points on basis of 2D delaunay rule. Thirdly, Kruskal algorithm is used to generate MSF(MST). This process is terminated when the weight of edge to be added is more than the slope threshold. Finally, the sub-MST in the MSF with the largest area is taken as DTM. The flowchart of the proposed method is shown in Figure 2. Here, it should be mentioned that the original points can be divided into many parts when its number is huge. Combination of MSTs with largest area in different parts will be the terrain. This will accelerate the computation process.

2.4 Threshold

When TIN-based method is used to construct the undirected graph, there is only one parameter in our algorithm which is the absolute slope threshold. Threshold estimation is really a difficult and should try to fit to more different situations. To provide a unbiased estimation, we adopt the following strategy. It is computed by randomly picking several points from the original LiDAR points, computing the absolute slopes between each picked point and its neighbours, and adding the maximum absolute slopes of these samples. Then, the average of the maximum absolute slopes is taken as threshold.

It is obvious that the DTM classification is highly related to different samples and samples’ size especially when the sample’s size is small and the original data contains complex structure. Furthermore, the relationship between samples and the threshold can not be effectively estimated. To overcome this obstacle, we adopt relatively more random samples to compute the threshold. Another minor reason is that the computation of threshold is less sensitive to sample’s size. Then, we empirically choose 100–1000 points to compute threshold. The use of the maximum absolute slope makes it very likely that the threshold is more than the slope of terrain and less than the slope between terrain and object or outlier. The threshold is unbiased because it is the average of random selection.

2.5 Computation complexity

In our method, the main operations are to generate the undirected graph and construct its MST using Kruskal algorithm. When the TIN is taken to generate undirected graph and the method “divide and conquer” is used to produce TIN, the computation complexity should be $O(NlogN)$, $N$ is the number of original LiDAR points http://en.wikipedia.org/wiki/Delaunay_triangulation. Threshold estimation takes $O(N_s)$ operation, $N_s$ is the number of samples which should be much less than $N$. Generation of MST takes $O(ElogN)$ computation, $E$, $N$ are the edge and vertex number. $E$ should be several times of $N$. Totally, the computational complexity is about $O(NlogN)$.

3 Experiments and results

This part attempts to verify the effectiveness of MST-based algorithm by employing it in different scenarios. For comparison, Terrasolid and morphological approach in [5] are also applied here. In MST-based method, the undirected graph for each tested data-set is constructed according to 2D delaunay rule. The default parameters are used in Terrasolid. The parameters of morphological method are cell size ($c$), base of the exponential window ($b$), increment step for windows ($k$), slope ($s$), initial threshold ($d{h}_0$) and maximum threshold ($d{h}_{max}$). Table 1 lists the corresponding parameters for all experimental datas. The parameter $c$ is set to be half of average distance among points for each data-set. The Parameters $b$, $d{h}_0$ and $d{h}_{max}$ are set by referencing the values in [5]. $k=0,1,2,...,M$ and $2b^M+1$ must be greater than the maximum building size. Units of parameters $c$, $b$, $d{h}_0$ and $d{h}_{max}$ are all meter. Parameter $s$ is determined by referencing the absolute slope threshold calculated for MST-based method.

Table 1. Parameters of morphological method.

Parameter	Data-Set 1	Data-Set 2	Data-Set 3	Data-Set 4	Data-Set 5
c	1	1.5	0.3	1	1
b	2	2	2	2	2
k	0,1,2,...,4	0,1,2,...,5	0,1,2,...,6	0,1,2,...,4	0,1,2,...,5
s	0.4	0.35	0.2	0.3	0.35
$d{h}_0$	0.25	0.2	0.25	0.2	0.25
$d{h}_{max}$	2.5	80	2.5	12	2.5

3.1 Effectiveness

As mentioned in the introduction, the superior characteristic of our method to others is that it is robust to outliers including positive and negative outliers. To validate this characteristic, we use a synthetic data-set by randomly introducing some positive and negative outliers to a LiDAR raw data. Figure 3(a) and (b) show the synthetic data containing 35,043 points. Obviously, this data contains lots of outliers and complex structures. As we all known, morphological method and Terrasolid consider the lowest point within local region as terrain, and further discriminate other points based on which. However, these lowest points may also be negative outliers, rather than terrain points. In this situation, these two methods both get bad results. Figure 3(c) and (d) are the detected result of morphological method with 13,878 remaining points. For Terrasolid, the detected terrain points is 1,308 shown in Figure 3(e). Figure 3(f) shows our result containing 19,477 points. In Figure 3, blue points denote the terrain points and other points including green, yellow and red ones represent outliers or objects. According to Figure 3(c) and (d), we can see that morphological method fails to distinguish and remove negative outliers, and many terrain points are taken as object incorrectly. By comparison, MST-based method is rather robust to both positive and negative outliers and Terrasolid obtain worst result. The absolute slope threshold calculated is 0.395 for this data-set. In addition, the white areas correspond to non-terrain points being removed. The similar situation will be appeared in the following experiments. In practical applications, the DTMs of these white regions should be estimated. A simple method is to perform linear interpolations. How to do interpolation will not be discussed.

Figure 3. Data-set 1 of downtown area containing many positive and negative outliers. Blue points represent the terrains and others denote outliers or objects. Images a and b show the original raw points, c and d show the result of morphological method (a and c are the top views and b and d are the side view showing the outliers). Images e and f show the detected terrains using Terrasolid and our method, respectively.

3.2 Comparisons with other approaches

To further illustrate the validity of our method, a DTM with 46,151 points is used. Data-set 2 contains a number of small hills without any man-made structures. Terrasolid, morphology and MST are used to detect DTM from this data. Figure 4 shows the comparisons among them, in which, point’s color denotes its elevation. Yellow points denote the upland and purple points represent low-lying land. For morphology and Terrasolid, lots of terrain points are wrongly taken as object and the number of remaining points are 36,494 and 34,333, respectively. That is, they wrongly classify some terrain points as object. However, this can not be clearly found from the Figure 4 due to the points are very dense. In our result, there are 46,132 terrain points. Only 19 points are wrongly classified. The absolute slope threshold is 0.367 for data-set 2.

Figure 4. Data-set 2: DTM of mountains area. The colors of points denote theirs heights. Purple points are low-lying terrains and yellow points represent upland. (a) is original DTM, (b) is result of morphological method, (c) is the detected DTM of Terrasolid, (d) is our result.

Figure 5. Data-set 3: Raw data-set of downtown area. Purple points represent terrains and others denote objects. (a) original raw points, (b), (c), (d) detected terrains of morphological method, Terrasolid and MST-based method. Three results are close.

Data-set 3 has 34,533 points and it is the raw LiDAR data containing complex structures such as tall buildings, small buildings and few trees shown in Figure 5(a). Figure 5(b) shows the detected DTM using morphological method, Figure 5(c) shows the detected DTM using Terrasolid and (d) is result of MST-based method. Visually, three methods obtain comparable results for this data-set. In fact, Terrasolid and morphology wrongly take some terrains as object. The details can be found in next section. In Figure 5, the purple points represent terrain and other points denote the objects. The absolute slope threshold used in our method is 0.168.

Figure 6(a) is data-set 4 in flat regions containing a causeway and a small ridge. Although ridge and causeway are obviously higher than other terrain points, yet there are slow changes in the height. Thus, the causeway should belong to terrain. However, in another viewpoint, causeway belongs to man-made structure. As shown in 6(b) and 6(c), most of the causeway is taken as object by morphology and Terrasolid. MST-based method obtains result in 6(d), in which, the ridge and causeway are all retained due to they are all naturally connected with terrain points. In this sense, the ability of other two methods is relatively better than MST-based method. The absolute slope threshold calculated is 0.328 for this data-set.

Figure 6. Data-set 4: Raw data-set of rural region containing ridge and causeway. Purple points represent terrains and others denote objects. (a) is original raw points, (b) is result of morphological method, (c) is the detected DTM of Terrasolid, (d) is the of MST-based method.

Figure 7. Data-set 5: Raw data-set of downtown area containing bridge. Purple points shows terrains and others represent objects. (a) original raw points. (b), (c), (d) detected ground points of morphological method, Terresolid and MST-based method.

Table 2. Statistics of classification errors.

Filtered	Reference		Error
	Ter	Obj	Type I	Type II	Total
Terrain (Ter)	$a$	$b$	${\displaystyle \frac{c}{a+c}}$	${\displaystyle \frac{b}{b+d}}$	${\displaystyle \frac{b+c}{a+b+c+d}}$
Object (Obj)	$c$	$d$

Data-set 5 has 35,485 points and it is a downtown LiDAR raw data-set (see Figure 7(a)). The points on bridge are naturally connected with the terrain points and the MST-based method cannot discriminate them. Figure 7(b) shows the filtered result of morphology with 22,682 points taken as DTM. The points in the middle bottom are wrongly considered as outliers or object due to some small objects including trees and electric wires are existed in them. Then, we can know that morphology is sensitive to outliers and objects. Figure 7(c) shows the detected terrain points of Terrasolid with 20,260 points. This result still is not satisfying. Figure 7(d) is the result of MST-based method. For our method, most of points on objects including buildings, trees are detected, yet the points on bridge can not be discriminated. The absolute slope threshold of this data is 0.335.

Table 3. Parameters of morphological method.

Data-Set	Method	Filtered	Reference		Error			Time(ms)
			Ter	Obj	Type I	Type II	Total
1	Morphology	Ter	13,002	876	32.57%	5.58%	20.42%	4598
		Obj	6279	14,886
	Terrasolid	Ter	973	335	94.95%	2.13%	53.20%	5
		Obj	18,308	15,427
	MST	Ter	19,011	466	1.40%	2.96%	2.10%	657
		Obj	270	15,296
2	Morphology	Ter	36,494	0	20.92%	$\backslash$	20.92%	686
		Obj	9657	0
	Terrasolid	Ter	34,333	0	25.61%	$\backslash$	25.61%	30
		Obj	11,818	0
	MST	Ter	46,132	0	0.04%	$\backslash$	0.04%	296
		Obj	19	0
3	Morphology	Ter	12,022	7	15.61%	0.03%	6.46%	9395
		Obj	2224	20,280
	Terrasolid	Ter	11,619	5	18.44%	0.02%	7.62%	20
		Obj	2627	20,282
	MST	Ter	14,152	170	0.66%	0.84%	0.76%	168
		Obj	94	20,117
4	Morphology	Ter	21,960	74	24.07%	1.03%	19.47%	994
		Obj	6960	7125
	Terrasolid	Ter	22,427	5	22.45%	0.07%	17.99%	20
		Obj	6493	7194
	MST	Ter	28,698	1,697	0.77%	23.57%	5.31%	234
		Obj	222	5502
5	Morphology	Ter	22,653	29	19.65%	0.40%	15.69%	938
		Obj	5539	7264
	Terrasolid	Ter	20,250	10	28.17%	0.14%	22.41%	20
		Obj	7942	7283
	MST	Ter	28,153	1,692	0.14%	23.20%	4.88%	281
		Obj	39	5601

To further evaluate the ability of our method, we will perform quantitative evaluation for all tested data-set.

3.3 Quantitative evaluation

Our quantitative evaluation criteria are shown in Table 2 which references the work in [19, 20].

In Table 2, $a$ is the correctly detected number of terrain points, $b$ is the number of object points that have been incorrectly identified as terrain (contribute to Type II error), $c$ is the number of terrain points that have been incorrectly identified as object (contribute to Type I error), $d$ is the number of correctly identified object.

The statistics of errors of above-mentioned data-sets are as Table 3. From the statistics, we can see that the type I errors of morphological method and Terrasolid are relatively large, which indicates that these two methods tend to consider terrains as objects. For data-set 4 and 5, the type II errors of our method is very big. This is because bridges are naturally connected with the terrain points. At the same time, we can see that the computation of our method is less than that of morphology and is more than that of Terrasolid. Especially for data-set 3, the cost of morphology method is rather higher. It is due to more iterations are needed to do classification for the urban point data-set. By comparison, the the Terrasolid runs rather fast (the times are approximately estimated since we cannot exactly capture the beginning and ending of the programme). Then, MST-based method should be accelerated.

Figure 8. Data-set 5: Post processing. Purple points shows terrains and others represent objects. (a) is the detected terrain points of MST-based method, (b) is the post processing result of (a) using morphological opening operation. (c) is the corresponding result of TIN-based method.

According to the figures and tables, our method can obtain comparable results in general situations and perform better when data-set contains lots of outliers. The limitation lies in it cannot discriminate the object points when they are naturally connected with terrain points. Thereby, post processing should be done.

3.4 Post processing

To further detect man-made object from the results of MST-based method, morphological operation or TIN-based approach is used to perform further processing in this article. Figure 8(a) is the initial detected terrain of data-set 5, 8(b) is the post processing result using morphological opening operation and 8(c) is the corresponding result of TIN-based method. Obviously, the bridge in 8(a) is removed. For data in Figure 8(b), the error of Type I and Type II are 1.62% and 0.02%. The corresponding values are 0.71% and 0.01%. That is, only a fewer object points are wrongly classified into terrain. This character makes our method be a effectively preprocess way in filtering.

4 Conlusion

In this article, we propose a new DTM detection method from LiDAR raw points based on minimum spanning tree. The distinctive characteristic of the proposed method is that it is robust in the scenarios when raw data-set contains lots of negative outliers. Other merits are that this method is simple and only one parameter is to be set. Through experiments, its effectiveness is demonstrated. However, this method is a bit limited due to it cannot correctly discriminate object points when they are naturally connected with terrains. Then, post processing may be needed. Another deficiency is that the MST with the largest area is taken as the terrain points in this article, which is true in most situations. Yet, the terrain points may be other MST or combination of some MSTs in a few scenarios. Our future work is to develop the corresponding techniques to strengthen this algorithm.

Acknowledgments

This study was supported by the Natural Science Foundation Fund of China under contract 41371426 and the Natural Science Foundation Fund of Hubei Province under contract 2012FFB04202.