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

Abstract

The invention discloses a kind of virtual array DOA estimation method based on L-type array, comprise the following steps：(1) based on motion immovability matter, by the submatrix Z of L-type array_x,Z_yTranslation obtains virtual array Z_x',Z_y', the rotational invariance of two subsignals is formd due to the motion immovability of submatrix, the signal of virtual submatrix is equivalent to L-type submatrix Z_x,Z_yInput signal is multiplied by twiddle factor and obtains respectively；(2) output of 4 submatrixs is merged, is constituted output signal Z (t) of virtual array；(3) signal subspace and noise subspace can be described with the feature decomposition of the covariance matrix of array output, and cross correlation process is carried out to array output signal Z (t), obtain R_zz, carry out Eigenvalues Decomposition and obtain signal subspace；(4) twiddle factor is solved by linear operation, is that can obtain signal wave up to direction by its diagonal element.The present invention need not calculate spectral function, and direction of arrival is solved indirectly without search peak, reduce complexity；Reduce equipment complexity and cost；With positioning precision higher.

Description

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

Technical field

L-type array is based on the invention belongs to Array Signal Processing arrival direction estimation technical field, more particularly to one kind Virtual array DOA estimation method.

Background technology

Target Bearing Estimation (DOA) estimates it is the letter such as sonar, radar, radio communication, medical imaging, microphone array column processing One important branch of number process field, the basic problem that DOA is solved is to determine multiple senses in a certain space field simultaneously The locus of the echo signal of interest, i.e. each echo signal reach the deflection of sensor array.Based on conventional wave beam shape Direction-finding method into scanning has intrinsic restriction, is limited by " Rayleigh criterion ", and the resolution ratio of estimation depends on array length, Only when the separation angle between two information sources in space more than array aperture it is reciprocal when, can just be resolved.In order to improve basic matrix Resolution ratio, when signal wavelength lambda must be, the general aperture length for only increasing basic matrix increases array element m or increase battle array spacing d, But increasing element number of array can improve equipment complexity and cost, and increase array element spacing will cause secondary maximum, for actual feelings again The limitation of condition, the yardstick of array is also impossible to be made very big, therefore relies solely on increase array aperture and reaches raising resolution ratio Way, be difficult to be applicable in Practical Project.In order to overcome this limitation, a kind of virtual array DOA based on L-type battle array is devised Method of estimation, its basic thought is the virtual expansion for obtaining array by optimized algorithm under the array case of limited dimensions, from And improve resolution ratio.

The content of the invention

The present invention is in view of the shortcomings of the prior art, it is proposed that a kind of virtual array DOA estimation method based on L-type battle array.

Technical scheme is as follows：

A kind of virtual array DOA estimation method based on L-type array, wherein, comprise the following steps：

Step 1：Construction L-type array, determines the signal model of array received；

Acoustic pressure time-domain signal is received by 2M-1 array element in L-type array, this L-type array is the uniform of M by array number in x-axis Linear array Z_xIt is the even linear array Z of M with array number in y-axis_yConstitute, wherein, 2M-1 is L-type array elements number, and M is not less than 2 Integer, d is array element spacing.Assuming that space has K information source to incide on array, its 2-d direction finding is Respectively k-th elevation angle and azimuth of signal source.

Assuming that it is K to incide the signal number on this array, then the signal point for being received by M array element respectively in x-axis, y-axis It is not such as following formula (1) and formula (2)：

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

Y (t)=A_ys(t)+n(t) (2)

S (t) is signal source matrix in formula, and n (t) is noise matrix, A_x, A_y∈C^M×K, respectively in L-type array x-axis, y-axis Direction matrix, be represented by：

Step 2：Output signal matrix Z is obtained by ESPRIT construction virtual array；

By ESPRIT by L gusts of submatrix Z_x,Z_yIt is submatrix Z to carry out virtual expansion_x',Z_y', due to the shifting of submatrix Consistency forms two rotational invariances of submatrix signal, 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, by formula (1) output signal of virtual submatrix is first obtained with formula (2), is then merged four submatrix outputs, constitute whole array Output signal matrix z (t) such as following formula (3)：

WhereinAssuming that the direction of arrival of each information source is different, thenRow Line independent between vector,

And

Wherein, matrix φ_x, φ_yIt is the diagonal matrix of K × K, its diagonal element is signal respectively in Z_x, Z_yIt is any on array Phase delay between array element, diag represents diagonal matrix, i.e., the element in addition to leading diagonal is zero square formation.

