A method and system for DOA estimation in a limited area based on atomic norm minimization

Abstract

This application provides a method and system for estimating DOA in a limited area based on atomic norm minimization. The method includes: limiting the DOA range of the incoming wave and converting the limited DOA range into a spatial frequency range; designing a corresponding set of atoms 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; solving the SDP problem to obtain an estimate of the spatial frequency; and converting the obtained spatial frequency estimate into an estimated DOA value. The advantages of this application are: compared with grid-based algorithms, DOA estimation is more accurate because it is not affected by grid mismatch; compared with the original CA-ANM algorithm, the method of this application is consistent with its performance in terms of accuracy, but with significantly reduced computation time.

Description

Technical Field

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

[0002] Estimating the direction of arrival (DOA) of underwater targets is a hot topic in underwater acoustic array signal processing. To overcome the Rayleigh limit and achieve super-resolution DOA estimation, subspace methods, such as multiple signal classification (MUSIC) and estimation of signal via rotational invariance techniques (ESPRIT), have been widely applied to DOA estimation. Although subspace methods significantly improve the resolution of DOA estimation, they require an accurate estimation of the signal covariance matrix, making them inapplicable to scenarios with a limited number of snapshots. In contrast, algorithms based on compressed sensing have become a research hotspot due to their ability to achieve high-resolution DOA estimation with a small number of snapshots. Early compressed sensing algorithms were based on grid dictionaries, where the target's DOA is divided into a discrete grid. The steering vectors corresponding to each grid point form an overcomplete dictionary. Since the target's DOA falls on a limited grid, it is sparse, which transforms the DOA estimation problem into a sparse recovery problem. In actual engineering, before performing DOA estimation, a certain prior knowledge of the DOA range of the target is usually available, so DOA estimation only needs to be performed within a limited range. In fact, when there is prior information about the target, the computational complexity can be reduced by limiting the DOA search range to the target area. This situation is referred to as DOA estimation in a limited area in the present invention. Yang et al. designed a null-stuck matrix filter so that the DOA of the incoming signal is limited to the target area, and proposed a DOA estimation method based on sparse spectrum fitting. On this basis, Zhang et al. achieved the limitation of the DOA search range by limiting the dictionary range to the target area, and implemented DOA estimation based on the sparse Bayesian learning (SBL) algorithm, thereby improving computational efficiency.

[0003] The aforementioned sparse methods assume that the target DOAs all fall on a predefined grid. However, the actual target DOAs may fall outside the grid, resulting in a grid mismatch problem. Subsequent research has proposed a series of compensatory methods to mitigate the impact of grid mismatch, collectively referred to as off-grid algorithms. However, these methods still assume the existence of a grid, resulting in inaccurate estimation of continuously valued DOAs. Gridless algorithms, such as atomic norm minimization (ANM), have garnered widespread attention in recent years due to their ability to obtain sparse solutions on a continuous dictionary. Because the ANM problem involves optimizing an infinite number of parameters, it cannot be solved directly. To this end, Candès et al. derived the dual problem of ANM and converted the dual atomic norm constraints into trigonometric polynomial inequality constraints. Using the bounded lemma, the trigonometric polynomial inequality constraints can be converted into semidefinite constraints, thereby transforming the dual problem into a semidefinite programming (SDP) problem, which can 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. In order to limit the DOA search range to the target area to reduce the complexity of the algorithm, analogous to the method of limiting the discrete dictionary to the target area in the grid algorithm, the DOA search range can be limited by limiting the constructed atomic set to the target area. However, this will cause the corresponding dual atomic norm constraint to not satisfy the form of the full-band trigonometric polynomial inequality constraint, resulting in the bounded lemma being unusable. This means that 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 scope of the ANM algorithm in recent years, these algorithms are only applicable to corresponding specific scenarios and are not yet universal, and therefore cannot be applied to the limited area DOA estimation problem considered in the present invention.

[0005] In addition, although the SDP problem obtained by 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 to the above method, the designing of a corresponding atom set based on a limited spatial frequency range includes:

[0018] K narrowband far-field signal sources are respectively emitted from the k The received signal x arriving at a certain N-element uniform linear array is expressed as:

[0019]

