An improved variational bayesian sparse learning outlier azimuth estimation method

Abstract

The application provides an improved sparse learning out-of-grid direction-of-arrival estimation method based on variational Bayesian. The method is characterized in that: the original data received by a hydrophone array is preprocessed, a real value transformation is used to convert a vectorized covariance matrix signal in a complex number field to a real number field, and the idea of variational sparse Bayesian learning and grid evolution is combined to make the grid evolve from an initial uniform grid to a non-uniform grid adaptively in an iteration process. The evolution process includes grid updating and grid fission. The evolved grid points are gradually close to the real source position through the alternately iterative grid updating process and grid fission process. Compared with the traditional compressed sensing method, the method has higher DOA estimation accuracy, reduces the operation complexity, optimizes the operation efficiency, improves the resolution capacity of the source, and has higher application value in actual engineering, especially in the case of few snapshots and low signal-to-noise ratio.

Description

Technical Field

[0001] This invention belongs to the field of array signal processing technology, specifically relating to an improved variational Bayesian sparse learning off-grid orientation estimation method. Background Technology

[0002] As an important branch of source parameter estimation, Direction-of-Arrival (DOA) estimation has always received attention from scholars and has been widely applied in radar, sonar, and mobile communications. Over the past few decades, subspace-based DOA estimation algorithms have become well-known, ushering in a high-resolution era for DOA estimation. The most representative subspace-based algorithms are those based on multiple signal classification (MUSIC) and those based on the estimation signal parameter via rotational invariance technique (ESPRIT). However, these subspace-based algorithms rely on a large number of snapshots and a high signal-to-noise ratio (SNR) to obtain accurate covariance. This leads to performance degradation or failure at low SNR or with few snapshots, and these algorithms are also ineffective for coherent sources.

[0003] In recent years, with the rapid development of compressed sensing technology, sparse signal recovery (SSR) has become a very popular DOA estimation technique. Sparse reconstruction algorithms can effectively overcome the problems of traditional subspace-based algorithms, maintaining good estimation performance even under low signal-to-noise ratio and limited snapshot conditions, while also improving robustness to noise. The L1-SVD method first applied SSR to DOA estimation, relaxing a non-convex L0 norm sparsity recovery problem into a convex L1 norm minimization problem, using L1 norm and singular value decomposition to enhance sparsity and reduce complexity. Although this method offers high estimation accuracy, it is computationally complex and computationally intensive. To address these issues, Bayesian inference was introduced. With the continuous development of compressed sensing theory, Sparse Bayesian Learning (SBL) is considered to have the same global convergence as L1-SVD-like convex optimization methods, and possesses significantly better computational efficiency than convex optimization algorithms. Traditional DOA estimation methods based on sparse Bayesian learning (SBL) improve resolution and enhance robustness against interference. However, the estimation performance of SBL-based algorithms is affected by the degree of spatial discretization. High discretization leads to high computational complexity, while low discretization results in reduced estimation accuracy due to off-grid errors, i.e., grid mismatch. To address this issue, Yang et al. proposed an off-grid sparse Bayesian inference (OGSBI) algorithm, which constructs an off-grid model using linear approximation. However, while OGSBI can effectively achieve off-grid DOA estimation, its performance under coarse grid conditions is unsatisfactory. Therefore, Dai et al. proposed a root-off-grid sparse Bayesian learning (RootSBL) algorithm, which uses discrete grid points as dynamic parameters and iteratively updates the model by solving polynomials. This method reduces computational complexity, improves computational efficiency, and also enhances estimation accuracy under coarse grid conditions. As is well known, the first-order Taylor approximation is sufficiently accurate when the mesh size is small. However, when the mesh size is not small enough to save computational costs, the estimation accuracy will be significantly biased. To address this, a second-order Taylor approximation has been introduced into off-grid DOA estimation. Higher-order Taylor approximations obviously reduce the approximation error, but the computational complexity remains high. Therefore, if a new DOA estimation method can improve both estimation accuracy and computational efficiency, it will greatly increase the practical engineering application value of this method.

[0004] This research, supported by the Shandong Provincial Natural Science Foundation (ZR2017MF024), addresses the aforementioned problems. This patent proposes an improved variational Bayesian sparse learning method for off-grid location estimation, combining real-valued transformation and grid evolution concepts. Compared to traditional compressed sensing methods, this method not only improves estimation accuracy and resolution but also reduces computational complexity and time, demonstrating significant engineering application value. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by proposing an improved variational Bayesian sparse learning off-grid location estimation method. Its key feature is that this patented method combines variational sparse Bayesian learning with grid evolution, adaptively evolving the grid from an initial uniform grid to a non-uniform grid during the iteration process, ultimately achieving source location estimation using a smaller number of grid cells.

