A meshless coherent signal DOA estimation method based on coprime array

Abstract

This invention discloses a meshless DOA estimation method for coherent signals based on coprime matrices. The method first averages the diagonal elements of the covariance matrix of the received data from the array, then reconstructs the Toeplitz matrix based on this average, and finally reconstructs the low-rank matrix using its Toeplitz structure. Next, the trace norm of the positive semi-definite matrix is ​​used to relax the non-convex low-rank matrix reconstruction problem, and the rank recovery matrix is ​​obtained using the convex optimization CVX toolbox. Finally, the DOA is estimated using the ESPRIT algorithm. Compared with traditional methods, this method does not require mesh generation, has stronger decoherence capabilities, and can accurately estimate the azimuth of coherent signals even under conditions of low signal-to-noise ratio, few snapshots, and small signal incident angle intervals. It has significant application value in underwater moorings or mobile observation platforms for ocean observation.

Description

Technical Field

[0001] This invention relates to a target detection and processing method for underwater moorings and mobile observation platforms for ocean observation. It is a gridless DOA estimation method based on a coprime array. Specifically, it is a method for estimating the incident azimuth of a target based on an underwater coprime hydrophone array. This method does not require grid search for the incident direction and directly calculates the estimated azimuth value. Background Technology

[0002] Over the past decade, coprime arrays have attracted significant attention in the field of DOA estimation because they can be extended into virtual uniform linear arrays with higher degrees of freedom, while suffering less coupling effects between elements. Numerous methods for DOA estimation have been proposed in coprime array research; however, most methods only address independent signals. Due to the complexity of real-world underwater acoustic environments, the signals received by the array are typically coherent, leading to rank loss in the array's received data covariance matrix. This severely degrades the estimation performance of these DOA estimation methods, rendering them ineffective for coherent signals.

[0003] Since the 1980s, research on DOA estimation algorithms for coherent signals has attracted considerable attention from scholars, resulting in a wealth of research findings and literature. Among these, most DOA estimation methods for coherent signals using coprime arrays are grid search methods. These methods require dividing the DOA value range into grids, with the number of grids far exceeding the actual number of sources, undoubtedly increasing computational costs. Furthermore, since the true DOA values ​​may not be precisely distributed across the grid points, grid-based methods typically suffer from grid errors, limiting the accuracy of the estimation. If a method could be developed that, by processing the covariance matrix of the array's received data, obtain accurate DOA estimates for coherent signals without requiring grid division of the DOA value range, it would significantly enhance the practical application value of such methods in engineering.

[0004] This patent proposes a meshless DOA estimation method for coherent signals based on coprime arrays. It achieves DOA estimation of coherent signals in the continuous domain of DOA values, saving the computational cost of mesh generation and search. It can obtain the azimuth estimation of coherent signals under conditions of low signal-to-noise ratio, few snapshots, and small incident angle intervals, and has important application value in practical engineering. Summary of the Invention

[0005] The purpose of this invention is to address the limitations of mesh-based methods in terms of estimation accuracy due to mesh errors, and the inability of the Covariance Matrix Reconstruction (CMRA) method to be applied to DOA estimation of coherent signals. This invention proposes a meshless DOA estimation method for coherent signals based on a coprime matrix. The key features are: First, the method averages the diagonal elements of the array-received data covariance matrix, then reconstructs a Toeplitz matrix based on this average, and utilizes its Toeplitz structure to reconstruct a low-rank matrix. Next, the non-convex low-rank matrix reconstruction problem is relaxed by using the trace norm of the positive semi-definite matrix, and the rank recovery matrix is ​​obtained by solving this problem using the convex optimization CVX toolbox. Finally, DOA estimation is performed using the ESPRIT algorithm. Compared to traditional methods, this invention eliminates the need for mesh generation, has stronger decoherence capabilities, and can accurately estimate the azimuth of coherent signals even under conditions of low signal-to-noise ratio, few snapshots, and small signal incident angle intervals, making it of significant application value in practical engineering.

[0006] The signal processing steps of the meshless coherent signal DOA estimation method based on coprime matrix of the present invention include the following steps:

