J. Cent. South Univ. (2017) 24: 1826-1836

DOI: https://doi.org/10.1007/s11771-017-3591-9

A novel SMC-PHD filter based on particle compensation

XU Cong-an(徐从安)¹, HE You(何友)¹, YANG Fu-cheng(杨富程)^{1, 2}, JIAN Tao(简涛)¹,

WANG Hai-peng(王海鹏)¹, LI Tian-mei(李天梅)³

1. Research Institute of Information Fusion, Naval Aeronautical and Astronautical University, Yantai 264001, China;

2. School of Electronics and Computer Science University of Southampton, SO17 1BJ, United Kingdom;

3. 302 Unit, Second Artillery Engineering University, Xi’an 710025, China

Central South University Press and Springer-Verlag GmbH Germany 2017

Abstract: As a typical implementation of the probability hypothesis density (PHD) filter, sequential Monte Carlo PHD (SMC-PHD) is widely employed in highly nonlinear systems. However, the particle impoverishment problem introduced by the resampling step, together with the high computational burden problem, may lead to performance degradation and restrain the use of SMC-PHD filter in practical applications. In this work, a novel SMC-PHD filter based on particle compensation is proposed to solve above problems. Firstly, according to a comprehensive analysis on the particle impoverishment problem, a new particle generating mechanism is developed to compensate the particles. Then, all the particles are integrated into the SMC-PHD filter framework. Simulation results demonstrate that, in comparison with the SMC-PHD filter, proposed PC-SMC-PHD filter is capable of overcoming the particle impoverishment problem, as well as improving the processing rate for a certain tracking accuracy in different scenarios.

Key words: random finite set (RFS); probability hypothesis density (PHD); particle filter (PF); particle impoverishment; particle compensation; multi-target tracking (MTT)

1 Introduction

With the increasing complexity of the modern war environments, multi-target tracking (MTT) is facing various technical challenges during the past decades. Generally, MTT can be viewed as a process which jointly determines the number and states of targets from a series of discrete measurements in conditions of uncertain association, uncertain detection, and clutters [1]. Many classical solutions have been studied in literatures, such as the joint probabilistic data association (JPDA) filter [2] and the multiple hypothesis tracking (MHT) filter [3], which both include a matching strategy that aims at optimizing the association between targets and measurements. In recent years, a range of alternative solutions based on random finite set (RFS), which avoid immediate data associations required in classical solutions, has been proven to be an effective way for the MTT problem [4-6]. However, in most practical applications, the multi-target Bayes filter is difficult to implement due to the multiple high-dimensional space integrals. To alleviate this intractability, the probability hypothesis density (PHD) filter [5] has been proposed and generated substantial interest.

The implementations of the PHD filter mainly include Gaussian mixture PHD (GM-PHD) filter [7], sequential Monte Carlo PHD (SMC-PHD) filter [8] and their modified versions [9-13]. The GM-PHD is constrained due to the linear Gaussian assumption. Compared to the GM implementation, the SMC implementation, which approximates the PHD by a set of random weighted particles, does not require any assumptions, and is applicable for most general situations. However, the SMC-PHD filter also suffers from some problems in practical applications.

The first one is the high computational burden problem. The computational burden of SMC-PHD filter highly grows with the number of targets and measurements, which limits the use for real-time applications and has naturally become one of the biggest shortages for the SMC-PHD filter [14]. In order to solve the problem, gating technique is widely used to increase the filtering speed [14-17]. The principle of gating technique in existing works is similar, i.e., classify either the measurements [14-16] or the particles [17] for different types of targets based on the differences between state estimations and measurements. The second one is the particle impoverishment problem. As a typical particle implementation of the PHD, like particle filter (PF) [18], SMC-PHD filter also faces the particle impoverishment problem. To the best of the authors’ knowledge, aiming at overcoming this drawback and improving the estimation accuracy, some approaches have been put forward, such as regularized particle filter (RPF) [19], Markov chain Monte Carlo (MCMC) algorithm [20], genetic algorithm (GA) [21] and distributed resampling algorithm [22]. However, these algorithms may either improve the computation complexity, or are mainly employed for some special applications. And they cannot be used in PHD filter, in other words, there are not specific algorithms for solving the particle impoverishment problem of SMC-PHD filter.

