DOA estimation method for nested hydrophone array with preposed auxiliary array elements

Abstract

The invention discloses a DOA estimation method for a nested hydrophone array with preposed auxiliary array elements. According to the scheme, the influence of amplitude and phase errors on the DOA estimation precision of the second-level nested array can be reduced, auxiliary array elements are added to the front end of the nested hydrophone array to form a new array, new array receiving data are divided into two uniform sub-array data, the amplitude and phase errors of the two uniform arrays are corrected through phase error vector and energy normalization, and the DOA estimation precision of the second-level nested array is improved. Reordering and combining the corrected data and calculating a covariance matrix of the corrected data; a virtual array is obtained through differential operation of hydrophone positions, after redundancy is removed and reordering is conducted, optimized filling of holes is conducted through nuclear norm minimization, and DOA estimation of incident signals is conducted in combination with backward smoothing. Simulation experiments show that the method provided by the invention reduces the influence of amplitude-phase errors on the DOA estimation performance of the nested array, improves the estimation precision, and has a relatively high 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 method for estimating DOA of a nested hydrophone array with auxiliary array elements in front. Specifically, an auxiliary array element is added to the front end of the nested hydrophone array to form a new array, the new array receiving data is divided into two uniform sub-array data, the phase error vector and energy normalization are used to correct the amplitude and phase errors of the two uniform arrays, the corrected data are reordered and combined, and the covariance matrix of the corrected data is calculated; a virtual array is obtained by differential operation of the hydrophone position, redundancy is removed and reordered, and the holes are filled by minimizing the nuclear norm, and a smoothing method is combined 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. 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 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. To address this problem, the phase information in the covariance matrix and energy normalization are used to correct the amplitude and phase errors.

[0004] In addition, the uniform linear array has limited degrees of freedom (DOF). A uniform linear array with N array elements can estimate at most N-1 signal sources. Moreover, in actual engineering applications, increasing the number of array elements will lead to a sharp increase in hardware equipment and costs. To address this problem, the idea of ​​nested arrays was proposed. However, when performing phase error correction on arrays with uneven array element positions, such as nested arrays, the calculated phase error is prone to deviation due to the uneven array element positions, which in turn causes a decrease in estimation performance.

[0005] This patent proposes a nested hydrophone array DOA estimation method with auxiliary array elements in front. This patented method adds auxiliary array elements at the front end of the nested hydrophone array to form a new array, divides the received data of the new array into two uniform sub-array data, uses phase error vector and energy normalization to realize the correction of amplitude and phase errors of the two uniform arrays, reorders and combines the corrected data, and calculates the covariance matrix of the corrected data; obtains a virtual array through differential operation of the hydrophone position, removes redundancy and reorders, fills the holes by minimizing the nuclear norm, and combines the smoothing method to perform DOA estimation. While being able to complete the underdetermined DOA estimation when the number of sources is greater than the number of array elements, it can well eliminate the influence of amplitude and phase errors on the estimation results, improve the estimation accuracy, and at the same time fill the holes, increase the degree of freedom, and have higher practical engineering application value. Summary of the invention

[0006] This patent aims to solve the problem that the estimation performance of the traditional nested array DOA estimation method is significantly reduced due to the influence of amplitude and phase errors in actual engineering, and proposes a nested hydrophone array DOA estimation method with auxiliary array elements in front. The characteristics are as follows: the patent method adds auxiliary array elements at the front end of the nested hydrophone array to form a new array, divides the received data of the new array into two uniform sub-array data, uses the phase error vector and energy normalization to realize the correction of the amplitude and phase errors of the two uniform arrays, reorders and combines the corrected data, and calculates the covariance matrix of the corrected data; obtains a virtual array through the differential operation of the hydrophone position, removes redundancy and reorders, fills the holes by minimizing the nuclear norm, and performs DOA estimation in combination with backward smoothing. While being able to complete the underdetermined DOA estimation when the number of sources is greater than the number of array elements, the patent method can well eliminate the influence of amplitude and phase errors on the estimation results, improve the estimation accuracy, and at the same time fill the holes, increase the degree of freedom, and has higher practical engineering application value.