[0020] Where A(θ)=[a(θ 1 ),a(θ 2 ),…,a(θ k ),…,a(θ K )],θ k ∈[-90°,90°],k=1,2,…,K, is the array steering vector, the superscript T is the vector transpose; f c is the carrier frequency; is the array element spacing, c is the wave velocity; s is the amplitude of the signal source; n is the noise; j represents the complex unit;

[0021] Design Atom Set for:

[0022]

[0023] As an improvement to 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] Among them, ∈ is the noise tolerance; is the optimization variable; is the atomic norm of x; ||·||2 represents the 2-norm;

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

[0028]

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

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

[0031]

[0032] Among them, q is the dual variable; represents the real part of the inner product; Represents the dual atomic norm of q:

[0033]

[0034] Among them, the dual trigonometric polynomial q k represents the kth element of q;

[0035] Using an ideal bandpass filter Filter q, the ideal bandpass filter is h after delay n n ; Construct an ideal bandpass filter operator H°q represents the ideal bandpass filter operator acting on q, q H represents the output of the ideal bandpass filter operator after acting on q;

[0036] The new filter operator H0 obtained by phase shifting H has a frequency band of [-B / 2, B / 2]: when 1 / B is an integer, it is directly downsampled by 1 / B; when 1 / B is not an integer, it is interpolated with an analog low-pass filter to obtain an analog signal, and then resampled according to the Nyquist sampling rate of B;

[0037] The filter matrix obtained by truncating H0 according to the set criteria is The output signal is It is the mapping operator.

[0038] As an improvement to the above method, the setting criteria are:

[0039] Take the delayed ideal filter h n The fourth to fifth side lobes in front of the main lobe center corresponding to index i0 and i N-1 The fourth to fifth posterior side lobes.

[0040] As an improvement to the above method, the atomic norm minimization problem is converted into an SDP problem based on a mapping operator, including:

[0041] The SDP questions are:

[0042]

[0043] Where 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 the optimization variable;

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

[0045]

[0046] The objective function is:

[0047]

[0048] Among them, τ is the regularization parameter.

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

[0050] Consider the atomic norm minimization problem:

[0051]

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

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

[0054]

[0055] is the dual cone, defined as

[0056] Use the following formula to solve for X and W:

[0057]

[0058] Where U is a matrix composed of eigenvectors whose eigenvalues ​​are not 0, 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] Use the following formula to solve for u:

[0060]

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

[0062] The limited memory method in the pseudo-Newton method is used to solve the above equation, where the Hessian matrix Initialized as:

[0063]

[0064] The calculation formula is:

[0065]

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

[0067] The spatial frequency range conversion module is used to limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range;

[0068] The SDP problem conversion module is used to design the corresponding atomic set and the corresponding atomic norm minimization problem based on the limited spatial frequency range, design the corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator;

[0069] The SDP problem solving module is used to solve the SDP problem and obtain the estimation of spatial frequency;

[0070] The DOA estimation module is used to convert the obtained spatial frequency estimation into an estimated value of DOA.

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

[0072] The present invention studies the problem of DOA estimation with limited orientation, 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. 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 gridded 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 is consistent with its performance in terms of accuracy, but the operation time is greatly reduced, which reflects the superiority of the fast algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0074] Figure 2 The figure shows a flow chart of the DOA estimation method for a limited area based on atomic norm minimization;

[0075] Figure 3(a) shows the number of array elements N m=50 for the number of successful estimations under different signal-to-noise ratios VSSNR;

[0076] Figure 3(b) shows the number of array elements N m = 50 under different signal-to-noise ratios;

[0077] Figure 3(c) shows the number of array elements N m =100 for the number of successful estimations under different signal-to-noise ratios VSSNR;

[0078] Figure 3(d) shows the number of array elements N m =100 under different signal-to-noise ratios;

[0079] Figure 4(a) shows the RMSE and computation time for different numbers of array elements, RMSEV.S.Number of array elements;

[0080] Figure 4(b) shows the RMSE and computation time for different numbers of array elements, and the computation time vs. the number of array elements.

[0081] Figure 5(a) shows the RMSE and operation time for different DOA search ranges,RMSEv.s.angle range;

[0082] Figure 5(b) shows the RMSE and computation time for different DOA search ranges, and computation time vs. angle range. DETAILED DESCRIPTION

