DIRECTION OF ARRIVAL ESTIMATION METHOD AND SYSTEM FOR SPARSE ARRAY BASED ON VANDERMONDE DECOMPOSITION RECONSTRUCTION

Abstract

Provided are a direction of arrival (DOA) estimation method and system for a sparse array based on Vandermonde decomposition reconstruction, relating to the technical field of array signal processing. The method includes: constructing a covariance matrix completion optimization model based on sparse array signals and uniform linear array signals, and performing Vandermonde decomposition by using characteristics of a uniform linear array; introducing a nuclear norm to optimize a rank function in the model, and updating the covariance matrix completion optimization model; and introducing an auxiliary variable to transform the model into a solvable optimization problem, and solving the problem by an alternating direction multiplier method to obtain an optimal estimation value. The DOA is estimated using a root multiple signal classification algorithm.

Claims

1. A method for identifying signal sources based on direction of arrival (DOA) estimation for a sparse array based on Vandermonde decomposition reconstruction, which is implemented in a system for identifying signal sources comprising a processor and a memory storing instructions to be executed by the processor to implement the method, wherein, the method comprises: acquiring sparse array signals by a sparse sensor array, wherein the sparse array signals are composed of a plurality of far-field narrow-band and uncorrelated signals entering the sparse sensor array from any directions; acquiring uniform linear array signals, wherein the uniform linear array signals are signals received by a presumed uniform linear array, and the presumed uniform linear array is obtained by transforming the sparse sensor array through an interpolation method; constructing a covariance matrix completion optimization model according to the sparse array signals and the uniform linear array signals, wherein the covariance matrix completion optimization model is a matrix completion model established by using characteristics of Vandermonde decomposition of a covariance matrix of the presumed uniform linear array; wherein constructing the covariance matrix completion optimization model comprises: computing a covariance matrix of the sparse array signals and a covariance matrix of the uniform linear array signals; decomposing the covariance matrix of the presumed uniform linear array into a Vandermonde matrix and a corresponding coefficient vector by using the characteristics of the Vandermonde decomposition; and constructing the covariance matrix completion optimization model according to the Vandermonde matrix and the coefficient vector; wherein a formula expression of the covariance matrix completion optimization model is as follows: $\underset{U,u}{\min }{.Math.}_{r=1}^R\mathrm{rank}\left(\left[\left[U\right]\right]\right)+{.Math.\left(u\right){}^H-{\hat{R}}_S.Math.}_F^2$ $\mathrm{subject}\mathrm{to}{\mathrm{UU}}^H=\left(u\right)$ wherein U=diag([.sub.1, . . . , .sub.L]) is a ||L-dimensional matrix, subject to is a feasible region for limiting an optimization variable, (u) is a covariance matrix and is a Hermitian-Toeplitz matrix, and a first column of the covariance matrix is u, .Math..sub.F represents Frobenius norm, is a Hankel matrix transformation operator, represents a wavelength of incident signals, R indicates a total number of columns, r is an r.sup.-th column, is a column-extraction operator, represents a ||||-dimensional compression matrix only comprising 0 and 1, is a sample covariance matrix, H represents conjugate transpose, is a Vandermonde matrix, and is a square root of power; introducing a nuclear norm for the covariance matrix completion optimization model to replace a rank function in the covariance matrix completion optimization model to obtain an updated covariance matrix completion optimization model; introducing an auxiliary variable for the updated covariance matrix completion optimization model, and transforming the updated covariance matrix completion optimization model into an equivalent form to obtain a solvable optimization problem; wherein introducing the auxiliary variable for the updated covariance matrix completion optimization model to perform equivalent form transformation on the updated covariance matrix completion optimization model to obtain the solvable optimization problem comprises: a formula expression of the solvable optimization problem is as follows: $\underset{U,V,u,B_r}{\min }{.Math.}_{r=1}^R{.Math.B_r.Math.}_{*}+{.Math.\left(u\right){}^H-{\hat{R}}_S.Math.}_F^2;$ wherein subject to UV.sup.H=(u), U=V, [[U]]=B.sub.r, r {1, . . . , R}; solving the solvable optimization problem by an alternating direction multiplier method to obtain an optimal estimation value of the covariance matrix of the sparse array signals; performing DOA estimation by using a root multiple signal classification algorithm according to the optimal estimation value of the covariance matrix of the sparse array signals; and identifying signal sources based on the DOA estimation.

2. The DOA estimation method for the sparse array based on Vandermonde decomposition reconstruction according to claim 1, wherein a formula expression of the sparse array signals is as follows: $x_{}\left(t\right)=\underset{l=1}{\overset{L}{.Math.}}{a}_{}\left({}_l\right)s_l\left(t\right)+n_{}\left(t\right)=A_{}\left(\right)s\left(t\right)+n_{}\left(t\right)$ wherein represents a summation symbol, s(t)=[s.sub.1(t), s.sub.2(t), . . . , s.sub.L(t)].sup.T represents an L-dimensional signal vector, l(1, L); t represents sampling time; (.Math.).sup.T represents a transpose operation; (t) represents a ||-dimensional independent and identically distributed IS complex-valued additive Gaussian noise vector; A.sub.s()=[a.sub.s(.sub.1), a.sub.s(.sub.2), . . . , a.sub.S(.sub.1)] represents a ||L-dimensional array manifold matrix; and a.sub.s(.sub.l) is a ||-dimensional steering vector with an angle of .sub.l.

