A Frequency Selection Method Based on the Pole Characteristics

Abstract

Due to the heavy jamming band of high frequency, frequency selecting strategies are serious issues for the system designed to achieve its best performance. Pole is independent of the direction and polarization of the incident wave, but the residue corresponding to the pole is related to the direction and polarization of the incident wave. And the value of residue is proportional to the value of the pole. This paper chooses the frequency which can maximize the residue in the high-frequency band as the optimal frequency for accurate extraction. The simulation result of a large number of ship targets shows remarkable rise in average recognition rate by using this method, compared with the average recognition rate of randomly selected frequency.

1. Introduction

With the extensive application of high-frequency radar, radar frequency selection has become an important research topic [1–5]. The propagation characteristics of the high-frequency electromagnetic waves are not only a function of the distance but also a function of frequency [4]. In the case of the same distance, the higher the frequency, the faster the attenuation of electromagnetic waves. The frequency selection will seriously affect the performance of the high-frequency surface wave radar on target detection. It is necessary to optimize the selection of operating frequency of the high-frequency radar in order to get excellent detection performance of high-frequency radar. This paper uses pole characteristic to select frequency, explore criteria for frequency selection optimization, and get the selected frequency optimization method. Finally, the result of simulation for a variety of ship target detection verified the correctness and validity of the method.

2. Frequency Optimization Criterion

To get an optimal operating frequency, the first thing to consider is the frequency optimization criterion. A high-frequency radar operating frequency can be obtained through scientific quantitative calculation in accordance with this criterion so that the high-frequency radar can achieve optimal detection performance. Radar is able to find the goal, that is to say, the average correct recognition rate to meet certain requirements. Paper [6] puts forward that radar working frequency optimization selection is investigated for an optimization selection method based on the minimum classification error rule in order to improve airplane target identification performance and that MER can change to the principle of average correct recognition rate maximum. Optimal operating frequency is the one which makes the average correct identification rate maximum. This also can be called as the principle of average correct recognition rate maximum. The average correct recognition rate can be got from dividing the number of correctly identified test by the total number of randomized trials.

Pole is independent of the direction and polarization of the incident wave, but the residue corresponding to the pole is related to the direction of polarization of the incident wave. And the value of residue is proportional to the value of the pole. Considering the above two factors, select the frequency which corresponds to the maximum value of residue as the optimal frequency.

3. Calculation of Residue

Matrix Pencil Method [7] is proposed on the basis of Pencil of Function [8] and has a higher computational efficiency. It only needs two Hankel matrices constructed by time-domain data sequence of the target. Then seeking the generalized eigenvalue of these two Hankel matrix we can calculate the poles of target. Because Matrix Pencil Method needs not to calculate the root of the characteristic equation as Prony method, Matrix Pencil Method has high antinoise performance. To improve the anti-noise performance of Matrix Pencil Method, Sarkar first applied Matrix Pencil Method in matrix singular value decomposition [9] and gave a method to decide the number of poles, through the matrix low rank approximation which could well inhibit the effects of noise, and the efficiency of the algorithm is improved greatly. But Matrix Pencil Method depends on the selection of relevant parameters; Sarkar studies the beam parameters selection [10] and the method for determining the number of target poles [9]. Matrix Pencil Method has high computational efficiency in certain SNR and good operation steadiness, which is very suitable for the extraction of the target pole.

From previous research we can see that the properties of Matrix Pencil Method are superior than others among a large number of pole extraction methods [6]. But there are some disadvantages, such as the unknown number of poles that will to a large extent affect the matrix method of pole extraction results. We use the improved the overall minimum square Matrix Pencil Method to extract the pole of complex targets [5, 9].

3.1. Principles of Matrix Beam Method

In 1971, after a lot of electromagnetic pulses tests, Baum found that, when objects are excited by electromagnetic pulse, their late response is a series of superposition of attenuation sinusoidal oscillation, and that attenuation factor and oscillation frequency composed the complex natural resonant frequency, which is the pole of target. Baum puts forward the theory of singularity expansion method on the basis of the above phenomena [11, 12].