In this work, a novel SMC-PHD filter based on particle compensation is proposed, where the novelty lies in the way of developing a new particle generating mechanism to compensate particles and maintain the diversity of particles. The main contributions of this paper are threefold. Firstly, a comprehensive analysis on the particle impoverishment problem of SMC-PHD filter is present. Secondly, a new particle generating mechanism is developed, where the particles are not only drawn from the proposal density but also from the compensatory probability density. Thirdly, based on the particle generating mechanism, a particle compensation SMC-PHD (PC-SMC-PHD) filter is proposed.

2 Background and problem formulation

2.1 PHD filter

In the MTT scenario, due to time-varying numbers of targets and the inherent defects existing in the detector, the states of the targets and the sources of measurements are unknown. Naturally, the uncertainty in both multi-target measurements and multi-target states can be modeled by RFS. Suppose at time k there are n_k targets with states x_k,1, …, x_k,nk in the state space X. Similarly, the measurements modeled by z_k,1, …, z_k,mk take values in the measurement space Z. Then the multi-target states and measurements set are represented by

                  (1)

                   (2)

where F(X) and F(Z) are the finite subsets of X and Z, respectively.

According to Bayesian estimation theory, the multi-target posterior probability density can be given by the optimal Bayes recursion as follows [4]:

      (3)

 (4)

where is the multi-target probability density, is the state transition density and is the multi-target likelihood.

The multi-target Bayes recursion involves multiple integrals and the computational complexity grows exponentially with the number of targets. As a practical suboptimal alternative, the PHD filter was proposed to alleviate the computational burden under FISST framework. The PHD is a multi-peak function in the state space; the peak number and the peak position are corresponding to the number of the targets and the expected value of the target states, respectively. Assume that denote the PHD at time k-1, then the prediction and update equations are given by [5]

       (5)

   (6)

whereandare the target survival probability and target detection probability, respectively; denotes the intensity of the spontaneously birth target RFS and denotes the intensity of the spawned target RFS; is the intensity of the clutter RFS, λ_k is the number of clutters which obeys the Poisson distribution, and c_k(·) is the probability density of the clutters.

2.2 SMC-PHD filter

As an approximate implementation of PHD filter, SMC-PHD filter ⁸ is considered as a promising filter for MTT problems. Given the setof particles and their associated weights approximating the PHD at time k-1:

                  (7)

where δ(·) is the Dirac delta function. Then, the SMC-PHD filter procedure can be briefly described as follows.

1) Prediction step

Considering the survival and newborn targets, the predicted PHD is given by

             (8)

where J_k denotes the number of new particles.

is drawn according to the proposal densities

      (9)

where q_k and p_k are the proposal distribution and the probability density of newborn target. And the predicted weights are defined as follows:

 (10)

2) Update step

The updated PHD can be written as

                (11)

where the weights are updated according to

       (12)

and

     (13)

3) Resampling step

Compute the total particle mass:

                (14)

To avoid the particle degeneracy, we resample from the set to obtain and rescale the weights by to obtain a new one

4) State estimation step

The expected target number N_k can be computed by

                             (15)

whereis the integer nearest to In this work, K-means clustering algorithm is used to determine the centers of clusters, which forms the particle approximation of the PHD function. These centers represent the peaks of the PHD function, which can be used as target state estimates.

2.3 Analysis on particle impoverishment problem of SMC-PHD filter

As previously mentioned, to avoid particle degeneracy and prevent exponential growth of the particle size, the resampling step is recommended in standard SMC-PHD filter, which obtains a new particle set by duplicating particles with larger weights and deleting particles with smaller weights. However, through a comprehensive analysis on particle evolution after resampling, some defects can be found. To show the particle evolution process, a vivid illustration in Fig. 1 is given.

Fig. 1 A vivid illustration of particle evolution process for SMC-PHD filter

In Fig. 1, true target positions and corresponding particles after resampling are present at five different times. The black circles denote true target positions; the gray circles denote the particles and the bigger circles mean bigger weights. As seen from Fig. 1, there are a mass of different particles around the true target at time k-2, while at time k+1, the remaining particles are only two types because of the so called “particle collapse”. That is to say, the diversity of particles gets worse as time increases. Especially, at time k-2, no particle surrounds the true target any more, which implies that the SMC-PHD filter loses the target.