[0083] The technical solution of this application is described in detail below with reference to the accompanying drawings.

[0084] In order to make full use of the prior 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 atomic set does not satisfy the form of trigonometric polynomial inequality constraints of the full band by bandpass filtering and downsampling the dual variables, converts the ANM problem into a programmable SDP problem, and proposes a limited area ANM method (CA-ANM). The present invention also proposes a fast method (FCA-ANM) that can be used to solve the CA-ANM problem, which improves the computational speed through the original dual interior point method. Theoretical analysis and simulation results show that compared with the original ANM algorithm, the CA-ANM method provided by this application has lower computational complexity and shorter computation time; the fast solution method FCA-ANM provided by this application further improves the computational speed under the condition that the estimation accuracy is consistent with the CA-ANM method solved based on the CVX toolkit.

[0085] The present application provides a method and system for estimating DOA in a limited area based on atomic norm minimization. The signal processed is a far-field plane wave signal emitted by a target and received by a uniform linear array (ULA). The schematic diagram is shown in FIG. Figure 1 shown.

[0086] Example 1

[0087] like Figure 2 As shown in the figure, the present application provides a limited area DOA estimation method based on atomic norm minimization. First, the DOA range of the incoming wave is limited, and then the limited DOA range is converted into a spatial frequency range, and the corresponding mapping operator is designed. Next, we design the corresponding atomic set and the corresponding ANM problem based on the limited spatial frequency range. Using a mapping operator, we convert the ANM problem into an SDP problem. We then use the CVX solver package or a fast algorithm based on the primal-dual interior point method to solve the SDP problem and obtain a spatial frequency estimate. Finally, we convert the resulting spatial frequency estimate into a DOA estimate, achieving DOA estimation within the limited region.

[0088] Assume that there are K narrowband far-field signal sources in space from directions θ k Arrives at a certain N-element uniform linear array. 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 array element spacing, c is the wave velocity, 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 of the target DOA range. For example, in the multi-frame DOA estimation problem, the target DOA range can be determined by the DOA estimation results of the first few frames. In the estimation of subsequent frames, the target DOA can be considered to exist only within the determined range. When interference exists, the spatial matrix filter can be used to make the incoming signal exist only within the pre-defined range. Assuming that the DOA range of the target is known to be [θ1,θ2], that is, Based on this feature, the target DOA can be searched only within the range of [θ1,θ2], so the algorithm complexity is relatively lower than the case where the DOA range is unknown.

[0091] When the range of the target's DOA is not restricted, in order to estimate the DOA from the received signal x in equation (1), the atomic set can be constructed: And define the atomic norm of x for

[0092]

[0093] Where t>0 is the optimization variable, c k >0 is the optimization variable, a k For atoms, Indicates a k Belong to the set The DOA of the target can be estimated by solving the following atomic norm minimization problem:

[0094]

[0095] When noise does not exist, the above formula is converted to

[0096]

[0097] in, is the optimization variable.

[0098] Since 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 optimization problem (4) is equivalent to the following semidefinite programming problem

[0099]

[0100] Where u0 represents the 0th element of vector u, and Toep(·) is the Toeplitz operator, which converts vector u into the Toeplitz matrix Toep(u). It is easy to prove that the Vandermonde decomposition of the solved Toep(u) yields

[0101]

[0102] where p k >0, is a real number greater than 0. The solution is This is the target DOA.

[0103] Next, we consider the case where we have a priori information about the target DOA range. In the dictionary-based grid parameter estimation framework, this prior information can be achieved by restricting the dictionary, that is, making the discrete division Corresponding to the gridless framework, that is, constructing a new atomic 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 the dual problem of (8) as

[0108]

[0109] Where q is the dual variable, represents the dual atomic norm of q, which can be calculated as

[0110]

[0111] in In order to facilitate the subsequent discussion, 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) we have

[0112]

[0113] The dual trigonometric polynomial is where q k Represents the kth element of q.

[0114] therefore Equivalent to

[0115]