[0007] This patented method uses M hydrophones to form a signal receiving system according to a two-level nested array. The element spacing of the first-level subarray is It is composed of M1 hydrophones; the array element spacing of the second-level subarray is d2=(M1+1)d1, and it is composed of M2 hydrophones, where M1+M2=M, and λ is the wavelength of the signal; the first array element position of the nested array is regarded as 0, and P=3 corrected auxiliary array elements are placed at -d1, -2d1 and -(M1+2)d1; the number of array elements in the new array is N=M+P; the physical position of the array elements is recorded as D=[-(M1+2),-2,-1,0,1,…,M1-1,M1,2M1+1,…,(M2-1)(M1+1)+M1] T *d1; the received data of N hydrophone arrays is recorded as x(t)=[x1(t),x2(t),…,x N (t)] T , the steps of the DOA estimation method of the hydrophone array are as follows:

[0008] Step 1: Perform Hilbert transform on the received data of N 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: Divide the new array into two uniform linear arrays. Subarray 1 consists of the first-level subarray and two auxiliary array elements located at -d1 and -2d1, with the number of array elements being M1+2. Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and -(M1+2)d1, with the number of array elements being M2+2.

[0010] Step 3: Since the received data vector x i (t)i=1,2,…,N and the array element position have a one-to-one correspondence, so the received data of the new array of N array elements is Receive data divided into two uniform linear arrays and According to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 1 Where L is the number of snapshots, [·] H represents the conjugate transpose of the matrix; according to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 2

[0011] Step 4: 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 the first two array elements are precisely calibrated auxiliary array elements, let m represents the row number, n represents the column number;

[0012] Step 5: Extract the first element of the first diagonal line next to the main diagonal 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:

[0013] Step 6: Multiply the received data of the M1+2 hydrophones of subarray 1 by the corresponding phase correction factor. The correction factor is in the form of To achieve phase correction, the corrected data is recorded as

[0014] Step 7: 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

[0015] Step 8: Phase-corrected data of M1+2 hydrophones in subarray 1 Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data of subarray 1 can be expressed as X1(t);

[0016] Step 9: Data covariance sub-matrix of sub-matrix 2 Repeat the calculation process of step 4, step 5, step 6, step 7 and step 8 to obtain the corrected data of subarray 2, which can be expressed as X2(t);

[0017] Step 10: According to the position relationship of each array element in the two uniform linear arrays, the vectors in X1(t) and X2(t) are calculated according to the formula

[0018] Sort and reorganize in the array to obtain the received data matrix X(t) after removing the amplitude and phase errors of the new array; according to the formula Calculate the covariance matrix of the data matrix X(t);

[0019] Step 11: Vectorize the covariance matrix Get virtual array element receiving data The position of the virtual array element is the difference set S of the original array element position D ={(D i -D j )×d1,1≤i,j≤N} are arranged in ascending order according to the position of the virtual array element, denoted by d u =[-5-(M2-1)(M1+1)-M1,-2-(M2-1)(M1+1)-M1,…,5+(M2-1)(M1+1)+M1] T d1; u The number of repetitions of each element in the array is recorded as w(d u ), d u and w(d u ) are all 2[(M2-1)(M1+1)+M1]+7-dimensional row vectors;

[0020] Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, when D i -D j =d uk When k = 1, 2, ..., 2[(M2-1)(M1+1)+M1]+7, it can be Find out w k (d u ) elements, of which the lth can be represented as r l (diff k ), the data corresponding to the kth virtual array element is at this time There are holes in the , the positions of the holes are ±(4+(M2-1)(M1+1)+M1)d1 and ±(3+(M2-1)(M1+1)+M1)d1, fill the holes and get

