Parallel hydrophone array DOA estimation method based on extended covariance and propagation operator

Abstract

The invention discloses a parallel hydrophone array DOA estimation method based on extended covariance and a propagation operator. The method comprises the following steps: firstly, constructing an expanded array receiving data matrix by using a sensor array receiving data matrix, then defining a propagation matrix and expanding the propagation matrix, further screening array manifolds of the expanded array receiving data matrix and constructing a linear relation matrix, and then combining an expanded propagation operator matrix to obtain a linear relation matrix; and finally, realizing joint estimation of the azimuth angle and the pitch angle of the incident signal by utilizing linear segmentation operation. Compared with a traditional DOA estimation method, the DOA estimation method has the advantages that array aperture characteristics are fully mined by fusing autocorrelation and cross covariance information of double parallel linear arrays, the information utilization rate is increased, DOA estimation precision and estimation performance of the parallel hydrophone array are effectively improved, parameter matching and an EVD process of a cross-correlation matrix are not needed, and the DOA estimation precision and estimation performance of the parallel hydrophone array are effectively improved. And spectral peak search is avoided, the calculation complexity is reduced, and the method has a relatively good engineering application value.

Description

Technical Field

[0001] The present invention belongs to the field of sensor array signal processing technology in the field of signal and information processing, and relates to a method for estimating the direction of arrival of an incident signal based on a parallel hydrophone array. Specifically, the DOA estimation of the parallel hydrophone array is achieved through an extended array received data class matrix and a propagation operator matrix. Background Art

[0002] Direction of Arrival (DOA) estimation is a research hotspot in array signal processing and is widely used in communications, radar detection, target positioning, and other fields. By processing the array's received signals, position information such as the azimuth and elevation angle of the incident signal source can be obtained. Due to their structural symmetry and translational invariance, dual parallel linear arrays exhibit significant advantages in two-dimensional direction of arrival estimation. Over the past few decades, researchers at home and abroad have proposed numerous two-dimensional DOA estimation methods, including two-dimensional multi-signal classification algorithms, signal parameter estimation algorithms based on rotational invariance techniques, and propagation operator algorithms. However, some of these methods suffer from problems such as spectral peak searching, high computational complexity, and the inability to achieve two-dimensional angle matching.

[0003] To address the above issues, Zhang Xiaofei et al. proposed a two-dimensional DOA estimation method (DOA-Matrix method) for dual parallel linear arrays. They constructed a DOA matrix using the properties of a partial covariance matrix and performed eigenvalue decomposition (EVD) on the DOA matrix to achieve two-dimensional DOA estimation. This method avoids spectral peak search, but ignores the information of some autocorrelation matrices and cross-covariance matrices, which reduces the estimation accuracy to a certain extent. Based on the traditional DOA-Matrix method, Dai Xiangrui et al. constructed an extended DOA matrix (EDOA-Matrix method) containing all autocorrelation and cross-correlation information. By obtaining eigenvalues ​​and eigenvectors from the extended DOA matrix, a more accurate DOA estimation can be obtained. However, this method inevitably increases the computational complexity. Sun Gaoli et al. proposed a two-dimensional DOA estimation method based on an extended covariance matrix of dual parallel linear arrays. This method constructs an extended covariance matrix based on the array received data, reorganizes the signal subspace of the extended matrix and constructs a relationship matrix. At the same time, the extended array manifold matrix is ​​screened and combined with the rotation invariance technology (ESPRIT) to achieve joint estimation of azimuth and elevation angles. Although this method further improves the DOA estimation accuracy, the block reorganization of the signal subspace also increases the computational complexity. This patent addresses the problems of high computational complexity of existing two-dimensional DOA estimation methods for dual parallel linear arrays and the inability to automatically match azimuth and elevation angles. A parallel hydrophone array DOA estimation method based on extended covariance and propagation operators is proposed. This method effectively improves the DOA estimation accuracy of the parallel hydrophone array, achieves automatic matching of two-dimensional angles, and reduces computational complexity. Summary of the Invention

