CN115015831B - Large-scale array target azimuth estimation method under combined influence of impulse noise and non-uniform noise - Google Patents

Abstract

The application provides a large-scale array target azimuth estimation method under the combined influence of impulse noise and non-uniform noise, and belongs to the field of underwater acoustic array signal processing. The application builds a variation sparse Bayesian signal model based on an underwater sound array, builds an S-t distributed impulse noise model, adopts Bernoulli distribution to detect the existence probability of impulse noise, the designed model is more close to real application, and the background noise simulates non-uniform Gaussian noise by non-uniform Gaussian distribution, thereby providing a high-precision direction-of-arrival estimation method under the mixed condition of non-uniform noise and impulse noise. Compared with the existing similar wave large-direction estimation method, the method has the advantages of higher estimation precision and stronger adaptability.

Description

Large-scale array target azimuth estimation method under combined influence of impulse noise and non-uniform noise

Technical Field

The invention relates to a large-scale array target azimuth estimation method under the combined influence of impulse noise and non-uniform noise, and belongs to the field of underwater acoustic array signal processing.

Background

The DOA estimation is a popular research direction in array signal processing, and in the field of underwater acoustic detection, underwater acoustic signals are observed or received by an array formed by a plurality of hydrophones, so that the DOA estimation of an underwater target is realized. The traditional beam forming methods such as CBF, MVDR and the like have limited resolution, are not enough to meet the requirement of high precision, and subspace azimuth estimation methods such as MUSIC, ESPRIT and the like can obtain azimuth estimation with high resolution, but require the number of information sources and the like as prior information.

Because the ocean environment is more complex and changeable, the interference of sea water flow velocity, high pressure, high waves and the like can be encountered when the array operation is carried out under water, the environment where the array is positioned cannot be directly assumed to be Gaussian noise environment, non-uniform Gaussian noise can exist among array elements of the array, impulse noise exists sometimes, and the direction finding precision of the traditional method is reduced due to the complex noise environment. Although sparse estimation methods for impulse noise exist at present, the situation that non-uniform noise and impulse noise exist at the same time cannot be processed.

Disclosure of Invention

The invention aims to provide a large-scale array target azimuth estimation method under the combined influence of impulse noise and non-uniform noise.

The purpose of the invention is realized in the following way: the method comprises the following steps:

step one: building an underwater uniform long-line array, uniformly distributing the array at half-wavelength intervals during distribution, and comprehensively building an incident signal model received by each array element module:

wherein Y is an observation signal, The method is characterized in that the method is an overcomplete matrix, X is a sparse representation matrix of an incident signal, Z is a pulse noise judgment matrix, V is a Gaussian noise matrix, 1 _M×L is a matrix with all values of 1, and E is a pulse noise matrix;

Step two: constructing prior distribution of each variable in the signal model, and carrying out hierarchical prior distribution on the array received signals:

Constructing a hierarchical prior distribution of the array received signals at the first moment is as follows:

Wherein, H, η represents a hidden variable matrix for controlling impulse noise, β represents a non-uniform noise variance vector, a subscript (·) _m,l represents a matrix element of an mth array element at a first moment, a subscript (·) _l represents a matrix vector of a first moment, a subscript (·) _m represents a matrix vector of an mth array element, for example, Y _m,l represents data received by the mth array element at the first moment, Y _l represents received data of the entire array at the first moment, and Y _l at all moments together form a data matrix Y; the meanings of the variables with the same subscript and case distinction are the same;

defining the noise variance beta as a hierarchical Gamma distribution:

Wherein g _m and h _m respectively represent a shape parameter and an inverse scale parameter of the noise variance beta _m at the m-th array element;

Defining the desired signal as a complex gaussian distribution:

where α is the variance matrix of the desired signal, and subscript (·) _n represents the matrix vector at the nth scanning orientation, e.g., α _n is the desired signal variance at the nth scanning orientation;

Respectively constructing layered Gamma distribution for the variance matrix alpha and the hidden variable matrix elements eta _l and H _l at the first moment:

Wherein a _n,ξ_m,l/2,c_m,l and gamma _m,l respectively represent shape parameters of corresponding distribution, b _n,ξ_m,l/2,d_m,l and rho _m,l respectively represent inverse scale parameters of corresponding distribution, ζ _l is variance vector of constraint hidden variable matrix element eta _l, and alpha is variance matrix of desired signal;

