2D-DOA estimation method and system for acoustic vector sensor of any array structure

Abstract

The invention provides an acoustic vector sensor 2D-DOA estimation method and system of any array structure, and relates to the array sensor direction finding technology. According to the method, the covariance matrix is obtained, characteristic decomposition is performed on the covariance matrix, the signal subspace is solved, rough estimation is performed on the two-dimensional DOA in the vector information according to the signal subspace by adopting the ESPRIT algorithm, and then phase compensation is performed on spatial ambiguity by utilizing the rough estimation result. The 2D-DOA accurate estimation of the acoustic vector sensor of any array structure is realized; the spatial layout of any array is realized by using the vector sensor, and the array element spacing of the sensor is not limited by half wavelength any more; according to the method, the information of a vector part is utilized, the information of a space phase is also utilized, the estimation precision is remarkably improved, spectral peak search and additional two-dimensional angle pairing calculation do not need to be carried out on a covariance matrix of a received signal, and the calculation complexity is remarkably reduced.

Description

Technical Field

[0001] The invention relates to an array sensor direction finding technology, and in particular to a 2D-DOA estimation method and system for an acoustic vector sensor with an arbitrary array structure. Background Art

[0002] MIMO (Multiple Input Multiple Output) radar is a radar system that draws on MIMO technology in the field of communications. It was first proposed by Bliss and Forsythe in 2003. The concept of MIMO radar marks a major development in radar technology and has aroused great interest from many scholars and research institutions at home and abroad. Since then, radar technology has entered a new stage of development. Unlike traditional single-input single-output (SISO) radar or phased array radar, MIMO radar transmits mutually orthogonal waveforms through multiple active antennas, so that the radar system has spatial and temporal diversity at the same time, and creates a virtual aperture, which improves the angular resolution and enables the radar system to more accurately determine the location of the target. Due to the use of multiple independent signal paths, the anti-interference ability of MIMO radar is also greatly increased. With the rapid development of MIMO radar, it has shown its unique advantages in the fields of air surveillance, ground surveillance, maritime surveillance, missile guidance, etc.

[0003] In the field of multiple-input multiple-output (MIMO) radar, the direction finding task has always been a hot research topic. At present, many excellent angle estimation methods have emerged, such as the multiple signal classification (MUSIC) method, the rotational invariance technique (ESPRIT), the propagator method (PM) method, and the tensor-based (PARAFAC) method. Among them, the MUSIC algorithm can handle the high-resolution angle estimation problem under arbitrary array geometry well, but it requires high-complexity peak search and high computational cost. The ESPRIT algorithm utilizes the invariance of MIMO radar angle estimation and does not require spectral peak search, but it needs to construct the covariance matrix of the received signal and its eigendecomposition; while the PM algorithm uses the internal relationship of the received data to calculate the propagation matrix, and then constructs the signal subspace or noise subspace, without the need for eigenvalue decomposition. Therefore, the PM algorithm has lower computational complexity than the ESPRIT algorithm, but the PM algorithm has poor estimation accuracy under low signal-to-noise ratio conditions and may even fail. The PARAFAC algorithm is also used in MIMO radar systems for multi-target localization based on trilinear decomposition. It can be applied to arrays of arbitrary geometric shapes, but iterations lead to high complexity and high computational cost.

[0004] At present, the direction-finding tasks of co-located MIMO radars are mostly concentrated on one-dimensional direction-finding using transmit (Tx) and receive (Rx) arrays configured with scalar sensors. If two-dimensional (2D) direction-finding is to be achieved, nonlinear sensor arrays must be used, such as L-shaped arrays, circular arrays, and uniform rectangular arrays (URA). However, although traditional sensor arrays can perform angle estimation, their spatial position layout is very fixed, and it is difficult to achieve such a fixed layout in some specific scenarios. On the other hand, since the direction-finding algorithm of co-located MIMO radars relies on the spatial phase characteristics of the scalar sensor array, the array element spacing between sensors must be less than (or equal to) half a wavelength, otherwise it will cause angle ambiguity.

[0005] In the prior art, vector sensor arrays can solve the above problems well, which can realize arbitrary array layout and make the sensor element spacing not limited by half wavelength. However, the current research on vector sensor arrays mainly focuses on the vector component part, ignoring the original spatial phase information, and the vector components are likely to be coupled with each other, which leads to poor accuracy of the predicted estimation results. Summary of the invention

[0006] The main purpose of the present invention is to provide a 2D-DOA estimation method and system for an acoustic vector sensor with an arbitrary array structure, so as to solve the technical problem that the prior art only focuses on the information of the vector component part and ignores the original spatial phase information, resulting in poor accuracy of the predicted estimation results.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is: a 2D-DOA estimation method of an acoustic vector sensor with an arbitrary array structure, comprising the following steps:

[0008] S1: The receiving end obtains the covariance matrix, performs eigendecomposition on the covariance matrix, and solves the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors;

[0009] S2: Using the signal subspace and the ESPRIT algorithm, a rough estimate of 2D-DOA is obtained;

[0010] S3: Simultaneously use the signal subspace to perform matrix transformation and solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, obtain the ambiguous spatial phase difference matrix;

[0011] S4: Perform phase compensation on the ambiguous spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA;

[0012] Here, 2D-DOA represents two-dimensional angle of arrival.

[0013] In the preferred embodiment, the specific steps of step S1 are as follows:

[0014] S101: Based on the fact that the array positions of the transmitting array element and the receiving array element are set arbitrarily, and the transmitting array element has M transmitting sensors and the receiving array element has N speed receiving sensors, the coordinates of the first row of the transmitting array and the receiving array are used as reference array elements, and the arrival time differences relative to the mth transmitting sensor and the nth speed receiving sensor are respectively:

[0015] d t,m,k =x tm cos(φ tk )sin(θ tk )+y tm sin(φ tk )sin(θ tk )+z tm cos(θ tk ) (1);

[0016] d r,n,k =x rn cos(φ rk )sin(θ rk )+yrn sin(φ rk )sin(θ rk )+z rn cos(θ rk ) (2);

[0017] In the formula, (x tm ,y tm ,z tm ) is the position of the mth transmitting sensor, (x rn ,y rn ,z rn ) is the position of the nth velocity receiving sensor, θ tk ,φ tk are the azimuth and elevation of the 2D-DOD of the kth target, θ rk ,φ rk are the azimuth and elevation of the 2D-DOA of the kth target, respectively; 2D-DOD represents the two-dimensional wave departure angle;

[0018] S102: Under the narrowband assumption, the transmitting sensor array spatial response vector a t (θ tk ,φ tk ) and the receiving sensor array spatial response vector a r (θ rk ,φ rk ) are:

[0019]

[0020] Where, d t,m,k is the arrival time difference of the mth transmitting sensor relative to the reference sensor array element, d r,n,k is the arrival time difference of the nth velocity receiving sensor relative to the reference sensor array element, λ is the wavelength of the transmitted signal;

[0021] S103: The output formula of the matching filter of the receiving end speed receiving sensor is:

[0022]

[0023] or:

[0024] Y(t)=[A t ⊙A r ⊙B r ]s(t)+n(t) (6);

[0025] In the formula, A t =[a t (θ t1 ,φ t1 ),a t (θt2 ,φ t2 ),...,a t (θ tK ,φ tK )]∈C M×K is the transmit array spatial response matrix, A r =[a r (θ r1 ,φ r1 ),a r (θ r2 ,φ r2 ),...,a r (θ rK ,φ rK )]∈C N×K is the spatial response matrix of the receiving array, s(t)=[s 1 (t),s 2 (t),…,s K (t)] T ∈C K×L , f k is the Doppler shift, β k is the amplitude, n(t)∈C MN×L is the received Gaussian white noise matrix, ⊙ represents the Khatri-Rao product, represents the Kroc inner product, L is the number of snapshots; B r =[b r1 ,b r2 ,b r3 ] T ∈C 3×K is the polarization response matrix, and:

[0026] b r1 =[sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ] T ,

[0027] b r2 =[sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,sinθ rK sinφ rK ] T ,

[0028] b r3 =[cosθ r1 ,cosθr2 ,...,cosθ rK ] T ;

[0029] S104: Define matrix A according to step S103:

[0030] A=A t ⊙A r ⊙B r (7);

[0031] The obtained signal model is:

[0032] Y = AS + N (8);

[0033] Calculate the covariance matrix R of the received signal Y for:

[0034] R Y =YY H / L (9);

[0035] In the formula, the superscript H is the conjugate transpose, and L is the number of snapshots;

[0036] S105: Solve the signal subspace, and calculate the received signal covariance matrix R in equation (9): Y Perform eigendecomposition to obtain the signal subspace. The eigendecomposition formula is:

[0037]

[0038] In the formula, is a K×K diagonal matrix whose diagonal elements contain K large eigenvalues, is a matrix composed of eigenvectors corresponding to K larger eigenvalues, The diagonal elements of contain 3MN-K smaller eigenvalues, is the matrix of eigenvectors corresponding to the remaining smaller eigenvalues, and can be regarded as a signal subspace, It can be regarded as a noise subspace;

[0039] Signal subspace The same as the subspace spanned by the spatial response vector matrix A, that is:

[0040]

[0041] In the formula, the matrix T∈C K×K , is a non-singular full-rank matrix.

[0042] In the preferred embodiment, the specific steps of step S2 are as follows:

[0043] S201: Set the selectivity matrix, the formula is:

[0044] J 1 =[I M ⊙I N ⊙I 3，1 ] (12a);

[0045] J 2 =[I M ⊙I N ⊙I 3,2 ] (12b);

[0046] J 3 =[I M ⊙I N ⊙I 3,3 ] (12c);

[0047] In the formula, 1 M is the M×M identity matrix, Ι N is the N×N identity matrix, Ι 3,k is the kth row of the 3×3 identity matrix;

[0048] S202: The selectivity matrix is ​​multiplied by the left side of the signal subspace matrix, and let Right now:

[0049]

[0050]

[0051] In the formula, B r,k Represented as matrix B r The kth row of , where:

[0052] Φ x =diag(sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ),

[0053] Φ y =diag(sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,sinθ rK sinφ rK ),

[0054] Φ z =diag(cosθ r1 ,cosθr2 ,...,cosθ rK );

[0055] S203: According to formula (13), the following matrix can be constructed:

[0056]

[0057] In the formula,

[0058] At the same time: let Φ xy =Φ z Φ xyz ,but:

[0059]

[0060] Replace A in formula (13c) 0 Substituting into formula (15), we can get:

[0061]

[0062] In the formula,

[0063] From formula (16), we can get: in represents the generalized inverse of a matrix, [·] -1 It means to find the inverse of the matrix;

[0064] S204: Perform eigendecomposition on the output matrix of step S203 to obtain a rough estimate of 2D-DOA. Specifically, Perform eigendecomposition to obtain eigenvalues and the eigenvector Right now:

[0065]

[0066] Then the rough estimate of the azimuth and elevation angles of 2D-DOA can be expressed as:

[0067]

[0068] In the formula, abs() means taking the amplitude, and angle() means taking the phase.

[0069] In the preferred solution, to solve the spatial ambiguity, step S3 specifically includes:

[0070] S301: Define the following spatial selectivity matrix, expressed as:

[0071] P n =[I M ⊙I N，n ⊙I3 ], n=1, ..., N (21);

[0072] In the formula, 1 N,n is the nth row of the N×N identity matrix;

[0073] S302: Multiply the selectivity matrix by the left side of the signal subspace matrix, that is:

[0074]

[0075] In the formula, is the spatial response vector A of the receiving array r The first line of is the spatial response vector A of the receiving array r The second line of

[0076] S303: and A r,2 =A r,1 Φ r1 , substituting into formula (22) we get:

[0077]

[0078] at this time:

[0079]

[0080] In the formula, u rk =cosφ rk sinθ rk , v rk = sinφ rk sinθ rk , w rk =cosφ rk , diag() means converting the row matrix into a diagonal matrix, let x r1 -x r2 =m 12 ,y r1 -y r2 =n 12 , z r1 -z r2 = l 12 ,but:

[0081]

[0082] S304: From formula (23), we can get Then substitute T in formula (18) to obtain:

[0083]

[0084] Similarly, using the selectivity matrix P 3Multiplying both sides of equation (11) at the same time, we can get And A r,3 =A r,2 Φ r2 ,Right now:

[0085]

[0086] at this time m 23 =x r2 -x r3 , n 23 =y r2 -y r3 , l 23 =z r2 -z r3 , similarly, substituting T in formula (18) into formula (27) yields:

[0087]

[0088] Same as above, using the selectivity matrix P 4,...,N ,Will The simulation formula (27)(28)(29) is obtained;

[0089] S305: Take The diagonal elements of Then put these N-1 row vectors into the initial phase matrix GG∈C N-1×K You can get:

[0090]

[0091] The initial phase matrix GG is phased to obtain the initial phase difference matrix GGp:

[0092]

[0093] That is: GGp=EP h (31);

[0094] in:

[0095]

[0096] In the preferred embodiment, step S4 specifically includes:

[0097] S401: First, using the rough estimation results in formulas (19) and (20), that is, obtaining the rough estimation results of the azimuth and elevation angles of the 2D-DOA in step S3, construct the following unambiguous and low-precision Poynting vector matrix: The expression is:

[0098]

[0099] S402: Use the unambiguous low-precision Poynting vector matrix By using formula (31), the unambiguous and low-precision phase difference matrix is ​​obtained: Through this unambiguous low-precision phase difference matrix Get blur Solve it, that is:

[0100]

[0101] S403: Blur By compensating to the periodically blurred phase difference matrix GGp through equation (34), we can obtain the unambiguous and high-precision phase difference matrix GGpt, namely:

[0102]

[0103] S404: Using the least squares principle, estimate the unambiguous and high-precision Poynting vector matrix Right now:

[0104]

[0105] Thus, the azimuth and elevation angles of 2D-DOA are obtained, that is, the azimuth and elevation angle estimation formulas of 2D-DOA are:

[0106]

[0107] The preferred solution also includes: using the root mean square error for evaluation, and the calculation formula is:

[0108]

[0109] In the formula, is θ rk The estimate, is φ rk Estimates.

[0110] A 2D-DOA estimation system for an acoustic vector sensor of an arbitrary array structure, comprising:

[0111] The signal subspace module is used for obtaining the covariance matrix at the receiving end, and performing eigendecomposition on the covariance matrix to solve the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors;

[0112] A rough estimation module is used to obtain a rough estimate of 2D-DOA using the signal subspace and the ESPRIT algorithm;

[0113] The spatial ambiguity module is used to simultaneously utilize the signal subspace and solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, to obtain the ambiguous spatial phase difference matrix;

[0114] An accurate estimation module is used to perform phase compensation on the spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA;

[0115] Here, 2D-DOA represents two-dimensional angle of arrival.

[0116] The present invention provides a 2D-DOA estimation method and system for an acoustic vector sensor of an arbitrary array structure. The method obtains a covariance matrix, performs eigendecomposition on the covariance matrix, solves a signal subspace, and then roughly estimates the two-dimensional DOA in the vector information based on the signal subspace and using an ESPRIT algorithm. The rough estimation result is then used to perform phase compensation for spatial ambiguity, thereby achieving accurate estimation of the 2D-DOA of an acoustic vector sensor of an arbitrary array structure. A vector sensor is used to realize the spatial layout of an arbitrary array, so that the array element spacing of the sensor is no longer limited by a half-wavelength. Not only the information of the vector part but also the information of the spatial phase is used, thereby significantly improving the estimation accuracy. There is no need to perform a spectrum peak search and additional two-dimensional angle pairing calculation on the covariance matrix of the received signal, thereby significantly reducing the computational complexity. BRIEF DESCRIPTION OF THE DRAWINGS

[0117] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0118] Figure 1 is a flow chart of the estimation method of the present invention;

[0119] Figure 2 It is a schematic diagram of a single-base MIMO radar system applicable to the present invention;

[0120] Figure 3 It is a Monte Carlo test simulation scatter plot of the present invention;

[0121] Figure 4 The present invention and Figure 3 The ESPRIT algorithm simulation scatter plot is used for comparison;

[0122] Figure 5 The present invention compares the RMSE under different signal-to-noise ratios SNR r picture;

[0123] Figure 6 is the RMSE of the present invention comparing different sampling numbers t picture;

[0124] Figure 7 is the RMSE of the present invention comparing different numbers of array elementsr picture. DETAILED DESCRIPTION

[0125] Example 1