3. The DOA estimation method for the sparse array based on Vandermonde decomposition reconstruction according to claim 2, wherein a formula expression of the uniform linear array signals is as follows: $x_{}\left(t\right)={.Math.}_{l=1}^L{a}_{}\left({}_l\right)s_l\left(t\right)+n_{}\left(t\right)=A_{}\left(\right)s\left(t\right)+n_{}\left(t\right)$ wherein (t) is a ||-dimensional noise vector; ()=[(.sub.1), (.sub.2), . . . , (.sub.L)] represents a ||L-dimensional array manifold matrix; and (.sub.l) represents a ||-dimensional steering vector with an angle of .sub.l.

4. The DOA estimation method for the sparse array based on Vandermonde decomposition reconstruction according to claim 3, wherein computing the covariance matrix of the sparse array signals comprises: computing the covariance matrix of the sparse array signals according to a formula $R_{}=\left\{x_{}\left(t\right)x^{\frac{H}{}}\left(t\right)\right\}=A_{}\left(\right){\mathrm{PA}}^{\frac{H}{}}\left(\right)+{}^{\frac{2}{n}}I$ wherein {.Math.} represents mathematical expectation, $P=\mathrm{diag}\left(\left[{}^{\frac{2}{1}},{}^{\frac{2}{2}},.Math.,{}^{\frac{2}{L}}\right]\right)$ represents a signal covariance matrix, diag(.Math.) represents a diagonal matrix generated by taking elements of one vector as diagonal elements; [[.sub.l.sup.2]] ${}^{\frac{2}{l}}$ represents power of an l.sup.-th signal, ${}^{\frac{2}{n}}$ is a noise term, and H represents conjugate transpose.

5. The DOA estimation method for the sparse array based on Vandermonde decomposition reconstruction according to claim 4, wherein computing the covariance matrix of the uniform linear array signals comprises: computing the covariance matrix of the uniform linear array signals according to a formula $R_{}=\left\{x_{}\left(t\right)x^{\frac{H}{}}\left(t\right)\right\}=A_{}\left(\right){\mathrm{PA}}^{\frac{H}{}}\left(\right)+{}^{\frac{2}{n}}I=R_s+{}^{\frac{2}{n}}I;$ wherein represents a noise-free covariance matrix of .

6. The DOA estimation method for the sparse array based on Vandermonde decomposition reconstruction according to claim 5, wherein the root multiple signal classification algorithm is as follows: $P_{\mathrm{MUSIC}}\left(\right)=\frac{1}{a^H\left(\right)U_N{U}_N^{H}a\left(\right)},$ wherein U.sub.N represents a noise sub-space of (u), and a() represents a steering vector.

7. (canceled)

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0023] To describe the technical solutions of the present disclosure more clearly, the following briefly introduces the accompanying drawings required for describing the embodiments. Apparently, the accompanying drawings in the following description show merely some embodiments of the present disclosure, and those of ordinary skill in the art may still derive other drawings from these accompanying drawings without creative efforts.

[0024] FIG. 1 is a flowchart of a DOA estimation method for a sparse array based on Vandermonde decomposition reconstruction according to an embodiment of the present disclosure;

[0025] FIG. 2 is a curve of RMSE (root mean square error) according to an embodiment of the present disclosure changing with a parameter R;

[0026] FIG. 3 is a curve of RMSE according to an embodiment of the present disclosure changing with a parameter ;

[0027] FIG. 4 is a schematic diagram of a spatial spectrum of a test algorithm when an incident angle is 0.8 and 0.8 according to an embodiment of the present disclosure;

[0028] FIG. 5A to FIG. 5D are schematic diagrams of spatial spectra of a test algorithm when the number of signal sources according to an embodiment of the present disclosure is 9, where FIG. 5A is a schematic diagram of a spatial spectrum of an SSR algorithm; FIG. 5B is a schematic diagram of a spatial spectrum of a VAI algorithm; FIG. 5C is a schematic diagram of a spatial spectrum of an SNCR algorithm; and FIG. 5D is a schematic diagram of a spatial spectrum of an algorithm provided by the present disclosure;

[0029] FIG. 6A to FIG. 6D are schematic diagrams of spatial spectra of a test algorithm when the number of signal sources according to an embodiment of the present disclosure is 12; FIG. 6A is a schematic diagram of a spatial spectrum of an SSR algorithm; FIG. 6B is a schematic diagram of a spatial spectrum of a VAI algorithm; FIG. 6C is a schematic diagram of a spatial spectrum of an SNCR algorithm; and FIG. 6D is a schematic diagram of a spatial spectrum of an algorithm provided by the present disclosure;

[0030] FIG. 7 is a structural diagram of a DOA estimation system for a sparse array based on Vandermonde decomposition reconstruction according to an embodiment of the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0031] The following clearly and completely describes the technical solution in the embodiments of the present disclosure with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are merely a part rather than all of the embodiments of the present disclosure. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present disclosure without creative efforts shall fall within the scope of protection of the present disclosure.

[0032] In order to make the objectives, features and advantages of the present disclosure more clearly, the present disclosure is further described in detail below with reference to the embodiments.

Embodiment 1

[0033] As shown in FIG. 1, a DOA estimation method for a sparse array based on Vandermonde decomposition reconstruction, including the following steps: [0034] Step 101: acquiring array received signals, where the array received signals include sparse array signals, and uniform linear array signals, the sparse array signals are composed of a plurality of far-field narrow-band and uncorrelated signals entering a sparse array from any directions, the uniform linear array signals are signals received by an presumed uniform linear array, and the presumed uniform linear array is obtained by transforming the sparse array through an interpolation method; [0035] Step 102: constructing a covariance matrix completion optimization model according to the sparse array signals and the uniform linear array signals in the array received signals, where the covariance matrix completion optimization model is a matrix completion model established by using characteristics of Vandermonde decomposition of a covariance matrix of the uniform linear array; [0036] Step 103: introducing a nuclear norm for the covariance matrix completion optimization model to replace a rank function in the covariance matrix completion optimization model to obtain an updated covariance matrix completion optimization model; [0037] Step 104: introducing an auxiliary variable for the updated covariance matrix completion optimization model, and transforming the updated covariance matrix completion optimization model into an equivalent form to obtain a solvable optimization problem; [0038] Step 105: solving the solvable optimization problem by an alternating direction multiplier method to obtain an optimal estimation value of a covariance matrix of the sparse array signals; and [0039] Step 106: performing DOA estimation by using a root multiple signal classification algorithm according to the optimal estimation value of the covariance matrix of the sparse array signals.

