CN116718980A - Robust sparse Bayesian two-dimensional arrival azimuth estimation method - Google Patents

Abstract

The invention discloses a robust sparse Bayesian two-dimensional arrival azimuth estimation method. The method is characterized in that: the method converts the two-dimensional DOA estimation problem into two one-dimensional DOA estimation problems by introducing a new auxiliary angle. The auxiliary angles and the pitch angles corresponding to the auxiliary angles are estimated respectively by solving two sparse reconstruction problems containing amplitude-phase errors, and then the azimuth angle estimation and the automatic angle matching process are completed according to the relation among the three angles. Aiming at the sparse reconstruction problem containing amplitude-phase errors, an expected maximum algorithm is adopted to deduce estimated expressions of all unknown parameters, a new spatial spectrum function representation method is obtained, and finally angle information is obtained. The method can effectively complete the automatic angle matching process, and improves the azimuth estimation precision and the angle resolution capability, so that the method has higher application value in actual engineering.

Description

Robust sparse Bayesian two-dimensional arrival azimuth estimation method

Technical Field

The invention belongs to the technical field of sensor array signal processing in the field of signal and information processing, relates to a robust sparse Bayesian two-dimensional arrival azimuth estimation method, and in particular relates to a method capable of realizing two-dimensional arrival azimuth estimation based on a sparse Bayesian method in the case that amplitude errors and phase errors exist in sensor array elements.

Background

The direction-of-arrival (DOA) estimation of a sensor array is an important research direction in other fields such as radar, sonar, and mobile communication, and thus research of the direction has attracted many scholars' interests. Currently, research on uniform linear arrays (uniform liner array, ULA) is mature, and many high-resolution DOA estimation methods have been proposed to improve accuracy and resolution, such as subspace-like methods, maximum likelihood methods, sparse bayesian learning-like methods, and the like. Among these methods, the sparse bayesian learning method has better adaptive performance under the conditions of a limited snapshot number, a low signal-to-noise ratio (SNR) and a spatially adjacent signal than other conventional methods by exploring the spatial sparsity of the incident signal. However, for ULA, it can only provide one-dimensional angle information, and it is more reasonable to find the two-dimensional direction of the source, i.e. azimuth and elevation, in practical application. In recent years, two-dimensional DOA estimation research has been increasing, and many array geometries for two-dimensional estimation such as circular arrays, rectangular arrays, cross arrays, etc. have been proposed by students. Compared with these two-dimensional arrays, the L-shaped array is simpler in structure, easy to realize and better in estimation performance, and has attracted much attention.

While the several types of azimuth estimation methods described above all have good DOA estimation performance, they rely on an ideal array manifold. However, in practical application scenarios, the array manifold is often affected by unknown amplitude phase errors. Without array manifold calibration, the performance of the azimuth estimation may be greatly reduced. Therefore, the DOA estimation method under the condition of array amplitude-phase error is researched and has theoretical significance and practical value. Existing methods of amplitude and phase error calibration fall into two general categories, active correction and self-correction. Generally, the active correction method estimates the amplitude and phase error by placing a calibration source accurately knowing the arrival direction in space and corrects the received data in normal operation, but it is difficult to implement because it is difficult to secure the existence of the calibration source in practical use. In contrast to the first method, the self-correction method can directly estimate the amplitude and phase errors of the array when in operation without placing a calibration source. Although such methods generally employ iterative methods, the computational effort is large, but their practical engineering application value is greater.

In order to solve the problems of low azimuth estimation precision and poor angle resolution of the existing method caused by unknown amplitude and phase errors of the array, the patent provides a robust sparse Bayesian two-dimensional DOA estimation method under the condition of the amplitude and phase errors by combining an L-shaped array. Compared with the traditional DOA estimation method, the method improves the azimuth estimation precision and the angle resolution capability, can effectively complete the angle automatic matching process, and has better engineering application value

Disclosure of Invention

The invention aims at overcoming the defects of the prior art, and provides a robust sparse Bayesian two-dimensional arrival azimuth estimation method. The method is characterized in that: the patent method converts 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-phase errors, an auxiliary angle and a pitch angle corresponding to the auxiliary angle are estimated respectively, and then according to the relation between the pitch angle and the auxiliary angle, azimuth angle estimation and automatic angle matching processes are completed simultaneously, so that the influence of the amplitude-phase errors on two sparse reconstruction is reduced, and the information source azimuth estimation with higher precision is realized.

