DOA estimation method for nested hydrophone array with auxiliary array elements arranged at two ends

Abstract

The invention discloses a DOA estimation method for a nested hydrophone array with auxiliary array elements arranged at two ends. According to the technical scheme, the influence of the amplitude-phase error on the DOA estimation performance under the secondary nested array is reduced. According to the scheme, auxiliary array elements are added at the front end, the middle part and the rear end of the nested hydrophone array to form a new array, the new array receiving data is divided into two uniform sub-array data, and the amplitude-phase error correction of the two uniform arrays is realized by using the phase error vector and energy normalization, so that the accuracy of the amplitude-phase error correction 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 the positions of hydrophones, and DOA estimation is carried out in combination with backward smoothing after redundancy removal and reordering. Simulation experiments show that the method provided by the invention greatly 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 placed at both ends. Specifically, it is a method for adding auxiliary array elements at the front end, middle part and rear end of the nested hydrophone array to form a new array, dividing the received data of the new array into two uniform sub-array data, using phase error vector and energy normalization to realize the correction of amplitude and phase errors of the two uniform arrays, re-arranging and combining the corrected data and calculating the covariance matrix of the corrected data; obtaining a virtual array through differential operation of hydrophone positions, removing redundancy and re-arranging, and combining with backward smoothing 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 DOA estimation method for a nested hydrophone array with auxiliary array elements placed at both ends. The patented method adds auxiliary array elements at the front, middle and back ends 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 it, 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 placed at both ends. It is characterized in that: the patent method adds auxiliary array elements at the front, middle and back ends 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 it, and combines it with backward smoothing to perform DOA estimation. While being able to complete the underdetermined DOA estimation when the number of signal sources is greater than the number of array elements, this 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 have 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, (M1+1)d1 and [M2(M1+1)+M1]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 three auxiliary array elements located at -d1, M1d1 and (M1+1)d1, with the number of array elements being M1+3; Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and [M2(M1+1)+M1]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 and last 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 diagonal element 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:

[0013] Step 6: Multiply the received data of the M1+3 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+3 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 =[-1-(M2(M1+1)+M1),-(M2(M1+1)+M1),…,1+M2(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(M1+1)+M1]+3-dimensional row vectors;

[0020] Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, Representation Matrix The element in row i and column j in i -D j =d uk When k=1,2,…,2[M2(M1+1)+M1]+3, we can find w k (d u ) elements, of which the lth can be represented by R l (diff k ),

[0021] Step 13: Divide the difference synthesis array 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 ii=1,2,…,n, represents the equivalent received signal of the ith subarray; the equivalent received signals of the n subarrays are arranged in columns to form a matrix Z=[Z1;Z2;…;Z n ];

[0022] Step 14: The equivalent covariance matrix of the i-th sub-matrix is ​​R i ; Take the average of the equivalent covariance matrices of the n sub-matrices and get

[0023] Step 15: 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;

[0024] Step 16: Make θ take values ​​uniformly 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.

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

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

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

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

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

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

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

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

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

[0034] Step 14: The equivalent covariance matrix R i The calculation formula is

[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 degree of freedom: The patented method of adding auxiliary array elements will not produce holes when the virtual array is differentially expanded, and the degree of freedom is greatly improved.

[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 This is 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, 4d1 and 15d1, 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 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 three auxiliary array elements located at -d1, M1d1 and (M1+1)d1, with the number of array elements being M1+3; Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and [M2(M1+1)+M1]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 and last 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 diagonal element 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:

[0052] Step 6: Multiply the received data of the M1+3 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+3 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 =[-1-(M2(M1+1)+M1),-(M2(M1+1)+M1),…,1+M2(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(M1+1)+M1]+3-dimensional row vectors;

[0059] Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, Representation Matrix The element in row i and column j in i -D j =d uk When k=1,2,…,2[M2(M1+1)+M1]+3, we can find w k (d u ) elements, of which the lth can be represented by R l (diff k ),

[0060] Step 13: Divide the difference synthesis array 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 i=1,2,…,n, represents the equivalent received signal of the ith subarray; the equivalent received signals of the n subarrays are arranged in columns to form a matrix Z=[Z1;Z2;…;Z n ];

[0061] Step 14: The equivalent covariance matrix of the i-th sub-matrix is ​​R i ; Take the average of the equivalent covariance matrices of the n sub-matrices and get

[0062] Step 15: 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;

[0063] Step 16: Make θ take values ​​uniformly 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.

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

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

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

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

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

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

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

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

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

[0073] Step 14: The equivalent covariance matrix R i The calculation formula is

[0074] The DOA estimation results of the patent method and the method without sub-array division are obtained by simulating the patent method based on the above conditions through MATLAB simulation software. Figure 2 As shown, analysis Figure 1 It can be seen that the estimated angles of both methods are concentrated around 30°.

[0075] 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 4 As shown, a performance analysis comparison diagram of the method without sub-array division 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 sub-array is M1=3, the number of array elements of the second-level sub-array is M2=3, the first array element of the first-level sub-array is located at position 0, P=3 accurately corrected auxiliary array elements are added at -d1, 4d1 and 15d1, 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.

[0076] analyze Figure 3 It can be seen that as the signal-to-noise ratio increases, the success probability of the patented method is always higher than that of the method without sub-array division; the patented method reaches 100% when the signal-to-noise ratio is 0dB, and has good estimation performance. When the signal-to-noise ratio is 10dB, the success probability of the patented method is about 20.333% higher than that of the method without sub-array division.

[0077] analyze Figure 4It can be seen that as the signal-to-noise ratio increases, the RMSE of the patented method is always lower than that of the method without sub-array division; the RMSE of the patented method is 0.236° when the signal-to-noise ratio is -5dB, and slowly decreases as the signal-to-noise ratio increases. When the signal-to-noise ratio is 10dB, the RMSE of the patented method is about 0.334° lower than that of the method without sub-array division.

[0078] 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 method without sub-array division 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 sub-array is M1=3, the number of array elements of the second-level sub-array is M2=3, the first array element of the first-level sub-array is located at position 0, P=3 accurately corrected auxiliary array elements are added at -d1, 4d1 and 15d1, 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 0dB, the number of snapshots starts from 100 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.

[0079] analyze Figure 5 It can be seen that as the number of snapshots increases, the success probability of the patented method is always higher than that of the method without sub-array division; when the number of snapshots is 250, the success probability of the patented method reaches 100%.

[0080] analyze Figure 6 It can be seen that with the increase of the signal-to-noise ratio, the RMSE of the patented method is always lower than that of the method without sub-array division; when the number of snapshots is 50, the RMSE of the patented method is about 0.185°, and slowly decreases with the increase of the number of snapshots.

[0081] 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 auxiliary array elements placed at both ends, 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, (M1+1)d1 and [M2(M1+1)+M1]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 three auxiliary array elements located at -d1, M1d1 and (M1+1)d1, with the number of array elements being M1+3; Subarray 2 consists of the second-level subarray and two auxiliary array elements located at -d1 and [M2(M1+1)+M1]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 and last array elements are precisely calibrated auxiliary array elements, let m represents the row number, n represents the column number; Step 5: Extract the first diagonal element next to the main diagonal and because Available At the same time, for the first diagonal element, each element plus And order i=1,2,…,M1+2, then, 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+3 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 i=1,2,…,M1+3, by the formula Get the amplitude correction vector Step 8: Phase-corrected data of M1+3 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 =[-1-(M2(M1+1)+M1),-(M2(M1+1)+M1),…,1+M2(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(M1+1)+M1]+3-dimensional row vectors; Step 12: Let i = 1, 2, ..., N, j = 1, 2, ..., N, Representation Matrix The element in row i and column j in i -D j =d uk When k=1,2,…,2[M2(M1+1)+M1]+3, we can find w k (d u ) elements, of which the lth can be represented by R l (diff k ), Step 13: Divide the difference synthesis array 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 i=1,2,…,n, represents the equivalent received signal of the ith subarray; the equivalent received signals of the n subarrays are arranged in columns to form a matrix Z=[Z1;Z2;…;Z n ]; Step 14: The equivalent covariance matrix of the i-th sub-matrix is ​​R i ; Take the average of the equivalent covariance matrices of the n sub-matrices and get Step 15: 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 16: Make θ take values ​​uniformly 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 auxiliary array elements placed at both ends 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. The method for estimating DOA of a nested hydrophone array with auxiliary array elements placed at both ends as claimed in 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 auxiliary array elements placed at both ends as claimed in 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 auxiliary array elements placed at both ends as claimed in 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 auxiliary array elements placed at both ends 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 method for estimating DOA of a nested hydrophone array with auxiliary array elements placed at both ends as claimed in 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 auxiliary array elements placed at both ends as claimed in claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is i=1,2,…,N,Z j represents the equivalent received signal of the jth sub-array.

9. The method for estimating DOA of a nested hydrophone array with auxiliary array elements placed at both ends as claimed in claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is i=1,2,…,N,Z j represents the equivalent received signal of the jth sub-array.

10. The method for estimating DOA of a nested hydrophone array with auxiliary array elements placed at both ends as claimed in claim 1, characterized in that: Equivalent covariance matrix R i The calculation formula is i=1,2,…,N.

Images