The late transient response in time domain of target can be represented by a finite sum of damped sinusoids: 𝑦(𝑡)=^𝑀_𝑖=1𝑅_𝑖𝑒^𝑠𝑖𝑡, (1) where 𝑅_𝑖=𝛼_𝑖+𝑗𝛽_𝑖, 𝑆_𝑖=𝜎_𝑖+𝑗𝜔_𝑖; 𝑅_𝑖 are complex residue; 𝑆_𝑖 is pole; 𝜎_𝑖, 𝜔_𝑖 is, respectively, attenuation factor and resonance pulsation. The discrete form after sampling is 𝑦(𝑛)=^𝑀_𝑖=1𝑅_𝑖𝑧^𝑛_𝑖, (2) where 𝑛=1,2,…,𝑁; 𝑁 are sampling points. To define the matrix disaggregatedly 𝐘=𝐲_𝐿,…,𝐲₁,𝐲₀,𝐘₁=𝐲_𝐿,…,𝐲₂,𝐲₁,𝐘₂=𝐲_𝐿−1,…,𝐲₁,𝐲₀, (3)(4)(5) where 𝐲_𝑙=[𝑦(𝑙),𝑦(𝑙+1),…,𝑦(𝑙+𝑁−𝐿−1)]^𝑇𝑙=1,2,…,𝐿; 𝐿 is called beam parameters.

Equations (4) and (5), can be written as 𝐘₁=𝐙_𝐿𝐑𝐙₀𝐙_𝑅,𝐘₂=𝐙_𝐿𝐑𝐙_𝑅, (6)(7) where, 𝐙_𝐿=𝑧11⋯1₁𝑧₂𝑧_𝑀𝑧⋮⋱⋮_{1𝑁−𝐿−1}𝑧₂^{𝑁−𝐿−1}⋯𝑧_𝑀^{𝑁−𝐿−1},𝐙_𝑅=𝑧₁^𝐿−1𝑧₁^𝐿−2𝑧⋯1_2𝐿−1𝑧₂^𝐿−21𝑧⋮⋱⋮_𝑀𝐿−1𝑧_𝑀^𝐿−2,𝑅⋯1𝐑=diag₁,𝑅₂,…,𝑅_𝑀,𝐙₀=diag𝑧₁,𝑧₂,…,𝑧_𝑀. (8)

Define the matrix pencil 𝐘₁−𝜆𝐘₂=𝐙_𝐿𝐑𝐙₀−𝜆𝐈𝐙_𝑅, (9) If beam parameters 𝐿 satisfie 𝑀≤𝐿≤𝑁−𝑀, then the rank of matrix pencil is 𝑀. If 𝜆=𝑧_𝑖, then the rank of matrix pencil decreased to 𝑀−1. In other words, {𝑧_𝑖} is generalized eigenvalue of matrix {𝑌₁,𝑌₂}.𝐘₁𝑟_𝑖=𝑧_𝑖𝐘₂𝑟_𝑖, (10)𝑟_𝑖 are eigenvectors corresponding to eigenvalue 𝑧_𝑖. Equation (10) can also be written as 𝐘^†₂𝐘₁−𝑧_𝑖𝐈𝑟_𝑖=0, (11) where superscript “†” represents the generalized inverse of the matrix. It can be seen from (11) that it is possible to obtain the required target pole when eigenvalue of 𝐘^†₁𝐘₂ is a known value.

And then, target complex residue can be obtained through the least squares method: 𝐱=𝐀^𝑇𝐀⁻¹𝐀^𝑇𝐛, (12) where 𝑧𝐀=11⋯1₁𝑧₂𝑧_𝑀𝑧⋮⋱⋮_1𝑁−1𝑧₂^𝑁−1⋯𝑧_𝑀^𝑁−1,𝐱=𝑅₁,𝑅₂,…,𝑅_𝑀^𝑇,𝐛=𝑦₀,𝑦₁,…,𝑦_𝑁−1^𝑇. (13)

3.2. The Choice of Pole Number

Advanced transient response of fully conducting targets can expand into an infinite number of decaying sinusoidal vibrations and forms. And formula (1) constrained the complex attenuation frequency on the 𝑀. It is One of the important questions about extracting pole of complex targets is how to select the appropriate 𝑀 [13]. If the 𝑀 value is selected too small, it will miss the target’s true pole, and the pole will have a large deviation; if the 𝑀 value is selected too large, it will not only produce false pole but also will lead to true pole position moved.