[0021] Step 13: Divide the difference matrix into n=(M 2 / 4+M / 2+5) sub-arrays, the number of elements in each sub-array is (M 2 / 4+M / 2+5), from the first (M 2 / 4+M / 2-i+6) to [(M 2 -2) / 2+M-i+11] array elements are the i-th sub-matrix, and the row vector Middle (M 2 / 4+M / 2-i+6) to [(M 2-2) / 2+M-i+11] elements constitute a column vector, denoted as Z i , represents the equivalent received signal of the i-th sub-array, where i = 1, 2, ..., n; the equivalent received signals of the n sub-arrays are arranged in columns to form a matrix Z = [Z1; Z2; ...; Z n ];

[0022] Step 14: According to the formula Optimizing the matrix Z can yield the matrix F, where ||·|| * represents the nuclear norm of the matrix, represents the Hadamard product, the matrix Its expression is F is the received data of each sub-array after the hole optimization filling;

[0023] Step 15: Use F i Construct the equivalent covariance matrix of the i-th submatrix as R i , where F i represents the i-th column of matrix F; averaging the equivalent covariance matrices of n sub-matrices, we get

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

[0025] Step 17: Make θ uniformly take values ​​in the range of [-π / 2,π / 2], calculate the spatial spectrum P(θ) of the signal; find the peak value and obtain the estimated value of the incident angle of the incident signal.

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

[0027] 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;

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

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

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

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

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

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

[0034] Step 15: The equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-matrix. In step 15, the equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-array.

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

[0036] 1. When there are amplitude and phase errors in the received data of the array elements, the estimation performance is better: an array with uneven element positions is divided into two uniform sub-arrays. The effect after amplitude and phase error correction is better, the estimation accuracy is significantly improved, and it has higher practical engineering application value.

[0037] 2. Improved degrees of freedom: The nuclear norm is used to minimize the filling of holes, making full use of the degree of freedom of the nested matrix.

[0038] 3. It can realize underdetermined DOA estimation when the number of signal sources is greater than the number of array elements. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 This is an array model diagram of the signal processing method of this patent;

[0040] Figure 2 This is the spatial spectrum of the patented signal processing method;

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

[0042] 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;

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

[0044] 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

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

[0046] Embodiment 1: The number of array elements of the first-level subarray is M1=3, the number of array elements of the second-level subarray is M2=3, the first array element of the first-level subarray is located at position 0, P=3 accurately calibrated auxiliary array elements are added at -d1, -2d1 and -5d1, the number of array elements of the new array is M=M1+M2+P=9, an incoherent signal is incident on the array from a direction of 30°, the array element spacing is d=0.03, the number of snapshots is L=300, 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:

[0047] Step 1: Perform Hilbert transform on the received data of N 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;

[0048] Step 2: Divide the new array into two uniform linear arrays. Subarray 1 consists of the first-level subarray and two auxiliary array elements located at -d1 and -2d1, with the number of array elements being M1+2. Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and -(M1+2)d1, with the number of array elements being M2+2.

[0049] Step 3: Since the received data vector x i (t)i=1,2,…,N and the array element position have a one-to-one correspondence, so the received data of the new array of N array elements is Receive data divided into two uniform linear arrays and According to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 1 Where L is the number of snapshots, [·] H represents the conjugate transpose of the matrix; according to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 2

[0050] Step 4: 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 the first two array elements are precisely calibrated auxiliary array elements, let m represents the row number, n represents the column number;

[0051] Step 5: Extract the first element of the first diagonal line next to the main diagonal 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:

[0052] Step 6: Multiply the received data of the M1+2 hydrophones of subarray 1 by the corresponding phase correction factor. The correction factor is in the form of To achieve phase correction, the corrected data is recorded as

[0053] Step 7: 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

[0054] Step 8: Phase-corrected data of M1+2 hydrophones in subarray 1 Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data of subarray 1 can be expressed as X1(t);

[0055] Step 9: Data covariance sub-matrix of sub-matrix 2 Repeat the calculation process of step 4, step 5, step 6, step 7 and step 8 to obtain the corrected data of subarray 2, which can be expressed as X2(t);

[0056] Step 10: According to the position relationship of each array element in the two uniform linear arrays, the vectors in X1(t) and X2(t) are calculated according to the formula

