CN104678350A - TLS-ESPRTT algorithm-based 2D DOA estimation in large scale MIMO system - Google Patents

Abstract

The invention discloses a method for estimating a 2D DOA (2-Dimensional Direction of Arrival) by utilizing a TLE ESPRIT (Total Least Square-Estimating Signal Parameter via Rotational in Variance Techniques) algorithm in a large scale MIMO system. The method comprises the following steps: dividing an area array into two partially-overlapped sub-arrays according to a received signal model of the large scale MIMO system, so that the two partially-overlapped sub-arrays meet the requirement for the array shift invariance of the ESPRIT algorithm; calculating a signal subspace of each sub-array according to a 1D ESPRIT (1-Dimensional ESPRIT) algorithm, and estimating elevation information according to a TLS rule; estimating an azimuth by utilizing a direction matrix of one of the sub-arrays, and automatically matching the azimuth with the elevation. According to the method, the problem about 2-dimensional DOA estimation in the large scale MIMO system is effectively solved by adopting the thought of the TLS-ESPRIT algorithm; the azimuth can be matched with the elevation automatically; the computing complexity is relatively low.

Description

Estimate based on the 2D DOA of TLS-ESPRTT algorithm in extensive mimo system

Technical field

The present invention relates generally to future mobile field.

Background technology

Following mobile data services amount will exponentially increase, and traditional MIMO technology can not meet such demand.In the end of the year 2010, Bell Laboratory scientist Thomas L.Marzetta proposes extensive MIMO (Large Scale MIMO), is also the concept of Massive MIMO.Extensive MIMO as a kind of gordian technique of future broadband wireless communication systems, at the antenna that base station end configured number is huge, with degree of depth digging utilization Spatial Dimension unlimited resources, elevator system spectrum efficiency and power.Extensive mimo system adopts face battle array, cylindrical array or special-shaped battle array in base station end usually.In extensive mimo system, because the multiple antennas of base station configures, cause base station to obtain channel condition information (Channel State Information, CSI) very difficult, the precoding channel feedback that places one's entire reliance upon in traditional base station is no longer applicable.2D DOA is estimated as channel precoding and provides accurate spatial information (si), and the cost volume that realizes of optimum downlink precoding is significantly reduced.

In numerous extensive MIMO base station antenna configuration, URA is relatively simple, can make full use of the spatial information that horizontal dimensions and vertical dimensions provide simultaneously, therefore, be more suitable for applying in extensive mimo system.Under URA configuration, Received signal strength comprises vertical dimensions and horizontal dimensions information simultaneously, and 2D Mutual coupling is exactly to realize the position angle of incoming wave signal level dimension and the elevation estimate of vertical dimension.

Mutual coupling is as one of the important research content of Array Signal Processing, and at sonar, radar, the fields such as wireless communication system have important using value.At present, existing 1D DOA estimation theory has developed to obtain relative maturity, 1D DOA estimation method is as multiple signal classification (Multiple Signal Classification, MUSIC) algorithm, invariable rotary Subspace algorithm (Estimating Signal Parameter via Rotational In variance Techniques, ESPRIT) algorithm application is relatively more extensive, but is difficult to obtain practical application in two-dimensional estimation.Compare with 1D DOA, 2D DOA can state signal space position better, therefore studies 2D DOA and has more practical significance.When utilizing classical MUSIC algorithm to carry out two-dimensional estimation, need to carry out spectrum peak search, operand is large, and consumes a large amount of time.Two dimensional ESPRIT algorithm solves this problem preferably, but needing to match to final result just can obtain effective angle estimation.On the one hand, due to antenna number increase and 2D DOA estimate in the increase of parameter dimension, cause the increase of estimation procedure complexity.If still adopt traditional subspace class algorithm such as 2D MUSIC, 2D ESPRIT algorithm to estimate angle, the problems such as the large and angle pairing of algorithm complexity, calculated amount will be faced, make it can not be used for actual direction-finding system; On the other hand, due to the interference of noise, itself there is error in sampled signal.Sampled signal is directly carried out computing as signal itself by existing most of DOA algorithm for estimating, does not meet actual direction finding environment.

Summary of the invention

The present invention is directed to the deficiency that current 2D DOA estimation method exists, propose the 2D DOA estimation method based on TLS-ESPRIT algorithm in extensive mimo system, the thought of TLS-ESPRIT algorithm is adopted effectively to solve arrival direction estimation problem in extensive mimo system, not only automatically can realize position angle and elevation angle pairing, and computation complexity is relatively low.

