A robust hydrophone array direction of arrival estimation method

Abstract

The application discloses a kind of robust hydrophone array azimuth estimation method.This patent scheme introduces amplitude-phase error correction, utilizes phase error vector and energy normalization, realizes the correction of amplitude-phase error.Covariance matrix is reconstructed using the corrected data, and spatial spectrum is calculated using the reconstructed covariance matrix, then the first order differential function of spectrum function is calculated, and the spectrum peak search is carried out on the differential spatial spectrum curve to estimate the incident signal direction of arrival.The simulation experiment and lake test experiment results show that the method greatly reduces the influence of amplitude-phase error on DOA estimation performance, improves the estimation accuracy of direction, and has low computational complexity.This patent scheme has important application value in submarine and ocean mobile observation platform and other ocean instruments, and has high practical engineering application value.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of ocean acoustic measurement and detection equipment, and is a robust hydrophone array direction of arrival (DOA) estimation method. Specifically, the present application is a method for correcting amplitude and phase errors using phase error vector and energy normalization, reconstructing the covariance matrix using the corrected data, performing eigenvalue decomposition on the covariance matrix, using the noise subspace to construct the first-order difference function of the spectral function, and performing DOA estimation. BACKGROUND

[0002] The total area of the ocean is about 360 million square kilometers, with an average depth of about 3795 meters. Only 5% of the ocean floor has been explored by humans today. In order to fully utilize the rich marine resources, more and more marine special instruments are used for underwater information collection, processing and transmission. Among them, since the submersible buoy can detect multiple data, it can realize multi-parameter synchronous detection, and can be placed in the ocean for a long time to continuously acquire data. It is widely used in mobile observation platforms, offshore platform-based observation stations, etc., and plays an important role in marine scientific research. However, the corrosion of seawater and the attachment of marine plankton can cause the submersible buoy to be damaged, thereby producing errors and affecting the detection accuracy. Under this background, it is of great significance to study a robust target parameter estimation method and apply it to marine special instruments such as submersible buoys and marine mobile observation platforms.

[0003] In array signal processing, direction of arrival (DOA) estimation is an important branch and is widely used in radar, sonar, wireless communication, underwater target detection, and biomedical engineering. The uniform linear array (ULA) is simple in structure and is a widely used traditional array structure. Traditional DOA estimation algorithms mainly include the multiple signal classification (MUSIC) method and the estimation of signal parameters via rotational invariance techniques (ESPRIT) method. In most cases, these classic algorithms can basically achieve high accuracy and meet daily needs. In practical engineering applications, the submersible buoy in seawater for a long time will produce amplitude and phase errors due to seawater erosion or plankton attachment, thereby affecting the array flow pattern and causing the DOA estimation performance to decrease significantly or even fail completely.

[0004] To solve the above problems, there are two types of solutions: active correction algorithm and self-correction algorithm. The active correction algorithm requires an auxiliary source with known direction information, which has a small amount of calculation but requires high accuracy of the auxiliary source. The self-correction algorithm uses an optimization function to realize joint estimation of error parameters and DOA, but the amount of calculation is large and the accuracy cannot be guaranteed. The amplitude and phase correction algorithm commonly used at present, such as the iterative signal modeling (ISM) algorithm, is only applicable to the case where the positions of the auxiliary array elements are adjacent, and the effect is not good in practical engineering applications.

[0005] This patent proposes a robust DOA estimation method for uniform linear arrays. This amplitude and phase error correction algorithm utilizes the phase error vector and energy normalization to correct amplitude and phase errors. The covariance matrix is ​​reconstructed using the corrected data, and the first-order difference function of the spectral function is calculated to achieve DOA estimation. While not significantly increasing computational complexity, it improves estimation accuracy, greatly reduces the impact of amplitude and phase errors, and the placement of auxiliary array elements is not limited to specific locations, thus possessing greater practical engineering application value. Summary of the Invention

