1.   引言

波达角(Direction Of Arrival, DOA)估计是阵列信号处理中的一个重要研究课题，在无线通信、雷达和医学成像等许多领域中有广泛的应用。在阵列结构方面，平面阵的应用最为广泛。在过去的几十年中，国内外学者已经提出了许多基于平面阵列结构的DOA估计算法，包括面阵^[1,2]、L型阵列^[3]和平行线阵^[4]等。平面阵列通常由几个均匀间隔的线性子阵列构成，具有有限的自由度，例如含$M$个阵元的L型均匀线阵阵列，最多可以估计$M - 1$个入射信号角度。

近年来，稀疏阵列由于能够有效提高阵列的自由度而受到广泛关注，例如最小冗余阵列^[5]、嵌套阵列^[6]和互质阵列^[7]等。与传统阵列相比，稀疏阵列可以在保证性能的前提下充分地减少阵元数，或在阵元数相同的情况下，拥有更大的阵列孔径、更低的旁瓣级，通过对阵元位置和加权的解算改善测向算法的精度、分辨率和自由度。在稀疏嵌套阵列方面，文献[8]和文献[9]将嵌套阵列从1维DOA估计推广到2维DOA估计，提出了一种包含两个均匀间隔线性子阵的2维嵌套阵列。在L型互质阵列方面，有基于迭代最小化和离网格稀疏学习^[10,11]。在平行互质阵列方面，文献[12]首次提出了利用平行互质阵列的互协方差矩阵构造1维向量，通过稀疏重构和最小二乘法得到相互匹配的DOA估计，提升了阵列的自由度，但是在小快拍情况下精度较差，且算法复杂度较高。文献[13]和文献[14]针对复杂度较高问题，将2维DOA估计问题转化为1维表示，利用互协方差矩阵和压缩感知方法进行角度估计，一定程度上提高了估计精度。上述算法主要利用了平行互质阵列的互协方差矩阵，文献[15]和文献[16]使用平行互质阵列的协方差矩阵，采用1维DOA估计结合功率进行匹配实现俯仰角和方位角估计，但是容易出现失配现象。

可以看出，现有平行互质阵列DOA估计算法只利用了阵列的协方差矩阵或互协方差矩阵，需要网格搜索和匹配，存在着计算复杂度较高、算法精度不足、容易出现失配等问题。针对现有算法的不足，本文提出了一种新的基于平行互质线阵2维DOA估计算法，利用两个线阵的协方差矩阵和互协方差矩阵构造新的DOA估计矩阵。在此基础上，利用SVD和ESPRIT算法，根据特征值和特征向量得到相匹配的方位角和俯仰角。与现有算法相比，本文算法充分地利用了自相关和互相关矩阵信息，可以估计更多的信源数，精度更高。同时由于阵元孔径的扩展，算法的分辨能力较高，计算复杂度较低，且在低信噪比和小快拍的情况下性能较好。

文中符号说明：${( \cdot )^{\rm{T}}}$, ${( \cdot )^*}$, ${( \cdot )^{\rm{H}}}$, ${( \cdot )^{ - 1}}$和${( \cdot )^ + }$分别表示矩阵转置、共轭、共轭转置、求逆和求伪逆；${\rm{diag}}(v)$表示以$v$为主对角线元素的对角矩阵；${\rm{vec}}( \cdot )$表示矩阵拉伸；$\otimes$表示Khatri-Rao积；${I}$表示单位矩阵；$\arg ( \cdot )$表示取相位角。

2.   信号模型

本文的阵列模型采用平行互质阵列结构，由两个相互平行的扩展互质阵列^[17]组成，如图1所示，包含子阵1和子阵2。子阵1有$L$个真实物理阵元，由两个不重合均匀线阵交叉构成：阵元间距为$Nd$的线阵，阵元个数为$2M - 1$；阵元间距为$Md$的线阵，阵元个数为$N$。可见，$L{\rm{ = 2}}M - 1 + N$。其中，$M$和$N$是互质的整数，$d$为基本阵元间距，$d = \lambda /2$, $\lambda$为信号的入射波长。这里，用$p = \{ {p_0}, {p_1}, ··· ,{p_{L - 1}}\}$代表其中每个阵元的排列位置。