[0040] In order to deeply analyze and understand the structural characteristics of array covariance matrix, an optimization model needs to be established. Through this model, the relationship between various sensors in the array and the signal propagation characteristics can be well understood. By the optimization model, the performance of the array can be further optimized, and the accuracy and efficiency of the signal processing are improved.

[0041] In a layout design of array sensors, Nyquist sampling positions are a key concept. The positions refer to the positions of a group of sensor positions with a fixed interval, and the interval is usually determined by a half wavelength of the incident signal. Such a sampling mode can ensure the uniform sampling of the signals in the space, thus avoiding the occurrence of aliasing phenomenon. By accurately setting the Nyquist sampling positions, it can be ensured that the array can effectively capture the details of the signals, thus improving the accuracy and reliability of signal processing.

[0042] Specifically, the Nyquist sampling positions represent the positions of a group of sensors with a fixed interval d, and d is the half-wavelength of the incident signal. may be represented as follows:

[00001]$\begin{array}{cc}=\left\{v_{1}d,v_{2}d,.Math.,v_{.Math..Math.}d\right\}=\left\{0,d,2d,.Math.,\left(.Math..Math.-1\right)d\right\}.& \left(1\right)\end{array}$

[0043] in the formula, ={u.sub.1 d, u.sub.2 d, . . . , d} is a subset of , whereu.sub.1=0.

[0044] In some embodiments, the execution of Step 101 to Step 102 may specifically include the following steps: [0045] acquiring array received signals, where the array received signals include sparse array signals, and uniform linear array signals, the sparse array signals are composed of multiple far-field narrow-band and uncorrelated signals entering a sparse array from any directions, the uniform linear array signals are signals received by an presumed uniform linear array, and the presumed uniform linear array is obtained by transforming the sparse array through an interpolation method; [0046] computing covariance matrixes of the sparse array signals and the uniform linear array signals according to the sparse array signals and the uniform linear array signals in the array received signals; [0047] tranforming the covariance matrix of the uniform linear array into a Vandermonde matrix and a corresponding coefficient vector by using characteristics of the Vandermonde decomposition; and [0048] constructing the covariance matrix completion optimization model according to the Vandermonde matrix and the coefficient vector.

[0049] Specifically, when acquiring the signals, it is assumed that there are L far-field narrow-band and uncorrelated signals incident on the sparse array from directions =[.sub.1, .sub.2, . . . , .sub.L], and the sparse array signals in the array received signals can be modeled as follows:

[00002]$\begin{array}{cc}{x}_{}\left(t\right)=\underset{l=1}{\overset{L}{.Math.}}{a}_{}\left({}_l\right)s_l\left(t\right)+n_{}\left(t\right)=A_{}\left(\right)s\left(t\right)+n_{}\left(t\right)& \left(2\right)\end{array}$

[0050] where represents a summation symbol, s(t)=[s.sub.1(t), s.sub.2(t), . . . , s.sub.L(t)].sup.T represents an L-dimensional signal vector, l(1, L); t represents sampling time; (.Math.).sup.T represents a transpose operation; (t) represents a ||-dimensional independent and identically distributed complex-valued additive Gaussian noise vector, which has a mean of zero and a covariance matrix of

[00003]${}_n^{2}I;{}_n^{2}I$

is noise power, ()=[(.sub.1), (.sub.2), . . . , (.sub.L)] represents a ||L-dimensional array manifold matrix; and a(.sub.l) is a ||-dimensional steering vector with an angle of .sub.l, which is represented as follows:

[00004]$\begin{array}{cc}{a}_{}\left({}_l\right)={\left[1,e^{-j\frac{2}{}{u}_{2}dsin{}_l},.Math.,e^{-j\frac{2}{}{u}_{.Math..Math.}dsin{}_l}\right]}^T& \left(3\right)\end{array}$

[0051] where represents a wavelength of the incident signals, j is an imaginary unit, and e is the base of a natural logarithm.

[0052] According to the acquired sparse array S, the covariance matrix is represented as follows:

[00005]$\begin{array}{cc}{R}_{}=\left\{x_{}\left(t\right)x_{}^H\left(t\right)\right\}=A_{}\left(\right)P{A}_{}^H\left(\right)+{}_n^{2}I& \left(4\right)\end{array}$

[0053] where{.Math.} represents mathematical expectation,

[00006]$P=\mathrm{diag}\left(\left[{}_1^2,{}_2^2,.Math.,{}_L^2\right]\right)$

represents a signal covariance matrix, diag(.Math.) represents a diagonal matrix generated by taking elements of one vector as diagonal elements;

[00007]${}_l^2$

represents power of an l.sup.th signal,

[00008]${}_n^{2}I$

is a noise term, and n represents conjugate transpose.

[0054] However, in practice, the covariance matrix R of the sparse array can only be estimated through the sample covariance matrix, which is expressed as follows:

[00009]$\begin{array}{cc}{\widehat{R}}_{}=\frac{1}{K}{.Math.}_{t=1}^K{x}_{}\left(t\right)x_{}^H\left(t\right)& \left(5\right)\end{array}$

[0055] where K represents the number of time domain samples (snapshot number).