[0007] Step 1: Represent the received data of the L-element coprime array consisting of two subarrays as x(t) = As(t) + n(t), where subarray 1 consists of M elements with an element spacing of Nd, and subarray 2 consists of N elements with an element spacing of Md. M and N are a pair of coprime numbers, d is half the signal wavelength λ, L = M + N - 1, and x(t) = [x1(t), x2(t), ..., x L (t)] T The array receives data matrix, A = [a(θ1), a(θ2), ..., a(θ...]. K ] is an array manifold matrix, Indicates the array steering vector. This represents the set of index numbers of the array elements in a coprime array. θ k Let represent the azimuth angle of the incident signal, k = 1, 2, ..., K, where K represents the number of incident signals, and s(t) = [s1(t), s2(t), ..., s...]. K (t)] T Let n(t) be the signal vector received by the array, where n(t) = [n1(t), n2(t), ..., n L (t)] T Let t = 1, 2, ..., T be the Gaussian white noise data matrix received by the array, where t = 1, 2, ..., T represents the sampling index of the array data, and the superscript "T" in the formula indicates matrix transpose;

[0008] Step 2: According to the formula Calculate the estimated value of the covariance matrix of the received data x(t) of the coprime array, where the superscript "H" indicates the conjugate transpose of the matrix;

[0009] Step 3: According to the formula and r(l)=r * Find the covariance matrix (-l), where l = 0, 1, ..., L-1. The average value of the diagonal elements of the lower triangular part is used to construct the Toplitz matrix.

[0010] Where Γ(u) is the Toplitz matrix constructed from the signal vectors received by the array, u=[u1,u2,…,u L ] T For the first column of Γ(u), N n The Toplitz matrix is ​​constructed for the noise vector received by the array, with the superscript "*" indicating the complex conjugate.

[0011] Step 4: According to the formula Calculate matrix W, and calculate...

[0012] Step 5: Set up an L×Q dimensional matrix G, where Q = MN - N + 1, and the l-th row of matrix G contains only the set The element in the middle is set to 1 at the corresponding position, and all other positions are set to 0;

[0013] Step Six: Minimize the problem using the CVX Convex Optimization Toolbox Solving for stΓ(u)≥0 yields the Q×Q dimensional new covariance matrix Γ(u) after rank recovery, where λ>0 is a user-defined regularization parameter;

[0014] Step 7: According to the formula Perform eigenvalue decomposition on the new covariance matrix Γ(u), Λ s Let U be the diagonal matrix corresponding to the K largest eigenvalues ​​of the new covariance matrix Γ(u). s For its corresponding signal subspace, Λ n U is the diagonal matrix corresponding to the remaining small eigenvalues ​​of the new covariance matrix Γ(u). n This is its corresponding noise subspace;

[0015] Step 8: For the signal subspace U s Divided into two parts Among them, U s1 From the signal subspace U s The first Q-1 rows constitute U s2 From the signal subspace U s The last Q-1 rows constitute the whole;

[0016] Step 9: According to the formula Calculate matrix F;

[0017] Step 10: Perform eigenvalue decomposition on matrix F to obtain K large eigenvalues.

[0018] Step 11: According to the formula Calculate the DOA estimate of the coherent signal.

[0019] Compared with the prior art, the present invention, employing the above technical solution, has the following technical effects:

[0020] (1) The method in this invention obtains the covariance matrix of the received data from the coprime array. The average of the diagonal elements is taken, and the Toeplitz matrix is ​​constructed using the average to achieve decoherence of the coherent signal;

[0021] (2) This invention utilizes the Toeplitz structure of the reconstructed matrix to perform low-rank matrix reconstruction, and relaxes the problem of non-convex low-rank matrix reconstruction by using the trace norm of the positive semi-definite matrix, thereby improving the estimation performance of coherent signals to a certain extent.

[0022] (3) This invention is a meshless algorithm that does not require meshing of the DOA value range, thus reducing computational cost compared to mesh-based algorithms; Attached Figure Description

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

[0024] Figure 2 The diagram shows the DOA estimation results of the signal processing method of this patent for coherent signals.

[0025] Figure 3 The curve showing the relationship between the root mean square error and the signal-to-noise ratio of the signal processing method of this patent is shown.

[0026] Figure 4 The curve showing the relationship between the resolution success probability and the signal-to-noise ratio of the signal processing method of this patent.

[0027] Figure 5 The curve showing the relationship between the root mean square error and the number of snapshots in the signal processing method of this patent.

[0028] Figure 6 The curve showing the relationship between the success rate of resolution and the number of snapshots in the signal processing method of this patent.

