Sequential Bearings-Only-Tracking Initiation with Particle Filtering Method

Abstract

The tracking initiation problem is examined in the context of autonomous bearings-only-tracking (BOT) of a single appearing/disappearing target in the presence of clutter measurements. In general, this problem suffers from a combinatorial explosion in the number of potential tracks resulted from the uncertainty in the linkage between the target and the measurement (a.k.a the data association problem). In addition, the nonlinear measurements lead to a non-Gaussian posterior probability density function (pdf) in the optimal Bayesian sequential estimation framework. The consequence of this nonlinear/non-Gaussian context is the absence of a closed-form solution. This paper models the linkage uncertainty and the nonlinear/non-Gaussian estimation problem jointly with solid Bayesian formalism. A particle filtering (PF) algorithm is derived for estimating the model’s parameters in a sequential manner. Numerical results show that the proposed solution provides a significant benefit over the most commonly used methods, IPDA and IMMPDA. The posterior Cramér-Rao bounds are also involved for performance evaluation.

1. Introduction

The problem of angle/bearings-only tracking (BOT) has been studied over many years due to its tremendous importance in a variety of practical communication and signal processing applications, such as localization in wireless sensor networks [1], submarine tracking using passive sonar [2], aircraft surveillance using radars [3], and sonar-based robotic navigation [4, 5]. This paper is concerned with radar and sonar tracking. Specifically, the focus of this paper lies in a fundamental problem in many sonar and radar tracking tasks, called track initiation. This problem consists in linking sets of point observations from different time steps to fit a desired model without any a priori track estimates. In general, track initiation approaches can be categorized into two schemes: sequential and batch schemes [6–8]. The sequential scheme processes a sequence of measurements received during consecutive radar/sonar scans sequentially one at a time. The sequential methods are widely used in radar and sonar tracking problems. For the batch technique approach, 𝑁 past scans of measurements are treated simultaneously to determine feasible target trajectories (see Figure 2). The batch techniques suffer from a heavy computational load and a slow process, while they find applications in image processing and tracking in heavy clutter background [6]. This paper only concerns sequential techniques.

In the context of radar and sonar applications, the target signal, if present, is often corrupted by additive noise and comes with other measurements unrelated to the target, which originate from thermal noise and various forms of clutter such as terrain and clouds [9]. So the issue of measurement origin uncertainty has to be addressed before tracking. This is called the data association problem and has been addressed in the context of multitarget tracking (MTT) in clutter [10, 11].

The tracking initiation problem considered here has its own challenges not only due to the measurement origin uncertainty, but also the uncertainty in the target presence. In fact, for tracking initiation, it is prerequisite to determine whether the target is present before dealing with data association and state filtering. Moreover, in BOT the measurement function is highly nonlinear, which results in a non-Gaussian posterior probability density function (pdf) in the optimal estimate of the target. The consequence of this nonlinear/non-Gaussian context is the absence of a closed-form solution [12–14]. As BOT initiation involves such interactive problems, that is, the measurement origin uncertainty, the target presence uncertainty, and the nonlinear/non-Gaussian posteriors in optimal estimation, it is therefore much more difficult than any individual-related tasks such as state filtering, target detection, and data association. It is noted that, besides the simulated annealing [15, 16] and the Ants algorithm [17], there are few reports on such BOT initiation problem, not to mention sequential BOT initiation.

This paper proposes an algorithm for sequential BOT initiation. The measurement origin uncertainty, target presence uncertainty, and the nonlinear non-Gaussian factors are handled jointly within a Bayesian sequential estimation framework. Based on such Bayesian formalism, a sequential Monte Carlo (SMC), a.k.a particle filtering (PF), algorithm is derived. Performance of the proposed approach is evaluated via numerical simulations and related methods are involved for performance comparison.

The remainder of this paper is organized as follows. Section 2 formulates the sequential BOT initiation problem. Section 3 presents the Bayesian sequential estimation framework for the problem under consideration. Section 4 derives the PF algorithm under the framework developed above. In Section 5, simulations are performed to evaluate the performance of the proposed algorithm. Finally, Section 6 concludes the paper.

2. Problem Formulation

2.1. Dynamic Model of the Target

