TWO-DIMENSIONAL DIRECTION-OF-ARRIVAL ESTIMATION METHOD FOR COPRIME SURFACE ARRAY BASED ON VIRTUAL DOMAIN TENSOR FILLING

Abstract

Disclosed in the present invention is a two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling, which mainly solves the problems of the loss of multi-dimensional signal structural information and the inability to fully utilize virtual domain statistics in the existing method. The steps thereof are as follows: constructing a coprime surface array; modeling a tensor of a received signal of the coprime surface array; constructing an augmented non-continuous virtual surface array based on cross-correlation tensor transformation of the coprime surface array; deriving a virtual domain tensor based on mirror extension of the non-continuous virtual surface array; dispersing contiguous missing elements by reconstructing the virtual domain tensor; filling the virtual domain tensor based on the minimization of a tensor kernel norm; and decomposing a filled virtual domain tensor to obtain a direction-of-arrival estimation result.

Description

BACKGROUND

Technical Field

The present invention belongs to the technical field of array signal processing and relates to a statistical signal processing technology based on sparse array virtual domain second-order high-dimensional statistics, in particular to a two-dimensional direction-of-arrival estimation method for coprime surface array based on virtual domain tensor filling, which can be used for target positioning.

Description of Related Art

As a sparse array with a systematic structure, a coprime array has the advantages of a large aperture, and high resolution. It can break through the performance bottleneck of traditional uniform array direction-of-arrival estimation in estimation performance and cost overhead. Since sparsely arranged array elements of the coprime array do not meet a Nyquist sampling rate, in order to realize a Nyquist matching direction-of-arrival estimation, a common practice is to calculate second-order statistics of a received signal of the coprime array to construct an augmented non-contiguous virtual array, and extract a continuous part therefrom to realize Nyquist matching processing based on a virtual domain second-order equivalent signal. Further, in order to make full use of all the non-contiguous virtual array elements, an existing method fills a non-contiguous virtual array to improve the performance of the direction-of-arrival estimation. However, the above method usually expresses the received signal as a vector, and derives the virtual domain second-order equivalent signal by vectorizing a covariance matrix of the received signal; in a scenario where a coprime surface array is deployed, since the received signal of the coprime surface array covers two-dimensional spatial information, this vectorized signal processing method destroys an original spatial information structure of the received signal of the coprime surface array, which will cause serious performance loss.

In order to preserve the structured information of the multi-dimensional received signal, a tensor, as a multi-dimensional data type, has been used in the field of array signal processing to characterize a received signal covering the multi-dimensional spatial information, and to perform feature analysis and effective information extraction, so as to achieve a high-precision and high-resolution direction-of-arrival estimation. However, when it relates to the statistic processing of the virtual field tensor of the coprime surface array, the augmented multi-dimensional non-contiguous virtual array would have holes in the whole piece, resulting in the corresponding virtual field tensor with contiguous missing elements. The traditional tensor filling method applied to image inpainting is premised on a random distribution of the missing elements in the tensor, so it cannot effectively fill the virtual field tensor. Therefore, for a virtual domain tensor model of the coprime surface array, how to effectively use all the non-continuous virtual domain tensor statistics information to achieve high-precision and high-resolution two-dimensional direction-of-arrival estimation is still an urgent problem to be solved.

SUMMARY

The purpose of the present invention is to propose a two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling in order to solve the problems of loss of multi-dimensional signal structure information and inability to fully utilize virtual domain statistics in existing methods. It provides a feasible idea and an effective solution to realize the high-precision and high-resolution two-dimensional direction-of-arrival estimation of Nyquist matching by making full use of all the non-continuous virtual domain tensor statistics information corresponding to the coprime surface array.

The purpose of the present invention is to realize through the following technical solutions: a two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling, wherein the method comprises the following steps:

(1) using 4M_xM_y+N_xN_y−1 physical antenna array elements by a receiving end, and performing constructing according to a structure of a coprime surface array, wherein M_x, N_x and M_y, N_y are a pair of coprime integers respectively; decomposing the coprime surface array into two sparse uniform sub-surface arrays ₁ and ₂, wherein ₁ contains 2M_x×2M_y antenna array elements, array element spacings in an x axial direction and a y axial direction are respectively N_xd and N_yd, ₂ includes N_x×N_y antenna array elements, array element spacings in the x axial direction and the y axial direction are respectively M_xd and M_yd, and the unit interval d is taken as half of the wavelength λ of an incident narrowband signal, that is, d=λ/2;

(2) assuming that there are K far-field narrowband uncorrelated signal sources from {(θ₁, φ₁), (θ₂, φ₂), . . . , (θ_K, φ_K)} directions, θ_k and φ_k are respectively an azimuth angle and an elevation angle of a kth incident signal source, k=1,2, . . . , K, utilizing a three-dimensional tensor ∈^2Mx×2My×T to express T sampling snapshot signals of a sparse uniform sub-surface array ₁ as follows:

$X_{{\mathbb{P}}_1}=\sum _{k=1}^K{a}_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ s_k+N_{{\mathbb{P}}_1},$

wherein s_k=[s_k,1, s_k,2, . . . , s_k,T]^T is a multi-snapshot sampling signal waveform corresponding to the kth incident signal source, [⋅]^T represents a transpose operation, ∘ represents the outer product of a vector, N₁ is a noise tensor independent of each signal source, (θ_k, φ_k) and (θ_k, φ_k) are respectively steering vectors of ₁ in an x axial direction and a y axial direction, correspond to a signal source with an incoming wave direction of (θ_k, φ_k), and are expressed as follows:

$a_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi x_{{\mathbb{P}}_1}^{\left(2\right)}{\mu }_k},\cdots ,e^{-j\pi x_{{\mathbb{P}}_1}^{\left(2M_x\right)}{\mu }_k}\right]}^T,a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi y_{{\mathbb{P}}_1}^{\left(2\right)}{v}_k},\cdots ,e^{-j\pi y_{{\mathbb{P}}_1}^{\left(2M_x\right)}{v}_k}\right]}^T,$