the bernoulli distribution is constructed for the noise occurrence state vector z _l at the first time:

Wherein pi _l represents the occurrence probability vector of the state vector z _l, the corresponding lower case letter is a specific element corresponding to the position of the subscript, and the occurrence probability pi _l is constructed into a layered Beta distribution as follows:

wherein, p _m,l and q _m,l are respectively Beta distribution parameters obeyed by the m-th array element at the first moment;

Step three: parameter initialization: setting initial iteration iter=1, initializing maxiter,N,a₀,b₀,c₀,d₀,g₀,h₀,γ₀,ρ₀,z_m,l,maxiter as iteration maximum times, N as the number of grids divided in the whole space azimuth range, a ₀ as the shape parameter of the distribution variance of the desired signal, b ₀ as the inverse scale parameter of the distribution variance of the desired signal, c ₀、g₀ and gamma ₀ as the shape parameter of the noise variance distribution, d ₀、h₀ and rho ₀ as the inverse scale parameter of the noise variance distribution, and p ₀ and q ₀ as the control occurrence probability pi _l;

Step four: iterative calculation, namely obtaining an iterative formula of each parameter according to the distribution in the second step, and then updating;

Step five: judging whether an iteration condition is met, if so, jumping out of iteration and outputting mu _x, and if not, continuing to carry out the step four;

Step six: spectral peak searching, and carrying out azimuth estimation according to the restored signal:

Step seven: and outputting an azimuth estimation result P _VB (theta).

The invention also includes such structural features:

1. The fourth step specifically comprises:

Updating

Updating

Update a _n:

update b _n:

Update c _m,l:

Update d _m,l:

Updating

Updating

Update g _m,l:

update h _m,l:

update gamma _m,l:

Update ρ _m,l:

Updating the average value of each variable:

Let iter=iter+1.

2. Updating varianceThe update variable < D _Δ > at the time is:

D_Δ＝diag(z_l⊙β_l+(1_M×1-z_l)⊙η_l⊙H_l)

wherein diag (·) is a diagonal operation.

3. The iteration termination condition of the fifth step is as follows:

iter≥maxiter

wherein, tol is a set termination threshold, and the value is 0.001; α is the desired signal variance, defined as:

where < · > represents the desire for the variable to be solved for.

Compared with the prior art, the invention has the beneficial effects that: the method fully considers the existence of impulse noise, the influence of impulse noise and the influence of non-uniform noise, utilizes the distribution characteristic of Bernoulli distribution control noise to construct a variable decibel leaf model, enables the model to be more similar to an actual marine environment, realizes accurate azimuth estimation, has ultrahigh separate capability, multi-target resolving power and the condition of resisting the simultaneous existence of non-uniform noise and impulse noise, and is suitable for complex environment arrival azimuth estimation.

Drawings

FIG. 1 is a diagram of an array receiver signal system according to the present invention;

FIG. 2 is a model diagram constructed in accordance with the present invention;

FIG. 3 is a flow chart of the direction estimation of the present invention;

Fig. 4 (a) - (b) are diagrams comparing the present invention (proposed method) with two other methods (conventional sparse method, impulse noise based sparse method): fig. 4 (a) root mean square error (Root Mean Square Error, RMSE) results for each method; fig. 4 (b) shows the detection probability results of different methods.

Detailed Description

The invention is described in further detail below with reference to the drawings and the detailed description.

The operation steps of the invention are as follows, in combination with the accompanying drawings:

(1) As shown in FIG. 1, an M-element underwater uniform long line array is built, the array elements are uniformly distributed at half-wavelength intervals during distribution, each array element module receives signals, and a received signal model is built as follows

Wherein the method comprises the steps ofData matrix received for hydrophone array,/>For array flow pattern matrix,/>For the sampled expected signal matrix,/>State matrix appearing for impulse noise at each moment of each channel of array,/>For the non-uniform noise matrix of each array element of the hydrophone array, the addition is matrix dot product operation, 1 _M×L is the full 1 matrix,/>The pulse noise matrix is a pulse noise matrix at each array element of the hydrophone array, L is the length of the array receiving signal, and N is the search azimuth space grid number.

