A robust sparse bayesian two-dimensional bearing of arrival estimation method

Abstract

The application discloses a kind of robust sparse bayesian two-dimensional boresight angle estimation method.It is characterized in that: the patent method is by introducing a new auxiliary angle, two-dimensional DOA estimation problem is converted into two one-dimensional DOA estimation problems.Through solving two sparse reconstruction problems containing amplitude and phase errors, the auxiliary angle and its corresponding pitch angle are estimated respectively, and according to the relationship between the three angles, the azimuth angle estimation and angle automatic matching process are completed simultaneously.For the sparse reconstruction problem containing amplitude and phase errors, the expectation maximum algorithm is used to derive the estimation expression of all unknown parameters, obtain a new spatial spectrum function representation method, and finally obtain the angle information.The patent method can effectively complete the angle automatic matching process, and improve the azimuth estimation accuracy and angle resolution capability, so it has high application value in practical engineering.

Description

Technical Field

[0001] This invention belongs to the field of sensor array signal processing technology in the field of signal and information processing. It relates to a robust sparse Bayesian two-dimensional azimuth-arrival estimation method. Specifically, it is a method that can robustly achieve two-dimensional azimuth-arrival estimation based on the sparse Bayesian method when there are amplitude errors and phase errors in the sensor array elements. Background Technology

[0002] Direction-of-arrival (DOA) estimation for sensor arrays is an important research area in radar, sonar, mobile communications, and other fields, attracting considerable interest from scholars. Currently, research on uniform liner arrays (ULAs) is relatively mature, and many high-resolution DOA estimation methods have been proposed to improve accuracy and resolution, such as subspace methods, maximum likelihood methods, and sparse Bayesian learning methods. Among these methods, by exploring the spatial sparsity of the incident signal, sparse Bayesian learning methods exhibit better adaptive performance than other traditional methods under conditions of limited snapshot numbers, low signal-to-noise ratio (SNR), and spatially adjacent signals. However, for ULAs, it can only provide one-dimensional angular information. In practical applications, finding the two-dimensional direction of the signal source, i.e., azimuth and elevation angles, is more reasonable. In recent years, research on two-dimensional DOA estimation has gradually increased, and scholars have proposed many array geometries for two-dimensional estimation, such as circular arrays, rectangular arrays, and cross arrays. Compared to these two-dimensional arrays, L-shaped arrays have a simpler structure, are easier to implement, and have better estimated performance, which has attracted widespread attention.

[0003] While the aforementioned azimuth estimation methods all exhibit good DOA estimation performance, they rely on an ideal array manifold. However, in practical applications, the array manifold is often affected by unknown amplitude and phase errors. Without array manifold calibration, the performance of azimuth estimation can be significantly degraded. Therefore, researching DOA estimation methods in the presence of array amplitude and phase errors has both theoretical and practical value. Existing amplitude and phase error calibration methods are broadly classified into two categories: active calibration and self-calibration. Generally, active calibration methods estimate amplitude and phase errors by placing a calibration source in space with accurately known arrival directions and then correcting the received data during normal operation. However, the presence of a calibration source is difficult to guarantee in practical applications, making implementation challenging. Compared to the first method, self-calibration methods can directly estimate the amplitude and phase errors of the array during operation without placing a calibration source. Although these methods typically employ iterative methods, resulting in high computational costs, their practical engineering application value is greater.

[0004] To address the issues of low azimuth estimation accuracy and poor angle resolution caused by unknown amplitude and phase errors in arrays, this patent proposes a robust sparse Bayesian two-dimensional DOA estimation method for L-shaped arrays with amplitude and phase errors. Compared to traditional DOA estimation methods, this method improves azimuth estimation accuracy and angle resolution, and can effectively perform automatic angle matching, 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 a robust sparse Bayesian two-dimensional DOA estimation method. Its key feature is that this patented method transforms the two-dimensional DOA estimation problem into two one-dimensional DOA estimation problems by introducing a new auxiliary angle. By solving two sparse reconstruction problems containing amplitude and phase errors, the auxiliary angle and its corresponding elevation angle are estimated respectively. Then, based on the relationship between the elevation and auxiliary angles, the azimuth estimation and automatic angle matching processes are simultaneously completed, reducing the impact of amplitude and phase errors on the two sparse reconstructions and achieving higher-precision source azimuth estimation.

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

