CN104407335B - DOA estimation method of 3-axis cross array - Google Patents

Abstract

The invention brings forward an SLS-NC-ESPRI algorithm based DOA estimation method of a 3-axis cross array, and is applied to estimation of the DOA of a strict second-order non-circular (NC) signal. Different from a conventional DOA estimation algorithm, the method brought forward by the invention has the following advantages: first of all, the by use of the NC characteristic of the signal, the virtual space of the array is expanded, the estimation precision is improved, and a detectable source number is enabled to be increased; secondly, the overlapping structure of sub-array configuration is utilized, the overlapping of signal subspace errors is taken into consideration, and shift-invariant equations in three cross array axes are effectively solved by use of a structure least square (SLS) method; and finally, through a mode of increasing constraint conditions, it is ensured that the three axial directions has approximately the same characteristic vectors, the rank defect problem introduced when arrival signals have the same projection at a certain axial direction is solved, the algorithmic validity is ensured, and a more accurate DOA estimated value is provided.

Description

A kind of doa method of estimation of 3 axle crossed arrays

Technical field

The present invention relates to array signal process technique field, more particularly, to a kind of doa method of estimation of 3 axle crossed arrays.

Background technology

Incoming wave signal doa estimates in radar, sonar, the multiple fields such as mobile communication and biomedical imaging obtain extensively Application.Numerous technology for solving high-resolution doa estimation problem are suggested, and wherein more classical technology has: Extremum search, root of a polynomial, matrix rotation is constant etc..As a kind of constant method of matrix rotation, esprit algorithm can be effective Ground goes out the doa of signal using the invariable rotary property calculation of signal subspace, and can provide analytic solutions.Generally, esprit Algorithm utilizes least square (least squares, ls) or total least square (total least squares, tls) to solve Doa estimate, has this feature of identical error because ls and tls method does not account for signal subspace lap, makes The doa estimate error that must obtain is larger.For the special construction information preferably being had using model, determined linear equation Algorithm is suggested the overall minimum norm of structuring (structured total least norm, stln) it is assumed that signal model is Ax ≈ b, there may be error in wherein a and b, stln remains the affine structure of a, for example, the structure such as Top profit thatch, Hunk that Type.Because the invariable rotary equation of esprit algorithm does not have the affine structure required for stln, therefore, stln is not suitable for Esprit (estimating signal parameters via rotational invariance techniques) algorithm Solution.Similar with stln method, sls algorithm also have ignored the quadratic term of the expansion of remaining matrix, can by emulation experiment To find, when being related to the configuration of overlapping subarrays, sls is more accurate than the doa estimate that ls and tls obtains.esprit Algorithm is initially that estimation range is [0, π] in order to determine one-dimensional orientation problem proposition.Doa in order to solve two dimension estimates to ask Topic, plane antenna structure is utilized, the array structure such as example circular, rectangle and chiasma type, and these array structures can be combined Estimate azimuth and the angle of pitch of signal.Due to the special construction of crossed array, compared to other multidimensional uniform array structures, It can provide the angular resolution of a larger aperture and Geng Gao.

In order to improve the performance of doa estimator further, miss except considering that signal subspace lap has identical This feature of difference, also uses the characteristic in time domain for the signal, i.e. nc characteristic.Not rounded signal such as bpsk, oqpsk, pam and ask adjust Signal processed, is all widely used in many Modern Communication System.In addition, when i/q lacks of proper care, circle signal also can be changed into non- How circle signal, so not rounded signal is very universal in actual applications, utilize the not rounded characteristic of signal, improves doa and estimates The performance of gauge has higher researching value.By using the not rounded characteristic of received signal, some are improved empty based on son Between doa estimate to be suggested, such as nc-music, nc-root-music, nc-esprit and nc-unitary esprit, Nc-esprit scheduling algorithm, this kind of algorithm can make array aperture double, and improves the precision that doa estimates.

Content of the invention

In order to overcome the problems referred to above, using signal subspace lap, there is identical error disturbance and arriving signal Not rounded characteristic, the present invention proposes a kind of 3 axle crossed array doa methods of estimation based on sls-nc-esprit algorithm.Due to rotation The signal subspace turning lap in constant equation (shift-invariant equations, sies) has identical by mistake Difference disturbance, sls algorithm can provide more accurate doa to estimate.

The present invention is achieved through the following technical solutions:

A kind of doa method of estimation of the 3 axle crossed arrays based on sls-nc-esprit algorithm, described 3 axle crossed arrays by Three along x-axis, y-axis, the linear sub-arrays of z-axis are formed with initial point for geometric center, all of bay be all independent and Equidistant, between the array element of 3 axles, distance is δ, and array number is m=m_x+m_y+m_z, have d far field narrow band signal to reach array, wherein m_x,m_y,m_zIt is illustrated respectively in x, y, the element number of array comprising of the subarray in z-axis direction；It is characterized in that: described inclusion is following Step:

1) x (t) represents reception sample data, using signal not rounded characteristic, obtains the measurement data vector expanded $x^{\left(\mathrm{nc}\right)}\left(t\right)=\left[\begin{array}{c}\mathrm{as}\left(t\right)\\ {\mathrm{\π}}_m{a}^{*}{s}^{*}\left(t\right)\end{array}\right]+\left[\begin{array}{c}n\left(t\right)\\ {\mathrm{\π}}_m{n}^{*}\left(t\right)\end{array}\right]=a^{\left(\mathrm{nc}\right)}s\left(t\right)+n^{\left(\mathrm{nc}\right)}\left(t\right),$ Wherein π_mSwitching matrix for m × m, this exchange The anti-diagonal element of matrix is 1, and other elements are 0, a is signal guide vector matrix, and s (t) is signal phasor, and n (t) is to make an uproar Sound；

2) calculate the sample covariance matrix of expansionWherein, x^(nc)Represent a 2m × n's Receiving data matrix, by n sampling snapshot data x^(nc)(t_n), 1≤n≤n composition；

3) set $\begin{array}{c}{j}_{\mathrm{\ξ}1}=[i_{m_{\mathrm{\ξ}}-1},0_{(m_{\mathrm{\ξ}}-1)\×1}]\\ j_{\mathrm{\ξ}2}=[0_{(m_{\mathrm{\ξ}}-1)\×1},i_{m_{\mathrm{\ξ}}-1}]\end{array},$ ξ ∈ { x, y, z }, now 3 axle crossed arrays there is the selection of Maximum overlap degree Matrix is as follows:

$k_{x1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{j}_{x1}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x2}{\mathrm{\π}}_{m_x}\end{array}\right]$

$k_{x2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{j}_{x2}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x1}{\mathrm{\π}}_{m_x}\end{array}\right]$

$k_{y1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y1}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_x}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_z}& 0_{(m_x-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y2}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]$

$k_{y2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y2}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_x}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_z}& 0_{(m_x-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y1}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]$

$k_{z1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_{z1}& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_z-1}{j}_{z2}{\mathrm{\π}}_{m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right]$

$k_{z2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_{z2}& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_z-1}{j}_{z1}{\mathrm{\π}}_{m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right];$

