Abstract
Channel state information (CSI) is important to improve the performance of wireless transmission. However, the problems of high propagation path loss, multipath, and frequency selective fading make it difficult to obtain the CSI in broadband millimeter-wave (mmWave) system. Based on the inherent multidimensional structure of mmWave multipath channels and the correlation between channel dimensions, mmWave multiple input multiple output (MIMO) channels are modelled as high-order parallel profiles with linear dependence (PARALIND) model in this paper, and a new PARALIND-based channel estimation algorithm is proposed for broadband mmWave system. Due to the structural property of PARALIND model, the proposed algorithm firstly separates the multipath channels of different scatterers by PARALIND decomposition and then estimates the channel parameters from the factor matrices decomposed from the model based on their structures. Meanwhile, the performance of mmWave channel estimation is analysed theoretically. A necessary condition for channel parameter estimation is given based on the uniqueness principle of PARALIND model. Simulation results show that the proposed algorithm performs better than traditional compressive sensing-based channel estimation algorithms.
1. Introduction
New generation wireless communications have tremendous requirements in high data rates and large connections. However, the severe shortage of spectrum resources has gradually become one of the main problems limiting its rapid development [1]. Millimeter wave (mmWave) has been introduced as one of the key technologies due to its advantage of massive available bandwidth, providing a foundation for satisfying the rapidly increasing requirements of wireless communications [2]. However, compared to the conventional microwave frequencies, inherent high signal attenuation and propagation loss are suffered in wireless communication when mmWave frequency is exploited. To compensate for such a loss, multiple input multiple output (MIMO) technique can be used to enhance the performance of signal transmission. Therefore, mmWave MIMO is a promising key technology in future wireless communication systems [3].
Channel state information (CSI) is essential for mmWave MIMO systems to improve the performance of wireless transmission. In wireless communication systems, channel information is usually obtained by channel estimation technologies. However, comparing to traditional wireless channels, the complexity, high propagation path loss, multipath, and Doppler of mmWave channels make it difficult in CSI estimation [4]. Meanwhile, lots of channel parameters are needed to be estimated due to the large number of antennas used in the system, leading to much more training overhead. Therefore, it is challenging to estimate the CSI in mmWave MIMO systems.
Since the path attenuation of mmWave channel is more serious than that of traditional microwave channel, the mmWave channel impulse response can be approximated by the paths with major energy. Therefore, mmWave channel is approximately sparse in angle or/and time domains and can be represented by the parameters of limited paths [5]. Based on this observation, mmWave channel can be uniformly quantized in angle domain by using the overcompleted dictionaries of angles of departure (AoDs) and angles of arrival (AoAs) [6]. By using the angle dictionaries as the sparse basis, mmWave channels can be formulated as a virtual channel model with sparse structure. Then, channel estimation can be formulated as a sparse recovery problem based on compressive sensing (CS) theory and can be solved by CS-based sparse recovery algorithms, such as orthogonal matching pursuit (OMP) algorithm. Different from [6], nonuniformly quantized dictionaries of AoD and AoA are adopted to formulate the virtual sparse version of mmWave channels in [7]. Nonuniformly quantization of angles can decrease the correlation between atomics in angle dictionaries, which may improve the performance of channel estimation in CS-based recovery algorithm. However, due to the use of nonuniform quantization method, the intervals between some quantized angles are relatively large, which may lead to more quantization error when the original AoDs and AoAs of paths fall into these intervals. To reduce the computation complexity, [8] proposed an adaptive CS-based algorithm for channel parameters estimation based on the virtual channel model of [6]. Similarly, [9] proposed a low-complexity channel estimation algorithm for mmWave system based on the sparsity of mmWave channels. However, different from [8], the sparsity representation of mmWave channels is formulated along two dimensions of AoA and AoD in [9], and the channels are estimated by separable compressive sensing method with low computation complexity.
The bandwidth of wireless transmission is increased significantly to meet the requirements of high rate and capacity in new generation mobile systems. Therefore, the problem of channel estimation for broadband mmWave systems has been studied recently. Based on the sparse structure of mmWave MIMO channel, [10] formulated the estimation of frequency selective mmWave channel as a sparse recovery problem based on CS theory and proposed a CS-based algorithm to estimate the channel in time domain, which reduces the pilot overhead. Since frequency selective fading is usually suffered in broadband wireless system, orthogonal frequency division multiplexing (OFDM) technology is often used in downlink transmission. [11] proposed a two-step channel estimation in time domain for broadband mmWave OFDM system. The algorithm firstly estimates the equivalent channels including beamforming, array response, and pulse shaping by least square method. Then, the channels are recovered by OMP algorithm with low complexity. Different from [10, 11], [12] proposed a hybrid algorithm to estimate channel for mmWave MIMO OFDM system by exploiting the sparse structure of the channel in angle domain, which can reduce power leaking in downlink transmission. [13] considers the downlink channel estimation and beamforming in mmWave OFDM system jointly. The beamforming training sequence is designed to construct measurement matrix for channel compression based on the sparse scattering characteristics of mmWave channels. By using the given measurement matrix, the CS-based algorithm can estimate the channel with low training overhead.
Above studies on mmWave MIMO channel estimation mainly exploit the sparse structure of channels in angle or time domain, and the channels are estimated by sparse recovery algorithm based on CS theory. In fact, the original mmWave channels are not strictly sparse in general. The CS-based algorithms need to transform the original channel matrix to the virtual sparse channel by using sparse dictionaries. The parameters of virtual channel are then recovered for channel estimation by exploiting its sparse structure. However, the sparse dictionaries are usually generated by the quantized angles or time slots, so that quantization error may be suffered between the quantized and the real parameters of propagation paths, which may decrease the performance of channel estimation. Meanwhile, previous studies usually focus on the scenarios of mmWave channels with single scattering path. However, mmWave signals may suffer from scattering multipath propagation during transmission in some scenarios, which brings challenges to the estimation of channel parameters. To solve these problems, this paper investigates the channel estimation problem of broadband mmWave system and proposes a new channel estimation algorithm based on high-order PARALIND model. The main contributions of this paper are list as follows: (1)The broadband mmWave channels are modelled as high-order PARALIND by exploiting the inherent multidimensional structure and the correlation between different dimensions of multipath channels. Due to the structural property of PARALIND model, the multipath channels of different scatterers can be separated by PARALIND decomposition. Then, the channel parameters can be estimated from the factor matrices of the model based on their structure with low computation complexity(2)The basic of proposed algorithm is not the sparse but multidimensional structure of the mmWave channels. Angle dictionaries are not used in channel modelling, so that no quantization error is suffered during channel estimation. Therefore, the accuracy of the proposed algorithm in channel estimation is higher than that of traditional CS-based algorithm(3)A necessary condition for channel parameter estimation is given based on the uniqueness property of PARALIND model. The performance of mmWave channel estimation is analysed theoretically, which can provide some theoretical supports for channel parameters estimation in broadband mmWave system
The rest of this paper is structured as follows. Section 2 presents the system model. Some preliminaries for PARALIND model are introduced in Section 3. In Section 4, we propose a new channel estimation algorithm based on PARALIND model. Section 5 presents some discussion on the proposed algorithm. An identifiability condition for channel parameters estimation is also given in this section. Section 6 gives simulation results for performance evaluation. The conclusions are drawn in the last section.
2. System Model
Consider the broadband mmWave MIMO downlink transmission system, where and antennas are used at transmitter and receiver. To overcome the impact of frequency selective channel, OFDM technology is used, where the band is divided into subcarrier for signal transmission. In order to improve the performance of downlink transmission, training signals are used to estimate the channel. The received signals at the th subcarrier is where is training signals transmitted at the th subcarrier. is precoding matrix, and is combining matrix. denotes addictive noise matrix at the th subcarrier, which follows complex Gaussian distribution with . is channel matrix of the th subcarrier between transmitter and receiver in frequency domain. Consider the scenario that scattering multipath is involved in mmWave channel, where scatterers existed between transmitter and receiver, and paths of signals are propagated from the th scatterer [14]. Let . The mmWave channel with delay spread can be formulated as follows. where is the fading coefficient of the th path from the th scatterer, following the distribution of . denotes the impacts from the pulse shaping filter, matching filter, etc. at time slot. is the array steering vector at transmitter, where is the AoD of the path at the th scatterer. is the array steering vector at receiver, where is the AoA of the th path at the th scatterer. Assume that uniform linear array is exploited at the transmitter and receiver, and can be formulated as where is the wavelength of signal and is the distance between adjacent antennas. By using points DFT on , the channel matrix in frequency domain at the th subcarrier is where
Define matrices , , . Let , , , where stands for a diagonal matrix with elements on its diagonal, and stands for a block diagonal matrix constructed by the submatrices and . Considering the signals from all paths jointly, Eq. (4) can be formulated in compact form as follows.
Assume that the training symbols with equal and unit power are used for channel estimation. Then, the received signal at the th subcarrier can be formulated as
According to (6), the matrices have block structure, and linearly dependent columns are involved in . Based on this observation, high-order PARALIND model can be exploited to analyse the multipath channels of broadband mmWave system. First, we present some basic of this model.
3. Preliminaries
3.1. PARALIND Model
PARALIND is a kind of high-order tensor model [15], which is an extension of the traditional parallel factor (PARAFAC) model. Consider a three-order tensor . In tensor analysis, three-order tensor can be decomposed into three factor matrices. If the uniqueness condition of tensor decomposition is satisfied, three-order tensor can be uniquely determined from the factor matrices, so that the factor matrices are usually used to represent high-order tensor for simplicity. Assume that are the factor matrices of three-order tensor, and linearly dependent columns are included in . A linearly dependent matrix is used to formulate the linearly dependent columns in the matrix , such as where is constructed by the independent columns in . is the linearly dependent matrix which represents the correlation between columns in . For example, assume that . Then, the and are
As we know, three-order tensor can be formulated as a cube. However, it is difficult to analyse a cube data directly in data processing. For simplicity, the cube data is usually converted to its matrix form by slicing it along one dimension [15]. By slicing along the third dimension, the th slicing matrix of based on PARALIND model is where represents a diagonal matrix with the elements of the th row of factor matrix on its diagonal. According to the definition of PARALIND model [15], the slicing form of in the third dimension is where denotes the Khatri-Rao product operator. Compared with the PARAFAC model, PARALIND model introduces a linearly dependent matrix, which can not only show the correlation structure of the vectors in the factor matrix but also represent the relationship between the corresponding data vectors of different factor matrices. More importantly, the high-order tensor with linear correlated columns can be effectively analysed by PARALIND model due to the linearly dependent matrix, which lays a foundation for the application of this model. The uniqueness of PARALIND model is analysed in the following section.
3.2. Uniqueness of PARALIND Model
In this section, we analyse the uniqueness properties of PARALIND model and present some results. First, we give the definition of rank and rank of a matrix.
Definition 1 (see [16]). Consider a matrix . The rank of is the maximal number such that any columns of is linearly independent.
Definition 2 (see [17]). Consider a matrix . Partition into blocks, such as , where the submatrix is formulated as , and . The rank of is the maximal number such that any set of submatrices of yields a set of linearly independent columns.
According to Definition 2, the rank of depends on the partitioned matrices. Specially, if , we have . Based on the definitions of rank and rank, the uniqueness condition of PARALIND model is given in [17].
Theorem 3 (see [17]). Assume that the three-order tensor can be modelled as PARALIND, of which the factor matrices are . The linearly dependent matrix of is , which is formulated as
Partition and according to the structure of , such as , where , , . Assume that the following condition is satisfied.
Then, if another triple matrices can be obtained from the PARALIND decomposition of , the satisfy the following equations where is permutation matrix and is nonsingular diagonal matrix. are block permutation matrices, and are block nonsingular diagonal matrices, with the same partitioned structure as and , respectively.
According to Theorem 3, if Eq. (13) is satisfied, structural uniqueness is involved in the decomposition of PARALIND model. That is, the factor matrix can be uniquely determined from PARALIND decomposition. The factor matrices and decomposed from PARALIND model have the block structure, and their column space can be uniquely determined. Based on the uniqueness results, PARALIND model has been applied in many areas. For example, PARALIND model has been used in the quantitative analysis of the compounds bifenthrin and tetramethrin [18] and has been demonstrated to be more suitable than PARAFAC model in the analysis of these compounds due to the overlap of the spectral peaks of several compounds. PARALIND model has also been applied in the quantitative analysis of flow injection data (FIA) and fluorescence spectrum data in [15, 19]. In this paper, the channel estimation problem of mmWave MIMO system is represented as the decomposition problem of PARALIND model, and a new algorithm is designed for channel estimation.
4. PARALIND-Based Channel Estimation
4.1. Problem Formulation
According to the system model, the received signals at the th subcarrier is
By arranging the elements of into vectors, such as , (15) has the following formulation. where is composed of the diagonal elements in and . In Eq. (16), we use the property of Khatri-Rao product, that is, , where denotes a vector constructed by the diagonal elements of . Let , , and . The received signals at all subcarriers can be formulated as
According to the structure of array steering manifold, can be formulated as , where
Then, (17) can be formulated as
Compared (19) to (11), it is shown that the received signals of mmWave MIMO system can be modelled as PARALIND, where , , and are three factor matrices, and is the linearly dependent matrix. Based on PARALIND model, we propose a new algorithm for channel estimation.
4.2. PARALIND-Based Channel Estimation Algorithm
In this section, a channel estimation algorithm for mmWave MIMO system is proposed based on the PARALIND model. For simplicity, let . (19) can be simplified as
According to (20), the received signals include three dimensions of transmit, receive, and frequency and can be modelled as PARALIND, where and are factor matrices and linearly dependent matrix. is the slicing matrix of the model in transmit dimension. Based on this observation, trilinear alternating least square (TALS) method can be used for PARALIND decomposition. First, the estimation of can be obtained by solving the following problem.
Then, the least square solution of is
Second, by arranging the elements of , the slicing matrix of PARALIND model in frequency dimension can be formulated as where is the slicing noise matrix in this dimension. The estimation of can be obtained by solving the following problem.
Then the least square solution of is
Similarly, the slicing matrix of PARALIND model in receive dimension can be obtained by arranging the elements of , such as where is the slicing noise matrix in this dimension. The estimation of can be obtained by solving the following problem.
Then, the least square solution of is
At last, the vectorized version of can be formulated as where . The estimation of can be obtained by solving the following problem.
Therefore, the least square solution of (30) is
The estimation of can be obtained by arranging . Repeat the above steps until convergence achieved. The PARALIND decomposition algorithm is summarized in Algorithm 1.
| ||||||||||||||||||
According to Theorem 3, matrix can be determined by PARALIND decomposition, so that the estimation of can be obtained from , formulated as
Then, can be obtained from . Similarly, can also be calculated from in least square sense, such as
According to the uniqueness property of PARALIND model, the matrix obtained by PARALIND decomposition is structurally related to the original . Based on this observation, we can uniquely determine the estimation of from . First, let and , where are the th columns of and , respectively. Partition and according to the structure of , such as and , where and , and . According to Theorem 3, and have the same column space, such as where stands for a column space constructed by vectors and . Define nonsingular matrices . (34) can be written as
Construct two submatrices and from , such as where stands for a submatrix composed of the elements of from the th row to the th row. According to (35), we have the following equations.
Let , . Note that is a Vandermonde matrix, such as where . Therefore, we have where . Then, and can be formulated as
Define , according to (40), can be written as
According to the definition of eigenvalue decomposition, is the eigenvalue of . Therefore, can be calculated by the eigenvalue decomposition of , such as where denotes the eigenvalue decomposition. Since , the can be estimated from , and then can be estimated from (38). Note that the channel parameters of multipath from each scatter can be estimated separately by using PARALIND-based algorithm, which can reduce the computation complexity in channel estimation. Then let . The estimation of can be obtained by solving the following problem, given and
The least square solution of (43) is
Denote . The estimation of channel fading coefficients at the th subcarrier can be obtained by . Lastly, the mmWave channel can be estimated by using (6). The proposed PARALIND-based channel estimation algorithm is summarized in Algorithm 2.
| ||||||||||||||||||||||||
5. Discussion
According to (6), the channel estimation problem is equivalent to identify the parameters in matrices . In this section, the identifiability conditions for channel parameter estimation are analysed. Without loss of generality, noise is ignored in the following analysis, so that the received signal can be simplified as where and . According to Algorithm 2, channel parameters can be estimated only if the factor matrices are uniquely determined from PARALIND decomposition. Therefore, the uniqueness of PARALIND decomposition is the key issue in channel parameters estimation. Based on this observation, we present an essential condition of channel estimation based on the uniqueness property of PARALIND model. According to Theorem 3, the uniqueness condition of (45) is
Based on (46), we present the identifiability results for and by discussing the value of and . First, we present the following theorem.
Theorem 4. Assume that is a Vandermonde matrix, such as
Let , where is a matrix with full column rank. Then, the rank of is .
Proof. Consider the following two cases.
(1)Since is a Vandermonde matrix, so that is full column rank. Note that is a matrix with full column rank. According to the property of rank, . Therefore, is also full column rank, and all columns of are linearly independent. According to the definition of rank, .
(2)Assume that is a submatrix constructed by any set of vectors in . Define matrix , we have . Note that is a Vandermonde matrix, so that . According to the property of rank, , which means that any vectors in are linear independent. Meanwhile, define matrix , where is a submatrix constructed by any set of vectors in , we have . Since is also a Vandermonde matrix, so that . Then, we have , which means that there exists column vectors of being linearly dependent. According to the definition of rank, .
Based on the former results, we have .
Since is a Vandermonde matrix, also has Vandermonde structure. Assume that . According to [20], . Since , we have based on Theorem 4.
The rank of depends on the partitioned structure of . Since , partition according to the structure of , such as , where . Then, can also be formulated as partitioned matrix, such as , where . Sort from largest to smallest according to the number of columns in the submatrix, such as , and assume that the AoAs of multipath signals are different. Define two submatrices of and , formulated as and , where . Two cases are considered in determining the rank of . (1)
Note that is a Vandermonde matrix, so that is full column rank, and the columns in are linearly independent. According to Theorem 4, the columns in are also linearly independent. Since consists of submatrices of with the most columns, the number of columns in any matrix consisting of submatrices of is less than . Therefore, the columns of any submatrices of are linearly independent. (2)
In this condition, is a matrix with full row rank. Based on the definition of matrix rank and the structure of , is also full row rank. Therefore, the columns of are linearly dependent.
According to the former analysis, if , the columns of any submatrices in are linearly independent. Note that is discrete, that is, and . Therefore, if , the columns of any submatrices in are linearly independent, and there are submatrices containing linearly dependent columns. According to the definition of rank, the rank of is .
Lastly, we discuss the rank of . Similarly, the rank of depends on its partitioned structure, so that we partition to submatrices according to the structure of , such as , where .. Note that is a random matrix whose columns are drawn independently from an absolutely continuous distribution. Based on the results in [21], are not only full rank but also full rank. Similar to the results of , if , the columns of any submatrices in are linearly independent, and there are submatrices containing linearly dependent columns. Then, the rank of is .
Based on the above results, the identifiability condition of PARALIND-based channel parameters estimation is related to the value of . In mmWave communication scenarios, the number of scatterers is usually less than that of antennas, such that . Then, (46) can be simplified as
According to Eq. (48), a necessary condition of and for channel parameter estimation can be derived. Take for example. Note that . It means that if and only if , the channel parameters can be determined from PARALIND-based algorithm. To satisfy this condition, should be larger than . Similarly, the value of should also be larger than . Therefore, a necessary condition of PARALIND-based channel parameters estimation is . It implies that the identifiability of channel parameters is related to the scatterers with most paths in mmWave channel. Specially, if only single path is generated from each scatterer, the value of and can be reduced to two, which can reduce the complexity and overhead in channel estimation effectively.
6. Numerical Results
In this section, we present some numerical results to illustrate the performance of the proposed channel estimation algorithm. The system is assumed to be operated at 28 GHz carrier frequency with a bandwidth of 100 MHz. Both the transmitter and receiver employ uniform linear arrays. The AoAs/AoDs of each path is randomly selected and uniformly distributed in [0,]. The fading coefficient of each path, the elements of precoding, and combining matrices are selected following the distribution . First, we evaluate the normalized mean square error (NMSE) performance of the proposed algorithm in channel estimation. The NMSE of the estimated channel is formulated as where is the expectation operator. Besides, the traditional CS-based mmWave channel estimation algorithm is also simulated just for comparison [6]. Figure 1 depicts the NMSE performance of algorithms under different SNR. The number of transmit and receive antennas is assumed to be 16. The length of training sequence and the number of combining vectors is 16. In order to evaluate the performance of the proposed algorithm under different numbers of subcarriers, we also present the NMSE performance curves when the number of subcarriers is 16, 32, and 64. It shows that the proposed PARALIND-based channel estimation algorithm performs better than the CS-based algorithm. As we discussed, the proposed algorithm exploits the multidimensional structure of mmWave channels in channel modelling. Therefore, compared to CS-based algorithm, the proposed algorithm is not suffered from the quantization error during the channel estimation. Furthermore, according to the basic of CS theory, the performance of CS-based recovery algorithm is related to the correlation of columns in measurement matrix [22]. However, when scattering multipath is involved in mmWave channel, the correlation between vectors in measurement matrix is high, resulting in the performance degradation of CS-based algorithm. Meanwhile, as the number of subcarriers increase, the performance of proposed PARALIND-based algorithm is improved. Equation (19) shows that the dimension of factor matrix in frequency is represented by the number of subcarriers. According to the uniqueness condition of PARALIND decomposition, the increase of the dimensions of factor matrix may improve the accuracy of trilinear decomposition, so that the performance of channel estimation is improved.
Figure 2 shows the NMSE performance of the proposed algorithm when different lengths of pilot sequence are used at the transmitter, where is set to be 16, 32, 64, and is 16. Meanwhile, Figure 3 evaluates the performance of proposed algorithm when different numbers of combining vectors are used at the receiver, and is set to be 16, 32, 64, and is 16. In these two simulations, both 16 antennas are used at transmitter and receiver. The number of subcarriers is 32. It shows that the performance of the proposed algorithm become better when the number of pilots and combining vectors increase. The reason for these results is that when more pilots and combining vectors are used, more samples can be obtained for channel estimation at the receiver, which may improve the performance of channel estimation. However, although increasing the length of pilot symbols and the number of combining vectors can improve the performance of channel estimation, it may also lead to the increase of pilot overhead and processing complexity at the transmitter and receiver, respectively. Therefore, a tradeoff between pilot overhead and system complexity is involved in channel estimation. For example, the pilot overhead can be reduced by increasing the complexity of receiver, and vice versa.
Figures 4 and 5 show the performance curves of the proposed algorithm when different numbers of antennas are used at the transmitter and receiver, where the number of transmit and receive antennas is 16, 32, and 64, respectively. We find that the performance of the proposed algorithm become worse when more antennas are used in the system. According to Algorithm 2, the parameters of paths are estimated to reconstruct array manifold for channel estimation. Due to the error in parameter estimation, the accuracy of reconstructed array manifold matrix may be decreased when the number of antennas increases, which leads to decreasing the performance of channel estimation.
Finally, we evaluate downlink transmission performance of the system using the channels estimated by the proposed algorithm. In this simulation, achievable spectral efficiency (ASE) is used as the performance metric [11], formulated as where and are precoding and combining matrices. According to [11], and can be designed by the singular value decomposition (SVD) of . Let , where and are matrices constructed by left singular value vectors and right singular value vectors of . is the singular value matrix of . The precoding matrix can be constructed by the first vectors of , and the combining matrix can be constructed by the first vectors of . In this simulation, the number of transmit and receive antennas are 16. The number of combining vectors and subcarriers are 16 and 32. The length of pilot sequence is chosen as 16, 32, and 64. We also present the ASE performance by using the optimal precoder, which is calculated from the perfect downlink channels. Figure 6 shows the ASE performance curves under different SNR. We find that the ASE performance of system using the estimated channel is close to that using the perfect channel, especially in the high SNR scenarios. It also shows that the performance of downlink transmission can be improved by using more pilot symbols.
7. Conclusions
In this paper, the estimation problem of mmWave channel is studied. According to the high-order dimension structure of broadband mmWave channel, the multipath channels are modelled as PARALIND, and a channel estimation algorithm is proposed based on PARALIND model decomposition. The proposed algorithm can separate the scattering multipath generated by different scatterers and estimate the channel parameters of each path. Meanwhile, we analyse the identifiability conditions of channel parameters based on the uniqueness property of PARALIND model and present a necessary condition for channel parameters estimation. Simulation results guarantee the effectiveness of the proposed algorithm.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported in part by Aeronautical Science Foundation of China (201901052001), Natural Science Foundation of China (61871237, 92067101), the Jiangsu Key Research and Development Program (BE2020084-3), China Postdoctoral Science Foundation (2012M521106), and Natural Science Foundation of China (61101104).