Robust hydrophone array orientation estimation method

Abstract

The invention discloses a robust hydrophone array orientation estimation method. According to the scheme, amplitude-phase error correction is introduced, and phase error vector and energy normalization are utilized to realize amplitude-phase error correction. The method comprises the following steps: correcting data, reconstructing a covariance matrix by using the corrected data, calculating a spatial spectrum by using the reconstructed covariance matrix, then calculating a first-order difference function of a spectral function, and performing spectral peak search on a differentiated spatial spectrum curve to estimate the direction of an incident signal. Simulation experiment and lake test experiment results show that the method greatly reduces the influence of amplitude-phase errors on DOA estimation performance, improves the orientation estimation precision, and is low in calculation complexity. The scheme of the patent has important application value in ocean instruments such as subsurface buoys and ocean mobile observation platforms, and has high practical engineering application value.

Description

Technical Field

[0001] The present invention belongs to the technical field of marine acoustic measurement and detection equipment, and is a robust hydrophone array azimuth estimation method. Specifically, it is a method that utilizes phase error vector and energy normalization to realize the correction of amplitude and phase errors, and then utilizes the corrected data to reconstruct the covariance matrix, performs eigendecomposition on it, and utilizes the noise subspace to construct the first-order difference function of the spectral function to perform DOA estimation. Background Art

[0002] The total area of ​​the ocean is about 360 million square kilometers, with an average water depth of about 3,795 meters. Only 5% of the seabed has been explored by humans today. In order to make full use of the rich marine resources, more and more marine-specific instruments are used to collect, process and transmit underwater information. Among them, because the buoy can detect a variety of data, it can realize multi-parameter synchronous detection, and can be placed in the ocean for a long time to continuously obtain data. It is widely used in mobile observation platforms, offshore platform-based observation platforms and other scenarios, and plays an important role in marine scientific research. However, the corrosion of seawater and the attachment of marine plankton can cause damage to the buoy, which in turn causes errors and affects the detection accuracy. In this context, it is of great significance to study robust target parameter estimation methods and apply them to marine-specific instruments such as buoys and marine mobile observation platforms.

[0003] In array signal processing, direction of arrival estimation is an important branch and is widely used in many fields such as radar, sonar, wireless communication, underwater target detection and biomedical engineering. The uniform linear array has a simple structure and is a widely used traditional array structure. Traditional DOA estimation algorithms mainly include the MUSIC method and ESPRIT method based on eigendecomposition. In most cases, these classic algorithms can basically achieve a high degree of accuracy and meet daily needs. In actual engineering applications, submerged buoys that have been in seawater for a long time will produce amplitude and phase errors due to seawater erosion or plankton attachment, which will affect the array flow pattern, resulting in a significant decrease in DOA estimation performance or even complete failure.

[0004] There are two types of solutions to the above problems: active correction algorithm and self-correction algorithm. The active correction algorithm requires an auxiliary source with known azimuth information. Although the amount of calculation is small, the accuracy of the auxiliary source is high. The self-correction algorithm uses an optimization function to achieve joint estimation of error parameters and DOA, but the amount of calculation is large and the accuracy cannot be guaranteed. Currently, commonly used amplitude and phase correction algorithms, such as the ISM algorithm, are only applicable to situations where the auxiliary array elements are adjacent to each other, and the effect is not good in actual engineering applications.

[0005] This patent proposes a robust uniform linear array DOA estimation method. The amplitude and phase error correction algorithm uses the phase error vector and energy normalization to achieve the correction of the amplitude and phase error. The covariance matrix is ​​reconstructed using the corrected data, and the first-order difference function of the spectrum function is calculated to achieve DOA estimation. While the amount of calculation is not significantly increased, the estimation accuracy is improved, the impact of the amplitude and phase error is greatly reduced, and the placement of auxiliary array elements is not limited to special positions, which has higher practical engineering application value. Summary of the invention

[0006] This patent proposes a robust hydrophone array azimuth estimation method to address the problem that the estimation performance of the traditional uniform hydrophone array azimuth estimation method is significantly reduced due to the influence of amplitude and phase errors in actual engineering. The characteristics are: the patent method uses energy normalization and phase error vector to realize the correction of amplitude and phase errors, and then reconstructs the covariance matrix using the corrected data, performs eigendecomposition on it to obtain the noise subspace, constructs the first-order difference function of the spectral function as a new "spectral function", and performs DOA estimation through spectral peak search. The patent method has low computational complexity, flexible auxiliary array element settings, can significantly improve estimation accuracy, and has important engineering value in practical applications.