The present invention adopts URA, to receive different directions incoming wave in base station end.According to the Received signal strength model of extensive mimo system, first face battle array is divided into two partly overlapping submatrixs, makes it meet the array motion immovability of ESPRIT algorithm requirement; Then the signal subspace of subarray is calculated according to 1D ESPRIT algorithm; Recycling TLS criterion, under noise disturbance condition, estimates accurate elevation information; Direction matrix finally by a wherein submatrix completes position angle and estimates, and automatically and the elevation angle match.

The present invention includes following steps:

Step one: base station adopts M × N to tie up URA, and consideration face battle array is positioned at XOZ plane, and supposing has K information source to incide this array, and array received signal model is as follows:

X＝AS+N

Wherein X is that array received signal phasor is tieed up in MN × 1, and S is K × 1 dimension space signal phasor, and N is MN × 1 white Gaussian noise vector, and A is that MN × K ties up direction matrix, and concrete form is as follows:

$A=\left[\begin{array}{cccc}a\left(u_1\right)\⊗a\left(v_1\right)& a\left(u_2\right)\⊗a\left(v_2\right)& ...& a\left(u_K\right)\⊗a\left(v_K\right)\end{array}\right]$

Wherein represent gram labor Roc inner product, a (u _i) be the direction vector of X-direction, a (v _i) be the direction vector of Z-direction.

Step 2: consider that base station end aerial array is positioned at XOZ plane, prolong Z-direction and M × N dimension URA is divided into two submatrixs with motion immovability, because the direction matrix of submatrix 1 and submatrix 2 differs a twiddle factor Φ, and this twiddle factor only comprises elevation information, specific as follows:

Φ＝diag[exp(-2πd cosθ ₁/λ) exp(-2πd cosθ ₂/λ) ... exp(-2πd cosθ _K/λ)]

Wherein d is array element distance, and λ is velocity of wave, θ _i, i=1,2,3...K represent the elevation angle of i-th user.

Step 3, signal subspace E according to 1D ESPRIT algorithm computing array _s, definition selection matrix J ₁, J ₂, therefrom select the signal subspace that two submatrixs are corresponding detailed process is as follows:

Definition

Then have:

$\begin{array}{cc}{J}_1=P_1\⊗I_{8\×8}& J_2=I_{8\×8}{\⊗P}_2\end{array}$

Then the signal subspace of submatrix 1 and submatrix 2 correspondence is:

There is noise disturbance in step 4, consideration, according to TLS criterion, estimates elevation information, be specifically expressed as follows:

According to TLS criterion, definition:

$E_{12}={\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right]}^H\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right]$

Right carry out Eigenvalues Decomposition:

$E_{12}^H{E}_{12}=\mathrm{E\Λ}{E}^H$

And E is resolved into the submatrix of K × K, that is:

$E=\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]$

Calculate and Eigenvalues Decomposition is carried out to Ψ:

Ψ＝HΩH ^-1＝TΦT ^-1

Eigenvalue λ in Ω _i, the diagonal entry in i=1,2...K homography Φ, H=T, calculates elevation estimate value:

θ _i＝arccos{angle(λ _i)λ/2πd}

Wherein d is array element distance, and λ is velocity of wave, θ _i, i=1,2,3...K represent the elevation angle of i-th user, () ^hrepresent the conjugate transpose asking for matrix, () ^-1represent and ask for inverse of a matrix.

The direction matrix of step 5, submatrix 1 can by signal subspace with the product representation of non-singular matrix H, utilize the column vector of its direction matrix to estimate position angle, concrete formula is as follows:

${\hat{A}}_1={\hat{E}}_{1}H$

in the i-th, i=1,2...K column vector meet:

Wherein represent that i-th arranges front 1:M (N-2) individual element, represent rear N-1:M (N-1) the individual element of the i-th row.

Order position angle φ can be estimated according to above formula:

φ _i＝arccos(angle(b _i)λ/(2πd sinθ _i))

Contemplated by the invention the noise that sampled signal exists, utilize TLS criterion, estimate the accurate elevation angle; Complete position angle by direction matrix column vector invariable rotary relation to estimate.The present invention only utilizes 1D DOA algorithm for estimating to realize 2D angle estimation, significantly reduces computation complexity, and can realize angle automatic matching.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of 2D DOA method in the extensive mimo system that proposes of the present invention;

Fig. 2 is extensive MIMO 3D model schematic;

Fig. 3 is that the URA that base station end adopts divides schematic diagram;

Fig. 4 is the arrival direction estimation method flow diagram based on TLS-ESPRIT algorithm in the extensive mimo system that proposes of the present invention.

Embodiment

In order to make the object, technical solutions and advantages of the present invention clearly, below in conjunction with accompanying drawing, the present invention is described in further detail:

Figure 1 shows that the FB(flow block) of 2D DOA method in the extensive mimo system that the present invention proposes, base station end configuration M × N ties up URA, wherein M represents URA horizontal dimension antenna number, N represents URA vertical dimension antenna number, ensure that each transmit antennas can both simultaneously processing horizontal dimension and the information in vertical dimensions, it is that K is individual that user side sends signal.According to the Received signal strength model of extensive mimo system, first face battle array is divided into two partly overlapping submatrixs, is calculated the signal subspace of subarray by 1D ESPRIT algorithm; Then utilize TLS criterion, estimate elevation information; Finally utilize the direction matrix of a wherein submatrix to estimate to complete position angle, and realize automatic matching.

Fig. 2 is extensive MIMO 3D model schematic, and base station height is h, and user is highly h _u, base station end configuration is positioned at the M capable N row URA of X-Z plane, and array element distance is d, and sending signal number is K.For each send signal have a position angle φ and elevation angle theta corresponding with it, system model is as follows:

X＝AS+N

Wherein X is that array received signal phasor is tieed up in MN × 1, and S is K × 1 dimension space signal phasor, and N is MN × 1 white Gaussian noise vector, and A is that MN × K ties up direction matrix, under concrete form:

$A=\left[\begin{array}{cccc}a\left(u_1\right)\⊗a\left(v_1\right)& a\left(u_2\right)\⊗a\left(v_2\right)& ...& a\left(u_K\right)\⊗a\left(v_K\right)\end{array}\right]$

a(u _i)＝[1 exp(-j2πdsinθ _icosφ _i/λ) ... exp(-j(M-1)2πdsinθ _icosφ _i/λ)] ^T

a(v _i)＝[1 exp(-j2πdcosθ _i/λ) ... exp(-j(N-1)2πdcosθ _i/λ)] ^T

Wherein represent gram labor Roc inner product, a (u _i) be the direction vector of X-direction, a (v _i) be the direction vector of Z-direction, θ _iand φ _irepresent respectively _ithe elevation angle of individual user and position angle, i=1,2,3...K, d are array element distance, and λ is velocity of wave, () ^trepresent and ask for transpose of a matrix.

Fig. 3 is that the URA that base station end adopts divides schematic diagram, consider that base station end aerial array is positioned at XOZ plane, prolong Z-direction and planar array is divided into two submatrixs with motion immovability, two submatrix array numbers are identical, the direction matrix of submatrix 1 and submatrix 2 differs a twiddle factor Φ, this twiddle factor only comprises elevation information θ, and concrete form is as follows:

A ₁Φ＝A ₂

Φ＝diag[exp(-2πd cosθ ₁/λ) exp(-2πd cosθ ₂/λ) ... exp(-2πd cosθ _K/λ)]

Wherein A ₁, Α ₂be respectively the direction matrix of submatrix 1 and submatrix 2, d is array element distance, and λ is velocity of wave, θ _i, i=1,2,3...K represent the elevation angle of i-th user.

Fig. 4 is the arrival direction estimation method flow diagram based on TLS-ESPRIT algorithm in extensive mimo system, and specific implementation step is as follows:

Step 41, as shown in Figure 3, the URA being positioned at X-Z face is divided into two partly overlapping submatrixs along Z axis, the direction matrix of submatrix 1 and submatrix 2 differs a twiddle factor Φ, and this twiddle factor only comprises elevation information θ.

Step 42, according to 1D ESPRIT algorithm, the covariance matrix of first computing array output:

$R_{\mathrm{XX}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{XX}}^H$

Wherein N is the fast umber of beats of signal sampling, then to R _xXcarry out Eigenvalues Decomposition, obtain signal subspace E _s:

$R_{\mathrm{XX}}=E_x{\mathrm{\Λ}}_s{E}_s^H+{\mathrm{\σ}}^{2}I$

Wherein σ ²for noise variance, () ^hrepresent the conjugate transpose asking for matrix.

Step 43, definition selection matrix J ₁, J ₂signal subspace E _sin select signal subspace corresponding to two submatrixs specifically be expressed as follows:

Definition:

Then have:

$\begin{array}{cc}{J}_1=P_1\⊗I_{8\×8}& J_2=I_{8\×8}{\⊗P}_2\end{array}$

Wherein represent gram labor Roc inner product.The signal subspace that so submatrix is corresponding is:

$\begin{array}{cc}{\hat{E}}_1=J_1{E}_s& {\hat{E}}_2=J_2{E}_s\end{array}$

Step 44, estimate the elevation angle according to TLS criterion, detailed process is as follows:

Under A, ideal conditions, the signal subspace E of submatrix 1 and submatrix 2 ₁, E ₂with corresponding direction matrix A ₁, Α ₂open into similar subspace, so the non-singular matrix T of existence anduniquess, following formula set up:

A ₁＝E ₁T A ₂＝E ₂T

According to conclusion A in Fig. 3 ₁Φ=A ₂, can E be obtained ₁t Φ T ^-1=E ₂; Definition Ψ=T Φ T ^-1, then E is had ₁Ψ=E ₂.Wherein, () ^-1represent and ask for inverse of a matrix.

There is error in B, actual samples signal, defences the signal subspace that difference solves jointly all there is noise disturbance:

${\hat{E}}_1\mathrm{\Ψ}\≈{\hat{E}}_2$

Utilize TLS criterion to solve Ψ more more accurate than traditional E SPRTI algorithm, be specifically expressed as follows:

Definition:

$E_{12}={\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right]}^H\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right]$

Right carry out Eigenvalues Decomposition:

$E_{12}^H{E}_{12}=\mathrm{E\Λ}{E}^H$

And E is resolved into the submatrix of K × K, namely

$E=\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]$

Calculate and Eigenvalues Decomposition is carried out to Ψ:

Ψ＝HΩH ^-1＝TΦT ^-1

Eigenvalue λ in Ω _i, the diagonal entry in i=1,2...K homography Φ, H=T, calculates elevation estimate value:

θ _i＝arccos{angle(λ _i)λ/2πd}

Wherein d is array element distance, and λ is velocity of wave, () ^hrepresent the conjugate transpose asking for matrix, () ^-1represent and ask for inverse of a matrix.

Step 45, utilizes the direction matrix of submatrix 1 to estimate position angle, and the direction matrix of submatrix 1 correspondence by the product representation of signal subspace and above-mentioned non-singular matrix H, specifically can be expressed as follows:

The direction matrix A of submatrix 1 ₁:

A ₁＝[α ₁α ₂… α _i… α _K] _M(N-1)×K

Wherein α _iconcrete form be: ${\mathrm{\α}}_i={\left[\begin{array}{c}1\\ \exp (-j2\mathrm{\π}d\cos {\mathrm{\θ}}_1/\mathrm{\λ})\\ .\\ .\\ .\\ \exp (-j2\mathrm{\π}(N-2)d\cos {\mathrm{\θ}}_1/\mathrm{\λ})\\ \exp (-j2\mathrm{\π}d\sin {\mathrm{\θ}}_1\cos {\mathrm{\φ}}_1/\mathrm{\λ})\\ \exp (-j2\mathrm{\πd}(\sin {\mathrm{\θ}}_1{\cos \mathrm{\φ}}_1+\cos {\mathrm{\θ}}_1)/\mathrm{\λ})\\ .\\ .\\ .\\ \exp (-j2\mathrm{\πd}(\sin {\mathrm{\θ}}_1{\cos \mathrm{\φ}}_1+(N-2)c{\mathrm{os\θ}}_1)/\mathrm{\λ})\\ .\\ .\\ .\\ .\\ .\\ .\\ \exp (-j2\mathrm{\πd}(M-1)\sin {\mathrm{\θ}}_1{\cos \mathrm{\φ}}_1/\mathrm{\λ})\\ \exp (-j2\mathrm{\πd}((M-1)\sin {\mathrm{\θ}}_1\cos {\mathrm{\φ}}_1+{\cos \mathrm{\θ}}_1)/\mathrm{\λ})\\ .\\ .\\ .\\ \exp (-j2\mathrm{\πd}((M-1)\sin {\mathrm{\θ}}_1{\cos \mathrm{\φ}}_1+(N-2){\cos \mathrm{\θ}}_1)/\mathrm{\λ})\end{array}\right]}_{M(N-1)\×1}$

