Robust DOA estimation method based on parallel hydrophone array

Abstract

The invention discloses a robust DOA estimation method based on a parallel hydrophone array. Two auxiliary array elements are respectively arranged in a subarray 1 and a subarray 2 of the parallel sensor array, and the auxiliary array elements are accurately corrected. The method comprises the following steps: estimating a phase correction vector of a sensor array by utilizing the property of a covariance matrix of a parallel sensor array, then calculating the average energy of collected signals, calculating an amplitude correction vector of the sensor array by utilizing the average energy, and performing amplitude phase correction on array receiving data by utilizing the correction vector to obtain an amplitude phase correction result; and expanding the received data matrix after amplitude-phase correction, generating a new covariance matrix, performing block processing on the covariance matrix, constructing two new matrixes, calculating eigenvalues, and estimating an azimuth angle and a pitch angle by using the eigenvalues. According to the method, the DOA estimation robustness, estimation precision and estimation performance of the parallel sensor array are effectively improved, and the method has 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 robust azimuth estimation method based on a parallel hydrophone array. Specifically, an auxiliary array element method is used to correct the phase error of a subarray, and an energy normalization method is used to correct the amplitude error of the subarray. A received data matrix after amplitude and phase error correction is expanded to construct a new covariance matrix, which is processed in blocks to construct two new matrices, eigenvalues ​​are calculated by eigendecomposition, and the DOA estimation of the parallel hydrophone array is realized by using the obtained method. Background Art

[0002] Direction of Arrival (DOA) estimation is one of the research hotspots in the field of array signal processing. It is widely used in the fields of communication, radar detection, target positioning, etc. By processing the array received signal, the position information such as the azimuth angle and pitch angle of the incident source is obtained. In the past few decades, domestic and foreign researchers have proposed many two-dimensional DOA estimation methods, including two-dimensional multi-signal classification algorithms, signal parameter estimation algorithms based on rotational invariance technology, and propagation operator algorithms. Some methods have problems such as spectral peak search, large amount of calculation, and inability to achieve two-dimensional angle matching. In addition, almost all DOA estimation algorithms are based on the accurate knowledge of the array flow pattern. In actual engineering applications, the influence of factors such as the processing technology level and the complex environment around the experiment may cause array errors (amplitude and phase errors, mutual coupling between array elements, and array element position errors, etc.), resulting in a certain deviation or disturbance between the actual array flow pattern and the ideal array flow pattern, which deteriorates the performance of most existing DOA estimation algorithms or even makes them fail.

[0003] In response to the above problems, Sun Xinyu et al. proposed a PM algorithm that does not require spectral peak search and has better performance than the propagation operator algorithm. However, the algorithm operates in the complex domain and the calculation is still relatively complex. Regarding the amplitude and phase error correction algorithm of two-dimensional arrays, Lu Zukun et al. proposed using L-type arrays and uniform circular arrays to self-correct the amplitude and phase errors, and obtained better estimated parameters. However, research on amplitude and phase error correction methods for parallel arrays has been lacking. This patent proposes a robust DOA estimation method based on a parallel hydrophone array to address the problems of amplitude and phase error correction, high computational complexity, and inability to automatically match azimuth and elevation angles in parallel sensor arrays. The auxiliary array element method is used, that is, two auxiliary array elements are assumed in subarray 1 and subarray 2 respectively. The amplitude and phase errors of the auxiliary array elements are precisely known, and the phase of the array elements is corrected based on the properties of the covariance matrix. The error is then calculated, and the average energy of each array element in the data matrix after phase correction is obtained. The amplitude correction vector is obtained by energy normalization. The received data matrix after amplitude and phase error correction is expanded to generate a new covariance matrix. The covariance matrix is ​​processed in blocks, two new matrices are constructed and the eigenvalues ​​are calculated to realize DOA estimation of the parallel hydrophone array. Compared with the estimation results before the amplitude and phase errors are corrected, the method of this patent effectively improves the DOA estimation accuracy and estimation performance of the parallel hydrophone array. At the same time, the Euler transform is used to convert the complex number operations in the PM algorithm into real number operations, which reduces the calculation complexity and realizes automatic matching of two-dimensional angles. Summary of the invention