(2) The construction of the variational inference model is shown in fig. 2.

(2.1) Constructing a hierarchical prior distribution of the array received signals at the first moment as follows:

wherein, H, η represents a hidden variable matrix for controlling impulse noise, β represents a non-uniform noise variance vector, a subscript (·) _m,l represents a matrix element of an mth array element at a first moment, a subscript (·) _l represents a matrix vector of a first moment, a subscript (·) _m represents a matrix vector of an mth array element, for example, Y _m,l represents data received by the mth array element at the first moment, Y _l represents received data of the entire array at the first moment, and Y _l at all moments together form a data matrix Y; similarly, the meanings of the variables with the same subscript and case distinction are the same;

defining the noise variance beta as a hierarchical Gamma distribution:

Wherein g _m and h _m respectively represent a shape parameter and an inverse scale parameter of the noise variance beta _m at the m-th array element;

Defining the desired signal as a complex gaussian distribution:

where α is the variance matrix of the desired signal, and subscript (·) _n represents the matrix vector at the nth scanning orientation, e.g., α _n is the desired signal variance at the nth scanning orientation;

Respectively constructing layered Gamma distribution for the variance matrix alpha and the hidden variable matrix elements eta _l and H _l at the first moment:

Wherein a _n,ξ_m,l/2,c_m,l and gamma _m,l respectively represent shape parameters of corresponding distribution, b _n,ξ_m,l/2,d_m,l and rho _m,l respectively represent inverse scale parameters of corresponding distribution, ζ _l is variance vector of constraint hidden variable matrix element eta _l, and alpha is variance matrix of desired signal;

the bernoulli distribution is constructed for the noise occurrence state vector z _l at the first time:

Wherein pi _l represents the occurrence probability vector of the state vector z _l, the corresponding lower case letter is a specific element corresponding to the position of the subscript, and the occurrence probability pi _l is constructed into a layered Beta distribution as follows:

Wherein, p _m,l and q _m,l are respectively Beta distribution parameters obeyed by the mth array element at the first moment.

(2.2) Defining a variable set omega= { X, alpha, H, eta, zeta, beta, Z and pi } by using a variation inference theory, and substituting the distribution matrix constructed in (2.1) into the following formula in sequence to solve the posterior probability of each variable:

Wherein q (·) represents the posterior probability of the variable (·), ln (·) represents taking the logarithm, < · > represents taking the expectation, p (·|·) represents the probability of the element therein, which can be obtained from (2.1), _q(Ω_i +.q) represents calculating the part of the set that does not contain the variable · and const represents a constant term.

(3) The mean and variance of each variable are calculated from the probability distribution of each variable obtained in (2).

(3.1) Parameter initialization: setting initial iteration iter=1, initializing maxiter,N,a₀,b₀,c₀,d₀,g₀,h₀,γ₀,ρ₀,z_m,l,maxiter as iteration maximum times, N as the number of grids divided in the whole space azimuth range, a ₀ as the shape parameter of the distribution variance of the desired signal, b ₀ as the inverse scale parameter of the distribution variance of the desired signal, c ₀、g₀ and gamma ₀ as the shape parameter of the noise variance distribution, d ₀、h₀ and rho ₀ as the inverse scale parameter of the noise variance distribution, and p ₀ and q ₀ as the control occurrence probability pi _l; then sequentially carrying out iterative updating on each parameter according to the probability distribution of each variable obtained in the step (2);

(3.2) updating the variance of the desired signal in turn And mean/>

Updating

Wherein ,D_Δ＝diag(z_l⊙β_l+(1_M×1-z_l)⊙η_l⊙H_l),D_α＝diag(α_l),diag(·) denotes a diagonal operation;

Updating

(3.3) Based on the results of (3.2), updating the respective distribution parameters, including:

Update a _n:

update b _n:

Update c _m,l:

Update d _m,l:

Updating

Updating

Update g _m,l:

update h _m,l:

update gamma _m,l:

Update ρ _m,l:

(3.4) based on the results of (3.3), updating the mean value of each variable, comprising:

The average value result of each variable is the variable result of iterative estimation, such as alpha= < alpha _n >, and the other is the same.