From the particle evolution process, it is clear that resampling introduces the problem of loss of diversity among the particles, referred as particle impoverishment problem, which may lead to poor tracking performance. As a consequence, some improvement should be done to overcome the revealed defect and enhance the performance of SMC-PHD filter.

3 Particle compensation based SMC-PHD filter

To solve the particle impoverishment problem, a novel SMC-PHD filter based on particle compensation (called PC-SMC-PHD) is proposed in this work. Based on the new particle generating mechanism, the PC-SMC-PHD maintains the diversity as while as keeps the effectiveness of particles.

3.1 Prediction step (new particle generating mechanism)

The excessive copying of the particles with high weight is the main reason for impoverishment phenomenon. For the purpose of improving the diversity of the particles, new particle generating mechanism is used to compensate the particles. Note that unlike the standard SMC-PHD filter, here particles are not only drawn from the proposal density but also from the compensatory probability density based on the one step prediction of the estimations.

At time k-1, suppose that the union of the two particle sets

 (16)

is the random sample representation of , i.e.,

 (17)

where and respectively are the sample particle sets based on the standard method and the compensatory method; and respectively are the standard PHD and the compensatory PHD, which correspond to the standard case and compensatory case.

In the standard case, sample the predicted particles and compute their weights to approximate For

                        (18)

  (19)

where we adoptfor simplicity. For

                            (20)

                   (21)

Then,is approximated by the particle set i.e.,

           (22)

While in the compensatory case, the idea is to use the multi-target estimation setto place the compensatory particles and obtain . The set can be given by

                            (23)

For each , a set of n_b new particles are generated by the compensatory probability density , which can be given by

            (24)

where

                  (25)

The total number of new particles generated in this way is

                                (26)

The weights of the particles are uniform, i.e.,

                            (27)

Then, is approximated by the particle set , i.e.,

           (28)

It should be noted thatincludes accurate velocity estimation, so the compensatory method by generating new particles can track the targets with much less particles than the standard method, which needs prior velocity information (inaccurate normally).

3.2 Update step

After generating the new particles, consider the update step of the PHD filter. In the standard case, can be approximated by the weighted sum of particles:

             (29)

where

    (30)

and

                  (31)

Similarly, is approximated by the weighted particle set, i.e.,

                (32)

where

                           (33)

3.3 Resampling step

In order to eliminate the samples with small weights and multiple the samples with large weights, the particle setapproximatingis next resampled times. Compute the total particle mass:

                           (34)

where denotes the approximation of the expected target number. After resampling, is approximated by

                   (35)

where .

The particle setapproximating can be also resampled (actually rescale particle weights). Compute the total particle mass:

                              (36)

wherecorresponds to the expected estimation number. After resampling, is approximated by

                  (37)

where

It is worth noting that the expected target number output by PC-SMC-PHD filter is , rather than the sum of and . State estimation is performed efficiently and accurately using the same procedure as standard SMC-PHD filter. For clarity, pseudo-code for the PC-SMC-PHD filter is given below.

Algorithm 1: Pseudo-code for PC-SMC-PHD filter

1 Input:

21) Particle set approximating

3 2) Particle set approximating

4 3) Multi-target estimation set at time k-1:

5 Union of input particle sets

6 for do

7 Draw

8 Compute weight

9 end for

10 for do

11 Draw

12 Compute predicted weight

13 end for

14

15for do

16 for do

17

18 Draw

19 Compute predicted weight

20 end for

21 end for

22 Compute updated weight using Eq. (30)

23 Compute updated weight using Eq. (33)

24 Compute the expected target numberusing Eq. (34)

25 Compute the expected estimation number using Eq. (36)

26

27 Resample times fromto obtain , with

28 Resampletimes fromto obtain , with

29 Output:

30 1) Expected target number

31 2) Particle setsand

4 Simulation results

In this section, we present simulation results to demonstrate the performance of the PC-SMC-PHD filter using a linear motion model and a nonlinear turn model, respectively. Scenario I considers targets with the linear model and Scenario II considers targets with the nonlinear turn model. All simulations are given in the same experimental environment (MATLAB 2009a, Pentium E6700 3.2 GHz processor with 2 GB of RAM).

4.1 Scenario I: Linear model

Assume that there is an unknown and time-varying number of targets moving in the two-dimensional scenariofor a period of 100 time steps. The sampling interval of the sensor is T=1 s. The true target trajectories are shown in Fig. 2.

