阅读说明：本技术 一种基于稀疏贝叶斯的低复杂度信号波达方向估计方法 (Low-complexity signal direction-of-arrival estimation method based on sparse Bayes ) 是由 王一凡 王芳芳 于 2020-12-09 设计创作，主要内容包括：本发明公开了一种基于稀疏贝叶斯的低复杂度信号波达方向估计方法,旨在解决现有技术中MMV模型在考虑时间相关问题时计算复杂度高、计算效率低的技术问题。其包括：对入射信号进行扩展,利用马尔科夫概率先验模型获得扩展后的入射信号的先验；根据扩展后的入射信号的先验获得GAMP算法的输入函数和输出函数,进行GAMP迭代,获得恢复的扩展后的入射信号,并根据恢复的扩展后的入射信号获得信号波达方向。本发明能够在考虑时间相关因素的前提下降低DOA估计的计算复杂度,提高计算效率,准确估计的信号波达方向。(The invention discloses a sparse Bayes-based low-complexity signal direction-of-arrival estimation method, and aims to solve the technical problems of high computation complexity and low computation efficiency of an MMV (multi-media video) model in the prior art when a time correlation problem is considered. It includes: expanding the incident signal, and acquiring the prior of the expanded incident signal by using a Markov probability prior model; and obtaining an input function and an output function of the GAMP algorithm according to the prior of the expanded incident signal, performing GAMP iteration to obtain a recovered expanded incident signal, and obtaining a signal arrival direction according to the recovered expanded incident signal. The DOA estimation method and the DOA estimation device can reduce the calculation complexity of DOA estimation under the premise of considering time related factors, improve the calculation efficiency and accurately estimate the signal arrival direction.)

1. A low-complexity signal direction-of-arrival estimation method based on sparse Bayes is characterized by comprising the following steps:

expanding the incident signal based on spatial domain meshing;

obtaining the prior of the expanded incident signal by using a Markov probability prior model;

obtaining an input function and an output function of the GAMP algorithm according to the expanded incident signal prior;

performing GAMP iteration on the snapshot of the expanded incident signal by using an input function and an output function of a GAMP algorithm to obtain a recovered expanded incident signal;

and obtaining the arrival direction of the signal according to the recovered expanded incident signal.

2. The sparse Bayesian-based low-complexity signal direction-of-arrival estimation method according to claim 1, wherein the specific steps of incident signal expansion are as follows:

setting an incident signal as x, wherein a receiving array consists of M uniform linear array antennas, and when the receiving array receives T snapshots of the incident signal, an array output model is obtained:

y＝A(θ)x+e

wherein y is a receiving signal of the receiving array, A (theta) is an array flow pattern matrix, and e is white Gaussian noise received by the receiving array;

carrying out grid division on the spatial domain of the receiving array by using an exhaustion method to obtain an ultra-complete angle setWherein the content of the first and second substances,denotes the nth overcomplete angle, H is the overcomplete set of anglesThe number of angles in (1), 2, …, H;

based onAnd A (theta) obtaining an expanded array flow pattern matrix:

wherein A (theta) _ spark represents the expanded array flow pattern matrix;

obtaining an ultra-complete array output model according to the expanded array flow pattern matrix:

wherein the content of the first and second substances,which represents the received signal after the spreading,representing the expanded incident signal.

3. The sparse bayes-based low-complexity signal direction-of-arrival estimation method according to claim 2, wherein the spatial domain of said receiving array is [ -90 °,90 ° ].

4. The sparse Bayesian-based low-complexity signal direction-of-arrival estimation method according to claim 2, wherein the expanded incident signal is obtained a priori by the following steps:

modeling the expanded incident signal by using a Markov probability prior model to obtain the relation between the nth expanded incident signal at the tth snapshot and the t-1 st snapshot:

wherein the content of the first and second substances,the nth spread incident signal representing the t-th snapshot,the nth spread incident signal representing the t-1 snapshot,to representAndthe relation between beta is a time correlation coefficient, and beta belongs to (-1,1) and gamma_nThe prior error of the nth expanded incident signal, T ═ 1,2, …, T;