[0126] This embodiment considers a single-base MIMO radar system composed of multiple transmitting sensors and multiple speed receiving sensors. Assume that the array positions of the transmitting array element and the receiving array element are arbitrary, and the transmitting array element has M transmitting sensors, the receiving array element has N speed receiving sensors, and the position of the mth transmitting sensor of the transmitting array element is (x tm ,y tm ,z tm ), the position of the nth velocity receiving sensor of the receiving array is (x rn ,y rn ,z rn ).like Figure 1 As shown in the figure, assuming that there are K targets in the far field of the MIMO radar, the two-dimensional departure direction angle 2D-DOD of the kth target relative to the transmitting array is (θ tk ,φ tk ), the two-dimensional receiving direction angle 2D-DOA relative to the receiving array is (θ rk ,φ rk ).

[0127] The transmitting antenna used in this embodiment transmits orthogonal waveform signals with the same bandwidth and center frequency, and the wavelength is λ. These signals are reflected by the target and then received by the speed sensor.

[0128] like Figure 1-7 As shown, a 2D-DOA estimation method for an acoustic vector sensor with an arbitrary array structure includes the following steps:

[0129] S1: The receiving end obtains the covariance matrix, performs eigendecomposition on the covariance matrix, and solves the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors;

[0130] S2: Using the signal subspace and the ESPRIT algorithm, a rough estimate of 2D-DOA is obtained;

[0131] S3: Simultaneously use the signal subspace to solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, obtain the ambiguous spatial phase difference matrix;

[0132] S4: Perform phase compensation on the ambiguous spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA;

[0133] Here, 2D-DOA represents two-dimensional angle of arrival.

[0134] This embodiment obtains the covariance matrix, performs eigendecomposition on the covariance matrix, solves the signal subspace, and then roughly estimates the two-dimensional DOA in the vector information based on the signal subspace and uses the ESPRIT algorithm. The rough estimation result is then used to perform phase compensation for spatial ambiguity, thereby achieving accurate estimation of the 2D-DOA of the acoustic vector sensor of any array structure. This embodiment not only utilizes the vector sensor to realize the spatial layout of any array, but also makes the array element spacing of the sensor no longer limited by half a wavelength; it not only utilizes the information of the vector part, but also utilizes the information of the spatial phase, which significantly improves the estimation accuracy. Moreover, the algorithm does not need to perform spectral peak search and additional two-dimensional angle pairing calculation on the covariance matrix of the received signal, which significantly reduces the computational complexity, improves the accuracy of the estimation, and enhances the robustness and applicability of the system.

[0135] Since the array element spacing is greater than half a wavelength, there will be phase ambiguity in the spatial phase difference matrix, and the spatial phase difference matrix is ​​the information in the spatial domain, which means that there is ambiguity in the spatial domain, namely, spatial ambiguity. In view of this feature, this embodiment not only utilizes the information of the vector part to make a rough estimate first, but also utilizes the information of the spatial domain to make a refined estimate, thereby improving the accuracy of the algorithm.

[0136] In the preferred embodiment, the specific steps of step S1 are as follows:

[0137] S101: Based on the fact that the array positions of the transmitting array element and the receiving array element are set arbitrarily, and the transmitting array element has M transmitting sensors and the receiving array element has N speed receiving sensors, the coordinates of the first row of the transmitting array and the receiving array are used as reference array elements, and the arrival time differences relative to the mth transmitting sensor and the nth speed receiving sensor are respectively:

[0138] d t,m,k =x tm cos(φ tk )sin(θ tk )+y tm sin(φ tk )sin(θ tk )+z tm cos(θ tk ) (1);

[0139] d r,n,k =x rn cos(φ rk )sin(θ rk )+y rn sin(φ rk )sin(θ rk )+z rn cos(θ rk ) (2);

[0140] In the formula, (x tm ,y tm ,z tm ) is the position of the mth transmitting sensor, (x rn ,y rn ,z rn ) is the position of the nth velocity receiving sensor, θ tk ,φ tk are the azimuth and elevation of the 2D-DOD of the kth target, θ rk ,φ rk are the azimuth and elevation of the 2D-DOA of the kth target, respectively; where 2D-DOD represents the two-dimensional wave departure angle;

[0141] S102: Under the narrowband assumption, the transmitting sensor array spatial response vector a t (θ tk ,φ tk ) and the receiving sensor array spatial response vector a r (θ rk ,φ rk ) are:

[0142]

[0143] Where, d t,m,k is the arrival time difference of the mth transmitting sensor relative to the reference sensor array element, d r,n,k is the arrival time difference of the nth velocity receiving sensor relative to the reference sensor array element, and λ is the wavelength of the transmitted signal.

[0144] S103: Assume that the distance between the target and the single-base radar system is far enough, so the signal reflected by the target can be regarded as a plane wave. The output formula of the matched filter of the receiving end velocity receiving sensor is:

[0145]

[0146] or:

[0147] Y(t)=[A t ⊙A r ⊙B r ]s(t)+n(t) (6);

[0148] In the formula, A t =[a t (θ t1 ,φ t1 ),a t (θ t2 ,φ t2 ),...,a t (θtK ,φ tK )]∈C M×K is the transmit array spatial response matrix, A r =[a r (θ r1 ,φ r1 ),a r (θ r2 ,φ r2 ),...,a r (θ rK ,φ rK )]∈C N×K is the spatial response matrix of the receiving array, s(t)=[s 1 (t),s 2 (t),…,s K (t)] T ∈C K×L , f k is the Doppler shift, β k is the amplitude, n(t)∈C MN×L is the received Gaussian white noise matrix, ⊙ represents the Khatri-Rao product, represents the Kroc inner product, L is the number of snapshots; B r =[b r1 ,b r2 ,b r3 ] T ∈C 3×K is the polarization response matrix, and