As formula (3) z (t) includes x-axis direction homogenous linear submatrix Z_xOutput signal x (t)_、Y-axis direction homogenous linear submatrix Z_yOutput signal y (t), Z_xThe virtual submatrix Z that translation is obtained_x'Output signal x'(t), Z_yThe virtual submatrix Z that translation is obtained_y' Output signal y'(t), the noise that each array received is arrived is identical, virtual submatrix Z_x'、Z_y'All for array number is the even linear array of M.

Step 3：Correlation matrix R is obtained from array output signal matrix Z_z,

The feature decomposition of the covariance matrix that signal subspace and noise subspace can export Z with array is obtained, such as formula (4) It is shown：

R in formula_sIt is the autocorrelation matrix of signal, σ²It is noise variance, I is unit matrix, the E [] in formula (4), ()^HPoint Mathematic expectaion, conjugate transposition computing are not expressed as.

Step 4：By correlation matrix R_zFeature decomposition is done, signal number is estimated；

Array correlation matrix R_zTwo spaces can be divided into, i.e.,The characteristic value of K Corresponding characteristic vector E_s=[s₁,s₂,...s_k] composition signal subspace, there is a non-singular matrix T of K × K and meetAnd due to the shifting invariant feature E of array_sCan be analyzed to 4 parts, E_x,E_y,E_x',E_y'∈C^M×K, corresponding submatrix Row are respectively Z_x, Z_y, Z_x', Z_y', as shown in formula (5),

Step 5：Construction φ_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 φ_yIt is similar matrix, possesses Identical characteristic value, and its characteristic value is twiddle factor φ_x, φ_yDiagonal element.

Step 6：Least square method solves twiddle factor φ_x, φ_y, calculate direction of arrival

Twiddle factor φ is solved with least square method_x, φ_yAs shown in formula (7), the direction of arrival of signal just can be therefrom drawn.

It is E_xPseudoinverse,It is E_yPseudoinverse, Eigenvalues Decomposition is carried out to F and is obtained Obtain simultaneouslyEstimate u_k, Eigenvalues Decomposition is carried out to H and is obtainedObtain simultaneouslyEstimate v_k。Can be by Formula (8) is estimated：

The invention has the advantages that：(1) compared with traditional coherent signal subspace algorithm, the present invention need not calculate spectrum Function, direction of arrival is solved without search peak indirectly, reduces complexity；(2) in the case where number of sensors determines, lead to Virtual array increase array aperture is crossed, equipment complexity and cost is reduced；(3) with linear operation direct solution sound source two dimension The estimate of DOA, with positioning precision higher.

Brief description of the drawings

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

Fig. 2 is the virtual array Z for translating along the y-axis direction_x'；

Fig. 3 is the virtual array Z for translating along the x-axis direction_y'。

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

Fig. 5 is MUSIC algorithms estimation performance (SNR=20dB) in the case of three information sources

Fig. 6 is two-dimensional virtual array algorithm estimation performance (SNR=20dB) in the case of three information sources

Specific embodiment

The invention will be further described with reference to the accompanying drawings and examples：

Step 1：Construction L-type array, determines the signal model of array received；

As shown in figure 1, receiving acoustic pressure time-domain signal by 2M-1 array element in L-type array, this L-type array is by array element in x-axis Number is the even linear array Z of M_xIt is the even linear array Z of M with array number in y-axis_yConstitute, wherein, 2M-1 is L-type array elements number, M It is the integer not less than 2, d is array element spacing.Assuming that space has K information source to incide on array, its 2-d direction finding is Respectively k-th elevation angle and azimuth of signal source.

Assuming that it is K to incide the signal number on this array, then the signal point for being received by M array element respectively in x-axis, y-axis It is not such as following formula (1) and formula (2)：

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

Y (t)=A_ys(t)+n(t) (2)

S (t) is signal source matrix in formula, and n (t) is noise matrix, A_x, A_y∈C^M×K, respectively in L-type array x-axis, y-axis Direction matrix, be represented by：

Step 2：Output signal matrix Z is obtained by ESPRIT construction virtual array；

As shown in Figure 2 and Figure 3, by ESPRIT by L gusts of submatrix Z_x,Z_yIt is submatrix Z to carry out virtual expansion_x', Z_y', because the motion immovability of submatrix forms two rotational invariances of submatrix signal, 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, the output signal of virtual submatrix is first obtained by formula (1) and formula (2), then closed four submatrix outputs And, constitute output signal matrix z (t) such as following formula (3) of whole array：