obtaining the influence of the previous snapshot according to the relation between the previous and the next snapshots of the expanded incident signalAnd the influence of the latter snapshotWherein the content of the first and second substances,representing incident signalsThe prior probability of (a) being,represents the mean value of t-1 snapshot forward transmission of the nth expanded incident signal t snapshot,represents the variance of t-1 snapshot forward delivery at the time of the nth expanded incident signal t snapshot,representing incident signalsThe prior probability of (a) being,represents the mean value of backward transfer of the nth expanded incident signal t +1 during snapshot,the variance of the backward transmission of the nth expanded incident signal t +1 during the snapshot is represented;

according toAndobtaining an expanded incident signalA priori of (a):

5. the sparse Bayesian-based low-complexity signal direction-of-arrival estimation method according to claim 4, wherein the nth expanded incident signal of the tth snapshotThe expression of (a) is as follows:

wherein the content of the first and second substances,

6. the sparse Bayesian-based low-complexity signal direction-of-arrival estimation method of claim 4, wherein an expression of an input function of the GAMP algorithm is as follows:

wherein, g_sRepresenting an input function, q^(t)An approximation of the received signal of the t-th snapshot without noise effects is represented,denotes q^(t)Noise variance of y^(t)Received signal, σ, representing the t-th snapshot reception matrix²Denotes y^(t)The noise variance of (2);

the expression of the output function of the GAMP algorithm is as follows:

wherein, g_xThe output function is represented by a function of the output,an approximation of the nth spread incident signal representing the t-th snapshot,to representThe noise variance of (2).

7. The sparse Bayesian-based low-complexity signal direction-of-arrival estimation method according to claim 6, wherein the operation of obtaining the recovered expanded incident signal is as follows:

initializing expanded incident signal for tth snapshot

Input function and output function pair using GAMP algorithmGAMP iteration is carried out, and the mean value and the variance of the expanded incident signals are updated;

updating iteration parameters by using an EM algorithm in each iteration process, and performing iteration convergence judgment by using an iteration convergence condition;

obtaining a recovered expanded incident signal according to the mean and variance of the expanded incident signal satisfying the iterative convergence judgment;

wherein the iteration convergence condition comprises a GAMP iteration convergence condition and an EM iteration convergence condition.

8. The method according to claim 7, wherein the maximum iteration number of the GAMP algorithm is G, and the iteration process of each iteration is as follows:

obtaining the expanded incident signal of the t-th snapshot in the ith iteration according to the output of the (i-1) th iterationVariance of (2)Mean value X^i(t)Prior error gammaⁱAnd the noise variance (σ) of the received signal²)ⁱWherein i ∈ [1, G ]]；

Let S ═ A (theta) _ spark-²By usingAnd X^i(t)Calculating an approximation q of the received signal of the t-th snapshot without noise effect in the ith iteration^i(t)：

Wherein, g_s ^(i-1)(t)An input function representing the t-snapshot of the i-1 th iteration;

using gammaⁱ、(σ²)ⁱAnd q is^i(t)Updating input function g of t-snapshot of ith iteration_s ^i(t)；

Using updated input function g_s ^i(t)Respectively calculating the approximate values r of the expanded incident signals of the t-th snapshot^i(t)And its noise variance

Wherein A (theta) _ spark^TDenotes the transpose of a (θ) _ spark, in order to input the damping coefficient,s^(i-1)(t)value, S, representing the i-1 th iteration^TDenotes the transpose of S, (g)_s ^i(t)) ' representing an input function g_s ^i(t)A derivative of (a);

using gammaⁱ、r^i(t)And τ_r ^i(t)Updating output function g of t snapshot of ith iteration_x ^i(t)；

According to the updated output function g_x ^i(t)Updating the variance and mean of the expanded incident signal:

wherein the content of the first and second substances,is shown asVariance of the expanded incident signal output for i iterations, (g)_x ^i(t)) ' representing the output function g_x ^{i (t)}Derivative of (A), X^(i+1)(t)Represents the mean of the expanded incident signal output for the ith iteration,in order to output the damping coefficient,

iterative convergence Condition on X with GAMP^(i+1)(t)Andperforming GAMP iterative convergence judgment, finishing iteration when GAMP iterative convergence conditions are met, outputting the mean value and the variance of the expanded incident signals of the ith iteration, and entering the next step when the GAMP iterative convergence conditions are not met;