[0057] Sort and reorganize in the array to obtain the received data matrix X(t) after removing the amplitude and phase errors of the new array; according to the formula Calculate the covariance matrix of the data matrix X(t);

[0058] Step 11: Vectorize the covariance matrix Get virtual array element receiving data The position of the virtual array element is the difference set S of the original array element position D ={(D i -D j )×d1,1≤i,j≤N} are arranged in ascending order according to the position of the virtual array element, denoted by d u =[-5-(M2-1)(M1+1)-M1,-2-(M2-1)(M1+1)-M1,…,5+(M2-1)(M1+1)+M1] T d1; u The number of repetitions of each element in the array is recorded as w(d u ), d u and w(d u ) are all 2[(M2-1)(M1+1)+M1]+7-dimensional row vectors;

[0059] Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, when D i -D j =d uk When k = 1, 2, ..., 2[(M2-1)(M1+1)+M1]+7, it can be Find out w k (d u ) elements, of which the lth can be represented as r l (diff k ), the data corresponding to the kth virtual array element is at this time There are holes in the , the positions of the holes are ±(4+(M2-1)(M1+1)+M1)d1 and ±(3+(M2-1)(M1+1)+M1)d1, fill the holes and get

[0060] Step 13: Divide the difference matrix into n=(M 2 / 4+M / 2+5) sub-arrays, the number of elements in each sub-array is (M 2 / 4+M / 2+5), from the first (M 2 / 4+M / 2-i+6) to [(M 2 -2) / 2+M-i+11] array element is the i-th sub-matrix, and the (M-th 2 / 4+M / 2-i+6) to [(M 2 -2) / 2+M-i+11] elements constitute a column vector, denoted as Z i , represents the equivalent received signal of the i-th sub-array, where i = 1, 2, ..., n; the equivalent received signals of the n sub-arrays are arranged in columns to form a matrix Z = [Z1; Z2; ...; Z n ];

[0061] Step 14: According to the formula Optimizing the matrix Z can yield the matrix F, where ||·|| * represents the nuclear norm of the matrix, represents the Hadamard product, the matrix Its expression is F is the received data of each sub-array after the hole optimization filling;

[0062] Step 15: Use F i Construct the equivalent covariance matrix of the i-th submatrix as R i , where F i represents the i-th column of matrix F; averaging the equivalent covariance matrices of n sub-matrices, we get

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

[0064] Step 17: Make θ uniformly take values ​​in the range of [-π / 2,π / 2], calculate the spatial spectrum P(θ) of the signal; find the peak value and obtain the estimated value of the incident angle of the incident signal.

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

[0066] 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;

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

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

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

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

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

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

[0073] Step 15: The equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-matrix. In step 15, the equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-array.

[0074] Step 15: The equivalent covariance matrix R i The calculation formula is R i =(F i F i H )i=1,2,…,N.

[0075] The patent method based on the above conditions is simulated by MATLAB simulation software, and the DOA estimation results of the patent method, the hole removal method and the non-sub-array division method are as follows: Figure 2 As shown, analysis Figure 1 It can be seen that the estimated angles of the two methods are concentrated around 30°, among which the estimated angle of the method in the present patent is closer to the incident angle of the signal, and can accurately estimate the angle of the incident signal when there are amplitude and phase errors in the received data.

[0076] Second 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 3 , Figure 4As shown, a performance analysis comparison diagram of the hole removal method and the non-divided subarray method is also given. The application conditions of the method of the present invention are as follows: the number of array elements of the first-level subarray is M1=3, the number of array elements of the second-level subarray is M2=3, the first array element of the first-level subarray is located at position 0, P=3 accurately corrected auxiliary array elements are added at -d1, -2d1 and -5d1, the number of array elements of the new array is M=M1+M2+P=9, an incoherent signal is incident on the array from the direction of 30, the array element spacing d=0.03, the number of snapshots L=300, the signal-to-noise ratio starts from -5dB, and increases 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.