The DOA estimation method process of the invention comprises the following steps:

step one, adopting a uniform L-shaped array with the array element number of 2M-1, wherein subarrays on a y axis and a z axis are uniform linear arrays containing M array elements, and the sampling length is T; the received data of the m-th array element of the two subarrays at the time t are respectively expressed as y _m (t) and z _m (T), wherein t=1, 2, …, T, m=1, 2, …, M;

step two, respectively arranging subarray receiving data of M array elements on the y axis and the z axis into a vector form y (t) = [ y ] ₁ (t),…,y _M (t)]，z(t)＝[z ₁ (t),…,z _M (t)]According to the formulaConstructing a sample covariance matrix, and then taking diagonal elements of the cross correlation matrix to obtain r=diag (R _zy ) Wherein the superscript "H" indicates that the matrix is subjected to conjugate transpose operation, and diag (·) indicates that diagonal elements are taken;

step three, the space area is formed by minus 90 degrees and 90 degrees]Evenly dividing into N parts to obtain an angle gridSetting the off-grid overcomplete sparse dictionary set to Φ (β) =a+bdiag (β), β= [ β ] ₁ ,…,β _N ]Representing a grid error vector, array manifold matrix +.>Andwherein->Is a direction vector of the array, and can be written as Is->Is used as a first derivative of (a), in the above equation, d represents the array element pitch, and λ represents the wavelength of the incident signal;

initializing parameters, b=d=f=0.01, a=c=e=b+1, δ [ ] ⁰ )＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，Φ ⁽⁰⁾ (β)＝A ⁽⁰⁾ ＝A，B ⁽⁰⁾ =b, i=0 andfor the set coarse grid, the values τ and Iter of the iteration termination condition are set _max ；

Step five, respectively obtaining the average mu of the sparse signal vector delta according to the sparse Bayesian inference _t Sum covariance matrix Σ _s Update formulas for other parameter vectors:Σ _s ＝[α ₀ Ω ^H (β,γ)Ω(β,γ)+diag(α _s )] ^-1 wherein-> Representing alpha _s Is the nth value of (1), wherein Σ _s,n,n Representation of sigma _s N elements of row n, mu _n,t Representation mu _t N-th value of (c);γ＝H ^-1 d， wherein->I _M Representing an M x M dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein In->Representing the real part;

step six, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max Stop iteration, pass->And alpha _s Carrying out spectrum peak search to obtain an estimated value of the auxiliary angle eta; otherwise, returning to the step five again to continue the loop iteration;

step seven, representing the recombination of the received data on the y axis and the z axis asRe-representing sparse dictionary set->Wherein-> Represents the solved auxiliary angle eta _k Is used for the estimation of the (c),

step eight, initializing parameters, b=d=f=0.01, a=c=e=b+1, s # ⁰ )＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，Andfor the set coarse grid, the values τ and Iter of the iteration termination condition are set _max ；

Step nine, respectively obtaining the average value of the sparse signal vector S according to the sparse Bayesian inferenceSum covariance matrixUpdate formulas for other parameter vectors: /> Wherein-> Representing alpha _s N-th value of (2), wherein->Representation->N elements of row n of +.>Representation->N-th value of (c); wherein->I _2M Representing a 2M x 2M dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein

Step ten, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max Stopping iteration and obtaining the estimated alpha _s Dividing the system into K blocks, and obtaining an estimated value +_of a pitch angle theta corresponding to each auxiliary angle by carrying out spectral peak search on each block>Otherwise, returning to the step nine again to continue the loop iteration;

step eleven, according to the relation between trianglesAnd obtaining an estimated value of the azimuth angle corresponding to each pitch angle, and realizing automatic matching.

Compared with the prior art, the technical scheme provided by the invention has the following technical effects:

(1) Aiming at the sparse reconstruction problem containing amplitude-phase errors, the method adopts an Expected Maximum (EM) algorithm to deduce the estimated expression of all unknown parameters, obtains a new space spectrum function expression method, and finally obtains angle information through spectrum peak search. The process reduces the influence of the amplitude-phase error on the estimation performance, improves the azimuth estimation precision and obviously improves the resolution success rate.