[0149] b r1 =[sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ] T ,

[0150] b r2 =[sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,sinθ rK sinφ rK ] T ,

[0151] b r3 =[cosθ r1 ,cosθ r2 ,...,cosθ rK ] T ;

[0152] S104: Define matrix A according to step S103:

[0153] A=A t ⊙A r ⊙B r (7); then the obtained signal model is:

[0154] Y = AS + N (8);

[0155] According to the above signal model, the covariance matrix R of the received signal can be calculated Y for:

[0156] R Y =YY H / L (9);

[0157] Where, superscript H is the conjugate transpose, and L is the number of snapshots.

[0158] S105: Solve the signal subspace, and calculate the received signal covariance matrix R in equation (9): Y Perform eigendecomposition to obtain the signal subspace. The eigendecomposition formula is:

[0159]

[0160] In the formula, is a K×K diagonal matrix whose diagonal elements contain K large eigenvalues, is a matrix composed of eigenvectors corresponding to K larger eigenvalues, The diagonal elements of contain 3MN-K smaller eigenvalues, is the matrix of eigenvectors corresponding to the remaining smaller eigenvalues, and can be regarded as a signal subspace, It can be regarded as a noise subspace;

[0161] Signal subspace The same as the subspace spanned by the spatial response vector matrix A, that is:

[0162]

[0163] In the formula, the matrix T∈C K×K , is a non-singular full-rank matrix.

[0164] In this embodiment, step S1 provides the necessary signal model and mathematical basis for the subsequent angle of arrival estimation by calculating the arrival time difference, determining the spatial response vector, defining the receiving end output formula, calculating the covariance matrix and solving the signal subspace by eigendecomposition, etc. After the signal subspace is obtained, a rough estimate of 2D-DOA can be performed.

[0165] In the preferred embodiment, the specific steps of step S2 are as follows:

[0166] S201: Set the selectivity matrix, the formula is:

[0167] J 1 =[I M ⊙I N ⊙I 3，1 ] (12a);

[0168] J 2 =[I M ⊙I N ⊙I 3，2 ] (12b);

[0169] J 3 =[I M ⊙I N ⊙I 3，3 ] (12c);

[0170] In the formula, 1 M is the M×M identity matrix, Ι N is the N×N identity matrix, Ι 3,k is the kth row of the 3×3 identity matrix;

[0171] S202: The selectivity matrix is ​​multiplied by the left side of the signal subspace matrix, and A is set 0 =A t ⊙A r ,Right now:

[0172]

[0173] In the formula, B r,k Represented as matrix B r The kth row of , where:

[0174] Φ x =diag(sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ),

[0175] Φ y =diag(sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,sinθ rK sinφ rK ),

[0176] Φ z =diag(cosθ r1 ,cosθ r2 ,...,cosθ rK );

[0177] S203: According to formula (13), the following matrix can be constructed:

[0178]

[0179] In the formula,

[0180] At the same time: let Φ xy =Φ z Φ xyz ,but:

[0181]

[0182] Replace A in formula (13c) 0 Substituting into formula (15), we can get:

[0183]

[0184] In the formula,

[0185] From formula (16), we can get: in represents the generalized inverse of a matrix, [·] -1 It means to find the inverse of the matrix;

[0186] S204: Perform eigendecomposition on the output matrix of step S203 to obtain a rough estimate of 2D-DOA. Specifically, Perform eigendecomposition to obtain eigenvalues and the eigenvector Right now:

[0187]

[0188] Then the rough estimate of the azimuth and elevation angles of 2D-DOA can be expressed as:

[0189]

[0190] In the formula, abs() means taking the amplitude, and angle() means taking the phase.

[0191] In the above steps, a rough estimation of 2D-DOA is achieved by constructing a selectivity matrix, applying the matrix to the signal subspace, constructing and processing a new matrix, and performing eigendecomposition.

[0192] In the preferred solution, to solve the spatial ambiguity, step S3 specifically includes:

[0193] S301: Define the following spatial selectivity matrix, expressed as:

[0194] P n =[I M ☉I N，n ☉I 3 ],n=1,...,N (21)

[0195] In the formula, 1 N,n is the nth row of the N×N identity matrix;

[0196] S302: The spatial selectivity matrix is ​​multiplied by the left side of the signal subspace matrix, that is:

[0197]

[0198] In the formula, is the spatial response vector A of the receiving array r The first line of is the spatial response vector A of the receiving array r The second line of

[0199] S303: and A r,2 =A r,1 Φ r1 , substituting into formula (22) we get:

[0200]

[0201] at this time:

[0202]

[0203] In the formula, u rk =cosφ rk sinθ rk , v rk = sinφ rk sinθ rk , w rk =cosφ rk , diag() means converting the row matrix into a diagonal matrix, let x r1 -x r2 =m 12 ,y r1 -y r2 =n 12 , z r1 -z r2 = l 12 ,but:

[0204]

[0205] S304: From formula (23), we can get Then substitute T in formula (18) to obtain:

[0206]

[0207] Similarly, using the selectivity matrix P 3 Multiplying both sides of equation (11) at the same time, we can get And A r 3 =A r 2 Φ r2 ,Right now:

[0208]

[0209] at this time m 23 =x r2 -x r3 , n 23 =y r2 -y r3 , l 23 =z r2 -z r3 , similarly, substituting T in formula (18) into formula (27) yields:

[0210]

[0211] Same as above, using the selectivity matrix P 4,...,N ,Will The simulation formula (27)(28)(29) is obtained;

[0212] S305: Take The diagonal elements of are used to construct the row vector Then put these N-1 row vectors into the initial phase matrix GG∈C N-1×K You can get:

[0213]

[0214] The initial phase matrix GG is phased to obtain the initial phase difference matrix GGp:

[0215]

[0216] That is: GGp=EP h (31);

[0217] in:

[0218]

[0219] At this time, when the phase is obtained using the built-in angle function in MATLAB, the phase values ​​greater than [-π,π] will be periodically mapped to the range of [-π,π]. Therefore, there will be phase ambiguity in the phase difference matrix GGp of the received adjacent sensors. This phase ambiguity will cause errors in the DOA estimation results. Therefore, the following phase compensation process is required to obtain accurate estimation results.

[0220] In the preferred embodiment, step S4 specifically includes:

[0221] S401: First, using the rough estimation results in formulas (19) and (20), that is, obtaining the rough estimation results of the azimuth and elevation angles of the 2D-DOA in step S3, construct the following unambiguous and low-precision Poynting vector matrix: The expression is:

[0222]

[0223] S402: Use the unambiguous low-precision Poynting vector matrix By using formula (31), the unambiguous and low-precision phase difference matrix is ​​obtained: Through this unambiguous low-precision phase difference matrix Get blur Solve it, that is:

[0224]

[0225] S403: Blur By compensating to the periodically blurred phase difference matrix GGp through equation (34), we can obtain the unambiguous and high-precision phase difference matrix GGpt, namely:

[0226]

[0227] S404: Using the least squares principle, estimate the unambiguous and high-precision Poynting vector matrix Right now:

[0228]

[0229] Thus, the azimuth and elevation angles of 2D-DOA are obtained, that is, the azimuth and elevation angle estimation formulas of 2D-DOA are:

[0230]

[0231] The usefulness of the scheme of the present invention is illustrated by a specific example. In the monostatic MIMO radar system, the number of transmitting sensors M=10, the number of velocity receiving sensors N=10, and the wavelength λ of the antenna operation are set to 1. Assume that there are K=3 targets within the range of the monostatic MIMO radar system, and the 2D-DOA of the three targets is θ r =(20°,36°,65°),φ r =(5°, 15°, 25°). Under the conditions of sampling number L = 100 and signal-to-noise ratio SNR = 12dB, the angle estimation is performed using the method proposed in the scheme of the present invention. The scatter plot results of 100 Monte Carlo test simulations are shown in the attached figure. Figure 2 As shown in the figure, in order to highlight the reliability of the solution of the present invention, the angle estimation is compared with the method of obtaining 2D-DOA by using a velocity receiving sensor through an ESPRIT-Like algorithm, as shown in the attached figure. Figure 3 As shown, it can be seen Figure 2 Shows more Figure 3 These findings indicate that the scheme of the present invention can effectively estimate the receiving direction angle 2D-DOA of the target signal relative to the velocity receiving sensor array, and the performance of the estimation algorithm is also significantly improved.

[0232] The preferred solution also includes: using the root mean square error (RMSE) for evaluation, and the calculation formula is:

[0233]

[0234] In the formula, is θ rk The estimate, is φ rk Estimates.

[0235] like Figure 4 As shown in the figure, the algorithm proposed in this embodiment compares the RMSE under different signal-to-noise ratios SNR under the conditions of sampling number L=100, number of transmitting sensors M=8, and number of speed receiving sensors N=8. r Figure. In order to highlight the reliability of the solution of the present invention, the ESPRIT_Like algorithm is used for comparison. At this point, it can be clearly seen that no matter how the signal-to-noise ratio SNR increases, the accuracy of the ESPRIT_Like algorithm is not very high. However, the algorithm proposed in the present invention, as the signal-to-noise ratio SNR increases, the estimation accuracy will become higher and higher, and is much higher than the ESPRIT_Like algorithm.

[0236] like Figure 5As shown in the figure, the algorithm proposed in the present invention compares the RMSE under different sampling sample numbers L under the conditions of the number of transmitting sensors M = 8, the number of speed receiving sensors N = 8, and the signal-to-noise ratio SNR = 5dB. r At this point, it can be clearly seen that as the number of samples L increases, the speed at which the accuracy of the ESPRIT_Like algorithm increases becomes slower and slower. However, the accuracy of the algorithm proposed in the present invention will become higher and higher, and is always higher than that of the ESPRIT_Like algorithm.

[0237] like Figure 6 As shown in the figure, the algorithm proposed in the present invention compares the RMSE of different array element numbers under the conditions of sampling number L = 100 and signal-to-noise ratio SNR = 10dB. r At this time, as before, the algorithm proposed by the present invention will have higher and higher estimation accuracy as the number of array elements increases, and will always be higher than the ESPRIT_Like algorithm.

[0238] This embodiment has the following beneficial effects:

[0239] 1. Different from the traditional angle estimation based on scalar sensor array, the use of vector sensor array can realize arbitrary layout of sensor spatial position, such as placing sensors on the surface of the aircraft to reduce radar cross section and enhance the stealth and maneuverability of the aircraft, and the spacing between sensor array elements is no longer limited by half wavelength, reducing sensor position error and mutual coupling effect, and improving the accuracy of the radar system.

[0240] 2. Different from the previous angle estimation based on electromagnetic vector sensor array, the speed receiving sensor is used. The hardware cost of this sensor is low, and it only has three vector channels. The calculation complexity is also low. Its low cost and high speed advantages can be used on a large scale.

[0241] 3. The embodiment not only utilizes the information of the vector part, but also fully utilizes the spatial phase information of the sensor array, which significantly improves the estimation accuracy; the proposed algorithm does not require spectral peak search and additional angle pairing calculations on the covariance matrix of the received signal, which significantly reduces the computational complexity.