4) sample covariance matrix to expansionDo svd decomposition and obtain signal subspace estimationNot rounded is believed Number invariable rotary equation be represented by:

$k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\hat{u}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}=k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\hat{u}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\&element;\{x,y,z\},$

Least square method solution is carried out to above formula, obtainsξ ∈ { x, y, z } is as γ_ξInitial estimation, signal space Estimate initial value be

5) improved signal subspace is estimated to be expressed as:Obtain residual matrix $r_{\mathrm{\ξ}}({\tilde{u}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})=k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{u}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\&element;\{x,y,z\};$

6) utilize the array matrix γ of 3 axle bifurcated arrays_ξ, ξ ∈ { x, y, z } have the characteristics that identical characteristic vector set up Matrix f₁,f₂And f₃:

f₁=γ_xγ_y-γ_yγ_x=0_d×d

f₂=γ_yγ_z-γ_zγ_y=0_d×d

f₃=γ_zγ_x-γ_xγ_z=0_d×d；

7) for kth time iteration, remaining matrix, f₁,f₂And f₃Can be expressed as:

$r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}$

f_1k=γ_xkγ_yk-γ_ykγ_xk

f_2k=γ_ykγ_zk-γ_zkγ_yk

f_3k=γ_zkγ_xk-γ_xkγ_zk；

8) remaining matrix, f in the iteration of kth+1 time₁,f₂,f₃Can be written as:

$\begin{array}{c}r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}})\\ \≈r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})+k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}}+k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξk}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\end{array}$

f_1(k+1)=(γ_xk+δγ_xk)(γ_yk+δγ_yk)-(γ_yk+δγ_yk)(γ_xk+δγ_xk)

≈f_1k+γ_xkδγ_yk+δγ_xkγ_yk-γ_ykδγ_xk-δγ_ykγ_xk

f_2(k+1)=(γ_yk+δγ_yk)(γ_zk+δγ_zk)-(γ_zk+δγ_zk)(γ_yk+δγ_yk)

≈f_2k+γ_ykδγ_zk+δγ_ykγ_zk-γ_zkδγ_yk-δγ_zkγ_yk

f_3(k+1)=(γ_zk+δγ_zk)(γ_xk+δγ_xk)-(γ_xk+δγ_xk)(γ_zk+δγ_zk)

≈f_3k+γ_zkδγ_xk+δγ_zkγ_xk-γ_xkδγ_zk-δγ_xkγ_zk；

9) set up cost function:

$\{{\widehat{\mathrm{\γ}}}_x,{\widehat{\mathrm{\γ}}}_y,{\widehat{\mathrm{\γ}}}_z\}=\underset{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}{\arg \min }\left|\right|{\left[\begin{array}{c}\mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{xk}})\right\}\\ \mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{yk}})\right\}\\ \mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{zk}})\right\}\\ \mathrm{vec}\left\{f_{1k}\right\}\\ \mathrm{vec}\left\{f_{2k}\right\}\\ \mathrm{vec}\left\{f_{3k}\right\}\\ \widetilde{\mathrm{\κ}}\·\mathrm{vec}\left\{\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]+h_k\·\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]\left|\right|}_f,$

Wherein,It is the error matrix of the signal subspace estimated in kth step, γ_ξk= γ_ξk-1+δγ_ξk-1, ξ ∈ { x, y, z }, vec { } are vectorization functions,