[0077] analyze Figure 3 It can be seen that with the increase of the signal-to-noise ratio, the success probability of the patented method is slightly higher than the success probability of the hole removal method, and much higher than the method without sub-array division; the patented method reaches 100% when the signal-to-noise ratio is 5dB, and has good estimation performance.

[0078] analyze Figure 4 It can be seen that with the increase of the signal-to-noise ratio, the RMSE of the patented method is slightly lower than the hole removal method, and much lower than the method without sub-array division; the RMSE of the patented method drops to about 0.355° when the signal-to-noise ratio reaches 0dB, and slowly decreases with the increase of the signal-to-noise ratio.

[0079] The third embodiment: The relationship between the success probability and the number of snapshots and the relationship between the root mean square error (RMSE) and the number of snapshots of the signal processing method of the patent are studied respectively, and the effect diagram is as follows: Figure 5 , Figure 6 As shown, a performance analysis comparison diagram of the hole removal method and the non-dividing subarray method is also given. The application conditions of the method of the present invention are as follows: the number of array elements of the first-level subarray is M1=3, the number of array elements of the second-level subarray is M2=3, the first array element of the first-level subarray is located at position 0, P=3 accurately corrected auxiliary array elements are added at -d1, -2d1 and -5d1, the number of array elements of the new array is M=M1+M2+P=9, an incoherent signal is incident on the array from the direction of 30, the array element spacing is d=0.03, the signal-to-noise ratio is 10dB, the number of snapshots starts from 150 and increases to 500 with a step size of 50, 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.

[0080] analyze Figure 5It can be seen that as the number of snapshots increases, the success probability of the patented method is slightly higher than that of the hole removal method, and much higher than that of the method without sub-array division; when the number of snapshots is 200, compared with the hole removal method, the success probability of the patented method is increased by 2%, and compared with the method without sub-array division, the success probability of the patented method is increased by 59%.

[0081] analyze Figure 6 It can be seen that with the increase of the signal-to-noise ratio, the RMSE of the patented method is slightly lower than that of the hole removal method and much lower than that of the non-sub-array division method; when the number of snapshots is 200, compared with the hole removal method, the RMSE of the patented method is reduced by 0.045°, and compared with the non-sub-array division method, the RMSE of the patented method is reduced by 1.359°.