At time k, let denotes the state vector, where x_k and y_k are the states of target position, and and are the states of target velocity. The single-target transition model is given by

                           (38)

where the target transition matrix F_k is given by

                         (39)

Fig. 2 True target trajectories

and the process noise v_k-1 follows a Gaussian distribution with zero-mean and covariance Q, which is given by

            (40)

where σ_v=0.1 denotes standard deviation of the process noise. The probability of survival is set to P_s=0.98. The bearing θ_k and range r_k measurement of a target are generated by a radar sensor located at [0, 0]^T and are given by

                    (41)

where w_1,k and w_2,k are zero-mean Gaussian white noise with respective standard deviations σ_θ=0.01 and σ_r=5. The probability of detection is P_D=0.95. The birth process follows a Poisson RFS with intensity:

                   (42)

where Clutter is uniformly distributed in surveillance region; the Poisson rate for the number of clutter measurements per time step is set to be λ=5. In Scenario I, 500 particles are used to represent one target, J_k=1000 particles are used to represent the newborn targets and n_b=200 particles are compensated for every estimation.

Figure 3 shows estimates of target positions in a trial by standard SMC-PHD filter and PC-SMC-PHD filter. As seen from Fig. 3(a), particle impoverishment problem of SMC-PHD filter discussed before leads to target loss. While as seen from Fig. 3(b), the proposed PC-SMC-PHD filter avoids the problem and can track all the targets effectively.

Fig. 3 Estimates of target positions in a trial:

Figure 4 shows estimates of target number averaged over 50 Monte Carlo trials. It can be seen that the PC-SMC-PHD filter has higher accuracy than the SMC-PHD filter, which underestimates the target number in most time. In other words, by generating new particles, the proposed filter increases the particle diversity and significantly improves the performance of number estimation.

To evaluate the average performance of the two filters, 50 Monte Carlo trials are performed and the OSPA distance [23] is introduced, which jointly captures differences in cardinality and individual elements between two finite sets. Let X and Y denote the sets of state vectors, then the OSPA is defined by

                (43)

where the parameters are set to be p=2 and c=30. If m>n,

Fig. 4 Estimates of target number averaged over 50 Monte Carlo trials

Figure 5 shows the average OSPA distance estimated by the two filters. It can be seen that the OSPA distance of the PC-SMC-PHD is much smaller than that of SMC-PHD, and the SMC-PHD suffers the particle impoverishment problem as time increases.

Fig. 5 Average OSPA distance estimated by SMC-PHD and PC-SMC-PHD

For further evaluation of the two filters, Fig. 6 shows the average OSPA distance versus the number of particles (λ=5). As seen, the average OSPA distance declines with the increase of the number of particles for the both filters, and the average OSPA distance of PC-SMC-PHD is much smaller with the same number of particles. That is because the proposed filter avoids the particle impoverishment problem by generating new particles. Figure 7 shows the average OSPA distance versus the clutter density. As seen, there is an overall increase of average OSPA distance with the increment of the clutter density. And it appears that the PC-SMC-PHD filter has much better performance than that of the standard SMC-PHD filter.

The average simulation time and average OSPA distance are listed in Table 1. As seen, the computational loads of two filters grow with the increase of the number of particles. From Table 1, when J_k=5000, average simulation time and average OSPA distance of SMC-PHD are respectively 278.20 s and 16.32, and when J_k=500, those of PC-SMC-PHD are respectively 136.86 s and 13.05. That is to say, compared with SMC-PHD filter, even when the average OSPA distance is smaller, the simulation time of the PC-SMC-PHD filter is highly reduced by (278.20-138.86)/278.20=50.81%. This indicates that the PC-SMC-PHD filter has faster processing rate for a certain accuracy of average OSPA distance.

Fig. 6 Average OSPA distance versus number of particles (λ=5)

Fig. 7 Average OSPA distance versus clutter density (J_k=1000)

Table 1 Average simulation time and average OSPA distance (λ=5)

4.2 Scenario II: Nonlinear model

In Scenario II, each target moves according to the nonlinear turn model. At time k, denotes the state vector, where Ω_k is the turn rate. Then, the target transition matrix with unknown Ω_k is given by

 (44)

And covariance Q is given by

 (45)

where The birth process can be given by Eq. (42), where , and