ItsAssuming that the direction of arrival of each information source is different, thenRow arrow Line independent between amount,

And

Wherein, matrix φ_x, φ_yIt is the diagonal matrix of K × K, its diagonal element is signal respectively in Z_x, Z_yIt is any on array Phase delay between array element, diag represents diagonal matrix, i.e., the element in addition to leading diagonal is zero square formation.

As formula (3) z (t) includes x-axis direction homogenous linear submatrix Z_xOutput signal x (t)_、Y-axis direction homogenous linear submatrix Z_yOutput signal y (t), Z_xThe virtual submatrix Z that translation is obtained_x'Output signal x'(t), Z_yThe virtual submatrix Z that translation is obtained_y' Output signal y'(t), the noise that each array received is arrived is identical, virtual submatrix Z_x'、Z_y'All for array number is the even linear array of M.

Step 3：Correlation matrix R is obtained from array output signal matrix Z_z,

The feature decomposition of the covariance matrix that signal subspace and noise subspace can export Z with array is obtained, such as formula (4) It is shown：

R in formula_sIt is the autocorrelation matrix of signal, σ²It is noise variance, I is unit matrix, the E [] in formula (4), ()^HPoint Mathematic expectaion, conjugate transposition computing are not expressed as.

Step 4：By correlation matrix R_zFeature decomposition is done, signal number is estimated；

Array correlation matrix R_zTwo spaces can be divided into, i.e.,The characteristic value of K Corresponding characteristic vector E_s=[s₁,s₂,...s_k] composition signal subspace, there is a non-singular matrix T of K × K and meetAnd due to the shifting invariant feature E of array_sCan be analyzed to 4 parts, E_x,E_y,E_x',E_y'∈C^M×K, corresponding submatrix Row are respectively Z_x, Z_y, Z_x', Z_y', as shown in formula (5),

Step 5：Construction φ_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 φ_yIt is similar matrix, possesses Identical characteristic value, and its characteristic value is twiddle factor φ_x, φ_yDiagonal element.Step 6：Least square method solve rotation because Sub- φ_x, φ_y, calculate direction of arrival

Twiddle factor φ is solved with least square method_x, φ_yAs shown in formula (7), the direction of arrival of signal just can be therefrom drawn.

It is E_xPseudoinverse,It is E_yPseudoinverse, Eigenvalues Decomposition is carried out to F and is obtained Obtain simultaneouslyEstimate u_k, Eigenvalues Decomposition is carried out to H and is obtainedObtain simultaneouslyEstimate v_k。Can be by Formula (8) is estimated：

Step 7：The operation simulation result of the virtual array DOA estimation method based on L-type battle array；

Simulated conditions are：The L-type array longitudinal axis, the transverse axis for using all are 8 even linear arrays of omnidirectional's array element composition, and sound exists Spread speed in air is c=340 m/s, and frequency of source is f=3000 Hz, and array element spacing takes d=λ/2, i.e. d=5.7 Cm, signal to noise ratio snr=20 dB.Contrast estimates performance, analogue simulation three based on MUSIC algorithms with the DOA of virtual array algorithm Sound source, actual direction of arrival is (10 °, 15 °), (30 °, 25 °), (50 °, 35 °), as shown in figure 4, in figure stain to be sound source true Real incident angle.The direction of arrival of sound source is estimated by searching for peak-to-peak value with MUSIC algorithms, as shown in figure 5, can only estimate Go out 1 direction of arrival of sound source, DOA is (28 °, 24 °), and secondary lobe is more, and positioning precision is relatively obscured, it is impossible to sound source is recognized accurately Accurate location.By the relation between sound source position and microphone, with the virtual array DOA estimation method based on L-type battle array Directly calculate three ripples of sound source and reach estimate for (10.01 °, 15.07 °), (30 °, 25.01 °), (49.98 °, 34.99 °), As shown in fig. 6, stain is the DOA of estimation in figure, compared with real direction of arrival, the DOA being calculated is between ± 1 ° Change, error is smaller, and the DOA estimation method positioning precision based on L-type battle array virtual array algorithm is of a relatively high.

As can be seen that compared with MUSIC algorithms, virtual array algorithm need not calculate spectral function in result, without search Peak value solves direction of arrival indirectly, reduces complexity；In the case where number of sensors determines, battle array is increased by virtual array Row aperture, reduces equipment complexity and cost；With the estimate of linear operation direct solution sound source two dimension DOA, with compared with Positioning precision high.