and wherein {, , . . . , }{, , . . . , } represent respectively actual positions of physical antenna elements of the sparse uniform sub-surface array ₁ in the x axial direction and y axial direction, and =0, =0, μ_k=sin(φ_k)cos(θ_k), ν_k=sin(φ_kk)sin(θ_k), j=√{square root over (−1)};

expressing the T sampled snapshot signals of the sparse uniform sub-surface array ₂ by another three-dimensional tensor ∈^Nx×Ny×T as follows:

$X_{{\mathbb{P}}_2}=\sum _{k=1}^K{a}_x^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)\circ s_k+N_{{\mathbb{P}}_2},$

wherein a noise tensor independent of each signal source, (θ_k, ϕ_k) and (″_k, φ_k) are respectively steering vectors of ₂ in an x axial direction and a y axial direction, correspond to a signal source with an incoming wave direction of (θ_k, φ_k), and are expressed as follows:

$a_x^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi x_{{\mathbb{P}}_2}^{\left(2\right)}{\mu }_k},\cdots ,e^{-j\pi x_{{\mathbb{P}}_2}^{\left(N_x\right)}{\mu }_k}\right]}^T,a_y^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi y_{{\mathbb{P}}_2}^{\left(2\right)}{v}_k},\cdots ,e^{-j\pi y_{{\mathbb{P}}_2}^{\left(N_x\right)}{v}_k}\right]}^T,$

and wherein {, , . . . , }{, , . . . , } represent respectively actual positions of physical antenna elements of the sparse uniform sub-surface array ₂ in the x axial direction and y axial direction, and =0, =0;

obtaining a second-order cross-correlation tensor ∈^{2Mx×2My×Nx×Ny} by solving cross-correlation statistic of three-dimensional tensors and :

$R_{{\mathbb{P}}_1{\mathbb{P}}_2}=E\left[<X_{{\mathbb{P}}_1},X_{{\mathbb{P}}_2}^{*}{>}_3\right]=\sum _{k=1}^K{\sigma }_k^2{a}_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_x^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\circ a_y^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)+N_{{\mathbb{P}}_1{\mathbb{P}}_2},\mathrm{wherein}{\sigma }_k^2=E\left[s_k{s}_k^{*}\right]$

represents power of a kth incident signal source, =E[<,>₃] represents the cross-correlation noise tensor, <⋅,⋅>_r represents a tensor contraction operation of the two tensors along a rth dimension, E[⋅] represents a mathematical expectation operation, and (⋅)* represents a conjugation operation; the cross-correlation noise tensor only has an element with a value σ_n² in the (1, 1, 1, 1)th position, wherein σ_n² represents the noise power, and elements in other positions have the same value 0;

(3) defining dimension sets J₁={1, 3}, J₂={2,4}, and obtaining a virtual domain signal U_W∈^2MxNx×2MyNy by performing a tensor transformation of dimension merging on the cross-correlation tensor :

$U_W\stackrel{\Delta }{=}{R}_{{\mathbb{P}}_1{\mathbb{P}}_{2_{\left\{J_1,J_2\right\}}}}=\sum _{k=1}^K{\sigma }_k^2{b}_x\left(k\right)\circ b_y\left(k\right),$

wherein by respectively forming a difference set array on the exponential term,

$b_x\left(k\right)=a_x^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\otimes a_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\mathrm{and}{b}_y\left(k\right)=a_y^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\otimes a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)$

construct a two-dimensional augmented virtual surface array along the x axial direction and the y axial direction, ⊗ represents the Kronecker product; therefore, U_W corresponds to a non-continuous virtual surface array W of size J_Wx×J_Wy, J_Wx=3M_xN_x−M_x−N_x+1, J_Wy=3M_yN_y−M_y−N_y+1, and the non-continuous virtual surface array W contains holes in an entire row and an entire column, that is, missing elements;

(4) constructing a virtual surface array W that mirrors the non-continuous virtual surface array W about the coordinate axis, and superimposing the W and W on a third dimension into a three-dimensional non-continuous virtual cubic array of size , wherein and =2; correspondingly, rearranging elements in the conjugate transposed signal U_W * of the virtual domain signal U_W to correspond to positions of virtual array elements in W, so as to obtain a virtual domain signal U_Wcorresponding to the virtual surface array W; superimposing U_W and U_Win the third dimension to obtain the virtual domain tensor corresponding to the non-continuous virtual cubic array , which is represented as:

$U_{\mathbb{Q}}=\sum _{k=1}^K{\sigma }_k^2{\tilde{b}}_x\left(k\right)\circ {\tilde{b}}_y\left(k\right)\circ c_k,$

wherein {tilde over (b)}_x(k) and {tilde over (b)}_y(k) are respectively steering vectors of the non-continuous virtual cubic array on the x axis and the y axis, and correspond to the signal source with an incoming wave direction (θ_k, φ_k); due to existence of the holes in , {tilde over (b)}_x(k) and {tilde over (b)}_y(k) respectively correspond to elements in the hole positions in in the x axial and y axial directions which are set to be zero,

$c_k={\left[1,e^{-j\pi \left(-\left(M_x{M}_y+M_x+M_y\right){\mu }_k-\left(N_x{N}_y+N_x+N_y\right)v_k\right)}\right]}^T$

represents a mirror transformation factor vector corresponding to W and W; since the non-continuous virtual surface array W contains the holes in the entire row and the entire column, the non-continuous virtual cubic array obtained by superimposing W with a mirror image part thereof W contains contiguous holes, corresponds to virtual field tensor of the non-continuous virtual cubic array , and thus contains contiguous missing elements;