[0007] In the technical solution of the present invention, the buoy detector uses a uniform linear array composed of M hydrophones to detect targets, and the array element spacing of the detection system is Where λ is the wavelength of the signal; the received data of the M hydrophone arrays is recorded as x(t)=[x1(t),x2(t),…,x M (t)] T , the steps of the azimuth estimation method of the hydrophone array are as follows:

[0008] Step 1: Perform Hilbert transform on the received data of M hydrophones to convert them into complex signals. The complex signal form of the hydrophone array received signal x(t) is: in represents the array flow matrix when there is amplitude and phase error, t=1,2,…,L represents the sampling sequence number, represents the incident signal vector, represents the complex noise vector of the hydrophone array;

[0009] Step 2: According to the formula Calculate the covariance matrix of the complex signal received by the hydrophone array, where L represents the number of snapshots, [·] H represents the conjugate transpose of a matrix;

[0010] Step 3: Assume that the P array elements of the uniform linear array are auxiliary array elements for precise correction, and P ≥ K + 1, K represents the number of incident signals, and extract the covariance matrix of subarray 1 Phase information Where angle(·) represents the extracted phase information; the phase error is calculated using the upper triangular elements of the phase information matrix. The upper triangular matrix composed of the upper triangular elements of the phase information matrix is ​​recorded as in Represents the phase error. Since it is assumed that the first P array elements are accurately calibrated auxiliary array elements, let m represents the row number, n represents the column number;

[0011] Step 4: Extract the first P-1 elements of the first diagonal line next to the main diagonal and because p=1,2,…,P, we can get At the same time, for the first diagonal element, each element plus And order i=1,2,…,M-1, then, Combining linear algebra, we get because The above formula can be simplified to The phase error vector is obtained from the above equation:

[0012] Step 5: Multiply the received data of the M hydrophones of the buoy detector by the corresponding phase correction factor. The correction factor is in the form of To achieve phase correction, the corrected data is recorded as

[0013] Step 6: According to the calculation formula of discrete signal energy Calculate the phase-corrected data The average energy per hydrophone By formula Get the amplitude correction vector

[0014] Step 7: Phase-corrected data from the M hydrophones of the buoy detector Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data can be expressed as X(t); calculate the covariance matrix of the received data matrix after removing the amplitude and phase errors

[0015] Step 8: Make θ uniformly distributed in the range of [-π / 2,π / 2], and use the covariance matrix estimate in step 7 Calculate the spatial spectrum P(θ);

[0016] Step 9: According to the formula

[0017]

[0018] Calculate the first-order difference function y(θ) of the spectrum function P(θ), use it as the new spatial spectrum, take the absolute value of y(θ), find the peak, average the two peaks with the smallest difference, and get the estimated value of the incident angle of the incident signal.

[0019] The calculation method of the spatial spectrum P(θ) in step eight is as follows:

[0020] A1: Covariance matrix Perform eigenvalue decomposition and sort the eigenvalues ​​from large to small, dividing them into signal subspace and noise subspace. N It is composed of the eigenvectors corresponding to the K+1th to Mth eigenvalues;

[0021] A2: Let θ take values ​​uniformly in the range of [-π / 2,π / 2], and calculate the spatial spectrum of the signal according to the formula:

[0022] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0023] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0024] The calculation formula of the spatial spectrum P(θ) in step 8 is: where ||·||2 represents the 2-norm of the vector.

[0025] The calculation formula of the spatial spectrum P(θ) in step 8 is: The eigenvalue λ i Sort from largest to smallest, q i is i The corresponding feature vector.

[0026] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0027] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0028] 1. It can overcome the influence of amplitude and phase errors of array element receiving data, provide good estimation performance, and have higher practical engineering application value.

[0029] 2. Low computational complexity: In the array receiving data, the phase error vector and amplitude error vector are extracted by extracting phase information and energy normalization to correct the amplitude and phase errors, which reduces the computational complexity and saves time costs.

[0030] 3. Compared with the ISM method, the selection of auxiliary array elements is not limited to the front end of the array element, the position selection is more flexible, and the estimation effect is better in the lake test experiment, which has higher practical application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 It is the normalized spatial spectrum of the signal processing method of this patent;

[0032] Figure 2 It is the relationship curve between the success probability and the average deviation of amplitude and phase errors of the patented signal processing method;

[0033] Figure 3 It is the relationship curve between the root mean square error and the average deviation of the amplitude and phase errors of the signal processing method of this patent;