Claims (1)

1. a kind of virtual array DOA estimation method based on L-type array, it is characterised in that：The method is comprised the following steps：

Step 1：Construction L-type array, determines the signal model of array received；

Acoustic pressure time-domain signal is received by 2M-1 array element in L-type array, this L-type array is the even linear array of M by array number in x-axis Z_xIt is the even linear array Z of M with array number in y-axis_yConstitute, wherein, 2M-1 is L-type array elements number, and M is whole not less than 2 Number, d is array element spacing；Assuming that space has K information source to incide on array, its 2-d direction finding isθ_k,Respectively k-th elevation angle and azimuth of signal source；

Assuming that it is K to incide the signal number on this array, then the signal for being received by M array element respectively in x-axis, y-axis is respectively Such as following formula (1) and formula (2)：

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

Y (t)=A_ys(t)+n(t) (2)

S (t) is signal source matrix in formula, and n (t) is noise matrix, A_x, A_y∈C^M×K, the side respectively in L-type array x-axis, y-axis To matrix, it is represented by：

Step 2：Output signal matrix Z is obtained by ESPRIT construction virtual array；In the feelings that number of sensors determines Under condition, array aperture is increased by virtual array, reduce equipment complexity and cost；

By ESPRIT by L gusts of submatrix Z_x,Z_yIt is submatrix Z to carry out virtual expansion_x',Z_y', because the shifting of submatrix is constant Property forms two rotational invariances of submatrix signal, i.e. Z_x' submatrix signal be actual submatrix Z_xInput signal be multiplied by rotation Factor φ_xObtain, Z_y'Submatrix signal be actual submatrix Z_yInput signal be multiplied by twiddle factor φ_yObtain, by formula (1) With the output signal that formula (2) first obtains virtual submatrix, then four submatrix outputs are merged, constituted the defeated of whole array Go out signal matrix z (t) such as following 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)$

WhereinAssuming that the direction of arrival of each information source is different, thenColumn vector it Between Line independent,

And

Wherein, matrix φ_x, φ_yIt is the diagonal matrix of K × K, its diagonal element is signal respectively in Z_x, Z_yAny array element on array Between phase delay, diag represents diagonal matrix, i.e., the element in addition to leading diagonal is zero square formation；

As formula (3) z (t) includes x-axis direction homogenous linear submatrix Z_xOutput signal x (t), y-axis direction homogenous linear submatrix Z_y's Output signal y (t), Z_xThe virtual submatrix Z that translation is obtained_x'Output signal x'(t), Z_yThe virtual submatrix Z that translation is obtained_y'It is defeated Go out signal y'(t), the noise that each array received is arrived is identical, virtual submatrix Z_x'、Z_y'All for array number is the even linear array of M；

Step 3：Correlation matrix R is obtained from array output signal matrix Z_z,

The feature decomposition of the covariance matrix that signal subspace and noise subspace can export Z with array is obtained, such as formula (4) institute Show：

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

R in formula_sIt is the autocorrelation matrix of signal, σ²It is noise variance, I is unit matrix, the E [] in formula (4), ()^HDifference table It is shown as mathematic expectaion, conjugate transposition computing；

Step 4：By correlation matrix R_zFeature decomposition is done, signal number is estimated；

Array correlation matrix R_zTwo spaces can be divided into, i.e.,The characteristic value of K is corresponding Characteristic vector E_s=[s₁,s₂,...s_k] composition signal subspace, there is a non-singular matrix T of K × K and meetAnd And due to the shifting invariant feature E of 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', as shown in formula (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：Construction φ_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 φ_yIt is similar matrix, possesses identical Characteristic value, and its characteristic value be twiddle factor φ_x, φ_yDiagonal element；

Step 6：Least square method solves twiddle factor φ_x, φ_y, calculate direction of arrival

Twiddle factor φ is solved with least square method_x, φ_yAs shown in formula (7), the direction of arrival of signal just can be therefrom drawn；

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

It is E_xPseudoinverse,It is E_yPseudoinverse, Eigenvalues Decomposition is carried out to F and is obtained Obtain simultaneouslysinθ_kEstimate u_k, Eigenvalues Decomposition is carried out to H and is obtained Obtain simultaneouslysinθ_kEstimate v_k；θ_k,Can be estimated by formula (8)：