[0242] Example 2

[0243] In combination with Example 1, a 2D-DOA estimation system of an acoustic vector sensor with an arbitrary array structure is proposed, comprising:

[0244] The signal subspace module is used for obtaining the covariance matrix at the receiving end, and performing eigendecomposition on the covariance matrix to solve the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors;

[0245] A rough estimation module is used to obtain a rough estimate of 2D-DOA using the signal subspace and the ESPRIT algorithm;

[0246] The spatial ambiguity module is used to simultaneously utilize the signal subspace and solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, to obtain the ambiguous spatial phase difference matrix;

[0247] An accurate estimation module is used to perform phase compensation on the ambiguous spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA;

[0248] Here, 2D-DOA represents two-dimensional angle of arrival.

[0249] The above embodiments are only preferred technical solutions of the present invention and should not be regarded as limiting the present invention. The protection scope of the present invention shall be the technical solutions recorded in the claims, including equivalent replacement solutions of the technical features in the technical solutions recorded in the claims. That is, equivalent replacement improvements within this scope are also within the protection scope of the present invention.

Claims

1. A 2D-DOA estimation method for an acoustic vector sensor of arbitrary array structure, characterized in that: The following steps are involved: S1: The receiving end obtains the covariance matrix, performs eigendecomposition on the covariance matrix, and solves the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors; S2: Using the signal subspace and the ESPRIT algorithm, a rough estimate of 2D-DOA is obtained; S3: Simultaneously use the signal subspace to perform matrix transformation and solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, obtain the ambiguous spatial phase difference matrix; S4: Perform phase compensation on the ambiguous spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA; Here, 2D-DOA represents two-dimensional angle of arrival.

2. According to the method for estimating 2D-DOA of an acoustic vector sensor with an arbitrary array structure according to claim 1, it is characterized in that: The specific steps of step S1 are as follows: S101: Based on the fact that the array positions of the transmitting array element and the receiving array element are set arbitrarily, and the transmitting array element has M transmitting sensors and the receiving array element has N speed receiving sensors, the coordinates of the first row of the transmitting array and the receiving array are used as reference array elements, and the arrival time differences relative to the mth transmitting sensor and the nth speed receiving sensor are respectively: d t,m,k =x tm cos(φ tk )sin(θ tk )+y tm sin(φ tk )sin(θ tk )+z tm cos(θ tk ) (1); d r,n,k =x rn cos(φ rk )sin(θ rk )+y rn sin(φ rk )sin(θ rk )+z rn cos(θ rk ) (2); In the formula, (x tm ,y tm ,z tm ) is the position of the mth transmitting sensor, (x rn ,y rn ,z rn ) is the position of the nth velocity receiving sensor, θ tk ,φ tk are the azimuth and elevation of the 2D-DOD of the kth target, θ rk ,φ rk are the azimuth and elevation of the 2D-DOA of the kth target, respectively; where 2D-DOD represents the two-dimensional wave departure angle; S102: Under the narrowband assumption, the transmitting sensor array spatial response vector a t (θ tk ,φ tk ) and the receiving sensor array spatial response vector a r (θ rk ,φ rk ) are: Where, d t,m,k is the arrival time difference of the mth transmitting sensor relative to the reference sensor array element, d r,n,k is the arrival time difference of the nth velocity receiving sensor relative to the reference sensor array element, λ is the wavelength of the transmitted signal; S103: The output formula of the matching filter of the receiving end speed receiving sensor is: or: Y(t)=[A t ☉A r ☉B r ]s(t)+n(t) (6); In the formula, A t =[a t (θ t1 ,φ t1 ),a t (θ t2 ,φ t2 ),...,a t (θ tK ,φ tK )]∈C M×K is the transmit array spatial response matrix, A r =[a r (θ r1 ,φ r1 ),a r (θ r2 ,φ r2 ),...,a r (θ rK ,φ rK )]∈C N×K is the spatial response matrix of the receiving array, s(t)=[s1(t),s2(t),…,s K (t)] T ∈C K ×L, f k is the Doppler shift, β k is the amplitude, n(t)∈C MN×L is the received Gaussian white noise matrix, ⊙ represents the Khatri-Rao product, represents the Kroc inner product, L is the number of snapshots; B r =[b r1 ,b r2 ,b r3 ] T ∈C 3×K is the polarization response matrix, and b r1 =[sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ] T , b r2 =[sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,sinθ rK sinφ rK ] T , b r3 =[cosθ r1 ,cosθ r2 ,...,cosθ rK ] T ; S104: Define matrix A according to step S103: A=A t ☉A r ☉B r (7); The obtained signal model is: Y = AS + N (8); Calculate the covariance matrix R of the received signal Y for: R Y =YY H / L (9); In the formula, the superscript H is the conjugate transpose, and L is the number of snapshots; S105: Solve the signal subspace, and calculate the received signal covariance matrix R in equation (9): Y Perform eigendecomposition to obtain the signal subspace. The eigendecomposition formula is: In the formula, is a K×K diagonal matrix whose diagonal elements contain K large eigenvalues, is a matrix composed of eigenvectors corresponding to K larger eigenvalues, The diagonal elements of contain 3MN-K smaller eigenvalues, is the matrix of eigenvectors corresponding to the remaining smaller eigenvalues, and can be regarded as a signal subspace, It can be regarded as a noise subspace; Signal subspace The same as the subspace spanned by the spatial response vector matrix A, that is: In the formula, the matrix T∈C K×K , is a non-singular full-rank matrix.