By α _iconcrete form can find out, A ₁in i-th arrange front 1:M (N-2) individual element and rear N-1:M (N-1) individual element meets invariable rotary relation:

Utilize the signal subspace of submatrix 1 direction matrix corresponding to submatrix is gone out with the product estimation of non-singular matrix H:

${\hat{A}}_1={\hat{E}}_{1}H$

Step 46, utilizes the invariable rotary relation of middle column vector is by twiddle factor solve, estimate position angle φ, concrete formula is as follows:

φ _i＝arccos(angle(b _i)λ/(2πd sinθ _i))

Wherein d is array element distance, and λ is velocity of wave, φ _i, i=1,2,3...K represent the position angle of i-th user, () ⁺represent the pseudoinverse asking for matrix.

Be more than that the preferred embodiments of the present invention are described in detail, should be appreciated that preferred embodiment only in order to the present invention is described, instead of in order to limit the scope of the invention.

Claims (5)

1. in extensive mimo system based on the 2D DOA estimation method of TLS-ESPRIT algorithm, it is characterized in that, comprise step: base station end adopts M × N to tie up uniform rectangular face battle array (URA), according to the Received signal strength model of extensive mimo system, first face battle array is divided into two partly overlapping submatrixs, makes it meet the array motion immovability of ESPRIT algorithm requirement; And then the signal subspace of subarray is calculated according to 1D ESPRIT algorithm; Then utilize TLS criterion, realize elevation estimate; The direction matrix of a wherein submatrix is finally utilized to estimate corresponding position angle.

2. 2D DOA estimation method according to claim 1, is characterized in that, described two partly overlapping submatrixs are submatrix 1 and submatrix 2, submatrix 1 direction matrix A ₁with submatrix 2 direction matrix A ₂meet relation: A ₁Φ=A ₂, wherein include elevation information in matrix Φ.

3. 2D DOA estimation method according to claim 1, it is characterized in that, the described signal subspace according to 1D ESPRIT algorithm calculating subarray is: according to 1D ESPRIT algorithm, the covariance matrix of pair array sampled data carries out Eigenvalues Decomposition, obtains signal subspace E _s, from E _sin select the signal subspace corresponding with submatrix 1 and submatrix 2 with that is: wherein J ₁, J ₂represent corresponding selection matrix respectively.

4. 2D DOA estimation method according to claim 3, is characterized in that, describedly utilizes TLS criterion, and realizing elevation estimate is: according to TLS criterion, estimates accurate elevation information, definition $E_{12}={\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right]}^H\left[\begin{array}{cc}{\hat{E}}_1& {\hat{E}}_2\end{array}\right],$ Right carry out Eigenvalues Decomposition and obtain eigenmatrix E, and E is resolved into the submatrix of K × K, that is:

$E=\left[\begin{array}{cc}{E}_{11}& E_{12}\\ E_{21}& E_{22}\end{array}\right]$

Calculate and Eigenvalues Decomposition is carried out to Ψ:

Ψ＝HΩH ^-1＝TΦT ^-1

Eigenvalue λ in Ω _i, the diagonal entry in i=1,2...K homography Φ, H=T, calculates elevation estimate value:

θ _i＝arccos{angle(λ _i)λ/2πd}

Wherein d is array element distance, and λ is velocity of wave, θ _irepresent the elevation angle of i-th signal, () ^hrepresent the conjugate transpose asking for matrix, () ^-1represent and ask for inverse of a matrix.

5. 2DDOA method of estimation according to claim 1, is characterized in that, the direction matrix of a described utilization wherein submatrix estimates that corresponding position angle concrete grammar is: the direction matrix of submatrix 1 can by its signal subspace with the product representation of non-singular matrix H: the invariable rotary relation of utilization orientation rectangular array vector can estimate corresponding position angle φ:

φ _i＝arccos(angle(b _i)λ/(2πdsinθ _i))

Wherein $b_i=\left[{{\widehat{\mathrm{\α}}}_i]}_{1:M(N-2)}\right[{\widehat{\mathrm{\α}}}_i{]}_{N-1:M(N-1)},[{{\widehat{\mathrm{\α}}}_i]}_{1:M(N-2)}$ Represent i-th arranges front 1:M (N-2) individual element, represent rear N-1:M (N-1) the individual element of the i-th row, d is array element distance, and λ is velocity of wave, θ _irepresent the elevation angle of i-th signal, φ _irepresent the _ithe position angle of individual signal, () ⁺represent the pseudoinverse asking for matrix.