[0006] This patent addresses the problem of significantly degraded estimation performance of traditional uniform hydrophone array azimuth estimation methods in practical engineering due to amplitude and phase errors. It proposes a robust hydrophone array azimuth estimation method. The key features are: the method utilizes energy normalization and phase error vectors to correct amplitude and phase errors; then, it reconstructs the covariance matrix using the corrected data and performs eigenvalue decomposition to obtain the noise subspace; a first-order difference function of the spectral function is constructed as a new "spectral function"; and DOA estimation is performed through spectral peak search. This patented method has low computational complexity, flexible auxiliary array element settings, and can significantly improve estimation accuracy, possessing significant engineering value in practical applications.

[0007] In the technical solution of this invention, the underwater glider detector uses a uniform linear array of M hydrophones for target detection, and the element spacing of the detection system is... Where λ is the wavelength of the signal; the received data from the M hydrophone array is denoted as x(t) = [x1(t), x2(t), ..., x M (t)] T The steps for azimuth estimation of a hydrophone array are as follows:

[0008] Step 1: Perform a Hilbert transform on the received data from the M hydrophones to convert it into a complex signal. The complex signal form of the hydrophone array received signal x(t) is as follows: in This represents the array manifold matrix when amplitude and phase errors exist, where t = 1, 2, ..., L represents the sampling sequence number. Represents the incident signal vector. Represents the noise complex 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 elements of the uniform linear array are auxiliary elements for precise correction, and P ≥ K + 1, where K represents the number of incident signals. Extract the subarray covariance matrix. Phase information where angle(·) represents extracting phase information; phase errors are calculated by using upper triangular elements of the phase information matrix, and an upper triangular matrix composed of the upper triangular elements of the phase information matrix is denoted as where represents a phase error, and since the first P elements are assumed to be accurately corrected auxiliary elements, let m represents a row index, and n represents a column index;

[0011] Step four: extracting the first P-1 elements of the first diagonal line adjacent to the main diagonal and Since p=1, 2, …, P, the following can be obtained At the same time, for the elements of the first diagonal line, each element is added by and let i=1, 2, …, M-1, then, by combining linear algebra, the following can be obtained Since the above formula can be simplified as the phase error vector is solved from the above formula

[0012] Step five: multiplying the received data of the M hydrophones of the submarine marker detector by the corresponding phase correction factor, and the form of the correction factor is to realize the correction of the phase, and the corrected data is denoted as

[0013] Step six: calculating the data after the phase correction according to the calculation formula of the discrete signal energy the average energy of each hydrophone in the data after the phase correction is obtained by the formula the amplitude correction vector is obtained

[0014] Step seven: multiplying the data of the M hydrophones of the submarine marker detector after the phase correction by the corresponding amplitude correction factor to realize the correction of the amplitude; the corrected data can be represented as X(t); the covariance matrix of the received data matrix after removing the amplitude and phase errors is calculated

[0015] Step eight: making θ take values uniformly in the range of [-π / 2, π / 2], and using the covariance matrix estimation value in step seven the spatial spectrum P(θ) is calculated;

[0016] Step nine: according to the formula

[0017]

[0018] ​​The first difference function y (theta) of the spectrum function P (theta) is calculated, and the new space spectrum is obtained.

[0019] The calculation method of the space spectrum P (theta) in step eight is calculated as follows:

[0020] A1: eigenvalue decomposition is performed on the covariance matrix The eigenvalues are sorted in descending order, and the signal subspace and noise subspace are divided into two parts, wherein the noise subspace U N The eigenvectors corresponding to the first K+1 to M eigenvalues are composed of the eigenvectors corresponding to the first K+1 to M eigenvalues;

[0021] A2: make theta in the range of [-pi / 2, pi / 2] uniform, and calculate the spatial spectrum of the signal according to the formula

[0022] The calculation formula of the space spectrum P (theta) in step eight is:

[0023] The calculation formula of the space spectrum P (theta) in step eight is:

[0024] The calculation formula of the space spectrum P (theta) in step eight is: Where ||·||2 represents the 2-norm of the vector.

[0025] The calculation formula of the space spectrum P (theta) in step eight is: The eigenvalues lambda i are sorted in descending order, q i is the eigenvector corresponding to lambda i .

[0026] The calculation formula of the space spectrum P (theta) in step eight is:

[0027] Compared with the prior art, the above technical scheme has the following technical effects:

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