[0056] In order to meet the Nyquist spatial sampling criterion (that is, the sampling frequency must be at least twice the highest frequency component of the signal) while remaining the aperture of the sparse array unchanged, the Nyquist space filling process needs to be implemented to generate an presumed uniform linear array. In this process, by interpolating the sparse array into a continuous uniform linear array , the problem of lacking actual sampling data at the position - is solved, where includes all elements of . In theory, the received signal of the presumed uniform linear array U may be modeled as follows:

[00010]$\begin{array}{cc}{x}_{}\left(t\right)={.Math.}_{l=1}^L{a}_{}\left({}_l\right)s_l\left(t\right)+n_{}\left(t\right)=A_{}\left(\right)s\left(t\right)+n_{}\left(t\right)& \left(6\right)\end{array}$

[0057] where (t) is a ||-dimensional noise vector; ()=[(.sub.1), (.sub.2), . . . , (.sub.1)] represents a ||L-dimensional array manifold matrix; and (.sub.1) represents a ||-dimensional steering vector with an angle of (.sub.l), which is represented as follows:

[00011]$\begin{array}{cc}{a}_{}\left({}_l\right)={\left[1,e^{-j\frac{2}{}{v}_2\mathrm{dsin}{}_l},.Math.,e^{-j\frac{2}{}{v}_{.Math..Math.}\mathrm{dsin}{}_l}\right]}^T& \left(7\right)\end{array}$

[0058] where a received signal of the sparse array S may be represented as follows through x:

[00012]$\begin{array}{cc}{x}_{}\left(t\right)=x_{}\left(t\right)& \left(8\right)\end{array}$

[0059] where is a ||||-dimensional compression matrix including only 0 and 1, and when the i.sup.-th element of overlaps with the j.sup.-th element of , .sub.i,j=1, which represents that an element of the i.sup.-th row and j.sup.-th column of is 1, and other all elements are zero. For example, for a sparse array .sub.example={0, d, 3d, 6d}, and a uniform linear array .sub.example={0, d, 2d, 3d, 4d, 5d, 6d}, the following can be obtained:

[00013]$\begin{array}{cc}{}_{\mathrm{example}}=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]& \left(9\right)\end{array}$

[0060] The covariance matrix of the presumed uniform linear array U is represented as follows:

[00014]$\begin{array}{cc}{R}_{}=\left\{x_{}\left(t\right)x_{}^H\left(t\right)\right\}=A_{}\left(\right){\mathrm{PA}}_{}^H\left(\right)+{}_n^{2}I=R_s+{}_n^{2}I& \left(10\right)\end{array}$

[0061] where represents a noise-free covariance matrix of .

[0062] ||||-dimensional sparse array covariance matrix may be augmented to the ||||-dimensional presumed uniform linear array covariance matrix as follows:

[00015]$\begin{array}{cc}\stackrel{\_}{R_{}}=R_{}={}^H{R}_{}& \left(11\right)\end{array}$

[0063] where =.sup.H is a ||||-dimensional diagonal matrix. It may be noted that for all id belonging to -, the (i+1).sup.-th row and the (i+1).sup.-th column of the augmented covariance matrix are zero vectors. Therefore, the problem of filling the assumed sensors at the Nyquist sampling positions can be transformed into a problem of completing the covariance matrix , thus reconstructing the covariance matrix .

[0064] Hankel matrix transformation operator represents that an n1-dimensional vector x is mapped into (nm+1) x m-dimensional Hankel matrix [x], that is, [x].sub.i,j=x.sub.i+j1, which represents that the (i+j1).sup.-th element of x is equal to an element at the i.sup.-th row and the j.sup.-th column of [x]. An adjoint operator * of represents a (nm+1)m-dimensional matrix X is mapped into an n-dimensional vector *[X], that is [*X].sub.k=.sub.i+j1=kX.sub.i,j, which represents that the k.sup.-th element of *[X] is equal to a sum of elements in X, a sum of whose row subscript and column subscript is equal to k+1. The column-extraction operator represents that a r.sup.th column of the matrix is extracted, that is [X]=[X.sub.:,r], X.sub.:,r represents the r.sup.th column of the matrix X. An adjoint operator of represents that an n-dimensional vector x is mapped into an nm-dimensional matrix, which is defined as follows:

[00016]$\begin{array}{cc}{\left[r_{*}^{}x\right]}_{:,k}=\{\begin{array}{cc}x,& k=r\\ 0,& kr\end{array},r\left\{1,.Math.,m\right\}& \left(12\right)\end{array}$

[0065] where r {1, . . . , m} represents that a value of r is respectively 1, . . . , m.

[00017]$z_l=e^{-j\frac{2}{}\mathrm{dsin}{}_l},$

where l=1, . . . , L. In theory, the noise-free array covariance matrix can be decomposed as follows:

[00018]$\begin{array}{cc}{R}_s={\left[\begin{array}{ccc}1& .Math.& 1\\ z_1& .Math.& z_L\\ .Math.& .Math.& .Math.\\ z_1^{.Math..Math.-1}& .Math.& z_L^{.Math..Math.-1}\end{array}\right]\left[\begin{array}{ccc}{}_1^2& & \\ & & \\ & & {}_L^2\end{array}\right]\left[\begin{array}{ccc}1& .Math.& 1\\ z_1& .Math.& z_L\\ .Math.& .Math.& .Math.\\ z_1^{.Math..Math.-1}& .Math.& z_L^{.Math..Math.-1}\end{array}\right]}^H=A_{}{\left\{\mathrm{diag}\left(\left[{}_1,.Math.,{}_L\right]\right)\right\}}^2{A}_{}^H={\mathrm{UU}}^H& \left(13\right)\end{array}$

[0066] where U=diag ([.sub.1, . . . , .sub.L]) is a ||L-dimensional matrix. The decomposition is called Vandermonde decomposition of a positive semi-definite Toeplitz matrix. As is a Vandermonde matrix, a l.sup.-th column of U correspond to a product of the l.sup.-th column of and one coefficient

[00019]${}_l^2.$

Therefore, in this embodiment, U is called to have a Vandermonde structure.

[0067] Based on the above matrix decomposition, the noise term

[00020]${}_n^{2}I$

is eliminated, the matrix is reconstructed, which involves finding a ||L-dimensional matrix with a Vandermonde structure, thus making =UU.sup.H. However, it is difficult to achieve the direct constraint of adding the Vandermonde structure. x is enabled to represent an NR-dimensional matrix without zero column, and apparently, there is

[00021]${.Math.}_{r=1}^R\mathrm{rank}\left(\left[\left[X\right]\right]\right)R,$

where rank(.Math.) represents the rank of the matrix. If and only if X has the Vandermonde structure,

[00022]${.Math.}_{r=1}^R\mathrm{rank}\left(\left[\left[X\right]\right]\right)=R.$

Therefore, minimizing

[00023]${.Math.}_{r=1}^R\mathrm{rank}\left(\left[\left[U\right]\right]\right)$

can encourage to have the Vandermonde structure. In addition, is a Hermitian-Toeplitz matrix, and available information in the augmented covariance is used as reference information in the reconstruction process. Therefore, the covariance matrix of the presumed uniform linear array can be reconstructed by solving the following optimization problem:

[00024]$\begin{array}{cc}\underset{U,u,}{\min }{.Math.}_{r=1}^R\mathrm{rank}\left(\left[\left[U\right]\right]\right)+{.Math.\left(u\right){}^H-{\hat{R}}_{}.Math.}_F^2& \left(14\right)\end{array}$$\mathrm{subject}\mathrm{to}{\mathrm{UU}}^H=T\left(u\right)$

[0068] U=diag([.sub.1, . . . , .sub.L]) is a ||L-dimensional matrix, subject to is used to limit a feasible region for limiting an optimization variable, (u) is a covariance matrix, which is a Hermitian-Toeplitz matrix, and a first column of the covariance matrix is u, .Math..sub.F represents Frobenius norm; a covariance fitting error is used to tolerate the noise term, and a regularization parameter is used to weigh the covariance fitting error and a rank function; the is a Hankel matrix transformation operator, represents a wavelength of an incident signal, R indicates a total number of columns, is a column-extraction operator, represents a ||||-dimensional compression matrix only including 0 and 1, and is as shown in a formula (5).

[0069] The execution of Step 103 to Step 104 may specifically performed as follows:

[0070] As the rank function makes the minimum problem difficult to be solved, a nuclear norm is introduced as a convex relaxation of the rank. The optimization problem (14) may be reformulated as follows:

[00025]$\begin{array}{cc}\underset{U,u,}{\min }{.Math.}_{r=1}^R{.Math.\left[\left[U\right]\right].Math.}_{*}+{.Math.\left(u\right){}^H-{\hat{R}}_{}.Math.}_F^2& \left(15\right)\end{array}$$\mathrm{subject}\mathrm{to}{\mathrm{UU}}^H=T\left(u\right)$

[0071] where .Math., represents the nuclear norm.

[0072] In practice, the number L of the signal sources is generally unknown. Therefore, the algorithm in this embodiment needs to preset a parameter R to determine the size of the matrix U. In order to solve the optimization problem conveniently, an auxiliary variable is introduced and (15) is rewritten as the following equivalent form:

[00026]$\begin{array}{cc}\underset{U,V,u,B_r}{\min }\underset{r=1}{\overset{R}{.Math.}}{.Math.B_r.Math.}_{*}+{.Math.\left(u\right){}^H-{\hat{R}}_{}.Math.}_F^2& \left(16\right)\end{array}$$\mathrm{subject}\mathrm{to}{\mathrm{UV}}^H=\left(u\right),U=V,\left[\left[U\right]\right]=B_r,r\left\{1,.Math.,R\right\}.$

[0073] In some embodiments, Step 105 may specifically performed as follows:

[0074] An alternating direction multiplier method is used to solve the optimization problem, where an augmented Lagrangian function of an optimization problem (16) is as follows:

[00027]$\begin{array}{cc}{}_{}\left(U,V,u,B_r,C,D,E_r\right)=\underset{r=1}{\overset{R}{.Math.}}{.Math.B_r.Math.}_{*}+{.Math.\left(u\right){}^H-{\hat{R}}_{}.Math.}_F^2+.Math.C,{\mathrm{UV}}^H-\left(u\right).Math.+\frac{}{2}{.Math.{\mathrm{UV}}^H-\left(u\right).Math.}_F^2+.Math.D,U-V.Math.+\frac{}{2}{.Math.U-V.Math.}_F^2+\underset{r=1}{\overset{R}{.Math.}}.Math.E_r,U-B_r.Math.+\frac{}{2}{.Math.R\left[\left[U\right]\right]-B_r.Math.}_F^2& \left(17\right)\end{array}$

[0075] in formula (17), represents a Lagrangian function, C, D and E.sub.r are Lagrangian multipliers, in <X, Y>=tr(X.sup.HY), tr(.Math.) denote traces of the matrix; and u is a penalty parameter.

[0076] Next, the variable Vis initialized in this embodiment according to the following ways.

[0077] First, is subjected to feature value decomposition, that is =QQ.sup.1. Then, the variable is initialized as V.sup.0=Q.sub.:,1:R, and Q.sub.:,1:R represents taking the first column to the R-th column of Q.

[0078] The alternating direction multiplier method employs the following iterative solution to solve the optimization problem (16):

[00028]$\begin{array}{cc}\underset{U}{\min }{}_{}\left(U,V^k,u^k,B_r^k,C^k,D^k,E_r^k\right)& \left(18\right)\end{array}$$\begin{array}{cc}\underset{V}{\min }{}_{}\left(U^{k+1},V,u^k,B_r^k,C^k,D^k,E_r^k\right)& \left(19\right)\end{array}$$\begin{array}{cc}\underset{U}{\min }{}_{}\left(U^{k+1},V^{k+1},u,B_r^k,C^k,D^k,E_r^k\right)& \left(20\right)\end{array}$$\begin{array}{cc}\underset{B_r}{\min }{}_{}\left(U^{k+1},V^{k+1},u^{k+1},B_r,C^k,D^k,E_r^k\right)& \left(21\right)\end{array}$$\begin{array}{cc}{C}^{k+1}=C^k+{}^k\left[{U^{k+1}\left(V^{k+1}\right)}^H-\left(u^{k+1}\right)\right]& \left(22\right)\end{array}$$\begin{array}{cc}{D}^{k+1}=D^k+{}^k\left(U^{k+1}-V^{k+1}\right)& \left(23\right)\end{array}$$\begin{array}{cc}{E}_r^{k+1}=E_r^k+{}^k\left(\left[\left[U^{k+1}\right]\right]-B_r^{k+1}\right),r\left\{1,.Math.,R\right\}.& \left(24\right)\end{array}$

[0079] (1) U is updated, which is specifically as follows:

[0080] In order to update U, a subproblem (18) needs to be solved, a form of which may be represented as follows:

[00029]$\begin{array}{cc}\underset{U}{\min }{.Math.{U\left(V^k\right)}^H-\left(u^k\right)+\frac{1}{{}^k}{C}^k.Math.}_F^2+{.Math.}_{r=1}^R{.Math.R\left[\left[U\right]\right]-B_r^k+\frac{1}{{}^k}{E}_r^k.Math.}_F^2+{.Math.U-V^k+\frac{1}{{}^k}{D}^k.Math.}_F^2.& \left(25\right)\end{array}$

[0081] As a target function of the optimization problem (25) is a convex function, the conjugate matrix U* of U is defined, a partial derivative of the target function about U* is 0, thus obtaining the following linear equation to solve the problem:

[00030]$\begin{array}{cc}{U\left(V^k\right)}^H{V}^k+\underset{r=1}{\overset{R}{.Math.}}{Q}_r^{*}\left[R^{*}\left[R\left[\left[U\right]\right]\right]\right]+U=\left(u^k\right)V^k-\frac{1}{{}^k}{C}^k{V}^k+V^k-\frac{1}{{}^k}{D}^k+{.Math.}_{r=1}^R{Q}_r^{*}\left[R^{*}\left[B_r^k-\frac{1}{{}^k}{E}_r^k\right]\right]& \left(26\right)\end{array}$

[0082] According to the definition of the operators and , it can be obtained that:

[00031]$\begin{array}{cc}*\left[\left[x\right]\right]=x.& \left(27\right)\end{array}$

[0083] In formula (27), is a vector, a k.sup.-th element of which is the number of elements on a k.sup.-th back-diagonal of the Hankel matrix [x], and represents Hadamard product. According to the definition of the operators and the following formula can be obtained:

[00032]$\begin{array}{cc}{.Math.}_r^L{Q}_r^{*}\left[{}^{*}\left[\left[Q_r\left[U\right]\right]\right]\right]=TU;& \left(28\right)\end{array}$

[0084] where T is a matrix, each column of which is . Therefore, the formula (26) may be represented as follows:

[00033]$\begin{array}{cc}{U_{i,:}\left(V^k\right)}^H{V}^k+U_{i,:}\mathrm{diag}\left(T_{i,:}\right)+U=Y_{i,:};& \left(29\right)\end{array}$

[0085] where Y is a right part of the expression (29). Therefore, a closed-form solution of each row of U can be obtained as follows:

[00034]$\begin{array}{cc}{U}_{i,:}^{k+1}={Y_{i,:}\left[{\left(V^k\right)}^H{V}^k+\mathrm{diag}\left(T_{i,:}\right)+I\right]}^{-1}.& \left(30\right)\end{array}$

[0086] Where I represents an RR-dimensional identity matrix, and (.Math.).sup.1 represents an inverse operation of the matrix.

[0087] (2) V is updated, which may be specifically as follows:

[0088] In order to update V, a subproblem (19) can be represented as follows:

[00035]$\begin{array}{cc}\underset{V}{\min }{.Math.U^{k+1}{V}^H-\left(u^k\right)+\frac{1}{{}^k}{C}^k.Math.}_F^2+{.Math.U^{k+1}-V+\frac{1}{{}^k}{D}^k.Math.}_F^2.& \left(31\right)\end{array}$

[0089] Similar to the update of U, a closed-form solution of V is represented as follows:

[00036]$\begin{array}{cc}{V}^{k+1}={\left\{\left[{\left(\left(u^k\right)\right)}^H+I-\frac{1}{{}^k}{\left(C^k\right)}^H\right]U^{k+1}+\frac{1}{{}^k}{D}^k\right\}\left[{\left(U^{k+1}\right)}^H{U}^{k+1}+I\right]}^{-1}.& \left(32\right)\end{array}$

[0090] (3) u is updated, which is specifically as follows:

[0091] In order to update u, a subproblem (20) can be represented as follows:

[00037]$\begin{array}{cc}\underset{u}{\min }{.Math.\left(u\right){}^H-{\hat{R}}_{}.Math.}_F^2+\frac{{}^k}{2}{.Math.\left(u\right)-{U^{k+1}\left(V^{k+1}\right)}^H-\frac{1}{}{C}^k.Math.}_F^2.& \left(33\right)\end{array}$

[0092] As a target function of the optimization problem (33) is a convex function, its partial derivative about u* is 0, thus obtaining a solution of the subproblem as follows:

[00038]$\begin{array}{cc}u={\left(\mathrm{diag}\left(d\right)+\frac{{}^k}{2}D\right)}^{-1}\left(r+\frac{{}^k}{2}z\right).& \left(34\right)\end{array}$

[0093] Specifically, variables involved in the formula (34) are defined as follows:

[00039]$\begin{array}{cc}{d}_i=\underset{j=1}{\stackrel{.Math..Math.-i+1}{.Math.}}\left[{\left({}^{H}J\right)}_{j,j+i-1}+{\left({}^{H}J\right)}_{j+i-1,j}^{*}\right]& \left(35\right)\end{array}$$z_i=\underset{j=1}{\stackrel{.Math..Math.-i+1}{.Math.}}\left(Z_{j,j+i-1}+Z_{j+i-1,j}^{*}\right)$$r_i=\underset{j=1}{\stackrel{.Math..Math.-i+1}{.Math.}}\left[{\left({}^H{\hat{R}}_S\right)}_{j,j+i-1}+{\left({}^H{\hat{R}}_{}\right)}_{j+i-1,j}^{*}\right]$$i\left\{1,2,.Math.,.Math..Math.\right\}$$\begin{array}{cc}D=\mathrm{diag}\left(\left[.Math..Math.,2\left(.Math..Math.-1\right),2\left(.Math..Math.-2\right),.Math.,2\right]\right),& \left(36\right)\end{array}$

[0094] in the formula, d, z and r are ||-dimensional vectors, d.sub.i, z.sub.i and r.sub.i represent the i-th elements of d, z and r

[00040]$Z={U^{k+1}\left(V^{k+1}\right)}^H+\frac{1}{}{C}^k,$

J represents a ||||-dimensional matrix with all elements of 1.

[0095] (4) B.sub.r is updated, which may be specifically as follows:

[0096] Finally, in order to update B.sub.r, a solution of B.sub.r can be obtained by solving the following problems:

[00041]$\begin{array}{cc}\underset{B_r}{\min }{.Math.B_r.Math.}_{*}+\frac{{}^k}{2}{.Math.\left[\left[U^{k+1}\right]\right]+\frac{1}{{}^k}{E}_r^k-B_r.Math.}_F^2.& \left(37\right)\end{array}$

[0097] A solution of the problem (37) is represented as follows through a soft threshold operator;

[00042]$\begin{array}{cc}{B}_r^{k+1}={}_{\frac{1}{{}^k}}\left(\left[\left[U^{k+1}\right]\right]+\frac{1}{{}^k}{E}_r^k\right).& \left(38\right)\end{array}$

[0098] In such an iterative algorithm, in order to accelerate convergence, a strategy of .sup.k+1=.sup.k is adopted in this embodiment to gradually increase the parameter u, and a reconstructed covariance matrix (u) is obtained after completing iteration.

[0099] In some embodiments, Step 106 may be performed as follows:

[0100] DOA estimation is carried out by using a multiple signal classification algorithm.

[0101] Specifically, the DOA estimation is carried out using a root multiple signal classification algorithm (Root-MUSIC) according to the reconstructed covariance matrix (u). A spatial spectrum of multiple signal classification (MUSIC) is as follows:

[00043]$\begin{array}{cc}{P}_{\mathrm{Music}}\left(\right)=\frac{1}{a^H\left(\right)U_N{U}_N^{H}a\left(\right)};& \left(39\right)\end{array}$

[0102] In the formula, U.sub.N represents a noise sub-space of (u); a() represents a steering vector. By searching the maximum {circumflex over (L)} peak values in P.sub.MUSIC, a DOA estimation result can be obtained. The number {circumflex over (L)} of the signal sources may be estimated using Akaike Information Criterion (AIC), the minimum description length (MDL), or a Gerschgorin disk estimation (GDE) method.

[0103] In conclusion, the whole technical solution of the present disclosure is as follows:

[0104] An input is sparse array received signals

[00044]${\left\{x_s\left(t\right)\right\}}_{t=1}^k;$

an output is a DOA estimation result {circumflex over ()}.sub.l, l=1, 2, . . . , L; and the variables are initialized as U, V, u, C, D, E.sub.r, .sub.r {1, . . . , R}, , .sup.0, , k=0.

[0105] When k<100, the following operations are performed: updating U according to the formula (30); updating V according to the formula (32); updating u according to the formula (34); for r=1, . . . , R, updating B, according to the formula (38); for r=1, . . . , R, updating E, according to the formula (24); updating C, D and k=k+1 according to the formula (22) and the formula (23); and then executing the root-MUSIC algorithm on (u) for DOA estimation.

[0106] Specifically, an embodiment is provided to illustrate the method involved in the present disclosure as follows:

[0107] Considering a coprime array composed of seven physical sensors, the positions of the sensors are ={0,3d, 5d, 6d, 9d, 10d, 12d}. In the following simulation test, for the algorithm provided by the present disclosure, the parameters are set as follows: .sup.0=110.sup.3 and =1.8, which are heuristically selected based on practical experience. As the number of the signal sources is unknown, the parameter R is also set as the maximum value R 13. The snapshot number is set as K=500.

[0108] The estimation accuracy of a test algorithm is evaluated by root mean square error (RMSE), which is defined as follows:

[00045]$\begin{array}{cc}\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{NL}}{.Math.}_{n=1}^N{.Math.}_{l=1}^L{\left({\hat{}}_{n,l}-{}_l\right)}^2};& \left(40\right)\end{array}$

