CN104730491A - Virtual array DOA estimation method based on L type array - Google Patents

Abstract

The invention discloses a virtual array DOA estimation method based on an L type array. The method comprises the steps that 1, based on the shift invariant property, a subarray Zx and a subarray Zy of the L type array horizontally shift to obtain a virtual array Zx' and a virtual array Zy', rotation invariance of two sub signals is formed due to the shift invariant property of the subarrays, and the virtual subarray signals are equal to the L type subarray Zx and L type subarray Zy input signals multiplied by a twiddle factor respectively; 2, the output of the four subarrays are combined to form a virtual array output signal Z (t); 3, the signal subspace and the noise subspace can be described by decomposing the features of covariance matrixes output by the array, mutual correlation processing is carried out on the array output signal Z(t) to obtain Rzz, and eigenvalue decomposition is carried out to obtain signal subspaces; 4, the twiddle factor is solved through linear operation, and the signal wave arrival direction can be obtained through the diagonal element of the twiddle factor. According to the virtual array DOA estimation method, no spectral function needs to be calculated, the phenomenon that the wave arrival direction is indirectly calculated without searching for the peak value is avoided, the complexity is lowered, the equipment complexity and cost are reduced, and the positioning precision is high.

Description

A kind of virtual array DOA estimation method based on L-type battle array

Technical field

The invention belongs to Array Signal Processing arrival direction estimation technical field, in particular a kind of Virtual array DOA estimation method based on L-type array.

Background technology

Target Bearing Estimation (DOA) estimation is an important branch in the signal transacting fields such as sonar, radar, radio communication, medical imaging, microphone array column processing, simultaneously the basic problem that DOA solves determines to be in the locus of multiple interested echo signal in a certain space field, and namely each echo signal arrives the deflection of sensor array.Direction-finding method based on conventional beamformer scanning has intrinsic restriction, and be subject to the restriction of " Rayleigh criterion ", the resolution of estimation depends on array length, only when the separation angle in space between two information sources is greater than array aperture reciprocal, just can be resolved.In order to improve the resolution of basic matrix, when signal wavelength lambda must be, generally only has the aperture length increasing basic matrix, namely increase array element m or increase battle array spacing d, but increasing element number of array can improve equipment complexity and cost, increases array element distance and will cause secondary maximum again, for the restriction of actual conditions, the yardstick of array also can not do very large, therefore only relies on increase array aperture to reach and puies forward high-resolution way, be difficult to be suitable in Practical Project.In order to overcome this limitation, devise a kind of virtual array DOA estimation method based on L-type battle array, its basic thought is under the array case of limited dimensions, is obtained the virtual expansion of array by optimized algorithm, thus improves resolution.

Summary of the invention

The present invention is directed to the deficiencies in the prior art, propose a kind of virtual array DOA estimation method based on L-type battle array.

Technical scheme of the present invention is as follows:

Based on a virtual array DOA estimation method for L-type array, wherein, comprise the following steps:

Step 1: structure L-type array, determines the signal model of array received;

Receive acoustic pressure time-domain signal by 2M-1 array element in L-type array, this L-type array is the even linear array Z of M by array number in x-axis _xwith the even linear array Z that array number in y-axis is M _yform, wherein, 2M-1 is L-type array elements number, M be not less than 2 integer, d is array element distance.Hypothesis space has K information source to incide on array, and its 2-d direction finding is be respectively the elevation angle and the position angle of a kth signal source.

Suppose that the signal number incided on this array is K, then signal x-axis, y-axis received by M array element is respectively respectively as shown in the formula (1) and formula (2):

x(t)＝A _xs(t)+n(t) (1)

y(t)＝A _ys(t)+n(t) (2)

In formula, s (t) is signal source matrix, and n (t) is noise matrix, A _x, A _y∈ C ^{m × K}, be respectively the direction matrix in L-type array x-axis, y-axis, can be expressed as:

Step 2: obtain output signal matrix Z by ESPRIT constructing virtual array;

By the submatrix Z of ESPRIT by L battle array _x, Z _ycarrying out virtual expansion is submatrix Z _x', Z _y', the motion immovability due to submatrix defines the rotational invariance of two submatrix signals, i.e. Z _x'submatrix signal be actual submatrix Z _xinput signal be multiplied by twiddle factor φ _xobtain, Z _y'submatrix signal be actual submatrix Z _yinput signal be multiplied by twiddle factor φ _yobtain, first obtained the output signal of virtual submatrix by formula (1) and formula (2), then four submatrixs are exported and merged, form output signal matrix z (t) of whole array as shown in the formula (3):