(5) designing a translation window of size P_x×P_y×2 to select a sub-tensor of the virtual domain tensor , which contains elements of which indices are (1: P_x−1), (1: P_y−1) and (1:2) respectively in three dimensions of ; then, translating the translation window by one element in turn along the x axial direction and the y axial direction, dividing into L_x×L_y sub-tensors, expressed as , s_x=1, 2, . . . , L_x, s_y=1, 2, . . . , L_y, wherein a value range of the translation window size is as follows:

$2\le P_x\le J_{{\mathbb{Q}}_x}-1,2\le P_y\le J_{{\mathbb{Q}}_y}-1,$

and L_x, L_y, P_x, P_y satisfy the following relationship:

$P_x+L_x-1=J_{{\mathbb{Q}}_x},P_y+L_y-1=J_{{\mathbb{Q}}_y},$

superimposing the sub-tensors with the same index subscript s_y in the fourth dimension to obtain L_y four-dimensional tensors with P_x×P_y×2×L_x dimensions; further, superimposing the L_y four-dimensional tensors in a fifth dimension to obtain a five-dimensional virtual domain tensor ∈^{Px×Py×2×Lx×Ly}, wherein the five-dimensional virtual domain tensor contains spatial angle information in the x axial direction and y axial direction, spatial mirror transformation information, and spatial translation information in the x axial direction and the y axial direction; defining the dimension sets K₁={1, 2}, K₂={3}, K₃={4,5}, and merging through the dimensions of to obtain a three-dimensional reconstructed virtual domain tensor ∈^{PxPy×LxLy×2}:

$K_{\mathbb{Q}}\stackrel{△}{=}{T}_{{\mathbb{Q}}_{\left\{K_1,K_2,K_3\right\}}}$

wherein the three dimensions of respectively represent the spatial angle information, the spatial translation information and the spatial mirror transformation information, thus, the contiguous missing elements in the original virtual domain tensor are randomly distributed to the three spatial dimensions contained in ;

(6) designing a virtual domain tensor filling optimization problem based on tensor kernel norm minimization:

$\underset{{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}}{\min }{{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}}_{*}s.t.P_{\stackrel{\_}{\Omega }}\left({\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}\right)=P_{\stackrel{\_}{\Omega }}\left(K_{\mathbb{Q}}\right),$

wherein the optimization variable ∈^{PxPy×LxLy×2} is the filled virtual domain tensor corresponding to the virtual uniform cubic array , ∥⋅∥_* represents the tensor kernel norm, Ω represents a position index set of the non-missing elements in , P_Ω(⋅) represents the mapping of the tensor on Ω; since the kernel norm is a convex function, the virtual domain tensor filling problem based on the minimization of the tensor kernel norm is a solvable convex optimization problem, and wherein the convex optimization problem is solved to obtain ;

(7) expressing the filled virtual domain tensor as follows:

${\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}=\sum _{k=1}^K{\sigma }_k^2{P}_k\circ q_k\circ c_k,\mathrm{wherein}{p}_k=d_y\left(k\right)\otimes d_x\left(k\right),q_k=g_y\left(k\right)\otimes g_x\left(k\right)\mathrm{are}\mathrm{the}\mathrm{space}\mathrm{factors}\mathrm{of}{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}},d_x\left(k\right)={\left[e^{-j\pi \left(-M_x{N}_x+M_x\right){\mu }_k},e^{-j\pi \left(-M_x{N}_x+M_x+1\right){\mu }_k},\dots ,e^{-j\pi \left(2M_x{N}_x-N_x\right){\mu }_k}\right]}^T,d_y\left(k\right)={\left[e^{-j\pi \left(-M_y{N}_y+M_y\right)v_k},e^{-j\pi \left(-M_y{N}_y+M_y+1\right)v_k},\dots ,e^{-j\pi \left(2M_y{N}_y-N_y\right)v_k}\right]}^T,$

respectively represent the steering vectors of the virtual uniform cubic array along the x axial direction and y axial directions,

$g_x\left(k\right)={\left[1,e^{-j\pi {\mu }_k},\dots ,e^{-j\pi \left(L_x-1\right){\mu }_k}\right]}^T,g_y\left(k\right)={\left[1,e^{-j\pi v_k},\dots ,e^{-j\pi \left(L_y-1\right)v_k}\right]}^T,$

are respectively spatial translation factor vectors corresponding to the x axial direction and y axial direction in the process of intercepting the sub-tensor by the translation window; performing canonical polyadic decomposition on the filled virtual domain tensor to obtain the estimated values of three factor vectors p_k, q_k and c_k, which are expressed as {circumflex over (p)}_k, {circumflex over (q)}_k and ĉ_k; and extracting angle parameters contained in exponential terms of {circumflex over (p)}_k and {circumflex over (q)}_k to obtain a two-dimensional direction-of-arrival estimation result ({circumflex over (θ)}_k, {circumflex over (φ)}_k).

Further, the coprime surface array structure described in step (1) is specifically described as follows: constructing a pair of sparse uniform sub-surface arrays ₁ and ₂ in a plane coordinate system xoy, wherein ₁ contains 2M_x×2M_y antenna elements, and the array element spacings in the x axial direction and the y axial direction are respectively N_xd and N_yd, and position coordinate thereof on xoy is {(N_xdm_x, N_ydm_y), m_x=0, 1, . . . , 2M_x−1,m_y=0,1, . . . , 2M_y−1}; ₂ contains N_x×N_y antenna array elements, the array element spacings in the x axial direction and the y axial direction are respectively M_xd and M_yd, and position coordinate thereof on xoy is {(M_xdn_x, M_ydn_y), n_x=0, 1, . . . , N−1,n_y=0, 1, . . . , N_y−1}; M_x and N_x, and M_y and N_y are respectively a pair of coprime integers; and combining the sub-arrays of ₁ and ₂ in the way that the array elements at (0, 0) position in the coordinate system overlap to obtain a coprime surface array that actually contains 4M_xM_y+N_xN_y−1 physical antenna array elements.