3. According to the method for estimating 2D-DOA of an acoustic vector sensor with an arbitrary array structure as described in claim 1, it is characterized in that: The specific steps of step S2 are as follows: S201: Set the selectivity matrix, the formula is: J1=[I M ☉I N ☉I 3，1 ] (12a); J2=[I M ☉I N ☉I 3，2 ] (12b); J3=[I M ☉I N ☉I 3，3 ] (12c); In the formula, 1 M is the M×M identity matrix, Ι N is the N×N identity matrix, Ι 3,k is the kth row of the 3×3 identity matrix; S202: Multiply the selectivity matrix by the left side of the signal subspace matrix, and let A0 = A t ⊙A r ,Right now: In the formula, B r ,k is represented by the matrix B r The kth row of , where: Φ x =diag(sinθ r1 cosφ r1 ,sinθ r2 cosφ r2 ,...,sinθ rK cosφ rK ), F y =diag(sinθ r1 sinφ r1 ,sinθ r2 sinφ r2 ,...,synth rK sinφ rK ), F z =diag(cosθ r1 ,cosθ r2 ,...,cosθ rK ); S203: According to formula (13), the following matrix can be constructed: In the formula, At the same time: let Φ xy =Φ z Φ xyz ,but: Substituting A0 in formula (13c) into formula (15), we can obtain: In the formula, From formula (16), we can get: in represents the generalized inverse of a matrix, [·] -1 It means to find the inverse of the matrix; S204: Perform eigendecomposition on the output matrix of step S203 to obtain a rough estimate of 2D-DOA. Specifically, Perform eigendecomposition to obtain eigenvalues and the eigenvector Right now: Then the rough estimate of the azimuth and elevation angles of 2D-DOA can be expressed as: In the formula, abs() means taking the amplitude, and angle() means taking the phase.

4. According to the method for estimating 2D-DOA of an acoustic vector sensor with an arbitrary array structure as described in claim 1, it is characterized in that: To solve the spatial ambiguity, step S3 specifically includes: S301: Define the following spatial selectivity matrix, expressed as: P n =[I M ☉I N,n ☉I3],n=1,...,N (21); In the formula, 1 N,n is the nth row of the N×N identity matrix; S302: Multiply the selectivity matrix by the left side of the signal subspace matrix, that is: In the formula, is the spatial response vector A of the receiving array r The first line of is the spatial response vector A of the receiving array r The second line of S303: and A r,2 =A r,1 Φ r1 , substituting into formula (22) we get: at this time: In the formula, u rk =cosφ rk sinθ rk , v rk = sinφ rk sinθ rk , w rk =cosφ rk , diag() means converting the row matrix into a diagonal matrix, let x r1 -x r2 =m 12 ,y r1 -y r2 =n 12 , z r1 -z r2 = l 12 ,but: S304: From formula (23), we can get Then substitute T in formula (18) to obtain: Similarly, by multiplying both sides of equation (11) with the selectivity matrix P3, we can get And A r,3 =A r,2 Φ r2 ,Right now: at this time m 23 =x r2 -x r3 , n 23 =y r2 -y r3 , l 23 =z r2 -z r3 , similarly, substituting T in formula (18) into formula (27) yields: Same as above, using the selectivity matrix P 4,...,N ,Will The simulation formula (27)(28)(29) is obtained; S305: Take The diagonal elements of are used to construct the row vector Then put these N-1 row vectors into the initial phase matrix GG∈C N-1×K You can get: The initial phase matrix GG is phased to obtain the initial phase difference matrix GGp: That is: GGp=EP h (31); in:

5. According to the method for estimating 2D-DOA of an acoustic vector sensor with an arbitrary array structure as described in claim 1, it is characterized in that: Step S4 specifically includes: S401: First, using the rough estimation results in formulas (19) and (20), that is, obtaining the rough estimation results of the azimuth and elevation angles of the 2D-DOA in step S3, construct the following unambiguous and low-precision Poynting vector matrix: The expression is: S402: Use the unambiguous low-precision Poynting vector matrix By using formula (31), the unambiguous and low-precision phase difference matrix is ​​obtained: Through this unambiguous low-precision phase difference matrix Get blur Solve it, that is: S403: Blur By compensating to the periodically blurred phase difference matrix GGp through equation (34), the unambiguous and high-precision phase difference matrix GGpt is obtained, namely: S404: Using the least squares principle, estimate the unambiguous and high-precision Poynting vector matrix Right now: Thus, the azimuth and elevation angles of 2D-DOA are obtained, that is, the azimuth and elevation angle estimation formulas of 2D-DOA are:

6. According to the method for estimating 2D-DOA of an acoustic vector sensor with an arbitrary array structure as described in claim 1, it is characterized in that: Also includes: The root mean square error is used for evaluation and is calculated as: In the formula, is θ rk The estimate, is φ rk Estimates.

7. A 2D-DOA estimation system for acoustic vector sensors of arbitrary array structure, characterized in that: include: The signal subspace module is used for obtaining the covariance matrix at the receiving end, and performing eigendecomposition on the covariance matrix to solve the signal subspace; wherein the receiving end adopts a vector sensor array, and the vector sensor array includes a plurality of speed receiving sensors; A rough estimation module is used to obtain a rough estimate of 2D-DOA using the signal subspace and the ESPRIT algorithm; The spatial ambiguity module is used to simultaneously utilize the signal subspace and solve the spatial ambiguity through an algorithm similar to ESPRIT, that is, to obtain the ambiguous spatial phase difference matrix; An accurate estimation module is used to perform phase compensation on the spatial phase difference matrix according to the rough estimation result of step 2 to obtain an accurate estimation of 2D-DOA; Here, 2D-DOA represents two-dimensional angle of arrival.

Images