$z\left(t\right)={[x\left(t\right),x^{\′}\left(t\right)y^{\′}\left(t\right)]}^T=\stackrel{\‾}{A}s\left(t\right)+n\left(t\right)$

Wherein suppose that the direction of arrival of each information source is different, then column vector between Line independent,

And

Wherein, matrix φ _x, φ _yfor the diagonal matrix of K × K, its diagonal element is that signal is respectively at Z _x, Z _yphase delay on array between any array element, diag represents diagonal matrix, and the element namely except principal diagonal is the square formation of zero.

X-axis direction homogenous linear submatrix Z is comprised such as formula (3) z (t) _xoutput signal x (t) _,y-axis direction homogenous linear submatrix Z _youtput signal y (t), Z _xthe virtual submatrix Z that translation obtains _x'output signal x'(t), Z _ythe virtual submatrix Z that translation obtains _y'output signal y'(t), the noise that each array received arrives is identical, virtual submatrix Z _x', Z _y'the even linear array of to be all array number be M.

Step 3: obtain correlation matrix R from array output signal matrix Z _z,

Signal subspace and noise subspace can obtain, shown in (4) by the feature decomposition that array exports the covariance matrix of Z:

$R_z=E\left[z\left(t\right)z^H\left(t\right)\right]=\stackrel{\‾}{A}{R}_s{\stackrel{\‾}{A}}^H+{\mathrm{\σ}}^{2}I---\left(4\right)$

R in formula _sfor the autocorrelation matrix of signal, σ ²for noise variance, I is unit matrix, the E [.] in formula (4), (.) ^hbe expressed as mathematical expectation, conjugate transpose operation.

Step 4: by correlation matrix R _zdo feature decomposition, estimated signal number;

Array correlation matrix R _ztwo spaces can be divided into, namely the eigenwert characteristic of correspondence vector E of K _s=[s ₁, s ₂... s _k] composition signal subspace, the non-singular matrix T that there is a K × K meets and move invariant feature E due to array _scan be analyzed to 4 parts, E _x, E _y, E _x', E _y'∈ C ^{m × K}, corresponding subarray is respectively Z _x, Z _y, Z _x', Z _y', shown in (5),

$E_s=\left[\begin{array}{c}{E}_x\\ E_y\\ E_{x^{\′}}\\ E_{y^{\′}}\end{array}\right]=\left[\begin{array}{c}{A}_{x}T\\ A_{y}T\\ A_x{\mathrm{\φ}}_{x}T\\ A_y{\mathrm{\φ}}_{y}T\end{array}\right]---\left(5\right)$

Step 5: structure φ _x, φ _ysimilar matrix F, H;

Formula (6) can be derived by formula (5):

E _x'＝E _xT ^-1φ _xT＝E _xT E _y'＝E _yT ^-1φ _yT＝E _yH (6)

Wherein, F=T ^-1φ _xt, H=T ^-1φ _yt, T are non-singular matrix, therefore F and φ _x, H and φ _yfor similar matrix, have identical eigenwert, and its eigenwert is twiddle factor φ _x, φ _ydiagonal element.

Step 6: least square method solves twiddle factor φ _x, φ _y, calculate direction of arrival

Twiddle factor φ is solved by least square method _x, φ _yshown in (7), just therefrom can draw the direction of arrival of signal.

$\begin{array}{cc}\hat{F}=E_x^{+}{E}_{x^{\′}}& \hat{H}=E_y^{+}{E}_y\end{array}----\left(7\right)$

for E _xpseudoinverse, for E _ypseudoinverse, Eigenvalues Decomposition is carried out to F and obtains obtain simultaneously the u of estimated value _k, Eigenvalues Decomposition is carried out to H and obtains obtain simultaneously the v of estimated value _k. can be estimated by formula (8):

Tool of the present invention has the following advantages: (1), compared with traditional coherent signal subspace algorithm, the present invention does not need to calculate spectral function, indirectly solves direction of arrival, reduce complexity without the need to search peak; (2) when number of sensors is determined, increase array aperture by virtual array, decrease equipment complexity and cost; (3) use the estimated value of linear operation direct solution sound source two dimension DOA, there is higher positioning precision.

Accompanying drawing explanation

Fig. 1 is L-type array structure schematic diagram;

Fig. 2 is the virtual array Z of translation along the y-axis direction _x';