Pole number of simple targets within a certain frequency band can be easy to get; but for the aircraft and ship target, the pole number 𝑀 contains in the signal is unknown. For this issue, this paper uses the selection method based on time-domain signal reconstruction error minimized to decide the number of poles [5, 14].

In order to get better results of pole extraction, we can apply the singular value decomposition technique to pole extraction problem with matrix beam methods. Let us suppose formula (1) is accurate; if we ignore deviation and noise which come from modeling and calculations, matrix 𝐘 will become a full rank matrix.

By singular value of formula (3), we can get this: 𝐘=𝐔𝐔^′Σ00Σ^′𝐕^𝐻𝐕^′𝐻,Σ=diag𝜎₁,𝜎₂,…,𝜎_𝑀,Σ^′=𝜎_𝑀+10⋯00𝜎_𝑀+20⋮⋱⋮00⋯0_{(𝑁−𝐿−𝑀)×(𝐿+1−𝑀)}. (14)𝜎₁,𝜎₂,…,𝜎_𝑀 is the biggest 𝑀 singular values of 𝐘. If there is no noise, 𝜎_𝑀+1,𝜎_𝑀+2,… will be zero. Take low rank of 𝐘 to approximation of matrix 𝐘: 𝐘=𝐔Σ𝐕^𝐇. (15) And we can remove the last column of 𝐘, then we will get matrix 𝐘₁, and if we remove the first column, we will get matrix 𝐘₂ like this 𝐘=𝐘₁𝐲_𝟎=𝐲_𝐋𝐘₂. (16)𝐘₁, 𝐘₂ singular values can be decomposed as 𝐘₁=𝐔Σ𝐕^𝐇_𝟏,𝐘₂=𝐔Σ𝐕^𝐇_𝟐. (17) removes the last column of matrix 𝐕^𝐇 will get 𝐕^𝐇_𝟏, and if remove the first column, we’ll get 𝐕^𝐇_𝟐.

Theorem [15] If matrix 𝐘’s rank greater than or equal 𝑀, then there exists the only matrix of the same dimensions 𝐘, and the rank is equal to 𝑀, which makes the norm of ‖𝐘−𝐘‖_𝐹 Frobenius minimum. The extent of best approximation is described by ‖𝐘−𝐘‖_𝐹=∑^𝑁−𝐿_𝑖=𝑀+1𝜎_𝑖, where ‖⋅‖_𝐹 indicates Frobenius norm.

Theorem explains that the use of 𝐘 is the optimal approximation to 𝐘 under F-norm in theory. We will get 𝐘₁ and 𝐘₂ through 𝐘; then the eigenvalue of 𝐘^†₁𝐘₂ will be got. This method can effectively inhibit the effects of noise on pole extraction. This singular value decomposition technique applied to a method of matrix beam method is called total least square matrix beam method [5].

Theorem in the assumption of formula (1) is accurate and on the condition of 𝑀 is known. In the actual pole extraction, 𝑀 is unknown, and it cannot be determined only through the distribution of singular values 𝜎₁,𝜎₂,…,𝜎_𝑀,𝜎_𝑀+1,𝜎_𝑀+2,…. Because of the distribution of singular value, neither 𝑀 is rendered when you select a value, then the singular value is far smaller than the preceding singular value, nor 𝑀 appears when you select a value, followed by singular value of size are essentially the same. Therefore you need to focus on pole number 𝑀 of selection method.

It can judged by pole extraction results to decide the parameters selection is good or not. We will use the deviations between theoretical and extraction as criteria to constrain parameters selection and to get the best parameters 𝑀. But the most target theoretical pole is unknown, so it needs to give a reasonable and actionable criterion.

If the target number of pole 𝑀 is known, pole 𝑧_𝑖 can be achieved by matrix beam method, and the complex residue of target 𝑅_𝑖 can be got by arithmetic expression (4). It also can reconstruct the time domain signal 𝑦_rec(𝑘) of target, then we will use the error between the reconstructed domain signal 𝑦_rec(𝑘) and domain signal 𝑦_cal(𝑘)𝑆(𝑀)=^𝑁_𝑘=1𝑦_rec(𝑘)−𝑦_cal(𝑘)². (18) As a criterion min_𝑀^𝑁_𝑘=1𝑦_rec(𝑘)−𝑦_cal(𝑘)² (19) Choosing the appropriate parameter 𝑀_opt.