[0004] This patent addresses the shortcomings of existing two-dimensional DOA estimation methods for dual parallel linear arrays and proposes a parallel hydrophone array DOA estimation method based on extended covariance and propagation operators. Based on the data matrix received by the sensor array, the invention constructs an extended array reception data class matrix, screens the corresponding array manifold, and then calculates the extended propagation matrix. By leveraging the relationship between the extended propagation matrix and the array manifold, the invention utilizes linear segmentation operations to obtain a joint estimate of the azimuth and elevation angles of the incident signal. This patented method achieves automatic matching of azimuth and elevation angles, improves the estimation accuracy of two-dimensional DOA, avoids spectral peak searching, reduces computational complexity, and has good engineering application value.

[0005] A parallel hydrophone array DOA estimation method based on extended covariance and propagation operator includes the following steps:

[0006] 1. A method for estimating the direction of arrival (DOA) of a parallel hydrophone array based on extended covariance and a propagation operator, characterized in that the parallel hydrophone array comprises two uniform linear subarrays 1 and 2, each subarray having M array elements, a spacing d between array elements in subarray 1 and subarray 2, and a spacing d between subarray 1 and subarray 2, where d = λ / 2, and λ represents the wavelength of the signal. The method comprises the following steps:

[0007] Step 1: The received data of subarray 1 and subarray 2 are expressed as Y(t)=A y s(t)+N y (t), Z(t)=A z s(t)+N z (t), where is the direction matrix of sub-array 1, K represents the number of incident signals, The superscript "T" indicates the matrix transpose operation, and the signal label k = 1, 2, ..., K, A z is the direction matrix of sub-matrix 2, A z =A y U, is the rotation factor matrix, diag{·} represents the diagonal matrix, θ k 、 are the elevation and azimuth angles of the kth incident signal, Represents the signal data vector received by the array at time t, N y (t) and N z (t) are the additive white Gaussian noise data vectors of sub-matrix 1 and sub-matrix 2 at time t;

[0008] Step 2: Calculate the covariance matrix R yy =E[Y(t)Y H (t)]=A y R s A y H +I M σ 2 ,

[0009] where R s =E[s(t)s H (t)] is the covariance matrix of the signal, σ 2 is the noise power of the array, where [·] H represents the conjugate transpose of the matrix, E[·] represents the calculation of the statistical mean, I M represents the identity matrix of rank M;

[0010] Step 3: Calculate the estimated value of the covariance matrix in step 2 Where L is the number of snapshots;

[0011] Step 4: Perform eigendecomposition and average the obtained MK smaller eigenvalues ​​to obtain the estimated value of the noise power right After removing the noise, we get

[0012] Step 5: Construct a new expanded covariance matrix in Where [·] T Represents the transpose of a matrix;

[0013] Step 6: Find the autocorrelation matrix and get right Divide into blocks, in The dimension is 3M×K, The dimension is 3M×(3M-K);

[0014] Step 7: Construct the matrix Further obtain P e =[I K ,P] H , I K represents the identity matrix of rank K;

[0015] Step 8: Construct the selection matrix I u1 =[O 2M×M ,I 2M ], I u2 =[I 2M ,O 2M×M ], where I3 represents the identity matrix of rank 3, I M-1 represents the identity matrix of rank M-1, I 2M represents the identity matrix of rank 2M, O (M-1)×1 represents an all-zero matrix of dimension (M-1)×1, O 2M×M represents an all-zero matrix of dimension 2M×M, represents the Kronecker product;

[0016] Step 9: Construct the matrix (·) + represents the generalized inverse of a matrix;

[0017] Step 10: Perform eigenvalue decomposition on the matrix π to obtain the diagonal element eigenvalues ​​φ1, φ2,…, φ K , the eigenvector in