Fig. 3 is the virtual array Z of translation along the x-axis direction _y'.

Fig. 4 is the actual direction of arrival of Fig. 4 sound source (/ °).

Fig. 5 is MUSIC algorithm estimated performance (SNR=20dB) in three information source situations

Fig. 6 is two-dimensional virtual array algorithm estimated performance (SNR=20dB) in three information source situations

Embodiment

Below in conjunction with drawings and Examples, the invention will be further described:

Step 1: structure L-type array, determines the signal model of array received;

As shown in Figure 1, receive acoustic pressure time-domain signal by 2M-1 array element in L-type array, this L-type array is the even linear array Z of M by array number in x-axis _xwith the even linear array Z that array number in y-axis is M _yform, wherein, 2M-1 is L-type array elements number, M be not less than 2 integer, d is array element distance.Hypothesis space has K information source to incide on array, and its 2-d direction finding is be respectively the elevation angle and the position angle of a kth signal source.

Suppose that the signal number incided on this array is K, then signal x-axis, y-axis received by M array element is respectively respectively as shown in the formula (1) and formula (2):

x(t)＝A _xs(t)+n(t) (1)

y(t)＝A _ys(t)+n(t) (2)

In formula, s (t) is signal source matrix, and n (t) is noise matrix, A _x, A _y∈ C ^{m × K}, be respectively the direction matrix in L-type array x-axis, y-axis, can be expressed as:

Step 2: obtain output signal matrix Z by ESPRIT constructing virtual array;

As shown in Figure 2 and Figure 3, by the submatrix Z of ESPRIT by L battle array _x, Z _ycarrying out virtual expansion is submatrix Z _x', Z _y', the motion immovability due to submatrix defines the rotational invariance of two submatrix signals, i.e. Z _x'submatrix signal be actual submatrix Z _xinput signal be multiplied by twiddle factor φ _xobtain, Z _y'submatrix signal be actual submatrix Z _yinput signal be multiplied by twiddle factor φ _yobtain, first obtained the output signal of virtual submatrix by formula (1) and formula (2), then four submatrixs are exported and merged, form output signal matrix z (t) of whole array as shown in the formula (3):

$z\left(t\right)={[x\left(t\right),x^{\′}\left(t\right)y^{\′}\left(t\right)]}^T=\stackrel{\‾}{A}s\left(t\right)+n\left(t\right)$

Its suppose that the direction of arrival of each information source is different, then column vector between Line independent,

And

Wherein, matrix φ _x, φ _yfor the diagonal matrix of K × K, its diagonal element is that signal is respectively at Z _x, Z _yphase delay on array between any array element, diag represents diagonal matrix, and the element namely except principal diagonal is the square formation of zero.

X-axis direction homogenous linear submatrix Z is comprised such as formula (3) z (t) _xoutput signal x (t) _,y-axis direction homogenous linear submatrix Z _youtput signal y (t), Z _xthe virtual submatrix Z that translation obtains _x'output signal x'(t), Z _ythe virtual submatrix Z that translation obtains _y'output signal y'(t), the noise that each array received arrives is identical, virtual submatrix Z _x', Z _y'the even linear array of to be all array number be M.

Step 3: obtain correlation matrix R from array output signal matrix Z _z,

Signal subspace and noise subspace can obtain, shown in (4) by the feature decomposition that array exports the covariance matrix of Z:

$R_z=E\left[z\left(t\right)z^H\left(t\right)\right]=\stackrel{\‾}{A}{R}_s{\stackrel{\‾}{A}}^H+{\mathrm{\σ}}^{2}I---\left(4\right)$

R in formula _sfor the autocorrelation matrix of signal, σ ²for noise variance, I is unit matrix, the E [.] in formula (4), (.) ^hbe expressed as mathematical expectation, conjugate transpose operation.

Step 4: by correlation matrix R _zdo feature decomposition, estimated signal number;

Array correlation matrix R _ztwo spaces can be divided into, namely the eigenwert characteristic of correspondence vector E of K _s=[s ₁, s ₂... s _k] composition signal subspace, the non-singular matrix T that there is a K × K meets and move invariant feature E due to array _scan be analyzed to 4 parts, E _x, E _y, E _x', E _y'∈ C ^{m × K}, corresponding subarray is respectively Z _x, Z _y, Z _x', Z _y', shown in (5),