图 1  平行互质阵列几何模型

下载: 全尺寸图片 幻灯片

在图1中，$\alpha$和$\beta$分别为辐射源与$X$轴和$Y$轴的夹角，满足$\cos \alpha = \sin \varphi \sin \theta$, $\cos \beta = \cos \varphi \sin \theta$。其中，$\varphi$和$\theta$为辐射源的方位角和俯仰角。当有$K$个远场窄带信号入射时，则子阵1第$l$个阵元接收信号模型可表示为

其中，${s_k}(t)$代表第$k$个入射信号的幅度，${n_l}(t)$为噪声。

则整个阵列的接收信号模型为

其中，${{{z}}_1}(t)$为子阵1的接收数据，${{{z}}_2}(t)$为子阵2的接收数据。方向矩阵${{A}}{\rm{ = [}}{{{a}}_1},{{{a}}_2}, ···,{{{a}}_K}{\rm{]}}$的维度为$L \times K$，包含入射信号的方位角和俯仰角信息，响应矢量${{{a}}_k} = [1,{{\rm{e}}^{ - {\rm{j2}}\pi /\lambda {p_1}\cos {\alpha _k}}},{{\rm{e}}^{ - {\rm{j2}}\pi /\lambda {p_2}\cos {\alpha _k}}}, ··· ,$${{\rm{e}}^{ - {\rm{j2}}\pi /\lambda {p_{L - 1}}\cos {\alpha _k}}}]^{\rm{T}}$对应第$k$个入射信号。${{s}}(t)$为$K$个非合作的窄带信号的幅度，${{{n}}_1}(t)$和${{{n}}_2}(t)$表示零均值加性高斯白噪声，且噪声和信号不相关。${{{z}}_2}(t)$中${{\varPhi}} = {\rm{diag}}({{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _1}}},{{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _2}}}, ··· , {{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _K}}})$为$K \times K$的对角矩阵。这里假设所有的入射信号不重合，即${\alpha _1} \ne {\alpha _2} \ne ··· \ne {\alpha _K}$和${\beta _1} \ne {\beta _2} \ne ··· \ne {\beta _K}$。

由于本文采用扩展互质子阵列组成平行阵列，通过虚拟域信号构造，单个子阵物理阵元为$L$($L{\rm{ = 2}}M - 1 + N$)个时，连续虚拟阵元个数为$2MN + 2M - 1$，其自由度可以达到$MN + M{\rm{ - }}1$，突破了传统阵列物理阵元对天线孔径的限制，有效的扩展了天线的孔径和增加了DOA估计的自由度。

3.   基于平行互质虚拟阵列的低复杂度2维DOA估计算法

前文已指出，现有的2维DOA估计算法存在着计算复杂度较高、算法精度不足、容易出现失配等问题。针对这一问题，本文重点研究了2维虚拟阵列，提出了一种适用于平行互质虚拟阵列的2维DOA估计算法。

3.1   扩展矩阵构造

根据接收的入射信号数据，先构造应用于2维虚拟阵列估计的DOA扩展矩阵。子阵1的接收信号${{{z}}_1} = {{{z}}_1}(t)$(对特定时刻$t$)的协方差矩阵为

其中，${{R}} = {\rm{diag}}({\rho _1}{\rm{, }}{\rho _2}, ··· ,{\rho _K})$, ${\rho _k}$$(k = 1, 2, ··· ,K)$为第$k$个入射信号的功率，$\sigma _n^2$代表着噪声的功率，${{I}}$为$K \times K$的单位矩阵。

将${{{R}}_{\rm{b}}}$向量化后可得