[0029] 2. Low calculation complexity: in array receiving data, the phase error vector and amplitude error vector are extracted by extracting phase information and energy normalization, etc. The calculation complexity is angle, and the time cost is saved.

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

[0031] Figure 1 The normalized spatial spectrum diagram of the signal processing method of the patent;

[0032] Figure 2 The relationship curve between the success probability of the signal processing method of the patent and the average bias value of the amplitude and phase error;

[0033] Figure 3 The relationship curve between the root mean square error of the signal processing method of the patent and the average bias value of the amplitude and phase error;

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

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

[0036] Figure 6 The relationship curve between the success probability of the signal processing method of the patent and the number of snapshots;

[0037] Figure 7 The spatial spectrum history diagram of the lake test experiment of the signal processing method of the patent; DETAILED DESCRIPTION

[0038] The present application will be further described in conjunction with embodiments, drawings:

[0039] The first embodiment: a uniform linear array with M=4 array elements, P=2 auxiliary array elements, and N=M+P=6 total array elements is used, a non-coherent 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 phase error is subject to a normal distribution with a mean of 0, and the amplitude error is subject to a uniform distribution of [-0.6, 0.6]. With the above conditions, the specific implementation process is as follows:

[0040] Step one: Hilbert transform is performed on the received data of the M hydrophones to convert them into complex signals, and the complex signal form of the hydrophone array received signal x(t) is wherein represents the array flow matrix in the presence of amplitude and phase errors, t=1, 2, …, L represents the sampling sequence number, represents the incident signal vector, represents the noise complex 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 elements of the uniform linear array are auxiliary elements for precise correction, and P ≥ K + 1, where K represents the number of incident signals. Extract the subarray covariance matrix. Phase information Where angle(·) represents the extraction of phase information; the phase error is calculated using the upper triangular elements of the phase information matrix, and the upper triangular matrix formed by the upper triangular elements of the phase information matrix is ​​denoted as . in Let represent the phase error. Since it is assumed that the first P elements are precisely corrected auxiliary elements, let . m represents the row number, and n represents the column number;

[0043] Step 4: Extract the first P-1 elements of the first diagonal line immediately adjacent to the main diagonal. because Given p = 1, 2, ..., P, we can obtain... At the same time, for each element on the first diagonal, add And order If i = 1, 2, ..., M-1, then, by Combining linear algebra, we get because The above formula can be simplified to The phase error vector can be obtained from the above equation.

[0044] Step 5: Multiply the received data from the M hydrophones of the underwater glider by the corresponding phase correction factor. The correction factor is in the form of... Phase correction is achieved, and the corrected data is denoted as...

[0045] Step Six: Calculate the energy of the discrete signal using the formula... Calculate the phase-corrected data Average energy of each hydrophone Through formula Obtain the amplitude correction vector

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

[0047] Step eight: make θ take value uniformly in [-π / 2, π / 2], and use the covariance matrix estimation value in step seven Calculate the spatial spectrum P(θ);

[0048] Step nine: according to the formula

[0049]

[0050] Calculate the first difference function y(θ) of the spectrum function P(θ), take it as a new spatial spectrum, take the absolute value of y(θ), find the peak value, and take the average of the two peak values with the smallest difference to obtain the estimation value of the incident angle of the incident signal.

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

[0052] A1: perform eigenvalue decomposition on the covariance matrix , sort it in descending order of eigenvalue, and divide it into signal subspace and noise subspace two parts, wherein the noise subspace U N is composed of the eigenvectors corresponding to the K+1 to M eigenvalues;