$E_s=\left[\begin{array}{c}{E}_x\\ E_y\\ E_{x^{\′}}\\ E_{y^{\′}}\end{array}\right]=\left[\begin{array}{c}{A}_{x}T\\ A_{y}T\\ A_x{\mathrm{\φ}}_{x}T\\ A_y{\mathrm{\φ}}_{y}T\end{array}\right]---\left(5\right)$

Step 5: structure φ _x, φ _ysimilar matrix F, H;

Formula (6) can be derived by formula (5):

E _x'＝E _xT ^-1φ _xT＝E _xT E _y'＝E _yT ^-1φ _yT＝E _yH (6)

Wherein, F=T ^-1φ _xt, H=T ^-1φ _yt, T are non-singular matrix, therefore F and φ _x, H and φ _yfor similar matrix, have identical eigenwert, and its eigenwert is twiddle factor φ _x, φ _ydiagonal element.Step 6: least square method solves twiddle factor φ _x, φ _y, calculate direction of arrival

Twiddle factor φ is solved by least square method _x, φ _yshown in (7), just therefrom can draw the direction of arrival of signal.

$\begin{array}{cc}\hat{F}=E_x^{+}{E}_{x^{\′}}& \hat{H}=E_y^{+}{E}_y\end{array}----\left(7\right)$

for E _xpseudoinverse, for E _ypseudoinverse, Eigenvalues Decomposition is carried out to F and obtains obtain simultaneously the u of estimated value _k, Eigenvalues Decomposition is carried out to H and obtains obtain simultaneously the v of estimated value _k. can be estimated by formula (8):

Step 7: based on the Dynamic simulation result of the virtual array DOA estimation method of L-type battle array;

Simulated conditions is: the L-type array longitudinal axis of employing, transverse axis are all the even linear array of 8 omnidirectional's array element compositions, the aerial velocity of propagation of sound is c=340 m/s, and frequency of source is f=3000 Hz, and array element distance gets d=λ/2, i.e. d=5.7 cm, signal to noise ratio snr=20 dB.Contrast the DOA estimated performance based on MUSIC algorithm and virtual array algorithm, analogue simulation three sound sources, actual direction of arrival is (10 °, 15 °), (30 °, 25 °), (50 °, 35 °), as shown in Figure 4, in figure, stain is the real incident angle of sound source.Use MUSIC algorithm to estimate the direction of arrival of sound source by search peak-to-peak value, as shown in Figure 5, the direction of arrival of 1 sound source can only be estimated, DOA is (28 °, 24 °), and secondary lobe is more, positioning precision is fuzzyyer, accurately can not identify the accurate location of sound source.By the relation between sound source position and microphone, the ripple using the virtual array DOA estimation method based on L-type battle array directly to calculate three sound sources reaches estimated value and is (10.01 °, 15.07 °), (30 °, 25.01 °), (49.98 °, 34.99 °), as shown in Figure 6, in figure, stain is the DOA of estimation, compared with real direction of arrival, the DOA calculated changes between ± 1 °, error is less, and the DOA estimation method positioning precision based on L-type battle array virtual array algorithm is relatively high.

Can find out in result, compared with MUSIC algorithm, virtual array algorithm does not need to calculate spectral function, indirectly solves direction of arrival, reduce complexity without the need to search peak; When number of sensors is determined, increase array aperture by virtual array, decrease equipment complexity and cost; Use the estimated value of linear operation direct solution sound source two dimension DOA, there is higher positioning precision.

Claims (1)

1., based on a virtual array DOA estimation method for L-type array, it is characterized in that: the method comprises the following steps:

Step 1: structure L-type array, determines the signal model of array received;

Receive acoustic pressure time-domain signal by 2M-1 array element in L-type array, this L-type array is the even linear array Z of M by array number in x-axis _xwith the even linear array Z that array number in y-axis is M _yform, wherein, 2M-1 is L-type array elements number, M be not less than 2 integer, d is array element distance; Hypothesis space has K information source to incide on array, and its 2-d direction finding is be respectively the elevation angle and the position angle of a kth signal source;

Suppose that the signal number incided on this array is K, then signal x-axis, y-axis received by M array element is respectively respectively as shown in the formula (1) and formula (2):

x(t)＝A _xs(t)+n(t) (1)

y(t)＝A _ys(t)+n(t) (2)

In formula, s (t) is signal source matrix, and n (t) is noise matrix, A _x, A _y∈ C ^{m × K}, be respectively the direction matrix in L-type array x-axis, y-axis, can be expressed as:

Step 2: obtain output signal matrix Z by ESPRIT constructing virtual array;