其中，${{B}} = [{{{a}}_1} \otimes {{a}}_1^*,{{{a}}_2} \otimes {{a}}_2^*, ··· ,{{{a}}_K} \otimes {{a}}_K^*]$，可以将其看做为虚拟线阵的方向矢量。${{r}} = {[{\rho _1},{\rho _2}, ··· ,{\rho _K}]^{\rm{T}}}$为信号功率，相当于虚拟线阵中的信号，${{{I}}_{{e}}} = {\rm{vec}}({{I}})$，则${{{v}}_1}$可以作为1维非均匀线阵中单次快拍接收的信号。

通过对式(3)中的矩阵进行特征值分解，可以从中估计出噪声功率。去除${{v}}_1^{}$中的噪声影响，可得

接着，考虑到${{{v}}_1}$中存在重复的元素，需要剔除重复元素，并截取连续的虚拟阵元，可以得到$2MN +$$2M - 1$个位于$\{ - (MN + M - 1), - (MN + M - 2),$$··· , - 1,{\rm{ 0, 1}}, ··· ,(MN + M - 2),(MN + M - 1)\}$的阵元位置上的元素，即

其中，${\bar{ B}} \!=\! [{{\bar{ B}}_1}{\rm{, }}{{\bar{ B}}_2}, ··· ,{{\bar{ B}}_K}]$, ${{\bar{ B}}_k} = [{{\rm{e}}^{{\rm{j}}\pi (MN + M - 1)\cos {\alpha _k}}},$${{\rm{e}}^{{\rm{j}}\pi \left( {MN \!+\! M \!-\! 2} \right)\!\cos\! {\alpha _k}}}\!, ··· ,\!0,{{\rm{e}}^{ - {\rm{j}}\pi \!\cos \!{\alpha _k}}}\!, ···s ,\!{{\rm{e}}^{ - {\rm{j}}\pi (MN \!+\! M \!-\! 1)\!\cos\! {\alpha _k}}}]$。${\bar{ B}}$为$(2MN + 2M - 1) \times K$的矩阵，${{r}}$为$K \times 1$的矩阵。

通过对向量${{\bar{ v}}_1}$中的元素进行重新排列可得

可见，矩阵${{{V}}_1}$的尺寸为$(MN + M) \times (MN + M)$。可以将${{{V}}_1}$表示为

其中，${{W}}$为对角矩阵，和信号的功率相关。${{D}}$中包含了入射信号与X轴的夹角信息。

同样，可以通过接收的信号，定义子阵1和子阵2的互协方差矩阵为

其中，${{\varPhi}} = {\rm{diag}}({{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _1}}},{{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _2}}},··· ,{{\rm{e}}^{ - {\rm{j}}\pi \cos {\beta _K}}})$，含有入射信号与Y 轴的夹角信息。由于两个子阵的噪声不相关，因此互协方差矩阵没有噪声的干扰项。对${{{R}}_{\rm{c}}}$进行向量化可得

同理，剔除${{v}}_2^{}$中的重复元素，并取其中连续的虚拟阵元，可以得到

与式(7)类似，构造互协方差矩阵，即对${{\bar{ v}}_2}$中的元素重新排列，可得

由式(8)和式(12)，构造DOA估计的扩展矩阵

对${{{R}}_{\rm{m}}}$分析可知，${\bar{ D}}$中包含了入射信号与X轴和Y轴的夹角信息所有信息，即${\bar{ D}}$为所求矩阵。通过对${{{R}}_{\rm{m}}}$进行奇异值分解，得到的信号子空间，即为${\bar{ D}}$张成的空间。

3.2   2维DOA估计

在3.1节中，通过子阵1的协方差矩阵及其与子阵2的互协方差矩阵，构造出DOA估计矩阵${{{R}}_{\rm{m}}}$。接下来，本文利用SVD和ESPRIT算法，由特征值和特征向量得到相匹配的方位角和俯仰角。

首先，对DOA估计矩阵${{{R}}_{\rm{m}}}$进行奇异值分解可得

其中，${{{U}}_1}$为信号子空间，${{{U}}_2}$为噪声子空间，则