Consider a discrete-time dynamic model 𝐱_𝑘+1=𝐟_𝑘𝐱_𝑘,𝐯_𝑘, (1) which is used to describe the target movement in a two-dimensional (2D) 𝑥-𝑦 plane. 𝑘 denotes the discrete time index, 𝐯_𝑘 the process noise, and 𝐱_𝑘 the target state vector. The target state consists of 2D position and velocity elements, which is defined as below: 𝐱_𝑘=𝑥_𝑘𝑥̇_𝑘𝑦_𝑘𝑦̇_𝑘^𝑇, (2) where (𝑥,𝑦) and (𝑥̇,𝑦)̇ denote the 2D position and velocity, respectively. The superscript 𝑇 denotes matrix transposition. If the movement pattern of the target is a priori known, then 𝐟 will have a specific form. For example, if the target moves with a near constant velocity, then 𝐟 can be modeled to be 𝐟_𝑘𝐱_𝑘,𝐯_𝑘=𝐅𝐱_𝑘−1+𝐯_𝑘, (3) where 𝐅=^𝐅𝐬0_0𝐅𝐬,  𝐅_𝐬=^1T₀₁, and T denotes the sampling period of the measurements. For more details on dynamic models for tracking problems, see [18].

2.2. Markov Chain Model for Uncertainty of Target Existence

In tracking initiation problem, the target presence is not a definitive issue either before the track under consideration is confirmed or after it is terminated. To take into account the uncertainty of the target presence, we introduce binary valued variable 𝜆_𝑘∈{0,1} and let 𝜆_𝑘=1 represent the target being present in the surveillance region at time 𝑘, and 𝜆_𝑘=0 denotes the opposite. Evolution of 𝜆 during the surveillance time is modeled by a Markov chain, which is characterized by the following transition probabilities: 𝑃_𝑏≜𝑃𝜆_𝑘=1∣𝜆_𝑘−1,𝑃=0_𝑑≜𝑃𝜆_𝑘=0∣𝜆_𝑘−1=1, (4)(5) and the initial probability of target existence, that is, 𝑝(𝜆₁). 𝑃_𝑏 and 𝑃_𝑑, have specific physical meanings. 𝑃_𝑏 is actually the probability of a new target’s birth or appearing at a time step and 𝑃_𝑑 is the probability of the target’s death or disappeance. The transitional probability matrix is given by Π=1−𝑃_𝑏𝑃_𝑏𝑃_𝑑1−𝑃_𝑑, (6) where 1−𝑃_𝑑 and 1−𝑃_𝑏 are equivalent to probabilities of target staying alive and remaining absent, respectively. 𝑃_𝑏, 𝑃_𝑑, and 𝑝(𝜆₁) are assumed a priori known.

2.3. Measurement Model

Let 𝐙_𝑘≜(𝐳₁,𝐳₂,…,𝐳_𝑘) denote the measurements received up to and including the time step 𝑘. Here, 𝐳_𝑘 represents the measurement received at time 𝑘. It is composed of only clutters (if the target is absent currently) or clutters along with the target originated measurement element (if the target is present at this time step). Specifically, 𝐳_𝑘≜𝛽_𝑘,𝑧_1,𝑘,𝑧_2,𝑘,…,𝑧_{𝑚𝑘,𝑘},if𝜆_𝑘𝑧=1_1,𝑘,𝑧_2,𝑘,…,𝑧_{𝑚𝑘,𝑘},if𝜆_𝑘=0, (7) where 𝛽_𝑘 denotes the target originated measurement, 𝑧_1,𝑘,𝑧_2,𝑘,…,𝑧_{𝑚𝑘,𝑘} measurement elements due to clutters, and 𝑚_𝑘 the number of clutters. Note that elements in 𝐳_𝑘 are not ordinal. For BOT, the relationship between 𝛽_𝑘 and 𝐱_𝑘 is specified to be 𝛽_𝑘𝑥=arctan^𝑘−𝑥_𝑜,𝑘𝑦_𝑘−𝑦_𝑜,𝑘+𝐰_𝑘, (8) where 𝐰_𝑘 is the measurement noise with known pdf and is independent of 𝐯_𝑘 in (1) and (𝑥_𝑜,𝑘,𝑦_𝑜,𝑘) the sensor’s position.