[0109] where N represents the number of Monte Carlo tests, L represents the number of the incident signals, {circumflex over ()}.sub.n,l represents an l.sup.-th estimation angle of the n.sup.-th test, and .sub.l represents a true angle of the l.sup.-th incident signal.

[0110] The influence of the parameter and the parameter R on the estimation accuracy of the technical method of the present disclosure are tested at first. It is set that the incident signal has an incident angle of 0 and a signal noise ratio (SNR) of 10 dB, and 5000 Monte Carlo tests are carried out. A curve of the RMSE changing with the parameter R is shown in FIG. 2. It can be observed that even if R is set to be greater than the number of the actual signal sources, the estimation accuracy of the technical method provided by the present disclosure is just slightly decreased. Then, 500 Monte Carlo tests are carried out respectively with respect to the signal noise ratios being set as 10 dB, 0 dB, 10 dB, and 20 dB. Curves of RMSE changing with the parameter under different signal noise ratios are shown in FIG. 3. Apparently, for different SNRs, there is an optimal parameter which makes the RMSE minimum. In addition, the higher the SNR, the greater the parameter set for the RMSE to obtain the minimum value.

[0111] Next, the estimation performance of the technical method provided by the present disclosure is compared with that of a sparse signal reconstruction (SSR) algorithm, a virtual array-based interpolation (VAI) algorithm and a structured Nyquist correlation reconstruction (SNCR) algorithm. The resolution performance of the test algorithms is compared at first. It is ssumed that two signal sources respectively have incident angles 0.8 and 0.8 and thus have a small angular interval, and have signal noise ratios of 0 dB. As shown in FIG. 4, spatial spectra of the test algorithms are displayed. Other algorithms cannot distinguish these two closely spaced signal sources. In contrast, the technical method of the present disclosure can distinguish the two closely spaced signal sources. The comparison result indicates that the technical method provided by the present disclosure has better resolution performance.