$h_k=\left[\begin{array}{cccc}{i}_d\&circletimes;\left(k_{x1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(m_y-1)d\×d^2}& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^t\&circletimes;k_{x1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{x2}^{\left(\mathrm{nc}\right)}\\ 0_{2(m_x-1)d\×d^2}& i_d\&circletimes;\left(k_{y1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;k_{y1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{y2}^{\left(\mathrm{nc}\right)}\\ 0_{2(m_x-1)d\×d^2}& 0_{2(m_y-1)d\×d^2}& i_d\&circletimes;\left(k_{z1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& {\mathrm{\γ}}_{\mathrm{zk}}^t\&circletimes;k_{z1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{z2}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}& i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×d^2}& 0_{d^2\×2\mathrm{md}}\\ 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{zk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{zk}}& i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×2\mathrm{md}}\\ i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{xk}}& 0_{d^2\×2\mathrm{md}}\\ 0_{2\mathrm{md}\×d^2}& 0_{2\mathrm{md}\×d^2}& 0_{2\mathrm{md}\×d^2}& \widetilde{\mathrm{\κ}}{i}_{2\mathrm{md}}\end{array}\right];$

10) judge whether cost function meetsξ ∈ { x, y, z }, if met, recognizes For converging to optimal solution, stop iteration, wherein ε is the worst error scope allowing；If be unsatisfactory for, making k=k+1, returning Step 9；

11) by the optimal solution obtainingIt is evd to decompose, obtain μ_xp,μ_yp,μ_zp, p=1 ..., the optimal estimation of d ValueP=1 ..., d, wherein,

${\mathrm{\μ}}_x({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_y({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\sin \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_z({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right),$

λ₀For wavelength, θ_i,φ_iRepresent the angle of pitch and the azimuth of i-th signal respectively；

According to d θ_i,φ_iObtain the estimate of the doa of d information source

The invention has the beneficial effects as follows: for 3 axle crossed arrays, when two or more arriving signals are in x-axis, y-axis, z-axis Have during identical projection it may appear that rank defect problem, now ls, tls and sls method will lose efficacy.In order to solve the problems, such as rank defect, this Invention ensure that in the form of increasing constraints 3 direction of principal axis have approximately uniform characteristic vector, thus ensure that algorithm Validity it is provided that more accurate doa estimate.Because the former efficiently utilizes the not rounded characteristic of signal, the present invention Propose can be provided than sls-esprit algorithm more based on 3 axle crossed array doa methods of estimation of sls-nc-esprit algorithm Good performance.

Brief description

Fig. 1 is the structural representation of 3 axle crossed arrays；

Fig. 2 is the 3 axle crossed array doa method of estimation flow processs based on sls-nc-e esprit algorithm that the present invention provides Figure；

Fig. 3 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and signal to noise ratio (θ₁=θ₂,φ₁≠φ₂, n=300)；

Fig. 4 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and signal to noise ratio (θ₁=θ₂,φ₁≠φ₂, n=20)；

Fig. 5 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and fast umber of beats (θ₁=θ₂,φ₁≠φ₂)；

Fig. 6 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and signal to noise ratio (θ₁≠θ₂,φ₁≠φ₂, n=300)；

Fig. 7 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and signal to noise ratio (θ₁≠θ₂,φ₁≠φ₂, n=20)；

Fig. 8 is the graph of relation of the rmse of spatial frequency estimate with traditional algorithm for the method for the present invention and fast umber of beats (θ₁≠θ₂,φ₁≠φ₂)；

Fig. 9 is the relation of the method for the present invention and the rmse of spatial frequency estimate and azimuthal spacings φ of traditional algorithm Curve map；

Figure 10 is that the method for the present invention is bent with the relation that the elevation angle is spaced θ with the rmse of spatial frequency estimate of traditional algorithm Line chart.

Specific embodiment

The present invention is further described for explanation and specific embodiment below in conjunction with the accompanying drawings.

3 axle cross array structures as shown in Figure 1, by three along x-axis, y-axis, the linear sub-arrays of z-axis (with initial point are Geometric center) composition.Assume that all of bay is all independent and equidistant, between the array element of 3 axles apart from δ be identical. The array number assuming 3 axle crossed arrays is m, has d far field narrow band signal arrival array, m=m here_x+m_y+m_z, wherein m_x,m_y, m_zIt is illustrated respectively in x, y, the element number of array comprising of the subarray in z-axis direction.Array received signal modeling is:

X (t)=as (t)+n (t) (1)

Wherein s (t)=[s₁(t),...,s_d(t)]^tFor signal phasor, n (t)=[n₁(t),...,n_m(t)]^tComprise additivity Sensor noise；For signal guide vector matrix, comprise d signal guide vectorI=1 ..., d, wherein θ_i,φ_iRepresent the angle of pitch and the azimuth of i-th signal respectively；Information source number d commonly assumes that It is previously known or using information theory criterion, it can be estimated.The steering vector of 3 axle crossed arrays is as follows:

Wherein,

Vector p_n, n=1 ..., m represent aerial position, μ_iCan be expressed as

Wherein:

$\begin{array}{c}{\mathrm{\μ}}_x({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)\\ {\mathrm{\μ}}_y({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\sin \left({\mathrm{\φ}}_i\right)\sin \left({\mathrm{\θ}}_i\right)\\ {\mathrm{\μ}}_z({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}\cos \left({\mathrm{\φ}}_i\right)\end{array}---\left(5\right)$

Here λ₀For wavelength it is assumed that noise is white Gaussian noise process, its average is zero, and covariance matrix is It is incoherent with signal s (t).

Assume that signal is strict not rounded signal that is to say, that in planisphere, all possible states of signal can represent On one wire, such signal is widely present in actual applications, for example, the signal such as bpsk, oqpsk and pam.Signal is sweared Amount is decomposed into:

S (t)=ψ s₀(t) (6)

Whereinφ_i, i=1 ..., d are the initial phase of source signal, s₀T () is real-valued signal arrow Amount, signal s_iT the not rounded coefficient of () is defined as:

${\mathrm{\γ}}_{\mathrm{sl}}=\frac{e\left\{s_l^2\left(t\right)\right\}}{e\left\{{\left|s_l\left(t\right)\right|}^2\right\}},l=1,...,d---\left(7\right)$

According to aforementioned it is assumed that the covariance matrix of receipt signal vector x (t) can be expressed as:

$r_x=e\left\{x\left(t\right)x{\left(t\right)}^h\right\}=a{r}_s{a}^h+{\mathrm{\σ}}_n^2{i}_m---\left(8\right)$

Wherein,It is the noise power of the reception of each array antenna, r_s=e { s (t) s (t)^hIt is signal covariance square Battle array, to r_xCarry out svd and decompose and can obtain:

r_x=u σ u^h(9)

Wherein σ=diag { λ₁,...,λ_m, λ₁,…,λ_mFor receipt signal covariance matrix r_xCharacteristic value, and meetU=[u₁,...,u_m], i=1 ..., m is corresponding characteristic vector.HaveWherein σ_s=diag { λ₁,...,λ_dIt is dominant eigenvalue, u_s=[u₁,...,u_d] it is signal Subspace, u_n=[u_d+1,...,u_m] it is noise subspace, according to formula (7), can obtain:

$a{r}_s{a}^h+{\mathrm{\σ}}_n^2{i}_m=u_s{\mathrm{\σ}}_s{u}_s^h+{\mathrm{\σ}}_n^2{u}_n{u}_n^h---\left(10\right)$

WillUse in (10), obtain:

$u_s=a{r}_s{a}^h{u}_s{[{\mathrm{\σ}}_s-{\mathrm{\σ}}_n^2{i}_m]}^{-1}=\mathrm{at}---\left(11\right)$

Wherein,It is nonsingular matrix, u_sFor signal subspace.Formula (11) is esprit The key point that high-resolution doa estimate can be given of algorithm.

Assume that esprit algorithm uses the subarray with Maximum overlap, now selection matrix is:

$\begin{array}{c}{j}_{\mathrm{\ξ}1}=[i_{m_{\mathrm{\ξ}}-1},0_{(m_{\mathrm{\ξ}}-1)\×1}]\\ j_{\mathrm{\ξ}2}=[0_{(m_{\mathrm{\ξ}}-1)\×1},i_{m_{\mathrm{\ξ}}-1}]\end{array},\mathrm{\ξ}\&element;\{x,y,z\}---\left(12\right)$

Now the selection matrix of 3 axle crossed arrays is as follows:

$\begin{array}{c}{k}_{\mathrm{xl}}=[j_{\mathrm{xl}},0_{(m_x-1)\×m_y},0_{(m_x-1)\×m_z}],l=\mathrm{1,2}\\ k_{\mathrm{yl}}=[0_{(m_y-1)\×m_x},j_{\mathrm{yl}},0_{(m_y-1)\×m_z}],l=\mathrm{1,2}\\ k_{\mathrm{zl}}=[0_{(m_z-1)\×m_x},0_{(m_z-1)\×m_y},j_{\mathrm{zl}}],l=\mathrm{1,2}\end{array}---\left(13\right)$

Because all subarrays have the characteristic of invariable rotary, can get

k_ξ1aφ_ξ=k_ξ2a,ξ∈{x,y,z} (14)

WhereinCan be obtained according to formula (11):

k_ξ1u_sγ_ξ=k_ξ2u_s(15)

Wherein

γ_ξ=t^-1φ_ξt (16)

For the not rounded information using signal, widely used linear processing techniques, the measurement data vector x of definition expansion^(nc)T () is:

Wherein π_mSwitching matrix for m × m, its anti-diagonal element is 1, and other elements are 0.Formula (1) and (6) are substituted into (17), obtain:

$\begin{array}{c}{x}^{\left(\mathrm{nc}\right)}\left(t\right)=\left[\begin{array}{c}\mathrm{as}\left(t\right)\\ {\mathrm{\π}}_m{a}^{*}{s}^{*}\left(t\right)\end{array}\right]+\left[\begin{array}{c}n\left(t\right)\\ {\mathrm{\π}}_m{n}^{*}\left(t\right)\end{array}\right]\\ =a^{\left(\mathrm{nc}\right)}s\left(t\right)+n^{\left(\mathrm{nc}\right)}\left(t\right)\end{array}---\left(18\right)$

Wherein a^(nc)The array guiding matrix extending can be considered as, compared with former steering vector, it makes array virtual aperture It is changed into original twice, n^(nc)T () represents the noise vector of extension.It should be pointed out that array virtual aperture doubles to carry The resolution ratio that high doa estimates, and make maximum can detect that information source number doubles.

If x^(nc)Represent the receiving data matrix of a 2m × n, by n sampling snapshot data x^(nc)(t_n), 1≤n≤n group Become.Similar to esprit algorithm, the structure of selection matrix and expansion sampling matrix x here^(nc)Definition have relation, and these Matrix form (19) is to propose first in the present invention.

$\begin{array}{c}{k}_{x1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{j}_{x1}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x2}{\mathrm{\π}}_{m_x}\end{array}\right]\\ k_{x2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{j}_{x2}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x1}{\mathrm{\π}}_{m_x}\end{array}\right]\\ k_{y1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y1}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_x}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_z}& 0_{(m_x-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y1}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]\\ k_{y2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y2}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& j_{z1}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y1}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]\\ k_{z1}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_{z1}& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_z-1}{j}_{z2}{\mathrm{\π}}_{m_z}{}_{}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right]\\ k_{z2}^{\left(\mathrm{nc}\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_2& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}& {\mathrm{\π}}_{m_z-1}{j}_{z1}{\mathrm{\π}}_{m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right]\end{array}---\left(19\right)$

Each subarray due to 3 axle crossed arrays is invariable rotary, therefore same a^(nc)There is invariable rotary special Property, therefore meet following formula:

$k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{a}^{\left(\mathrm{nc}\right)}{\mathrm{\φ}}_{\mathrm{\ξ}}=k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{a}^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\&element;\{x,y,z\}---\left(20\right)$

According to formula (17), by reception data matrix x^(nc)Carry out svd decomposition, extract d maximum singular value and correspond to Characteristic vector, composition signal subspace estimationNow a^(nc)WithThere is approximately uniform column space, that is,Wherein t^(nc)∈c^d×dIt is nonsingular matrix.Similar with esprit, the sies of not rounded signal can represent For:

$k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\hat{u}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}=k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\hat{u}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\&element;\{x,y,z\}---\left(21\right)$

The modal method of solution formula (21) is the method based on ls, can be calculated γ using ls method_ξEstimate beξ ∈ { x, y, z }, however, ls algorithm is assumedIt is known and error free, Error exists only inIn, this situation obviously and in practical application does not correspond.Different from the algorithm of ls, tls Method is assumedWithAll there is error, be present in by minimizingWithInterferenceWithThus obtaining doa estimate more accurate than ls.But, only the array element in two subarrays does not have In the case of overlap, tls can provide the doa estimate of optimum.Conversely, working as along x, each linear array in y and z-axis direction When two subarrays have the array element of overlap, sls method due to considering the physical relationship existing in structure between subarray, because This, with respect to tls method, can improve the accuracy of estimate.So present invention discusses solving not rounded using sls algorithm The sies of signal, simulation results show, in the case of degree of overlapping maximum between subarray, sls algorithm substantially increases pitching Angle and azimuthal accuracy of estimation.

BecauseOnly open into the signal subspace comprising noise it is assumed that there is little disturbing in signal subspace DynamicImproved signal subspace is estimated to be expressed asDefining residual matrix is:

$r_{\mathrm{\ξ}}({\tilde{u}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})=k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_s^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{u}}_s^{\left(\mathrm{nc}\right)},\mathrm{\ξ}\&element;\{x,y,z\}---\left(22\right)$

Sls algorithm is makingFrobenius norm minimum on the basis of, keep simultaneously's Frobenius norm is as little as possible, in sum, sets up following cost function, is minimized by making table below reach formula, enters one The precision that step raising doa estimates:

$\underset{\mathrm{\δ}{\hat{u}}_s^{\left(\mathrm{nc}\right)},\mathrm{\δ}{\hat{u}}_s^{\left(\mathrm{nc}\right)}}{\min }{\left|\right|\left[\begin{array}{c}{r}_{\mathrm{\ξ}}({\tilde{u}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}})\\ \mathrm{\κ}\·\mathrm{\δ}{\hat{u}}_s^{\left(\mathrm{nc}\right)}\end{array}\right]\left|\right|}_f,\mathrm{\ξ}\&element;\{x,y,z\}---\left(23\right)$

WhereinIt is adjustmentAnd residual matrixBetween relation Regulatory factor.If the value of α is more than 1, in an iterative processEach element range of disturbance with respect toIn element big, but existing document verifies that sls algorithm is for the selection of α and insensitive.DefinitionWithIn (i, j) individual element be respectively r_ijAnd e_ij, formula (23) can be written as:

$\underset{\mathrm{\δ}{\hat{u}}_s^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}}}{\min }\underset{j=1}{\overset{d}{\mathrm{\σ}}}(\frac{\mathrm{\α}}{2(m-3)d}\underset{i=1}{\overset{2(m-3)}{\mathrm{\σ}}}{\left|r_{\mathrm{ij}}\right|}^2+\frac{1}{2\mathrm{md}}\underset{i=1}{\overset{2m}{\mathrm{\σ}}}{\left|e_{\mathrm{ij}}\right|}^2)---\left(24\right)$

Because the subarray on each axle of 3 axle crossed arrays is linear, when the signal of two or more arrival exists When identical projection is had on axle ξ, ξ ∈ { x, y, z }, matrixL=1,2 be likely to occur rank defect or The situation of approximate rank defect, therefore, ls, tls and sls algorithm will be unable to provide correct spatial frequency to estimate on specific ξ axle Value.But, because matrix γ_ξ, ξ ∈ { x, y, z } has identical characteristic vector composition matrix t, thus only exists one effectively Feasible solution.According to the knowledge of matrix analysis, the sufficient and necessary condition that matrix a and b has identical set of eigenvectors is ab =ba.Using this characteristic, the method for the present invention can be by with matrix γ_ξ, the relation between ξ ∈ { x, y, z } carries further The performance of high algorithm, defines matrix f₁,f₂And f₃, set up following equation:

f₁=γ_xγ_y-γ_yγ_x=0_d×d

f₂=γ_yγ_z-γ_zγ_y=0_d×d(25)

f₃=γ_zγ_x-γ_xγ_z=0_d×d

Formula (25) ensure that matrix γ_ξ, ξ ∈ { x, y, z } has identical set of eigenvectors, can solve when two or more The signal reaching has the rank defect problem leading to during identical projection on axle ξ, ξ ∈ { x, y, z }, thus the precision of the doa estimating It is improved.

By being iterated to formula (23) solving, minimize f simultaneously₁,f₂,f₃Frobenius norm be obtained with this Invent the improved sls-nc-esprit algorithm being adopted to solve the angle of pitch and the orientation angular estimation based on 3 axle crossed arrays Problem.

In the iteration of kth time, residual matrixf₁,f₂And f₃Can be expressed as:

$\begin{array}{c}r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ f_{1k}={\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{xk}}\\ f_{2k}={\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{yk}}\\ f_{3k}={\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{zk}}\end{array}---\left(26\right)$

Wherein ξ ∈ { x, y, z }. therefore, residual matrix in+1 iteration of kthf₁,f₂And f₃For:

$\begin{array}{c}r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})=r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}})\\ \≈r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})+k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}}+k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}{\mathrm{\γ}}_{\mathrm{\ξk}}-k_{\mathrm{\ξ}2}^{\left(\mathrm{nc}\right)}\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ f_{1(k+1)}=({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}})({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}})-({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}})({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}})\\ \≈f_{1k}+{\mathrm{\γ}}_{\mathrm{xk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}-\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{xk}}\\ f_{2(k+1)}=({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}})({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}})-({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}})({\mathrm{\γ}}_{\mathrm{yk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}})\\ \≈f_{2k}+{\mathrm{\γ}}_{\mathrm{yk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}-\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{yk}}\\ f_{3(k+1)}=({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}})({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}})-({\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}})({\mathrm{\γ}}_{\mathrm{zk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}})\\ \≈f_{3k}+{\mathrm{\γ}}_{\mathrm{zk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}-\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}{\mathrm{\γ}}_{\mathrm{zk}}\end{array}---\left(27\right)$