(3.5) Update iter:

iter＝iter+1

(4) Judging whether the iteration termination condition is satisfied, wherein the iteration termination condition is as follows

iter≥maxiter

Wherein, tol is a set termination threshold, and the value is 0.001; α is the expected signal variance, and when one of the above iteration termination conditions is met, the iteration is skipped and μ _x is output, otherwise the iteration of (3.2) - (3.5) is continued;

(5): spectral peak searching, and carrying out azimuth estimation according to the restored signal:

wherein, |·| ₁ represents a matrix-norm operation, || _∞ representation matrix is an infinite norm operation of (1).

(6): And outputting an azimuth estimation result P _VB (theta).

Fig. 3 is an overall flow of position estimation. The invention fully considers the existence and non-existence of impulse noise and the influence of the impulse noise on the size and the non-uniform noise, and the constructed model is more fit with the actual ocean environment, so that a better estimation result can be obtained.

Simulation study of the invention:

Simulation conditions:

Two single-frequency pulse signals are used as incident signals, the incident directions are-30 degrees and-15 degrees respectively, the snapshot number is 100, one half of snapshots are added with non-uniform noise, the noise variance of each array element is randomly changed between [0.1,5], the other half of snapshots are added with S alpha S pulse noise, so that the generalized signal to noise ratio (GSNR) is changed between [ -10 and 20], and the traditional sparse method, the sparse method based on the pulse noise and the method provided by the invention are subjected to comparative analysis.

Fig. 4 (a) is a plot of root mean square error (Root Mean Square Error, RMSE) as a function of GSNR for each method under an impulse noise environment. Comparing the curves of the methods, the traditional sparse method can be found to have serious failure; although the sparse method based on impulse noise also has a descending trend, when non-uniform noise and impulse noise alternate, a jump phenomenon exists, so that an estimation result is unstable; the method provided by the invention has the advantages of lowest RMSE, most stability, capability of adapting to noise change conditions under different SNR, good estimation result and minimum deviation.

As shown in fig. 4 (b), the probability curve of success of detection when each method changes with GSNR under impulse noise environment is defined as success of detection within 1 ° of target deviation. Comparing the curves of the methods, the probability of success of detection of the traditional sparse method can be found to rise slowly along with the increase of SNR, that is, the method can not find the target azimuth accurately in time under the background of impulse noise and non-uniform noise; the probability of success of detection is increased relatively based on the sparse method of impulse noise, but when non-uniform noise and impulse noise alternate, a jump phenomenon exists, so that an estimation result is unstable; the method provided by the invention has the highest detection success probability, is the most stable, and has the strongest estimated function on the target.

In summary, the application provides a large-scale array target azimuth estimation method under the combined influence of impulse noise and non-uniform noise, and belongs to the field of underwater acoustic array signal processing. According to the application, a variational sparse Bayesian signal modeling based on the underwater acoustic array is constructed, an S-t distribution impulse noise model is established, the presence or absence of impulse noise is detected by adopting Bernoulli distribution, the model is more close to reality, a non-uniform Gaussian noise model with non-uniform Gaussian distribution is established, and high-precision DOA estimation under the mixed condition of non-uniform noise and impulse noise is realized. The method is more close to the actual application, and the background noise simulates the non-uniform Gaussian noise by using the non-uniform Gaussian distribution, so that a high-precision direction-of-arrival estimation method under the mixed condition of the non-uniform noise and the impulse noise is provided. Compared with the existing similar wave large-direction estimation method, the method has the advantages of higher estimation precision and stronger adaptability. The effectiveness and feasibility of the application are verified by simulation research results.

Claims (4)

1. A large-scale array target azimuth estimation method under the combined influence of impulse noise and non-uniform noise is characterized by comprising the following steps:

step one: building an underwater uniform long-line array, uniformly distributing the array at half-wavelength intervals during distribution, and comprehensively building an incident signal model received by each array element module:

wherein Y is an observation signal, The method is characterized in that the method is an overcomplete matrix, X is a sparse representation matrix of an incident signal, Z is a pulse noise judgment matrix, V is a Gaussian noise matrix, 1 _M×L is a matrix with all values of 1, and E is a pulse noise matrix;