[0112] Finally, the performance of the test methods in the estimation of multiple signal sources is compared. It is considered that the nine signal sources are uniformly distributed in the range of [50, 50] and the signal noise ratio is 0 dB. Spatial spectra of the test algorithms are shown in FIG. 5A to FIG. 5D. SSR produces multiple pseudo peaks in the spatial spectrum, while the other three algorithms can accurately estimate DOA of the nine signal sources. Afterwards, the number of the signal sources is increased to L=12 to evaluate the maximum achievable degree of freedom of the technical method provided by the present disclosure. Spatial spectra of the test algorithms are shown in FIG. 6A to FIG. 6D. Apparently, except the SSR algorithm, all algorithms have successfully identified all twelve signal sources. However, the VAI algorithm produces several estimation results that obviously deviate from true directions, while the SNCR algorithm and the technical method provided by the present disclosure achieve higher estimation accuracy.

Embodiment 2

[0113] As shown in FIG. 7, this embodiment provides a DOA estimation system for a sparse array based on Vandermonde decomposition reconstruction, including: [0114] a signal acquisition module 701, configured to acquire array received signals, where the array received signals include sparse array signals, and uniform linear array signals, the sparse array signals are composed of a plurality of far-field narrow-band and uncorrelated signals entering a sparse array from any directions, the uniform linear array signals are signals received by an presumed uniform linear array, and the presumed uniform linear array is obtained by transforming the sparse array through an interpolation method; [0115] a model construction module 702, configured to construct a covariance matrix completion optimization model according to the sparse array signals and the uniform linear array signals in the array received signals, where the covariance matrix completion optimization model is a matrix completion model established by using characteristics of Vandermonde decomposition of a covariance matrix of the uniform linear array; [0116] a first model optimization model 703, configured to introduce a nuclear norm for the covariance matrix completion optimization model to replace a rank function in the covariance matrix completion optimization model, thus obtaining an updated covariance matrix completion optimization model; [0117] a second model optimization module 704, configured to introduce an auxiliary variable for the updated covariance matrix completion optimization model, and transform the updated covariance matrix completion optimization model into an equivalent form to obtain a solvable optimization problem; [0118] a computing module 705, configured to solve the solvable optimization problem by an alternating direction multiplier method to obtain an optimal estimation value of a covariance matrix of the sparse array signals; and [0119] a DOA estimation module 706, configured to perform DOA estimation by using a root multiple signal classification algorithm according to the optimal estimation value of the covariance matrix of the sparse array signals.

[0120] In conclusion, the present disclosure has the following beneficial effects:

[0121] In the prior art, the sparse signal reconstruction (SSR) algorithm has the problems of basis mismatch, and poor estimation performance when the number of signal sources exceeds the number of sensors in the sparse array, and the virtual array-based interpolation (VAI) algorithm and the structured Nyquist correlation reconstruction (SNCR) algorithm may not be able to distinguish different incident signals for the signal sources with small angular intervals. When the angular interval between the incident signals is small and the number of the signal sources exceeds the number of the sensors in the sparse array, the method and system provided by the present disclosure can achieve better spatial resolution and maximize using the extra degree of freedom provided by the sparse array.

[0122] The technical features of the above embodiments can be combined at will. In order to make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, it should be considered that these combinations of technical features fall within the scope recorded in this specification provided that these combinations of technical features do not have any conflicts.

[0123] Specific examples are used herein for illustration of the principles and embodiments of the present disclosure. The description of the embodiments is merely used to help illustrate the method and its core principles of the present disclosure. In addition, a person of ordinary skill in the art can make various modifications in terms of specific embodiments and scope of application in accordance with the teachings of the present disclosure. In conclusion, the content of this specification shall not be construed as a limitation to the present disclosure.