Notice that last approximate expression in (27) is ignored quadratic term and obtained.According to matrix analysis, for known square Battle arrayHave:

$\mathrm{vec}\left\{c_1{c}_2{c}_3\right\}=(c_3^t\&circletimes;c_1)\mathrm{vec}\left\{c_2\right\}---\left(28\right)$

Vectored calculations are applied in formula (27), obtain:

$\begin{array}{c}\mathrm{vec}\{r({\tilde{u}}_{s(k+1)}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξ}(k+1)}\}\≈\mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{\ξk}})\right\}\\ +[i_d\&circletimes;\left(k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right){\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}\right)]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}}\right\}]\\ +[{\mathrm{\γ}}_{\mathrm{\ξk}}^t\&circletimes;k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{\mathrm{\ξ}1}^{\left(\mathrm{nc}\right)}]\mathrm{vec}\left\{\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\\ \mathrm{vec}\left\{f_{1(k+1)}\right\}\≈\mathrm{vec}\left\{f_{1k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ +[i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{xk}}-{\mathrm{\γ}}_{\mathrm{xk}}\&circletimes;i_d]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{f_{2(k+1)}\right\}\≈\mathrm{vec}\left\{f_{2k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{zk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{zk}}]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ +[i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}\&circletimes;i_d]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{f_{3(k+1)}\right\}\≈\mathrm{vec}\left\{f_{3k}\right\}\\ +[{\mathrm{\γ}}_{\mathrm{xk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{xk}}\}\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ +[i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{zk}}-{\mathrm{\γ}}_{\mathrm{zk}}\&circletimes;i_d]\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\end{array}---\left(29\right)$

DefinitionIt is the error in kth time iteration with respect to the signal subspace initial value estimated Matrix, then the signal subspace estimate after improvingIt is represented by:

${\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}={\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{\tilde{u}}_{s,k}^{\left(\mathrm{nc}\right)}---\left(30\right)$

According to above iterative formula, nc-sls-esprit can be expressed with following formula:

$\{{\widehat{\mathrm{\γ}}}_x,{\widehat{\mathrm{\γ}}}_y,{\widehat{\mathrm{\γ}}}_z\}=\underset{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\δ}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}}{\arg \min }\left|\right|{\left[\begin{array}{c}\mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{xk}})\right\}\\ \mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{yk}})\right\}\\ \mathrm{vec}\left\{r({\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{u}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)},{\mathrm{\γ}}_{\mathrm{zk}})\right\}\\ \mathrm{vec}\left\{f_{1k}\right\}\\ \mathrm{vec}\left\{f_{2k}\right\}\\ \mathrm{vec}\left\{f_{3k}\right\}\\ \widetilde{\mathrm{\κ}}\·\mathrm{vec}\left\{\mathrm{\δ}{\tilde{u}}_{s,k}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]+h_k\·\left[\begin{array}{c}\mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{xk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{yk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{zk}}\right\}\\ \mathrm{vec}\left\{\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right\}\end{array}\right]\left|\right|}_f---\left(31\right)$