(2) According to the invention, a new auxiliary angle is introduced, the two-dimensional DOA estimation problem is converted into two one-dimensional DOA estimation problems, the auxiliary angle and the corresponding pitch angle are respectively estimated by solving two sparse reconstruction problems containing amplitude-phase errors, and then the corresponding azimuth angle is solved according to the relation between the pitch angle and the auxiliary angle, so that the angle is automatically matched, and the method has higher practical engineering application value.

Drawings

FIG. 1 is a diagram of an L-type array model of the signal processing method of the present patent;

FIG. 2 is a graph of spatial spectrum estimation at different signal-to-noise ratios for the signal processing method of the present patent;

FIG. 3 is a spatial spectrum estimation diagram of the signal processing method of the present patent and other methods;

FIG. 4 is a graph showing the relation between the root mean square error and the signal to noise ratio in the signal processing method of the present patent;

FIG. 5 is a graph showing the relation between root mean square error and snapshot number in the signal processing method of the present patent;

FIG. 6 is a graph showing the relation between the root mean square error and the standard deviation coefficient of the amplitude and phase error in the signal processing method of the present patent;

FIG. 7 is a graph showing the relationship between the probability of successful resolution and the signal-to-noise ratio of the signal processing method of the present patent;

FIG. 8 is a graph showing the relationship between the probability of successful resolution and the standard deviation coefficient of amplitude and phase error in the signal processing method of the present patent;

Detailed Description

The invention will now be further described with reference to examples, figures:

example 1: FIG. 1 shows a model of an L-array in the present invention. Firstly, we construct a uniform L-shaped array with array elements of 2M-1=15, subarrays on the y axis and the z axis are uniform linear arrays with M=8 array elements, the array spacing d is half-wavelength 0.03M, and three incoherent signals are distributed from (theta ₁ ,φ ₁ )＝(14.2°,40.3°)，(θ ₂ ,φ ₂ )＝(30.5°,15.6°)，(θ ₃ ,φ ₃ ) The direction incident receiving array of = (60.8 °,80.7 °), the signal-to-noise ratio is set to 0dB and 15dB, respectively, and the snapshot count is set to t=100. By adopting the conditions, the specific implementation process is as follows:

step one, adopting a uniform L-shaped array with the array element number of 2M-1=15, wherein subarrays on the y axis and the z axis are uniform linear arrays with M=8 array elements, and the sampling length is T=100; the receiving data of the m-th array element of the two subarrays at the time t are respectivelyDenoted as y _m (t) and z _m (t), wherein t=1, 2, …,100, m=1, 2, …,8;

step two, respectively arranging subarray receiving data of M=8 array elements on the y axis and the z axis into a vector form y (t) = [ y ] ₁ (t),…,y ₈ (t)]，z(t)＝[z ₁ (t),…,z ₈ (t)]According to the formulaConstructing a sample covariance matrix, and then taking diagonal elements of the cross correlation matrix to obtain r=diag (R _zy ) Wherein the superscript "H" indicates that the matrix is subjected to conjugate transpose operation, and diag (·) indicates that diagonal elements are taken;

step three, the space area is formed by minus 90 degrees and 90 degrees]Evenly dividing into N=91 parts to obtain an angle gridSetting the off-grid overcomplete sparse dictionary set to Φ (β) =a+bdiag (β), β= [ β ] ₁ ,…,β ₉₁ ]Representing a grid error vector, array manifold matrix +.>Andwherein->Is a direction vector of the array, and can be written as Is->Is used as a first derivative of (a), in the above equation, d represents the array element pitch, and λ represents the wavelength of the incident signal;

initializing parameters, b=d=f=0.01, a=c=e=b+1, δ [ ] ⁰ )＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，Φ ⁽⁰⁾ (β)＝A ⁽⁰⁾ ＝A，B ⁽⁰⁾ =b, i=0 andfor the set coarse mesh, the value τ=10 of the iteration termination condition is set ^-3 And Iter _max ＝300；