based on EM algorithm, utilizing X output by ith iteration^(i+1)(t)Andrespectively calculating prior error gamma in the (i + 1) th iterationⁱ⁺¹Sum noise variance (σ)²)ⁱ⁺¹The calculation formula is as follows:

(σ²)ⁱ⁺¹＝||y^(t)-A(θ)_sparse·X^(i+1)(t)||²

using EM iterative convergence conditions on X^(i+1)(t)Andperforming EM iterative convergence judgment, and outputting the expanded incident signal of the (i + 1) th iteration when the EM iterative convergence condition is satisfiedAnd (5) finishing iteration by means of the mean value and the variance, and continuing the iteration when the EM iteration convergence condition is not met.

9. The sparse bayes-based low-complexity signal direction of arrival estimation method according to claim 8, wherein the GAMP iteration convergence condition is:

or the iteration number i reaches the maximum iteration number G of the GAMP algorithm, wherein epsilon_gampRepresenting the normalized tolerance parameter of the GAMP algorithm iteration.

10. The sparse Bayesian-based low-complexity signal direction-of-arrival estimation method according to claim 8, wherein the EM iteration convergence condition is:

or the iteration number i reaches the maximum iteration number K of the EM algorithm, wherein epsilon_emNormalized tolerance parameters representing iterations of the EM algorithm.

Technical Field

The invention relates to a low-complexity signal direction of arrival estimation method based on sparse Bayes, belonging to the technical field of signal processing.

Background

The Direction of Arrival (DOA) of a signal is the Direction of incoming waves of incoming and outgoing signals estimated by using received data of an antenna array in a noisy environment, the basic principle is that the spatial angle of the signal arriving at the array is estimated by using the phase difference existing between the received data of different array elements of a spatial array, and the method can be applied to the fields of radio communication, radar, sonar, navigation, seismic detection, biomedicine and the like and has important significance.

The SBL algorithm which has emerged in recent years is originally proposed by Tipping et al in 2001 and later, and is introduced into the field of sparse Signal recovery/compressed sensing, SBL is initially applied to a model of a Single Measurement Vector (SMV), and is gradually expanded to Multiple Measurement Vectors (MMV), so that the SBL algorithm has the advantages of high resolution and low calculation amount, but the DOA model used in MMV is ideal, and does not consider more practical factors. Zhang Chilean et al in 2011 applied SBL to time-related direction finding scenes, proposed a time-related Sparse Bayesian Learning (TSBL) algorithm, introduced a hyper-parameter to control the time correlation between snapshots, so that although the time-related problem can be solved, the introduced hyper-parameter will lead to the increase of the calculated amount, so that the efficiency becomes very low, if on large-scale problems such as more array elements, accurate direction finding cannot be realized, and further algorithm optimization must be performed.

Disclosure of Invention

In order to solve the problems of high computation complexity and low computation efficiency of an MMV model in the prior art when the time correlation problem is considered, the invention provides a low-complexity signal direction-of-arrival estimation method based on sparse Bayes.

In order to solve the technical problems, the invention adopts the following technical means:

the invention provides a low-complexity signal direction of arrival estimation method based on sparse Bayes, which comprises the following steps:

expanding the incident signal based on spatial domain meshing;

obtaining the prior of the expanded incident signal by using a Markov probability prior model;

obtaining an input function and an output function of the GAMP algorithm according to the expanded incident signal prior;

performing GAMP iteration on the snapshot of the expanded incident signal by using an input function and an output function of a GAMP algorithm to obtain a recovered expanded incident signal;

and obtaining the arrival direction of the signal according to the recovered expanded incident signal.

Further, the specific steps of incident signal spreading are as follows:

setting an incident signal as x, wherein a receiving array consists of M uniform linear array antennas, and when the receiving array receives T snapshots of the incident signal, an array output model is obtained:

y＝A(θ)x+e (1)

wherein y is a receiving signal of the receiving array, A (theta) is an array flow pattern matrix, and e is white Gaussian noise received by the receiving array;

carrying out grid division on the spatial domain of the receiving array by using an exhaustion method to obtain an ultra-complete angle setWherein the content of the first and second substances,represents the nth overcomplete angle, and H is the overcomplete set of anglesThe number of angles in (1), 2, …, H;

based onAnd A (theta) obtaining an expanded array flow pattern matrix:

wherein A (theta) _ spark represents the expanded array flow pattern matrix;