[0004] This patent aims at the shortcomings of the existing DOA estimation methods of parallel hydrophone arrays with amplitude and phase errors, and proposes a robust DOA estimation method based on parallel hydrophone arrays. Based on the properties of the covariance matrix, this patented method studies the auxiliary array element method based on the covariance matrix to correct the phase error of the array element. According to the data matrix after phase correction, the amplitude error of the array element is corrected by energy normalization. The received data matrix after amplitude and phase error correction is expanded, a new covariance matrix is ​​generated and processed in blocks, and two new matrices are constructed to calculate the eigenvalues ​​through eigendecomposition, so as to obtain accurate estimates of the azimuth and elevation angles. This patented method effectively improves the DOA estimation accuracy and estimation performance of the parallel hydrophone array, and reduces the computational complexity.

[0005] The parallel hydrophone array of the present invention is composed 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, d = λ / 2, λ represents the wavelength of the signal, and the first two array elements of subarray 1 and subarray 2 are respectively assumed to be auxiliary array elements. A robust DOA estimation method based on a parallel hydrophone array includes the following steps:

[0006] Step 1: The received data of subarray 1 and subarray 2 are represented as X(t) and Y(t), respectively, where t represents the sampling index; the received data of the two subarrays are Hilbert transformed to obtain the complex signal received by the array:

[0007] Step 2: Use the complex signal of sub-array 1 to calculate the covariance matrix of sub-array 1: Calculate the covariance matrix of sub-array 2 using the complex signal of sub-array 2 In the formula [·] H represents the conjugate transpose of the matrix, L represents the number of snapshots;

[0008] Step 3: Extract the covariance matrix R of sub-matrix 1 xx The phase information can be obtained as follows: Where angle(·) represents the phase angle of the calculated covariance matrix;

[0009] Step 4: Extract the matrix angle (R xx ), and use it to construct the upper triangular matrix T1.

[0010] Step 5: Since the first two array elements of subarray 1 are precisely calibrated auxiliary array elements, let h 12 = ψ1, subtract ψ1 from all elements in the first subdiagonal above the main diagonal of matrix T1, and we get α i =h i(i+1) -ψ1,i=1,2,...,M-1, construct the equation system The phase error vector is obtained This vector can be used to construct the phase error vector of subarray 1: The phase error vector can be used to obtain the phase correction matrix Multiply it by Obtain the received data matrix X'(t) after phase correction;

[0011] Step 6: According to the formula Calculate the average energy of each array element in the phase-corrected data matrix X'(t) i=1,2,...,M, through the formula Calculate the amplitude correction factor of the hydrophone and get the amplitude correction vector of subarray 1 Use this vector to generate the amplitude correction matrix Multiply it by X'(t) to get the received data matrix after amplitude and phase error correction of subarray 1

[0012] Step 7: Replace the covariance matrix in step 3 with R yyRepeat steps 3, 4, 5, and 6 to get the received data matrix after amplitude and phase error correction of subarray 2.

[0013] Step 8: Exploitation and structure Let x c (t) = [x(t) + x * (t)] / 2,x s (t) = [x(t) - x * (t)] / 2j, where (·) * Represents the conjugate of the matrix, construct Calculate the covariance matrix

[0014] Step 9: R r Divide into blocks to get R1 and R2, R r =[R1 R2], Further we get P = R1 -1 R2, construction

[0015] Step 10: Construct the matrix J1 = [I M I M 0 M 0 M ],J2=[0 M 0 M -I M I M ],

[0016]

[0017] Step 11: According to the constructed P c , J1, J2, J3, J4, calculate the matrices P1 and P2, where P1 = (J1P c ) + J2P c , P2=(J3P c ) + J4P c , perform eigendecomposition on P1 and obtain the kth diagonal element γ of the eigenvalue matrix k , perform eigendecomposition on P2 and obtain the kth diagonal element ρ of the eigenvalue matrix k , where k = 1, 2, ..., K, K represents the number of incident signals, where (·) + represents the generalized inverse of a matrix;

[0018] Step 12: According to the formula The estimated value of the pitch angle can be further calculated and the estimated value of the azimuth

[0019] Compared with the prior art, the technical solution of this patent has the following technical effects:

[0020] (1) The present invention effectively solves the influence of amplitude and phase errors in array elements on DOA estimation, and improves the DOA estimation accuracy of parallel sensor arrays when amplitude and phase errors exist.