Step five, respectively obtaining the average mu of the sparse signal vector delta according to the sparse Bayesian inference _t Sum covariance matrix Σ _s Update formulas for other parameter vectors:Σ _s ＝[α ₀ Ω ^H (β,γ)Ω(β,γ)+diag(α _s )] ^-1 wherein-> Representing alpha _s Is the nth value of (1), wherein Σ _s,n,n Representation of sigma _s N elements of row n, mu _n,t Representation mu _t N-th value of (c);γ＝H ^-1 d， wherein->I ₈ Representing an 8 x 8 dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein In->Representing the real part;

step six, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max Stop iteration by +.>And alpha _s Carrying out spectrum peak search to obtain an estimated value of the auxiliary angle eta; otherwise, returning to the step five again to continue the loop iteration;

step seven, representing the recombination of the received data on the y axis and the z axis asRe-representing sparse dictionary set->Wherein-> Represents the solved auxiliary angle eta _k Is used for the estimation of the (c),

step eight, initializing parameters, b=d=f=0.01, a=c=e=b+1, s ⁽⁰⁾ ＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，Andfor the set coarse mesh, the value τ=10 of the iteration termination condition is set ^-3 And Iter _max ＝300；

Step nine, respectively obtaining the average value of the sparse signal vector S according to the sparse Bayesian inferenceSum covariance matrixUpdate formulas for other parameter vectors: /> Wherein-> Representing alpha _s N-th value of (2), wherein->Representation->N elements of row n of +.>Representation->N-th value of (c);γ＝H ^-1 d， wherein->I _2×8 Representing a 16 x 16 dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein

Step ten, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max =300, stopping iteration, and estimating α _s Dividing the model into K=3 blocks, and obtaining an estimated value of a pitch angle theta corresponding to each auxiliary angle by carrying out spectrum peak search on each block>Otherwise, returning to the step nine again to continue the loop iteration;

step eleven, according to the relation between trianglesAnd obtaining an estimated value of the azimuth angle corresponding to each pitch angle, and realizing automatic matching.

The estimation result of 150 Monte Carlo experiments performed by the method under the conditions of 0dB signal to noise ratio and 15dB signal to noise ratio is shown in figure 2; as can be seen from fig. 2, the DOA estimation result of the method of the present patent is close to the real direction of the incident wave when the signal-to-noise ratio is 15dB; when the signal-to-noise ratio is reduced to 0dB, the approximate position of the source can be effectively estimated although the estimation accuracy is reduced; when the signal-to-noise ratio is 10dB, the spatial spectrums estimated by the method and the 2D-MUSIC method, the OMP method and the OGSBI method are shown in figure 3, and compared with other three algorithms, the method reduces the influence of the amplitude-phase error on the azimuth estimation, and is closest to the real information source position, and the estimation performance of the OMP algorithm is worst.

Example 2: researching the relation between the estimation precision and the signal-to-noise ratio, the snapshot number and the standard deviation coefficient of the amplitude-phase error of the signal processing method, and displaying the variation result of the root mean square error along with the signal-to-noise ratio, wherein the obtained effect diagram is shown in figure 4; the root mean square error is displayed as a function of the snapshot number, as shown in fig. 5; and displaying the variation of the root mean square error with the standard deviation coefficient of the amplitude phase error, as shown in fig. 6. The algorithm of the invention is applied under the following conditions:

by taking a uniform L-shaped array with array element number of 2M-1=15 as an example, the array spacing is 0.03M of half wavelength, and three incoherent signal slaves (theta ₁ ,φ ₁ )＝(14.2°,40.3°)，(θ ₂ ,φ ₂ )＝(30.5°,15.6°)，(θ ₃ ,φ ₃ ) The direction incidence receiving array with the angle of 60.8 degrees and 80.7 degrees is firstly provided with the snapshot number T=100, the standard deviation coefficient of the amplitude and phase error is 0.5, the signal to noise ratio is changed from-5 dB, and the step length is increased to 15dB by 2 dB; then keeping the signal-to-noise ratio at 10dB and the standard deviation coefficient of the amplitude and phase error at 0.5 unchanged, and increasing the number of changed snapshots from 20 to 200 by taking 20 as a step length; and finally setting the snapshot number T=100, changing the standard deviation coefficient of the amplitude-phase error from 5 with the signal-to-noise ratio of 10dB unchanged, respectively carrying out 300 independent Monte Carlo experiments with the step length of 0.5 reduced to 0.5, carrying out simulation by a MATLAB software system, and observing the simulation result. The simulation results are shown in fig. 4, 5 and 6.