[0029] Figure 7 The curve showing the relationship between the root mean square error and the angle interval of the signal processing method of this patent is shown.

[0030] Figure 8The curve showing the relationship between the resolution success probability and the angular interval of the signal processing method of this patent is shown. Detailed Implementation

[0031] The present invention will now be further described in conjunction with the embodiments and accompanying drawings:

[0032] Example 1: Figure 1 The coprime array model in this invention is presented. First, a 10-element coprime array is constructed. Subarray 1 consists of 6 elements with an element spacing of 5d, and subarray 2 consists of 5 elements with an element spacing of 6d, where d is half the signal wavelength, 0.5m. Two coherent signals are incident on the coprime array from directions θ1 = -4° and θ2 = 6°, respectively, with a signal-to-noise ratio of 0dB and a snapshot number T = 200. Azimuth estimation is performed using the above conditions, and the specific implementation process is as follows:

[0033] Step 1: Represent the received data of a 10-element coprime array consisting of two subarrays as x(t) = As(t) + n(t), where subarray 1 consists of 6 elements with an element spacing of 5d, and subarray 2 consists of 5 elements with an element spacing of 6d, where d is half the signal wavelength (0.5m), L = M + N - 1 = 10, and x(t) = [x1(t), x2(t), ..., x 10 (t)] T Let A be the array receiving data matrix, and A = [a(θ1), a(θ2)] be the array manifold matrix. Indicates the array steering vector. This represents the index of the element position in a coprime array. θ k Let represent the azimuth angle of the incident signal, k = 1, 2, ..., K, where K = 2 represents the number of incident signals, and s(t) = [s1(t), s2(t)]. T Let n(t) be the signal vector received by the array, where n(t) = [n1(t), n2(t), ..., n 10 (t)] T Let t = 1, 2, ..., 200 be the Gaussian white noise data matrix received by the array, and let t = 1, 2, ..., 200 represent the sampling index of the array data. The superscript "T" in the formula indicates matrix transpose.

[0034] Step 2: According to the formula Calculate the estimated value of the covariance matrix of the received data x(t) of the coprime array, where the superscript "H" indicates the conjugate transpose of the matrix;

[0035] Step 3: According to the formula and r(l)=r * Find the covariance matrix of (-l), l = 0, 1, ..., 9. The average value of the diagonal elements of the lower triangular part is used to construct the Toplitz matrix.

[0036] Where Γ(u) is the Toplitz matrix constructed from the signal vectors received by the array, u=[u1,u2,…,u 10 ] T For the first column of Γ(u), N n The Toplitz matrix is ​​constructed for the noise vector received by the array, with the superscript "*" indicating the complex conjugate.

[0037] Step 4: According to the formula Calculate matrix W, and calculate...

[0038] Step 5: Set up an L×Q dimensional matrix G, where Q = 26, and the l-th row of matrix G contains only sets. The element in the middle is set to 1 at the corresponding position, and all other positions are set to 0;

[0039] Step Six: Minimize the problem using the CVX Convex Optimization Toolbox Solving for stΓ(u)≥0 yields the rank-restored 26×26-dimensional new covariance matrix Γ(u), where λ>0 is a user-defined regularization parameter, and in this patent, λ=0.25 is set.

[0040] Step 7: According to the formula Perform eigenvalue decomposition on the new covariance matrix Γ(u), Λ s Let U be the diagonal matrix corresponding to the K largest eigenvalues ​​of the new covariance matrix Γ(u). s For its corresponding signal subspace, Λ n U is the diagonal matrix corresponding to the remaining small eigenvalues ​​of the new covariance matrix Γ(u). n This is its corresponding noise subspace;

[0041] Step 8: For the signal subspace U s Divided into two parts Among them, U s1 From the signal subspace U s The first 25 lines constitute U s2 From the signal subspace U s The last 25 lines constitute the whole;

[0042] Step 9: According to the formula Calculate matrix F;

[0043] Step 10: Perform eigenvalue decomposition on matrix F to obtain two large eigenvalues.

[0044] Step 11: According to the formula The DOA estimate of the coherent signal was calculated, and the DOA estimation result of this patented method was obtained by conducting 50 independent experiments. Figure 2As shown, from Figure 2 The simulation results show that the angle estimated by the method of this patent is basically consistent with the actual angle.