[0034] Figure 4 This is the relationship curve between the success probability and the signal-to-noise ratio of the patented signal processing method;

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

[0036] Figure 6 is a curve showing the relationship between the success probability and the number of snapshots of the patented signal processing method;

[0037] Figure 7 This is the spatial spectrum process diagram of the lake test experiment of the patented signal processing method; DETAILED DESCRIPTION

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

[0039] Example 1: A uniform linear array with M=4 elements, P=2 auxiliary elements, and N=M+P=6 total elements is used. An incoherent signal is incident on the array from a direction of 15°. The element spacing d=0.03, the number of snapshots L=400, the signal-to-noise ratio is 10dB, and the phase error obeys the variance The normal distribution has a mean of 0, and the amplitude error follows a uniform distribution of [-0.6, 0.6]. Using the above conditions, the specific implementation process is:

[0040] Step 1: Perform Hilbert transform on the received data of M hydrophones to convert them into complex signals. The complex signal form of the hydrophone array received signal x(t) is: in represents the array flow matrix when there is amplitude and phase error, t=1,2,…,L represents the sampling sequence number, represents the incident signal vector, represents the complex noise vector of the hydrophone array;

[0041] Step 2: According to the formula Calculate the covariance matrix of the complex signal received by the hydrophone array, where L represents the number of snapshots, [·] H represents the conjugate transpose of a matrix;

[0042] Step 3: Assume that the P array elements of the uniform linear array are auxiliary array elements for precise correction, and P ≥ K + 1, K represents the number of incident signals, and extract the covariance matrix of subarray 1 Phase information Where angle(·) represents the extracted phase information; the phase error is calculated using the upper triangular elements of the phase information matrix. The upper triangular matrix composed of the upper triangular elements of the phase information matrix is ​​recorded as in Represents the phase error. Since it is assumed that the first P array elements are accurately calibrated auxiliary array elements, let m represents the row number, n represents the column number;

[0043] Step 4: Extract the first P-1 elements of the first diagonal line next to the main diagonal and because p=1,2,…,P, we can get At the same time, for the first diagonal element, each element plus And order i=1,2,…,M-1, then, Combining linear algebra, we get because The above formula can be simplified to The phase error vector is obtained from the above equation:

[0044] Step 5: Multiply the received data of the M hydrophones of the buoy detector by the corresponding phase correction factor. The correction factor is in the form of To achieve phase correction, the corrected data is recorded as

[0045] Step 6: According to the calculation formula of discrete signal energy Calculate the phase-corrected data The average energy per hydrophone By formula Get the amplitude correction vector

[0046] Step 7: Phase-corrected data from the M hydrophones of the buoy detector Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data can be expressed as X(t); calculate the covariance matrix of the received data matrix after removing the amplitude and phase errors

[0047] Step 8: Make θ uniformly distributed in the range of [-π / 2,π / 2], and use the covariance matrix estimate in step 7 Calculate the spatial spectrum P(θ);

[0048] Step 9: According to the formula

[0049]

[0050] Calculate the first-order difference function y(θ) of the spectrum function P(θ), use it as the new spatial spectrum, take the absolute value of y(θ), find the peak, average the two peaks with the smallest difference, and get the estimated value of the incident angle of the incident signal.

[0051] The calculation method of the spatial spectrum P(θ) in step eight is as follows:

[0052] A1: Covariance matrix Perform eigenvalue decomposition and sort the eigenvalues ​​from large to small, dividing them into signal subspace and noise subspace. N It is composed of the eigenvectors corresponding to the K+1th to Mth eigenvalues;

[0053] A2: Let θ take values ​​uniformly in the range of [-π / 2,π / 2], and calculate the spatial spectrum of the signal according to the formula:

[0054] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0055] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0056] The calculation formula of the spatial spectrum P(θ) in step 8 is: where ||·||2 represents the 2-norm of the vector.

[0057] The calculation formula of the spatial spectrum P(θ) in step 8 is: The eigenvalue λ i Sort from largest to smallest, q i is i The corresponding feature vector.

[0058] The calculation formula of the spatial spectrum P(θ) in step 8 is:

[0059] The patent method based on the above conditions is simulated by MATLAB simulation software, and the angle estimation results of the patent method, ISM method and MUSIC method are as follows: Figure 1 As shown, analysis Figure 1It can be seen that the angles estimated by the three methods are all around 15°, but the energy spectrum of the MUSIC method is wider.