Cause 𝑀 can meet 𝑀∈𝐍,1≤𝑀≤min{𝐿,𝑁−𝐿} (20)

And usually target scattering data length 𝑁 is not too long; the value of 𝐿 is generally around 𝑁/3. So the value of 𝑀 is an integer from 1 to 𝑁/3. We can directly choose all parameters 𝑀 within the interval range to calculate 𝑆(𝑀) select minimal parameter 𝑆(𝑀) which is corresponding to 𝑀 as the most appropriate parameters 𝑀_opt. Figure 1 shows that the ship target error 𝑆(𝑀) and pole number determine the results of process variation of 𝑆(𝑀) with parameter 𝑀. It can be seen that as 𝑀 is increasing, errors it first shows shock reducing, then minimum is reached, and at last it will be larger again.

Figure 1 

The error is changing with the value of 𝑀 in different.

We can see from Figure 1 that in this paper, the best number of pole is 𝑀 = 46. This approach can reduce the impact because models are inaccurate, and there are deviations in the data and the noise, extracting the precise targets pole.

4. Simulation

4.1. The Rationality of the Residue Selected Frequency

To prove the rationality of the residue selected frequency, we need to verify that the residue is a main contribution to RCS. Because the residue is obtained from the late response, we need to compare the contribution of the early and late response to the RCS for the first step.

4.1.1. How to Recover RCS by Using Late Response Only

Target time-domain transient response and frequency response are the Fourier transform for each other. If we can recover the time-domain transient response, the frequency domain response of the target can be obtained naturally. Target time-domain transient response is composed of early and late responses, which can be written as 𝑦_all(𝑡)=𝑦_ear(𝑡)+𝑦_lat(𝑡), (21) where 𝑦_ear(𝑡), 𝑦_lat(𝑡) represent the early and late responses. And 𝑦_lat(𝑡) is the 𝑦(𝑡) in formula 𝑦(𝑡)=∑^𝑀_𝑖=1𝑅_𝑖𝑒^𝑠𝑖𝑡; we can rewrite formula 𝑦(𝑡)=∑^𝑀_𝑖=1𝑅_𝑖𝑒^𝑠𝑖𝑡 as 𝑦_lat(𝑡)=^𝑀_𝑖=1𝑅_𝑖𝑒^𝑠𝑖𝑡. (22)

In theory, no matter what form of target transient response 𝑦_all(𝑡) is in formula 𝑦_all(𝑡)=𝑦_ear(𝑡)+𝑦_lat(𝑡), it can be expanded into an infinite sum number of complex exponential 𝑦_all(𝑡)=^∞_𝑖=1𝑅_𝑖𝑒𝑠̃_𝑖𝑡, (23) where 𝑅_𝑖 is the coefficient of expansion corresponding to exponential term 𝑒𝑠̃_𝑖𝑡. Basis function set {𝑠₁̃,𝑠₂̃,…} definitely contains the target pole set {𝑠₁,𝑠₂,…,𝑠_𝑀}. Although using only the pole set {𝑠₁,𝑠₂,…,𝑠_𝑀} to expand the transient response 𝑦_all(𝑡) may exist certain error, but that error is admissible. 𝑦_all(𝑡)≈^𝑀_𝑖=1𝑅_𝑖𝑒^𝑠𝑖𝑡, (24) where 𝑅_𝑖 can be called the new residue corresponding 𝑆_𝑖 to different from 𝑅_𝑖.

4.1.2. The Comparison Simulation

Here we get the original RCS by using software FEMAP to model the ship target and FEKO to calculate the RCS. The target we choose 49.05 meters long, 12.40 meters high, and 56.39 meters wingspan airplane, the model of airplane, is shown in Figure 2. The comparison result is shown in Figures 3 and 4, red line represents recovered RCS only by late time response and blue line shows the original RCS, which proves that the late response is a very important contribution to RCS. The above conclusion implies that the residue calculated from late response is also very important to RCS. Next, we will test and verify the above conclusion by estimating the contribution of residue to RCS.

Figure 2 

Airplane model using FEMAP (𝐿 × 𝑊 × 𝐻 (m): 49.05 × 56.39 × 12.40).

(a) Forward

(b) Backward

Figure 3 

The comparison between original RCS and the RCS recovered only by late time response in frequency domain. (a) The angle between the incident and the bow is 0°. (b) The angle between the incident and the bow is 180°.

