CN104375133B - Estimation method for space two-dimensional DOA - Google Patents

Abstract

The invention discloses an estimation method for a space two-dimensional DOA and belongs to the field of the radar technology. A cost function is constructed through an obtained guiding vector, reverse solution is conducted on direction and pitch two-dimensional angles through the maximum zero point, a Vandermonde structure that all array elements corresponding to the two-dimensional angles receive data is considered in the construction of the guiding vector, the defect that the angle is determined through a single feature value in an ESPRIT method is overcome, and the estimation accuracy is improved.

Description

A kind of evaluation method of space two-dimensional DOA

Technical field

The present invention relates to Radar Technology field, particularly to a kind of evaluation method of space two-dimensional DOA.

Background technology

In field of radar, determine that DOA (Direction of Arrival, direction of arrival) is always the important class of research Topic.

In existing technology, conventional has maximum likelihood and MUSIC (Multiple Signal Classification, Multiple Signal Classification) and ESPRIT (Estimating Signal Parameters via Rotational Invariance Techniques, by ESPRIT estimating signal parameter) method, wherein ESPRIT By calculating the method for closed solutions it is possible to obtain azimuth and two important parameters of the angle of pitch of information source, thus completing to DOA Estimation, be not required to as maximum likelihood and MUSIC method will to significantly reduce related data through scanning for spectral peak Amount of calculation and amount of storage.

During realizing the present invention, inventor finds that prior art at least has problems with：

As wherein more outstanding method, ESPRIT method needs fractal dimension calculation and parameter pairing in computing, for example It is azimuth and angle of pitch timing really carrying out two-dimentional computing, need separately to be processed above-mentioned two parameter, now Mistake in computation occurs, leads to not DOA is accurately estimated.

Content of the invention

In order to solve problem of the prior art, the invention provides a kind of evaluation method of space two-dimensional DOA, as Fig. 1 institute Show, described inclusion：

Step one, disposes receiving array, includes 2M omnidirectional's electromagnetic transducer in described receiving array, and array contains two Subarray, each submatrix contain M sensor, spacing d between two wherein adjacent sensors less than wavelength X two/ One；

Step 2, is received by described receiving array and includes the information source of P narrowband target the first of the very first time Signal value x₁(t)=A₁s(t)+n₁T (), described receiving array is pressed the first coordinate η=(η_x,η_y) translated, after being translated Receiving array, by described translation after receiving array receive secondary signal value x in the second time for the described information source₂(t) =A₂s(t)+n₂(t)；

Step 3, described first signal value and described secondary signal value are entered ranks storehouse, form high dimensional data

And obtain the piecemeal correlation matrix of described x (t),

After Eigenvalues Decomposition, averaged power spectrum is carried out to minimal eigenvalue, determine noiseApproximationAnd to described R_xCarry out denoising computing, obtain the piecemeal correlation matrix after denoising,

Step 4, obtains random initial mask U (0)=[u₁(0)u₂(0)...u_P(0) l], is made to be cycle-index, construction the One cost function

WhereinC_i1Associate matrix each other, U^H(l-1) with U (l-1) associate matrix each other, to described first Cost function carries out feature decomposition, obtains fisrt feature breakdown

WhereinFor corresponding eigenvalue matrix,For corresponding eigenvectors matrix, if front P special greatly The corresponding characteristic vector of value indicative constitutes V (l-1), constructs the second cost function

WhereinC_i1Associate matrix each other, V^H(l-1) with V (l-1) associate matrix each other, to described second Cost function carries out feature decomposition, obtains second feature breakdown

By the first cost function described in looping construct and described second cost function, and carry out feature decomposition, until full Foot

||U(L)U^H(L)-U(L-1)U^H(L-1)||_F≤ε

When stop circulation, wherein said ε is default threshold value, according to the l stopping during circulation, determines matrix U (L), fixed Adopted signal subspace G=U (L)；

Step 5, constructs piecemeal dimensionality reduction matrix J, and described signal subspace G is on diagonal in described piecemeal dimensionality reduction matrix J Element, by described piecemeal dimensionality reduction matrix J to the piecemeal correlation matrix C after described denoising_xCarry out dimension-reduction treatment, obtain

Described

Step 6, to describedCarry out Eigenvalues Decomposition, obtain

Wherein U_SAnd U_NCorresponding signal space and spatial noise, U_S1And U_S2It is same dimension matrix, to U_S1Generalized inverse matrix With U_S2ProductCarry out feature decomposition, obtain

Construction steering vectorAnd then obtain the estimation of described steering vectorWhereinFor G^HGeneralized inverse matrix；

Step 7, chooses described steering vectorIn m row p column element be a_mp, constructed fuction ρ₁(μ)、ρ₂(μ) as follows：

Wherein,If vectorial z_m=(x_m,y_m),x_mAnd y_mIt is in the 1st subarray The transverse and longitudinal coordinate of m sensor.If vectorial WithIt is the 2nd The transverse and longitudinal coordinate of m-th sensor in subarray.If vectorial

μ=(cos α_psinβ_p,sinα_psinβ_p), if ＜ is z_m, μ ＞ represents vectorial z_mWith the inner product of μ, ifInt represents floor operation, to described function ρ₁(μ)、ρ₂(μ) modulus value is carried out square, has after arrangement

Wherein vectorial z_n=(x_n,y_n),x_nAnd y_nIt is the horizontal stroke of n-th sensor in the 1st subarray Vertical coordinate.If vectorial WithIt is n-th in the 2nd subarray The transverse and longitudinal coordinate of sensor.

Described | ρ₁(μ)|²、|ρ₂(μ)|²Middle and make derivative value be zero respectively to μ derivation, obtain

Choose zero point δ of maximum respectively₁And δ₂, solve from zero point is counter

＜ (z_m-z_n), μ ＞=angle (δ₁),

Substitute into z_m-z_n, μ, have equation below group

IfCan solve from above formula formulaCan get azimuth and the angle of pitch Estimation is respectively

The beneficial effect brought of technical scheme that the present invention provides is：

Cost function is constructed by required steering vector, solves orientation and pitching two dimension angular, guiding by maximum zero point is counter The construction of vector considers the vandermonde structure of the corresponding all array element receiving datas of two dimension angular, only overcomes ESPRIT method The drawbacks of determine angle by single feature value, improves the accuracy of estimation.

Brief description

In order to be illustrated more clearly that technical scheme, the accompanying drawing of required use in embodiment being described below Be briefly described it should be apparent that, drawings in the following description are only some embodiments of the present invention, general for this area For logical technical staff, on the premise of not paying creative work, other accompanying drawings can also be obtained according to these accompanying drawings.

Fig. 1 is the invention provides a kind of flow chart of the evaluation method of space two-dimensional DOA；

Fig. 2 is that the mean square error of the angle estimation that the present invention provides is illustrated with signal to noise ratio change and corresponding Cramér-Rao lower bound Figure；

Fig. 3 is the mean square error of the angle estimation that the present invention provides with hits change and corresponding Cramér-Rao lower bound situation Schematic diagram；

Fig. 4 is the planisphere of the two dimension angular of 3 targets that the evaluation method that the present invention provides is estimated；

Fig. 5 is the planisphere of the two dimension angular of 3 targets that the ESPRIT algorithm that the present invention provides is estimated.

Specific embodiment

Structure and advantage for making the present invention are clearer, below in conjunction with accompanying drawing, the structure of the present invention are made further Description.

Embodiment one

The present embodiment provides a kind of evaluation method of space two-dimensional DOA, and described evaluation method includes：

Step one, disposes receiving array, includes 2M omnidirectional's electromagnetic transducer in described receiving array, and array contains two Subarray, each submatrix contain M sensor, spacing d between two wherein adjacent sensors less than wavelength X two/ One.

In force, comprise 2 subarrays in this receiving array, in each subarray, comprise m omnidirectional's electromagnetic transducer, And subarray can be randomly topologically structured, the value of m is natural number here.

Step 2, receives the first signal including the information source of P narrowband target in the very first time by receiving array Value x₁(t)=A₁s(t)+n₁T (), receiving array is pressed the first coordinate η=(η_x,η_y) translated, the reception battle array after being translated Row, the receiving array after translation receives secondary signal value x in the second time for the information source₂(t)=A₂s(t)+n₂(t).

In force, if two submatrix spacing are η=(η_x,η_y), if space has an information source, in this information source, it is disposed with P Narrowband target, the P narrowband target s_pIt is (α with respect to the two-dimentional arrival direction of receiving array_p,β_p), α_pAnd β_pThe side of representative respectively Parallactic angle and the angle of pitch, if m-th sensor (m=1,2 ..., 2M) coordinate in XOY face is (x in receiving array_m,y_m), if The noise of m-th sensor reception is n_m, then under the t time time sampling, total sound of the p narrowband target that sensor m receives Should be

Wherein a_mp=exp (j2 π Δ_mp/ λ), λ is wavelength, Δ_mpFor

Δ_mp=x_mcosα_psinβ_p+y_msinα_psinβ_p

The t time time sampling receipt signal of 1st submatrix be

x₁(t)=A₁s(t)+n₁(t),

The t time time sampling receipt signal of 2nd submatrix be

x₂(t)=A₂s(t)+n₂(t),

Wherein A₂=A₁Λ, Λ=diag (μ₁,μ₂,…,μ_P), diag represents diagonal matrix,

Step 3, the first signal value and secondary signal value are entered ranks storehouse, form high dimensional data

And obtain the piecemeal correlation matrix of x (t),

After Eigenvalues Decomposition, averaged power spectrum is carried out to minimal eigenvalue, determine noiseApproximationAnd to R_xCarry out denoising computing, obtain the piecemeal correlation matrix after denoising,

In force, above-mentioned first signal value x₁(t) and secondary signal value x₂T () is a series of discrete values, by One signal value and secondary signal value enter ranks storehouse will above-mentioned a series of discrete values as the element value in column vector, from And the vector x (t) of the only string getting, the element in this vector x (t) is x from the first row to last column₁(t)、x₂ Discrete values in (t).According to x (t), through computing R_x=E [x (t) x^H(t)], obtain correlation matrix R_x, it is to x (t) here With associate matrix x^HT () asks desired after carrying out product, in formula (1)Be respectively with to be estimated , due to there is A in azimuth and the corresponding noise figure of the angle of pitch₂=A₁The relation of Λ, A therefore therein₂All use A₁Correlation form Replace, I represents unit matrix.To R_xAfter carrying out Eigenvalues Decomposition, averaged power spectrum is carried out according to minimal eigenvalue, can get noiseApproximation, from R_xMiddle removal noise approximation, that obtain is exactly the correlation matrix C after denoising_x.

Step 4, obtains random initial mask U (0)=[u₁(0) u₂(0) ... u_P(0) l], is made to be cycle-index, structure Make the first cost function

WhereinC_i1Associate matrix each other, U^H(l-1) with U (l-1) associate matrix each other, to the first cost Function carries out feature decomposition, obtains fisrt feature breakdown

WhereinFor corresponding eigenvalue matrix,For corresponding eigenvectors matrix, if front P big The corresponding characteristic vector of eigenvalue constitutes V (l-1), constructs the second cost function

WhereinC_i1Associate matrix each other, V^H(l-1) with V (l-1) associate matrix each other, to the second cost Function carries out feature decomposition, obtains second feature breakdown

By looping construct first cost function and the second cost function, and carry out feature decomposition, until meeting

||U(L)U^H(L)-U(L-1)U^H(L-1)||_F≤ε

When stop circulation, wherein ε is default threshold value, according to the l stopping during circulation, determines matrix U (L), definition letter Work song space G=U (L).

In force, two cost functions are constructed to initial matrix U (0)WithRightObtain after carrying out Eigenvalues Decomposition To withRelated multinomial, rightObtain after carrying out feature decomposition withRelated multinomial, now easily sends out Existing formula (2) to formula (3), be from U (l-1) multinomial toPolynomial process, formula (4) to (5) is from V (l- 1) multinomial arrivesPolynomial process, therefore, formula (2) to formula (3), formula (4) to formula (5) are once (l-1) To the circulation of (l), often complete one cycle, the content in () adds 1, until when the numerical value satisfaction of (l) is not more than threshold value ε, root The numerical value of l when accordingly, determines matrix U (L), and then determines signal subspace G=U (L).

Step 5, constructs piecemeal dimensionality reduction matrix J, and signal subspace G is the element on diagonal in piecemeal dimensionality reduction matrix J, By piecemeal dimensionality reduction matrix J to the piecemeal correlation matrix C after denoising_xCarry out dimension-reduction treatment, obtain

In force, in order to signal subspace G is carried out with dimension-reduction treatment, special construction dimensionality reduction matrix J, this matrix is to angular moment Battle array, and each element is signal subspace G, according to formula (6), obtains denoising piecemeal correlation matrix

Step 6 is rightCarry out Eigenvalues Decomposition, obtain

Wherein U_SAnd U_NCorresponding signal space and spatial noise, U_S1And U_S2It is same dimension matrix, to U_S1Generalized inverse matrixWith U_S2ProductCarry out feature decomposition, obtain

Construction steering vectorAnd then obtain the estimation of steering vectorWherein For G^HGeneralized inverse matrix.

In force, to the denoising piecemeal correlation matrix having gotCarry out feature decomposition, obtain as formula (7) institute The multinomial showing, to therein to U_S1Generalized inverse matrix and U_S2ProductCarry out feature decomposition, obtain comprising T's Multinomial, further, obtains steering vectorAnd steering vector estimation

Step 7, chooses steering vectorIn m row p column element be a_mp, constructed fuction ρ₁(μ)、ρ₂(μ) as follows：

Wherein,If vectorial z_m=(x_m,y_m),x_mAnd y_mIt is m in the 1st subarray The transverse and longitudinal coordinate of individual sensor.If vectorial WithIt is the 2nd The transverse and longitudinal coordinate of m-th sensor in subarray.If vectorial μ=(cos α_psinβ_p,sinα_psinβ_p), if ＜ is z_m, μ ＞ represent to Amount z_mWith the inner product of μ, ifInt represents floor operation, to function ρ₁(μ)、ρ₂(μ) modulus value is carried out square, Have after arrangement

Wherein vectorial z_n=(x_n,y_n),x_nAnd y_nIt is the horizontal stroke of n-th sensor in the 1st subarray Vertical coordinate.If vectorial WithIt is n-th in the 2nd subarray The transverse and longitudinal coordinate of sensor.

| ρ₁(μ)|²、|ρ₂(μ)|²Middle and make derivative value be zero respectively to μ derivation, obtain

Choose zero point δ of maximum respectively₁And δ₂, solve from zero point is counter

＜ (z_m-z_n), μ ＞=angle (δ₁),

Substitute into z_m-z_n, μ, have equation below group

IfCan solve from above formula formulaCan get azimuth and the angle of pitch Estimation is respectively

In force, two function ρs related to the element in steering vector estimated matrix are constructed₁(μ)、ρ₂(μ), pass through To ρ₁(μ)、ρ₂(μ) a series of conversion, after derivation, obtain maximum zero point, from the anti-angle equation such as formula solving of this zero point (8) shown in, wherein, except Y_pBeyond the trigonometric function part representing, it is given value, therefore, by antitrigonometric function computing, that is, The estimated value of azimuth and the angle of pitch can be respectively obtained.

Following for the accuracy of checking this method, described in system model, in XOY face, arrangement has translation invariant Two sub- sensor arrays of structure, each submatrix contains 18 sensors, and adjacent sensors maximum spacing is less than half-wavelength, m= 18, three extraterrestrial targets respectively at (35 °, 30 °), (35 °, 40 °), it is average that (55 °, 30 °) each experiment is 100 Monte Carlos Result.

Test the angle root-mean-square error contrast that 1 the inventive method is estimated with 2D ESPRIT method

The mean square error that Fig. 2 and 3 is respectively angle estimation changes and corresponding Cramér-Rao lower bound with signal to noise ratio and hits (CRB) situation.In Fig. 1 hits be 500 times, signal to noise ratio change from 0dB to 50dB, in Fig. 2 signal to noise ratio be 10dB, hits from 100 to 1000 changes.It can be seen that in the case of signal to noise ratio and hits are less, the mean square error of angle estimation of the present invention is substantially better than 2D ESPRIT algorithm, with the increase of signal to noise ratio or hits, the inventive method and 2D ESPRIT algorithm performance convergence are simultaneously same When approach Cramér-Rao lower bound, be better than 2DESPRIT all the time in the whole mean square error emulating institute's extracting method in interval wider scope Algorithm.

Test the angle planisphere contrast that 2 the inventive method are estimated with 2D ESPRIT method

Figure 4 and 5 are respectively the planisphere of the two dimension angular of 3 targets of inventive algorithm and the estimation of 2D ESPRIT algorithm, Can intuitively find out that the space angle that inventive algorithm is estimated is more concentrated at being distributed near actual value from figure, therefore angle estimation Error less, in addition inventive algorithm is without separable dimension processing, without parameter pairing and overlapping in the one-dimensional component of two dimension angular When still can correctly estimate information source angle.

A kind of evaluation method of the space two-dimensional DOA proposing in the present embodiment, constructs cost letter by required steering vector Number, solves orientation and pitching two dimension angular by maximum zero point is counter, and it is corresponding all that the construction of steering vector considers two dimension angular The vandermonde structure of array element receiving data, the drawbacks of overcoming ESPRIT method and only determine angle by single feature value, improves The accuracy estimated.

It should be noted that：A kind of embodiment of the evaluation method of space two-dimensional DOA that above-described embodiment provides, is only used as Explanation in actual applications in this evaluation method, can also be according to actual needs and by above-mentioned evaluation method in other application field Use in scape, it implements process similar to above-described embodiment, repeats no more here.

The foregoing is only embodiments of the invention, not in order to limit the present invention, all in the spirit and principles in the present invention Within, any modification, equivalent substitution and improvement made etc., should be included within the scope of the present invention.

Claims (1)

1. a kind of evaluation method of space two-dimensional DOA is it is characterised in that include：

Step one, disposes receiving array, includes 2M omnidirectional's electromagnetic transducer in described receiving array, and array contains two submatrixs Row, each submatrix contains M sensor, and spacing d between two wherein adjacent sensors is less than 1/2nd of wavelength X；

Step 2, receives the first signal including the information source of P narrowband target in the very first time by described receiving array Value x₁(t)=A₁s(t)+n₁T (), described receiving array is pressed the first coordinate η=(η_x,η_y) translated, connecing after being translated Receive array, the receiving array after described translation receives secondary signal value x in the second time for the described information source₂(t)=A₂s (t)+n₂(t)；

Step 3, described first signal value and described secondary signal value are entered ranks storehouse, form high dimensional data

$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]=As\left(t\right)+n\left(t\right),$

And obtain the piecemeal correlation matrix of described x (t),

$R_x=\left[\begin{array}{cc}{A}_1{R}_s{A}_1^H+{\mathrm{\σ}}_1^{2}I& A_1{R}_s{\mathrm{\Λ}}^H{A}_1^H\\ A_1{\mathrm{\ΛR}}_s{A}_1^H& A_1{\mathrm{\ΛR}}_s{\mathrm{\Λ}}^H{A}_1^H+{\mathrm{\σ}}_2^{2}I\end{array}\right],$

After Eigenvalues Decomposition, averaged power spectrum is carried out to minimal eigenvalue, determine noiseApproximationAnd to described R_xCarry out denoising computing, obtain the piecemeal correlation matrix after denoising,

$C_x=R_x-\left[\begin{array}{cc}{\widehat{\mathrm{\σ}}}_1^{1}I& \\ & {\widehat{\mathrm{\σ}}}_2^{2}I\end{array}\right]=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right];$

Step 4, obtains random initial mask U (0)=[u₁(0) u₂(0) ... u_P(0) l], is made to be cycle-index, construction the One cost function