As can be seen from fig. 4, the root mean square error curves of the four methods gradually decrease with the increase of the signal-to-noise ratio, the root mean square error curves of the 2D-MUSIC and OMP methods are higher, the error is larger, and the method is superior to other methods in the signal-to-noise ratio range, and the root mean square error curve is the lowest and the error is the smallest. At a signal-to-noise ratio of 15dB, the root mean square error of the method is 0.4188, the estimated error is reduced by 0.9531 relative to the OGSBI method, the estimated error is reduced by 2.0574 relative to the 2D-MUSIC method, and the estimated error is reduced by 2.1952 relative to the OMP method. This demonstrates that the method of the present patent has better estimation performance in the presence of amplitude-phase errors, and especially has more obvious advantages in the case of high signal-to-noise ratio.

As can be seen from FIG. 5, the root mean square error of the four methods all showed a decreasing trend with increasing snapshot count. The root mean square error curve of the method is lower than that of the other three methods in the snapshot interval, and the root mean square error is 1.7596 only when the snapshot number is 20. This illustrates that the method of the present patent performs best under different snapshot count conditions than other methods.

As can be seen from fig. 6, the root mean square error of the four methods all increases with the increase of the amplitude and phase error, but the method of the present invention has small influence on the amplitude and phase error, the root mean square error curve is kept to be the lowest in the interval, the 2D-MUSIC method is more obviously influenced by the amplitude and phase error, and the root mean square error floats the most. This shows that compared with other methods, the method of the patent can effectively reduce the influence of amplitude and phase errors, and has the best estimation performance.

Example 3: the relation between the angle resolution capability and the signal-to-noise ratio and the standard deviation coefficient of the amplitude-to-phase error of the signal processing method of the patent is studied, and the result that the resolution success rate changes along with the signal-to-noise ratio and the standard deviation coefficient of the amplitude-to-phase error changes is displayed, so that effect diagrams are shown in figures 7 and 8. The algorithm of the invention is applied under the following conditions:

by taking a uniform L-shaped array with array element number of 2M-1=15 as an example, the array spacing is 0.03M of half wavelength, and three incoherent signal slaves (theta ₁ ,φ ₁ )＝(14.2°,40.3°)，(θ ₂ ,φ ₂ )＝(30.5°,15.6°)，(θ ₃ ,φ ₃ ) The direction incidence receiving array with the angle of 60.8 degrees and 80.7 degrees is firstly provided with the snapshot number T=100, the standard deviation coefficient of the amplitude and phase error is 0.5, the signal to noise ratio is changed from-5 dB, and the step length is increased to 15dB by 2 dB; then setting the snapshot number T=100, setting the signal-to-noise ratio as 15dB, changing the standard deviation coefficient of the amplitude-phase error from 5, reducing the step length to 0.5 by 0.5, respectively performing 300 independent Monte Carlo experiments, performing simulation by a MATLAB software system, and observing the simulation result. The simulation results are shown in fig. 7 and 8.

As can be seen from fig. 7, the probability of success in resolution of the method of this patent is higher than that of the other methods, both in the case of low signal-to-noise ratio and in the case of high signal-to-noise ratio. The probability of success of resolution is 100% when the signal-to-noise ratio is 9dB, but the probability of success of resolution is lower when the signal-to-noise ratio is 15dB, and the resolution capability is poor. This shows that the method has stronger angle resolution capability under different signal-to-noise ratios, and the robustness is better than that of other methods.

As can be seen from fig. 8, under the condition of larger amplitude-phase error, the resolution of other methods is lower, so that OMP and 2D-MUSIC methods can hardly be resolved successfully, and the resolution probability of the method can reach 66.7%. For the case of smaller amplitude-phase errors, the successful resolution probability of the method can reach 100%, and compared with other methods, the resolution capability is obviously better.

The specific examples described in this patent are offered by way of illustration only. Those skilled in the art may make various modifications, additions or substitutions to the described embodiments without departing from the invention or the scope thereof as defined in the accompanying claims.

Claims (1)