[0060] Second embodiment: The relationship between the success probability of the signal processing method of the present invention and the average deviation of the amplitude and phase errors and the root mean square error (RMSE) and the average deviation of the amplitude and phase errors is studied, and the effect diagram is as follows: Figure 2 , Figure 3 As shown; the performance analysis comparison diagram of the MUSIC method is also given. The conditions for the application of the algorithm in the present invention are as follows: a uniform linear array with an array element number M = 4, an auxiliary array element number P = 2, and a total array element number N = M + P = 6 is used, an incoherent signal is incident on the array from a direction of 15°, the array element spacing d = 0.03, the number of snapshots L = 400, the signal-to-noise ratio is 10dB, the average deviation of the amplitude and phase error increases from 0.07 to 0.21 with a step size of 0.02, and 500 independent Monte Carlo experiments are performed.

[0061] analyze Figure 2 It can be seen that with the increase of the average deviation of the amplitude and phase errors, the success probability of the patented method fluctuates around 94%, maintaining good estimation performance; the success probability of the MUSIC method continues to decrease, from 84.25% when the average deviation of the amplitude and phase errors is 0.07 to 39% when the average deviation of the amplitude and phase errors is 0.21.

[0062] analyze Figure 3 It can be seen that with the increase of the average bias of the amplitude and phase errors, the RMSE of the patented method fluctuates around 0.259°, maintaining good estimation performance; the success probability of the MUSIC method continues to decrease, from 0.337° when the average bias of the amplitude and phase errors is 0.07 to 0.966° when the average bias of the amplitude and phase errors is 0.21.

[0063] The third embodiment: The relationship between the success probability and the signal-to-noise ratio and the relationship between the root mean square error (RMSE) and the signal-to-noise ratio of the signal processing method of the patent are studied respectively, and the effect diagram is as follows: Figure 4 , Figure 5 As shown, the performance analysis comparison diagram of the ISM method is also given. The application conditions of the algorithm in the present invention are as follows: a uniform linear array with an array element number M = 4, an auxiliary array element number P = 2, and a total array element number N = M + P = 6 is used. An incoherent signal is incident on the array from a direction of 15°, the array element spacing d = 0.03, the number of snapshots L = 400, the signal-to-noise ratio is changed from -5dB to 20dB in steps of 5dB, and the phase error obeys the variance The normal distribution with mean 0 and the amplitude error obey the uniform distribution of [-0.6, 0.6] were used, and 500 independent Monte Carlo experiments were performed.

[0064] analyze Figure 4It can be seen that when the number of snapshots L = 400, the success probability of the two amplitude and phase error correction algorithms increases with the increase of the signal-to-noise ratio. The success probability of the patented method is always higher than that of the ISM method. When the signal-to-noise ratio is 5dB, the success probability of the patented method is 3% higher than that of the ISM method. When the signal-to-noise ratio is 10dB, the success probability of the resolution of the patented method reaches 100%.

[0065] analyze Figure 5 It can be seen that when the snapshot number L = 400, as the signal-to-noise ratio increases, the RMSE of the two amplitude and phase error correction algorithms continues to decrease with the increase of the signal-to-noise ratio. The RMSE of the patented method is always lower than that of the ISM method. When the signal-to-noise ratio is 0dB, the RMSE of the patented method is 0.14° lower than that of the ISM method.

[0066] The fourth embodiment: The relationship between the root mean square error (RMSE) and the number of snapshots of the signal processing method of the present patent is studied respectively, and the effect diagram is as follows: Figure 6 The application conditions of the algorithm in the present invention are as follows: a uniform linear array with M=4 array elements, P=2 auxiliary array elements, and N=M+P=6 total array elements is used. An incoherent signal is incident on the array from a direction of 15°. The array element spacing is d=0.03. The signal-to-noise ratio is 10dB. The number of snapshots is changed from L=50 and increases from 50 to 400 in steps of 50. The phase error obeys the variance The normal distribution with mean 0 and the amplitude error obey the uniform distribution of [-0.6, 0.6] were used, and 500 independent Monte Carlo experiments were performed.

[0067] analyze Figure 6 It can be seen that when the signal-to-noise ratio is 10dB, as the number of snapshots increases, the RMSE of the two amplitude and phase error correction algorithms decreases slightly when the value is small, and the RMSE of the patented method is smaller than that of the ISM method most of the time. When the number of snapshots is 100, the RMSE of the patented method is 0.18° lower than that of the ISM method.