WhereinFor making vectorOptimization independent of other column vectors regulation The factor, wherein h_kIt is defined as follows:

$h_k=\left[\begin{array}{cccc}{i}_d\&circletimes;\left(k_{x1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(m_y-1)d\×d^2}& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^t\&circletimes;k_{x1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{x2}^{\left(\mathrm{nc}\right)}\\ 0_{2(m_x-1)d\×d^2}& i_d\&circletimes;\left(k_{y1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;k_{y1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{y2}^{\left(\mathrm{nc}\right)}\\ 0_{2(m_x-1)d\×d^2}& 0_{2(m_y-1)d\×d^2}& i_d\&circletimes;\left(k_{z1}^{\left(\mathrm{nc}\right)}{\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\right)& {\mathrm{\γ}}_{\mathrm{zk}}^t\&circletimes;k_{z1}^{\left(\mathrm{nc}\right)}-i_d\&circletimes;k_{z2}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}& i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×d^2}& 0_{d^2\×2\mathrm{md}}\\ 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{zk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{zk}}& i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×2\mathrm{md}}\\ i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{yk}}-{\mathrm{\γ}}_{\mathrm{yk}}^t\&circletimes;i_d& 0_{d^2\×d^2}& {\mathrm{\γ}}_{\mathrm{xk}}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{\mathrm{xk}}& 0_{d^2\×2\mathrm{md}}\\ 0_{2\mathrm{md}\×d^2}& 0_{2\mathrm{md}\×d^2}& 0_{2\mathrm{md}\×d^2}& \widetilde{\mathrm{\κ}}{i}_{2\mathrm{md}}\end{array}\right]---\left(32\right)$

By the cost function in optimized-type (31), obtain γ_x,γ_yAnd γ_zOptimal estimationWithRight WithCarry out feature decomposition, be calculated spatial frequency μ_xp,μ_yp,μ_zp, p=1 ..., the optimal estimation value of dp =1 ..., d.Note in an iterative process, by carrying out deformation process to the vector obtaining, to recover by following formulaWith γ_ξk:

$\begin{array}{c}{\tilde{u}}_{s(k+1)}^{\left(\mathrm{nc}\right)}={\tilde{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}+\mathrm{\δ}{\hat{u}}_{\mathrm{sk}}^{\left(\mathrm{nc}\right)}\\ {\mathrm{\γ}}_{\mathrm{\ξ}(k+1)}={\mathrm{\γ}}_{\mathrm{\ξk}}+\mathrm{\δ}{\mathrm{\γ}}_{\mathrm{\ξk}},\mathrm{\ξ}\&element;\{x,y,z\}\end{array}---\left(33\right)$

The initial value that signal subspace is estimated is made to be the signal subspace estimate in svd decomposition, that is, By obtained using ls methodAs γ_ξInitial value.When cost function meets following constraintsIt is believed that improved sls-nc-esprit of the present invention during ξ ∈ { x, y, z } To optimal solution, wherein ε is the worst error scope allowing to algorithmic statement.

Accompanying drawing 2 is the 3 axle crossed array doa method of estimation streams based on sls-nc-e esprit algorithm that the present invention provides Cheng Tu.

Step 1: obtain the receipt signal vector x expanding using signal not rounded characteristic^(nc)(t), and calculate the sampling of expansion Covariance matrix

Step 2: decomposed using svdObtain signal subspace to estimateAssume that improved signal subspace is estimated Can be expressed as:Obtain residual matrix

Step 3: using the array matrix γ of 3 axle bifurcated arrays_ξ, ξ ∈ { x, y, z } has the characteristics that identical characteristic vector Set up matrix f₁,f₂And f₃As shown in formula (25)；

Step 4: cost function is set up by formula (31), when cost function meets following constraintsIt is believed that converging to optimal solution during ξ ∈ { x, y, z }, wherein ε is that the maximum allowing is missed Difference scope；

Step 5: obtained by formula (31)It is done with evd decompose, d information source be can be obtained by according to formula (5) The estimate of doa.

The present invention mainly discusses spatial frequency μ to incoming wave signal_xp,μ_yp,μ_zp, p=1 ..., the pushing away of the estimate of d Lead process.But independently carry out because the svd of each axle decomposes, so the estimation of their calculated spatial frequencys Value is not one-to-one.The marriage problem of spatial frequency is all discussed in many documents, because not being the core of the present invention Content, so do not describe in detail.Find angle of arrival and the angle of pitch of signal, i.e. θ by formula (5)_pWithP=1 ... d, be By spatial frequency μ_xp,μ_yp,μ_zpWell-determined, the spatial frequency estimate obtaining with the inventive method can be unique Determine the doa of incoming wave signal.

Consider that one 3 axle crossed arrays comprise m=3m array element, along x, each subarray of y and z-axis all has m array element.Root Estimate rmse according to formula (34)

${\mathrm{rmse}}_p=\sqrt{\frac{1}{m_c}\underset{i=1}{\overset{m_c}{\mathrm{\σ}}}\underset{r=x,y,z}{\mathrm{\σ}}{({\widehat{\mathrm{\μ}}}_{p,i}^{\left(r\right)}-{\mathrm{\μ}}_p^{\left(r\right)})}^2}---\left(34\right)$

Wherein Monte Carlo experiment number of times is m_c=1000,P=1 ..., d tests p-th for i ＆ lt The estimation space frequency estimation of signal, andP=1 ..., d is corresponding spatial frequency actual value.Institute It is assumed that there being the bpsk signal of two equal power in some examples, not rounded coefficient is 1, and initial phase is respectively 5 °, 20 °.Battle array Unit's sum is m=30, and the interval between the array element of three subarrays is all equal to be 0.45 λ₀, λ here₀For incoming wave signal wavelength. The error margin factor of sls algorithm is ε=10^-6.In order to verify sls-nc-esprit of the present invention in l-G simulation test The performance advantage of algorithm, by itself and esprit, sls-esprit, nc-esprit algorithm has carried out Performance comparision.Attached in order to study The impact to sls-esprit and sls-nc-esprit algorithm for the addition of constraints condition, also provide in simulation result not additional about The characteristic of sls-esprit and sls-nc-esprit of bundle condition, in analogous diagram, uses sls-esprit (no respectively Constraint) represent with sls-nc-esprit (no constraint).

θ is worked as in experiment 1₁=θ₂,φ₁≠φ₂When rmse characteristic

Work as θ₁=θ₂When, μ can be obtained by formula (5)_z1=μ_z2, that is, two signals have identical to project in z-axis.According to formula (21) Understand rank defect problem occurs as ξ=z, now ls and do not have the sls method of additional constraint condition all can lose efficacy, will lead to Esprit, nc-esprit and the sls-esprit poor performance not having additional constraint condition are it is assumed that the azimuth of two bpsk signals It is respectively 0 ° and 20 °, the angle of pitch is 60 °.

By accompanying drawing 3 it is found that being fixed as n=300 in fast umber of beats, the rmse of the spatial frequency estimate of distinct methods with The change of snr, the change with snr is big, does not have the sls-esprit algorithm performance of additional constraint condition to decline notable, esprit and The rmse of nc-esprit algorithm is also larger, because they cannot solve the problems, such as rank defect.Sls-nc- proposed by the present invention Esprit algorithm performance in all algorithms is best.Subtracting fewer snapshots is n=20, from accompanying drawing 4, proposed by the present invention The performance of sls-nc-esprit algorithm remains optimum.When fixing snr is -5db, accompanying drawing 5 illustrates each algorithm output rmse With the change of fast umber of beats, esprit and do not have additional constraint condition sls-esprit algorithm performance worst, when fast umber of beats increases, Performance is not almost improved, although nc-esprit and do not have the sls-nc-esprit algorithm performance of additional constraint condition to carry Rise, but still can not solve the problems, such as rank defect, sls-esprit and sls-nc-esprit algorithm proposed by the present invention have preferably Performance, but because sls-nc-esprit proposed by the present invention considers the not rounded characteristic of signal, rmse is minimum for its output.

θ is worked as in experiment 2₁≠θ₂,φ₁≠φ₂When rmse characteristic

Fixing fast umber of beats n=300, the azimuth coverage of two bpsk signals is (- 4 °, 4 °), elevation coverage be (18 °, 28°).From accompanying drawing 6, sls-nc-esprit algorithm proposed by the present invention is more excellent than other algorithm performances, particularly in noise Become apparent from than less than performance advantage in the case of -5db.Although now there is not rank defect problem, sls-esprit and the present invention The sls-nc-esprit algorithm proposing still has less output rmse than do not have an additional constraint condition based on sls algorithm. In fig. 7, subtracting fewer snapshots is n=20, now still it can be found that sls-nc-esprit performance proposed by the present invention is better than it His algorithm.Accompanying drawing 8 illustrates each algorithm and exports the change with fast umber of beats for the rmse, using the nc characteristic of signal nc-esprit and The performance of sls-nc-esprit algorithm is improved, but it is found that sls-nc-esprit proposed by the present invention Algorithm performance in all algorithms is optimum.

Experiment 3rmse is with the relation at azimuth and the interval variation of the angle of pitch

Fixing fast umber of beats n=300, signal to noise ratio snr=-5db, the doas of first aim is (2 °, 18 °), second mesh Target doas is (2 ° of+δ φ, 28 °), wherein δ φ ∈ (2 °, 10 °).Output rmse is illustrated with two aspects by accompanying drawing 9 The change curve at the interval between angle, as seen from the figure, sls-nc-esprit algorithm performance proposed by the present invention is better than other Algorithm.In fig. 10, fast umber of beats and signal to noise ratio adopt and accompanying drawing 9 identical parameter setting.The doas of fixing first aim For (- 4 °, 18 °), the doas of second target is (4 °, 18 ° of+δ θ), δ θ ∈ (2 °, 10 °), and observation finds to be less than when the angle of pitch When 5 °, the performance of sls-nc-esprit algorithm proposed by the present invention obtains being significantly better than other algorithms.Due to esprit and nc- Esprit have ignored the overlapping features in structure between subarray, and therefore they provide poor estimated accuracy, although sls- Esprit algorithm can make full use of the overlay structure between subarray, but due to not utilizing the nc characteristic of signal, therefore its characteristic ratio Sls-nc-esprit algorithm proposed by the present invention is poor, particularly in the case of 5 ° of δ θ ＜.

Above content is to further describe it is impossible to assert with reference to specific preferred embodiment is made for the present invention Being embodied as of the present invention is confined to these explanations.For general technical staff of the technical field of the invention, On the premise of present inventive concept, some simple deduction or replace can also be made, all should be considered as belonging to the present invention's Protection domain.

Claims (3)

1. a kind of doa method of estimation of the 3 axle crossed arrays based on sls-nc-esprit algorithm, described 3 axle crossed arrays are by three Individual along x-axis, y-axis, the linear sub-arrays of z-axis are formed with initial point for geometric center, and all of bay is all independence and waits Away from, between the array element of 3 axles, distance is δ, and array number is m=m_x+m_y+m_z, have d far field narrow band signal to reach array, wherein m_x,m_y,m_zIt is illustrated respectively in x, y, the element number of array comprising of the linear sub-arrays in z-axis direction；

It is characterized in that: described comprise the following steps:

1) x (t) represents reception sample data, using signal not rounded characteristic, obtains the measurement data vector expandedWherein π_mSwitching matrix for m × m, this exchange The anti-diagonal element of matrix is 1, and other elements are 0, a is signal guide vector matrix, and s (t) is signal phasor, and n (t) is to make an uproar Sound；

2) calculate the sample covariance matrix of expansionWherein, x^(nc)Represent the reception of a 2m × n Data matrix, by n sampling snapshot data x^(nc)(t_n), 1≤n≤n composition；

3) setξ ∈ { x, y, z }, now 3 axle crossed arrays there is the selection matrix of Maximum overlap degree such as Under:

$k_{x1}^{\left(nc\right)}=\left[\begin{array}{cccccc}{j}_{x1}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x2}{\mathrm{\π}}_{m_x}\end{array}\right]$

$k_{x2}^{\left(nc\right)}=\left[\begin{array}{cccccc}{j}_{x2}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_x}\\ 0_{(m_x-1)\×m_x}& 0_{(m_x-1)\×m_y}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_z}& 0_{(m_x-1)\×m_y}& {\mathrm{\π}}_{m_x-1}{j}_{x1}{\mathrm{\π}}_{m_x}\end{array}\right]$

$k_{y1}^{\left(nc\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y1}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_x}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_z}& 0_{(m_x-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y2}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]$

$k_{y2}^{\left(nc\right)}=\left[\begin{array}{cccccc}{0}_{(m_y-1)\×m_x}& j_{y2}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_z}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_x}\\ 0_{(m_y-1)\×m_x}& 0_{(m_y-1)\×m_y}& 0_{(m_y-1)\×m_z}& 0_{(m_x-1)\×m_z}& {\mathrm{\π}}_{m_y-1}{j}_{y1}{\mathrm{\π}}_{m_y}& 0_{(m_y-1)\×m_z}\end{array}\right]$

$k_{z1}^{\left(nc\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_{z1}& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_z-1}{j}_{z2}{\mathrm{\π}}_{m_2}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right]$