[0018] Step 11: Construction The eigenvalues ​​of the diagonal elements are obtained by eigendecomposition: γ1,γ2,…γ K ,in

[0019] Step 12: According to φ k ,γ k An estimate of the azimuth angle can be calculated and the estimated pitch angle

[0020] Compared with the existing technology, the technical solution of this patent has the following technical effects:

[0021] (1) The present invention constructs an extended covariance matrix based on the received data matrix, thereby improving the DOA estimation accuracy of the parallel sensor array.

[0022] (2) The method of the present invention does not require spectral peak search and signal subspace division, and has low computational complexity.

[0023] (3) The method of the present invention realizes automatic matching of azimuth angle and pitch angle. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 This is the actual scenario application diagram of this system;

[0025] Figure 2 It is the two-dimensional DOA estimation result of the patented signal processing method under two signal sources;

[0026] Figure 3 This is the relationship curve between the root mean square error of the azimuth angle and the signal-to-noise ratio of the signal processing method of this patent;

[0027] Figure 4 This is the relationship curve between the root mean square error of the pitch angle and the signal-to-noise ratio of the signal processing method of this patent;

[0028] Figure 5 This is the relationship curve between the root mean square error of the azimuth angle and the number of snapshots of the signal processing method of this patent.

[0029] Figure 6 This is the relationship curve between the root mean square error of the pitch angle and the number of snapshots of the signal processing method of this patent. DETAILED DESCRIPTION

[0030] The present invention will now be further described with reference to the embodiments and accompanying drawings:

[0031] The first embodiment: Figure 1The array model of the parallel hydrophone used in this invention is presented. The parallel hydrophone array consists of subarray 1 and subarray 2. Each subarray has 12 elements, i.e., M = 12. The spacing between adjacent elements and the two parallel linear arrays is 0.03 m. The number of signals is K = 2, the azimuth angle is (40°, 60°), the elevation angle is (3°, 20°), the signal-to-noise ratio is 10 dB, and the number of snapshots is L = 500.

[0032] The above conditions are used to estimate the direction of arrival. The specific implementation process is as follows:

[0033] Step 1: The received data of subarray 1 and subarray 2 are expressed as Y(t)=A y s(t)+N y (t), Z(t)=A z s(t)+N z (t), where is the direction matrix of sub-array 1, K=2 represents the number of incident signals, The superscript "T" indicates the matrix transposition operation, the signal number k = 1, 2, A z is the direction matrix of sub-matrix 2, A z =A y U, is the rotation factor matrix, diag{·} represents the diagonal matrix, θ k 、 are the elevation and azimuth angles of the kth incident signal, Represents the signal data vector received by the array at time t, N y (t) and N z (t) are the additive white Gaussian noise data vectors of sub-matrix 1 and sub-matrix 2 at time t;

[0034] Step 2: Calculate the covariance matrix R yy =E[Y(t)Y H (t)]=A y R s A y H +I 12 σ 2 ,

[0035] where R s =E[s(t)s H (t)] is the covariance matrix of the signal, σ 2 is the noise power of the array, where [·] H represents the conjugate transpose of the matrix, E[·] represents the calculation of the statistical mean, I 12 represents the identity matrix of rank 12;

[0036] Step 3: Calculate the estimated value of the covariance matrix in step 2 The number of snapshots is 500;

[0037] Step 4: Perform eigendecomposition and average the 10 smaller eigenvalues ​​to obtain an estimate of the noise power. right After removing the noise, we get

[0038] Step 5: Construct a new expanded covariance matrix in Where [·] T Represents the transpose of a matrix;

[0039] Step 6: Find the autocorrelation matrix and get right Divide into blocks, in The dimension is 36×2, The dimensions are 36×34;

[0040] Step 7: Construct the matrix Further obtain P e =[I2,P] H , I2 represents the identity matrix of rank 2;