[0006] The DOA estimation method of the present invention includes the following steps:

[0007] Step 1: Use a uniform hydrophone array with M array elements and a sampling length of T. The received real data of the m-th array element is represented as follows: Where m = 1, 2, ..., M represents the element number of the hydrophone array, for y m Perform a Hilbert transform to convert the data received by the array elements into complex data.

[0008] Step 2: Arrange the complex data of the M hydrophone array elements into a matrix form. According to the formula Construct the sample covariance matrix, and then vectorize it to obtain y = vec(R), where the superscript "H" indicates that the matrix is ​​conjugate transpose.

[0009] Step 3: Divide the spatial region [-90°, 90°] into N equal parts to obtain an angular grid. Set the overcomplete sparse dictionary set off the grid as Φ(β)=A+Bdiag(β), β=[β1,…,β N ] represents the mesh error vector, array manifold matrix and in The direction vector of the array can be specifically written as: ,and yes The first derivative of the above formula, where d represents the spacing between hydrophone array elements and λ represents the wavelength of the incident signal;

[0010] Step 4: To address the issue of unknown noise accuracy, the formula Φ(β)=[Φ(β) 1 is used. MAn extended, complete sparse dictionary set, in which e m This represents a vector where all elements except the m-th element are zero;

[0011] Step 5: Initialize parameters, ρ = 0.01, δ (0) =0,β (0) =0, Φ (0) (β)=A (0) =A,B (0) =B, i=0 and To set a coarse mesh, set the values ​​of the iteration termination condition τ and Iter. max ;

[0012] Step 6: Based on the variational sparse Bayesian inference, the update formulas for the mean μ and covariance matrix Σ of the sparse signal vector s, as well as the elements of the hyperparameter vector δ, are obtained respectively: μ = ΣΦ H (β)y,Σ=(Φ H (β)Φ(β)+Λ -1 ) -1 , Where Λ = diag(δ), and the hyperparameter δ is represented as δ = [δ1, δ2, ..., δ N+1 ] T , n = 1, 2, ..., N+1, μ n and Σ nn These are the nth element of μ and the nth element of the diagonal of the covariance matrix Σ, respectively.

[0013] Step 7: Update grid points. First, determine if matrix P = Re{B} T B * ⊙ μμ H + Σ Is it reversible? If so, then according to the formula β = P -1 v is used to update the grid points. If it is irreversible, then the formula β is used. i =v i / P ii To update the grid points, where 1 w =vec{I M}, Re{·} denotes taking the real part, ||·||2 denotes calculating the 2-norm of a matrix, tr(·) denotes calculating the trace of a matrix, and '⊙' denotes the Hadamard product. μ Let μ be a vector consisting of the first N elements of μ, and let μ0 be the last element of μ. Σ Let Σ be a matrix consisting of the first N rows and the first N columns of Σ, and let Υ be a vector consisting of the first N columns of the (N+1)th row of Σ.

[0014] Step 8: According to the formula Update grid points if the conditions are met. Then receive this new If the conditions are not met, then the original grid points are retained;

[0015] Step 9: Grid point fission, calculate each element s of s. n The average power is Select the grid point corresponding to the local maximum value among the first M-1 maximum average power Q(n) as the grid point that needs to be split;

[0016] Step 10: Determine whether the left and right intervals of the grid points to be split are less than a certain value, i.e., r(n-1) < θ. ε and r(n) < θ ε If the condition is met, the splitting of the grid point stops; if the condition is not met, the splitting of the grid point is performed. The splitting method is as follows: when the selected splitting grid point is the nth grid point, the position of the newly added grid point after splitting is the midpoint of the left and right grid spacings r(n-1) and r(n) of the grid point. Then, rows corresponding to this grid point with s, β and array streams A and B are added to adapt to the next iteration update.

[0017] Step 11: Determine if the termination condition of the iteration is met. If it is met... Or it reaches the maximum number of iterations, i ≥ Iter max Stop iteration and output. δ (i+1) Otherwise, return to step six and continue with the next iteration.

[0018] Step 12, according to and δ (i+1) Draw the spatial spectrum and find the angle values ​​corresponding to the spatial spectrum extrema. This is the estimated angle of incidence of the target.

[0019] To further improve the computational efficiency of the algorithm, further improvements were made based on the above method, proposing an improved variational Bayesian sparse learning off-grid azimuth estimation method. The method is characterized by: firstly, transforming the vectorized covariance matrix signal in the complex domain to the real domain through a new real-valued transformation; then, combining the ideas of variational sparse Bayesian learning and grid evolution, the grid is adaptively evolved from an initial uniform grid to a non-uniform grid during the iteration process, ultimately achieving the azimuth estimation of the source.

[0020] The DOA estimation method of the present invention includes the following steps:

[0021] Step 1: Use a uniform hydrophone array with M array elements and a sampling length of T. The received real data of the m-th array element is represented as follows: Where m = 1, 2, ..., M represents the element number of the hydrophone array, for y m Perform a Hilbert transform to convert the data received by the array elements into complex data.

[0022] Step 2: Arrange the complex data of the M hydrophone array elements into a matrix form. According to the formula Construct the sample covariance matrix, and then vectorize it to obtain y = vec(R), where the superscript "H" indicates that the matrix is ​​conjugate transpose.

[0023] Step 3: Divide the spatial region [-90°, 90°] into N equal parts to obtain an angular grid. Set the overcomplete sparse dictionary set off the grid as Φ(β)=A+Bdiag(β), β=[β1,…,β N ] represents the mesh error vector, array manifold matrix and in The direction vector of the array can be specifically written as: ,and yes The first derivative of the above formula, where d represents the spacing between hydrophone array elements and λ represents the wavelength of the incident signal;

[0024] Step 4: To address the issue of unknown noise accuracy, the formula Φ(β)=[Φ(β) 1 is used. M An extended, complete sparse dictionary set, in which e m This represents a vector where all elements except the m-th element are zero;

[0025] Step 5: Use the formula y = U H W -1 / 2 y and Φ(β)=U H W -1 / 2 Φ(β) performs a real-valued transformation on the vectorized sample covariance matrix and the overcomplete sparse dictionary set, transforming the array signal from the complex domain to the real domain, where... W -1 / 2 W -1 The Hermitian square root, U, is obtained by applying W. -1 / 2 A is obtained through singular value decomposition, i.e., UΛV H =svd(W -1 / 2 A);

[0026] Step 6: Initialize parameters, ρ = 0.01, δ (0) =0,β (0) =0, Φ (0) (β)=A (0) =A,B (0) =B, i=0 and To set a coarse mesh, set the values ​​of the iteration termination condition τ and Iter. max ;

[0027] Step 7: Based on the variational sparse Bayesian inference, the update formulas for the mean μ and covariance matrix Σ of the sparse signal vector s, as well as the elements of the hyperparameter vector δ, are obtained respectively: μ = ΣΦ H (β)y,Σ=(Φ H (β)Φ(β)+Λ -1 ) -1 , Where Λ = diag(δ), and the hyperparameter δ is represented as δ = [δ1, δ2, ..., δ N+1 ] T , n = 1, 2, ..., N+1, μ n and Σ nn These are the nth element of μ and the nth element of the diagonal of the covariance matrix Σ, respectively.

[0028] Step 8: Update grid points. First, determine if matrix P = Re{(B U ) T B U * ⊙ μμ H + Σ Is it reversible? If so, then according to the formula β = P -1 v is used to update the grid points. If it is irreversible, then the formula β is used. i =v i / P ii To update the grid points, where A U =U H W -1 / 2 A, B U =U H W -1 / 2 B, 1 U =U H W -1 / 2 vec{I M}, Re{·} denotes taking the real part, ||·||2 denotes calculating the 2-norm of a matrix, tr(·) denotes calculating the trace of a matrix, and '⊙' denotes the Hadamard product. μ Let μ be a vector consisting of the first N elements of μ, and let μ0 be the last element of μ. Σ Let Σ be a matrix consisting of the first N rows and the first N columns of Σ, and let Υ be a vector consisting of the first N columns of the (N+1)th row of Σ.

[0029] Step 9: According to the formula Update grid points if the conditions are met. Then receive this new If the conditions are not met, then the original grid points are retained;

[0030] Step 10: Grid point fission, calculate each element s of s. n The average power is Select the grid point corresponding to the local maximum value among the first M-1 maximum average power Q(n) as the grid point that needs to be split;

[0031] Step 11: Determine whether the left and right intervals of the grid points to be split are less than a certain value, i.e., r(n-1) < θ. ε and r(n) < θ ε If the condition is met, the splitting of the grid point stops; if the condition is not met, the splitting of the grid point is performed. The splitting method is as follows: when the selected splitting grid point is the nth grid point, the position of the newly added grid point after splitting is the midpoint of the left and right grid spacings r(n-1) and r(n) of the grid point. Then, rows corresponding to this grid point with s, β and array streams A and B are added to adapt to the next iteration update.

[0032] Step 12: Determine if the termination condition of the iteration is met. If it is met... Or it reaches the maximum number of iterations, i ≥ Iter max Stop iteration and output. δ (i+1) Otherwise, return to step seven and continue with the next iteration.

[0033] Step Thirteen, according to and δ (i+1) Draw the spatial spectrum and find the angle values ​​corresponding to the spatial spectrum extrema. This is the estimated angle of incidence of the target.

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

[0035] (1) The method of this invention combines the ideas of variational sparse Bayesian learning and grid evolution. During the iteration process, the grid adaptively evolves from an initial uniform grid to a non-uniform grid. The evolution process includes grid update and grid fission. During the alternating iteration of grid update and grid fission, the grid points gradually approach the true source location. This process not only improves the accuracy of azimuth estimation, but also significantly improves the resolution success rate due to the adaptive grid evolution. Furthermore, this evolution process effectively controls the number of grids, thus reducing the computational load of iteration and increasing computational efficiency.

[0036] (2) This invention transforms the original complex sparse signal reconstruction problem into a real-value problem by utilizing the special properties of the virtual steering vector of the linear array through real-value transformation, thereby reducing the amount of computation and computational cost, and thus having higher practical engineering application value. Attached Figure Description

[0037] Figure 1 This is a schematic diagram of the grid point evolution process of the signal processing method of this patent;

[0038] Figure 2 This describes the grid point update and fission process of the signal processing method in this patent.

[0039] Figure 3 This is the spatial spectrum of the signal processing method of this patent;

[0040] Figure 4 The curve showing the relationship between the estimation performance and signal-to-noise ratio of the signal processing method of this patent is shown.

[0041] Figure 5 The curve showing the relationship between the estimation performance of the signal processing method of this patent and the number of snapshots is shown.

[0042] Figure 6 This is a curve showing the relationship between the resolution success rate of the signal processing method of this patent and the distance between the signal sources at high signal-to-noise ratio.

[0043] Figure 7 This is a curve showing the relationship between the resolution success rate of the signal processing method of this patent and the distance between the signal sources at low signal-to-noise ratio.

[0044] Figure 8 This is a curve showing the relationship between the computation time and the number of snapshots for the signal processing method of this patent. Detailed Implementation

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

[0046] Example 1: Figure 1A schematic diagram of the grid point evolution process of the signal processing method in this invention is given. First, a uniform linear array with M=8 array elements is constructed, with an array spacing d equal to a half wavelength of 0.03m. Two incoherent signals are incident on the receiving array from directions θ=-12.234° and 4.565°, respectively. The signal-to-noise ratio is set to 10dB, and the number of snapshots is T=100. Under these conditions, the specific implementation process is as follows:

[0047] Step 1: Use a uniform hydrophone array with M=8 array elements and a sampling length of T=100. The received real data of the m-th array element is represented as follows: Where m = 1, 2, ..., 8 represents the element number of the hydrophone array, and for y m Perform a Hilbert transform to convert the data received by the array elements into complex data.

[0048] Step 2: Arrange the complex data of M=8 hydrophone array elements into a matrix form. According to the formula Construct the sample covariance matrix, and then vectorize it to obtain y = vec(R), where the superscript "H" indicates that the matrix is ​​conjugate transpose.

[0049] Step 3: Divide the spatial region [-90°, 90°] evenly into N = 10 parts to obtain an angular grid. Set the overcomplete sparse dictionary set off the grid as Φ(β)=A+Bdiag(β), β=[β1,…,β 10 ] represents the mesh error vector, array manifold matrix and in The direction vector of the array can be specifically written as: ,and yes The first derivative of the above formula, where d = 0.03 represents the spacing between hydrophone array elements and λ = 0.06 represents the wavelength of the incident signal;

[0050] Step 4: To address the issue of unknown noise accuracy, the formula Φ(β)=[Φ(β) 1 is used. M An extended, complete sparse dictionary set, in which e m This represents a vector where all elements except the m-th element are zero, except for the m-th element which is 1.

[0051] Step 5: Use the formula y = U H W -1 / 2 y and Φ(β)=U H W -1 / 2Φ(β) performs a real-valued transformation on the vectorized sample covariance matrix and the overcomplete sparse dictionary set, transforming the array signal from the complex domain to the real domain, where... W -1 / 2 W -1 The Hermitian square root, U, is obtained by applying W. -1 / 2 A is obtained through singular value decomposition, i.e., UΛV H =svd(W -1 / 2 A);

[0052] Step 6: Initialize parameters, ρ = 0.01, δ (0) =0,β (0) =0, Φ (0) (β)=A (0) =A,B (0) =B, i=0 and To set a coarse mesh, the iteration termination condition is set to τ = 10. -3 and Iter max =500;

[0053] Step 7: Based on the variational sparse Bayesian inference, the update formulas for the mean μ and covariance matrix Σ of the sparse signal vector s, as well as the elements of the hyperparameter vector δ, are obtained respectively: μ = ΣΦ H (β)y,Σ=(Φ H (β)Φ(β)+Λ -1 ) -1 , Where Λ = diag(δ), and the hyperparameter δ is represented as δ = [δ1, δ2, ..., δ N+1 ] T , μ n and Σ nn These are the nth element of μ and the nth element of the diagonal of the covariance matrix Σ, respectively.

[0054] Step 8: Update grid points. First, determine if matrix P = Re{(B U ) T B U * ⊙ μμ H + Σ Is it reversible? If so, then according to the formula β = P -1 v is used to update the grid points. If it is irreversible, then the formula β is used. i =v i / P ii To update the grid points, where A U =U H W -1 / 2 A, BU =U H W -1 / 2 B, 1 U =U H W -1 / 2 vec{I M}, Re{·} denotes taking the real part, ||·||2 denotes calculating the 2-norm of a matrix, tr(·) denotes calculating the trace of a matrix, and '⊙' denotes the Hadamard product. μ Let μ be a vector consisting of the first N elements of μ, and let μ0 be the last element of μ. Σ Let Σ be a matrix consisting of the first N rows and the first N columns, and let Υ be a vector consisting of the first N columns of the (N+1)th row of Σ.

[0055] Step 9: According to the formula Update grid points if the conditions are met. Then receive this new If the conditions are not met, then the original grid points are retained;

[0056] Step 10: Grid point fission, calculate each element s of s. n The average power is Select the grid point corresponding to the local maximum value among the first M-1 maximum average power Q(n) as the grid point that needs to be split;

[0057] Step 11: Determine whether the left and right intervals of the grid points to be split are less than a certain value, i.e., r(n-1) < θ. ε =1 and r(n) < θ ε When = 1, if the condition is met, the grid point stops splitting; if not, the grid point splits. The splitting method is as follows: when the selected splitting grid point is the nth grid point, the position of the newly added grid point after splitting is the midpoint of the left and right grid spacings r(n-1) and r(n) of the grid point. Then, rows corresponding to this grid point and the array streams A and B are added to adapt to the next iteration update.

[0058] Step 12: Determine if the termination condition of the iteration is met. If it is met... Or it reaches the maximum number of iterations, i ≥ Iter max =500, stop iteration, output δ (i+1) Otherwise, return to step seven and continue with the next iteration.

[0059] Step Thirteen, according to and δ (i+1) Draw the spatial spectrum and find the angle values ​​corresponding to the spatial spectrum extrema. This is the estimated angle of incidence of the target.

[0060] Figure 2 The graph illustrates the mesh update and fission process during iteration. The x-axis represents the iteration number, and the y-axis represents the position of the mesh points in each iteration. It can be seen that the mesh points are uniformly distributed from [-90°, 90°], with an initial mesh interval of 20°. As iterations proceed, the updated and fissuring mesh points increasingly approximate the actual source location and are densely distributed in the vicinity of that location. After approximately 48 iterations, mesh fission ceases, and only an update process with 46 meshes is performed. Finally, iteration stops after the 101st iteration, resulting in the final spatial spectrum. Figure 3 As shown. By Figure 3 It can be seen that the method of this patent can estimate the specific location of two incoherent signals relatively accurately, with good estimation accuracy and spatial resolution.

[0061] The second embodiment: This study investigates the relationship between the estimation performance of the signal processing method of this patent and the signal-to-noise ratio (SNR) and the number of snapshots, and displays the root mean square error (RMSE) as a function of the SNR. The resulting graph is shown below. Figure 4 As shown; and displaying the root mean square error as a function of the number of snapshots, as shown. Figure 5 As shown. The conditions for applying the algorithm in this invention are as follows:

[0062] We used eight uniform linear arrays as an example, with an array spacing d equal to half the wavelength of 0.03m. Two incoherent signals were incident on the receiving array from directions θ = -12.234° and 4.565°. Initially, the number of snapshots T = 100 was set. The signal-to-noise ratio (SNR) was changed from 0dB, increasing in 2dB steps to 16dB. Then, keeping the SNR constant at 10dB, the number of snapshots was changed from 20, increasing in 20 steps to 200. 500 independent Monte Carlo experiments were performed, and simulations were conducted using MATLAB software. The simulation results are as follows: Figure 4 , 5 As shown.

[0063] from Figure 4 As can be seen, comparing the L1-SVD method, OGSBI method, RootSBL method, Variational Sparse Bayesian Learning (VBL) method, Method 1 of this patent, and Method 2 of this patent, the root mean square error (RMSE) curves of all six methods gradually decrease with increasing signal-to-noise ratio (SNR). The L1-SVD method has the highest RMSE, indicating a larger error. Method 2 of this patent outperforms the other methods within this SNR range, exhibiting the lowest RMSE and the smallest error. At a low SNR of 0 dB, the RMSE of Method 2 of this patent is only 0.4143. This demonstrates that Method 2 of this patent has better estimation performance, especially under low SNR conditions, where its advantages are more pronounced.

[0064] Depend on Figure 5It can be seen that the root mean square error (RMSE) of all six methods decreases with the increase of the number of snapshots, and gradually stabilizes after the number of snapshots reaches 80. The RMSE curve of Method 2 of this patent is lower than that of the other five methods within this range of snapshot numbers, and the RMSE is only 0.2585 when the number of snapshots is 20. This indicates that Method 2 of this patent has the best estimation performance compared to the other methods under different snapshot number conditions.

[0065] The third embodiment: This study investigates the relationship between the resolution capability of the signal processing method of this patent and the incident source angular spacing, and displays the resolution success rate as a function of the incident angular spacing under different signal-to-noise ratios, as shown in the following figure. Figure 6 , 7 As shown. The conditions for applying the algorithm in this invention are as follows:

[0066] We used eight uniform linear arrays as an example, with an array spacing d equal to half a wavelength of 0.03 m. Two incoherent signals were incident on the receiving array from directions θ1 and θ2 = θ1 + Δ. θ1 was fixed, and Δ varied from 1° to 10° in 1° increments. The number of snapshots T = 100, the high signal-to-noise ratio (SNR) was 10 dB, and the low SNR was 0 dB. The position of the first signal source was fixed, and 500 independent Monte Carlo experiments were performed. The simulations were conducted using MATLAB software, and the simulation results were observed. The simulation results under the high SNR condition are as follows: Figure 6 As shown, the simulation results under low signal-to-noise ratio conditions are as follows: Figure 7 As shown.

[0067] Depend on Figure 6 It can be seen that under high signal-to-noise ratio (SNR) conditions, for different incident source grid spacings, Method 2 of this patent can achieve a 100% resolution success rate at 4°. In contrast, the OGSBI and RootSBL methods require an incident source grid spacing of 5° or higher to achieve a 100% resolution success rate. Under the same conditions, the L1-SVD method still cannot achieve 100% resolution success even if the incident source grid spacing is increased to 10°. Therefore, it is evident that Method 2 of this patent has the strongest resolution capability under high SNR conditions.

[0068] pass Figure 7 It is evident that, under low signal-to-noise ratio conditions, although the resolution success rate of all six methods cannot reach 100%, the resolution performance advantage of Method 2 in this patent becomes increasingly apparent as the grid spacing of the incident signal source increases. Combined with... Figure 6 and 7 This demonstrates that regardless of whether the signal-to-noise ratio is high or low, Method 2 of this patent has good resolution capabilities, especially at high signal-to-noise ratios, where it can distinguish between two relatively close signals.

[0069] Fourth embodiment: The relationship between the computation time and the number of snapshots in the signal processing method of this patent is studied, and the results of the computation time changing with the number of snapshots are displayed, as shown in the figure below. Figure 8 As shown. The conditions for applying the algorithm in this invention are as follows:

[0070] We used eight uniform linear arrays as an example, with an array spacing d equal to half the wavelength of 0.03m. Two incoherent signals were incident on the receiving array from directions θ = -12.234° and 4.565°, respectively. Maintaining a constant signal-to-noise ratio of 10dB, the number of snapshots was changed from 20 to 200 in increments of 20. 500 independent Monte Carlo experiments were performed, and simulations were conducted using MATLAB software. The simulation results are shown below. Figure 8 As shown.

[0071] Figure 8 The graph shows the computation time as a function of the number of snapshots. It can be seen that within this range of snapshot counts, the L1-SVD method has the longest computation time, while the proposed method 2 has the shortest computation time and the highest efficiency. Furthermore, the computation time of the other five methods increases with the number of snapshots, but the computational efficiency of the proposed method 2 is less affected by changes in the number of snapshots, exhibiting a stable trend.

[0072] The specific examples described in this patent 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. An improved variational Bayesian sparse learning method for off-grid location estimation, characterized in that: The orientation estimation method includes the following steps: Step 1: Use an array element count of The sampling length is A uniform hydrophone array, the array of which is the first The received real data of each array element is represented as follows: ,in Indicates the element number of the hydrophone array, for Perform a Hilbert transform to convert the data received by the array elements into complex data. ; Step 2, The complex data of the hydrophone array elements are arranged in matrix form. According to the formula Construct the sample covariance matrix, and then vectorize it to obtain... superscript This indicates the operation of conjugate transpose on a matrix; Step 3: Divide the spatial area Evenly divided into The angle grid is obtained. Set the overcomplete sparse dictionary set off from the grid as , Represents the mesh error vector, array manifold matrix and ,in , The direction vector of the array can be specifically written as: ,and , yes The first derivative of , in the above formula, Indicates the spacing between hydrophone array elements. Indicates the wavelength of the incident signal; Step 4: To address the issue of unknown noise accuracy, the formula is used. An extended, complete sparse dictionary set, in which , Then it means except for the first A vector in which 1 element is 1 and all other elements are zero; Step 5: Using the formula and A real-valued transformation is performed on the vectorized sample covariance matrix and the overcomplete sparse dictionary set to transform the array signal from the complex domain to the real domain. , express Hermitian square root, Then through the Singular value decomposition yields, i.e. ; Step 6: Initialize parameters. , , , , , as well as Set the value of the iteration termination condition for the coarse grid. and ; Step 7: Based on the variational sparse Bayesian inference, obtain the sparse signal vectors respectively. mean Covariance Matrix and hyperparameter vectors The update formulas for each element in the formula are as follows: , , ,in hyperparameters Represented as , , , and They are respectively The Each element and its covariance matrix The diagonal of Each element value; Step 8: Update grid points. First, determine the matrix. Is it reversible? If so, then according to the formula... Update the grid points; if it's irreversible, then follow the formula. To update the grid points, where , , , , Indicates taking the real part, This indicates finding the 2-norm of a matrix. This indicates finding the trace of a matrix. 'Represents the Hadamard product, for The former A vector consisting of n elements for The last element, For the reason The former line and front A matrix composed of columns Indicates by The Before the journey A vector consisting of column elements; Step 9: According to the formula Update grid points if the conditions are met. Then receive this new If the conditions are not met, then the original grid points are retained. Step 10: Grid point fission, calculation Each element The average power is ,in front Maximum average power Select the grid point location corresponding to the local maximum value as the grid point that needs to be split; Step 11: Determine whether the left and right spacing of the grid points to be split is less than a certain value, i.e. and If the condition is met, the fission of that grid point stops; otherwise, the fission of that grid point begins. The fission method is as follows: when the selected fission grid point is the th... If there are 1 grid point, then the position of the newly added grid point after fission is the left and right grid spacing of that grid point. and The midpoint, then add with , and array popularity and The row corresponding to this grid point is used to adapt to the next iteration update; Step 12: Determine if the termination condition of the iteration is met. If it is met... Or, the maximum number of iterations is reached. Stop iteration and output. , Otherwise, return to step seven and continue with the next iteration. Step Thirteen, according to and Draw the spatial spectrum and find the angle values ​​corresponding to the spatial spectrum extrema. This is the estimated angle of incidence of the target.

Images