${{{U}}_1}$为$2\left( {MN + M} \right) \times K$的矩阵，$K$为入射信号的个数。由于${{{U}}_1}$为信号子空间，这里存在一个$K \times K$的矩阵${{T}}$使得式(17)成立

由上述分析可知，信号子空间中包含着入射信号的方位角和俯仰角。其中，${{D}}$为方向矩阵，包含入射信号与$X$轴夹角信息$\alpha$; ${{\varPhi}}$为对角矩阵，包含入射信号与$Y$轴夹角信息$\beta$。则有

其中，${{{U}}_{11}}$和${{{U}}_{12}}$的维度均为$\left( {MN + M} \right) \times K$。

通过${{{U}}_{11}}$和${{{U}}_{12}}$的关系，构造矩阵${{F}}$，使得

可以得到${{\varPhi}}$中对应的特征值${\rm{\{ }}{\psi _k},k = 1,2, ··· ,K{\rm{\} }}$，对应入射的$K$个信号与Y 轴夹角信息

通过对式(20)进行特征值分解，其对应的特征向量构成的矩阵为${{T}}$，通过(18)式有

其中，$\hat {{D}}$为$(MN + M) \times K$的矩阵。

在文献[2]中，通过构造一个谱峰搜索函数进行角度搜索，进而得到入射信号与X轴的夹角，算法复杂度高，效率较低。本文借鉴旋转不变思想，将矩阵$\hat {{D}}$进行分块，${{{C}}_1}$取矩阵$\hat {{D}}$第1行到$MN + M - 1$行，${{{C}}_2}$取矩阵$\hat { {D}}$的第2行到$MN + M$行，则有

其中，${{\varPsi}} {\rm{ = diag(}}{{\rm{e}}^{{\rm{j}}\pi \cos {\alpha _1}}},{{\rm{e}}^{{\rm{j}}\pi \cos {\alpha _2}}}, ···,{{\rm{e}}^{{\rm{j}}\pi \cos {\alpha _K}}}{\rm{)}}$，为旋转不变因子，则有

通过对${{\varPsi}}$进行特征值分解，得到特征值${\gamma _k}$，则有

通过联立式(21)和式(25)，可得到相互匹配的方位角和俯仰角，即

3.3   算法步骤

基于以上理论分析，这里给出本文算法的具体步骤。

步骤 1　根据实际应用，通过有限的快拍数对协方差矩阵和互协方差矩阵进行估计，即

其中，$P$为接收信号的快拍数。

步骤 2　由${{\hat{ R}}_{\rm{b}}}$和${{\hat{ R}}_{\rm{c}}}$构造新的DOA估计矩阵${{{R}}_{\rm{m}}}$。

步骤 3　对矩阵${{{R}}_{\rm{m}}}$进行SVD，得到式(14)所示的信号子空间${{{U}}_1}$。

步骤 4　构造角度估计矩阵${{F}}$，求出辐射源角度信息${\beta _k}$。

步骤 5　借用ESPRIT思想，构造包含方位角信息的旋转因子${{\varPsi}}$，估计出辐射源角度信息${\alpha _k}$。

步骤 6　由式(26)和式(27)，得到入射信号的方位角和俯仰角。

上述算法步骤主要包括协方差和互协方差矩阵构建、EVD和SVD操作。经推导，算法复杂度约为$O\{ 2{L^2}P + 2{K^3} + {(2MN + 2M)^3}\}$，其中$L{\rm{ = 2}}M -$$1 + N$, P为采样快拍数，M和N为互质数，K为入射信号个数。

本文所提出的算法，利用互质虚拟阵列的思想，对于2维平行阵列进行扩展，能利用有限的阵元估计更多的入射角度，同时由于阵元孔径的扩展，对入射角度的分辨能力更好，并且在低信噪比和小快拍的情况下精确度相对较高，且算法复杂度更低。

4.   实验仿真

