Target direction finding method based on expectation maximization in multi-noise environment

Abstract

The invention discloses a target direction finding method based on expectation maximization in a multi-noise environment, and the method comprises the steps: constructing a noisy observation model y = diag (1 + e) a (w) x + n, e represents multiplicative noise, a (w) represents an array flow pattern vector, x represents signal intensity, and n represents additive noise; the expectation maximization algorithm and the optimization minimization algorithm are used for solving the frequency w and the intensity x of the signal for the noisy observation model. The method can jointly estimate the frequency and the intensity of the signal under the condition that additive noise and multiplicative noise exist; according to the method, the expected maximum value and the optimization minimization method are used for jointly estimating the frequency and the intensity, the estimation accuracy is high, and the convergence can be ensured.

Description

Technical Field

[0001] The invention belongs to the field of target direction finding methods, and in particular relates to a target direction finding method based on expectation maximization in a multiple noise environment. Background Art

[0002] Target direction finding is a very important problem in the field of signal processing, and has many applications in fields such as radar, sonar, wireless sensor networks, etc. In general, existing target direction finding methods can be divided according to different types of observations, such as arrival time, arrival time difference, received signal strength, or even a mixture of these.

[0003] In the problem of target direction finding using observations based only on range, most of the work focuses on additive noise, where the observation error is independent of the distance from the sensor to the target. Based on this problem, some scholars proposed a least squares approach based on range measurements (R-LS), which uses a semidefinite relaxation technique. Other scholars proposed a least squares approach based on squared range measurements (LS based on squared range measurements, SR-LS), which uses a bisection algorithm to find the global optimal solution, but both methods have high computational complexity.

[0004] In reality, the distance observation error will gradually increase with the increase of distance. Therefore, if the model only models the additive noise, it will lead to inaccurate positioning. The existence of multiplicative noise must also be considered. In order to solve the problem of target direction finding related to observation and distance in the presence of both additive and multiplicative noise, some scholars proposed a Newton-Raphson iterative algorithm to solve the approximate maximum likelihood estimation problem. For target tracking, the initial point is always very good, and the Newton-Raphson iterative algorithm has very good convergence. However, for target direction finding, when the initial point is not very good, the algorithm performs poorly. Summary of the invention

[0005] The object of the present invention is to provide a target direction finding method based on expectation maximization in a multiple noise environment, aiming at the target direction finding problem related to observation and distance in the presence of both additive and multiplicative noises, and the method has low computational complexity and high estimation accuracy.

[0006] To achieve the above object, the present invention provides the following technical solutions:

[0007] A target direction finding method based on expectation maximization in a multiple noise environment, the method comprising:

[0008] Construct a noisy observation model y = diag (1 + e) ​​a (w) x + n, where y ∈ C M×1 Noisy observation vector, M represents the number of elements in the uniform linear array, diag() represents the diagonal matrix, e∈C M×1 represents multiplicative noise, e={e1,...,e m , ..., e M}, w is the signal frequency, a(w) represents the array flow vector, x∈C 1×1 represents the signal strength, n represents the additive noise, n={n1,...,n m , ..., n M}, CN represents complex Gaussian distribution, μ e，m represents the mean, and denote the variance of multiplicative noise and additive noise respectively;

[0009] The frequency w and intensity x of the signal are solved for the noisy observation model using the expectation maximization algorithm and the optimization minimization algorithm, including:

[0010] Step 1: Set the initial point w of signal frequency and strength 0 、x 0 , parameter ε and the maximum number of iterations I max ;

[0011] Step 2, set t=0, a=0;

[0012] Step 3: Calculate the hidden variable h m =1+e m About the likelihood function p(h m |y m ;x k , w k ) and the expected μ post,m ,

[0013]

[0014]

[0015] in, j is an imaginary number, k represents the current iteration number;

[0016] Step 4: Restrict w to a finite discrete set Ω, maximizing the function max w |a H (w)λ| 2 w cAs output, H represents the conjugate transpose,

[0017] Step 5: Determine if a<3 or t≤=I max , if yes, execute the loop strategy, otherwise go to step 6;

[0018] Step 6, output w t , x t .

[0019] Preferably, the discrete set γ represents the oversampling ratio.

[0020] Preferably, the circulation strategy includes:

[0021] Step 5.1, maximize the function max w |a H (w)λ| 2 Use Newton's method to correct the frequency w and determine the objective function g after frequency correction t Less than before correction, if so, accept the current frequency correction result, otherwise the current frequency remains unchanged;

[0022] The objective function is:

[0023]

[0024] z represents the equivalent noise, represents the variance of the equivalent noise, △ represents the mean of the equivalent noise; Step 5.2, update

[0025]

[0026] Step 5.3, calculate the objective function g t , and judge whether t≥=1 or |g t -g t-1 |<ε, if so, then a=a+1 and execute step 5.4, otherwise a remains unchanged and execute step 5.5 directly;

[0027] Step 5.4, Update and μ post，m ;

[0028] Step 5.5, t=t+1.

[0029] Compared with the prior art, the present invention has the following beneficial effects:

[0030] The present invention can jointly estimate the signal frequency and strength in the presence of both additive and multiplicative noise;

[0031] The present invention utilizes the expected maximum value and the optimization minimization method to jointly estimate the frequency and intensity, and has high estimation accuracy and can ensure convergence;

[0032] The present invention utilizes the optimization maximization algorithm to execute the M step of the expectation maximization algorithm, with a small number of iterations, low computational complexity, and easy hardware implementation;

[0033] The present invention uses Newton's method to estimate the signal frequency, effectively reduces the impact of the grid point mismatch problem, and improves the frequency estimation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 This is the computational flow chart of the noisy observation model. DETAILED DESCRIPTION

[0035] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work belong to the protection scope of the present invention.

[0036] Additive noise generally refers to thermal noise, shot noise, etc., and their relationship with the signal is addition. Noise exists regardless of whether there is a signal or not. Multiplicative noise is generally caused by an imperfect channel. Their relationship with the signal is multiplication. It exists when the signal is there, and it does not exist when the signal is not there. In general communications, additive randomness is regarded as the background noise of the system; while multiplicative randomness is regarded as caused by the time-varying nature of the system (such as fading or Doppler) or nonlinearity.

[0037] For a uniform linear array consisting of M sensors, only one signal is of interest, and the number of snapshots of the signal T is also 1. The noisy observation model is modeled as shown in formula (1):

[0038] y=diag(1+e)a(w)x+n (1)

[0039] Where y∈C M×1 Noisy observation vector, M represents the number of elements in the uniform linear array, diag(1+e) represents a square matrix with 1+e elements in the diagonal position and 0 in other positions, e∈C M×1 represents multiplicative noise, e={e1,...,e m , ..., e M}, x∈C 1×1 represents the signal strength, n represents the multiplicative noise, n={n1,...,n m , ..., n M ), CN represents complex Gaussian distribution, μ e，m represents the mean, and represents the variance, a(w) represents the array flow vector, w is the signal frequency, j is an imaginary number, and T is the transpose.

[0040] For a noisy observation vector y, any element y m As shown in formula (2):

[0041]

[0042] By defining the equivalent noise z m :

[0043]

[0044] Change formula (2) into:

[0045]

[0046] In the present invention, the signal frequency and signal strength of the noisy observation model are solved, and the frequency and strength of the signal are used to perform target direction finding.

[0047] In the present invention, because e m and n m is a Gaussian variable that follows a complex Gaussian distribution,

[0048]

[0049] Then the probability density function is

[0050]

[0051] In the formula, △ represents the mean value, represents the variance of the equivalent noise,

[0052] The maximum likelihood estimation problem corresponding to formula (1) can be modeled as shown in formula (5):

[0053] min w,x ~g(w,x) (5)

[0054] in

[0055]

[0056] The above problem can be solved by gradient descent or Newton-Raphson method, where the step size can be chosen by using the backward line search technique in each iteration, but this leads to a slow convergence rate.

[0057] The present invention provides a new iterative method called R-EM-MM method, which combines the expectation maximization algorithm EM with the optimization minimization algorithm MM, making full use of the structure of the positioning problem to ensure convergence. The principle of the method is as follows:

[0058] 1+e m Set as hidden variable h m , as shown in formula (7):

[0059] h m =1+e m , m=1…,M.(7)

[0060] Then formula (2) becomes:

[0061]

[0062] For the kth iteration, let w k and x k denotes the estimated values ​​of w and x.

[0063] In the expectation maximization algorithm, the E-step determines the entire likelihood function log p(y, h; w, x) relative to the posterior probability density function p(y|h; x k , w k ), that is:

[0064] Q(w,x;w k , x k ) = E[log p(y, h; w, x)]

[0065] =E[log p(y|h;x,w)]+const (9)

[0066] in,

[0067]

[0068] * represents conjugate, R represents the real part operation, and const represents a set of constants that are independent of variables.

[0069] Taking the expectation of (9), the M-step determines the maximum Q(w, x; w k , x k ) is equivalent to formula (11):

[0070]

[0071] Notice

[0072]

[0073] in,

[0074]

[0075] Therefore p(h|y;x k , w k ) is also complex Gaussian distributed, and has

[0076] Let μ post，m and Indicates h m About p(h m |y m ;x k , w k ), then we have:

[0077]

[0078] Let|h m | 2 and h m If the coefficients in the two equations agree, then we have

[0079]

[0080]

[0081] Further:

[0082]

[0083]

[0084] because

[0085]

[0086] And combined with (11), M-step is equivalent to optimizing and minimizing g(w, x), as shown in formula (21):

[0087] min w，x g(w,x) (21)

[0088] in,

[0089]

[0090] This means

[0091]

[0092] therefore

[0093]

[0094] Substituting x into g(w,x), we have:

[0095]

[0096]

[0097] It can be seen It is unrelated to w, so the penalty function for w is as shown in formula (26):

[0098] max w |a H (w)λ| 2 (26)

[0099] in, H is the conjugate transpose.

[0100] In the present invention, the R-EM-MM method is used to solve the frequency and intensity of the signal in the noisy observation model, including the following steps:

[0101] Step 1: Set the initial point w of signal frequency and strength 0 、x 0 , parameter ε and the maximum number of iterations I max ;

[0102] Step 2, set t=0, a=0;

[0103] Step 3: Calculate the hidden variable h m About the likelihood function p(h m |y m ;x k , w k ) and the expected μ post,m ;

[0104]

[0105]

[0106] Step 4: Restrict w to a finite discrete set 0, 1, ..., (γM-1)}, γ represents the oversampling rate, which will maximize the function max w |a H (w)λ| 2 w c As output, w c The value w that maximizes the function among all grid point frequencies;

[0107] Step 5: Determine if a<3 or t≤=I max If yes, execute the loop strategy, otherwise output w t , x t ;

[0108] The cycle strategy includes:

[0109] (5-1) For the maximization function max w |a H (w)λ| 2 Use Newton's method to correct the frequency w and determine the target function g after frequency correction t Less than before correction, if so, accept the current frequency correction result, otherwise the current frequency remains unchanged;

[0110] (5-2) Update

[0111]

[0112] (5-3) Calculate the objective function g t , and judge whether t≥=1 or |g t -g t-1 |<ε, if so, then a=a+1 and execute step (5-4), otherwise execute step (5-5)

[0113] Among them, the objective function g t It is expressed as:

[0114]

[0115] (5-4) Update and μ post,m , execute step (5-5);

[0116] (5-5)t=t+1.

[0117] In the present invention, how to use Newton's method to correct the current frequency is a conventional technical means in the field, and those skilled in the art can set it according to actual conditions, which will not be described in detail here.

[0118] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A target direction finding method based on expectation maximization in a multiple noise environment, characterized in that: The method comprises: Construct a noisy observation model y = diag (1 + e) ​​a (w) x + n, where y ∈ C M×1 Noisy observation vector, M represents the number of elements in the uniform linear array, diag() represents the diagonal matrix, e∈C M×1 represents multiplicative noise, e={e1,...,e m , ..., e M }, w is the signal frequency, a(w) represents the array flow vector, x∈C 1×1 represents the signal strength, n represents the additive noise, n={n1,...,n m , ..., n M }, CN represents complex Gaussian distribution, μ e，m represents the mean, and denote the variance of multiplicative noise and additive noise respectively; The frequency w and intensity x of the signal are solved for the noisy observation model using the expectation maximization algorithm and the optimization minimization algorithm, including: Step 1: Set the initial point w of signal frequency and strength 0 、x 0 , parameter ε and the maximum number of iterations I max ; Step 2, set t=0, a=0; Step 3: Calculate the hidden variable h m =1+e m About the likelihood function p(h m |y m ;x k , w k ) and the expected μ post，m , in, j is an imaginary number, k represents the current iteration number; Step 4: Restrict w to a finite discrete set Ω, maximizing the function max w |a H (w)λ| 2 w c As output, H represents the conjugate transpose, Step 5: Determine if a<3 or t≤=I max , if yes, execute the loop strategy, otherwise go to step 6; Step 6, output w t , x t .

2. A target direction finding method based on expectation maximization in a multiple noise environment as claimed in claim 1, characterized in that: The discrete set γ represents the oversampling ratio.

3. The target direction finding method based on expectation maximization in a multiple noise environment as claimed in claim 1, characterized in that: The cycle strategy includes: Step 5.1, maximize the function max w |a H (w)λ| 2 Use Newton's method to correct the frequency w and determine the objective function g after frequency correction t Less than before correction, if so, accept the current frequency correction result, otherwise the current frequency remains unchanged; The objective function is: z represents the equivalent noise, represents the variance of equivalent noise, △ represents the mean of equivalent noise; Step 5.2, Update Step 5.3, calculate the objective function g t , and judge whether t≥=1 or |g t -g t-1 |<ε, if so, then a=a+1 and execute step 5.4, otherwise a remains unchanged and execute step 5.5 directly; Step 5.4, Update and μ post,m ; Step 5.5, t=t+1.

Images