1. A robust sparse Bayesian two-dimensional arrival azimuth estimation method is characterized in that: the two-dimensional direction of arrival estimation method comprises the following steps:

step one, adopting a uniform L-shaped array with the array element number of 2M-1, wherein subarrays on a y axis and a z axis are uniform linear arrays containing M array elements, and the sampling length is T; the received data of the m-th array element of the two subarrays at the time t are respectively expressed as y _m (t) and z _m (T), wherein t=1, 2, T, m=1, 2, M;

step two, respectively arranging subarray receiving data of M array elements on the y axis and the z axis into a vector form y (t) = [ y ] ₁ (t),,y _M (t)]，z(t)＝[z ₁ (t),,z _M (t)]According to the formulaConstructing a sample covariance matrix, and then taking diagonal elements of the cross correlation matrix to obtain r=diag (R _zy ) Wherein the superscript "H" indicates that the matrix is subjected to conjugate transpose operation, and diag (·) indicates that diagonal elements are taken;

step three, the space area is formed by minus 90 degrees and 90 degrees]Evenly dividing into N parts to obtain an angle gridSetting the off-grid overcomplete sparse dictionary set to Φ (β) =a+bdiag (β), β= [ β ] ₁ ,,β _N ]Representing a grid error vector, array manifold matrix +.>And->Wherein the method comprises the steps ofIs a direction vector of the array, and can be written asIs->Is used as a first derivative of (a), in the above equation, d represents the array element pitch, and λ represents the wavelength of the incident signal;

step four, initializing parameters, b=d=f=0.01, a=c=e=b+1, δ ⁽⁰⁾ ＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，Φ ⁽⁰⁾ (β)＝A ⁽⁰⁾ ＝A，B ⁽⁰⁾ =b, i=0 andfor the set coarse grid, the values τ and Iter of the iteration termination condition are set _max ；

Step five, obtaining sparsity according to the sparse Bayesian inferenceMean value mu of signal vector delta _t Sum covariance matrix Σ _s Update formulas for other parameter vectors:Σ _s ＝[α ₀ Ω ^H (β,γ)Ω(β,γ)+diag(α _s )] ^-1 wherein-> Representing alpha _s Is the nth value of (1), wherein Σ _s,n,n Representation of sigma _s N elements of row n, mu _n,t Representation mu _t N-th value of (c);γ＝H ^-1 d， wherein->I _M Representing an M x M dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein In->Representing the real part;

step six, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max Stop iteration, pass->And alpha _s Carrying out spectrum peak search to obtain an estimated value of the auxiliary angle eta; otherwise, returning to the step five again to continue the loop iteration;

step seven, representing the recombination of the received data on the y axis and the z axis asRe-representing sparse dictionary setsWherein-> Represents the solved auxiliary angle eta _k Is used for the estimation of the (c),

step eight, initializing parameters, b=d=f=0.01, a=c=e=b+1, s ⁽⁰⁾ ＝1，α _s ⁽⁰⁾ ＝0，α _γ ⁽⁰⁾ ＝0，α ₀ ⁽⁰⁾ ＝mean(var(r))，γ ⁽⁰⁾ ＝1，β ⁽⁰⁾ ＝0，i=0 and +.>For the set coarse grid, the values τ and Iter of the iteration termination condition are set _max ；

Step nine, respectively obtaining the average value of the sparse signal vector S according to the sparse Bayesian inferenceAnd covariance matrix->Update formulas for other parameter vectors: /> Wherein-> Representing alpha _s Of (c), whereinRepresentation->N elements of row n of +.>Representation->N-th value of (c); wherein->I _2M Representing a 2M x 2M dimensional identity matrix; alpha _γ Is +.about.the update of the mth value of (a)>β＝P ^-1 v, wherein

Step ten, judging whether the iteration termination condition is met, if soOr the maximum iteration number is reached, i.e. i is not less than Iter _max Stopping iteration and obtaining the estimated alpha _s Dividing the system into K blocks, and obtaining an estimated value +_of a pitch angle theta corresponding to each auxiliary angle by carrying out spectral peak search on each block>Otherwise, returning to the step eight again to continue the loop iteration;

step eleven, according to the relation between trianglesAnd obtaining an estimated value of the azimuth angle corresponding to each pitch angle, and realizing automatic matching.