Step two: constructing prior distribution of each variable in the signal model, and carrying out hierarchical prior distribution on the array received signals:

Constructing a hierarchical prior distribution of the array received signals at the first moment is as follows:

Wherein, H, η represents a hidden variable matrix for controlling impulse noise, β represents a non-uniform noise variance vector, a subscript (·) _m,l represents a matrix element of an mth array element at a first moment, a subscript (·) _l represents a matrix vector of a first moment, a subscript (·) _m represents a matrix vector of an mth array element, for example, Y _m,l represents data received by the mth array element at the first moment, Y _l represents received data of the entire array at the first moment, and Y _l at all moments together form a data matrix Y; the meanings of the variables with the same subscript and case distinction are the same;

defining the noise variance beta as a hierarchical Gamma distribution:

Wherein g _m and h _m respectively represent a shape parameter and an inverse scale parameter of the noise variance beta _m at the m-th array element;

Defining the desired signal as a complex gaussian distribution:

where α is the variance matrix of the desired signal, and subscript (·) _n represents the matrix vector at the nth scanning orientation, e.g., α _n is the desired signal variance at the nth scanning orientation;

Respectively constructing layered Gamma distribution for the variance matrix alpha and the hidden variable matrix elements eta _l and H _l at the first moment:

Wherein a _n,ξ_m,l/2,c_m,l and gamma _m,l respectively represent shape parameters of corresponding distribution, b _n,ξ_m,l/2,d_m,l and rho _m,l respectively represent inverse scale parameters of corresponding distribution, ζ ₁ is variance vector of constraint hidden variable matrix element eta _l, and alpha is variance matrix of desired signal;

the bernoulli distribution is constructed for the noise occurrence state vector z _l at the first time:

Wherein pi _l represents the occurrence probability vector of the state vector z _l, the corresponding lower case letter is a specific element corresponding to the position of the subscript, and the occurrence probability pi _l is constructed into a layered Beta distribution as follows:

wherein, p _m,l and q _m,l are respectively Beta distribution parameters obeyed by the m-th array element at the first moment;

Step three: parameter initialization: setting initial iteration iter=1, initializing maxiter,N,a₀,b₀,c₀,d₀,g₀,h₀,γ₀,ρ₀,z_m,l,maxiter as iteration maximum times, N as the number of grids divided in the whole space azimuth range, a ₀ as the shape parameter of the distribution variance of the desired signal, b ₀ as the inverse scale parameter of the distribution variance of the desired signal, c ₀、g₀ and gamma ₀ as the shape parameter of the noise variance distribution, d ₀、h₀ and rho ₀ as the inverse scale parameter of the noise variance distribution, and p ₀ and q ₀ as the control occurrence probability pi _l;

Step four: iterative calculation, namely obtaining an iterative formula of each parameter according to the distribution in the second step, and then updating;

Step five: judging whether an iteration condition is met, if so, jumping out of iteration and outputting mu _x, and if not, continuing to carry out the step four;

Step six: spectral peak searching, and carrying out azimuth estimation according to the restored signal:

Step seven: and outputting an azimuth estimation result P _VB (theta).

2. The method for estimating a target azimuth of a large-scale array under the combined influence of impulse noise and non-uniform noise according to claim 1, wherein the step four specifically comprises:

Updating

Updating

Update a _n:

update b _n:

Update c _m,l:

Update d _m,l:

Updating

Updating

Update g _m,l:

update h _m,l:

update gamma _m,l:

Update ρ _m,l:

Updating the average value of each variable:

Let iter=iter+1.

3. The method for estimating a target azimuth of a large-scale array under the combined influence of impulse noise and non-uniform noise according to claim 2, wherein the variance is updatedThe update variable < D _Δ > at the time is:

D_Δ＝diag(z_l⊙β_l+(1_M×1-z_l)⊙η_l⊙H_l)

wherein diag (·) is a diagonal operation.

4. The method for estimating a target azimuth of a large-scale array under the combined influence of impulse noise and non-uniform noise according to claim 1, wherein the iteration termination condition of the fifth step is:

iter≥maxiter

wherein, tol is a set termination threshold, and the value is 0.001; α is the desired signal variance, defined as:

where < · > represents the desire for the variable to be solved for.