[0021] (2) The method of the present invention can still obtain robust DOA estimation performance when the number of snapshots is small and the signal-to-noise ratio is low, has high angular resolution, and achieves matching of azimuth and elevation angles.

[0022] (3) The method of the present invention has strong robustness, can accurately estimate the azimuth information of the signal, does not require spectrum peak search, and has low computational complexity. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0025] Figure 3 It is the relationship curve between the success probability and the signal-to-noise ratio of the patented signal processing method;

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

[0027] Figure 5 is the relationship curve between the success probability and the number of snapshots of the patented signal processing method;

[0028] Figure 6 This is the relationship curve between the root mean square error and the number of snapshots of the patented signal processing method. DETAILED DESCRIPTION

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

[0030] First embodiment: Figure 1 The array model of the parallel hydrophone used in the present invention is given. The parallel hydrophone array consists of subarray 1 and subarray 2, each subarray has 6 array elements, that is, M=6, and the spacing between adjacent array elements and two parallel linear arrays in the array is 0.03m. Assume that the first two array elements of subarray 1 and subarray 2 are auxiliary array elements, the number of signals K=1, the azimuth angle is 10°, the elevation angle is 45°, the signal-to-noise ratio is 10dB, and the number of snapshots is L=500.

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

[0032] Step 1: The received data of subarray 1 and subarray 2 are represented as X(t) and Y(t), respectively, where t represents the sampling index; the received data of the two subarrays are Hilbert transformed to obtain the complex signal received by the array:

[0033] Step 2: Use the complex signal of sub-array 1 to calculate the covariance matrix of sub-array 1: Calculate the covariance matrix of sub-array 2 using the complex signal of sub-array 2 In the formula [·] H represents the conjugate transpose of the matrix, the number of snapshots L = 500;

[0034] Step 3: Extract the covariance matrix R of sub-matrix 1 xx The phase information can be obtained as follows: Where angle(·) represents the phase angle of the calculated covariance matrix;

[0035] Step 4: Extract the matrix angle (R xx ), and use it to construct the upper triangular matrix T1.

[0036] Step 5: Since the first two array elements of subarray 1 are precisely calibrated auxiliary array elements, let h 12 = ψ1, subtract ψ1 from all elements in the first subdiagonal above the main diagonal of matrix T1, and we get α i =h i(i+1) -ψ1,i=1,2,...,5, construct the equation system The phase error vector is obtained This vector can be used to construct the phase error vector of subarray 1: The phase error vector can be used to obtain the phase correction matrix Multiply it by Obtain the received data matrix X'(t) after phase correction;

[0037] Step 6: According to the formula Calculate the average energy of each array element in the phase-corrected data matrix X'(t) i=1,2,...,6, by the formula Calculate the amplitude correction factor of the hydrophone and get the amplitude correction vector of subarray 1 Use this vector to generate the amplitude correction matrix Multiply it by X'(t) to get the received data matrix after amplitude and phase error correction of subarray 1

[0038] Step 7: Replace the covariance matrix in step 3 with R yy Repeat steps 3, 4, 5, and 6 to get the received data matrix after amplitude and phase error correction of subarray 2.

[0039] Step 8: Exploitation and structure Let x c (t) = [x(t) + x * (t)] / 2,x s (t) = [x(t) - x * (t)] / 2j, where (·) * Represents the conjugate of the matrix, construct Calculate the covariance matrix

[0040] Step 9: R r Divide into blocks to get R1 and R2, R r =[R1 R2], Further we get P = R1 -1 R2, construction

[0041] Step 10: Construct the matrix J1 = [I6 I6 06 06], J2 = [06 06 -I6 I6],

[0042]

[0043] Step 11: According to the constructed P c , J1, J2, J3, J4, calculate the matrices P1 and P2, where P1 = (J1P c ) + J2P c , P2=(J3P c ) + J4P c , perform eigendecomposition on P1 and obtain the first diagonal element γ of the eigenvalue matrix k , perform eigendecomposition on P2 and obtain the first diagonal element ρ of the eigenvalue matrix k , where (·) + represents the generalized inverse of a matrix;

[0044] Step 12: According to the formula The estimated value of the pitch angle can be further calculated and the estimated value of the azimuth