In Scenario II, the state vector has higher dimensions, so we adopt more particles to track the targets. Here, 500 particles are used to represent one target, J_k=2000 particles are used to represent the newborn targets and n_b=500 particles are compensated for every estimation. Other parameters are the same as in Scenario I. The true target trajectories are shown in Fig. 8.

Figure 9 shows estimates of target positions in a trial by standard SMC-PHD filter and PC-SMC-PHD filter. As seen from Fig. 9, SMC-PHD filter suffers from the particle impoverishment problem discussed before and loses one target, while the proposed PC-SMC-PHD filter avoids the problem and can track all the targets effectively.

Fig. 8 True target trajectories

Fig. 9 Estimates of target positions in a trial:

Figure 10 shows estimates of target number averaged over 50 Monte Carlo trials. As seen, the PC-SMC-PHD filter has higher accuracy of target number estimation than the SMC-PHD filter. That is because the proposed filter increases the particle diversity by generating new particles, which can ensure the tracking continuity and improve the performance of number estimation.

Fig. 10 Estimates of target number averaged over 50 Monte Carlo trials

Figure 11 shows the average OSPA distance estimated by the two filters. It can be seen that the OSPA distance of the PC-SMC-PHD is much smaller than that of SMC-PHD, and the gap is widening as time increases.

Fig. 11 Average OSPA distance estimated by SMC-PHD and PC-SMC-PHD

Figure 12 shows the average OSPA distance versus the number of particles (λ=5). As seen, similar to Scenario I, the average OSPA distance declines with the increase of the number of particles for the both filters, and owing to employing the new particle generating mechanism, the PC-SMC-PHD has much smaller average OSPA distance with the same number of particles.

Fig. 12 Average OSPA distance versus number of particles (λ=5)

Figure 13 shows the average OSPA distance versus the clutter density. As expected, there is an overall increase of average OSPA distance with the increment of the clutter density. And the PC-SMC-PHD filter has much better performance than that of the standard SMC-PHD filter.

Fig. 13 Average OSPA distance versus clutter density (J_k=2000)

The average simulation time and average OSPA distance are listed in Table 2. As seen, the computational loads of two filters grow with the increase of the number of particles. Like in Scenario I, we observe these two statistics, i.e., average simulation time and average OSPA distance of SMC-PHD when J_k=5000 (453.22 s and 20.47, respectively), and those of PC-SMC-PHD when J_k=500 (325.02 s and 14.34, respectively). That is to say, compared with SMC-PHD filter, even when the average OSPA distance is smaller, the simulation time of the PC-SMC-PHD filter is highly reduced by (453.22- 315.02)/453.22=30.49%. This indicates that the PC- SMC-PHD filter has faster processing rate for a certain accuracy of average OSPA distance.

Table 2 Average simulation time and average OSPA distance

5 Conclusions

The particle impoverishment problem of SMC-PHD filter is studied and a novel SMC-PHD filter based on particle compensation is proposed. Aiming at maintaining the diversity and keeping the effectiveness of particles, some compensatory particles are added into the particle system by developing a new particle generating mechanism, and then all the particles are naturally integrated into the SMC-PHD filter framework (Prediction step, Update step and Resampling step). Simulation results demonstrate that the PC-SMC-PHD filter can overcome the particle impoverishment problem well in some different scenarios. Meanwhile, compared with the SMC-PHD filter, the PC-SMC-PHD filter has faster processing rate for a certain tracking accuracy.

References

[1] BAR-SHALOM Y, FORTMANN T E. Tracking and data Association [M]. San Diego, CA, USA: Academic, 1988: 15-18.

[2] FORTMANN T E, BAR-SHALOM Y, SCHEFFE, M. Sonar tracking of multiple targets using joint probabilistic data association [J]. IEEE Journal of Oceanic Engineering, 1983, 8(3): 173-184.

[3] BLACKMAN S. Multiple hypothesis tracking for multiple target tracking [J]. IEEE Transactions on Aerospace and Electronic Systems, 2004, 19(1): 5-18.

[4] MAHLER R. Statistical multisource-multitarget information fusion [M]. Norwood: Artech House, 2007: 20-25.