[0045] The second embodiment: The relationship between the root mean square error (RMSE) and the signal-to-noise ratio (SNR) and the probability of successful resolution and the SNR of the signal processing method of this patent were studied respectively, and the effect diagram is shown below. Figure 3 , Figure 4 As shown. The application conditions of the algorithm in this invention are as follows:

[0046] We used a 10-element coprime array, where subarray 1 consists of 6 elements with an element spacing of 5d, and subarray 2 consists of 5 elements with an element spacing of 6d, where d is half the signal wavelength, 0.5m. Two coherent signals were incident on the coprime array from θ1 = -4° and θ2 = 6°, respectively. The number of snapshots T = 200. The signal-to-noise ratio was changed from -15dB to 5dB in 5dB increments, and 200 independent Monte Carlo experiments were conducted. The simulation was performed using MATLAB software.

[0047] analyze Figure 3 It can be seen that when the number of snapshots N=200, the RMSE of the four algorithms decreases as the signal-to-noise ratio increases, and the RMSE of the algorithm of this patent is always lower than that of the other three algorithms. When the signal-to-noise ratio is 5dB, the RMSE of the algorithm of this patent is about 0.127°, which is about 0.060° higher than the forward smoothing algorithm, about 0.188° higher than the forward smoothing algorithm, and about 0.320° higher than the coprime array interpolation algorithm.

[0048] analyze Figure 4 It can be seen that when the number of snapshots N=200, the success rate of the four algorithms increases with the increase of the signal-to-noise ratio, and the success rate of the algorithm of this patent is always higher than that of the other three algorithms. When the signal-to-noise ratio is -15dB, the success rate of the algorithm of this patent is about 12.5%, while the success rate of the other three algorithms is almost zero. When the signal-to-noise ratio is -5dB, the success rate of the algorithm of this patent reaches 100%, while the success rates of the other algorithms are still low.

[0049] The third embodiment: The relationship between the root mean square error (RMSE) and the number of snapshots and the success probability of resolution and the number of snapshots in the signal processing method of this patent were studied respectively, and the effect diagram is shown below. Figure 5 , Figure 6 As shown. The application conditions of the algorithm in this invention are as follows:

[0050] We used a 10-element coprime array, where subarray 1 consists of 6 elements with an element spacing of 5d, and subarray 2 consists of 5 elements with an element spacing of 6d, where d is half the signal wavelength, 0.5m. Two coherent signals were incident on the coprime array from θ1 = -4° and θ2 = 6°, respectively, with a signal-to-noise ratio snr = 0dB. We changed the number of snapshots from 50 to 400 in steps of 50, and conducted 200 independent Monte Carlo experiments, which were simulated using MATLAB software.

[0051] analyze Figure 5 It can be seen that when the signal-to-noise ratio is 0dB, the RMSE of the four algorithms decreases with the increase of the number of snapshots, and the RMSE of the algorithm in this patent decreases the fastest, tending to stabilize when the number of snapshots is 150; when the number of snapshots is 400, the RMSE of the algorithm in this patent is about 0.141°, which is about 0.106° higher than the forward smoothing algorithm; about 0.310° higher than the forward smoothing algorithm; and about 0.369° higher than the coprime array interpolation algorithm.

[0052] analyze Figure 6 It can be seen that when the signal-to-noise ratio is 0dB, the success rate of the four algorithms increases with the increase of the number of snapshots. Moreover, the success rate of the algorithm in this patent is always higher than that of the other three algorithms. When the number of snapshots is 50, the success rate of the algorithm in this patent is about 93%, which is about 78.5% higher than the forward smoothing algorithm, about 86% higher than the forward smoothing algorithm, and about 69% higher than the coprime array interpolation algorithm.

[0053] Fourth embodiment: The relationship between the root mean square error (RMSE) and the angle interval, and the relationship between the success probability of resolution and the angle interval of the signal processing method of this patent are studied respectively, and the effect diagram is shown below. Figure 7 , Figure 8 As shown. The application conditions of the algorithm in this invention are as follows:

[0054] We used a 10-element coprime array, where subarray 1 consists of 6 elements with an element spacing of 5d, and subarray 2 consists of 5 elements with an element spacing of 6d, where d is half the signal wavelength, 0.5m. Two coherent signals were incident on the coprime array from θ1 = -4° and θ2 = 6°, respectively. The signal-to-noise ratio snr = 0dB, the number of snapshots was 200, and the angle interval was changed from 1° to 10° in increments of 1°. 200 independent Monte Carlo experiments were conducted and simulated using MATLAB software.