$\tilde{C}=\underset{i=1}{\overset{2}{\Σ}}\[C_{i1}^{H}U(l-1)U^H(l-1)C_{i1}\],$

WhereinC_i1Associate matrix each other, U^H(l-1) with U (l-1) associate matrix each other, to described first cost Function carries out feature decomposition, obtains fisrt feature breakdown

${\tilde{C}}_{eig}=\stackrel{\‾}{V}(l-1)\stackrel{\‾}{D}(l-1){\stackrel{\‾}{V}}^H(l-1),$

WhereinFor corresponding eigenvalue matrix,For corresponding eigenvectors matrix, if front P big eigenvalue Corresponding characteristic vector constitutes V (l-1), constructs the second cost function

WhereinC_i1Associate matrix each other, V^H(l-1) with V (l-1) associate matrix each other, to described second cost Function carries out feature decomposition, obtains second feature breakdown

By the first cost function described in looping construct and described second cost function, and carry out feature decomposition, until meeting

||U(L)U^H(L)-U(L-1)U^H(L-1)||_F≤ε

When stop circulation, wherein said ε is default threshold value, according to the l stopping during circulation, determines matrix U (L), definition letter Work song space G=U (L)；

Step 5, constructs piecemeal dimensionality reduction matrix J, described signal subspace G is the unit on diagonal in described piecemeal dimensionality reduction matrix J Element, by described piecemeal dimensionality reduction matrix J to the piecemeal correlation matrix C after described denoising_xCarry out dimension-reduction treatment, obtain

Described

Step 6, to describedCarry out Eigenvalues Decomposition, obtain

${\left({\tilde{C}}_x\right)}_{eig}=\left[\begin{array}{c}{U}_{S1}\\ U_{S2}\end{array}\right]{\mathrm{\Λ}}_S{\left[\begin{array}{c}{U}_{S1}\\ U_{S2}\end{array}\right]}^H+U_N{\mathrm{\Λ}}_N{U}_N^H,$

Wherein U_SAnd U_NCorresponding signal space and spatial noise, U_S1And U_S2It is same dimension matrix, to U_S1Generalized inverse matrixWith U_S2ProductCarry out feature decomposition, obtain

${\left(U_{S1}^c{U}_{S2}\right)}_{eig}=T^{-1}nT,$

Construction steering vectorAnd then obtain the estimation of described steering vectorWherein (G^H)^□ For G^HGeneralized inverse matrix；

Step 7, chooses described steering vectorIn m row p column element be a_mp, constructed fuction ρ₁(μ)、ρ₂(μ) as follows：

Wherein,If vectorial z_m=(x_m,y_m),x_mAnd y_mIt is m-th in the 1st subarray The transverse and longitudinal coordinate of sensor；If vectorial WithIt is the 2nd son The transverse and longitudinal coordinate of m-th sensor in array；If vectorial μ=(cos α_psinβ_p,sinα_psinβ_p), if ＜ is z_m, μ ＞ represent vector z_mWith the inner product of μ, ifInt represents floor operation, to described function ρ₁(μ)、ρ₂(μ) modulus value is put down Side, has after arrangement

Wherein vectorial z_n=(x_n, y_n),x_nAnd y_nIt is the transverse and longitudinal seat of n-th sensor in the 1st subarray Mark；If vectorial WithIt is n-th sensing in the 2nd subarray The transverse and longitudinal coordinate of device；Described | ρ₁(μ)|²、|ρ₂(μ)|²Middle and make derivative value be zero respectively to μ derivation, obtain

Choose zero point δ of maximum respectively₁And δ₂, solve from zero point is counter

＜ (z_m-z_n), μ ＞=angle (δ₁),

Substitute into z_m-z_n, μ, have equation below group

IfCan solve from above formulaCan get azimuth and the estimation of the angle of pitch divides It is not

${\widehat{\mathrm{\α}}}_p=\arctan (y_{p2}/y_{p1}),$

${\widehat{\mathrm{\β}}}_p=\arcsin \sqrt{y_{p1}^2+y_{p2}^2}.$