The measurement elements due to clutters are assumed to be uniformly distributed over the surveillance region, so the pdf 𝑝_𝑐(⋅)=1/𝑉, where 𝑉 denotes volume of the surveillance region. The number of clutters per scan, 𝑚_𝑘, is Poisson distributed with density parameter 𝑀_𝑐; that is, 𝑚_𝑘~Poission(𝑀_𝑐).

2.4. The BOT Initiation Problem

The information on target presence along with the target’s dynamic state (if target presence has been confirmed) is of interest here, while the only probe to get the desired information is the received measurements.

In Figure 1, we give an example showing the measurements generated in a BOT initiation case, where the target exists from the 8th second to the 32nd second. In the upper plot, the target-originated measurements are linked by a series of line segments, while the other points denote measurements generated due to clutters. The bottom plot represents the true scenario encountered in practice where there is no indication on the source of the measurements.

(a)

(b)

Figure 1 

Received BOT measurements versus discrete time index. The left plot indicates target originated measurements by crossing them by a line, while the right plot does not give any indication on the source of each measurement.

Figure 2 

The observer trajectory and the target trajectory in the simulation.

The task of track initiation just consists in taking sets of point observations from different time steps and linking together those observations that fit a desired model without any previous track estimates. As shown in Figure 1, we may draw many different lines to link many different sets of measurements, while there is only one true answer underlying these measurements. In general, such track initiation problem suffers from a combinatorial explosion in the number of potential tracks that must be evaluated.

As indicated by Figure 1, it is reasonable to input the measurements to the tracker algorithm in a batch mode, as latter observations can be used to benefit making decisions at previous time steps. That is why so many batch mode techniques are developed for track initiation, for example, the Hough-transform-based approach [6, 8, 19–21], the multiple-kd tree algorithm [22], and so on.

However, the problem of interest here is different. The focus is how to do BOT initiation in a sequential manner. We formulate this problem as follows. Given 𝑝(𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘), once 𝐳_𝑘+1 is observed, how to calculate 𝑝(𝐱_𝑘+1,𝜆_𝑘+1∣𝐙_𝑘+1) immediately, 𝑘=1,2,…. Actually, 𝑝(𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘) covers all the information about (𝐱_𝑘,𝜆_𝑘) that we are able to get based on 𝐙_𝑘. It is noted that, provided that sufficient statistics of 𝑝(𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘) can be obtained, there is no need to reload and process 𝐙_𝑘 for calculating 𝑝(𝐱_𝑘+1,𝜆_𝑘+1∣𝐙_𝑘+1) upon the arrival of 𝐳_𝑘+1.

3. Sequential Bayesian Estimation Framework for BOT Initiation

Here, the distribution of interest is the posterior 𝑝(𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘). We present a Bayesian estimation framework for computing this posterior distribution sequentially in a recursive manner. The recursion starts from an initial distribution (𝐱₀,𝜆₀)∼𝑝(𝐱₀,𝜆₀)∣𝐙₀, where 𝐙₀ is empty, since at the very beginning, there is no measurement element received. The successive posteriors are computed as follows: 𝑝𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘−1=𝑝𝐱_𝑘,𝜆_𝑘∣𝐱_𝑘−1,𝜆_𝑘−1𝑝𝐱_𝑘−1,𝜆_𝑘−1∣𝐙_𝑘−1×𝑑𝐱_𝑘−1,𝜆_𝑘−1𝑝𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘∝𝑝𝐳_𝑘∣𝐱_𝑘,𝜆_𝑘𝑝𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘−1. (9) Now let us examine individual parts that constitute the above integrals. First, consider the state transition part: 𝑝𝐱_𝑘,𝜆_𝑘∣𝐱_𝑘−1,𝜆_𝑘−1=𝑝𝐱_𝑘∣𝐱_𝑘−1,𝜆_𝑘−1,𝜆_𝑘𝑝𝜆_𝑘∣𝜆_𝑘−1. (10) Note that 𝑝(𝜆_𝑘∣𝜆_𝑘−1) is just determined by the birth/death Markov model (6). If 𝜆_𝑘=0, the target is absent, so 𝐱_𝑘 is undefined; otherwise, the pdf of 𝐱_𝑘 conditional on 𝐱_𝑘−1 and 𝜆_𝑘−1 is given by 𝑝𝐱_𝑘∣𝐱_𝑘−1,𝜆_𝑘−1,𝜆_𝑘=𝑝𝐱=1^𝑘∣𝐱_𝑘−1for𝜆_𝑘−1𝑝=1,_𝐵𝐱_𝑘for𝜆_𝑘−1=0, (11) where 𝑝_𝐵(⋅) denotes the initial pdf of the target on its appearance.