[0045] The patent method based on the above conditions is simulated by MATLAB simulation software, and the angle estimation result of the patent method is obtained as follows Figure 2 As shown in the figure, when the input azimuth angle is 10° and the pitch angle is 45°, the estimated angles are 10.71° and 45.02°, respectively. However, when the amplitude and phase error correction is not performed, the estimated angles are 14.15° and 42.28°, respectively. The azimuth of the signal cannot be accurately estimated. Figure 2 It can be seen that the experimental results show that the patented method can accurately locate the signal source.

[0046] The second embodiment: the relationship curve between the success probability and the signal-to-noise ratio of the signal processing method of the present patent, such as Figure 3 As shown, a comparison chart of the success probability of estimating the signal direction is also given when there is an amplitude and phase error and other conditions remain unchanged, as the signal-to-noise ratio increases. The application conditions of the method of the present invention are as follows:

[0047] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1 As shown, the number of array elements in each subarray is a uniform linear array of M=6, the array element spacing is half wavelength, that is, d=0.03m, the incident angle of the source is 10° in azimuth and 45° in elevation, the number of snapshots is L=500, 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, and the obtained success probability curve is shown as follows Figure 3 shown.

[0048] Depend on Figure 3 It can be seen that the success probability of the angle estimation method of the patent and when there are amplitude and phase errors increases with the increase of the signal-to-noise ratio. Overall, the performance of the method of the present invention is better than the DOA estimation result when the amplitude and phase errors are not corrected, and it also has a strong estimation performance under low signal-to-noise ratio.

[0049] Example 3: The relationship curve between the root mean square error and the signal-to-noise ratio of the signal processing method of this patent, the effect diagram obtained is as follows Figure 4 As shown. At the same time, a comparison chart of the performance analysis of the two methods is given when there is an amplitude and phase error and other conditions remain unchanged as the signal-to-noise ratio increases. The application conditions of the method of the present invention are as follows:

[0050] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1As shown, the number of array elements in each subarray is a uniform linear array of M=6, the array element spacing is half wavelength, that is, d=0.03m, the incident angle of the source is 10° in azimuth and 45° in elevation, the number of snapshots is L=500, 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, and the performance analysis curve obtained is as follows Figure 4 shown.

[0051] Depend on Figure 4 It can be seen that the root mean square error of the method of this patent and when there are amplitude and phase errors decreases with the increase of the signal-to-noise ratio. From the overall effect, the performance of the method of the present invention is better than the DOA estimation result when the amplitude and phase errors are not corrected, and the azimuth information of the signal can be accurately estimated.

[0052] Example 4: The relationship between the success probability and the number of snapshots of the signal processing method of the present invention is shown in the following figure: Figure 5 As shown. At the same time, a comparison chart of the success probability of estimating the signal direction is given when there is an amplitude and phase error and other conditions remain unchanged, as the number of snapshots increases. The application conditions of the method of the present invention are as follows:

[0053] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1 As shown, the number of array elements in each subarray is a uniform linear array of M=6, the array element spacing is half wavelength, that is, d=0.03m, the incident angle of the source is 10° in azimuth and 45° in elevation, the signal-to-noise ratio is 10dB, the number of snapshots changes from 100, and increases from 100 to 1000 at a step size of 100 each time, 1000 independent Monte Carlo experiments are performed, and simulation analysis is performed using MATLAB simulation software, and the obtained success probability curve is shown as follows Figure 5 shown.

[0054] Depend on Figure 5 It can be seen that the success probability of the angle estimation method of the patent and when there are amplitude and phase errors increases with the increase of the number of snapshots. Overall, the performance of the method of the present invention is better than the DOA estimation result when the amplitude and phase errors are not corrected, and it also has a strong estimation performance at a low number of snapshots.

[0055] Example 5: The relationship curve between the root mean square error and the number of snapshots of the signal processing method of the present patent, the effect diagram obtained is as follows Figure 6 As shown. At the same time, a comparison chart of the performance analysis of the two methods is given when there is an amplitude and phase error and other conditions remain unchanged as the number of snapshots increases. The application conditions of the method of the present invention are as follows:

[0056] The method of the present invention is based on the array model of double parallel linear arrays. Figure 1As shown, the number of array elements in each subarray is a uniform linear array of M=6, the array element spacing is half wavelength, that is, d=0.03m, the incident angle of the source is 10° in azimuth and 45° in elevation, the signal-to-noise ratio is 10dB, the number of snapshots changes from 100, and increases from 100 to 1000 at a step size of 100 each time, 1000 independent Monte Carlo experiments are performed, and simulation analysis is performed using MATLAB simulation software, and the performance analysis curve obtained is shown as follows Figure 6 shown.

[0057] Depend on Figure 6 It can be seen that the root mean square error of the method of the patent and the case where amplitude and phase errors exist decreases with the increase of the number of snapshots. From the overall effect, the performance of the method of the present invention is better than the DOA estimation result when amplitude and phase errors are not corrected, and the azimuth information of the signal can be accurately estimated.

[0058] The specific examples described herein are merely illustrative of the present invention. A person skilled in the art may make various modifications, supplements or substitute similar methods to the specific examples described, but this will not deviate from the present invention or exceed the scope defined by the appended claims.

Claims

1. A robust DOA estimation method based on a parallel hydrophone array, characterized by: The parallel hydrophone array is composed 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, d=λ / 2, λ represents the wavelength of the signal, and the first two array elements of subarray 1 and subarray 2 are auxiliary array elements. A robust DOA estimation method based on the parallel hydrophone array includes the following steps: Step 1: The received data of subarray 1 and subarray 2 are represented as X(t) and Y(t), respectively, where t represents the sampling index; the received data of the two subarrays are Hilbert transformed to obtain the complex signal received by the array: Step 2: Use the complex signal of sub-array 1 to calculate the covariance matrix of sub-array 1: Calculate the covariance matrix of sub-array 2 using the complex signal of sub-array 2 In the formula [·] H represents the conjugate transpose of the matrix, L represents the number of snapshots; Step 3: Extract the covariance matrix R of sub-matrix 1 xx The phase information can be obtained as follows: Where angle(·) represents the phase angle of the calculated covariance matrix; Step 4: Extract the matrix angle (R xx ), and use it to construct the upper triangular matrix T1. Step 5: Since the first two array elements of subarray 1 are precisely calibrated auxiliary array elements, let h 12 = ψ1, subtract ψ1 from all elements in the first subdiagonal above the main diagonal of matrix T1, and we get α i =h i(i+1) -ψ1,i=1,2,...,M-1, construct the equation system The phase error vector is obtained This vector can be used to construct the phase error vector of subarray 1: The phase error vector can be used to obtain the phase correction matrix Multiply it by Obtain the received data matrix X'(t) after phase correction; Step 6: According to the formula Calculate the average energy of each array element in the phase-corrected data matrix X'(t) i=1,2,...,M, by the formula Calculate the amplitude correction factor of the hydrophone and get the amplitude correction vector of subarray 1 Use this vector to generate the amplitude correction matrix Multiply it by X'(t) to get the received data matrix after amplitude and phase error correction of subarray 1 Step 7: Replace the covariance matrix in step 3 with R yy Repeat steps 3, 4, 5, and 6 to get the received data matrix after amplitude and phase error correction of subarray 2. Step 8: Exploitation and structure Let x c (t) = [x(t) + x * (t)] / 2, xs(t)=[x(t)-x*(t)] / 2j, where (·) * Represents the conjugate of the matrix, construct Calculate the covariance matrix Step 9: R r Divide into blocks to get R1 and R2, R r =[R1 R2], Further we get P = R1 -1 R2, construction Step 10: Construct the matrix J1 = [I M I M 0 M 0 M ],J2=[0 M 0 M -I M I M ], Step 11: Based on the constructed Pc, J1, J2, J3, J4, calculate the matrices P1 and P2, where P1 = (J1P c ) + J2P c , P2=(J3P c ) + J4P c , perform eigendecomposition on P1 and obtain the kth diagonal element γ of the eigenvalue matrix k , perform eigendecomposition on P2 and obtain the kth diagonal element ρ of the eigenvalue matrix k , where k = 1, 2, ..., K, K represents the number of incident signals, where (·) + represents the generalized inverse of a matrix; Step 12: According to the formula Two intermediate parameters can be estimated, and the estimated value of the pitch angle can be further calculated using these two intermediate parameters and the estimated value of the azimuth

Images