A robust DOA estimation method based on sparse Bayesian learning

Abstract

The application discloses a kind of robust DOA estimation methods based on sparse bayesian learning, belong to array signal processing field, the method of this application first introduces position error parameter, and determines the prior distribution of grid error and element position error, secondly with the aid of sparse bayesian learning model, establish joint probability density distribution function, finally with the aid of maximum expectation algorithm to each unknown parameter is iterated, obtains spatial spectrum diagram, estimates the direction of incident signal.This application method can better azimuth estimation accuracy and azimuth resolution, still have good performance under the condition of less snapshot, low signal-to-noise ratio, and the application method has good robustness to the array system that exists element position error, has greater engineering application value.

Description

A robust DOA estimation method based on sparse Bayesian learning

[0001] The present invention relates to the field of array signal processing, and in particular to a robust DOA estimation method based on sparse Bayesian learning. Specifically, it is a robust DOA estimation method based on sparse Bayesian learning under the condition of array element position error. Background Technology

[0002] Direction of Arrival (DOA) estimation has always been a hot topic in array signal processing. Its principle involves analyzing the characteristic information of signals received by sensor arrays using various methods, and it has wide applications in radar, sonar, microphones, and medicine. Typical subspace-based DOA estimation methods, such as the MUSIC method and the ESPRIT method, limit the maximum number of resolvable targets to the number of receiving array elements. Coprime arrays can overcome the limitation of the physical number of array elements on the maximum number of resolvable signals. Furthermore, coprime arrays have a larger array aperture than uniform linear arrays with the same number of elements, thus exhibiting superior DOA estimation performance.

[0003] For coprime array models, researchers have proposed methods such as Spatial Smoothing MUSIC (SS-MUSIC) and the Joint ESPRIT method. The SS-MUSIC method can detect more signal sources than sensors while retaining its high-resolution performance. However, the application of spatial smoothing techniques leads to the loss of half of the continuous DOFs, resulting in a significant decrease in array detection performance. Furthermore, when using differential operations on coprime arrays to form augmented virtual arrays, non-continuous virtual array elements are ignored during smoothing, and the information of the virtual array is not fully utilized. The Joint ESPRIT method decomposes the coprime array into two uniform linear subarrays and uses the ESPRIT method to estimate the direction of arrival (DOA) separately. The method finds the coincident estimates from the two subarray estimation results to determine the azimuth of the incident signal. This method has a much lower complexity than the SS-MUSIC method. However, by decomposing the coprime array into two subarrays for separate calculations, the number of estimable targets is reduced by at least 50% compared to traditional uniform arrays.

[0004] In recent years, with the deepening research on compressed sensing and sparse reconstruction methods, researchers have found that this method can be used for virtual degrees of freedom extended by coprime arrays. Therefore, compared with subspace-based algorithms, compressed sensing algorithms have better DOA estimation performance, leading to the proposal of numerous DOA estimation methods based on spatial sparsity. Researchers proposed a DOA estimation method based on Off-Grid Sparse Bayesian Learning (OGSBL) models. This method introduces an offset parameter and uses Sparse Bayesian Learning (SBL) to solve for the estimated offset, improving the orientation estimation performance in cases where the incident signal is off-grid. Subsequently, Dai et al. proposed the Root-OGSBL method, reducing the computational cost of the OGSBL method.

[0005] However, the aforementioned methods all ignore the influence of various error factors. In practical engineering, element position errors are inevitable, which greatly affect the estimation performance of various methods and may even cause them to fail. Therefore, with the support of the Shandong Provincial Natural Science Foundation (ZR2017MF024), this problem was studied, and a high-precision DOA estimation method based on sparse Bayesian learning was explored. This patent proposes a robust DOA estimation method based on sparse Bayesian learning, which can still achieve good orientation estimation performance even with a small number of snapshots and a low signal-to-noise ratio, greatly improving the practical engineering application value of this method. Summary of the Invention

[0006] To address the problem of performance degradation in various estimation methods due to element position errors in practical engineering, this invention proposes a robust DOA estimation method based on sparse Bayesian learning. The method is characterized by: firstly, introducing position error parameters and determining the prior distributions of grid errors and element position errors; secondly, establishing a joint probability density function using a sparse Bayesian learning model; and finally, iterating through the unknown parameters using the expectation-maximization algorithm to obtain the spatial spectrum. To implement the above DOA estimation method, this invention provides a robust DOA estimation method based on sparse Bayesian learning for the presence of element position errors, the process of which includes the following steps:

