Limited area DOA estimation method and system based on atom norm minimization

Abstract

The invention provides a limited area DOA estimation method and system based on atom norm minimization, and the method comprises the steps: limiting the DOA range of an incoming wave, and converting the limited DOA range into a spatial frequency range; designing a corresponding atom set and a corresponding atom norm minimization problem based on a limited spatial frequency range, designing a corresponding mapping operator, and converting the atom norm minimization problem into an SDP problem based on the mapping operator; the SDP problem is solved, and estimation of the spatial frequency is obtained; and converting the obtained spatial frequency estimation into an estimated value of the DOA. The method has the advantages that compared with a grid algorithm, the DOA estimation precision is higher due to the fact that the method is not affected by grid mismatch; compared with an original algorithm CA-ANM, the method is consistent with the original algorithm CA-ANM in precision and performance, but the operation time is greatly shortened.

Description

Technical Field

[0001] This application belongs to the field of underwater acoustic target detection, and specifically relates to a method and system for DOA estimation in a limited area based on atomic norm minimization. Background Art

[0002] The estimation of the direction of arrival (DOA) of underwater targets is one of the hot issues in underwater acoustic array signal processing. In order to break through the Rayleigh limit and achieve super-resolution DOA estimation, subspace methods represented by multiple signal classification (MUSIC) and estimating signal via rotational invariance techniques (ESPRIT) have been widely applied to the DOA estimation problem. Although subspace methods have greatly improved the resolution of DOA estimation, they require accurate estimation of the signal covariance matrix, resulting in the inability of such algorithms to be applied to scenarios with insufficient snapshot numbers. In contrast, algorithms based on compressive sensing have gradually become a research hotspot because they can achieve high-resolution DOA estimation under the condition of a small number of snapshots. Early compressive sensing algorithms were based on grid dictionaries, that is, the DOA of the target was divided into discrete grids, and the steering vectors corresponding to each grid point could form an over-complete dictionary. Since the target DOA only falls on a limited number of grids, that is, it has sparse characteristics, this makes the DOA estimation problem transformed into a sparse recovery problem. In actual engineering, there is usually a certain prior knowledge about the DOA range of the target before DOA estimation, so it is only necessary to perform DOA estimation within the limited range. In fact, in the presence of target prior information, the computational complexity can be reduced by restricting the DOA search range to the target area. In this invention, this situation is called DOA estimation in a limited area. Yang et al. designed a nulling matrix filter to limit the DOA of the incoming signal to the target area and gave a DOA estimation method based on sparse spectrum fitting. On this basis, Zhang et al. restricted the DOA search range by limiting the dictionary range to the target area and achieved DOA estimation based on the sparse Bayesian learning (SBL) algorithm, improving the computational efficiency.

[0003] The above sparse method assumes that the target DOAs all fall on the pre-defined grid. In reality, the target DOAs can fall outside the grid, resulting in the grid mismatch problem. Therefore, subsequent research has proposed a series of compensation methods to reduce the impact of grid mismatch. Such methods are collectively referred to as off-grid algorithms. However, these methods still assume the existence of the grid, resulting in the continuous-valued DOAs still being unable to be accurately estimated. The gridless algorithms represented by the atomic norm minimization (ANM) method have received extensive attention in recent years because they can obtain sparse solutions on a continuous dictionary. Since it involves the optimization of an infinite number of non-enumerable parameters, the ANM problem itself cannot be directly solved. Therefore, Candès et al. derived the dual problem of ANM and transformed the constraint of the dual atomic norm into a triangular polynomial inequality constraint. Using the bounded lemma, the triangular polynomial inequality constraint can be transformed into a semi-definite constraint, thus transforming the dual problem into a semi-definite programming (SDP) problem, which can then be solved by programming.

[0004] Although the bounded lemma provides a method for solving the ANM problem, it also limits the application scenarios of the ANM method. To limit the DOA search range to the target area to reduce the algorithm complexity, by analogy with the method of restricting the discrete dictionary to the target area in the grid-based algorithm, the DOA search range can be restricted by restricting the constructed atom set to the target area at this time. However, this will make the corresponding dual atomic norm constraint not satisfy the form of the triangular polynomial inequality constraint in the full frequency band, resulting in the inapplicability of the bounded lemma. This makes the gridless DOA estimation method based on ANM can only search for the target DOA in the entire area, resulting in computational redundancy. Although there have been some attempts to broaden the application range of the ANM algorithm in recent years, these algorithms are only applicable to the corresponding specific scenarios and do not have universality, so they cannot be applied to the DOA estimation problem in the limited area considered in the present invention.

[0005] In addition, although the SDP problem obtained through bounded lemma transformation is equivalent to the atomic norm minimization problem, due to its scale and other issues, existing general solvers, such as the CVX toolkit, have difficulty solving it efficiently. To improve the solution efficiency, Hansen et al. proposed a fast interior point method based on the primal-dual interior point method. By analytically solving the augmented Karush-Kuhn-Tucker (KKT) conditions and using the pseudo-Newton method L-BFGS to obtain the update step size, the semi-positive constraints are replaced by finite-point trigonometric polynomial inequality constraints, and the fast Fourier transform (FFT) is used to determine the constraints, which significantly improves the solution speed. However, the above algorithm has strict requirements on the form of the signal. Although Gao et al. extended the fast interior point method to the case of partial observation of MMV, its application scenarios are still limited. Summary of the Invention

[0006] The purpose of this application is to overcome the defect of low efficiency in solving SDP problems in the prior art.

[0007] To achieve the above objectives, this application proposes a limited area DOA estimation method based on atomic norm minimization, including:

[0008] Step 1: Limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range;

[0009] Step 2: Based on the limited spatial frequency range, design the corresponding atomic set and the corresponding atomic norm minimization problem, design the corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator;

[0010] Step 3: Solve the SDP problem to obtain an estimate of the spatial frequency;

[0011] Step 4: Convert the obtained spatial frequency estimate into an estimate of DOA.

[0012] As an improvement to the above method, step 1 includes:

[0013] The target DOA range is determined as [θ1,θ2] by the DOA estimation results of the previous set number of frames;

[0014] Map the incident angle θ to the spatial frequency f:

[0015]

[0016] Among them, f L =-sinθ2 / 2,f H =-sinθ1 / 2, representing the lower and upper bounds of the spatial angular frequency respectively.

[0017] As an improvement of the above method, designing the corresponding atom set based on the defined spatial frequency range includes:

[0018] The received signals \(x\) of \(K\) narrowband far - field signal sources arriving at a certain \(N\) - element uniform linear array from the incident angles \(\theta\) k are expressed as:

[0019]

[0020] where \(A(\Theta)=[a(\theta\) 1 ),a(\theta\) 2 ),\(\cdots,a(\theta\) k ),\(\cdots,a(\theta\) K )], \(\theta\) k \(\in[-90^{\circ},90^{\circ}],k = 1,2,\cdots,K\), is the array steering vector, and the superscript \(T\) is the vector transpose; \(f\) c is the carrier frequency; is the element spacing, \(c\) is the wave speed; \(s\) is the amplitude of the signal source; \(n\) is the noise; \(j\) represents the imaginary unit;

[0021] Design the atom set as:

[0022]

[0023] As an improvement of the above method, the atomic norm minimization problem is expressed as:

[0024] When there is noise, the corresponding atomic norm minimization problem is:

[0025]

[0026] where \(\epsilon\) is the noise tolerance; is the optimization variable; is the atomic norm of \(x\); \(\|\cdot\|_2\) represents the 2 - norm;

[0027] When there is no noise, the corresponding atomic norm minimization problem is:

[0028]

[0029] As an improvement of the above method, designing the corresponding mapping operator includes:

[0030] The dual problem of the atomic norm minimization problem is:

[0031]

[0032] where \(q\) is the dual variable; represents the real part of the inner product; Denote the dual atomic norm of \(q\):

[0033]

[0034] where the dual trigonometric polynomial \(q\) k denotes the \(k\)-th element of \(q\);

[0035] Filter \(q\) using an ideal band-pass filter The ideal band-pass filter becomes \(h\) after a delay of \(n\) n ; construct the ideal band-pass filter operator \(H\circ q\) represents the action of the ideal band-pass filter operator on \(q\), and \(q\) H denotes the output after the ideal band-pass filter operator acts on \(q\);

[0036] Obtain a new filtering operator \(H_0\) by phase-shifting \(H\), with the frequency band \([-B / 2, B / 2]\): When \(1 / B\) is an integer, directly downsample it by a factor of \(1 / B\); when \(1 / B\) is not an integer, interpolate it using an analog low-pass filter to obtain an analog signal, and then resample it at the Nyquist sampling rate of \(B\);

[0037] The filtering matrix obtained by truncating \(H_0\) according to the set criterion is The output signal is which is the mapping operator.

[0038] As an improvement of the above method, the set criterion is:

[0039] Take the ideal filter \(h\) with delay n The fourth to fifth side lobes and \(i\) in front of the index \(i_0\) corresponding to the center of the main lobe N-1 The fourth to fifth side lobes at the back.

[0040] As an improvement of the above method, converting the atomic norm minimization problem to an SDP problem based on the mapping operator includes:

[0041] The SDP problem is:

[0042]

[0043] where \(y\) and \(u\) are optimization variables; \(\text{Toep}(u)\) represents the Toeplitz operator; the superscript \(H\) represents the conjugate transpose; \(u_0\) represents the 0-th element of the vector \(u\); \(t > 0\) is an optimization variable;

[0044] When considering noise and truncation error, the SDP problem is:

[0045]

[0046] The objective function is:

[0047]

[0048] where τ is the regularization parameter.

[0049] As an improvement of the above method, solving the SDP problem includes:

[0050] Consider the atomic norm minimization problem:

[0051]

[0052] where M represents the dimension of the positive semi - definite matrix, Y = G H Y0, μ = [u T , vec(X) T , vec(W) T T , X, W represent the intermediate variables for optimization; Y represents a quantity with known value, that is, the received signal containing noise;

[0053] p * is the optimal value of the optimized objective function, defined as:

[0054]

[0055] is the dual cone, defined as

[0056] Solve for X and W using the following formula:

[0057]

[0058] where U is a matrix composed of eigenvectors with non - zero eigenvalues, and each column is an eigenvector; Σ is a diagonal matrix with eigenvalues on the diagonal; I represents the identity matrix; ξ represents the dual variable;

[0059] Solve for u using the following formula:

[0060]

[0061] where ψ t (U) is the objective function, defined as

[0062] Use the limited - memory method in the quasi - Newton method to implement the solution of the above formula, where the Hessian matrix is initialized as:

[0063] ​

[0064] The calculation formula is as follows:

[0065]

[0066] This application also provides a DOA estimation system based on atomic norm minimization, which is implemented based on the above method. The system includes:

[0067] A conversion spatial frequency range module for limiting the DOA range of the incoming wave and converting the limited DOA range into a spatial frequency range;

[0068] A conversion SDP problem module for designing a corresponding atom set and a corresponding atomic norm minimization problem based on the limited spatial frequency range, designing a corresponding mapping operator, and converting the atomic norm minimization problem into an SDP problem based on the mapping operator;

[0069] An SDP problem solving module for solving the SDP problem to obtain an estimate of the spatial frequency;

[0070] A DOA estimation module for converting the obtained spatial frequency estimate into an estimated value of the DOA.

[0071] Compared with the prior art, the advantages of this application are as follows:

[0072] The present invention studies the DOA estimation problem with limited azimuth and proposes a gridless DOA estimation method CA-ANM based on the atomic norm minimization method. In order to improve the efficiency of the CA-ANM algorithm, the present invention further proposes a fast solution algorithm FCA-ANM based on the primal-dual interior point method. The simulation results show that, compared with the original ANM algorithm, the proposed CA-ANM algorithm has a shorter operation time and can achieve more accurate DOA estimation in a noisy environment. Compared with the grid algorithm, the proposed algorithm has higher DOA estimation accuracy because it is not affected by grid mismatch. Compared with the original algorithm CA-ANM, the proposed fast algorithm FCA-ANM has the same performance in terms of accuracy, but the operation time is greatly reduced, reflecting the superiority of the fast algorithm. Description of the Drawings

[0073] Figure 1 Shown is a schematic diagram of a uniform linear array (ULA) receiving a far-field plane wave signal emitted by a target;

[0074] Figure 2 Shown is a flowchart of a DOA estimation method for a limited area based on atomic norm minimization;

[0075] As shown in Fig. 3(a), the number of array elements N mNumber of successful estimations vs. SNR at different SNRs when = 50;

[0076] Figure 3(b) shows the number of array elements N m RMSE at different SNRs when = 50;

[0077] Figure 3(c) shows the number of array elements N m Number of successful estimations vs. SNR at different SNRs when = 100;

[0078] Figure 3(d) shows the number of array elements N m RMSE at different SNRs when = 100;

[0079] Figure 4(a) shows the RMSE and operation time when the number of array elements is different, RMSE vs. the number of array elements;

[0080] Figure 4(b) shows the RMSE and operation time when the number of array elements is different, operation time vs. the number of array elements;

[0081] Figure 5(a) shows the RMSE and operation time when the DOA search range is different, RMSE vs. the angular range;

[0082] Figure 5(b) shows the RMSE and operation time when the DOA search range is different, operation time vs. the angular range. Detailed implementation manners

[0083] The technical solutions of the present application will be described in detail below with reference to the accompanying drawings.

[0084] In order to make full use of the prior knowledge of the target DOA range to reduce the computational complexity, the present invention mainly solves the DOA estimation problem in a limited area, and designs a gridless DOA estimation method based on the atomic norm minimization theory to avoid the influence of grid mismatch. The present invention solves the problem that the dual atomic norm constraint corresponding to the constructed atom set does not satisfy the form of the triangular polynomial inequality constraint of the full frequency band by performing band - pass filtering and down - sampling on the dual variable, converts the ANM problem into a solvable SDP problem, and proposes a limited - area ANM method (CA - ANM). The present invention also proposes a fast method (FCA - ANM) for solving the CA - ANM problem, which improves the operation speed through the primal - dual interior - point method. Theoretical analysis and simulation results show that, compared with the original ANM algorithm, the CA - ANM method provided by the present application has lower computational complexity and shorter operation time; the fast solution method FCA - ANM provided by the present application further improves the operation speed on the condition that the estimation accuracy is consistent with that of the CA - ANM method solved by the CVX toolbox.

[0085] A method and system for DOA estimation in a limited area based on atomic norm minimization provided by this application processes signals that are far-field plane wave signals emitted by a target received by a uniform linear array (ULA). The schematic diagram is as Figure 1 shown.

[0086] Embodiment 1

[0087] As Figure 2 shown, a method for DOA estimation in a limited area based on atomic norm minimization provided by this application first limits the DOA range of the incoming wave, then converts the limited DOA range into a spatial frequency range, and designs a corresponding mapping operator Next, based on the limited spatial frequency range, a corresponding atom set and the corresponding ANM problem are designed, and the ANM problem is converted into an SDP problem based on the mapping operator. Then, the CVX solver package or a fast algorithm based on the primal-dual interior point method is used to solve the above SDP problem, and the estimation of the spatial frequency is obtained. Finally, the obtained spatial frequency estimation is converted into the estimated value of the DOA, realizing the DOA estimation in the limited area.

[0088] Assume that K narrowband far-field signal sources existing in space arrive at a certain N-element uniform linear array from directions θ k respectively. The received signal x can be expressed as

[0089]

[0090] where A(Θ) = [a(θ 1 ), a(θ 2 ), …, a(θ K ), θ k ∈[-90°, 90°], k = 1, 2, …, K, is the array steering vector, f c is the carrier frequency, is the element spacing, c is the wave speed, s is the amplitude of the signal source, and n is the noise. The goal of DOA estimation is to recover the direction from the received signal x. In practical engineering, there is usually a priori knowledge about the range of the target DOA. For example, in the problem of multi-frame DOA estimation, the range of the target DOA can be determined through the DOA estimation results of the previous few frames, and it can be considered that the DOA of the target only exists within the determined range in the estimation of subsequent frames; when interference exists, the incoming wave signals can be made to only exist within the previously demarcated range by using a spatial matrix filter. Assume that the known DOA range of the target is [θ1, θ2], that is Based on this characteristic, the target DOA can be searched only within the range of [θ1, θ2], so the algorithm complexity is relatively lower compared to the case where the DOA range is unknown.

[0091] When the range of the DOA of the target is not restricted, in order to estimate the DOA from the received signal \(x\) in Equation (1), an atom set can be constructed and the atomic norm of \(x\) is defined as

[0092]

[0093] where \(t>0\) is the optimization variable, \(c\) k \(>0\) is the optimization variable, \(a\) k is an atom, denotes that \(a\) k belongs to the set The DOA of the target can be estimated by solving the following atomic norm minimization problem

[0094]

[0095] When there is no noise, the above equation is converted to

[0096]

[0097] where, is the optimization variable.

[0098] Since the optimization problems (3) and (4) involve the optimization of an infinite number of parameters, they cannot be solved directly. Candès et al. pointed out that the optimization problem (4) is equivalent to the following semidefinite programming problem

[0099]

[0100] where \(u_0\) represents the 0th element of the vector \(u\), and Toep(·) is the Toeplitz operator, whose role is to convert the vector \(u\) into the Toeplitz matrix Toep(\(u\)). It is easy to prove that by performing a Vandermonde decomposition on the obtained Toep(\(u\)), we get

[0101]

[0102] where \(p\) k \(>0\), is a real number greater than 0. The obtained is the target DOA.

[0103] Next, consider the case where there is a prior for the range of the target DOA. In the grid-based parameter estimation framework based on a dictionary, this prior information can be realized by restricting the dictionary, that is, making the discrete partition correspond to the gridless framework, that is, constructing a new atom set The corresponding atomic norm minimization problem can be proposed as

[0104]

[0105] where ∈ is the noise tolerance. When there is no noise, the corresponding atomic norm minimization problem is

[0106]

[0107] It is easy to obtain that the dual problem of (8) is

[0108]

[0109] where q is the dual variable, denotes the dual atomic norm of q and can be calculated as

[0110]

[0111] where denotes the real part of the inner product, and <·,·> denotes the inner product. For the convenience of subsequent discussion, here the incident angle θ is mapped to the spatial frequency f, that is where f L =-sinθ2 / 2, f H =-sinθ1 / 2, representing the lower and upper bounds of the spatial angular frequency respectively. Substituting them into (10) gives

[0112]

[0113] where the dual trigonometric polynomial where q k represents the k-th element of q.

[0114] Therefore is equivalent to

[0115]

[0116] Note that by using the bounded lemma, can be converted into a semi-definite constraint. However, only the case where f ∈ [f L , f H is restricted in (12), so the bounded lemma is not applicable. Note that the dual polynomial q(f) can actually be regarded as the discrete-time Fourier transform of a finite-support sequence q. Therefore, the bounded lemma can be used by processing q. An intuitive processing method is to filter the sequence q using a finite impulse response (FIR) band-pass filter to convert it into a band-pass signal with a passband of [f L , f H . At this time and Approximately equivalent. However, this method has very poor robustness and it is difficult to achieve accurate DOA estimation when the signal-to-noise ratio is low. In the present invention, first, an ideal band-pass filter is used to filter q. Let the ideal low-pass filter be After delaying it by n, it becomes h n , and an ideal band-pass filter operator can be constructed as where H represents the ideal band-pass filter operator, H°q represents applying the ideal band-pass filter operator to q, and q H represents the output after applying the ideal band-pass filter operator to q. At this time, At this time, since the output q H has an infinite dimension, it is impossible to convert the problem into an SDP problem of finite size. Suppose we hope to construct an M×M positive semi-definite matrix using the bounded lemma, then a dimension reduction operator needs to be designed such that

[0117]

[0118] where is the discrete-time Fourier transform operator. Simulation experiments prove that it is sufficient to design G by simply truncating q H . Specifically, suppose the index corresponding to the center of the main lobe of the delayed ideal filter h n is i n , then take the part of the fourth to fifth side lobes on the front side of i0 and the same fourth to fifth side lobes on the back side of i N-1 . Let the truncated h n be The obtained output is At this time, does not satisfy the conditions of LTI, but a better effect can be obtained than using an LTI FIR band-pass filter.

[0119] Although truncation is performed, the output by the above method usually has a higher dimension than q, resulting in a larger scale of the finally obtained SDP problem. In fact, since 's spectrum is restricted in [f L , f H , and its bandwidth B = f H - f L < 1, it is entirely possible to reduce the number of samples by reducing the sampling rate. Specifically, going back to q H before truncation, this signal is a strict band-pass signal, so it can be resampled to fill its frequency band . Specifically, first, phase shift q H to make it a low-pass signal The frequency band is [-B / 2, B / 2]. When 1 / B is an integer, it can be directly decimated by a factor of 1 / B; when 1 / B is not an integer, it can be interpolated with an analog low-pass filter to obtain an analog signal, and then resampled at the Nyquist sampling rate of B. The above process can be directly applied to H, let the new filtering operator obtained be H0, and let the filtering matrix obtained by truncating with a similar criterion as before be The output signal is At this time The dimension N0 of is usually smaller than q, approximately O(NB). Based on the above discussion, can be Approximately equivalent to From this, the SDP form of the dual problem can be obtained as

[0120]

[0121] At this time, the SDP of the corresponding primal problem is

[0122]

[0123] where y and u are optimization variables. By performing Vandermonde decomposition on the solved Toep(u), it can be written in the following form:

[0124]

[0125] where a 0 (f 0 ) = [1, exp(j2πf 0 ), …, exp(j2πf 0 (N0 - 1))] T . At this time, the frequency can be extracted from it According to the discussion above, the spatial frequency f can be obtained through and mapping the spatial frequency back to the incident angle θ can achieve DOA estimation.

[0126] When considering noise and truncation error, (15) can be modified to

[0127]

[0128] where ∈ is the noise tolerance.

[0129] In summary, the DOA estimation problem of the above restricted region problem can be solved by (17). The constraint on noise in (17) is a hard constraint, so there may be a problem of no feasible solution in some cases. For this reason, the objective function can be rewritten in the form of atomic norm regularization

[0130]

[0131] where τ is the regularization parameter. To efficiently solve (18), the present invention also provides a fast solution method based on the primal-dual interior point method.

[0132] The fast solution method proposed by the present invention mainly considers the following more general atomic norm minimization problem:

[0133]

[0134] where Y = G H Y0, μ = [u T , vec(X) T , vec(W) T [[ID=2f1]]] T , where p * is the optimal value of the optimized objective function, defined as is the dual cone, defined as Denote the dual variable of μ as λ = [ξ T , vec(S) T , vec(V) T T , where ξ and S represent dual variables, and X and W represent intermediate variables in the optimization. Y represents a quantity with a known value in the fast solution algorithm and can be understood as the received signal containing noise in the context described above. The Lagrange dual of the optimization problem (19) can be written in the following form:

[0135]

[0136] where is 's dual cone. Gao et al. derived

[0137]

[0138] From this, the corresponding dual problem can be obtained as

[0139]

[0140] The dual gap is

[0141] η(μ, λ) = f(μ) - g(λ). (23)

[0142] Since the original problem (19) satisfies the KKT condition, thus

[0143] η(μ * , λ * ) = p​* -d * = 0. (24)

[0144] It is easy to prove that the value of the dual function g(λ) is

[0145]

[0146] where is the eigenvalue decomposition of GG H , where U is a matrix composed of eigenvectors with non-zero eigenvalues, each column of which is an eigenvector, and Σ is a diagonal matrix with eigenvalues on the diagonal.

[0147] The objective of the present invention is to solve the primal problem (19) within the framework of the primal-dual interior point method. Specifically, the present invention uses the logarithmically homogeneous function F(μ) as the barrier function corresponding to the primal cone . F(μ) is defined as follows:

[0148]

[0149] The primal-dual interior point method solves (19) by iteratively solving the following problem

[0150]

[0151] where t > 0 increases in each iteration. It can be proved that when t → ∞, the solution of (27) converges to the solution of (19). The augmented KKT conditions corresponding to (27) are

[0152]

[0153] Here the gradient operator is defined as the Wirtinger derivative, which is defined as

[0154] From (28), we have

[0155]

[0156] and

[0157]

[0158] where Λ = Toep -1 (u)X. The primal-dual variables μ and λ can be solved from (30) and (29) respectively. Specifically, u can be obtained from (30)(c), X can be solved from (30)(b), and W can be obtained through (30)(a). The dual variable can be obtained from (29). In the problem studied by the present invention, the main technical difficulty comes from the solution of (30)(b). The present invention proposes the following method to efficiently solve (30)(b). Based on GGH The eigenvalue decomposition of

[0159] -τΛ=UΣU H (Toep(u)Λ - Y0). (31)

[0160] Multiplying both sides of (31) on the left by UU H , we can obtain

[0161] -τUU H Λ=UΣU H (Toep(u)Λ - Y0)= -τΛ. (32)

[0162] From this, we can get

[0163] UU H Λ=Λ. (33)

[0164] From (31), we can get

[0165] -τUΣ -1 U H Λ=UU H (Toep(u)Λ - Y0). (34)

[0166] Therefore, we have

[0167] UU H (Toep(u)+τUΣ -1 U H )Λ=UU H Y0. (35)

[0168] Noting that UU H Λ=Λ, so we have

[0169] UU H (Toep(u)+τUΣ -1 U H )UU H Λ=UU H Y0. (36)

[0170] Since U is a column full - rank matrix, so we have

[0171]

[0172] Combining with (33), we can get the following result:

[0173] Λ=U(U H Toep(u)U+τΣ -1 ) -1 U H Y0. (38)

[0174] It is worth noting that in (38), only the inverse of (U H Toep(u)U + τΣ -1 ) needs to be calculated, and its scale is significantly reduced compared to the matrix whose inverse needs to be calculated in (30)(b). It should be noted that usually only the observation Y = G H Y0 of Y0 can be obtained. To calculate U H Y0 that appears in (38), it can be deduced that

[0175] U H Y0 = Σ -1 U H GY. (39)

[0176] Based on the above discussion, when u is obtained, the original variables and dual variables can be analytically solved and written as functions of u, which are μ(u) and λ(u) respectively. The above results are summarized as follows:

[0177]

[0178] It should be noted that (40)(1), (40)(P1), and (40)(P2) can be solved efficiently. Specifically, Toep(u)U and Toep(u)Λ in (40)(1) and (40)(P1) can be solved by the fast Fourier transform, while Toep -1 (u)X in (40)(P2) can be solved quickly by the Levinson - Durbin algorithm.

[0179] Finally, consider the problem of solving u from (30)(c). It is easy to verify that the solution of (30)(c) is equivalent to the following optimization problem:

[0180]

[0181] where ψ t (U) is the objective function, defined as

[0182] Problem (41) can be solved by the Newton or quasi - Newton method. The present invention uses the Limited Memory BFGS (L - BFGS) algorithm in the quasi - Newton method to efficiently solve the search direction. The L - BFGS algorithm needs to initialize the Hessian matrix, which is usually initialized as the identity matrix Hansen et al. found that such an initialization would lead to slow convergence of the algorithm, so they proposed a heuristic diagonal matrix approximation of the Hessian matrix

[0183]

[0184] In the present invention, this initialization scheme is adopted, and the is calculated as follows:

[0185]

[0186] So far, the derivation of the main update steps of the algorithm has been completed. Based on the framework of the primal-dual interior point method, the complete algorithm is summarized in the algorithm in Table 1.

[0187] Table 1 Limited Region Atomic Norm Minimization Method Based on Primal-Dual Interior Point Method (FCA-ANM)

[0188]

[0189] Embodiment 2

[0190] This application also provides a limited region DOA estimation system based on atomic norm minimization, which is implemented based on the above method. The system includes:

[0191] A conversion space frequency range module for limiting the DOA range of the incoming wave and converting the limited DOA range into a space frequency range;

[0192] A conversion SDP problem module for designing a corresponding atom set and a corresponding atomic norm minimization problem based on the limited space frequency range, designing a corresponding mapping operator, and converting the atomic norm minimization problem into an SDP problem based on the mapping operator;

[0193] An SDP problem solving module for solving the SDP problem to obtain an estimate of the space frequency;

[0194] A DOA estimation module for converting the obtained space frequency estimate into an estimated value of the DOA.

[0195] Simulation experiment:

[0196] Suppose there are 3 far-field plane wave signals illuminating a uniform linear array with the number of array elements N m whose DOAs are 0°, 13.2° and -17.8° respectively. The above DOA settings are mainly to reflect the influence of grid mismatch of the grid method. Specifically, Θ is obtained by dividing at intervals of 1° dic , and at this time, there are DOAs of the target that do not fall on the grid, resulting in grid mismatch. Suppose the range of the target DOA is known to be [-30°, 30°] before detection. Let the received noise be circular Gaussian white noise obeying The element spacing is taken as half wavelength There is no aliasing at this time. This simulation is all based on MATLAB 2021b, and the optimization problems involved are all solved using the CVX toolbox, with the solver being SDPT3. The CPU of the device used in the simulation is 12th Gen Intel(R) Core(TM) i7-12700H 2.30GHz, and the onboard RAM is 16.0GB.

[0197] The following algorithms are mainly used in this simulation for DOA estimation: First, ignoring the prior of the DOA range, directly using the atomic norm minimization method described in (5) for DOA estimation, and excluding the DOA estimation values outside the defined region after estimation. This method is called the ANM method; Second, the grid method based on the dictionary constructed by discrete partitioning Θ dic To improve efficiency, the SBL method is used for solution, and this method is called the Constrained Area-SBL method, abbreviated as the CA-SBL method; Then, it is solved based on the atomic norm minimization method described in (17) proposed in the present invention, named the Constrained Area-ANM method, abbreviated as the CA-ANM method; Finally, the method of using the fast solution method based on the primal-dual interior point method proposed in the present invention to solve the CA-ANM is called the Fast CA-ANM method, abbreviated as the FCA-ANM method.

[0198] First, verify the performance of each algorithm in a noisy environment. Verify the performance of the algorithms when the number of array elements N m = 50 and N m = 100 respectively, and set the signal-to-noise ratio (SNR) to -20dB to 30dB and -20dB to 10dB. Conduct 100 Monte Carlo simulations for each SNR. The performance of the algorithm can be characterized from two aspects: the number of successful estimations and the root mean square error (RMSE) of the estimation. For a single Monte Carlo simulation, define the root square error (RSE) of the estimation as

[0199]

[0200] where N Tar = 3 is the number of targets, is the DOA estimation of the k-th target, and θ k is the true DOA of the k-th target. Define that the estimation is successful when RSE ≤ 5°. For the simulation rounds with successful estimations, define the RMSE of the estimation as

[0201]

[0202] wherein is the number of successful estimations, is the DOA estimation of the k-th target in the n-th successful estimation. When N m = 50, the curves of the number of successful estimations and RMSE of each algorithm versus SNR are plotted in Figs. 3(a) and 3(b) respectively. When N m = 100, the curves of the number of successful estimations and RMSE of each algorithm versus SNR are plotted in Figs. 3(c) and 3(d) respectively.

[0203] When N m = 50, it can be seen from Fig. 3(a) that when the signal-to-noise ratio is higher than 0 dB, all methods can successfully achieve DOA estimation. Among them, the proposed estimation methods with limited regions, CA-SBL, CA-ANM, and FCA-ANM, of the present invention have a slightly higher success rate than the original ANM method at low signal-to-noise ratios. This shows that limiting the DOA estimation range is beneficial to improving the robustness of the algorithm. It can be seen from Fig. 3(b) that the estimation accuracies of all methods are comparable at low signal-to-noise ratios. However, when the signal-to-noise ratio is higher than 0 dB, due to the influence of grid mismatch, the performance of the CA-SBL method no longer improves significantly, while the other meshless methods improve the estimation accuracy as the signal-to-noise ratio increases, demonstrating the superiority of meshless methods. Among the meshless methods, the performances of ANM, CA-ANM, and FCA-ANM are almost the same. When N m = 100, it can be seen from Figs. 3(c) and 3(d) that the trends of the number of successful estimations and RMSE of each algorithm versus SNR at this time are the same as those when N m = 50. But at the same SNR, when N m = 100, the number of successful estimations and accuracies of each algorithm are higher than those when N m = 50. Thus, it can be seen that increasing the number of array elements is beneficial to improving the robustness of the algorithm.

[0204] Next, the operation speeds of each algorithm are compared, and at this time, the SNR is fixed at 30 dB. First, the influence of the number of array elements on the operation speed is simulated. The number of array elements N m is set to range from 20 to 100. For each N m , 100 Monte Carlo simulations are performed and the operation times of each algorithm are recorded. The times recorded in each simulation are averaged as the actual operation time. It is still considered that the estimation is successful when the RSE of a single simulation ≤ 5°. The simulation results show that for each N m , the number of successful estimations of all algorithms is 100 times at this time. The curves of RMSE and operation time of each algorithm versus the number of array elements N m are plotted in Figs. 4(a) and 4(b) respectively.

[0205] As can be seen from Fig. 4(a), due to the influence of grid mismatch, the estimation accuracy of the grid-based CA-SBL algorithm is weaker than that of the gridless ANM, CA-ANM, and FCA-ANM methods. Among the gridless algorithms, the accuracy of the ANM algorithm for DOA estimation over the entire region is slightly better than that of the proposed DOA estimation algorithm for the limited region, while the fast algorithms FCA-ANM and CA-ANM perform basically the same, further demonstrating the reliability of the fast algorithms. Since the solution mechanism of the grid-based CA-SBL algorithm is different from that of the other gridless algorithms, for fairness, only the operation times of the three gridless algorithms are compared in Fig. 4(b). As can be seen from Fig. 4(b), as the number of array elements increases, the operation times of all algorithms show an upward trend. At the same time, for all numbers of array elements, the operation time of the fast algorithm FCA-ANM proposed in the present invention is significantly lower than that of other algorithms, demonstrating the high efficiency of the proposed fast algorithm. In addition, it can also be found that when the number of array elements is large, the operation time of the CA-ANM algorithm proposed in the present invention is significantly lower than that of the ANM algorithm. This is consistent with the previous discussion, that is, the dimension of the dual variable q in the ANM algorithm is N, while the dimension of the dual variable is approximately O(NB), where B < 1 is the bandwidth of the spatial frequency. Therefore, the computational complexity of the CA-ANM algorithm is less than that of the ANM algorithm.

[0206] Next, the influence of the DOA search range on the algorithm speed is simulated. Still fixing SNR = 30 dB, with 0° as the center, the angular range is set from 60° (i.e., [-30°, 30°]) to 180° (i.e., the full space), and 100 Monte Carlo simulations are performed. The simulation results show that for each search range, the number of successful estimations of all algorithms is 100 times at this time. The curves of the RMSE and operation time of each algorithm varying with the angular range are plotted in Figs. 5(a) and 5(b) respectively.

[0207] As can be seen from Fig. 5(a), due to the influence of grid mismatch, the estimation accuracy of the CA-SBL algorithm with grids is weaker than that of the gridless ANM, CA-ANM, and FCA-ANM methods. For different DOA search ranges, the accuracy of each algorithm is basically the same, indicating that restricting the DOA search range does not affect the accuracy of the algorithm. As can be seen from Fig. 5(b), the operation times of the proposed restricted region algorithms CA-ANM and FCA-ANM both increase with the expansion of the angle search range. When the DOA search range does not exceed 120°, the operation time of the proposed CA-ANM algorithm is less than that of the ANM algorithm, indicating that the operation speed can be improved by reducing the DOA search range, which is consistent with the previous analysis. When the DOA search range is large, at this time the bandwidth B of the spatial frequency is close to 1, and the reduction of the algorithm scale by the CA-ANM algorithm is not obvious. At the same time, the form of the corresponding optimization problem is more complex than that of the ANM algorithm, resulting in its operation time exceeding that of the ANM algorithm. In contrast, the operation speed of the proposed fast algorithm FCA-ANM is significantly better than that of the CA-ANM algorithm, and the operation time is always less than that of the ANM algorithm, demonstrating the efficiency of the proposed fast algorithm.

[0208] The present application can also provide a computer device, including: at least one processor, a memory, at least one network interface, and a user interface. Each component in the device is coupled together through a bus system. It can be understood that the bus system is used to realize the connection and communication between these components. In addition to the data bus, the bus system also includes a power bus, a control bus, and a status signal bus.

[0209] Among them, the user interface can include a display, a keyboard, or a pointing device. For example, a mouse, a trackball, a touchpad, or a touch screen, etc.

[0210] It can be understood that the memory in the disclosed embodiments of the present application can be a volatile memory or a non-volatile memory, or can include both volatile and non-volatile memories. Among them, the non-volatile memory can be a read-only memory (ROM), a programmable read-only memory (PROM), an erasable programmable read-only memory (EPROM), an electrically erasable programmable read-only memory (EEPROM), or a flash memory. The volatile memory can be a random access memory (RAM), which is used as an external cache. By way of example but not limitation, many forms of RAM are available, such as static random access memory (SRAM), dynamic random access memory (DRAM), synchronous dynamic random access memory (SDRAM), double data rate synchronous dynamic random access memory (DDR SDRAM), enhanced synchronous dynamic random access memory (ESDRAM), synchlink dynamic random access memory (SLDRAM), and direct rambus random access memory (DRRAM). The memories described herein are intended to include but not be limited to these and any other suitable types of memories.

[0211] In some embodiments, the memory stores the following elements, executable modules, or data structures, or subsets thereof, or extended sets thereof: an operating system and application programs.

[0212] Among them, the operating system includes various system programs, such as a framework layer, a core library layer, a driver layer, etc., and is used to implement various basic services and handle hardware-based tasks. The application programs include various application programs, such as a media player and a browser, etc., and are used to implement various application services. The program for implementing the method of the disclosed embodiments of the present application can be included in the application programs.

[0213] In the above-mentioned embodiments, by calling the programs or instructions stored in the memory, specifically, the programs or instructions stored in the application programs, the processor is configured to:

[0214] Execute the steps of the above method.

[0215] The above method can be applied to a processor or implemented by a processor. The processor may be an integrated circuit chip with the ability to process signals. During implementation, each step of the above method can be completed by the integrated logic circuit of the hardware in the processor or instructions in the form of software. The above-mentioned processor may be a general-purpose processor, a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components. It can implement or execute the various methods, steps, and logic block diagrams disclosed above. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc. Combining the steps of the above-disclosed method can be directly embodied as being executed and completed by a hardware decoding processor, or executed and completed by a combination of the hardware and software modules in the decoding processor. The software module may be located in a mature storage medium in the art such as a random access memory, a flash memory, a read-only memory, a programmable read-only memory, or an electrically erasable programmable memory, a register, etc. This storage medium is located in the memory, and the processor reads the information in the memory and combines its hardware to complete the steps of the above method.

[0216] It can be understood that these embodiments described in the present application can be implemented using hardware, software, firmware, middleware, microcode, or a combination thereof. For hardware implementation, the processing unit can be implemented in one or more application specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field-programmable gate arrays (FPGAs), general-purpose processors, controllers, microcontrollers, microprocessors, other electronic units for performing the functions described in the present application, or a combination thereof.

[0217] For software implementation, the technology of the present application can be implemented by executing the functional modules of the present application (such as procedures, functions, etc.). The software code can be stored in the memory and executed by the processor. The memory can be implemented inside or outside the processor. [[ID=!]]

[0218] The present application may also provide a non-volatile storage medium for storing a computer program. When the computer program is executed by a processor, the various steps in the above method embodiments can be implemented.

[0219] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application and not to limit them. Although the present application has been described in detail with reference to the embodiments, those of ordinary skill in the art should understand that any modification or equivalent replacement of the technical solutions of the present application does not depart from the spirit and scope of the technical solutions of the present application, and they should all be covered within the scope of the claims of the present application.

Claims

1. A method for estimating DOA in a limited area based on atomic norm minimization, comprising: Step 1: Limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range; Step 2: Design a corresponding atom set and the corresponding atomic norm minimization problem based on the limited spatial frequency range, design a corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator; Step 3: Solve the SDP problem to obtain an estimate of the spatial frequency; Step 4: Convert the obtained spatial frequency estimate into an estimate of the DOA.

2. The DOA estimation method for a limited area based on atomic norm minimization according to claim 1, wherein The said Step 1 includes: Determine the range of the target DOA as [θ1, θ2] through the DOA estimation results of the previous set number of frames; Map the incident angle θ to the spatial frequency f: where f L = -sinθ2 / 2, f H = -sinθ1 / 2, represent the lower and upper bounds of the spatial angular frequency respectively.

3. The method for DOA estimation in a limited area based on atomic norm minimization according to claim 2, wherein The said design of a corresponding atom set based on the limited spatial frequency range includes: K narrowband far-field signal sources respectively arrive from the incident angle θ k The received signal x of a certain N-element uniform linear array is expressed as: where \(A(\Theta)=[a(\theta 1 ),a(\theta 2 ),\cdots,a(\theta k ),\cdots,a(\theta K )]\), \(\theta k \in[-90^{\circ},90^{\circ}], k = 1,2,\cdots,K\), is the array steering vector, and the superscript \(T\) represents the vector transpose; \(f c is the carrier frequency; is the element spacing, \(c\) is the wave speed; \(s\) is the amplitude of the signal source; \(n\) is the noise; \(j\) represents the imaginary unit; Set of design atoms is as follows:

4. The method for estimating the DOA of a limited area based on atomic norm minimization according to claim 3, characterized in that, The atomic norm minimization problem is expressed as: When there is noise, the corresponding atomic norm minimization problem is: where ∈ is the noise tolerance; is the optimization variable; is the atomic norm of x; ||·||2 represents the 2-norm; When there is no noise, the corresponding atomic norm minimization problem is:

5. The method for estimating the DOA of a limited area based on atomic norm minimization according to claim 4, wherein The said design of a corresponding mapping operator includes: The dual problem of the atomic norm minimization problem is: where \(q\) is the dual variable; denotes the real part of the inner product;; denotes the dual atomic norm of \(q\): Among them, the dual trigonometric polynomial q k represents the k-th element of q; Using an ideal band-pass filter Filter q. The ideal band-pass filter has a delay of n and is denoted as h n ; Construct the ideal band-pass filter operator H°q means applying the ideal band-pass filter operator to q. q H represents the output after applying the ideal band-pass filter operator to q; A new filtering operator H0 obtained by phase-shifting H, with the frequency band [-B / 2, B / 2]: When 1 / B is an integer, directly downsample it by 1 / B times; when 1 / B is not an integer, interpolate it with an analog low-pass filter to obtain an analog signal, and then resample it at the Nyquist sampling rate of B; The filtering matrix obtained by truncating H0 according to the set criterion is The output signal is which is the mapping operator.

6. The DOA estimation method for a limited area based on atomic norm minimization according to claim 5, wherein The said setting criterion is: Take the delay ideal filter h n The fourth to fifth side lobes in front of the index i0 corresponding to the center of the main lobe and i N-1 The fourth to fifth side lobes at the back.

7. The DOA estimation method for a limited area based on atomic norm minimization according to claim 5, wherein The said conversion of the atomic norm minimization problem into an SDP problem based on the mapping operator includes: The SDP problem is: Wherein, y and u are optimization variables; Toep(u) represents the Toeplitz operator; the superscript H represents the conjugate transpose; u0 represents the 0th element of the vector u; t > 0 is an optimization variable; When considering noise and truncation error, the SDP problem is: The objective function is: Wherein, τ is the regularization parameter.

8. The method for DOA estimation in a limited area based on atomic norm minimization according to claim 7, wherein The said solution of the SDP problem includes: Consider the atomic norm minimization problem: Among them, M represents the dimension of the positive semi - definite matrix, Y = G H Y0, μ = [u T , vec(X) T , vec(W) T T , X, W represent the intermediate variables for optimization; Y represents a quantity with a known value, that is, the received signal containing noise;​ p * The optimal value of the optimized objective function, defined as: is the dual cone, defined as Solve to obtain X and W using the following formula: Wherein, U is a matrix composed of eigenvectors with non-zero eigenvalues, and each column is an eigenvector; Σ is a diagonal matrix with eigenvalues on the diagonal; I represents the identity matrix; ξ represents the dual variable; Solve to obtain u using the following formula: Among them, ψ t (U) is the objective function, defined as The above formula is solved using the limited memory method in the quasi-Newton method, where the Hessian matrix is initialized as: The calculation formula is as follows:

9. A DOA estimation system in a limited area based on atomic norm minimization, implemented based on any one of the methods described in claims 1-8, characterized in that The said system includes: A conversion spatial frequency range module, used to limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range; A conversion SDP problem module, used to design a corresponding atom set and the corresponding atomic norm minimization problem based on the limited spatial frequency range, design a corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator; An SDP problem solution module, used to solve the SDP problem to obtain an estimate of the spatial frequency; and A DOA estimation module, used to convert the obtained spatial frequency estimate into an estimate of the DOA.

Images