Now, we consider the calculation of the measurement likelihood. To begin with, introduce variable 𝑛_𝑘 to denote the number of measurement elements in 𝐳_𝑘. If the target is absent, that is, 𝜆_𝑘=0, then 𝑛_𝑘 just equals the number of measurements due to clutters, that is, 𝑚_𝑘; otherwise, 𝑛_𝑘=𝑚_𝑘+1. Now, we introduce the association variable 𝛾_𝑘. We use 𝛾_𝑘=𝑖,𝑖∈{1,…,𝑛_𝑘}, to denote the event that the 𝑖th element in 𝐳_𝑘 is target-originated, and the others are generated due to clutters. Additionally, 𝛾_𝑘=0 denotes that no elements in 𝐳_𝑘 are target-originated. Now, the measurement likelihood is specified to be 𝑝𝐳_𝑘∣𝐱_𝑘,𝜆_𝑘=^𝑛𝑘_𝑖=0𝑝𝐳_𝑘∣𝐱_𝑘,𝜆_𝑘,𝛾_𝑘=𝑖𝑝𝛾_𝑘=𝑖∣𝐱_𝑘,𝜆_𝑘, (12) where 𝑝𝐳_𝑘∣𝐱_𝑘,𝜆_𝑘,𝛾_𝑘=𝑝𝑧=𝑖_𝑖,𝑘∣𝐱_𝑘⋅𝑝_𝑐(⋅)^𝑚𝑘,if𝜆_𝑘𝑝=1,𝑖≠0,_𝑐(⋅)^𝑛𝑘,if𝜆_𝑘𝑝𝛾=0,𝑖=0,0,otherwise,_𝑘=𝑖∣𝐱_𝑘,𝜆_𝑘=1𝑛_𝑘,if𝜆_𝑘=1,𝑖=1,…,𝑛_𝑘,1,if𝜆_𝑘=0,𝑖=0,0,otherwise. (13)

4. Sequential Monte Carlo Implementation

In this section, SMC techniques are developed to approximately implement the Bayesian estimation framework presented in Section 3. The idea is to use the Monte Carlo sampling methodology to approximate the involved distributions by a set of weighted random samples. It is noted that convergence of SMC methods in dealing with nonlinear non-Gaussian Bayesian tracking problem has been proved [23–25].

To begin with, denote {𝐱^𝑛_𝑘−1,𝜆^𝑛_𝑘−1,𝑤^𝑛_𝑘−1}^𝑁_𝑛=1 as a random measure that approximates the posterior at 𝑘−1; namely, 𝑝(𝐱_𝑘−1,𝜆_𝑘−1∣𝐙_𝑘−1). Here, 𝑁 is the sample size; 𝑤^𝑛_𝑘−1 is the importance weight of the 𝑛th sample (𝐱^𝑛_𝑘−1,𝜆^𝑛_𝑘−1). In the following, an attempt is made to derive an evolution law for the current random set of particles to get a particle approximation to 𝑝(𝐱_𝑘,𝜆_𝑘∣𝐙_𝑘).

First, evolve the existence variable one time step further for each particle based on the birth/death Markov model (6). For example, if 𝜆^𝑛_𝑘−1=0, then the probabilities of 𝜆^𝑛_𝑘=1 and 𝜆^𝑛_𝑘=0 are 𝑃_𝑏 and 1−𝑃_𝑏, respectively, according to (4) and (6).

Next, we show how to generate state vector value 𝐱_𝑘 for each particle according to (10) and (11). Let us focus on particles that are associated with positive 𝜆_𝑘 values. They are active particles at time step 𝑘. These particles can be categorized into two classes:(i)newborn particles: this set of particles is characterized by the transition from 𝜆^𝑛_𝑘−1=0 to 𝜆^𝑛_𝑘=1;(ii)staying alive particles: this is a set of particles that continues to stay active with 𝜆^𝑛_𝑘−1=1 and 𝜆^𝑛_𝑘=1.According to (11), we draw state vector sample values 𝐱_𝑘 for the newborn particles from a given initial pdf 𝑝_𝐵(⋅) and draw state vector sample values for the staying alive particles from the dynamic transition prior density, which are determined by the target dynamic model (1). For the rest of particles that are associated with 𝜆_𝑘=0, just keep their state vector values undefined, since they are totally useless in the algorithm.