[0082] 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 method for estimating DOA of a nested hydrophone array with an auxiliary array element in front, characterized by: The buoy detector uses M hydrophones to form a signal receiving system according to the two-level nested array. The element spacing of the first-level subarray is It is composed of M1 hydrophones; the array element spacing of the second-level subarray is d2=(M1+1)d1, and it is composed of M2 hydrophones, where M1+M2=M, and λ is the wavelength of the signal; the first array element position of the nested array is regarded as 0, and P=3 corrected auxiliary array elements are placed at -d1, -2d1 and -(M1+2)d1; the number of array elements in the new array is N=M+P; the physical position of the array elements is recorded as D=[-(M1+2),-2,-1,0,1,…,M1-1,M1,2M1+1,…,(M2-1)(M1+1)+M1] T *d1; the received data of N hydrophone arrays is recorded as x(t)=[x1(t),x2(t),…,x N (t)] T , the steps of the DOA estimation method of the hydrophone array are as follows: Step 1: Perform Hilbert transform on the received data of N 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: Divide the new array into two uniform linear arrays. Subarray 1 consists of the first-level subarray and two auxiliary array elements located at -d1 and -2d1, with the number of array elements being M1+2. Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and -(M1+2)d1, with the number of array elements being M2+2. Step 3: Since the received data vector x i (t)i=1,2,…,N and the array element position have a one-to-one correspondence, so the received data of the new array of N array elements is Receive data divided into two uniform linear arrays and According to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 1 Where L is the number of snapshots, [·] H represents the conjugate transpose of the matrix; according to the formula Calculate the covariance matrix of the complex signal of the received data of subarray 2 Step 4: 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 the first two array elements are precisely calibrated auxiliary array elements, let m represents the row number, n represents the column number; Step 5: Extract the first element of the first diagonal line next to the main diagonal 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 6: Multiply the received data of the M1+2 hydrophones of subarray 1 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 7: 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 8: Phase-corrected data of M1+2 hydrophones in subarray 1 Multiply by the corresponding amplitude correction factor to achieve amplitude correction; the corrected data of subarray 1 can be expressed as X1(t); Step 9: Data covariance sub-matrix of sub-matrix 2 Repeat the calculation process of step 4, step 5, step 6, step 7 and step 8 to obtain the corrected data of subarray 2, which can be expressed as X2(t); Step 10: According to the position relationship of each array element in the two uniform linear arrays, the vectors in X1(t) and X2(t) are calculated according to the formula Sort and reorganize in the array to obtain the received data matrix X(t) after removing the amplitude and phase errors of the new array; according to the formula Calculate the covariance matrix of the data matrix X(t); Step 11: Vectorize the covariance matrix Get virtual array element receiving data The position of the virtual array element is the difference set S of the original array element position D ={(D i -D j )×d1,1≤i,j≤N} are arranged in ascending order according to the position of the virtual array element, denoted by d u =[-5-(M2-1)(M1+1)-M1,-2-(M2-1)(M1+1)-M1,…,5+(M2-1)(M1+1)+M1] T d1; u The number of repetitions of each element in the array is recorded as w(d u ), d u and w(d u ) are all 2[(M2-1)(M1+1)+M1]+7-dimensional row vectors; Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, when D i -D j =d uk When k = 1, 2, ..., 2[(M2-1)(M1+1)+M1]+7, it can be Find out w k (d u ) elements, of which the lth can be represented as r l (diff k ), the data corresponding to the kth virtual array element is at this time There are holes in the , the positions of the holes are ±(4+(M2-1)(M1+1)+M1)d1 and ±(3+(M2-1)(M1+1)+M1)d1, fill the holes and get Step 13: Divide the difference matrix into n=(M 2 / 4+M / 2+5) sub-arrays, the number of elements in each sub-array is (M 2 / 4+M / 2+5), from the first (M 2 / 4+M / 2-i+6) to [(M 2 -2) / 2+M-i+11] array elements are the i-th sub-matrix, and the row vector Middle (M 2 / 4+M / 2-i+6) to [(M 2 -2) / 2+M-i+11] elements constitute a column vector, denoted as Z i , represents the equivalent received signal of the i-th sub-array, where i = 1, 2, ..., n; the equivalent received signals of the n sub-arrays are arranged in columns to form a matrix Z = [Z1; Z2; ...; Z n ]; Step 14: According to the formula Optimizing the matrix Z can yield the matrix F, where ||·|| * represents the nuclear norm of the matrix, represents the Hadamard product, the matrix Its expression is F is the received data of each sub-array after the hole optimization filling; Step 15: Use F i Construct the equivalent covariance matrix of the i-th submatrix as R i , where F i represents the i-th column of matrix F; averaging the equivalent covariance matrices of n sub-matrices, we get Step 16: Covariance matrix R SS Perform eigenvalue decomposition and sort the eigenvalues ​​from large to small, dividing them into signal subspace and noise subspace. N It is composed of eigenvectors corresponding to the K+1th to nth eigenvalues; Step 17: Make θ uniformly take values ​​in the range of [-π / 2,π / 2], calculate the spatial spectrum P(θ) of the signal; find the peak value and obtain the estimated value of the incident angle of the incident signal.

2. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of 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 a uniform value in the range of [-π / 2,π / 2], and calculate the spatial spectrum of the signal according to the formula:

3. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

4. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

5. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is: where ||·||2 represents the 2-norm of the vector.

6. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of 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 method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: The calculation formula of spatial spectrum P(θ) is:

8. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-array.

9. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is F j represents the equivalent received signal of the jth sub-array.

10. The method for estimating DOA of a nested hydrophone array with an auxiliary array element in front of claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is R i =(F i F i H )i=1,2,…,N.

Images