(a) Forward

(b) Backward

Figure 4 

The comparison between original RCS and the RCS recovered only by late time response in time domain. (a) The angle between the incident and the bow is 0°. (b) The angle between the incident and the bow is 180°.

To make the simulation contain common complex goal, this time I use software FEMAP and FEKO to model one ship and calculate its RCS, the model is shown in Figure 5. Then calculating the residue according to the method I described in the paper, incident angle ranges from 0° to 180° at a 10° interval. The frequency selection is from 0.2 MHz to 20 MHz at a 0.2 MHz interval. The number of the pole is 46. The result is shown in Figure 6, We can see that the residue changes with frequency and angle. It shows peak appeared when angle is 90° and frequency is 9.6 MHz.

Figure 5 

Ship model using FEMAP Destroyer (𝐿 × 𝑊 × 𝐻 (m): 154.0 × 17.0 × 6.0).

Figure 6 

The ship target residue varies with frequency and angle where the 𝑥-axis is frequency and 𝑦-axis is residue. The red line represents the residue changing trend with frequency when angle is 0°, blue line represents 90°, and green line is 180°.

We can set that maximum residue as zero and then recover the RCS. Comparing the variation of RCS before and after. Figure 7 shows that when maximum residue was replaced by 0, RCS declined to 8 dB–10 dB; the above simulation proves the influence of residue on the value of RCS. The simulation results provide a theoretical basis for using residue to select frequency.

Figure 7 

Comparison of RCS between original and residue = 0.

4.2. Compare the Average Correct Recognition Rate

4.2.1. Nearest Neighbor Rule Classification

Modeling 𝐸 types different target by software FEMAP, discrete the high-frequency band into 𝐾 frequency points with a fixed interval, obtain the RCS by FEKO calculating in all angle, we can get a target database in high-frequency band: 𝐃=𝐃₁𝐃₂⋯𝐃_𝐸=𝐷₁₁𝐷₁₂⋯𝐷_1𝐾𝐷₂₁𝐷₂₂⋯𝐷_2𝐾𝐷⋮⋮⋱⋮_𝐸1𝐷_𝐸2⋯𝐷_𝐸𝐾. (25)

Extracted 𝑘 frequency points from 𝐾 frequency points, are equivalent to extracted 𝑘 column vectors from 𝐾 column vectors in above matrix, to constitute a submatrix: 𝐃=𝐃₁𝐃₂⋯𝐃_𝐸^𝑇=𝐷₁₁𝐷₁₂⋯𝐷_1𝑘𝐷₂₁𝐷₂₂⋯𝐷_2𝑘𝐷⋮⋮⋱⋮_𝐸1𝐷_𝐸2⋯𝐷_𝐸𝑘. (26)𝐃_𝑖 is a radar target feature vector; name 𝐃 as radar target feature space composite of 𝐸 types target feature vector.

Obtain the RCS on 𝑘 frequency points by actual measurement, write it as 𝐃=𝐷₁,𝐷₂,…𝐷_𝑘, (27) name it radar target measurements eigenvectors.

According to the existing research Ksienski and Lin show that the nearest neighbor classification method for the class of complex radar target aircraft has better recognition results [11]. The Euclidean distance between measuring characteristic feature and the 𝐸 types target in target feature space is 𝑑_𝑚=𝐃−𝐃_𝑚=^𝑘_𝑗=1(𝐷_𝑗−𝐷_𝑚𝑗)^21/2. (28) According to nearest neighbor classification, if 𝑑_𝑚<𝑑_𝑛,𝑚,𝑛=1,2,…𝐸;𝑛≠𝑚 (29) target is 𝐸 type.

4.2.2. Simulation

In order to study the effect of the maximum number residue frequency-selective method of ship target also verifies that the maximum residue frequency selection method can effectively improve the recognition effect. I do a large number of random tests to verify the correctness of the theory and the effectiveness of the method.

The simulation selected 6 typical classes of ship target as a target to be classified from the library of high-band ship target characteristics; the type of ship targets and basic dimensions are shown in Table 1, Including aircraft carriers, cruisers, destroyers, frigates, hunting submarines, missile boats, and different types of ships. Here added noise to calculated target RCS to simulate the measured data, cause many of them are secret we cannot get the measured data. Here we can name this RCS calculated RCS with noise to different from below. Assume that select the 𝑚 target, add Gaussian white noise with zero mean; mean square deviation is 𝜎, and signal-to-noise ratio is calculated using 1SNR=𝑘^𝑘_𝑗=1𝐷_𝑚𝑗𝜎². (30)