[0116] Note that using the bounded lemma we can Convert to a semi-positive constraint. However, in (12) only for f∈[f L ,f H ] is restricted to the case in [f], so the bounded lemma does not apply. Note that the dual polynomial q(f) can actually be regarded as the discrete-time Fourier transform of the finite support sequence q, so q can be processed to use the bounded lemma. An intuitive approach is to filter the sequence q using a finite impulse response (FIR) bandpass filter to transform it into a sequence with a passband of [f L ,f H ] bandpass signal at this time and Approximately equivalent. However, this method has poor robustness and is difficult to achieve accurate DOA estimation when the signal-to-noise ratio is low. In the present invention, q is first filtered using an ideal bandpass filter. Let the ideal low-pass filter be After a delay of n, it is h n , the ideal bandpass filter operator can be constructed as Where H represents the ideal bandpass filter operator, H°q represents the ideal bandpass filter operator acting on q, q H It represents the output of the ideal bandpass filter operator after acting on q. At this time, due to the output q H It has infinite dimensions, so it is impossible to convert the problem into a finite-sized SDP problem. Assuming that we want to use the bounded lemma to construct an M×M positive semidefinite matrix, we need to design a dimensionality reduction operator Make

[0117]

[0118] in is the discrete time Fourier transform operator. Simulation experiments show that by H It is sufficient to perform a simple truncation to design G. Specifically, assuming an ideal filter with time delay h n The index corresponding to the center of the main lobe is i n , then take the fourth to fifth side lobes in front of i0 and i N-1 The fourth to fifth side lobes can be used on the rear side. n for The output obtained is At this time The conditions for LTI are not met, but better results can be achieved than an FIR bandpass filter using LTI.

[0119] Despite the truncation, the output of the above method is Usually has a higher dimension than q, resulting in a larger SDP problem. In fact, due to The spectrum is limited to [f L ,f H ], its bandwidth B = f H -f L <1, so it is entirely possible to reduce the number of samples by lowering the sampling rate. Specifically, return to q before truncation H , this signal is strictly bandpass, so it can be resampled to fill the frequency band Specifically, first of all, H Phase shifting 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 downsampled by 1 / B; when 1 / B is not an integer, it can be interpolated by an analog low-pass filter to obtain an analog signal, and then resampled according to the Nyquist sampling rate of B. The above process can be directly performed on H, and the new filter operator is H0. Let the filter matrix obtained after truncation with the same criteria as before be The output signal is at this time The dimension N0 of is usually smaller than q, about O(NB). Based on the above discussion, we can Approximately equivalent to From this we can get the SDP form of the dual problem:

[0120]

[0121] The SDP of the original problem is:

[0122]

[0123] Where y and u are optimization variables. By performing Vandermonde decomposition on the solved Toep(u), it can be written as follows:

[0124]

[0125] where a 0 (f 0 )=[1,exp(j2πf 0 ),…,exp(j2πf 0 (N0-1))] T . Now we can extract the frequency According to the previous discussion, the spatial frequency f can be obtained by By converting the spatial frequency back to the incident angle θ, DOA estimation can be achieved.

[0126] When noise and truncation error are taken into account, (15) can be modified to

[0127]

[0128] Where ∈ is the noise tolerance.

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

[0130]

[0131] Where τ is the regularization parameter. In order 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 in this invention mainly considers the following more general atomic norm minimization problem:

[0133]

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

[0135]

[0136] in for The dual cone of . Gao et al. derived

[0137]

[0138] From this we can get the corresponding dual problem:

[0139]

[0140] The duality gap is

[0141] eta(μ,λ)=f(μ)-g(λ). (twenty three)

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

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

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

[0145]

[0146] in It's GG H Eigenvalue decomposition of , where U is a matrix consisting of eigenvectors whose eigenvalues ​​are not 0, each column of which is an eigenvector, and Σ is a diagonal matrix with eigenvalues ​​on the diagonal.

[0147] The goal of this invention is to solve the original problem (19) under the framework of the original dual interior point method. Specifically, this invention uses the logarithmic homogeneous function F(μ) as the original cone The corresponding barrier function 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 is incremented in each iteration. It can be shown that the solution of (27) converges to the solution of (19) when t→∞. The corresponding augmented KKT condition of (27) is

[0152]

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

[0154] From (28) we can get

[0155]

[0156] and

[0157]