[0041] Step 8: Construct the selection matrix I u1 =[O 24×12 ,I 24 ], I u2 =[I 24 ,O 24×12 ], where I3 represents the identity matrix of rank 3, I 11 represents the identity matrix of rank 11, I 24 represents the identity matrix of rank 24, O 11×1 represents an all-zero matrix of dimension 11×1, O 24×12 represents an all-zero matrix of dimension 24×12, represents the Kronecker product;

[0042] Step 9: Construct the matrix (·) + represents the generalized inverse of a matrix;

[0043] Step 10: Perform eigenvalue decomposition on the matrix π to obtain the diagonal element eigenvalues ​​φ1, φ2, and eigenvectors in

[0044] Step 11: Construction The eigenvalues ​​of the diagonal elements are obtained by eigendecomposition, γ1, γ2, where

[0045] Step 12: According to φ k ,γ k An estimate of the azimuth angle can be calculated and the estimated pitch angle

[0046] The patent method based on the above conditions is simulated by MATLAB simulation software, and the angle estimation results of the patent method are obtained as follows: Figure 2 As shown in Figure 2, the two-dimensional DOA estimation results of other methods (EDOA-Matrix method and ECDOA-Matrix method) are also given. When the input azimuth angle is (40°, 60°) and the pitch angle is (3°, 20°), the estimated angles are azimuth angle (3.005°, 19.994°) and pitch angle (40.014°, 59.922°). Figure 2 It can be seen that the experimental results show that the patented method can accurately locate the location of the signal source.

[0047] Example 2: The relationship curve between the root mean square error and the signal-to-noise ratio of the signal processing method of this patent is as follows: Figure 3 and Figure 4 As shown. A comparison chart of the performance analysis of other methods is also given. The application conditions of the method of the present invention are as follows:

[0048] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1 As shown in FIG, the number of array elements in each subarray is a uniform linear array of M=12, the array element spacing is half a wavelength, i.e., d=0.03m, the incident angle of the signal source is azimuth (40°, 60°), and the pitch angle is (3°, 20°), the snapshot number L=500 is used, the signal-to-noise ratio changes from 0dB and increases to 10dB with a step size of 2dB each time, 1000 independent Monte Carlo experiments are performed, and simulation analysis is performed using MATLAB simulation software. The performance analysis curve obtained is shown in FIG. Figure 3 and Figure 4 shown.

[0049] Depend on Figure 3 It can be seen that the root mean square errors of the azimuth angles of the three methods all decrease with the increase of the signal-to-noise ratio, among which the root mean square error of the method of the present invention is the smallest. From the overall effect, the performance of the method of the present invention is better than the DOA estimation results of the other two methods, and the azimuth angle of the signal can be accurately estimated.

[0050] Depend on Figure 4It can be seen that the root mean square errors of the pitch angles of the three methods all decrease with the increase of the signal-to-noise ratio, among which the root mean square error of the method of the present invention is the smallest. From the overall effect, the performance of the method of the present invention is better than the DOA estimation results of the other two methods, and can accurately estimate the pitch angle of the signal.

[0051] Example 3: The relationship curve between the root mean square error and the number of snapshots of the signal processing method of this patent is as follows: Figure 5 and Figure 6 As shown. A comparison chart of the performance analysis of other methods (EDOA-Matrix method and ECDOA-Matrix method) is also given. The application conditions of the method of the present invention are as follows:

[0052] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1 As shown in the figure, the number of array elements in each subarray is a uniform linear array of M=12, the array element spacing is half a wavelength, that is, d=0.03m, the incident angle of the source is (40°, 60°) in azimuth and (3°, 20°) in elevation, the signal-to-noise ratio is 10dB, the number of snapshots starts from 100 and increases to 500 in steps of 100 each time, 1000 independent Monte Carlo experiments are performed, and simulation analysis is performed using MATLAB simulation software. The performance analysis curve obtained is shown in the figure. Figure 5 and Figure 6 shown.