Further, for the cross-correlation tensor derivation described in step (2), in practice, obtaining by calculating the cross-correlation statistic of the tensors and to approximate, that is, sampling cross-correlation tensor ∈^{2Mx×2My×Nx×Ny}:

${\hat{R}}_{{\mathbb{P}}_1{\mathbb{P}}_2}=\frac{1}{T}{\langle X_{{\mathbb{P}}_1},X_{{\mathbb{P}}_2}^{*}\rangle }_3.$

Further, in step (7), performing the canonical polyadic decomposition on the filled virtual domain tensor to obtain the factor vectors {circumflex over (p)}_k, {circumflex over (q)}_k and ĉ_k, and, then extracting the parameters {circumflex over (μ)}_k=sin({circumflex over (φ)}_k)cos({circumflex over (θ)}_k) and {circumflex over (ν)}_k=sin(φ_k)sin(θ_k) from {circumflex over (p)}_k{circumflex over (q)}_k as follows:

${\hat{\mu }}_k=\left(\angle \left(\frac{{\hat{p}}_{k_{\left({\eta }_1+1\right)}}}{{\hat{p}}_{k_{\left({\eta }_1\right)}}}\right)+\angle \left(\frac{{\hat{q}}_{k_{({\eta }_2+10}}}{{\hat{q}}_{k_{\left({\eta }_2\right)}}}\right)\right)/2\pi ,{\hat{v}}_k=\left(\angle \left(\frac{{\hat{p}}_{k_{\left(P_x+{\delta }_1\right)}}}{{\hat{p}}_{k_{\left({\delta }_1\right)}}}\right)+\angle \left(\frac{{\hat{q}}_{k_{\left(P_y+{\delta }_2\right)}}}{{\hat{q}}_{k_{\left({\delta }_2\right)}}}\right)\right)/2\pi ,$

wherein ∠(⋅) represents the operation of taking the argument of a complex number, and a_(a) represents the ath element of a vector a; here, according to the Kronecker structure of {circumflex over (p)}_k and {circumflex over (q)}_k, η∈[1, P_xP_y−1] and η₂∈[1, L_xL_y−1] respectively satisfy mod(η₁, P_x)≠0 and mod(η₂, P_y)≠0, and δ₁∈[1, P_xP_y−P_x], δ₂∈[1, L_xL_y−L_x], mod(⋅) represents the operation of taking a remainder; according to the relationship between the parameter (μ_k, ν_k) and the two-dimensional direction-of-arrival (θ_k, φ_k), obtaining a closed-form solution of the two-dimensional direction-of-arrival estimation ({circumflex over (θ)}_k, {circumflex over (φ)}_k) as follows:

${\hat{\theta }}_k=\mathrm{arc}\tan \left(\frac{{\hat{v}}_k}{{\hat{\mu }}_k}\right),{\hat{\phi }}_k=\sqrt{{\hat{v}}_k^2+{\hat{\mu }}_k^2}.$

Compared with the prior art, the present invention has the following advantages:

(1) The present invention derives the augmented non-continuous virtual surface array based on the cross-correlation tensor, and utilizes the mirror image expansion of the non-continuous virtual surface array to construct the three-dimensional non-continuous virtual cubic array and its corresponding virtual domain tensor, which fully preserves the structure information of all non-contiguous virtual domain statistics of the coprime array;

(2) The present invention proposes a virtual domain tensor filling mechanism for the non-contiguous virtual arrays. The virtual field tensor is reconstructed to disperse its missing elements to meet the low-rank fallibility of the virtual domain tensor. The virtual domain tensor is effectively filled to achieve the high-precision and high-resolution two-dimensional direction-of-arrival estimation.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a general flow block diagram of the present invention.

FIG. 2 is a schematic diagram of the structure of a coprime array constructed by the present invention.

FIG. 3 is a schematic diagram of an augmented non-contiguous virtual surface array derived by the present invention.

FIG. 4 is a schematic diagram of a non-contiguous virtual cubic array derived by the present invention.

FIG. 5 is a performance comparison diagram of the direction-of-arrival estimation accuracy of the method proposed in the present invention under different signal-to-noise ratio conditions.

FIG. 6 is a performance comparison diagram of the direction-of-arrival estimation accuracy of the method proposed in the present invention under different sampling snapshot numbers.

DESCRIPTION OF THE EMBODIMENTS

The technical solution of the present invention will be described in further detail below with reference to the accompanying drawings.

The present invention proposes a two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling in order to solve the problems of loss of multi-dimensional signal structure information and inability to fully utilize virtual domain statistics in existing methods. Effective filling of contiguous missing elements of an original virtual domain tensor is used to realize a Nyquist-matched two-dimensional direction-of-arrival estimation of a coprime surface array. With reference to FIG. 1. the implementation steps of the present invention are as follows:

Step 1: constructing a coprime surface array. Constructing a coprime surface array using 4M_xM_y+N_xN_y−1 physical antenna array elements by a receiving end, as shown in FIG. 2: constructing a pair of sparse uniform sub-surface arrays ₁ and ₂ in a plane coordinate system xoy, wherein ₁ contains 2 M_x×2 M_y antenna elements, and the array element spacings in the x axial direction and the y axial direction are respectively N_xd and N_yd, and position coordinate thereof on xoy is {(N_xdm_x,N_ydm_y),m_x=0, 1, . . . , 2M_x−1,m_y=0, 1, . . . , 2M_y−1}; ₂ contains N_x×N_y antenna array elements, the array element spacings in the x axial direction and the y axial direction are respectively M_xd and M_yd, and position coordinate thereof on xoy is {(M_xdn_x, M_ydn_y),n_x=0, 1, . . . , N_x−1, n_y=0, 1, . . . , N_y−1}; M_x and N_x, and M_y and N_y are respectively a pair of coprime integers; taking the unit interval d as half of the wavelength λ of an incident narrowband signal, that is, d=λ/2; and combining the sub-arrays of ₁ and ₂ in the way that the array elements at (0,0) position in the coordinate system overlap to obtain a coprime surface array that actually contains 4M_xM_y+N_xN_y−1 physical antenna array elements.

Step 2: modeling a tensor of a received signal of the coprime surface array. Assuming that there are K far-field narrowband uncorrelated signal sources from {(θ₁, φ₁), (θ₂, φ₂), . . . , (θ_K, φ_K)} directions, θ_k and φ_k are respectively an azimuth angle and an elevation angle of a kth incident signal source, k=1, 2, . . . , K, superimposing the T sampling snapshot signals of the sparse uniform sub-surface array ₁ in the coprime surface array in the third dimension to obtain a three-dimensional tensor signal ∈^2Mx×2My×T, which is modeled as follows:

$X_{{\mathbb{P}}_1}=\sum _{k=1}^K{a}_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ s_k+N_{{\mathbb{P}}_1},$

wherein s_k=[s_k,1, s_k,2, . . . , s_k,T]^T is a multi-snapshot sampling signal waveform corresponding to the kth incident signal source, [⋅]^T represents a transpose operation, ∘ represents the outer product of a vector, is a noise tensor independent of each signal source, (θ_k, φ_k) and (θ_k, φ_k) are respectively steering vectors of ₁ in an x axial direction and a y axial direction, correspond to a signal source with an incoming wave direction of (θ_k, φ_k), and are expressed as follows:

$a_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi x_{{\mathbb{P}}_1}^{\left(2\right)}{\mu }_k},\dots ,e^{-j\pi x_{{\mathbb{P}}_1}^{\left(2M_x\right)}{\mu }_k}\right]}^T,a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi y_{{\mathbb{P}}_1}^{\left(2\right)}{v}_k},\dots ,e^{-j\pi y_{{\mathbb{P}}_1}^{\left(2M_y\right)}{v}_k}\right]}^T,$