[5] MAHLER R. Multi-target Bayes filtering via first-order multi-target moments [J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(4): 1152-1178.

[6] MAHLER R. PHD filters of higher order in target number [J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(4): 1523-1543.

[7] VO B N,MA W K.The Gaussian mixture probability hypothesis density filter [J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4091-4104.

[8] VO B N, SINGH S, DOUCET A. Sequential Monte Carlo methods for Bayesian multi-target filtering with random finite sets [J]. IEEE Transactions on Aerospace and Electronic Systems, 2005, 41(4): 1224-1245.

[9] SHI Zhi-guo, LIU Yong-kang, HONG Shao-hua, CHEN Ji-ming, SHEN Xue-min. POSE: Design of hardware-friendly particle-based observation selection PHD filter [J]. IEEE Transactions on Industrial Electronics, 2014, 61(4): 1944-1956.

[10] LIAN Feng, HAN Chong-zhao, LIU Wei-feng. Novel extension of the probability hypothesis density filter for tracking multi-target in clutter with an unknown intensity [J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46: 2066-2078.

[11] WU Xin-hui, HUANG Gao-ming, GAO Jun. Particle filters for probability hypothesis density filter with the presence of unknown measurement noise covariance [J]. Chinese Journal of Aeronautics, 2013, 26(6): 1517-1523.

[12] LI Tian, CORCHADO J M, SIYAU M F. A particle dyeing approach for track continuity for the SMC-PHD filter [C]// 17th International Conference on Information Fusion. Salamanca, Spain, 2014: 1-8.

[13] LI Wen-ling, JIA Ying-min. Gaussian mixture PHD filter for jump Markov models based on best-fitting Gaussian approximation [J]. Signal Processing, 2011, 91: 1036-1042.

[14] ZHENG Yun-mei, SHI Zhi-guo, LU Rong-xing. An efficient data-driven particle PHD filter for multi-target tracking [J]. IEEE Transactions on Industrial Informatics, 2012, 9(4): 2318-2326.

[15] SWAIN A, CLARK D. The single-group PHD filter: An analytic solution [C]// 17th International Conference on Information Fusion. Salamanca, Spain, 2014, 41-48.

[16] LI Tian-cheng, SUN Shu-dong, Sattar T P. High-speed Sigma-gating SMC-PHD filter [J]. Signal Processing, 2013, 93(3): 2586-2593.

[17] WANG Ya-dong, WU Jian-kang, KASSIM A. Data-driven probability hypothesis density filter for visual tracking [J]. IEEE Transactions on Circuits and Systems for Video Technology, 2008, 18(8): 1085-1095.

[18] CARVA H, MORAL P D. Optimal Nonlinear Filtering in GPS/INS Integration [J]. IEEE Transactions on Aerospace and Electronic Systems, 1997, 33(2): 835-849.

[19] ARULA M S, MASKELL S, GORDON N, CLAPP T. A tutorial on particle filters for online nonlinear/ non-Gaussian Bayesian tracking [J]. IEEE Transactions on Signal Processing, 2002, 50(2): 174-188.

[20] FEARNH P. Markov chain Monte Carlo, sufficient statistics and particle filters [J]. Journal of Computional and Graphical Statistics, 2002, 11(4): 848-862.

[21] YE L, WANG J L, ZHANG Q. Genetic resampling particle Filter [J]. Acta Automatica Sinica, 2007, 33(8): 885-897.

[22] BOLIC M, DJURIC P M, HONG S J. Resampling algorithms and architectures for distributed particle filters [J]. IEEE Transactions on Signal Processing, 2005, 53(7): 2442-2450.

[23] SCHU D, VO B T, VO B N. A consistent metric for performance evaluation of multi-object filters [J]. IEEE Transactions on Signal Processing, 2008 56(8): 3447-3457.

(Edited by YANG Bing)

Cite this article as: XU Cong-an, HE You, YANG Fu-cheng, JIAN Tao, WANG Hai-peng, LI Tian-mei. A novel SMC-PHD filter based on particle compensation [J]. Journal of Central South University, 2017, 24(8): 1826-1836. DOI: https://doi.org/10.1007/s11771-017-3591-9.

Foundation item: Projects(61671462, 61471383, 61671463, 61304103) supported by the National Natural Science Foundation of China; Project(ZR2012FQ004) supported by the Natural Science Foundation of Shandong Province, China

Received date: 2015-12-14; Accepted date: 2016-12-26

Corresponding author: HE You, Professor, PhD; Tel: +86-535-6635877; E-mail: xcatougao@163.com