Table 1 

Ship targets size.

Here does not consider the situation that the target is not included in the target class to be classified.

1000 randomized tests are conducted to study the changes of the correct identification rate in different signal-to-noise ratios with the different frequency selection method. Figure 8 shows average correct recognition rate comparison between maximum residue frequency selection method and random frequency selection method of 6 goals. Apparently, the average correct recognition rate of maximum residue selection frequency method is better than random frequency selection method. The former method is about six percent higher than the latter. The average correct recognition rate of maximum residue selection frequency method can be as high as 0.8.

Figure 8 

Average correct recognition rate changes with the signal-to-noise ratio by different frequency selection methods for 6 classes of target.

To make the simulation much more closer to real conditions, we need to consider that there may be deviations between the measurement RCS and calculated RCS from simulation. That deviations may come from the bug of software or the noise that affects the measurement of RCS. In this simulation, a larger amount of uniformly distributed random value is added to calculated RCS as that deviation. Taking into account that the RCS deviation may occur in the reality as 𝑥 dB, to 𝑚 target’s 𝑗 frequency, we add uniformly distributed random value in the range of [−Δ𝐷_𝑚𝑗,Δ𝐷_𝑚𝑗], satisfying 𝑥=10⋅log(Δ𝐷_𝑚𝑗/𝐷_𝑚𝑗). At the same time, Gaussian white noise we discuss in Figure 6 is still added (calculated RCS with noise). Figure 9 shows the recognition results of 6 ship targets at two operating frequencies 13.4 MHz and 15 MHz when taking consideration of the deviations between the measure RCS and calculated RCS.

(a) The measuring deviation is ±1 dB

(b) The measuring deviation is ±2 dB

Figure 9 

The average correct identification rate changes trend of 6 targets with different SNR when deviations exist. (a) The measuring deviation is ±1 dB. (b) The measuring deviation is ±2 dB.

From Figure 8, Compare the average correct recognition rate of the above method and the random frequency selection method when angle is 90°; average correct identification rate of my method is 82.76% when SNR is 13 dB, while the average correct identification rate of the random frequency selection mode is 61.33% when SNR is 13 dB, which verifies the effectiveness of the method.

In Figure 9, we can see that the average correct recognition rate declined quickly when measuring deviation becomes larger. When measuring deviation is ±1 dB, the average correct recognition rate can reach 0.87. When measuring deviation is ±2 dB, the average correct recognition rate can reach 0.68. Comparing Figure 9 with Figure 8 we can that the measuring deviation does have an influence on the average correct recognition rate. But it can be overcome by adding more frequency. The simulation of 4 frequencies 13.4 MHz, 15 MHz, 16.6 MHz, and 20 MHz is shown in Figure 10, from which we can see that the average correct recognition rate can reach 0.85.

Figure 10 

The measuring deviation is ±2 dB, and the frequencies are 4: 13.4 MHz, 15 MHz, 16.6 MHz, and 20 MHz.

From above comparison we can get the conclusion that the measuring deviation has a great influence on recognition rate. But this problem can be overcome by adding more frequency even when the measuring deviation exists, and the average recognition rate of the method adopted in this paper is significantly larger than random selected frequency method, which indicated that the robustness of the method is absolutely better than random selected frequency method. And the method is feasible.

5. Conclusion

This paper studies the frequency selection method in high-frequency band radar target recognition. Based on the average correct recognition rate maximization principle, considering the residues corresponding to pole is related to the incident wave direction and polarization; the value of residue is proportional to the value of the pole, a maximum residue frequency selection method was proposed, and satisfactory recognition result is obtained with the simulation experiment on multiclass ship targets. Considering the difference between real measured RCS, and software calculated RCS, we add measurement deviation in software calculated RCS, and the ship target recognition effect under the conditions of different measurement deviations is analyzed. This verifies the correctness and validity of the frequency optimization method.

Acknowledgments

This work was supported in part by the Harbin Institute of Technology (HIT) and in part by the National Natural Science Foundation of China (NSFC) under the Project no. 61171182.