Now we see, for each particle, say the 𝑛th, its associated state vector value 𝐱^𝑛_𝑘 and existence variable value 𝜆^𝑛_𝑘 are obtained, while the importance weight 𝑤^𝑛_𝑘 needs to be updated when the new measurement 𝐳_𝑘 arrives at time step 𝑘. We update the importance weights according to the solid Bayesian formalism given by (9). Given 𝐱^𝑛_𝑘 and 𝜆^𝑛_𝑘, the (unnormalized) importance weight 𝑤^𝑛_𝑘 is just equivalent to the likelihood of 𝐳_𝑘. It is calculated according to (12). A normalization step is processed as follows: 𝑤^𝑛_𝑘=𝑤^𝑛_𝑘∑^𝑁_𝑘=1𝑤^𝑛_𝑘, (14) to make sure that the summation of the importance weights equals 1. Then, the posterior probability of target existence, that is, 𝑝(𝜆_𝑘=1∣𝐙_𝑘), is computed as below 𝑃_𝑘=^𝑁_𝑛=1𝜆^𝑛_𝑘𝑤^𝑛_𝑘, (15) which satisfies 0≤𝑃_𝑘≤1. The target is declared to be present if 𝑃_𝑘 is bigger than a given threshold, 0.6 is used in our simulation test. If the target is determined to be present, its state is estimated as follows: 𝐱_𝑘∣𝑘=∑_{𝑛∈{1,…,𝑁}∶𝜆}^𝑛_𝑘=1𝐱^𝑛_𝑘𝑤^𝑛_𝑘∑_{𝑛∈{1,…,𝑁}∶𝜆}^𝑛_𝑘=1𝑤^𝑛_𝑘. (16) To avoid the problem of degeneracy of SMC algorithm, that is, avoiding the situation that all but one of the importance weights are close to zero, a resampling procedure is performed if the effective number of particles is less than a given threshold 𝑁_thr. An estimate of the effective number of particles is computed as follows [26]: 𝑁_eff=1∑^𝑁_𝑛=1𝑤^𝑛_𝑘². (17) The resampling procedure can be summarized as follows: draw particles from the current particle set with probabilities proportional to their weights and then replace the current particle set with this new one. For more discussions on resampling methods in context of PF, refer to [23, 27, 28]. The above SMC algorithm is initialized at 𝑘=1 by drawing samples 𝜆^𝑛₁, 𝑛=1,…,𝑁 in accordance with the initial target existence probability 𝜇₁. The initialization of the active particles’ state vector is the same as for the newborn particles described above.

5. Simulation Study

In this section, the performance of the proposed PF algorithm for sequential BOT initiation is evaluated via simulations. The estimation accuracy of the proposed algorithm is compared to the Posterior Cramér-Rao lower bounds (CRLB) [29]. Two generally accepted approaches for sequential track initiation, namely, integrated probabilistic data association (IPDA) [9] and integrated multiple model probabilistic data association (IMMPDA) [30, 31], are also involved for performance comparison.

5.1. Simulation Setting

The scenario to be investigated involves BOT of a target from an observer, as shown in Figure 1. The observer travels at a fixed speed of 10 m/s and executes 2 maneuvers during the period under surveillance. The scanning process lasts 40 seconds in total. During the first 8 seconds, the target does not emit any energy outward, so there is no target-originated elements in the received measurements. During the period from the 8th to the 32nd second, the target moves and emits energies outward, so the received measurements include the target-originated elements and the tracker algorithm is expected to be able to detect the target. After the 32nd seconds, generation of target-originated elements is terminated in the received measurements. This setting is used to test whether a tracker algorithm under examination can terminate the track in time. The target motion is simulated according to (3) subjected to an amount of process noise with 𝑞_𝑠=0.1. The target’s 2D position and velocity are initialized to be (380,200) and (12.25,−12.25), respectively. The other parameters used in this simulation are listed in Table 1. An example show of the received measurements versus discrete time index is presented in Figure 1.

Table 1 