在本节实验中，令$M$=3, $N$=5，则子阵1的真实阵元数为$2M - 1 + N$，即L=10，阵元位置为${\rm{\{ }}0{\rm{\;3\;5\;6\;9\;10\;12\;15\;20\;25\} }}$。考虑多个入射信号的情况，假设信号源个数$K{\rm{ = 11}}$，即入射信号大于实际阵元数。此时，传统平行线阵算法均已失效，而本文算法仍然可以有效地估计出入射角度，如图2所示。本文提出算法能很好地扩展阵列孔径，提高阵元利用率，最多可估计$MN + M - 1$个入射角度，远高于传统算法。

图 2  K=11时算法估计结果(SNR=10 dB, P=500)

下载: 全尺寸图片 幻灯片

为了进一步验证提出算法的高分辨性能，我们假设入射信号个数为2，入射角为$({10^\circ },{11^\circ })$和$({11^\circ },{12^\circ })$。此时，传统算法失效。本文算法的实验结果如图3所示，可以看出，虽然存在一定误差，但依然可以分辨出两个非常接近的入射角度。需要注意的是，由于角度较为接近，实验中需要较高的信噪比和快拍数。

图 3  高分辨率实验(K=2, SNR=20 dB, P=500)

下载: 全尺寸图片 幻灯片

为了衡量算法的估计精度，本文采用均方根误差(RMSE)准则，定义为

其中，$K$代表着入射信号的个数，$Q$代表蒙克卡罗实验次数，${\hat \varphi _{k,q}}$和${\hat \theta _{k,q}}$代表第$q$次实验的DOA估计值。

图4为快拍数对本文算法性能的影响，并与已有算法进行比较。其中，文献[4]基于平行线阵2维PM算法，文献[12]基于基于平行互质阵列算法。考虑现有算法的适用范围，这里对4个信源进行估计，入射角分别为$({10^\circ },{10^\circ })$, $({20^\circ },{20^\circ })$, $({30^\circ },{30^\circ })$和$({40^\circ },{40^\circ })$。可以看出，本文算法在少量快拍下也有较好的估计性能。特别当快拍数P=10时候，本文算法依然可以实现有效DOA估计，可适用于小快拍场景下工作。

图 4  不同快拍数算法性能对比(K=4, SNR=10 dB)

下载: 全尺寸图片 幻灯片

接下来，分析信噪比变化对图4中各种算法性能的影响。为了减小快拍数对不同算法的影响，实验中采用较大的快拍数，即P=1000。实验中入射信号情况与图4相同，实验结果如图5所示。可以看出，本文算法在低信噪比情况下依然具有较高的估计精度。

图 5  不同信噪比下的性能分析(K=4, P=1000)

下载: 全尺寸图片 幻灯片

最后，为了对不同DOA估计算法的复杂度进行比较分析，在相同硬件和软件条件下，进行1000次蒙特卡罗实验，统计各种算法的运行时间。实验中计算机的CPU为I7-8550U，内存为8 G。实验中入射信号情况与图3和图4相同，采样快拍数P = 500，信噪比SNR=20 dB。实验统计结果如表1所示，可以看出，本文算法的计算复杂度优于文献[12]中的平行互质阵列算法，但略高于平行线阵PM算法。

表 1  不同2维DOA估计算法运行时间(s)

估计算法	本文算法	文献[12]	文献[4]
运行时间	2.283	26.413	1.268

下载: 导出CSV 

| 显示表格

5.   结论

本文提出一种基于平行互质虚拟阵列的低复杂度2维DOA估计算法，将传统平行线阵与互质虚拟阵列相结合。本文采用扩展的DOA矩阵，在进行DOA估计时，通过奇异值分解并提取旋转不变因子，避免了传统算法的谱峰搜索，降低了算法复杂度，获得了自动匹配的DOA估计。同时，采用虚拟阵列扩展了阵元孔径，解决了传统DOA估计算法入射信源数小于真实物理阵元数问题。仿真实验结果表明，本文所提算法具有更高的分辨能力，可以分辨更多的辐射源信号，而且在低信噪比和小快拍情况下也优于传统DOA估计算法。