obtaining an ultra-complete array output model according to the expanded array flow pattern matrix:

wherein the content of the first and second substances,which represents the received signal after the spreading,representing the expanded incident signal.

Further, the spatial domain of the receiving array is [ -90 °,90 ° ].

Further, the step of a priori acquiring the expanded incident signal is as follows:

modeling the expanded incident signal by using a Markov probability prior model to obtain the relation between the nth expanded incident signal at the tth snapshot and the t-1 st snapshot:

wherein the content of the first and second substances,the nth spread incident signal representing the t-th snapshot,the nth spread incident signal representing the t-1 snapshot,to representAndin relation of between, beta isCorrelation coefficient between beta and gamma, beta is from (-1,1)_nThe prior error of the nth expanded incident signal, T ═ 1,2, …, T;

obtaining the influence of the previous snapshot according to the relation between the previous and the next snapshots of the expanded incident signalAnd the influence of the latter snapshotWherein the content of the first and second substances,representing incident signalsThe prior probability of (a) being,represents the mean value of t-1 snapshot forward transfer at the time of the nth expanded incident signal t snapshot,represents the variance of t-1 snapshot forward delivery at the time of the nth expanded incident signal t snapshot,representing incident signalsThe prior probability of (a) being,represents the mean value of backward transfer of the nth expanded incident signal t +1 during snapshot,representing the variance of the backward transmission of the nth expanded incident signal t +1 during the snapshot;

according toAndobtaining an expanded incident signalA priori of (a):

further, the nth expanded incident signal of the tth snapshotThe expression of (a) is as follows:

wherein the content of the first and second substances,

further, the expression of the input function of the GAMP algorithm is as follows:

wherein, g_sRepresenting an input function, q^(t)An approximation of the received signal of the t-th snapshot without noise effects is represented,denotes q^(t)Noise variance of y^(t)Received signal, σ, representing the t-th snapshot reception matrix²Denotes y^(t)The noise variance of (2);

the expression of the output function of the GAMP algorithm is as follows:

wherein, g_xThe output function is represented by a function of the output,an approximation of the nth spread incident signal representing the t-th snapshot,to representThe noise variance of (2).

Further, the specific operation of obtaining the recovered expanded incident signal is as follows:

initializing expanded incident signal for tth snapshot

Input function and output function pair using GAMP algorithmGAMP iteration is carried out, and the mean value and the variance of the expanded incident signals are updated;

updating iteration parameters by using an EM algorithm in each iteration process, and performing iteration convergence judgment by using an iteration convergence condition;

obtaining a recovered expanded incident signal according to the mean and variance of the expanded incident signal satisfying the iterative convergence judgment;

wherein the iteration convergence condition comprises a GAMP iteration convergence condition and an EM iteration convergence condition.

Further, if the maximum iteration number of the GAMP algorithm is G, the iteration process of each iteration is as follows:

obtaining the output of the i-1 th iterationExpanded incident signal of tth snapshotVariance of (2)Mean value X^i(t)Prior error gammaⁱAnd the noise variance (σ) of the received signal²)ⁱWherein i ∈ [1, G ]]；

Let S ═ A (theta) _ spark-²By usingAnd X^i(t)Calculating an approximation q of the received signal of the t-th snapshot without noise effect in the ith iteration^i(t)：

Wherein, g_s ^(i-1)(t)An input function representing the t-snapshot of the i-1 th iteration;

using gammaⁱ、(σ²)ⁱAnd q is^i(t)Updating input function g of t-snapshot of ith iteration_s ^i(t)；

Using updated input function g_s ^i(t)Respectively calculating the approximate values r of the expanded incident signals of the t-th snapshot^i(t)And its noise variance

Wherein A (theta) _ spark^TDenotes the transpose of a (θ) _ spark, in order to input the damping coefficient,s^(i-1)(t)value, S, representing the i-1 th iteration^TDenotes the inversion of S, (g)_s ^i(t)) ' representing an input function g_s ^i(t)A derivative of (a);

using gammaⁱ、r^i(t)Andupdating output function g of t snapshot of ith iteration_x ^i(t)；

According to the updated output function g_x ^i(t)Updating the variance and mean of the expanded incident signal:

wherein the content of the first and second substances,represents the variance of the expanded incident signal output for the ith iteration, (g)_x ^i(t)) ' representing the output function g_x ^i(t)Derivative of (A), X^(i+1)(t)Represents the mean of the expanded incident signal output for the ith iteration,in order to output the damping coefficient,

iterative convergence Condition on X with GAMP^(i+1)(t)Andperforming GAMP iterative convergence judgment, finishing iteration when GAMP iterative convergence conditions are met, outputting the mean value and the variance of the expanded incident signals of the ith iteration, and entering the next step when the GAMP iterative convergence conditions are not met;

based on EM algorithm, utilizing X output by ith iteration^(i+1)(t)Andrespectively calculating prior error gamma in the (i + 1) th iterationⁱ⁺¹Sum noise variance (σ)²)ⁱ⁺¹The calculation formula is as follows:

(σ²)ⁱ⁺¹＝||y^(t)-A(θ)_sparse·X^(i+1)(t)||² (15)

using EM iterative convergence conditions on X^(i+1)(t)Andand performing EM iterative convergence judgment, outputting the mean value and the variance of the expanded incident signal of the i +1 th iteration when the EM iterative convergence condition is met, ending the iteration, and continuing the iteration when the EM iterative convergence condition is not met.

Further, the GAMP iteration convergence condition is as follows:

or the iteration number i reaches the maximum iteration number G of the GAMP algorithm, wherein epsilon_gampRepresenting the normalized tolerance parameter of the GAMP algorithm iteration.

Further, the EM iteration convergence condition is:

or the iteration number i reaches the maximum iteration number K of the EM algorithm, wherein epsilon_emNormalized tolerance parameters representing iterations of the EM algorithm.

The following advantages can be obtained by adopting the technical means:

the invention provides a sparse Bayesian-based low-complexity signal direction-of-arrival estimation method, which is characterized in that a Markov probability prior model is used for modeling, time-related factors during multiple beat time of direction-of-arrival estimation are taken into consideration, and the accuracy of a DOA model under a large-block beat time-related MMV model is improved. After the prior of the incident signal is obtained, the GAMP idea is combined, GAMP iteration is carried out by introducing intermediate parameters, and the step of matrix inversion in the process of solving the posterior probability density under the traditional Bayes framework is replaced by iterative operation, so that the calculation complexity is greatly reduced, the calculation efficiency is improved, and the resolution of the estimation of the signal arrival direction is improved.

In the signal direction of arrival estimation process, the method not only considers non-ideal environments such as time correlation and the like, but also balances good performances such as high resolution, low complexity and the like, breaks through the limitation of the existing SBL algorithm in the aspects of precision, complexity, robustness and the like, and has very important significance for the research of modern application scenes.

Drawings

FIG. 1 is a flow chart of the steps of a low-complexity signal direction of arrival estimation method based on sparse Bayesian of the present invention.

FIG. 2 is a diagram illustrating message passing among multiple snapshot time snapshots in an embodiment of the invention.

Fig. 3 is a normalized power spectrum of the TMSBL algorithm and the method of the present invention under the condition of low signal-to-noise ratio in the embodiment of the present invention.

Fig. 4 is a normalized power spectrum of the TMSBL algorithm and the method of the present invention under the condition of high signal-to-noise ratio in the embodiment of the present invention.

FIG. 5 is a comparison of the method of the present invention and TMSBL algorithm in RMSE in an embodiment of the present invention.

FIG. 6 is a comparison diagram of the method of the present invention and TMSBL algorithm in CUP-time according to the embodiment of the present invention.

Detailed Description

The technical scheme of the invention is further explained by combining the accompanying drawings as follows:

the invention provides a low-complexity signal direction of arrival estimation method based on sparse Bayes, which specifically comprises the following steps as shown in FIG. 1:

step 1, expanding an incident signal based on space domain grid division;

step 2, acquiring the prior of the expanded incident signal by using a Markov probability prior model;

step 3, obtaining an input function and an output function of the GAMP algorithm according to the expanded incident signal prior;

step 4, carrying out GAMP iteration on the snapshot of the expanded incident signal by using an input function and an output function of a GAMP algorithm to obtain a recovered expanded incident signal;

and 5, acquiring a signal arrival direction according to the recovered and expanded incident signal.