and wherein {, , . . . , }{, , . . . , } represent respectively actual positions of physical antenna elements of the sparse uniform sub-surface array ₁ in the x axial direction and y axial direction, and =0, =0, 82_k=sin(φ_k)cos(θ_k), ν_k=sin(φ_k)sin(θ_k), j=√{square root over (−1)}.

Similarly, the received signals of the sparse uniform sub-surface array ₂ can be expressed by three-dimensional tensor ∈^Nx×Ny×T as follows:

$X_{{\mathbb{P}}_2}=\sum _{k=1}^K{a}_x^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)\circ s_k+N_{{\mathbb{P}}_2},$

wherein is a noise tensor independent of each signal source, (θ_k, φ_k) and (θ_k, φ_k) are respectively steering vectors of ₂ in an x axial direction and a y axial direction, correspond to a signal source with an incoming wave direction of (θ_k, φ_k), and are expressed as follows:

$a_x^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi x_{{\mathbb{P}}_2}^{\left(2\right)}{\mu }_k},\dots ,e^{-j\pi x_{{\mathbb{P}}_2}^{\left(N_x\right)}{\mu }_k}\right]}^T,a_y^{\left({\mathbb{P}}_2\right)}\left({\theta }_k,{\phi }_k\right)={\left[1,e^{-j\pi y_{{\mathbb{P}}_2}^{\left(2\right)}{v}_k},\dots ,e^{-j\pi y_{{\mathbb{P}}_2}^{\left(N_y\right)}{v}_k}\right]}^T,$

and wherein {, , . . . , }{, , . . . , }represent respectively actual positions of physical antenna elements of the sparse uniform sub-surface array ₂ in the x axial direction and y axial direction, and =0, =0.

obtaining a second-order cross-correlation tensor ∈^{2Mx×2My×Nx×Ny} by solving cross-correlation statistic of three-dimensional tensor signals and :

$\begin{array}{c}{R}_{{\mathbb{P}}_1{\mathbb{P}}_2}=E\left[{\langle X_{{\mathbb{P}}_1},X_{{\mathbb{P}}_2}^{*}\rangle }_3\right]\\ =\sum _{k=1}^K{\sigma }_k^2{a}_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\circ a_x^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\circ a_y^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)+\\ N_{{\mathbb{P}}_1{\mathbb{P}}_2}\end{array},$

wherein σ_k²=E[s_ks_k*] represents power of a kth incident signal source, =E[<,>₃] represents the cross-correlation noise tensor, <⋅,⋅>_r represents a tensor contraction operation of the two tensors along a rth dimension, E[⋅] represents a mathematical expectation operation, and (⋅)* represents a conjugation operation. Here, the cross-correlation noise tensor only has an element with a value σ_n² in the (1, 1, 1, 1)th position, wherein σ_n² represents the noise power, and elements in other positions have the same value 0. In practice, obtaining by calculating the cross-correlation statistic of the tensors and to approximate, that is, sampling cross-correlation tensor ∈^{2Mx×2My×Nx×Ny}:

${\hat{R}}_{{\mathbb{P}}_1{\mathbb{P}}_2}=\frac{1}{T}{\langle X_{{\mathbb{P}}_1},X_{{\mathbb{P}}_2}^{*}\rangle }_3;$

Step 3: constructing an augmented non-continuous virtual surface array based on the cross-correlation tensor transformation of the coprime surface array. Since the cross-correlation tensor contains the spatial information corresponding to the two sparse uniform sub-surface arrays and , by merging the dimensions representing the spatial information in the same direction in , the steering vectors corresponding to the two sparse uniform sub-surface arrays can form a difference set array in the exponential term so as to construct a two-dimensional augmented virtual surface array. Specifically, the first and third dimensions of the cross-correlation tensor (represented by the steering vectors (θ_k, φ_k) and

$a_x^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right))$