$k_{z2}^{\left(nc\right)}=\left[\begin{array}{cccccc}{0}_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& j_{z2}& 0_{(m_z-1)\×m_z}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\\ 0_{(m_z-1)\×m_x}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_z}& {\mathrm{\π}}_{m_z-1}{j}_{z1}{\mathrm{\π}}_{m_2}& 0_{(m_z-1)\×m_y}& 0_{(m_z-1)\×m_x}\end{array}\right];$

4) sample covariance matrix to expansionDo svd decomposition and obtain signal subspace estimationThe rotation of not rounded signal Turn constant equation to be represented by:

$k_{\mathrm{\ξ}1}^{\left(nc\right)}{\hat{u}}_s^{\left(nc\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}=k_{\mathrm{\ξ}2}^{\left(nc\right)}{\hat{u}}_s^{\left(nc\right)},\mathrm{\ξ}\&element;\{x,y,z\},$

Least square method solution is carried out to above formula, obtainsξ ∈ { x, y, z } is as γ_ξInitial estimation, signal space estimate Initial value beγ_ξNot rounded coefficient for not rounded signal；

5) improved signal subspace is estimated to be expressed as:Obtain residual matrixξ∈{x,y,z}；

6) utilize the array matrix γ of 3 axle bifurcated arrays_ξ, ξ ∈ { x, y, z } has the characteristics that identical characteristic vector sets up matrix f₁,f₂And f₃:

f₁=γ_xγ_y-γ_yγ_x=0_d×d

f₂=γ_yγ_z-γ_zγ_y=0_d×d

f₃=γ_zγ_x-γ_xγ_z=0_d×d；

7) for kth time iteration, remaining matrix, f₁,f₂And f₃Can be expressed as:

$r({\tilde{u}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{\mathrm{\ξ}k})=k_{\mathrm{\ξ}1}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}{\mathrm{\γ}}_{\mathrm{\ξ}}-k_{\mathrm{\ξ}2}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}$

f_1k=γ_xkγ_yk-γ_ykγ_xk

f_2k=γ_ykγ_zk-γ_zkγ_yk

f_3k=γ_zkγ_xk-γ_xkγ_zk；

8) remaining matrix, f in the iteration of kth+1 time₁,f₂,f₃Can be written as:

$\begin{array}{c}r({\tilde{u}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{\mathrm{\ξ}k})=r({\tilde{u}}_{sk}^{\left(nc\right)}+\mathrm{\δ}{\hat{u}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{\mathrm{\ξ}k}+{\mathrm{\δ\γ}}_{\mathrm{\ξ}k})\\ \≈r({\tilde{u}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{\mathrm{\ξ}k})+k_{\mathrm{\ξ}1}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}{\mathrm{\δ\γ}}_{\mathrm{\ξ}k}+k_{\mathrm{\ξ}1}^{\left(nc\right)}\mathrm{\δ}{\tilde{u}}_{sk}^{\left(nc\right)}{\mathrm{\γ}}_{\mathrm{\ξ}k}-k_{\mathrm{\ξ}2}^{\left(nc\right)}\mathrm{\δ}{\tilde{u}}_{sk}^{\left(nc\right)}\end{array}$

$\begin{array}{c}{f}_{1(k+1)}=({\mathrm{\γ}}_{xk}+{\mathrm{\δ\γ}}_{xk})({\mathrm{\γ}}_{yk}+{\mathrm{\δ\γ}}_{yk})-({\mathrm{\γ}}_{yk}+{\mathrm{\δ\γ}}_{yk})({\mathrm{\γ}}_{xk}+{\mathrm{\δ\γ}}_{xk})\\ \≈f_{1k}+{\mathrm{\γ}}_{xk}{\mathrm{\δ\γ}}_{yk}+{\mathrm{\δ\γ}}_{xk}{\mathrm{\γ}}_{yk}-{\mathrm{\γ}}_{yk}{\mathrm{\δ\γ}}_{xk}-{\mathrm{\δ\γ}}_{yk}{\mathrm{\γ}}_{xk}\end{array}$

$\begin{array}{c}{f}_{2(k+1)}=({\mathrm{\γ}}_{yk}+{\mathrm{\δ\γ}}_{yk})({\mathrm{\γ}}_{zk}+{\mathrm{\δ\γ}}_{zk})-({\mathrm{\γ}}_{zk}+{\mathrm{\δ\γ}}_{zk})({\mathrm{\γ}}_{yk}+{\mathrm{\δ\γ}}_{yk})\\ \≈f_{2k}+{\mathrm{\γ}}_{yk}{\mathrm{\δ\γ}}_{zk}+{\mathrm{\δ\γ}}_{yk}{\mathrm{\γ}}_{zk}-{\mathrm{\γ}}_{zk}{\mathrm{\δ\γ}}_{yk}-{\mathrm{\δ\γ}}_{zk}{\mathrm{\γ}}_{yk}\end{array}$

$\begin{array}{c}{f}_{3(k+1)}=({\mathrm{\γ}}_{zk}+{\mathrm{\δ\γ}}_{zk})({\mathrm{\γ}}_{xk}+{\mathrm{\δ\γ}}_{xk})-({\mathrm{\γ}}_{xk}+{\mathrm{\δ\γ}}_{xk})({\mathrm{\γ}}_{zk}+{\mathrm{\δ\γ}}_{zk})\\ \≈f_{3k}+{\mathrm{\γ}}_{zk}{\mathrm{\δ\γ}}_{xk}+{\mathrm{\δ\γ}}_{zk}{\mathrm{\γ}}_{xk}-{\mathrm{\γ}}_{xk}{\mathrm{\δ\γ}}_{zk}-{\mathrm{\δ\γ}}_{xk}{\mathrm{\γ}}_{zk}\end{array};$

9) set up cost function:

$\{{\widehat{\mathrm{\γ}}}_x,{\widehat{\mathrm{\γ}}}_y,{\widehat{\mathrm{\γ}}}_z\}=\underset{{\mathrm{\δ\γ}}_{\mathrm{\ξ}k},\mathrm{\δ}{\tilde{u}}_{sk}^{\left(nc\right)}}{\arg \min }\left|\right|\left[\begin{array}{c}vec\left\{r({\tilde{u}}_{sk}^{\left(nc\right)}+{\mathrm{\δu}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{xk})\right\}\\ vec\left\{r({\tilde{u}}_{sk}^{\left(nc\right)}+{\mathrm{\δu}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{yk})\right\}\\ vec\left\{r({\tilde{u}}_{sk}^{\left(nc\right)}+{\mathrm{\δu}}_{sk}^{\left(nc\right)},{\mathrm{\γ}}_{zk})\right\}\\ vec\left\{f_{1k}\right\}\\ vec\left\{f_{2k}\right\}\\ vec\left\{f_{3k}\right\}\\ \widetilde{\mathrm{\κ}}\·vec\left\{\mathrm{\δ}{\tilde{u}}_{s,k}^{\left(nc\right)}\right\}\end{array}\right]+h_k\·\left[\begin{array}{c}vec\left\{{\mathrm{\δ\γ}}_{xk}\right\}\\ vec\left\{{\mathrm{\δ\γ}}_{yk}\right\}\\ vec\left\{{\mathrm{\δ\γ}}_{zk}\right\}\\ vec\left\{\mathrm{\δ}{\hat{u}}_{sk}^{\left(nc\right)}\right\}\end{array}\right]|{|}_f,$

Wherein,It is the error matrix of the signal subspace estimated in kth step, γ_ξk=γ_ξk-1+δ γ_ξk-1, ξ ∈ { x, y, z }, vec { } are vectorization functions,Wherein, α is a numerical value；

$h_k=\left[\begin{array}{cccc}{i}_d\&circletimes;\left(k_{x1}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}\right)& 0_{2(m_y-1)d\×d^2}& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{xk}^t\&circletimes;k_{x1}^{\left(nc\right)}-i_d\&circletimes;k_{x2}^{\left(nc\right)}\\ 0_{2(m_x-1)d\×d^2}& i_d\&circletimes;\left(k_{y1}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}\right)& 0_{2(m_z-1)d\×d^2}& {\mathrm{\γ}}_{yk}^t\&circletimes;k_{y1}^{\left(nc\right)}-i_d\&circletimes;k_{y2}^{\left(nc\right)}\\ 0_{2(m_x-1)d\×d^2}& 0_{2(m_y-1)d\×d^2}& i_d\&circletimes;\left(k_{z1}^{\left(nc\right)}{\tilde{u}}_{sk}^{\left(nc\right)}\right)& {\mathrm{\γ}}_{zk}^t\&circletimes;k_{z1}^{\left(nc\right)}-i_d\&circletimes;k_{z2}^{\left(nc\right)}\\ {\mathrm{\γ}}_{yk}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{yk}& i_d\&circletimes;{\mathrm{\γ}}_{yk}-{\mathrm{\γ}}_{yk}^t\&circletimes;i_d& 0_{d^2\×d^2}& 0_{d^2\×2md}\\ 0_{d^2\×d^2}& {\mathrm{\γ}}_{zk}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{zk}& i_d\&circletimes;{\mathrm{\γ}}_{yk}-{\mathrm{\γ}}_{yk}^t\&circletimes;i_d& 0_{d^2\×2md}\\ i_d\&circletimes;{\mathrm{\γ}}_{yk}-{\mathrm{\γ}}_{yk}^t\&circletimes;i_d& 0_{d^2\×d^2}& {\mathrm{\γ}}_{xk}^t\&circletimes;i_d-i_d\&circletimes;{\mathrm{\γ}}_{xk}& 0_{d^2\×2md}\\ 0_{2md\×d^2}& 0_{2md\×d^2}& 0_{2md\×d^2}& \widetilde{\mathrm{\κ}}{i}_{2md}\end{array}\right];$

10) judge whether cost function meetsξ ∈ { x, y, z }, if met, thinks receipts Hold back optimal solution, stop iteration, wherein ε is the worst error scope allowing；If be unsatisfactory for, make k=k+1, return to step 9；

11) by the optimal solution obtainingIt is evd to decompose, obtain spatial frequency μ_xp,μ_yp,μ_zp, the optimum of p=1 ..., d estimates EvaluationP=1 ..., d, wherein,

${\mathrm{\μ}}_x({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}cos\left({\mathrm{\φ}}_i\right)sin\left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_y({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}sin\left({\mathrm{\φ}}_i\right)sin\left({\mathrm{\θ}}_i\right)$

${\mathrm{\μ}}_z({\mathrm{\θ}}_i,{\mathrm{\φ}}_i)=\frac{2\mathrm{\π}}{{\mathrm{\λ}}_0}\mathrm{\δ}cos\left({\mathrm{\φ}}_i\right),$

λ₀For wavelength, θ_i,φ_iRepresent the angle of pitch and the azimuth of i-th signal respectively；

12) according to d θ_i,φ_iObtain the estimate of the doa of d information source.

2. method according to claim 1 it is characterised in that: the not rounded characteristic using signal improves doa estimated accuracy And make to detect that information source number increases, using 3 axle bifurcated array array matrix γ_ξ, ξ ∈ { x, y, z } has identical characteristic vector Feature effectively solving rank defect problem, thus provide more accurate doa estimating.

3. method according to claim 1 it is characterised in that: described δ γ_ξk-1, ξ ∈ { x, y, z } andIteration at the beginning of Initial value be appoint to a value.