[0007] Step 1: The distance between the M array elements and the reference array element is d = [d1, ..., d2]. m ,…,d M ] T The positional error of each array element is Δ P =[Δ P1 ,…,Δ Pm ,…,Δ PM ] TThe received data is represented as Y(t)=A(θ,Δ P Y(t) + N(t), t = 1, 2, ..., T, where T represents the number of snapshots, Y(t) represents the data received by the array, s(t) represents the transmitted signal, and N(t) represents the noise signal received by the array, all having a mean of 0 and a covariance of σ. 2 A narrowband Gaussian distribution, A(θ,Δ P The manifold matrix of the array is represented as A(θ,Δ). P )=[a(θ1,Δ P ),…,a(θ k ,Δ P ),…,a(θ K ,Δ P [), where K is the number of signal sources.

[0008] Step 2: Divide the spatial angle range [-90°, 90°] evenly into N parts to obtain the grid set. The sparse signal model of the array received data Y(t) is established as Y(t) = Φ(Δt). P ,δ)X(t)+N(t), where X(t) is the zero-spread of the original signal s(t), which has values ​​only at grid points close to the incident angle and is zero at all other locations. The `diag(·)` operation expands a vector into a diagonal matrix. δ represents the grid error.

[0009] Step 3: Initialize the hyperparameters b, c, and e, setting them to be less than 5 × 10. -2 The determined value, the parameter signal precision γ that needs to be updated during initialization is 0. N×1 Noise accuracy α n Satisfy α n ∈[10 -2 [1], mesh error δ=0 N×1 Array element position error Δ P =0 M×1 The accuracy of the array element position error ρ = 0 M×1 Initialize the loop iteration factor l = 1;

[0010] Step 4: Calculate the mean and covariance of the posterior probability of the sparse signal X, which follow a pattern with a mean of μ. x The covariance is Σ x The Gaussian distribution of , where the covariance matrix Σ x =(α n Φ H (ΔP ,δ)Φ(Δ P ,δ)+Λ -1 ) -1 , The superscript "H" indicates the operation of taking the conjugate transpose of a matrix, α n Indicates noise precision, Λ = diag(γ), mean.

[0011] Step 5: Update the value of the signal precision γ, resulting in the expression: Where c is the hyperparameter of the gamma distribution, 1≤n≤N, μ x(:,t) Represents the mean matrix μ x The t-th column, Σ x(n,n) This represents the nth row and nth column of the covariance matrix;

[0012] Step 6: Update noise accuracy α n The value of is used to obtain the expression . Among them, Y t Let Y represent the t-th column of the received data Y, b be the hyperparameter of the gamma distribution, Re{·} denote the operation of taking the real part, and Tr{·} denote the operation of taking the trace of the matrix;

[0013] Step 7: Calculate the element position error Δ P The expression for the position error of the m-th array element is: in, m = 1, 2, ..., M This indicates that the m-th position is 1, and the remaining positions are 0;

[0014] Step 8: Calculate the positional error precision ρ of the array elements. The calculation expression for the m-th element is as follows: m = 1, ..., M, where e is a hyperparameter;

[0015] Step 9: Calculate parameters Where l represents the number of iterations. The mean matrix calculated in step 4. If the parameter κ satisfies the error precision ε or the maximum number of iterations maxIter, proceed to step 10. If neither condition is satisfied, then... l = l + 1, then re-enter step 4 and iterate again;

[0016] Step 10: Update the mesh error δ, using Estimate the grid error vector. in, This represents the grid spacing, where 2 ≤ n ≤ N. Represents the Hadamard product;

[0017] Step 11: Calculate the spatial spectrum This represents the row mean vector, and the superscript "*" indicates the conjugate operation of the vector;

[0018] Step 12: Update the spatial grid points using the grid error calculated in Step 10, i.e. At the same time, it corresponds one-to-one with the spatial spectrum in step 11, and the angle corresponding to the peak of the spatial spectrum is the estimated direction of arrival of the K signals.

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

[0020] (1) Compared with traditional uniform linear arrays, it can break through the limitation of the number of physical array elements on the maximum number of resolvable signals and can realize target detection with more than the number of array elements; at the same time, coprime arrays have a larger array aperture than uniform linear arrays with the same number of array elements, which improves the degree of freedom and has better DOA estimation performance.

[0021] (2) Compared with other methods, the method of the present invention can still obtain robust DOA estimation performance when the number of snapshots is small and the signal-to-noise ratio is low. The estimation accuracy is better than that of conventional methods. At the same time, the method of the present invention also has a high angle resolution capability.

[0022] (3) The method of the present invention effectively solves the problem of array element position error, greatly improves the robustness of the coprime array-based orientation estimation method, and has great application value in practical engineering. Attached Figure Description

[0023] Figure 1 shows the power spectrum of the simulation experiment of the method of this patent and other methods;

[0024] Figure 2 shows the relationship between the root mean square error and the signal-to-noise ratio of the processing method of this patent.

[0025] Figure 3 shows the relationship between the root mean square error and the number of snapshots in the processing method of this patent.

[0026] Figure 4 shows the relationship between the success rate of resolution and the signal-to-noise ratio of the processing method of this patent.

[0027] Figure 5 shows the relationship between the success rate of resolution and the number of snapshots in the processing method of this patent. Detailed Implementation

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

[0029] In the method of this invention, we constructed a coprime array with M=7 array elements. Signals were incident on the receiving array from two directions, -27.32° and 17.75°, respectively, and the number of snapshots was T=500. Under these conditions, the specific implementation process is as follows:

[0030] Step 1: The distance between the 7 array elements and the reference array element is d = [0, 3l, 5l, 6l, 9l, 10l, 12l]. T Where l represents the unit spacing, with a value of l = λ / 2, λ represents the wavelength, and the position error of each array element is Δ. P =[-0.12l,-0.19l,0.22l,0.09l,-0.11l,0.06l,-0.22l] T The received data is represented as Y(t)=A(θ,Δ P Let Y(t) + N(t), t = 1, 2, ..., 500, where Y(t) represents the array received data, s(t) represents the transmitted signal, and N(t) represents the array received noise signal, all having a mean of 0 and a covariance of σ. 2 A narrowband Gaussian distribution, A(θ,Δ P The manifold matrix of the array is represented as A(θ,Δ). P )=[a(θ1,Δ P ),…,a(θ k ,Δ P ),…,a(θ2,Δ P [), where K is the number of signal sources.

[0031] Step 2: Divide the spatial angle range [-90°, 90°] into 61 grid points with a step size of 3°, resulting in the following grid set. The sparse signal model of the array received data Y(t) is established as Y(t) = Φ(Δt). P ,δ)X(t)+N(t), where X(t) is the zero-extension of the original signal s(t), which has a value only near the incident angle and is zero at other locations. diag(·) represents the operation of converting a vector into a diagonal matrix. δ represents the grid error.

[0032] Step 3: Initialize the hyperparameters b, c, and e, setting b = c = e = 10 -3 Set the maximum number of iterations (maxIter) to 300 and the error precision (ε) to 10. -3 Initialization requires updating the signal precision γ = 0. N×1Noise accuracy α n =0.1, mesh error δ=0 N×1 Array element position error Δ P =0 M×1 The accuracy of the array element position error ρ = 0 M×1 Initialize the loop iteration factor l = 1;

[0033] Step 4: Calculate the mean and covariance of the posterior probability of the sparse signal X, which follow a pattern with a mean of μ. x The covariance is Σ x The Gaussian distribution of , where the covariance matrix Σ x =(α n Φ H (Δ P ,δ)Φ(Δ P ,δ)+Λ -1 ) -1 , The superscript "H" indicates the operation of taking the conjugate transpose of a matrix, α n The noise level is represented by Λ = diag(γ), and the mean is μ. x =α n Σ x Φ H (Δ P ,δ)Y,

[0034] Step 5: Update the value of the signal precision γ, resulting in the expression: Where c is the hyperparameter of the gamma distribution, 1≤n≤61, μ x(:,t) Represents the mean matrix μ x The t-th column, Σ x(n,n) Represents the covariance matrix Σ x The nth row and nth column;

[0035] Step 6: Update noise accuracy α n The value of is used to obtain the expression . Among them, Y t Let Y represent the t-th column of the received data Y, b be the hyperparameter of the gamma distribution, Re{·} denote the operation of taking the real part, and Tr{·} denote the operation of taking the trace of the matrix;

[0036] Step 7: Calculate the element position error Δ P The expression for the position error of the m-th array element is: m = 1, 2, ..., 7, where diag(·) represents extracting the diagonal elements of a matrix or transforming a vector into a diagonal matrix. This indicates that the m-th position is 1, and the remaining positions are 0;

[0037] Step 8: The expression for the array element position error accuracy ρ is as follows: m = 1, 2, ..., 7, where e is a hyperparameter;

[0038] Step 9: Calculate parameters Where l represents the number of iterations. The mean matrix calculated in step 4. If the parameter κ satisfies the error precision ε or the maximum number of iterations maxIter, proceed to step 10. If neither condition is satisfied, then... l = l + 1, then re-enter step 4 and iterate again;

[0039] Step 10: Update the mesh error δ using the formula Estimate the grid error vector. in, Represents the Hadamard product;

[0040] Step 11: Calculate the spatial spectrum This represents the row mean vector, and the superscript "*" indicates the conjugate operation of the vector;

[0041] Step 12: Update the spatial grid points using the grid error calculated in Step 10, i.e. Simultaneously, the spatial spectrum corresponds one-to-one with the spatial spectrum in step 11, and the angle corresponding to the peak value of the spatial spectrum is the estimated direction of arrival of the two signals. Combining the above steps, we can obtain the spatial spectrum images of this patent method and other methods, as shown in Figure 1.

[0042] As can be seen from Figure 1, both the method of this patent and other methods can accurately estimate the incident direction of the signal. Through the spatial spectrum image, we can see that the curves of the method of this patent and the MUSIC method are smooth in the direction independent of the incident angle of the signal. However, the power of the method of this patent in the direction independent of the incident angle of the signal is much lower than that of the MUSIC method. The SunFG method and the OGSBI method have large fluctuations in the direction independent of the incident angle of the signal, which may interfere with the estimation results. Therefore, the method of this patent has a more robust estimation performance. Through some detailed magnification, the method of this patent is closest to the actual incident angle of the signal. Therefore, the method of this patent has better DOA estimation accuracy.

[0043] We used a coprime array under the above conditions as an example to explore the relationship between the number of snapshots and the root mean square error of the angle. The signal was incident on the receiving array from a direction of -27.32°, the number of snapshots T = 500, and the signal-to-noise ratio was changed from -10dB to 10dB in steps of 2dB. 200 independent Monte Carlo experiments were conducted, and the signal-to-noise ratio versus root mean square error curve was obtained by MATLAB software, as shown in Figure 2.

[0044] As can be seen from Figure 2, the root mean square error of the four methods decreases with the increase of signal-to-noise ratio. Throughout the entire signal-to-noise ratio range, the root mean square error of the estimated angle of the present invention is always smaller than that of the other three methods. Obviously, the DOA estimation accuracy of the present invention is better than that of the other three methods.

[0045] We used a coprime array under the above conditions as an example to explore the relationship between the number of snapshots and the root mean square error of the angle. The signal was incident on the receiving array from a direction of -27.32° with a signal-to-noise ratio of 0dB. The number of snapshots was changed from 50 to 500 in step size. 200 independent Monte Carlo experiments were conducted and simulated using MATLAB software. The curve of the number of snapshots versus the root mean square error of the angle is shown in Figure 3.

[0046] As can be seen from Figure 3, the root mean square error of the four methods decreases with the increase of the number of snapshots. Throughout the entire snapshot number range, the root mean square error of the estimated angle of the present invention is always smaller than the root mean square error of the angle estimation of the other three methods. Obviously, the DOA estimation accuracy of the present invention is better than that of the other three methods.

[0047] We used a coprime array under the above conditions as an example to explore the relationship between signal-to-noise ratio (SNR) and angular resolution. The signal was incident on the receiving array from a direction of -27.32°, with a snapshot number T = 500. The SNR was changed from -10dB to 10dB in steps of 2dB, and 200 independent Monte Carlo experiments were conducted. The simulation was performed using MATLAB software. If the difference between the estimated angle value and the actual incident angle was less than 0.2°, the resolution was considered successful; otherwise, it was considered a failure. The relationship between SNR and angular resolution is shown in Figure 4.

[0048] As can be seen from Figure 4, the angle resolution probability estimated by the four methods increases with the increase of signal-to-noise ratio. Throughout the entire signal-to-noise ratio range, the angle resolution probability of the method of this invention is always higher than that of the other three methods. Obviously, the DOA estimation resolution capability of the method of this invention is better than that of the other three methods.

[0049] We used a coprime array under the above conditions as an example to explore the relationship between the number of snapshots and the angular resolution. The signal was incident on the receiving array from a direction of -27.32° with a signal-to-noise ratio of 0dB. The number of snapshots was changed from 50 to 500 in step size, and 200 independent Monte Carlo experiments were conducted. The simulation was performed using MATLAB software. If the difference between the estimated angle value and the actual incident angle is less than 0.2°, the resolution is considered successful; otherwise, it is considered a failure. The relationship between the number of snapshots and the angular resolution is shown in Figure 5.

[0050] As can be seen from Figure 5, the angle resolution probability estimated by the four methods increases with the increase of the number of snapshots. Throughout the entire snapshot number range, the angle resolution probability of the method of this invention is always higher than that of the other three methods. Obviously, the DOA estimation resolution capability of the method of this invention is better than that of the other three methods.

[0051] 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 robust DOA estimation method based on sparse Bayesian learning, characterized in that: The DOA estimation method includes the following steps: Step 1: The distance interval between the M array elements and the reference array element is... The positional error of each array element is The received data is represented as ， ,in, Indicates the number of snapshots. This indicates that the array is receiving data. Indicates the transmission of a signal. The noise signal received by the array has a mean of 0 and a covariance of . Narrow-band Gaussian distribution, The manifold matrix of the array is represented as ， The number of signal sources. Step 2: Determine the airspace angle range Evenly divided into The resulting grid set is... Establish an array to receive data The sparse signal model is ,in, It is the original signal The zero expansion only has values ​​at grid points close to the incident angle; all other locations are zero. ， This represents the operation of expanding a vector into a diagonal matrix. ， ， Indicates grid error, ， Step 3: Adjust hyperparameters ， ， Perform initialization, setting the maximum number of iterations (maxIter) and the error precision. Initialize parameters that need to be updated, including signal precision. Noise accuracy Grid error Array element position error sum of array element position error accuracy Set loop parameters Step 4: Calculate the mean and covariance of the posterior probability of the sparse signal X, which follow a mean of . covariance is The Gaussian distribution, where the covariance matrix is ​​. ， The superscript "H" indicates the operation of taking the conjugate transpose of a matrix. Indicates noise accuracy. mean , Step 5: Update signal accuracy The value of is used to obtain the expression . Where c is the hyperparameter of the gamma distribution. ， Represents the mean matrix The List, The first element of the covariance matrix represents the first element of the covariance matrix. Line number Column; Step 6: Update noise accuracy The value of is used to obtain the expression . ,in, Indicates that the array receives data The List, Let be the hyperparameters of the gamma distribution. This represents the trace operation of the matrix; Step 7: Calculate the position error of the array elements. The expression for the position error of the m-th array element is: ,in, ， ， ， ， ， ， ， ， Indicates the first One position is 1, and the rest are 0. Indicates the operation of taking the real part; Step 8: Calculate the accuracy of the array element position error. The expression for calculating its m-th element is: ， ,in, It is a hyperparameter; Step 9: Calculate the parameters ,in Indicates the number of iterations. It is the mean matrix calculated in step 4, if the parameter Meets error accuracy or If the maximum number of iterations (maxIter) is satisfied, proceed to step 10; if neither condition is satisfied, then... ， Re-enter step 4 for another iteration; Step 10: Update mesh error ,use Estimate the grid error vector. ,in, ， ， ， ， ， Indicates grid spacing. ， Represent the Hadamard product; Step 11: Calculate the spatial spectrum ， This represents the row mean vector, and the superscript "*" indicates the vector conjugate operation; Step 12: Update the spatial grid points using the grid error calculated in Step 10, i.e. At the same time, it corresponds one-to-one with the spatial spectrum in step 11, and the angle corresponding to the peak of the spatial spectrum is the estimated direction of arrival of the K signals.

Images