represent the spatial information of the x axial direction, and the second and fourth dimensions (represented by the steering vectors (θ_k, φ_k) and

$a_y^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right))$

represent the spatial information of the y axial direction; for this reason, defining the dimension sets J₁ ={1, 3}, J₂={2, 4}, and obtaining a virtual domain signal U_W∈^2MxNx×2MyNy by performing the tensor transformation of dimension merging on the cross-correlation tensor :

$U_W\stackrel{△}{=}{R}_{{\mathbb{P}}_1{\mathbb{P}}_{2_{\left\{J_1,J_2\right\}}}}=\sum _{k=1}^K{\sigma }_k^2{b}_x\left(k\right)\circ b_y\left(k\right),$

wherein the

$b_x\left(k\right)=a_x^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\otimes a_x^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)\mathrm{and}{b}_y\left(k\right)=a_y^{{\left({\mathbb{P}}_2\right)}^{*}}\left({\theta }_k,{\phi }_k\right)\otimes a_y^{\left({\mathbb{P}}_1\right)}\left({\theta }_k,{\phi }_k\right)$

are equivalent to the steering vectors of the non-continuous virtual surface array W on the x axis and the y axis, which corresponds to the signal source with the incoming wave direction (θ_k, φ_k), and ⊗ represents the Kronecker product. The non-contiguous virtual surface array W has a size of J_Wx×J_Wy, and contains the holes (that is, missing elements) in an entire row and an entire column, as shown in FIG. 3, J_Wx=3M_xN_x−M_x−N_x+1, J_Wy=3M_yN_y−M_y−N_y+1. Here, in order to simplify the derivation process, the cross-correlation noise tensor is omitted in the theoretical modeling step of U_W. However, in practice, since the sampled cross-correlation tensor is used to replace the theoretical cross-correlation tensor , is still contained in the statistical processing of virtual domain signals;

Step 4: deriving the virtual domain tensor based on the mirror expansion of the non-contiguous virtual surface array. Constructing a virtual surface array W mirroring the non-continuous virtual surface array W about the coordinate axis, and superimposing the W and W on the third dimension to form a three-dimensional non-continuous virtual cubic array of size , as shown in FIG. 4. Here, =J_Wx, =J_Wy, and =2. Correspondingly, arranging the elements in the conjugate transposed signal U_W* of the virtual domain signal U_W to correspond to the positions of the virtual array elements in W, so that the virtual domain signal U_Wcorresponding to the virtual surface array W can be obtained; and superimposing the U_W and U_Won the third dimension, so as to obtain the virtual domain tensor Ucorresponding to the non-contiguous virtual cubic array , which is expressed as:

$U_{\mathbb{Q}}=\sum _{k=1}^K{\sigma }_k^2{\tilde{b}}_x\left(k\right)\circ {\tilde{b}}_y\left(k\right)\circ c_k,$

wherein {tilde over (b)}_x(k) and {tilde over (b)}_y(k) are respectively steering vectors of the non-continuous virtual cubic array on the x axis and the y axis, and correspond to the signal source with an incoming wave direction (θ_k, φ_k); due to existence of the missing elements (holes) in , {tilde over (b)}_x(k) and {tilde over (b)}_y(k) respectively correspond to elements in the hole positions in in the x axial and y axial directions which are set to be zero,

$c_k={\left[1,e^{-j\pi \left(-\left(M_x{M}_y+M_x+M_y\right){\mu }_k-\left(N_x{N}_y+N_x+N_y\right)v_k\right)}\right]}^T$

represents a mirror transformation factor vector corresponding to W and W; since the non-continuous virtual surface array W contains the holes in the entire row and the entire column, the non-continuous virtual cubic array obtained by superimposing W with a mirror image part thereof W contains contiguous holes, corresponds to virtual field tensor of the non-continuous virtual cubic array , and thus contains contiguous missing elements;

Step 5: dispersing contiguous missing elements thereof by virtual domain tensor reconstruction. In order to construct a virtual uniform cubic array to realize the Nyquist matched signal processing, it is necessary to fill the contiguous missing elements in the virtual domain tensor , so as to correspond to a virtual uniform cubic array . However, the low-rank tensor filling technique is premised on the random distribution of missing elements in the tensor, and cannot effectively fill the virtual domain tensor with contiguous missing elements. For this reason, dispersing contiguous missing elements thereof by reconstructing virtual domain tensor . Specifically, designing a translation window of size P_x×P_y×2 to select a sub-tensor of the virtual domain tensor , which contains elements of which indices are (1: P_x−1), (1: P_y−1) and (1:2) respectively in three dimensions of ; then, translating the translation window by one element in turn along the x axial direction and the y axial direction, and dividing into L_x×L_y sub-tensors, expressed as , s_x=1, 2, . . . , L_x, s_y=1, 2, . . . , L_y. A value range of the translation window size is as follows:

$2\le P_x\le J_{{\mathbb{Q}}_x}-1,2\le P_y\le J_{{\mathbb{Q}}_y}-1,$

and L_x, L_y, P_x, P_y satisfy the following relationship:

$P_x+L_x-1=J_{{\mathbb{Q}}_x},P_y+L_y-1=J_{{\mathbb{Q}}_y},$

Superimposing the sub-tensors with the same index subscript s_y in the fourth dimension to obtain L_y four-dimensional tensors with P_x×P_y×2×L_x dimensions; further, superimposing the L_y four-dimensional tensors in a fifth dimension to obtain a five-dimensional virtual domain tensor ∈^{Px×Py×2×Lx×Ly}, wherein the five-dimensional virtual domain tensor contains spatial angle information in the x axial direction and y axial direction, spatial mirror transformation information, and spatial translation information in the x axial direction and the y axial direction; and, merging along the first and second dimensions representing the spatial angle information, and at the same time merging it along the fourth and fifth dimensions representing the spatial translation information, and retaining the third dimension representing the spatial mirror transformation information. Specifically, defining the dimension sets K₁={1,2}, K₂={3}, K₃={4,5} and merging through the dimensions of to obtain a three-dimensional reconstructed virtual domain tensor ∈^{PxPy×LxLy×2}:

$K_{\mathbb{Q}}\stackrel{△}{=}{T}_{{\mathbb{Q}}_{\left\{K_1,K_2,K_3\right\}}},$

wherein the three dimensions of respectively represent the spatial angle information, the spatial translation information and the spatial mirror transformation information, thus, the contiguous missing elements in the virtual domain tensor are randomly distributed to the three spatial dimensions contained in ;

Step 6: performing virtual domain tensor filling based on tensor kernel norm minimization. In order to fill the reconstructed virtual domain tensor , designing a virtual domain tensor filling optimization problem based on tensor kernel norm minimization as follows:

$\underset{{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}}{\min }{{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}}_{*}s.t.P_{\stackrel{\_}{\Omega }}\left({\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}\right)=P_{\stackrel{\_}{\Omega }}\left(K_{\mathbb{Q}}\right),$

wherein the optimization variable ∈^{PxPy×LxLy×2} is the filled virtual domain tensor, corresponding to the virtual uniform cubic array , ∥⋅∥_* represents the tensor kernel norm, Ω represents a position index set of the non-missing elements in , P_Ω(⋅) represents the mapping of the tensor on ∩. Since the kernel norm is a convex function, the virtual domain tensor filling problem based on the minimization of the tensor kernel norm is a solvable convex optimization problem. The convex optimization problem is solved to obtain ;

Step 7: decomposing the filled virtual domain tensor to obtain the direction-of-arrival estimation result. Expressing the filled virtual domain tensor as follows:

${\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}}=\sum _{k=1}^K{\sigma }_k^2{p}_k\circ q_k\circ c_k,\mathrm{wherein}{p}_k=d_y\left(k\right)\otimes d_x\left(k\right),q_k=g_y\left(k\right)\otimes g_x\left(k\right)\mathrm{are}\mathrm{the}\mathrm{space}\mathrm{factors}\mathrm{of}{\stackrel{\_}{K}}_{\stackrel{\_}{\mathbb{Q}}},d_x\left(k\right)={\left[e^{-j\pi \left(-M_x{N}_x+M_x\right){\mu }_k},e^{-j\pi \left(-M_x{N}_x+M_x+1\right){\mu }_k},\dots ,e^{-j\pi \left(2M_x{N}_x-N_x\right){\mu }_k}\right]}^T,d_y\left(k\right)={\left[e^{-j\pi \left(-M_y{N}_y+M_y\right)v_k},e^{-j\pi \left(-M_y{N}_y+M_y+1\right)v_k},\dots ,e^{-j\pi \left(2M_y{N}_y-N_y\right)v_k}\right]}^T,$

respectively represent the steering vectors of the virtual uniform cubic array along the x axial direction and y axial directions,

$g_x\left(k\right)={\left[1,e^{-j\pi {\mu }_k},\dots ,e^{-j\pi \left(L_x-1\right){\mu }_k}\right]}^T,g_y\left(k\right)={\left[1,e^{-j\pi v_k},\dots ,e^{-j\pi \left(L_y-1\right)v_k}\right]}^T,$

are respectively the spatial translation factor vectors corresponding to the x axial direction and y axial direction in the process of intercepting the sub-tensor by the translation window. Performing the canonical polyadic decomposition on the filled virtual domain tensor to obtain estimated values of the factor vectors p_k, q_k and c_k, which are represented as {circumflex over (p)}_k, {circumflex over (q)}_k and ĉ_k, and, then extracting the parameters {circumflex over (μ)}_k=sin({circumflex over (φ)}_k)cos({circumflex over (θ)}_k) and {circumflex over (ν)}_k=sin({circumflex over (φ)}_k)sin({circumflex over (θ)}_k) from {circumflex over (p)}_k and {circumflex over (q)}_k as follows:

${\hat{\mu }}_k=\left(\angle \left(\frac{{\hat{p}}_{k_{\left({\eta }_1+1\right)}}}{{\hat{p}}_{k_{\left({\eta }_1\right)}}}\right)+\angle \left(\frac{{\hat{q}}_{k_{\left({\eta }_2+1\right)}}}{{\hat{q}}_{k_{\left({\eta }_2\right)}}}\right)\right)/2\pi ,{\hat{v}}_k=\left(\angle \left(\frac{{\hat{p}}_{k_{\left(P_x+{\delta }_1\right)}}}{{\hat{p}}_{k_{\left({\delta }_1\right)}}}\right)+\angle \left(\frac{{\hat{q}}_{k_{\left(P_y+{\delta }_2\right)}}}{{\hat{q}}_{k_{\left({\delta }_2\right)}}}\right)\right)/2\pi ,$

wherein ∠(⋅) represents the operation of taking the argument of a complex number, and a_(a) represents the ath element of a vector a; here, according to the Kronecker structure of {circumflex over (p)}_k and {circumflex over (q)}_k, η₁∈[1, P_xP_y−1] and η₂∈[1, L_xL_y−1] respectively satisfy mod(η₁, P_x)≠0 and mod(η₂, P_y)≠0, and δ₁∈[1, P_xP_y−P_x], δ₂∈[1, L_xL_y×L_x], mod(⋅) represents the operation of taking a remainder. According to the relationship between the parameter (μ_k, ν_k) and the two-dimensional direction-of-arrival (θ_k, φ_k), obtaining a closed-form solution of the two-dimensional direction-of-arrival estimation ({circumflex over (θ)}_k, {circumflex over (φ)}_k) as follows:

${\hat{\theta }}_k=\mathrm{arc}\tan \left(\frac{{\hat{v}}_k}{{\hat{\mu }}_k}\right),{\hat{\phi }}_k=\sqrt{{\hat{v}}_k^2+{\hat{\mu }}_k^2}.$