In the embodiment of the invention, W narrow-band far-field signal sources with the wavelength of lambda are arranged, and the angle theta is used_wThe angle of (a) is transmitted to a receiving end, W belongs to {1,2, ·, W }, the incident signal is x, and a receiving array of the receiving end consists of M uniform linear array antennas.

When the receiving array receives T snapshots of the incident signal, an array output model is obtained:

y＝A(θ)x+e (16)

wherein y is the received signal of the receiving array, A (theta) is the array flow pattern matrix, e is the Gaussian white noise received by the receiving array, e satisfies the normal distribution, e-N (0, sigma)²I)，σ²Representing the noise variance of the received signal.

The expression of the array flow pattern matrix is as follows:

based on the theory of sparse representation, the incoming direction of the incident signal is sparse in the whole space domain, and the distinguishable space domain of the receiving array is [ -90 degrees, 90 degrees DEG ]]Carrying out grid division on the spatial domain of the receiving array by using an exhaustion method to obtain an over-complete angle setThe incident signal is expanded, wherein,denotes the nth overcomplete angle, H is the overcomplete set of anglesThe number of angles in (1), 2, …, H.

Based onAnd A (theta) obtaining an expanded array flow pattern matrix A (theta) _ spark:

because the incoming signal falls on the grid, an ultra-complete array output model is obtained according to the expanded array flow pattern matrix:

wherein the content of the first and second substances,which represents the received signal after the spreading,which represents the incident signal after the spreading,

the invention recovers the expanded incident signal by using the received signal and A (theta) _ sparkReusing the expanded incident signalThe DOA direction of the incident signal x is obtained.

In step 2, the expanded incident signal is acquired a priori by the following steps:

step 201, under a Bayesian framework, considering time-related factors, namely, the incident signal is influenced by the previous and next snapshots at the current moment, a Markov probability prior model is used for modeling a signal source, and the product of the two messages of the previous and next snapshots is combined into a message.

The invention introduces a time correlation coefficient beta, wherein the beta belongs to (-1,1), and the beta can represent the relation between the previous snapshot and the next snapshot of a signal. Nth expanded incident signal of tth snapshotCan be expressed as:

wherein the content of the first and second substances,the nth spread incident signal representing the t-1 snapshot,γ_nthe nth expanded incident signal has an a priori error of T ═ 1,2, …, T.

The relationship between the nth expanded incident signal at the tth snapshot and the t-1 st snapshot is obtained according to equation (19):

wherein the content of the first and second substances,to representAndthe relationship between them.

Under MMV model, orderAs an experience of the nth expanded incident signal of the tth snapshot, the message transmission among the multiple snapshot is shown in fig. 2.

Assuming that the prior of forward message delivery is a Gaussian distribution, and parameters eta and psi are introduced as the mean and variance, the effect of the previous snapshot isIf the prior check of message transmission is also a Gaussian distribution, and parameters upsilon and phi are introduced as the mean value and the variance, the influence of the later snapshot isWherein the content of the first and second substances,representing incident signalsThe prior probability of (a) being,represents the mean value of t-1 snapshot forward transmission of the nth expanded incident signal t snapshot,represents the variance of t-1 snapshot forward delivery at the time of the nth expanded incident signal t snapshot,representing incident signalsThe prior probability of (a) being,represents the mean value of backward transfer of the nth expanded incident signal t +1 at snapshot time t,and the variance of the backward transmission of the nth expanded incident signal t +1 at the snapshot time t is shown.

MergingAndobtaining the nth expanded incident signal of the t fast shootingA priori of (a):

after the prior of the incident signal is obtained, the invention combines the thought of GAMP and carries out subsequent calculation through intermediate parameters to replace the relevant step of matrix inversion,and the computational complexity is reduced. The invention introduces an intermediate parameter r as an approximate value of an expanded incident signal, and gives a noise variance tau of r_rIntroducing an intermediate parameter q as an approximation of A (theta) x, i.e. of the received signal without noise influence, giving the noise variance tau of q_q。

Substituting the mean value of forward message transfer prior obtained by formula (20) according to the Gaussian probability density function and convolution operationSum varianceExpression (c):

therefore:

wherein the content of the first and second substances,an approximation of the nth spread incident signal representing the t-1 snapshot,to representThe noise variance of (2).

The mean value of backward message transfer prior can be obtained by the same methodSum varianceExpression (c):

the expression of the input function can be derived from the definition of the input function in the GAMP algorithm:

wherein, g_sRepresenting an input function, q^(t)An approximation of the received signal of the t-th snapshot without noise effects is represented,denotes q^(t)Z is a (θ) x, y^(t)Received signal, σ, representing the t-th snapshot reception matrix²Denotes y^(t)The noise variance of (2).

According to the definition and prior of output function in GAMP algorithmAn expression of the output function can be derived:

wherein, g_xThe output function is represented by a function of the output,after nth expansion for showing the t-th snapshotAn approximation of the incident signal is made,to representThe noise variance of (2).

In the embodiment of the present invention, the specific operation of step 4 is as follows:

initializing expanded incident signal for tth snapshot

Input function and output function pair using GAMP algorithmGAMP iteration is carried out, and the mean value and the variance of the expanded incident signals are updated;

updating iteration parameters by using an EM algorithm in each iteration process, and performing iteration convergence judgment by using an iteration convergence condition;

obtaining a recovered expanded incident signal according to the mean and variance of the expanded incident signal satisfying the iterative convergence judgment;

the iteration convergence conditions comprise GAMP iteration convergence conditions and EM iteration convergence conditions, convergence judgment is carried out twice in each iteration process according to the 2 iteration convergence conditions, and iteration is terminated when any convergence judgment is passed.

The maximum iteration number of the GAMP algorithm is set as G, when the first iteration is carried out, numerical values such as the iteration number, the mean value, the variance, the prior error of the expanded incident signal, the noise variance of the received signal and the like need to be initialized, and the numerical values can be changed along with the increase of the iteration number until the iteration convergence condition is met, and the recovered expanded incident signal is output.

Taking the ith iteration as an example, the iteration process is specifically as follows:

(1) obtaining the output of the i-1 th iterationExpanded incident signal of the tth snapshotVariance of (2)Mean value X^i(t)Prior error gammaⁱAnd the noise variance (σ) of the received signal²)ⁱWherein i ∈ [1, G ]]。

(2) Let S ═ A (theta) _ spark-²And obtaining a calculation formula of the variance of the approximation of the received signal:

wherein the content of the first and second substances,representing the variance of the approximation of the received signal of the t snapshot in the ith iteration without noise contribution.

(3) By usingAnd X^i(t)Calculating an approximation q of the received signal of the t-th snapshot without noise effect in the ith iteration^i(t)：

Wherein, g_s ^(i-1)(t)Represents the input function for the t-snapshot of the i-1 th iteration.

(4) Will gammaⁱ、(σ²)ⁱAnd q is^i(t)The input function g of the t snapshot of the ith iteration is updated by substituting the formula (28) in an equal way_s ^{i (t)}。

(5) For convergence of GAMP iteration, input damping coefficient is introducedAnd output damping coefficientIn addition, an intermediate parameter s is introduced, and the expression of s is as follows:

wherein s is^i(t)The value representing the ith iteration, which has no practical physical significance, is simply a mathematical equation, s^{(i -1)(t)}Representing the value of the (i-1) th iteration.

(6) Using updated input function g_s ^i(t)Respectively calculating the approximate values r of the expanded incident signals of the t-th snapshot^i(t)And its noise variance

Wherein A (theta) _ spark^TDenotes the transposition of A (theta) _ spark, S^TDenotes the transpose of S, (g)_s ^i(t)) ' represents an input function g_s ^i(t)The derivative of (c).

(7) Will gammaⁱ、r^i(t)And τ_r ^i(t)The output function g of the t snapshot of the ith iteration is updated by substituting the equation (29)_x ^{i (t)}。

(8) According to the updated output function g_x ^i(t)Updating the variance and mean of the expanded incident signal:

wherein the content of the first and second substances,represents the variance of the expanded incident signal output for the ith iteration, (g)_x ^i(t)) ' representing the output function g_x ^i(t)Derivative of (A), X^(i+1)(t)Representing the mean of the expanded incident signal output for the ith iteration.

(9) Iterative convergence Condition on X with GAMP^(i+1)(t)Andperforming GAMP iterative convergence judgment, wherein GAMP iterative convergence conditions are as follows:

or the iteration number i reaches the maximum iteration number G of the GAMP algorithm, wherein epsilon_gampRepresenting the normalized tolerance parameter of the GAMP algorithm iteration.

When the GAMP iteration convergence condition is satisfied, ending the iteration and outputtingAnd X^(i+1)(t)And when the GAMP iterative convergence condition is not met, the next step is carried out.

(10) Updating iteration parameters based on EM algorithm (expectation maximization algorithm), and outputting X by using ith iteration^(i+1)(t)Andrespectively calculating prior error gamma in the (i + 1) th iterationⁱ⁺¹Sum noise variance (σ)²)ⁱ⁺¹The calculation formula is as follows:

(σ²)ⁱ⁺¹＝||y^(t)-A(θ)_sparse·X^(i+1)(t)||² (37)

(11) using EM iterative convergence conditions on X^(i+1)(t)Andperforming EM iterative convergence judgment, wherein the EM iterative convergence conditions are as follows:

or the iteration number i reaches the maximum iteration number K of the EM algorithm, wherein epsilon_emNormalized tolerance parameters representing iterations of the EM algorithm.

When the EM iteration convergence condition is met, performing the (i + 1) th iteration, performing only the steps (1) - (8) in the (i + 1) th iteration, ending the iteration, and outputting the mean value X of the expanded incident signal of the (i + 1) th iteration^(i+2)(t)Sum varianceAnd when the EM iteration convergence condition is not met, continuing the iteration.

In step 5, each row of the recovered expanded incident signal contains only W non-zero values, the positions of which correspond to the DOA direction of the incident signal x.

Several sets of simulation experiments are provided to verify the effect of the method of the invention:

simulation experiment 1:

adopt even linear array to accept incident signal, establish array element number M of receiving array and become 12, be equipped with 4 signal sources, its incoming wave angle does respectively: -30.15, -10.23, 40.74, 28.39. In an MMV model, the sampling fast beat number T is 100, and the method and the prior TMSBL algorithm are respectively utilized to estimate the signal direction of arrival, the normalized power spectrograms of the method and the prior TMSBL algorithm under different noise environments are shown in figures 3 and 4, and it can be seen from the figures that under the condition that the SNR (signal to noise ratio) is 10dB, the TMSBL algorithm which also takes time related factors into consideration cannot estimate, but the method can still estimate the signal direction of arrival more accurately; under the condition of high signal-to-noise ratio, namely the SNR is 40dB, the peak of the power spectrum of the TMSBL algorithm is still inaccurate, and a false peak exists, so that the method is suitable for the DOA model estimation in the non-ideal environment.

Simulation experiment 2:

let the number of array elements M of the receiving array be 10, there are 2 signal sources, the incoming wave directions are-11.21 and 0.74, respectively, and the snapshot number T is 40. The predetermined grid point range is [ -90 °,90 ° ], with a resolution of 1 °. The method and the conventional TMSBL algorithm are utilized to research the performance comparison condition of the RMSE along with the change of the signal-to-noise ratio (SNR), and the root mean square error is obtained by averaging 500 Monte Carlo experiments aiming at each group of simulation parameters.

FIG. 5 is a comparison of the method of the present invention and TMSBL algorithm in RMES, and it can be seen that, when the signal-to-noise ratio is low and the time correlation factor is considered, the error of the method of the present invention is much smaller than that of the TMSBL algorithm, so the method of the present invention has higher accuracy in non-ideal environment.

Simulation experiment 3:

under the condition that other simulation conditions are guaranteed to be the same as those of the simulation experiment 2, the signal-to-noise ratio (SNR) is set to be 30dB, the performance comparison condition of CUP time changing along with the snapshot number is researched by utilizing the method and the TMSBL algorithm, the result is shown in figure 6, compared with the TMSBL, under the condition of large snapshot, the time required by the method is far shorter than the TMSBL, the calculation complexity of the method under the condition of large snapshot number is lower, and the calculation efficiency is higher.

In conclusion, the method solves the problem of signal source time correlation in MMV direction-of-arrival estimation under a Bayesian framework, utilizes GAMP to replace a matrix inversion part in the traditional SBL algorithm, reduces the computational complexity, and can estimate DOA more quickly and accurately. Under the non-ideal environment such as considering time correlation, the method has the advantages of high resolution, low complexity, high accuracy, better robustness and the like.

The above description is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, several modifications and variations can be made without departing from the technical principle of the present invention, and these modifications and variations should also be regarded as the protection scope of the present invention.