[0068] Embodiment 5: When there is an amplitude and phase error in the lake test, the spatial spectrum history diagram obtained by the signal processing method of this patent is as follows Figure 7 As shown in the figure, the spatial spectrum process diagram of the ISM method is also given. The application conditions of the method of the present invention are as follows: the number of array elements M = 4, the number of auxiliary array elements P = 2. A uniform linear array with a total number of array elements N = M + P = 6, the auxiliary array elements are located at the front end, the far-field narrowband CW pulse signal is incident on the array from the direction of 18°, the source center frequency is 30kHz, the array element spacing d = 0.03, the signal is divided into 50 groups, and the number of snapshots in each group is L = 1000.

[0069] analyze Figure 7It can be seen that in the presence of amplitude and phase errors, the power of the patented method is almost entirely concentrated near the incident angle of the incident signal, and it can be clearly seen that the angle changes over time in a straight line; the energy of the ISM method is mainly concentrated between 9° and 26°, and the energy of the MUSIC algorithm is mainly concentrated between 6° and 32°. The actual incident angle of the incident signal cannot be observed from the history diagram. Therefore, the patented method can overcome the impact of amplitude and phase errors in actual engineering applications.

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

Claims

1. A robust hydrophone array azimuth estimation method, characterized by: The buoy detector uses a uniform linear array composed of M hydrophones to detect targets. The array element spacing of the detection system is Where λ is the wavelength of the signal; the received data of the M hydrophone arrays is recorded as x(t)=[x1(t),x2(t),…,x M (t)] T , the steps of the azimuth estimation method of the hydrophone array are as follows: Step 1: Perform Hilbert transform on the received data of M hydrophones to convert them into complex signals. The complex signal form of the hydrophone array received signal x(t) is: in represents the array flow matrix when there is amplitude and phase error, t=1,2,…,L represents the sampling sequence number, represents the incident signal vector, represents the complex noise vector of the hydrophone array; Step 2: According to the formula Calculate the covariance matrix of the complex signal received by the hydrophone array, where L represents the number of snapshots, [·] H represents the conjugate transpose of a matrix; Step 3: Assume that the P array elements of the uniform linear array are auxiliary array elements for precise correction, and P ≥ K + 1, K represents the number of incident signals, and extract the covariance matrix of subarray 1 Phase information Where angle(·) represents the extracted phase information; the phase error is calculated using the upper triangular elements of the phase information matrix. The upper triangular matrix composed of the upper triangular elements of the phase information matrix is ​​recorded as in Represents the phase error. Since it is assumed that the first P array elements are accurately calibrated auxiliary array elements, let m represents the row number, n represents the column number; Step 4: Extract the first P-1 elements of the first diagonal line next to the main diagonal and because Available At the same time, for the first diagonal element, each element plus And order Then, by Combining linear algebra, we get because The above formula can be simplified to The phase error vector is obtained from the above equation: Step 5: Multiply the received data of the M hydrophones of the buoy detector by the corresponding phase correction factor. The correction factor is in the form of To achieve phase correction, the corrected data is recorded as Step 6: According to the calculation formula of discrete signal energy Calculate the phase-corrected data The average energy per hydrophone By formula Get the amplitude correction vector Step 7: Phase-corrected data from the M hydrophones of the buoy detector Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data can be expressed as X(t); calculate the covariance matrix of the received data matrix after removing the amplitude and phase errors Step 8: Make θ uniformly distributed in the range of [-π / 2,π / 2], and use the covariance matrix estimate in step 7 Calculate the spatial spectrum P(θ); Step 9: According to the formula Calculate the first-order difference function y(θ) of the spectrum function P(θ), use it as the new spatial spectrum, take the absolute value of y(θ), find the peak, average the two peaks with the smallest difference, and get the estimated value of the incident angle of the incident signal.

2. A robust hydrophone array azimuth estimation method as claimed in claim 1, characterized in that: The calculation method of spatial spectrum P(θ) is as follows: A1: Covariance matrix Perform eigenvalue decomposition and sort the eigenvalues ​​from large to small, dividing them into signal subspace and noise subspace. N It is composed of the eigenvectors corresponding to the K+1th to Mth eigenvalues; A2: Let θ take values ​​uniformly in the range of [-π / 2,π / 2], and calculate the spatial spectrum of the signal according to the formula:

3. A robust hydrophone array azimuth estimation method as claimed in claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

4. The method for robust hydrophone array azimuth estimation according to claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

5. The robust hydrophone array azimuth estimation method according to claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is: where ||·||2 represents the 2-norm of the vector.

6. A robust hydrophone array azimuth estimation method as claimed in claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is: The eigenvalue λ i Sort from largest to smallest, q i is i The corresponding feature vector.

7. The robust hydrophone array azimuth estimation method according to claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

Images