The effect of the present invention will be further described below in conjunction with a simulation example.

Simulation example: the coprime array is used to receive an incident signal, and its parameters are selected as M_x=2, M_y=3, N_x=3, N_y =4, that is, the constructed coprime surface array contains 4M_xM_y+N_xN_y−1=35 physical array elements. The translation window size of the sub-tensor is 6×15×2. Assuming that there are 2 narrowband incident signals, the azimuth and elevation angles of the incident direction are respectively [30.6°, 25.6°] and [40.5°, 50.5°]. Comparing the two-dimensional direction-of-arrival estimation method of coprime surface array based on virtual domain tensor filling proposed by the present invention with the traditional Multiple Signal Classification (MUSIC) method and Tensor Multiple Signal Classification (Tensor MUSIC) method which only utilize the contiguous part of the virtual domain, under the condition of the number of sampling snapshots T=300, plotting performance comparison curves of Root Mean Square Error (RMSE) as a function of signal-to-noise ratio (SNR), as shown in FIG. 5; under the condition of SNR=0 dB, plotting performance comparison curves of RMSE as a function of the number of sampling snapshots T, as shown in FIG. 6.

It can be seen from the comparison results of FIG. 5 and FIG. 6 that no matter in different expected signal-to-noise ratio (SNR) scenarios or in different numbers of the sampling snapshots T scenarios, the method proposed in the present invention has a performance advantage in the direction-of-arrival estimation accuracy. Compared with the traditional MUSIC method, the method proposed in the present invention makes full use of the structural information of the received signal of the coprime surface array by constructing the virtual domain tensor, so as to have better performance of direction-of-arrival estimation; and, compared with the Tensor MUSIC method, the performance advantage of the method proposed in the present invention comes from the use of all non-contiguous virtual domain statistics information through the virtual domain tensor filling, while the Tensor MUSIC method only extracts the continuous part of the non-contiguous virtual array for the virtual domain signal processing, resulting in loss of the virtual domain statistics information.

In summary, the present invention realizes the random distribution of the contiguous missing elements through the virtual domain tensor reconstruction, and based on this, designs the virtual domain tensor filling method based on the tensor kernel norm minimization, and successfully utilizes all the non-continuous virtual domain statistics information, which realizes the high-precision two-dimensional direction-of-arrival estimation of the coprime surface array.

The above descriptions are only preferred embodiments of the present invention. Although the present invention has been disclosed above with preferred embodiments, they are not intended to limit the present invention. Any person skilled in the art, without departing from the scope of the technical solution of the present invention, can make many possible changes and modifications to the technical solution of the present invention by using the methods and technical contents disclosed above, or modify the technical solution of the present invention into equivalent examples. Therefore, any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention without departing from the contents of the technical solution of the present invention still falls within the protection scope of the technical solution of the present invention.

Claims

1. A two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling, wherein the method comprises the following steps: (1) using 4MxMy+NxNy−1 physical antenna array elements by a receiving end, and performing constructing according to a structure of a coprime surface array, wherein Mx, Nx and My, Ny are a pair of coprime integers respectively; decomposing the coprime surface array into two sparse uniform sub-surface arrays 1 and 2, wherein 1 contains 2Mx×2My antenna array elements, array element spacings in an x axial direction and a y axial direction are respectively Nxd and Nyd, 2 includes Nx×Ny antenna array elements, array element spacings in the x axial direction and the y axial direction are respectively Mxd and Myd, and an unit interval d is taken as half of wavelength λ of an incident narrowband signal;(2) if there are K far-field narrowband uncorrelated signal sources from {(θ1, φ1 ), (θ2, φ2), . . . , (θK, φK} directions, θk and φk are respectively an azimuth angle and an elevation angle of a kth incident signal source, k=1, 2, . . . , K, utilizing a three-dimensional tensor ∈Mx×2My'T to express T sampling snapshot signals of a sparse uniform sub-surface array 1 as follows:

2. The two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling according to claim 1, wherein the coprime surface array structure described in the step (1) is specifically described as follows: constructing a pair of sparse uniform sub-surface arrays 1 and 2 in a plane coordinate system xoy, wherein 1 contains 2Mx×2My antenna elements, and the array element spacings in the x axial direction and the y axial direction are respectively Nxd and Nyd , and position coordinate thereof on xoy is {(Nxdmx, Nydmy), mx=0, 1, . . . , 2Mx−1, my=0, 1, . . . , 2My−1 }; 2 contains Nx×Ny antenna array elements, the array element spacings in the x axial direction and the y axial direction are respectively Mxd and Myd, and position coordinate thereof on xoy is {(Mxdnx, Mydny), nx=0, 1, . . . , Nx−1, ny=0, 1, . . . , Ny−1}; Mx and Nx, and My and Ny are respectively a pair of coprime integers; and combining the sub-arrays of 1 and 2 in the way that the array elements at (0, 0) position in the coordinate system overlap to obtain a coprime surface array that actually contains 4MxMy+NxNy−1 physical antenna array elements.

3. The two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling according to claim 1, wherein for the cross-correlation tensor derivation described in the step (2), obtaining by calculating the cross-correlation statistic of the tensor and to approximate, that is, sampling cross-correlation tensor └2Mx×2My×Nx×Ny:

4. The two-dimensional direction-of-arrival estimation method for a coprime surface array based on virtual domain tensor filling according to claim 1, wherein in the step (7), performing the canonical polyadic decomposition on the filled virtual domain tensor to obtain the factor vectors {circumflex over (p)}k, {circumflex over (q)}k and ĉk, and, then extracting the parameters {circumflex over (μ)}k=sin({circumflex over (φ)}k)cos({circumflex over (θ)}k) and {circumflex over (ν)}k=sin(φk)sin(θk) from {circumflex over (p)}k{circumflex over (q)}k as follows: