1. A Bayesian robust beam forming method under a motion interference environment is characterized by comprising the following steps:

step 1: establishing an array received signal model;

step 2: establishing a Bayesian hierarchical probability model, wherein the specific contents are as follows:

describing the azimuthally varying characteristics of the interference + noise component by the Dirichlet process, i.e.

n_k～G

In the formula, the-representation obeys,is a discrete random variable and is supported in its support setThe sum of the probabilities of (c) is 1, δ (·) is a dirac function,π_f～Beta(1,γ)∝γ(1-π_f)^γ-1and oc represents proportional to, γ represents the scaling factor and obeys the gamma distribution as follows:

the f-th set of interference + noise components share a precision matrix, namely:

in the formula G₀Denotes the reference distribution, Λ_fExpressing the precision matrix corresponding to the f-th group of Gaussian distribution, and obeying the Wishart distribution

p(Λ_f)＝W(Λ_f|W,v)∞|Λ_f|^v-Netr{-W^-1Λ_f}

In the formula (I), the compound is shown in the specification,representing the mean matrix, v representing degrees of freedom;

introducing an allocation vector z_kWhich is a random variable extracted from a multi-valued distribution (the parameter set of the multi-valued distribution is { ω }_f}_f＝1,…,K) Namely:

in the formula z_kIs a K × 1 dimensional vector with only one element being 1;

definition 1[ z_k＝f]Denotes z_kThe f-th element in (1) is 1, i.e. n_kAnd if the class f is assigned, the likelihood function of the array observation data follows the following complex Gaussian distribution:

in the formulaIs thatA corresponding accuracy matrix;

the statistical properties of the desired signal amplitude are represented in a complex gaussian distribution:

wherein s ═ s₁,…,s_K]^TIs a vector of the amplitudes of the desired signal waveforms for the K snapshots,the desired signal accuracy;

desired signal accuracyObey the following gamma distribution

Wherein, gamma function is represented by gamma, a > 0 represents shape parameter, b > 0 represents expansion parameter;

the desired signal steering vector follows the following Watson distribution:

p(v)＝w(v|μ,λ)＝c_p(λ)exp(λ|μ^Hv|²)

in the formula (I), the compound is shown in the specification,is a regularization parameter, λ is a focusing parameter; c. C_pThe hypergeometric function in (λ) is defined as the following form:

the hyperparameters mu, lambda obey the following Watson-Gamma joint distribution

By collectionsAll unknown variables in the probability model are represented, and the statistical characteristics of the unknown variables are combined with the distribution of the observation data Y to obtain the following joint probability density:

and step 3: obtaining an updating formula of posterior distribution parameters of all hidden variables by using a variational inference method, wherein the concrete contents are as follows:

(1) initializing the probabilistic model parameters, i.e. let v ═ μ ═ m₀＝v_pWherein v is_pFor a preset desired signal steering vector, Λ_f＝I_N×N(f＝1,…,K)，π_f(f＝1,…,K)＝0.5，v＝N，a＝a₀＝b＝b₀＝c＝d＝10^-6，s＝1_K×1，W＝10⁺⁶·I_N×NIn which 1 is_K×1Is a K x 1 dimensional full 1 vector, I_N×NAn identity matrix of dimension NxN;

(2) f group interference + noise component Λ_fA posteriori probability q (Λ)_f) Subject to the Wishart distribution with its degrees of freedomSum-mean matrixThe update formulas of (a) and (b) are respectively:

in the formula (II)<q(z_k＝f)>＝φ_k,f，<·>Expressing the mathematical expectation operation; phi is a_k,fThe calculation formula of (2) is as follows:

in the formula (I), the compound is shown in the specification,

(3)π_fobey a beta distribution, i.e. q (pi)_f)＝Beta(π_f|c_f,d_f) Wherein the hyperparameter c_fAnd d_fThe update formulas of (a) and (b) are respectively:

in the formula (I), the compound is shown in the specification,

(4) the posterior distribution of the scaling factor gamma is the gamma distribution, i.e.Wherein the shape parameterAnd a parameter of expansion and contractionThe update formulas of (a) and (b) are respectively:

in the formula<ln(1-π_f)>＝Ψ(d_f)-Ψ(c_f+d_f) Ψ (. cndot.) is a digamma function;

(5) posterior probability q(s) of expected signal waveform amplitude in kth snapshot_k) Subject to complex Gaussian distributions, i.e.Mean value thereofAnd precisionThe iterative update formula of

(6) The posterior probability of the desired signal steering vector v obeys a Watson distribution, i.e.Wherein the hyperparameterAndare respectively a matrixAnd a principal eigenvector and a principal eigenvalue of, and

(7) the posterior distribution of μ is the Watson distribution, i.e.WhereinIs a positive definite matrix (beta mm)^H+vv^H) Is determined by the maximum characteristic value of the image,is its corresponding feature vector;

(8) the posterior distribution of λ is a Gamma distribution, i.e., q (λ) ═ Gam (λ | a)₁,b₁) Wherein the shape parameter a₁And a parameter of expansion and contraction b₁The update formulas of (a) and (b) are respectively:

in the above formula, the first and second carbon atoms are,calculated from the last iteration step;

(9) the posterior probability of the accuracy of the desired signal obeys the gamma distribution, i.e.Wherein the shape parameter a₂And a parameter of expansion and contraction b₂The update formulas of (a) and (b) are respectively:

a₂＝K+a

(10) iterating the processes (2) - (9) until a convergence condition is met;

and 4, step 4: converging and solving the f group interference + noise precision matrix and the expected signal guide vectorAnd v_optSubstitution formulaCalculating to obtain the optimal weight vector corresponding to the f-th array receiving dataAnd obtaining the optimal weight vectors corresponding to the received data of the other groups of arrays by adopting the same calculation method.

2. The bayesian robust beamforming method under the motion interference environment according to claim 1, wherein the array receiving signal model in step 1 is specifically as follows:

taking a linear array with an array element number N as an example, the received data at the kth sampling time can be represented as:

in the formulas_k，Each represents the array output at time k, the amplitude of the desired signal waveform and the interference + noise component, v represents the desired signal steering vector, and Y ═ Y₁,…,y_K]Representing the array received data set over the observation timeK is the number of sampled fast beats; under the line array signal model, the steering vector can be generally written as follows:

in the formula of₀Is the carrier wavelength, [ d₁,…，d_N]For physical coordinates of each array element, theta₀Is the spatial orientation of the desired signal.

3. The Bayesian robust beamforming method according to claim 1, wherein the convergence condition in step 3 isWhereinThe desired signal accuracy value calculated for the current step,the desired signal accuracy value calculated for the previous step.

4. A computer system, comprising: one or more processors, a computer readable storage medium, for storing one or more programs, wherein the one or more programs, when executed by the one or more processors, cause the one or more processors to implement the method of claim 1.

5. A computer-readable storage medium having stored thereon computer-executable instructions for, when executed, implementing the method of claim 1.

6. A computer program comprising computer executable instructions which when executed perform the method of claim 1.

Background

One effective measure for improving the interference suppression capability of the beamformer is to use an adaptive beamforming technique to adaptively generate a notch in the interference direction, thereby improving the Signal-to-interference-plus-noise ratio (SINR). However, since adaptive patterns typically have sharp notches, moving interferers tend to move out of the pattern "null," causing the performance of conventional robust adaptive beamforming algorithms to degrade dramatically in non-stationary interference environments. In a beam forming system for a motion interference source, because the orientation of an interference signal changes rapidly along with time, when the interference orientation is unknown a priori, an effective strategy for fusing beam output results on discrete time points to finally obtain efficient suppression of motion interference cannot be provided. Under the circumstance, how to effectively synthesize data at a plurality of discrete time points so as to realize beam forming aiming at a non-stationary interference environment and obtain stronger non-ideal signal environment adaptability than the existing method is a key problem to be solved by the invention.

The problem of motion interference suppression widely exists in the field of array beam forming, and can be regarded as an extension of the multi-model multi-observation problem in the field of array signal processing. The method for solving the problem of motion interference suppression in the multi-model multi-observation system is different from the static interference suppression method in the existing single-model system under the influence of model difference. J.R. Guerci (J.R. Guerci, Theory and application of fundamental matrix properties for robust adaptive beamforming, IEEE trans. Signal Process.47(4) (1999)977-985.) proposes a motion interference suppression algorithm based on the sampling covariance matrix tapering idea. The Covariance Matrix Taper (CMT) algorithm achieves the goal of enhancing beamformer robustness by broadening the notch at the interference azimuthIn (1). However, this algorithm only focuses on the problem of suppression of the motion interference source and is not robust to the desired signal steering vector error. The use of the various order derivatives of the received data as constraints for the beamforming algorithm, a.b.gershman et al (a.b.gershman, u.nickel, j.f).Adaptive beamforming algorithms with robustness against jammer motion,IEEE Trans.Signal Process.45(7)(1997)1878-1885.)(A.B.Gershman,E.Németh,J.F.The IEEE Trans. Signal Process.48(1) (2000)246 and 250) improves the traditional loaded sampling covariance matrix inversion algorithm and eigenvector projection algorithm and populates the algorithm to the motion interference environment. However, the two algorithms face the technical implementation difficulties that the diagonal loading factor is difficult to accurately calculate, and noise and signal subspace aliasing are difficult to realize under low SNR. S.A. Vorobyov et al (S.A. Vorobyov, A.B. Gershman, Z.Luo, N.Ma, Adaptive beamforming with joint robustness against interference mismatch vector and interference mismatch, IEEE Signal Process.Lett.11(2) (2004)108-111.) generalize the Worst Case algorithm to the non-stationary interference environment, however, the implementation of the algorithm faces difficulties such as that the covariance matrix of the array received data and the uncertain set of the desired Signal steering vector are unknown a priori. L.zhang et al (l.zhang, b.li, l.huang, t.kirubarajan, h.c.so, Robust dispersion less response beamforming against interference fast-moving interference, Signal process.140(2017)190-197.) propose a Robust beamforming algorithm based on the minimum deviation criterion, which achieves the purpose of eliminating interference energy from the beamformer output by constraining the average power in the interference dynamic azimuth set to be 0. The algorithm has the defects that the global convergence is difficult to guarantee and the algorithm is easy to disperse to a local optimal solution.

In view of the fact that the multistage probability structure of the Bayes method can eliminate coupling among different signal components, local features of the signal components can be well reserved, and accordingly high output SINR can be obtained. The method can eliminate the notch offset which may occur in the beam pattern of the traditional method by fully utilizing the time-varying characteristics of the interference azimuth, and simultaneously, the negative influence caused by the real-time variation of the interference signal array manifold is better avoided.

Disclosure of Invention

Technical problem to be solved

In order to avoid the defects of the prior art, the invention provides a Bayesian robust beam forming method under the motion interference environment.

Technical scheme

The unparameterized Bayes method is a Bayes combined reconstruction method for extracting interesting target information by comprehensively utilizing observation data obtained by different models, so that the method can be used for solving the problem of beam forming containing a plurality of observation models, namely motion interference suppression. The invention adopts a Dirichlet process to self-adaptively adjust a weight vector and track an interference azimuth according to an actual signal environment, and based on a local variation inference criterion, aims to solve the problem that the hidden variable posterior distribution caused by a non-conjugate distribution pair under a layered probability model is difficult to estimate.

The expected signal and the interference signal are assumed to be far-field narrow-band signals, which are incident on a uniform linear array composed of N array elements, and the spacing between the array elements is half of the wavelength of the incident signal. When the receiving signal of the linear array is used for self-adaptive beam forming, the technical scheme adopted by the invention comprises the following steps:

step 1: establishing an array received signal model;

step 2: establishing a Bayesian hierarchical probability model, wherein the specific contents are as follows:

describing the azimuthally varying characteristics of the interference + noise component by the Dirichlet process, i.e.

n_k～G

In the formula, the-representation obeys,is a discrete random variable and is supported in its support setThe sum of the probabilities of (c) is 1, δ (·) is a dirac function,π_f～Beta(1,γ)∝γ(1-π_f)^γ-1and oc represents proportional to, γ represents the scaling factor and obeys the gamma distribution as follows:

the f-th set of interference + noise components share a precision matrix, namely:

in the formula G₀Denotes the reference distribution, Λ_fExpressing the precision matrix corresponding to the f-th group of Gaussian distribution, and obeying the Wishart distribution

p(Λ_f)＝W(Λ_f|W，v)∝|Λ_f|^v-Netr{-W^-1Λ_f}

In the formula (I), the compound is shown in the specification,representing the mean matrix, v representing degrees of freedom;

introducing an allocation vector z_kWhich is a random variable extracted from a multi-valued distribution (the parameter set of the multi-valued distribution is { ω }_f}_f＝1,…,K) Namely:

in the formula z_kIs only oneA K × 1 dimensional vector with elements of 1;

definition 1[ z_k＝f]Denotes z_kThe f-th element in (1) is 1, i.e. n_kAnd if the class f is assigned, the likelihood function of the array observation data follows the following complex Gaussian distribution:

in the formulaIs thatA corresponding accuracy matrix;

the statistical properties of the desired signal amplitude are represented in a complex gaussian distribution:

wherein s ═ s₁,…,s_K]^TIs a vector of the amplitudes of the desired signal waveforms for the K snapshots,the desired signal accuracy;

desired signal accuracyObey the following gamma distribution

Wherein, gamma function is represented by gamma, a > 0 represents shape parameter, b > 0 represents expansion parameter;

the desired signal steering vector follows the following Watson distribution:

in the formula (I), the compound is shown in the specification,is a regularization parameter, λ is a focusing parameter; c. C_pThe hypergeometric function in (λ) is defined as the following form:

the hyperparameters mu, lambda obey the following Watson-Gamma joint distribution

By collectionsAll unknown variables in the probability model are represented, and the statistical characteristics of the unknown variables are combined with the distribution of the observation data Y to obtain the following joint probability density:

and step 3: obtaining an updating formula of posterior distribution parameters of all hidden variables by using a variational inference method, wherein the concrete contents are as follows:

(1) initializing the probabilistic model parameters, i.e. let v ═ μ ═ m₀＝v_pWherein v is_pFor a preset desired signal steering vector, Λ_f＝I_N×N(f＝1，…,K)，π_f(f＝1,…,K)＝0.5，v＝N，a＝a₀＝b＝b₀＝c＝d＝10^-6，s＝1_K×1，W＝10⁺⁶·I_N×NIn which 1 is_K×1Is a K x 1 dimensional full 1 vector, I_N×NAn identity matrix of dimension NxN;

(2) f group interference + noise component Λ_fA posteriori probability q (Λ)_f) Subject to the Wishart distribution with its degrees of freedomSum-mean matrixThe update formulas of (a) and (b) are respectively:

wherein, is defined as < q (z)_k×f)〉＝φ_k,fAnd (h) represents the mathematical expectation operation; phi is a_k,fThe calculation formula of (2) is as follows:

in the formula (I), the compound is shown in the specification,

(3)π_fobey a beta distribution, i.e. q (pi)_f)＝Beta(π_f|c_f,d_f) Wherein the hyperparameter c_fAnd d_fThe update formulas of (a) and (b) are respectively:

in the formula (I), the compound is shown in the specification,

(4) the posterior distribution of the scaling factor gamma is the gamma distribution, i.e.Wherein the shape parameterAnd a parameter of expansion and contractionThe update formulas of (a) and (b) are respectively:

in the formula < ln (1-pi)_f)〉＝Ψ(d_f)-Ψ(c_f+d_f) Ψ (. cndot.) is a digamma function;

(5) posterior probability q(s) of expected signal waveform amplitude in kth snapshot_k) Subject to complex Gaussian distributions, i.e.Mean value thereofAnd precisionThe iterative update formula of

(6) The posterior probability of the desired signal steering vector v obeys a Watson distribution, i.e.Wherein the hyperparameterAndare respectively a matrixAnd a principal eigenvector and a principal eigenvalue of, and

(7) the posterior distribution of μ is the Watson distribution, i.e.WhereinIs a positive definite matrix (beta mm)^H+vv^H) Is determined by the maximum characteristic value of the image,is its corresponding feature vector;

(8) the posterior distribution of λ is a Gamma distribution, i.e., q (λ) ═ Gam (λ | a)₁，b₁) Wherein the shape parameter a₁And a parameter of expansion and contraction b₁The update formulas of (a) and (b) are respectively:

in the above formula, the first and second carbon atoms are,calculated from the last iteration step;

(9) the posterior probability of the accuracy of the desired signal obeys the gamma distribution, i.e.Wherein the shape parameter a₂And a parameter of expansion and contraction b₂The update formulas of (a) and (b) are respectively:

a₂＝K+a

(10) iterating the processes (2) - (9) until a convergence condition is met;

and 4, step 4: converging and solving the f group interference + noise precision matrix and the expected signal guide vectorAnd v_optSubstitution formulaCalculating to obtain the optimal weight vector corresponding to the f-th array receiving dataAnd obtaining the optimal weight vectors corresponding to the received data of the other groups of arrays by adopting the same calculation method.

The further technical scheme of the invention is as follows: the array received signal model in the step 1 is specifically as follows:

taking a linear array with an array element number N as an example, the received data at the kth sampling time can be represented as:

in the formulaEach represents the array output at time k, the amplitude of the desired signal waveform and the interference + noise component, v represents the desired signal steering vector, and Y ═ Y₁,…,y_K]Representing the array received data set over the observation time, K being the number of sampled fast beats; under the line array signal model, the steering vector can be generally written as follows:

in the formula of₀Is the carrier wavelength, [ d₁，…,d_N]For physical coordinates of each array element, theta₀Is the spatial orientation of the desired signal.

The further technical scheme of the invention is as follows: the convergence condition in the step 3 isWhereinThe desired signal accuracy value calculated for the current step,the desired signal accuracy value calculated for the previous step.

A computer system, comprising: one or more processors, a computer readable storage medium, for storing one or more programs, which when executed by the one or more processors, cause the one or more processors to implement the above-described method.

A computer-readable storage medium having stored thereon computer-executable instructions for performing the above-described method when executed.

A computer program comprising computer executable instructions which when executed perform the method described above.

Advantageous effects

The existing theoretical achievements in the field of bayesian beamforming do not support an effective suppression of motion interference well yet. In the existing beam forming algorithm, a large distance exists between a static prior hypothesis and an array actual model about signal space domain distribution information contained in observation data, and the static prior hypothesis and the array actual model cannot be directly used for solving the beam forming problem of a plurality of observation models containing azimuth time-varying signals. However, in some typical application environments such as radar and sonar, although the interference incidence direction is unknown and changes rapidly with time, the azimuth characteristic of the interference incidence direction is stable within a time interval, and the utilization of the characteristic helps to simplify the implementation process of the beam forming method.

Therefore, aiming at the specific technical background, the invention firstly provides a Bayesian combined reconstruction method for extracting interesting target information by comprehensively utilizing observation data obtained by different models, namely a non-parametric Bayesian method-Dirichlet process, so as to enhance the adaptability to signal environments such as array manifold mismatch, small samples and the like, effectively utilize the intrinsic relation of common spatial structure among observation data at different moments, and further improve the precision of the data reconstruction method based on a single observation model. Then, the invention adopts an alternate iterative optimization method to realize the joint estimation of the filter parameters. The variational inference algorithm can be used for realizing synchronous optimization of multiple groups of parameters, and the effectiveness of the variational inference algorithm is verified in processing various complex hierarchical probability models including non-conjugate prior distribution. The advantage of the variation inference method in parameter estimation is represented by the reduction of local extremum solution brought by convex approximation, and under the condition of proper sample number and signal-to-noise ratio, the estimation precision better approaches to the corresponding theoretical lower bound on the basis of ensuring global convergence. Compared with the traditional probability sampling method, the variational inference method adopted by the invention has higher calculation efficiency.

The invention provides a Bayesian combined reconstruction method for extracting interesting target information by comprehensively utilizing observation data obtained by different models, namely a nonparametric Bayesian method-Dirichlet process, so as to enhance the adaptability to a motion interference environment and effectively improve the precision of a data reconstruction method based on a single observation model. In the parameter estimation process, the invention adopts an alternate iterative optimization method to realize the joint estimation of the filter parameters, the advantages of the invention in the parameter estimation aspect are that the local extremum solution brought by convex approximation is reduced, and under the condition of proper sample number and signal-to-noise ratio, the estimation precision better approaches to the corresponding theoretical lower bound on the basis of ensuring global convergence. The technical approach adopted by the invention can automatically cluster the observation data according to the structural characteristics of the observation data (the category attribute is identified by a mixed factor in the Dirichlet process), so that the intrinsic relation of a common spatial domain structure among the observation data at different moments can be effectively utilized, and the reconstruction precision of the Bayesian weight coefficient is improved.

Because the technical approach adopted by the invention can automatically cluster the observation data according to the structural characteristics of the observation data (the category attribute is identified by the mixed factor in the Dirichlet process), the obtained technical result also has a great reference value for relevant researches in the fields of multi-task learning, frequency estimation under the non-uniform sampling condition, target detection and tracking based on a multi-observation platform and the like.

The basic principles and embodiments of the present invention were verified by computer numerical simulations.

Drawings

The drawings are only for purposes of illustrating particular embodiments and are not to be construed as limiting the invention, wherein like reference numerals are used to designate like parts throughout.

Fig. 1 is a graph showing the variation of the output SINR of the array with the input SNR obtained by the method of the present invention and five other robust beamforming methods in an implementation example under the situation of mismatching of the steering vector of the desired signal caused by the DOA estimation error, and a comparison graph between the output SINR and the optimal output SINR curve.

Fig. 2 is a graph showing the variation of the array output SINR with the input SNR obtained by the method of the present invention and five other robust beamforming methods in an implementation example under the condition of the mismatch of the desired signal steering vector caused by coherent local scattering, and a comparison graph of the output SINR with the optimal output SINR curve.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

The technical scheme adopted by the invention for solving the technical problem can be divided into the following 5 steps:

the method comprises the following steps: using array element position coordinates as d₁,…,d_NThe N-element uniform linear array is used as a receiving array to receive far-field narrow-band expected and interference signals. Each sensor array element on the linear array converts the received physical signal into an electric signal and obtains a discrete time domain signal through an amplifying circuit and a data acquisition unit.

Step two: assuming that the fast array sampling rate is K, the array received data at each sampling time can be expressed as:

in the formula, y_kV and n_kAll are N × 1 dimensional vectors, which respectively represent the array output signal, the target component and the interference + noise component at the sampling time k. s_kThe amplitude of the desired signal waveform at sampling instant k. The target, interference and noise components are assumed to be independent of each other. v is a desired signal steering vector expressed asWherein λ₀For the wavelength of the incident signal, { d₁,d₂,…,d_NIs the position coordinate of the array element, theta₀DOA of desired signal, Y ═ Y₁,…,y_K]Representing the array received data set over the observation time

Step three: establishing a hierarchical probability framework according to the array received signal model established in the second step, wherein the hierarchical probability framework comprises the following sub-steps:

the first substep: describing the azimuthally varying characteristics of the interference + noise component by the Dirichlet process, i.e.

n_k～G

In the formula, the-representation obeys,is a discrete random variable and is supported in its support setThe sum of the probabilities of (c) is 1, δ (·) is a dirac function,π_f～Beta(1,γ)∝γ(1-π_f)^γ-1and oc represents proportional to, γ represents the scaling factor and obeys the gamma distribution as follows:

and a second substep: the f-th set of interference + noise components share a precision matrix, namely:

in the formula G₀Denotes the reference distribution, Λ_fExpressing the precision matrix corresponding to the f-th group of Gaussian distribution, and obeying the Wishart distribution

p(Λ_f)＝W(Λ_f|W，v)∝|Λ_f|^v-Netr{-W^-1Λ_f}

In the formula (I), the compound is shown in the specification,representing the mean matrix, v representing degrees of freedom;

and a third substep: introducing an allocation vector z_kWhich is a random variable extracted from a multi-valued distribution (the parameter set of the multi-valued distribution is { ω }_f}_{f＝1，…，K}) Namely:

in the formula z_kIs a K × 1 dimensional vector with only one element being 1;

and a fourth substep: definition 1[ z_k＝f]Denotes z_kThe f-th element in (1) is 1, i.e. n_kAnd if the class f is assigned, the likelihood function of the array observation data follows the following complex Gaussian distribution:

in the formulaIs thatA corresponding accuracy matrix;

and a fifth substep: the statistical properties of the desired signal amplitude are represented in a complex gaussian distribution:

wherein s ═ s₁,…,s_K]^TIs a vector of the amplitudes of the desired signal waveforms for the K snapshots,the desired signal accuracy;

and a sixth substep: desired signal accuracyObey the following gamma distribution

Wherein, gamma function is represented by gamma, a > 0 represents shape parameter, b > 0 represents expansion parameter;

and a seventh substep: the desired signal steering vector follows the following Watson distribution:

in the formula (I), the compound is shown in the specification,for the regularization parameter, λ is the focusing parameter. c. C_pThe hypergeometric function in (λ) is defined as the following form:

the hyperparameters mu, lambda obey the following Watson-Gamma joint distribution

And a substep eight: to aggregateAll unknown variables in the probability model are represented, and the statistical characteristics of the unknown variables are combined with the distribution of the observation data Y to obtain the following joint probability density:

step four: calculating posterior distribution parameters of each variable in theta by using a variational inference method, comprising the following substeps:

the first substep: initializing the probabilistic model parameters, i.e. let v ═ μ ═ m₀＝v_pWherein v is_pFor a preset desired signal steering vector, Λ_f＝I_N×N(f＝1,…,K)，π_f(f＝1,…,K)＝0.5，v＝N，a＝a₀＝b＝b₀＝c＝d＝10^-6，s＝1_K×1，W＝10⁺⁶·I_N×NIn which 1 is_K×1Is a K x 1 dimensional full 1 vector, I_N×NAn identity matrix of dimension NxN;

and a second substep: f group interference + noise component Λ_fA posteriori probability q (Λ)_f) Subject to the Wishart distribution with its degrees of freedomSum-mean matrixThe update formulas of (a) and (b) are respectively:

wherein, is defined as < q (z)_k＝f)〉＝φ_k,fAnd (h) represents the mathematical expectation operation. Phi is a_k,fThe calculation formula of (2) is as follows:

in the formula (I), the compound is shown in the specification,

and a third substep: pi_fObey a beta distribution, i.e. q (pi)_f)＝Beta(π_f|c_f,d_f) Wherein the hyperparameter c_fAnd d_fThe update formulas of (a) and (b) are respectively:

in the formula (I), the compound is shown in the specification,

and a fourth substep: the posterior distribution of the scaling factor gamma is the gamma distribution, i.e.Wherein the shape parameterAnd a parameter of expansion and contractionThe update formulas of (a) and (b) are respectively:

in the formula < ln (1-pi)_f)〉＝Ψ(d_f)-Ψ(c_f+d_f) Ψ (. cndot.) is a digamma function;

and a fifth substep: amplitude of desired signal waveform in k-th snapshotPosterior probability of degree q(s)_k) Subject to complex Gaussian distributions, i.e.Mean value thereofAnd precisionThe iterative update formula of

And a sixth substep: the posterior probability of the desired signal steering vector v obeys a Watson distribution, i.e.Wherein the hyperparameterAndare respectively a matrixAnd a principal eigenvector and a principal eigenvalue of, and

and a seventh substep: the posterior distribution of μ is the Watson distribution, i.e.WhereinIs a positive definite matrix (beta mm)^H+vv^H) Is determined by the maximum characteristic value of the image,is its corresponding feature vector;

and a substep eight: the posterior distribution of λ is a Gamma distribution, i.e., q (λ) ═ Gam (λ | a)₁,b₁) Wherein the shape parameter a₁And a parameter of expansion and contraction b₁The update formulas of (a) and (b) are respectively:

in the above formula, the first and second carbon atoms are,calculated from the last iteration step;

and a substep nine: the posterior probability of the accuracy of the desired signal obeys the gamma distribution, i.e.Wherein the shape parameter a₂And a parameter of expansion and contraction b₂The update formulas of (a) and (b) are respectively:

a₂＝K+a

and a substep of ten: iterating the substeps two to nine until a convergence condition is met, i.e.WhereinThe desired signal accuracy value calculated for the current step,calculating an expected signal accuracy value for the previous step;

step five: converging and solving the f group interference + noise precision matrix and the expected signal guide vectorAnd v_optSubstitution formulaCalculating to obtain the optimal weight vector corresponding to the f-th array receiving dataAnd obtaining the optimal weight vectors corresponding to the received data of the other groups of arrays by adopting the same calculation method.

And carrying out numerical simulation by using a computer to check the estimation performance of the method provided by the invention.

The simulation adopts a 10-element uniform linear array with half-wavelength array element spacing, the received noise of each array element obeys complex Gaussian distribution with the mean value of 0 and the variance of 1, and the complex Gaussian distribution is independent of each other. The total number of the incident signals of the array is 3, wherein the number of the expected signals is 1, the number of the equal-power motion interference signals is 2, and each incident signal waveform obeys complex Gaussian distribution. The Interference-TO-noise ratio (INR) is set TO 30dB, and the change rule of the incident DOA of the Interference signal along with the sampling time is theta₁(k) 30 ° +0.1 ° k and θ₂(k) 60 ° -0.1 ° k. The preset value of the desired signal incident DOA is 5 °. The number of Monte Carlo experiments is set as 100, and in each experiment, the position error of each array element obeys [ -0.05 lambda ]₀,0.05λ₀]A uniform distribution in the range, wherein₀Is the wavelength of the incident signal. The calculation formula of the SINR output by the kth sampling time array isWherein w_kAnd R_i+n(k) Weight vectors each representing a sampling instant kThe optimal weight vector can be obtained by using an Interference-plus-noise covariance (INC) matrix and the optimal array output SINR of the sampling time kAnd solving by substituting the formula.

1) When the DOA estimation error exists, the output SINR results of the method provided by the invention and the five existing stable beam forming methods are compared

Considering the mismatching situation of the direction vector of the expected signal caused by the DOA estimation error, namely, the error between the preset value and the real value of the DOA of the expected signal in each Monte Carlo experiment obeys [ -3 degrees, 3 degrees ]]Uniform distribution within the range. The array sample fast beat number is 50. The proposed method is utilized in combination with the five existing methods, namely CMT method (J.R. Guerci, the Theory and application of innovative matrix properties for ROBUST adaptive beamforming, IEEE trans. Signal Process.47(4) (1999) 977-.Adaptive beamforming algorithms with ROBUST aggregation of the messages from the beginning of the game, IEEE Trans. Signal Process.45(7) (1997) 1878-.An external performance of Adaptive beamforming in a source environment with a heated array and a moving interface source, IEEE Transmission Signal Process.48(1) (2000) 246-. FIG. 1 shows the beam forming results of six methods compared with the optimal SINRFigure (a). From the results shown in fig. 1, it can be seen that the proposed method has better interference suppression performance than the existing method due to higher INC matrix and expected signal steering vector estimation accuracy.

2) When coherent local scattering exists, the output SINR results of the method provided by the invention and the five existing robust beam forming methods are compared

The true desired signal steering vector may be expressed asWhereinA vector is directed for a preset desired signal,representing coherent scatter paths, η_iAnd i is 1,2,3 and 4, which are corresponding phases on each path. In each Monte Carlo experiment, the parameter { theta }_0i}，{η_iIndependently of the same distribution, where { theta } is_0iSubject to a Gaussian distribution with a mean of 3 and a variance of 1 °, (η)_iObey [0,2 pi ]]Uniformly distributed therein. Setting the fast beat number to be 50, selecting five methods of CMT, ROBUST LSMI, ROBUST EP, WORST CASE and MDDR as comparison objects of the method provided by the invention, and respectively calculating the beam output SINR obtained by each method under different input SNR conditions. Fig. 2 shows the output SINR and the optimum SINR versus input SNR curves of the above six methods. From the results shown in fig. 2, it can be seen that the proposed method has better interference suppression performance than other existing methods, and reflects that the idea of utilizing signal characteristics in the modeling process is advantageous, and the idea can effectively utilize the spatial domain non-stationarity of the motion interference.

While the invention has been described with reference to specific embodiments, the invention is not limited thereto, and various equivalent modifications or substitutions can be easily made by those skilled in the art within the technical scope of the present disclosure.