[0158] where Λ=Toep -1 (u)X. The original 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 from (30)(a). The dual variables 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] At the same time, multiply both sides of (31) by UU on the left H , we can get

[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, there is

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

[0168] Note 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, we have

[0171]

[0172] Combined with (33), the following results can be obtained:

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

[0174] It is worth mentioning that in (38) only (U H Toep(u)U+τΣ -1 ), the size of the matrix to be inverted is significantly reduced compared to that in (30)(b). It should be noted that usually only the observation Y=G of Y0 can be obtained. H Y0. In order to calculate U appearing in (38) H Y0, can be deduced

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

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

[0177]

[0178] It is worth noting 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 fast Fourier transform, while Toep in (40)(P2) -1 (u)X can be quickly solved using 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 Newton or pseudo-Newton method. The present invention uses the limited memory Broyden-Fletcher-Goldfarb-Shanno (Limited Memory BFGS, L-BFGS) algorithm in the pseudo-Newton method to achieve efficient solution of the search direction. The L-BFGS algorithm requires the Hessian matrix to be initialized, usually initialized to the identity matrix Hansen et al. found that such initialization would cause the algorithm to converge very slowly, so they proposed a heuristic diagonal matrix approximation of the Hessian matrix.

[0183]

[0184] The present invention follows this initialization scheme and The following calculations were performed:

[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. FCA-ANM based on the primal-dual interior point method.

[0188]

[0189] Example 2

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

[0191] The spatial frequency range conversion module is used to limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range;

[0192] The SDP problem conversion module is used to design the corresponding atomic set and the corresponding atomic norm minimization problem based on the limited spatial frequency range, design the corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator;

[0193] The SDP problem solving module is used to solve the SDP problem and obtain the estimation of spatial frequency;

[0194] The DOA estimation module is used to convert the obtained spatial frequency estimation into an estimated value of DOA.

[0195] Simulation experiment:

[0196] Assume that there are three far-field plane wave signals irradiating the array element number N m On the uniform linear array, the DOAs are 0°, 13.2° and -17.8° respectively. The above DOA settings are mainly to reflect the effect of grid mismatch in the grid method. Specifically, the Θ is obtained by dividing the grid into 1° intervals. dic , at this time, the target has a DOA that does not fall on the grid, resulting in grid mismatch. Assume that the target DOA range is known to be [-30°, 30°] before detection. Let the receiving noise be The array element spacing is half wavelength. Aliasing does not occur in this case. This simulation was performed using MATLAB 2021b. The optimization problems involved were solved using the CVX toolkit, using the SDPT3 solver. The simulation was performed using a 12th Gen Intel(R) Core(TM) i7-12700H 2.30GHz CPU with 16.0GB of onboard RAM.

[0197] This simulation mainly uses the following algorithms for DOA estimation: First, the prior knowledge of the DOA range is ignored and the atomic norm minimization method described in (5) is directly used to estimate the DOA. After the estimation, the DOA estimates outside the limited area are excluded. This method is called the ANM method. Second, the DOA estimation is based on the discrete partition Θ dic The constructed dictionary has a grid method. In order to improve the efficiency, the SBL method is used for solution, which is called the constrained area (Constrained Area)-SBL method, abbreviated as CA-SBL method; then the solution is based on the atomic norm minimization method described in (17) proposed by the present invention, which is named the constrained area (Constrained Area)-ANM method, abbreviated as CA-ANM method; finally, the method of solving CA-ANM using the fast solution method based on the primal-dual interior point method proposed by the present invention is called the fast CA-ANM method, abbreviated as FCA-ANM method.

[0198] First, verify the performance of each algorithm in a noisy environment. Verify the number of array elements N m =50 and N m = 100, setting the signal-to-noise ratio (SNR) to -20dB to 30dB and -20dB to 10dB. 100 Monte Carlo simulations are performed for each SNR. The performance of the algorithm can be characterized from two perspectives: the number of successful estimates and the root mean square error (RMSE) of the estimate. For each Monte Carlo simulation, the root square error (RSE) of the estimate is defined as

[0199]

[0200] where N Tar =3 is the target number, is the DOA estimate of the kth target, θ k is the true DOA of the kth target. The estimation is defined as successful when RSE ≤ 5°. For the simulation rounds in which the estimation is successful, the RMSE of the estimation is defined as

[0201]

[0202] in is the estimated number of successes, is the DOA estimate of the kth target in the nth successful estimation. m When N = 50, the curves of the number of successful estimations and RMSE of each algorithm versus SNR are plotted in Figure 3(a) and Figure 3(b), respectively. m When =100, the curves of the number of successful estimations and RMSE of each algorithm versus SNR are plotted in Figure 3(c) and Figure 3(d), respectively.

[0203] When N m =50, as can be seen from Figure 3(a), when the signal-to-noise ratio is higher than 0dB, all methods can successfully achieve DOA estimation, among which the limited area estimation methods CA-SBL, CA-ANM and FCA-ANM proposed in 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 estimation range of DOA is beneficial to improving the robustness of the algorithm. As can be seen from Figure 3(b), the estimation accuracy of all methods is comparable at low signal-to-noise ratios. However, when the signal-to-noise ratio is higher than 0dB, due to the influence of grid mismatch, the performance of the CA-SBL method is no longer significantly improved, while the other gridless methods improve the estimation accuracy as the signal-to-noise ratio increases, showing the superiority of the gridless method. Among the gridless methods, the performance of ANM, CA-ANM and FCA-ANM are almost the same. When N m = 100, as can be seen from Figure 3(c) and Figure 3(d), the trend of the number of successful estimations and RMSE of each algorithm with SNR changes is similar to that of N m =50. However, at the same SNR, N m =100, the number of successful estimations and the accuracy of each algorithm are higher than N m =50, it can be seen that increasing the number of array elements is beneficial to improving the robustness of the algorithm.

[0204] Next, we compare the operation speed of each algorithm, and fix SNR = 30dB. First, we simulate the effect of the number of array elements on the operation speed. m Set to 20 to 100. For each N m , conduct 100 Monte Carlo simulations and record the operation time of each algorithm, and average the time recorded in each simulation as the actual operation time. It is still considered that the estimation is successful when the RSE of a simulation is ≤ 5°. The simulation results show that for each N m , the number of successful estimations for all algorithms is 100. The RMSE and operation time of each algorithm are calculated with the number of array elements N. m The changing curves are plotted in Figure 4(a) and Figure 4(b) respectively.

[0205] As shown in Figure 4(a), due to the influence of mesh mismatch, the estimation accuracy of the meshed CA-SBL algorithm is weaker than that of the meshless ANM, CA-ANM, and FCA-ANM methods. Among the meshless algorithms, the ANM algorithm, which estimates DOA for the entire area, slightly outperforms the proposed DOA estimation algorithm for a limited area. The fast algorithms FCA-ANM and CA-ANM perform similarly, further demonstrating the reliability of the fast algorithms. Because the solution mechanism of the meshed CA-SBL algorithm differs from that of the other meshless algorithms, for the sake of fairness, only the computational time of the three meshless algorithms is compared in Figure 4(b). As shown in Figure 4(b), the computational time of all algorithms increases with the number of array elements. At the same time, for all array element numbers, the computational time of the fast algorithm FCA-ANM proposed in this paper is significantly lower than that of the other algorithms, demonstrating the efficiency of the proposed fast algorithm. Furthermore, it can be seen that when the number of array elements is large, the computational time of the proposed CA-ANM algorithm is significantly lower than that of the ANM algorithm. This is consistent with the previous discussion that the dimension of the dual variable q in the ANM algorithm is N, while the dimension of the dual variable in the CA-ANM algorithm is The dimension of is approximately O(NB), where B<1 is the bandwidth of the spatial frequency. Therefore, the computational complexity of the CA-ANM algorithm will be smaller than that of the ANM algorithm.

[0206] The following simulation examines the effect of the DOA search range on algorithm speed. Maintaining a fixed SNR of 30 dB, the algorithm performed 100 Monte Carlo simulations with the angle range set from 60° (i.e., [-30°, 30°]) to 180° (i.e., the entire search space) centered at 0°. The simulation results show that all algorithms achieve 100 successful estimations within each search range. Figures 5(a) and 5(b) plot the RMSE and computation time of each algorithm as a function of angle range, respectively.

[0207] As shown in Figure 5(a), due to the influence of grid mismatch, the estimation accuracy of the gridded CA-SBL algorithm is weaker than that of the gridless ANM, CA-ANM, and FCA-ANM methods. The accuracy of each algorithm is roughly consistent for different DOA search ranges, indicating that limiting the DOA search range does not affect the accuracy of the algorithm. Figure 5(b) shows that the computation time of the proposed restricted area algorithms, CA-ANM and FCA-ANM, increases with the expansion of the angular search range. When the DOA search range does not exceed 120°, the computation time of the proposed CA-ANM algorithm is less than that of the ANM algorithm, indicating that reducing the DOA search range can improve computational speed, which is consistent with previous analysis. When the DOA search range is large, the spatial frequency bandwidth B approaches 1, and the CA-ANM algorithm does not significantly reduce the algorithm size. Furthermore, the corresponding optimization problem is more complex than that of the ANM algorithm, resulting in a longer computation time. In comparison, the operation speed of the proposed fast algorithm FCA-ANM algorithm 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 high efficiency of the proposed fast algorithm.

[0208] The present application may also provide a computer device comprising: at least one processor, memory, at least one network interface, and a user interface. The various components in the device are coupled together via a bus system. It will be understood that the bus system is used to enable communication between these components. In addition to a data bus, the bus system also includes a power bus, a control bus, and a status signal bus.

[0209] The user interface may include a display, a keyboard, or a pointing device, such as a mouse, a trackball, a touchpad, or a touch screen.

[0210] It is understood that the memory in the embodiments disclosed in the present application may be a volatile memory or a non-volatile memory, or may include both volatile and non-volatile memories. Among them, the non-volatile memory may 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 may be a random access memory (RAM), which is used as an external cache. By way of example and not limitation, many forms of RAM are available, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate synchronous DRAM (DDRSDRAM), enhanced synchronous DRAM (ESDRAM), synchronous link DRAM (SLDRAM), and direct RAM bus RAM (DRRAM). The memories described herein are intended to include, but are not limited to, these and any other suitable types of memory.

[0211] In some embodiments, the memory stores the following elements, executable modules or data structures, or a subset or an extension thereof: an operating system and applications.

[0212] The operating system includes various system programs, such as the framework layer, core library layer, and driver layer, which are used to implement various basic services and handle hardware-based tasks. Application programs include various application programs, such as media players and browsers, which are used to implement various application services. The program that implements the method of the embodiment of the present disclosure can be included in the application program.

[0213] In the above embodiment, the processor may also call a program or instruction stored in the memory, specifically, a program or instruction stored in the application program, to:

[0214] Perform the steps of the above method.

[0215] The above method can be applied to or implemented by a processor. The processor may be an integrated circuit chip with signal processing capabilities. During implementation, each step of the above method can be completed by hardware integrated logic circuits in the processor or by software instructions. The above 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, or discrete hardware components. The above-disclosed methods, steps, and logic block diagrams can be implemented or executed. The general-purpose processor may be a microprocessor or any conventional processor. The steps of the above-disclosed method can be directly implemented and executed by a hardware decoding processor, or by a combination of hardware and software modules in the decoding processor. The software module can be located in a storage medium well-known in the art, such as random access memory, flash memory, read-only memory, programmable read-only memory, electrically erasable programmable memory, registers, etc. The storage medium is located in the memory, and the processor reads the information in the memory and, in conjunction with its hardware, completes the steps of the above method.

[0216] It is understood that the embodiments described herein may be implemented using hardware, software, firmware, middleware, microcode, or a combination thereof. For hardware implementation, the processing unit may 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, or other electronic units or combinations thereof for performing the functions described herein.

[0217] For software implementation, the technology of the present application can be implemented by executing the functional modules (e.g., procedures, functions, etc.) of the present application. The software code can be stored in a memory and executed by a processor. The memory can be implemented in the processor or external to the processor.

[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, each step in the above method embodiment can be implemented.

[0219] Finally, it should be noted that the above embodiments are intended only to illustrate the technical solutions of this application and are not intended to limit the scope of the present invention. Although this application has been described in detail with reference to the embodiments, it should be understood by those skilled in the art that modifications or equivalent substitutions to the technical solutions of this application do not depart from the spirit and scope of the technical solutions of this application and should be encompassed by the claims of this 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: 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; Step 3: Solve the SDP problem to obtain an estimate of the spatial frequency; Step 4: Convert the obtained spatial frequency estimate into an estimated value of DOA; The corresponding mapping operators of the design include: The dual problem of the atomic norm minimization problem is: ； in, is the dual variable; represents the real part of the inner product; express The dual atomic norm of : ； Among them, the dual trigonometric polynomial , express No. elements; Using an ideal bandpass filter right Filtering, ideal bandpass filter delay Later ; Construct an ideal bandpass filter operator , represents the ideal bandpass filter operator acting on , represents the ideal bandpass filter operator acting on The output after right The new filter operator obtained by phase shifting , the frequency band is :when When it is an integer, it is directly downsampled times; when When it is not an integer, it is interpolated by an analog low-pass filter to obtain an analog signal, and then Resample at the Nyquist sampling rate; right The filter matrix obtained after truncation according to the set criteria is: , the output signal is ; It is the mapping operator; The method of converting the atomic norm minimization problem into the SDP problem based on the mapping operator includes: The SDP questions are: ； in, and is the optimization variable; represents the Toeplitz operator; superscript H represents conjugate transpose; Represents a vector The 0th element of is the optimization variable; When considering noise and truncation error, the SDP problem is: ； The objective function is: ； in, is the regularization parameter; Solving the SDP problem includes: Consider the atomic norm minimization problem: ； in, , M represents the dimension of the positive semidefinite matrix, , , ; 、 Indicates the optimization of intermediate variables; Represents a quantity with a known value, namely the received signal containing noise; is the optimal value of the objective function after optimization, which is defined as: ； is the dual cone, defined as ; Solve using the following formula to get and : ； in, A matrix consisting of eigenvectors whose eigenvalues ​​are not 0, where each column is an eigenvector; is a diagonal matrix with eigenvalues ​​on the diagonal; represents the identity matrix; represents the dual variable; Solve using the following formula to get : ； in, is the objective function, defined as ; The limited memory method in the pseudo-Newton method is used to solve the above equation, where the Hessian matrix Initialized as: ； The calculation formula is: 。 2. The method for estimating DOA in a limited area based on atomic norm minimization according to claim 1, wherein: The step 1 comprises: The range of target DOA is determined by the DOA estimation results of the previous set number of frames. ; The incident angle Mapping to spatial frequencies : ； in, , represent the lower and upper bounds of the spatial angular frequency respectively.

3. The method for estimating DOA in a limited area based on atomic norm minimization according to claim 2, wherein: The designing of a corresponding atom set based on a limited spatial frequency range includes: The narrowband far-field signal sources are respectively Arrive at a Received signal of uniform linear array Expressed as: ； in, , , is the array steering vector, the superscript T Transpose of the vector; is the carrier frequency; is the array element spacing, is the wave speed; is the amplitude of the signal source; For noise; j Indicates plural units; Design Atom Set for: 。 4. The method for estimating DOA in a limited area based on atomic norm minimization according to claim 3, wherein: The atomic norm minimization problem is expressed as: When there is noise, the corresponding atomic norm minimization problem is: ； in, is the noise margin; is the optimization variable; for The atomic norm of represents the 2-norm; When there is no noise, the corresponding atomic norm minimization problem is: 。 5. The method for estimating DOA in a limited area based on atomic norm minimization according to claim 1, wherein: The setting criteria are: Take the delayed ideal filter The index corresponding to the center of the main lobe The fourth to fifth side lobes and The fourth to fifth posterior side lobes.

6. A limited area DOA estimation system based on atomic norm minimization, implemented based on the method of any one of claims 1 to 5, characterized in that: The system comprises: The spatial frequency range conversion module is used to limit the DOA range of the incoming wave and convert the limited DOA range into a spatial frequency range; The SDP problem conversion module is used to design the corresponding atomic set and the corresponding atomic norm minimization problem based on the limited spatial frequency range, design the corresponding mapping operator, and convert the atomic norm minimization problem into an SDP problem based on the mapping operator; A SDP problem solving module, used to solve the SDP problem and obtain an estimate of the spatial frequency; and The DOA estimation module is used to convert the obtained spatial frequency estimation into an estimated value of DOA.

Images