[0007] Step 1: Use a uniform L-shaped array with 2M-1 array elements. The subarrays on both the y-axis and z-axis are uniform linear arrays containing M elements, with a sampling length of T. The received data of the m-th element of each subarray at time t is represented as y... m (t) and z m (t), where t = 1, 2, ..., T, m = 1, 2, ..., M;

[0008] Step 2: Arrange the received data of the M subarray elements on the y-axis and z-axis into vector form y(t) = [y1(t), ..., y2(t)]. M z(t)],z(t)=[z1(t),…,z M (t)], according to the formula Construct the sample covariance matrix, and then take the diagonal elements of its cross-correlation matrix to obtain r = diag(R). zy ), where the superscript "H" indicates the conjugate transpose operation on the matrix, and diag(·) indicates taking the diagonal elements;

[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: yes The first derivative, In the above formula, d represents the spacing between array elements, and λ represents the wavelength of the incident signal;

[0010] Step 4: Initialize parameters: b = d = f = 0.01, a = c = e = b + 1, δ( 0 )=1,α s (0) =0, α γ (0) =0, α0 (0) =mean(var(r)), γ (0) =1,β (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 ;

[0011] Step 5: Based on the sparse Bayesian inference, obtain the mean μ of the sparse signal vector δ. t The covariance matrix Σ s And the update formulas for other parameter vectors: Σ s =[α0Ω H (β,γ)Ω(β,γ)+diag(α s )] -1 ,in Represents α s The nth value, where Σ s,n,n Represents Σ s The element in the nth row and nth column, μ n,t μ t The nth value; γ=H -1 d, in I M Represents an M×M dimensional identity matrix; α γ The update formula for the m-th value is β=P -1 v, where In the formula Indicates taking the real part;

[0012] Step 6: 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, via and α s Perform a spectral peak search to obtain an estimate of the auxiliary angle η; otherwise, return to step five and continue the iterative process.

[0013] Step 7: Reassemble the received data on the y-axis and z-axis into a single representation. Representing sparse dictionary sets in This represents the auxiliary angle η obtained from the solution. k The estimated value,

[0014] Step 8: Initialize parameters, b = d = f = 0.01, a = c = e = b + 1, S( 0 )=1,α s (0) =0, α γ (0) =0, α0 (0) =mean(var(r)), γ (0) =1,β (0) =0, as well as To set a coarse mesh, set the values ​​of the iteration termination condition τ and Iter. max ;

[0015] Step 9: Based on the sparse Bayesian inference, obtain the mean of the sparse signal vector S. Covariance Matrix And the update formulas for other parameter vectors: in Represents α s The nth value, where express The element in the nth row and nth column, express The nth value; in I 2M Represents a 2M×2M dimensional identity matrix; α γ The update formula for the m-th value is β=P -1 v, where

[0016] Step 10: 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 the iteration and obtain the estimated α. s Dividing the spectrum into K blocks, the estimated pitch angle θ corresponding to each auxiliary angle is obtained by performing a spectral peak search on each block. Otherwise, return to step nine and continue the loop iteration;

[0017] Step 11: Based on the relationship between the triangles The estimated value of the azimuth angle corresponding to each pitch angle is obtained to achieve automatic matching.

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

[0019] (1) The method of this invention addresses the sparse reconstruction problem containing amplitude and phase errors. It employs the expectation-maximization (EM) algorithm to derive estimation expressions for all unknown parameters, resulting in a new spatial spectral function representation. Then, the angle information is finally obtained through spectral peak search. This process reduces the impact of amplitude and phase errors on estimation performance, not only improving the accuracy of azimuth estimation but also significantly enhancing the resolution success rate.

[0020] (2) This invention introduces a new auxiliary angle, transforming the two-dimensional DOA estimation problem into two one-dimensional DOA estimation problems. By solving two sparse reconstruction problems containing amplitude and phase errors, the auxiliary angle and its corresponding pitch angle are estimated respectively. Then, based on the relationship between the pitch angle and the auxiliary angle, the corresponding azimuth angle is solved, realizing automatic angle matching, which has higher practical engineering application value. Attached Figure Description

[0021] Figure 1 This is a diagram of the L-shaped array model of the signal processing method of this patent.

[0022] Figure 2 This is a spatial spectrum estimation diagram of the signal processing method of this patent under different signal-to-noise ratios;

[0023] Figure 3 Spatial spectrum estimation diagram of the signal processing method of this patent compared with other methods;

[0024] Figure 4 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.

[0025] 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 is shown.

[0026] Figure 6This is the relationship curve between the root mean square error and the standard deviation coefficient of amplitude and phase error in the signal processing method of this patent.

[0027] Figure 7 The curve showing the relationship between the successful resolution probability and the signal-to-noise ratio of the signal processing method in this patent;

[0028] Figure 8 The curve showing the relationship between the successful resolution probability and the standard deviation coefficient of amplitude and phase error in the signal processing method of this patent. Detailed Implementation

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

[0030] Example 1: Figure 1 A model diagram of the L-shaped array in this invention is given. First, a uniform L-shaped array with 2M-1 = 15 array elements is constructed. The subarrays on the y-axis and z-axis are both uniform linear arrays containing M = 8 array elements. The array spacing d is half a wavelength of 0.03 m. Three incoherent signals are incident on the receiving array from directions (θ1,φ1) = (14.2°, 40.3°), (θ2,φ2) = (30.5°, 15.6°), and (θ3,φ3) = (60.8°, 80.7°), respectively. The signal-to-noise ratios are set to 0 dB and 15 dB, and the number of snapshots is set to T = 100. Under these conditions, the specific implementation process is as follows:

[0031] Step 1: Use a uniform L-shaped array with 2M-1=15 array elements. The subarrays on the y-axis and z-axis are both uniform linear arrays containing M=8 array elements, with a sampling length of T=100. The received data of the m-th element of each subarray at time t is represented as y... m (t) and z m (t), where t = 1, 2, ..., 100, m = 1, 2, ..., 8;

[0032] Step 2: Arrange the received data of the M=8 subarray elements on the y-axis and z-axis into vector forms y(t)=[y1(t),…,y8(t)] and z(t)=[z1(t),…,z8(t)], respectively, according to the formula... Construct the sample covariance matrix, and then take the diagonal elements of its cross-correlation matrix to obtain r = diag(R). zy ), where the superscript "H" indicates the conjugate transpose operation on the matrix, and diag(·) indicates taking the diagonal elements;

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

[0034] Step 4: Initialize parameters: b = d = f = 0.01, a = c = e = b + 1, δ( 0 )=1,α s (0) =0, α γ (0) =0, α0 (0) =mean(var(r)), γ (0) =1,β (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 =300;

[0035] Step 5: Based on the sparse Bayesian inference, obtain the mean μ of the sparse signal vector δ. t The covariance matrix Σ s And the update formulas for other parameter vectors: Σ s =[α0Ω H (β,γ)Ω(β,γ)+diag(α s )] -1 ,in Represents α s The nth value, where Σ s,n,n Represents Σ s The element in the nth row and nth column, μ n,t μ t The nth value; γ=H -1 d, in I8 represents an 8×8 identity matrix; α γ The update formula for the m-th value is β=P -1 v, where In the formula Indicates taking the real part;

[0036] Step 6: 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 =300, stop iteration, pass and α s Perform a spectral peak search to obtain an estimate of the auxiliary angle η; otherwise, return to step five and continue the iterative process.

[0037] Step 7: Reassemble the received data on the y-axis and z-axis into a single representation. Representing sparse dictionary sets in This represents the auxiliary angle η obtained from the solution. k The estimated value,

[0038] Step 8: Initialize parameters: b = d = f = 0.01, a = c = e = b + 1, S (0) =1, α s (0) =0, α γ (0) =0, α0 (0) =mean(var(r)), γ (0) =1,β (0) =0, as well as To set a coarse mesh, the iteration termination condition is set to τ = 10. -3 and Iter max =300;

[0039] Step 9: Based on the sparse Bayesian inference, obtain the mean of the sparse signal vector S. Covariance Matrix And the update formulas for other parameter vectors: in Represents α s The nth value, where express The element in the nth row and nth column, express The nth value; γ=H -1 d, in I 2×8Represents a 16×16 dimensional identity matrix; α γ The update formula for the m-th value is β=P -1 v, where

[0040] Step 10: 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 =300, stop iteration, and calculate the estimated α. s Divided into K=3 blocks, the estimated pitch angle θ corresponding to each auxiliary angle is obtained by searching for spectral peaks in each block. Otherwise, return to step nine and continue the loop iteration;

[0041] Step 11: Based on the relationship between the triangles The estimated value of the azimuth angle corresponding to each pitch angle is obtained to achieve automatic matching.

[0042] The estimation results of this patented method after 150 Monte Carlo experiments under signal-to-noise ratio (SNR) conditions of 0 dB and 15 dB are as follows: Figure 2 As shown; by Figure 2 It can be seen that when the signal-to-noise ratio (SNR) is 15 dB, the DOA estimation result of this patented method is close to the true direction of the incident wave; when the SNR decreases to 0 dB, although the estimation accuracy decreases, it can still effectively estimate the approximate location of the source; when the SNR is 10 dB, the spatial spectrum estimated by this patented method, the 2D-MUSIC method, the OMP method, and the OGSBI method is as follows: Figure 3 As shown, through Figure 3 It can be seen that, compared with the other three algorithms, the method of this patent reduces the impact of amplitude and phase errors on the position estimation and is closest to the true source position. The OMP algorithm has the worst estimation performance.

[0043] The second embodiment: This study investigates the relationship between the estimation accuracy of the signal processing method of this patent and the signal-to-noise ratio, snapshot number, and standard deviation coefficient of amplitude and phase error. The results of the root mean square error changing with the signal-to-noise ratio are shown in the following figure. Figure 4 As shown; the result of the root mean square error changing with the number of snapshots is displayed, as follows. Figure 5 As shown; and displaying the results of the root mean square error varying with the standard deviation coefficient of amplitude and phase error, as shown. Figure 6 As shown. The conditions for applying the algorithm in this invention are as follows:

[0044] We used a uniform L-shaped array with 2M-1=15 elements as an example, with an array spacing of 0.03m (half a wavelength). Three incoherent signals were incident on the receiving array from directions (θ1,φ1)=(14.2°,40.3°), (θ2,φ2)=(30.5°,15.6°), and (θ3,φ3)=(60.8°,80.7°). We initially set the snapshot number T=100 and the amplitude and phase error standard deviation coefficient to 0.5. We then changed the signal-to-noise ratio starting from -5dB. The step size was increased from 2dB to 15dB; then, keeping the signal-to-noise ratio at 10dB and the standard deviation coefficient of amplitude and phase error constant at 0.5, the number of snapshots was changed from 20 to 200 in step sizes of 20; finally, the number of snapshots T=100 was set, the signal-to-noise ratio remained constant at 10dB, and the standard deviation coefficient of amplitude and phase error was changed from 5 to 0.5 in step sizes of 0.5. 300 independent Monte Carlo experiments were performed, and simulations were conducted using MATLAB software. The simulation results were observed. The simulation results are as follows: Figure 4 , Figure 5 as well as Figure 6 As shown.

[0045] from Figure 4 As can be seen, the root mean square error (RMSE) curves of all four methods gradually decrease with increasing signal-to-noise ratio (SNR). The 2D-MUSIC and OMP methods have higher RMSE curves and larger errors. The method proposed in this patent is superior to the other methods within this SNR range, exhibiting the lowest RMSE curve and the smallest error. At an SNR of 15 dB, the RMSE of the proposed method is 0.4188, which is a reduction of 0.9531 compared to the OGSBI method, 2.0574 compared to the 2D-MUSIC method, and 2.1952 compared to the OMP method. This indicates that the proposed method has better estimation performance under the condition of amplitude and phase errors, especially at high SNRs, where its advantages are more pronounced.

[0046] Depend on Figure 5 It can be seen that the root mean square error (RMSE) of all four methods decreases with increasing snapshot count. The RMSE curve of this patented method is lower than that of the other three methods within this snapshot count range, and the RMSE is only 1.7596 when the snapshot count is 20. This indicates that this patented method has the best estimation performance compared to other methods under different snapshot count conditions.

[0047] Depend on Figure 6 It can be seen that the root mean square error (RMSE) of all four methods increases with the increase of amplitude and phase error. However, the method of this patent is less affected by the magnitude of amplitude and phase error, and its RMSE curve remains the lowest in this range. The 2D-MUSIC method is more significantly affected by amplitude and phase error, with the largest fluctuation in RMSE. This indicates that compared with other methods, the method of this patent can more effectively reduce the influence of amplitude and phase error and has the best estimation performance.

[0048] The third embodiment: This study investigates the relationship between the angle resolution capability of the signal processing method of this patent and the signal-to-noise ratio (SNR) and the standard deviation coefficient of the amplitude and phase error. It displays the results of the resolution success rate changing with the SNR and the standard deviation coefficient of the amplitude and phase error, as shown in the following figure. Figure 7 , 8 As shown. The conditions for applying the algorithm in this invention are as follows:

[0049] We used a uniform L-shaped array with 2M-1=15 elements as an example, with an array spacing of half a wavelength of 0.03m. Three incoherent signals were incident on the receiving array from directions (θ1,φ1)=(14.2°,40.3°), (θ2,φ2)=(30.5°,15.6°), and (θ3,φ3)=(60.8°,80.7°). First, we set the number of snapshots T=100 and the standard deviation coefficient of amplitude and phase error to 0.5. We changed the signal-to-noise ratio (SNR) starting from -5dB and increasing it to 15dB in 2dB steps. Then, we set the number of snapshots T=100, the SNR to 15dB, and changed the standard deviation coefficient of amplitude and phase error starting from 5 and decreasing it to 0.5 in 0.5 steps. We conducted 300 independent Monte Carlo experiments, simulated them using MATLAB software, and observed the simulation results. The simulation results are as follows: Figure 7 , Figure 8 As shown.

[0050] Depend on Figure 7 It can be seen that, regardless of whether the signal-to-noise ratio (SNR) is low or high, the success rate of the proposed method is higher than that of the other methods. The proposed method achieves a 100% success rate at an SNR of 9 dB, while the other three methods do not reach 100% at an SNR of 15 dB and exhibit lower success rates and poorer resolution. This indicates that the proposed method has stronger angle resolution capabilities under different SNR conditions, and its robustness is also superior to the other methods.

[0051] Depend on Figure 8 It is known that when the amplitude and phase errors are large, the success resolution of other methods is low, with OMP and 2D-MUSIC methods almost failing to resolve the issue. In contrast, the resolution probability of the method proposed in this patent can reach 66.7%. When the amplitude and phase errors are small, the success resolution probability of the method proposed in this paper can reach 100%, which is significantly better than the resolution capabilities of other methods.

[0052] 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. A robust sparse Bayesian two-dimensional wave direction-of-arrival estimation method, characterized in that: The two-dimensional direction-of-arrival estimation method includes the following steps: Step 1: Use an array element count of A uniform L-shaped array, its shaft and Subarrays on the axis all contain A uniform linear array with n elements, and a sampling length of n. ; the first of the two subarrays Individual elements The received data at each time point are represented as follows: and ,in , ; Step 2: Separately shaft and On the axis The subarray data received by each element is arranged in vector form. , According to the formula Construct the sample covariance matrix, and then take the diagonal elements of its cross-correlation matrix to obtain... The superscript " " indicates that the matrix is ​​transposed using the conjugate operation. This indicates retrieving the diagonal elements; 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: , yes The first derivative, This represents the diagonal first derivative matrix of the array direction vectors corresponding to each point on the grid. The above formula represents the diagonal matrix of the array direction vectors corresponding to each point on the grid. Indicates the spacing between array elements. Indicates the wavelength of the incident signal; Step 4: Initialize parameters. , , , , , , , , , , as well as Set the value of the iteration termination condition for the coarse grid. and ; Step 5: Based on the sparse Bayesian inference, obtain the sparse signal vectors respectively. mean Covariance Matrix And the update formulas for other parameter vectors: , ,in , ; express The There are values, among which express The OK List the elements, express The One value; ; , , ,in , express 3D identity matrix; The The update formula for each value is: ; ,in , In the formula The symbol indicates taking the real part, and the superscript "*" indicates taking the complex conjugate operation; Step 6: Determine if the termination condition of the iteration is met. If it is met... Or, the maximum number of iterations is reached. Stop iteration, via and Auxiliary angles are obtained by performing spectral peak search. If the estimated value is obtained, then return to step five and continue the iterative process. Step 7: shaft and The received data on the axis is reassembled and represented as ,in , Representing sparse dictionary sets ,in , , This represents the auxiliary angle obtained from the solution. The estimated value, , ; , , ; , ; Step 8: Initialize parameters. , , , , , , , , , , as well as Set the value of the iteration termination condition for the coarse grid. and ; Step 9: Based on the sparse Bayesian inference, obtain the sparse signal vectors respectively. mean Covariance Matrix And the update formulas for other parameter vectors: , ,in , ; express The There are values, among which express The OK List the elements, express The One value; ; , , ,in , express 3D identity matrix; The The update formula for each value is: ; ,in , ; Step 10: Determine if the termination condition of the iteration is met. If it is met... Or, the maximum number of iterations is reached. Stop iterating and obtain the estimated value. Divided into The pitch angle corresponding to each auxiliary angle is obtained by performing a spectral peak search on each block. The estimated value Otherwise, return to step eight and continue the loop iteration; Step 11: Based on the relationship between the triangles The estimated value of the azimuth angle corresponding to each pitch angle is obtained to achieve automatic matching.

Images