Parameter setting for simulation.

5.2. Performance Comparison on Target Detection

Detection performance of the involved algorithm is investigated via Monte Carlo simulations with 100 independent runs of each algorithm. The outputted probabilities of target existence at each time step are averaged and then plotted in Figure 3. As is shown, the proposed SMC algorithm beats both IPDA and IMMPDA remarkably in the rate of convergence.

Figure 3 

The true and the estimated probabilities of target existence versus time index.

Quantitative comparison is also conducted via hypothesis testing. Hypotheses 𝐻₀ and 𝐻₁ are defined to represent the events “target being absent” and “target being present,” respectively. At each time step 𝑘, four possible cases may happen for each algorithm:(1)𝐻₀ is true and the algorithm chooses 𝐻₀,(2)𝐻₀ is true and the algorithm chooses 𝐻₁,(3)𝐻₁ is true and the algorithm chooses 𝐻₀,(4)𝐻₁ is true and the algorithm chooses 𝐻₁.We select probability of false alarm, denoted by 𝑃^FA, and probability of detection, denoted by 𝑃^𝐷, as the performance metrics. Let 𝐶_𝑖 count the number of times of case 𝑖 appearance in the simulation test, 𝑖∈{1,2,3,4}. Then 𝑃^FA=𝐶₂/(𝐶₁+𝐶₂) and 𝑃^𝐷=𝐶₄/(𝐶₃+𝐶₄). So a good algorithm is expected to yield bigger 𝑃^𝐷 and smaller 𝑃^FA. 𝑃^𝐷 and 𝑃^FA are calculated over 100 independent runs of each algorithm and the final result is listed in Table 2. It is demonstrated again that the proposed SMC algorithm is superior to IPDA and IMMPDA in detection performance.

Table 2 

Detection performance comparison.

5.3. Examination on Estimation Accuracy

The posterior Cramer-Rao lower bounds (CRLB) [29] are used to indicate the best estimation accuracy an algorithm can achieve. The performance metric, namely, the root-mean square (RMS) position error is defined as follows: RMS_𝑘=∑^𝑀_𝑖=1𝑥^𝑖_𝑘−𝑥_𝑘²+𝑦^𝑖_𝑘−𝑦_𝑘²𝑀, (18) where (𝑥_𝑘,𝑦_𝑘) denotes the target’s 2D position at time step 𝑘 in the simulations, (𝑥^𝑖_𝑘,𝑦^𝑖_𝑘) the estimation of (𝑥_𝑘,𝑦_𝑘) outputted in the 𝑖th test of the algorithm under consideration, and 𝑀 the number of independent tests under consideration. Define 𝐽⁻¹_𝑥,𝑘, 𝐽⁻¹_𝑦,𝑘 as the 2D position elements of the inverse information matrix for the problem at hand; the corresponding CRLB for the metric given by (18) is calculated as below: CRLBRMS_𝑘=𝐽⁻¹_𝑥,𝑘+𝐽⁻¹_𝑦,𝑘. (19) Here, we only consider the scanning time interval between the 8th and the 32nd time steps, when the target is truly present and the proposed algorithm successfully detects it. Figure 4 shows the RMS position error of the proposed algorithm in comparison with the calculated CRLB.

Figure 4 

RMS position error in comparison with the posterior CRLB.

6. Conclusions

This paper addresses the problem of sequential BOT initiation in the context of sonar and radar applications. A sequential Bayesian estimation framework is developed for this problem. This theoretical framework addresses measurement origin uncertainty and target existence uncertainty jointly via solid Bayesian formalism. An SMC approximate algorithm is derived under this framework and its performance is evaluated via simulations. It is shown that the proposed algorithm provides a remarkable performance improvement in target detection, compared with the commonly used PDA based methods. The proposed algorithm also gives accurate estimation of the target’s state, as indicated by a comparison with the posterior CRLB.

Acknowledgments

The authors appreciate the insightful comments and suggestions provided by the anonymous reviewers. This work was supported by the National Natural Science Foundation (NSF) of China (Grants nos. 61302158, 61172166, 61100135, 61170065, and 61373017), Provincial Science and Technology Plan (NSF) of Jiangsu province (Grant no. BK20130869), and Natural Science Research Project for Colleges and Universities in Jiangsu province (Grant no. 13KJB520019).