[0053] Depend on Figure 5 It can be seen that the root mean square errors of the azimuth angles of the three methods all decrease with the increase of the number of snapshots, among which the root mean square error of the method of the present invention is the smallest. From the overall effect, the performance of the method of the present invention is better than the DOA estimation results of the other two methods, and the azimuth angle of the signal can be accurately estimated.

[0054] Depend on Figure 6 It can be seen that the root mean square errors of the pitch angles of the three methods all decrease with the increase of the number of snapshots, among which the root mean square error of the method of the present invention is the smallest. From the overall effect, the performance of the method of the present invention is better than the DOA estimation results of the other two methods, and it can accurately estimate the pitch angle of the signal.

[0055] The specific examples described herein are merely illustrative of the present invention. Persons skilled in the art may make various modifications, additions, or substitutions to the specific examples described herein without departing from the present invention or exceeding the scope of the appended claims.

Claims

1. A method for DOA estimation of a parallel hydrophone array based on extended covariance and propagation operators. The parallel hydrophone array consists of two uniform linear subarrays 1 and 2. The number of array elements in each subarray is M. The spacing between the array elements of subarray 1 and subarray 2 is d, and the spacing between subarray 1 and subarray 2 is also d, where d = λ / 2, where λ represents the wavelength of the signal. The method is characterized by: A parallel hydrophone array DOA estimation method based on extended covariance and propagation operator includes the following steps: Step 1: The received data of subarray 1 and subarray 2 are expressed as Y(t)=A y s(t)+N y (t), Z(t)=A z s(t)+N z (t), where is the direction matrix of sub-array 1, K represents the number of incident signals, The superscript "T" indicates the matrix transpose operation, and the signal label k = 1, 2, ..., K, A z is the direction matrix of sub-matrix 2, A z =A y U, is the rotation factor matrix, diag{·} represents the diagonal matrix, θ k 、 are the elevation and azimuth angles of the kth incident signal, Represents the signal data vector received by the array at time t, N y (t) and N z (t) are the additive white Gaussian noise data vectors of sub-matrix 1 and sub-matrix 2 at time t; Step 2: Calculate the covariance matrix R yy =E[Y(t)Y H (t)]=A y R s A y H +I M σ 2 , where R s =E[s(t)s H (t)] is the covariance matrix of the signal, σ 2 is the noise power of the array, where [·] H represents the conjugate transpose of the matrix, E[·] represents the calculation of the statistical mean, I M represents the identity matrix of rank M; Step 3: Calculate the estimated value of the covariance matrix in step 2 Where L is the number of snapshots; Step 4: Perform eigendecomposition and average the obtained MK smaller eigenvalues ​​to obtain the estimated value of the noise power right After removing the noise, we get Step 5: Construct a new expanded covariance matrix in Where [·] T Represents the transpose of a matrix; Step 6: Find the autocorrelation matrix and get right Divide into blocks, in The dimension is 3M×K, The dimension is 3M×(3M-K); Step 7: Construct the matrix Further obtain P e =[I K ,P] H , I K represents the identity matrix of rank K; Step 8: Construct the selection matrix I u1 =[O 2M×M ,I 2M ], I u2 =[I 2M ,O 2M×M ], where I3 represents the identity matrix of rank 3, I M-1 represents the identity matrix of rank M-1, I 2M represents the identity matrix of rank 2M, O (M-1)×1 represents an all-zero matrix of dimension (M-1)×1, O 2M×M represents an all-zero matrix of dimension 2M×M, represents the Kronecker product; Step 9: Construct the matrix (·) + represents the generalized inverse of a matrix; Step 10: Perform eigenvalue decomposition on the matrix π to obtain the diagonal element eigenvalues ​​φ1, φ2,…, φ K , the eigenvector in Step 11: Construction The eigenvalues ​​of the diagonal elements are obtained by eigendecomposition: γ1,γ2,…γ K ,in Step 12: According to φ k ,γ k An estimate of the azimuth angle can be calculated and the estimated pitch angle

Images