[0055] analyze Figure 7It can be seen that when the signal-to-noise ratio is 0dB and the number of snapshots N=200, the RMSE of the four algorithms decreases as the incident angle interval increases, and the RMSE of the algorithm in this patent decreases the fastest. When the incident angle interval is 10°, the RMSE of the algorithm in this patent is about 0.141°, which is about 0.250° higher than the forward smoothing algorithm, about 0.502° higher than the forward smoothing algorithm, and about 0.394° higher than the coprime array interpolation algorithm.

[0056] analyze Figure 8 It can be seen that when the signal-to-noise ratio is 0dB and the number of snapshots N=200, the success rate of the four algorithms increases with the increase of the incident angle interval. The success rate of the algorithm of this patent is always higher than that of the other three algorithms. When the incident angle interval is 6°, the success rate of the algorithm of this patent reaches 100% and tends to stabilize, while the success rates of the other three algorithms are still low.

[0057] The specific examples described herein are merely illustrative of the invention. Those skilled in the art to which this invention pertains may make various modifications, additions, or similar substitutions to the described specific examples without departing from the invention or exceeding the scope defined by the appended claims.

Claims

1. A meshless coherent signal DOA estimation method based on coprime arrays, characterized in that: The DOA estimation method comprises the following steps: Step 1: Represent the received data of the L-element coprime array consisting of two subarrays as x(t) = As(t) + n(t), where subarray 1 consists of M elements with an element spacing of Nd, and subarray 2 consists of N elements with an element spacing of Md. M and N are a pair of coprime numbers, d is half the signal wavelength λ, L = M + N - 1, and x(t) = [x1(t), x2(t), ..., x L (t)] T The array receives data matrix, A = [a(θ1), a(θ2), ..., a(θ...]. K ] is an array manifold matrix, Indicates the array steering vector. This represents the set of index numbers of the elements in a coprime array. θ k Let represent the azimuth angle of the incident signal, k = 1, 2, ..., K, where K represents the number of incident signals, and s(t) = [s1(t), s2(t), ..., s...]. K (t)] T Let n(t) be the incident signal vector received by the array, where n(t) = [n1(t), n2(t), ..., n L (t)] T Let t = 1, 2, ..., T represent the sampling index of the array data, where the superscript "T" indicates matrix transpose. Step two, compute the estimate of the covariance matrix of the received data x(t) according to the formula RHH= x(t)Hx(t) (1) where the superscript "H" denotes the conjugate transpose of a matrix. Step three, the average of the diagonal elements of the lower triangular part of the covariance matrix R is calculated according to the formula and r(l) = r * (-l), l = 0, 1, …, L-1 The average of the diagonal elements of the lower triangular part of the covariance matrix R is calculated according to the formula where Γ(u) is a Toeplitz matrix constructed by the array received signal vector, u = [u1, u2, …, u L ] T is the first column of Γ(u), N n is a Toeplitz matrix constructed by the array received noise vector, and the superscript "*" represents taking the complex conjugate. Step four, compute the matrix W, and compute compute the matrix W, and compute Step five, set L x Q dimensional matrix G, where Q = MN - N + 1, the first l row of matrix G only has the set corresponding position of the element is 1, and all other positions are 0; Step six, solve the minimization problem by CVX convex optimization toolbox s.t. Γ(u) ≥ 0, to obtain the new Q x Q dimensional covariance matrix Γ(u) after rank restoration, where λ > 0 is a user-defined regularization parameter; Step seven, according to the formula Eigenvalue decomposition is performed on the new covariance matrix Γ(u), Λ s is the diagonal matrix corresponding to the K largest eigenvalues of the new covariance matrix Γ(u), U s is the signal subspace corresponding thereto, Λ n is the diagonal matrix corresponding to the remaining small eigenvalues of the new covariance matrix Γ(u), U n is the noise subspace corresponding thereto; Step eight, signal subspace U s Divided into two parts where U s1 consists of the first Q - 1 rows of U s consists of the first Q - 1 rows of U s2 consists of the last Q - 1 rows of U s consists of the last Q - 1 rows of U Step nine, compute the matrix F according to the formula F = (A - B) Step ten, eigenvalue decomposition of matrix F to obtain K large eigenvalues Step eleven, calculate the DOA estimate of the coherent signal according to the formula DOA estimate of the coherent signal.

Images