[0053] A2: make θ take value uniformly in [-π / 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 eight is:

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

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

[0057] The calculation formula of the spatial spectrum P(θ) in step eight is: Sort the eigenvalues λ i in descending order, q i is the eigenvector corresponding to λ i .

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

[0059] The simulation of the method of the patent based on the above conditions is carried out by MATLAB simulation software, and the angle estimation results of the method of the patent, ISM method and MUSIC method are as shown in Figure 1 , and the analysis Figure 1So, the estimated angles of the three methods are all near 15°, but the energy spectrum of the MUSIC method is wider.

[0060] The second embodiment: study the relationship between the success probability of the signal processing method of the patent and the average bias and root mean square error (RMSE) of the amplitude and phase error, and obtain the effect diagram as shown in Figure 2 、 Figure 3 The performance analysis comparison diagram of the MUSIC method is also given. The conditions for applying the algorithm in the application 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, a non-coherent signal incident on the array from a direction of 15°, an array element spacing d=0.03, a snapshot number L=400, a signal-to-noise ratio of 10dB, an average bias of the amplitude and phase error increasing from 0.07 to 0.21 with a step of 0.02, and 500 independent Monte Carlo experiments are performed.

[0061] Analysis Figure 2 It can be seen that, with the increase of the average bias of the amplitude and phase error, the success probability of the method of the patent fluctuates around 94%, maintaining good estimation performance; the success probability of the MUSIC method continuously decreases from 84.25% when the average bias of the amplitude and phase error is 0.07 to 39% when the average bias of the amplitude and phase error is 0.21.

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

[0063] The third embodiment: study the relationship between the success probability of the signal processing method of the patent and the signal-to-noise ratio and the relationship between the root mean square error (RMSE) and the signal-to-noise ratio, and obtain the effect diagram as shown in Figure 4 、 Figure 5 The performance analysis comparison diagram of the ISM method is also given. The conditions for applying the algorithm in the application 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, a non-coherent signal incident on the array from a direction of 15°, an array element spacing d=0.03, a snapshot number L=400, a signal-to-noise ratio changing from -5dB to 20dB with a step of 5dB, a phase error following a normal distribution with a variance and a mean of 0, an amplitude error following a uniform distribution of [-0.6, 0.6], and 500 independent Monte Carlo experiments are performed.

[0064] Analysis 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, and the success probability of the method of the patent is always higher than that of the ISM method. When the signal-to-noise ratio is 5dB, the success probability of the method of the patent is increased by 3% than that of the ISM method. When the signal-to-noise ratio is 10dB, the resolution success probability of the method of the patent reaches 100%.

[0065] Analysis Figure 5 It can be seen that when the number of snapshots L=400, the RMSE of the two amplitude and phase error correction algorithms decreases with the increase of the signal-to-noise ratio, and the RMSE of the method of the patent is always lower than that of the ISM method. When the signal-to-noise ratio is 0dB, the RMSE of the method of the patent is decreased by 0.14° than that of the ISM method.

[0066] The fourth embodiment: the root mean square error (RMSE) of the signal processing method of the patent and the number of snapshots are studied respectively, and the effect diagram is shown as Figure 6 The performance analysis comparison diagram of the ISM method is also given. The algorithm application conditions in the application are as follows: a uniform linear array with array element number M=4, auxiliary array element number P=2, and total array element number N=M+P=6 is adopted, a non-coherent signal is incident to the array from the direction of 15°, the array element spacing d=0.03, the signal-to-noise ratio is 10dB, the number of snapshots is changed from L=50 to 400 with a step of 50, the phase error obeys the normal distribution with variance and mean of 0, and the amplitude error obeys the uniform distribution of [-0.6, 0.6], and 500 independent Monte Carlo experiments are carried out.

[0067] Analysis Figure 6 It can be seen that when the signal-to-noise ratio is 10dB, the RMSE of the two amplitude and phase error correction algorithms decreases slightly in the case of small value, and the RMSE of the method of the patent is mostly smaller than that of the ISM method. When the number of snapshots is 100, the RMSE of the method of the patent is decreased by 0.18° than that of the ISM method.

[0068] The fifth embodiment: when the lake test experiment exists amplitude and phase error, the spatial spectrum history diagram obtained by the signal processing method of the patent is shown as Figure 7 The spatial spectrum history diagram of the ISM method is also given. The application conditions of the method of the application are as follows: a uniform linear array with array element number M=4, auxiliary array element number P=2, and total array element number N=M+P=6 is adopted, the auxiliary array element is located at the front end, a far-field narrowband CW pulse signal is incident to the array from the direction of 18°, the center frequency of the signal source is 30kHz, the array element spacing d=0.03, and the signal is divided into 50 groups, and the number of snapshots of each group is L=1000.

[0069] Analysis Figure 7It can be seen that in the presence of amplitude and phase errors, the power of the method of the present patent is almost entirely concentrated near the incident angle of the incident signal, and it can be clearly seen that the change of the angle with time presents 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°, and the actual incident angle of the incident signal cannot be observed from the history diagram. Therefore, the method of the present patent can overcome the influence of amplitude and phase errors in actual engineering applications.

[0070] The specific examples described in the present patent are merely illustrative of the present application. Those skilled in the art to which the present application belongs can make various modifications, supplements or substitutions with similar ways to the described specific examples, but will not deviate from the present application or exceed the scope defined by the appended claims.

Claims

1. A robust hydrophone array bearing estimation method, characterized in that: The underwater target detector uses a uniform linear array 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 array is denoted as x(t)=[x1(t),x2(t),…,x M (t)] T The steps of the bearing estimation method of the hydrophone array are as follows: Step one: Hilbert transform is performed 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 wherein represents the array flow pattern matrix in the presence of amplitude and phase errors, t = 1, 2, …, L represents the sampling sequence number, represents the incident signal vector, represents the noise complex vector of the hydrophone array; Step two: Compute the covariance matrix of the hydrophone array received complex signals according to the formula where L denotes the number of snapshots, [·] H denotes the conjugate transpose of a matrix; Step three: assuming that the first P elements of the uniform linear array are precisely corrected auxiliary elements, and P≥K+1, K represents the number of incident signals, extract the phase information of the subarray 1 covariance matrix where angle(·) represents the extraction of phase information; the upper triangular elements of the phase information matrix are used to calculate the phase error, and the upper triangular matrix composed of the upper triangular elements of the phase information matrix is denoted as where the phase error is denoted as m represents the row label, and n represents the column label;​​ Step four: extract the first P-1 elements of the diagonal immediately above the main diagonal and Since The phase error vector can be obtained as Meanwhile, for the first diagonal elements, each element is added to Let Then, from Combining linear algebra, we have Since The above equation can be simplified as The phase error vector can be obtained from the above equation Step five: the received data of the M hydrophones of the submarine detector are multiplied by the corresponding phase correction factor, the form of the correction factor is The phase correction is realized, and the corrected data is denoted as Step six: Calculate the energy of the discrete signals according to the formula Calculate the data after phase correction Average energy of each hydrophone in the middle Get the amplitude correction vector by the formula Get the amplitude correction vector by the formula Step seven: the phase-corrected data of M hydrophones of the submarine marker detector multiplied by the corresponding amplitude correction factor to correct the amplitude; the corrected data can be represented as X(t); calculate the covariance matrix of the received data matrix after removing the amplitude and phase errors Step eight: θ is uniformly taken in the range [-π / 2, π / 2], and the covariance matrix estimate value in step seven is used The spatial spectrum P(θ) is calculated; Step nine: Calculate the first difference function y(θ) of the spectrum function P(θ) according to the formula Take the absolute value of y(θ) to find the peak value, and take the average of the two peaks with the smallest difference to obtain the estimated value of the incident angle of the incident signal.

2. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The calculation method of the spatial spectrum P(θ) is as follows: A1: performing eigenvalue decomposition on the covariance matrix performing eigenvalue decomposition, sorting in order of eigenvalue from large to small, dividing into signal subspace and noise subspace two parts, wherein the noise subspace U N consisting of eigenvectors corresponding to the K+1 to M-th eigenvalues A2: make θ take value uniformly in the range [-π / 2, π / 2], and calculate the spatial spectrum of the signal according to the formula 3. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The formula for calculating the spatial spectrum P(θ) is:

4. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The formula for calculating the spatial spectrum P(θ) is:

5. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The formula for computing the spatial spectrum P(0) is: where || · ||2denotes the 2-norm of a vector.

6. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The calculation formula of the spatial spectrum P(θ) is: The eigenvalue λ i are sorted from large to small, q i is the eigenvalue λ i corresponding eigenvector.

7. The robust hydrophone array direction of arrival estimation method of claim 1, wherein: The formula for calculating the spatial spectrum P(θ) is:

Images