By the submatrix Z of ESPRIT by L battle array _x, Z _ycarrying out virtual expansion is submatrix Z _x', Z _y', the motion immovability due to submatrix defines the rotational invariance of two submatrix signals, i.e. Z _x'submatrix signal be actual submatrix Z _xinput signal be multiplied by twiddle factor φ _xobtain, Z _y'submatrix signal be actual submatrix Z _yinput signal be multiplied by twiddle factor φ _yobtain, first obtained the output signal of virtual submatrix by formula (1) and formula (2), then four submatrixs are exported and merged, form output signal matrix z (t) of whole array as shown in the formula (3):

$z\left(t\right)=[x\left(t\right),y\left(t\right),x^{\′}\left(t\right),y^{\′}\left(t\right){]}^T=\stackrel{\‾}{A}s\left(t\right)+n\left(t\right)---\left(3\right)$

Wherein suppose that the direction of arrival of each information source is different, then column vector between Line independent, and

Wherein, matrix φ _x,φ _yfor the diagonal matrix of K × K, its diagonal element is that signal is respectively at Z _x, Z _yphase delay on array between any array element, diag represents diagonal matrix, and the element namely except principal diagonal is the square formation of zero;

X-axis direction homogenous linear submatrix Z is comprised such as formula (3) z (t) _xoutput signal x (t) _,y-axis direction homogenous linear submatrix Z _youtput signal y (t), Z _xthe virtual submatrix Z that translation obtains _x'output signal x'(t), Z _ythe virtual submatrix Z that translation obtains _y'output signal y'(t), the noise that each array received arrives is identical, virtual submatrix Z _x',z _y'the even linear array of to be all array number be M;

Step 3: obtain correlation matrix R from array output signal matrix Z _z,

Signal subspace and noise subspace can obtain, shown in (4) by the feature decomposition that array exports the covariance matrix of Z:

$R_z=E\left[z\left(t\right)z^H\left(t\right)\right]=\stackrel{\‾}{A}{R}_s{\stackrel{\‾}{A}}^H+{\mathrm{\σ}}^{2}I---\left(4\right)$

R in formula _sfor the autocorrelation matrix of signal, σ ²for noise variance, I is unit matrix, the E [.] in formula (4), (.) ^hbe expressed as mathematical expectation, conjugate transpose operation;

Step 4: by correlation matrix R _zdo feature decomposition, estimated signal number;

Array correlation matrix R _ztwo spaces can be divided into, namely the eigenwert characteristic of correspondence vector E of K _s=[s ₁, s ₂... s _k] composition signal subspace, the non-singular matrix T that there is a K × K meets , and move invariant feature E due to array _scan be analyzed to 4 parts, E _x, E _y, E _x', E _y'∈ C ^{m × K}, corresponding subarray is respectively Z _x, Z _y, Z _x', Z _y', shown in (5),

$E_s=\left[\begin{array}{c}{E}_x\\ E_y\\ E_{x^{\′}}\\ E_{y^{\′}}\end{array}\right]=\left[\begin{array}{c}{A}_{x}T\\ A_{y}T\\ A_x{\mathrm{\φ}}_{x}T\\ A_y{\mathrm{\φ}}_{y}T\end{array}\right]---\left(5\right)$

Step 5: structure φ _x, φ _ysimilar matrix F, H;

Formula (6) can be derived by formula (5):

E _x'＝E _xT ^-1φ _xT＝E _xT E _y'＝E _yT ^-1φ _yT＝E _yH (6)

Wherein, F=T ^-1φ _xt, H=T ^-1φ _yt, T are non-singular matrix, therefore F and φ _x, H and φ _yfor similar matrix, have identical eigenwert, and its eigenwert is twiddle factor φ _x, φ _ydiagonal element;

Step 6: least square method solves twiddle factor φ _x, φ _y, calculate direction of arrival

Twiddle factor φ is solved by least square method _x, φ _yshown in (7), just therefrom can draw the direction of arrival of signal; $\hat{F}=E_x^{+}{E}_{x^{\′}},\hat{H}=E_y^{+}{E}_{y^{\′}}---\left(7\right)$

for E _xpseudoinverse _, for E _ypseudoinverse, Eigenvalues Decomposition is carried out to F and obtains obtain simultaneously sin θ _kthe u of estimated value _k, Eigenvalues Decomposition is carried out to H and obtains obtain simultaneously sin θ _kthe v of estimated value _k; can be estimated by formula (8):