1.   引言

到达角(Direction-Of-Arrival, DOA)估计是信号处理的一个重要分支，广泛应用于雷达目标探测、语音通信、无源定位等领域^[1-4]。在实际应用中，通常需要通过2维DOA估计获得辐射源的方位角和俯仰角信息。不同的阵列结构会对DOA估计性能产生不同的影响，其中L型阵列结构简单，受到许多研究者的关注。

最初的2维DOA估计算法大多基于2维谱峰搜索^[5]，计算量大且估计精度不高。文献[6]利用ESPRIT思想分别进行方位角和俯仰角估计，存在角度配对问题。文献[7]通过解析优化获得方位角与仰角间的关系，避免了复杂的计算，但估计精度有所下降。文献[8]将三线性分解运用到2维DOA估计中，该算法无需谱峰搜索，且能够实现参数自动配对，但其中的交替迭代过程仍具有较大的计算量。文献[9]提出了一种自动配对的DOA估计算法，在保证精度的前提下进一步降低了运算复杂度。

然而，当存在相干信号时，上述方法性能出现恶化。针对相干信号，文献[10]将空间平滑(Spatial Smoothing, SS)方法扩展到2维阵列，但是需要复杂的2维谱峰搜索。文献[11]提出了一种互相关矩阵ESPRIT(Cross Correlation Matrix ESPRIT, CCM-ESPRIT)算法，将互相关矩阵划分成多个子阵来消除信号的相干性，但是该方法孔径损失较大。文献[12]在此基础上提出了互相关矩阵相关矢量(Cross Correlation Matrix Correlation Vector Method, CCM-CVM)算法，减少了孔径损失，但也降低了估计精度，该文献同时提出了一种名为CCM-MCVM的改进算法，略微提高了估计精度，但增加了算法复杂度。文献[11]和文献[12]都使用了同一种角度自动配对算法，这种配对算法在多个相干信源条件下极易导致配对错误。文献[13]将稀疏重构运用到L型阵列DOA估计中，利用正交匹配追踪(OMP)法分别估计独立信号和相干信号，有效减小了孔径损失问题，但同时也引入了网格失配误差。

文献[14]针对相干信号DOA估计提出了一种主奇异矢量分析法(Principal-singular-vector Utilization for Modal Analysis, PUMA)，具有估计精度高、复杂度低的优点。受文献[14]启发，本文将PUMA算法扩展到L型阵列，提出了L-PUMA算法及其改进算法L-MPUMA，进一步提高了估计精度，减少了孔径损失。本文使用${\left( \cdot \right)^{\rm{T}}}$, ${\left( \cdot \right)^{\rm{H}}}$, ${\left( \cdot \right)^{ - 1}}$和${\left( \cdot \right)^{\rm{ + }}}$分别表示转置、共轭转置、逆运算和伪逆运算，${\rm{E}}\left\{ \cdot \right\}$, ${\rm{vec}}\left( \cdot \right)$, ${\rm{tr}}\left( \cdot \right)$, ${\rm{sum}}\left( \cdot \right)$, ${\rm{diag}}\left( \cdot \right)$和${\rm{blkdiag}}\left( \cdot \right)$分别表示期望、矢量化操作、矩阵的迹、求和、对角化和分块对角化操作，${{{I}}_m}$表示$m \times m$维单位矩阵，${{{{\textit{0}}}}_{m \times n}}$表示$m \times n$维零矩阵，${{{J}}_m}$表示$m \times m$维置换矩阵，$\otimes$表示Kronecker积。

2.   信号模型

假设接收阵列是由$x - z$平面内沿$x$轴和$z$轴的$M$元均匀线阵组成的L型阵列，共用原点位置的参考阵元${z_0}$，同一轴上相邻阵元间距$d$等于半波$\lambda /2$，总共有阵元$2M - 1$个，阵列结构如图1所示。假设3维空间内有$K$个窄带信号入射到L型阵列上，满足$K < M$。前${N_u}$个信号是独立信号，后$K - {N_u}$个相干信号分为${N_c}$组，组间信号相互独立，独立信号与相干信号组数之和$N = {N_u} + {N_c}$，第$k$个信号的仰角为${\theta _k}$，方位角为${\phi _k}$。

图 1  L型阵列几何结构

下载: 全尺寸图片 幻灯片

$x$轴和$z$轴的接收数据可分别表示为

式中${{S}}(t) = {\left[ {{s_1}(t),{s_2}(t), ··· ,{s_K}(t)} \right]^{\rm{T}}} \in {{\mathbb{C}}^{K \times 1}}$为入射信号矢量，${{{N}}_1}(t)$和${{{N}}_{_2}}(t)$表示$M$维零均值高斯白噪声，噪声功率为${\sigma ^2}$，且相互独立。${{{A}}_{{x}}} = \left[ {{a}}({\theta _1}, {\phi _1}),{{a}}({\theta _2},{\phi _2}), ··· ,{{a}}({\theta _K},{\phi _K}) \right] \in {{\mathbb{C}}^{M \times K}}$为$x$轴阵列流形，${{{A}}_{{z}}} = \left[ {{{a}}({\theta _1}),{{a}}({\theta _2}), ··· ,{{a}}({\theta _K})} \right] \in {{\mathbb{C}}^{M \times K}}$为$z$轴阵列流形。其中${{a}}({\theta _i},{\phi _i}) = \left[ 1,{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\sin ({\theta _i}){\rm{cos}}({\phi _i})/\lambda }}, ··· , {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}\left( {M - 1} \right)d\sin ({\theta _i}){\rm{cos}}({\phi _i})/\lambda }} \right]^{\rm{T}}$, ${{a}}({\theta _i}) = \left[ 1,{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\cos ({\theta _i})/\lambda }}, ··· ,{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}\left( {M - 1} \right)d\cos ({\theta _i})/\lambda }} \right]^{\rm{T}}$。

3.   L-PUMA算法

L-PUMA算法分为3个步骤，首先是降噪处理，然后是仰角和方位角估计，最后是角度配对。

3.1   互协方差矩阵降噪

在L型阵列中，除了原点处的噪声，$x$轴和$z$轴的噪声矢量是相互独立的，其噪声协方差为0，而$x$轴和$z$轴又共用原点处的参考阵元，因此原点位置的噪声协方差为${\sigma ^2}$。噪声协方差矩阵可表示为

易得，接收数据的互协方差矩阵为

式中${\varOmega _{{\rm{S}}1}}$为大奇异值，${\varOmega _{{\rm{n}}1}}$为小奇异值，${{{U}}_{{\rm{S}}1}}$, ${{{V}}_{{\rm{S}}1}}$分别为大奇异值对应的左、右奇异矢量矩阵，${{{U}}_{{\rm{n}}1}}$, ${{{V}}_{{\rm{n}}1}}$分别为小奇异值对应的左、右奇异矢量矩阵。同理，接收数据的自协方差矩阵

式中${\varOmega _{{\rm{S}}2}}$为大特征值，${\varOmega _{{\rm{n}}2}}$为小特征值，${{{U}}_{{\rm{S}}2}}$和${{{U}}_{{\rm{n}}2}}$分别为大特征值和小特征值对应的特征矢量矩阵。

许多文献选择使用互协方差来做降噪处理^[11,12]，但没有考虑到互协方差共同参考阵元而引入的噪声矩阵${{{Q}}_{\rm{0}}}$，实际上${{{Q}}_{\rm{0}}}$会比${\sigma ^2}{{{I}}_M}$带来更大的估计误差。将式(5)对应的自协方差矩阵左乘${{{U}}^{\rm{H}}_{{\rm{n}}2}}$，右乘${{{U}}_{{\rm{n}}2}}$，有

同时减去${\sigma ^2}{{{I}}_{M - K}}$，有${{{U}}^{\rm{H}}_{{\rm{n}}2}}{{{A}}_{{z}}}{{{R}}_{\rm{S}}}{{{A}}^{\rm{H}}_{{z}}}{{{U}}_{{\rm{n}}2}} = {{0}}$，即${{{U}}^{\rm{H}}_{{\rm{n}}2}}{{{A}}_{{z}}} = 0$。

同理，将式(4)对应的互协方差矩阵左乘${{{U}}^{\rm{H}}_{{\rm{n}}1}}$，右乘${{{V}}_{{\rm{n}}1}}$，有

将式(7)第1和第3行相减，有${{{U}}^{\rm{H}}_{{\rm{n}}1}}{{{A}}_{{z}}}{{{R}}_{\rm{S}}}{{{A}}^{\rm{H}}_{{x}}}{{{V}}_{{\rm{n}}1}} = {\varOmega _{{\rm{n}}1}} - {{{U}}^{\rm{H}}_{{\rm{n}}1}}{{{Q}}_{\rm{0}}}{{{V}}_{{\rm{n}}1}} \ne 0$，所以有${{{U}}^{\rm{H}}_{{\rm{n}}1}}{{{A}}_{{z}}} \ne 0$。可知阵列流形${{{A}}_{{z}}}$与噪声子空间${{{U}}_{{\rm{n}}2}}$正交，而与噪声子空间${{{U}}_{{\rm{n}}1}}$不正交。

在实际情况中，快拍数$L$为有限值，此时互协方差矩阵近似表示为

式中，$L$为快拍数，${{Z}}$和${{X}}$均为$M \times L$维矢量，${{{Q}}_{\rm{1}}}\left( {{\sigma ^2},L} \right)$是${\sigma ^2}$和$L$的函数，包含了有限快拍产生的误差以及信号噪声间的互相关成分。同理，自协方差矩阵近似表示为

式中${{{Q}}_{\rm{2}}}\left( {{\sigma ^2},L} \right)$包含了相同的噪声成分。易知，当快拍数有限时，噪声子空间${{{U}}_{{\rm{n}}1}}$, ${{{U}}_{{\rm{n}}2}}$与阵列流形${{{A}}_{{z}}}$均不正交。若选取互协方差矩阵进行后续的DOA估计，则需要减少误差的影响，其中${{{Q}}_{\rm{1}}}$和${{{Q}}_{\rm{2}}}$为随机误差，无法消除，而${{{Q}}_{\rm{0}}}$只需要估计出噪声功率即可进行消除。

根据式(5)对${{{R}}_Z}$进行特征值分解得到噪声功率的估计值${\hat \sigma ^2} = {\rm{sum}}\left( {{\rm{diag}}\left( {{\varOmega _{{\rm{n}}2}}} \right)} \right)/\left( {M - N} \right)$，最后可得降噪后的互协方差矩阵${{R}} = {{\hat{ R}}_{Z,X}} - {{\hat{ Q}}_{\rm{0}}}$，其中${{\hat{ Q}}_{\rm{0}}}$为${{{Q}}_{\rm{0}}}$的估计值。

3.2   仰角和方位角估计

L-PUMA算法需要预先知道相干信号的组数$N$和总数$K$，可以根据文献[15]所提的改进盖氏圆算法获得，在估计出$N$和$K$之后，进行仰角和方位角的估计。对${{R}}$进行奇异值分解，得到

${\varOmega _{\rm{S}}}$为$N$个大奇异值，${\varOmega _{\rm{n}}}$为$M - N$个小奇异值，${\varOmega _{\rm{S}}} = {\rm{diag}}\left( {{\lambda _1},{\lambda _2}, ··· ,{\lambda _N}} \right)$, ${\varOmega _{\rm{n}}} = {\rm{diag}}\left( {{\lambda _{N + 1}},{\lambda _{N + 2}}, ···, {\lambda _M}} \right)$。${{{U}}_{\rm{S}}}$和${{{V}}_{\rm{S}}}$为大奇异值对应的左、右主奇异矢量矩阵，${{{U}}_{\rm{n}}}$和${{{V}}_{\rm{n}}}$为小奇异值对应的噪声奇异矢量矩阵。不同于非相干条件，此时界定信号子空间和噪声子空间的不再是信号个数$K$，而是信号组数$N$。由文献[16]可知，相干条件下信号子空间和阵列流形间的关系可以表示为

式中，${{{G}}_1}$, ${{{G}}_2}$均为$K \times N$维变换矩阵，其秩为$N$。

先对仰角$\theta$进行估计，构造以${{z_i} = {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\cos ({\theta _i})/\lambda }},} \;\; {i = 1,2, ··· ,K}$为根的$K$次方程

将式(13)展开成

不妨令${b_0} = 1$，只需要确定多项式系数${{b}} = \left[ {b_1},{b_2}, ··· , {b_K} \right]^{\rm{T}}$，即可通过方程求根得到信号的仰角。构造$M \times \left( {M - K} \right)$维Toeplitz矩阵${{B}}$

对于$i = 1,2, \ ··· \ ,K$，有

因此，${{{B}}^{\rm{H}}}{{{A}}_{{z}}} = {{{\textit{0}}}}_{\left( {M - K} \right) \times K}$，可知${{B}}$是阵列流形${{{A}}_{{z}}}$的正交补。又因为 ${{{U}}_{\rm{S}}} = {{{A}}_{{z}}}{{{G}}_1}$，所以有

将${{{U}}_{\rm{S}}}$写成${{{U}}_{\rm{S}}} = \left[ {{{{u}}_1},{{{u}}_2}, ··· ,{{{u}}_N}} \right]$，其中${{{u}}_i}$为第$i$个主奇异矢量，有

式中

${\left[ {{{{u}}_i}} \right]_K}$表示矢量${{{u}}_i}$的第$K$个元素。对式(17)进行矢量化

式中${{F}} = {\left[ {{{{F}}_1}^{\rm{T}},{{{F}}_2}^{\rm{T}},··· ,{{{F}}_N}^{\rm{T}}} \right]^{\rm{T}}}$, ${{g}} = \left[ {{{g}}_1}^{\rm{T}},{{{g}}_2}^{\rm{T}}, ··· , {{{g}}_N}^{\rm{T}} \right]^{\rm{T}}$。在误差条件下，式(21)可以表示为

${\hat{ F}}$和${\hat{ g}}$是${{F}}$和${{g}}$的测量值。用加权最小二乘法(Weighted Least Squares, WLS)来求解${{b}}$，即

式中

${{\hat{ u}}_{\rm{S}}} = {\rm{vec}}\left( {{{{\hat{ U}}}_{\rm{S}}}} \right)$为矢量化后的主特征矢量。容易得到式(23)的解为

然而${{W}}$又是${{b}}$的函数，因此不能直接用式(26)求解${{b}}$，需要对${{W}}$作简化。文献[14]给出了${{W}}$的近似解

式中

至此${\hat{ b}}$可以通过交替迭代求得。首先通过式(22)的最小二乘解${{\hat{ b}}_0} = {{\hat{ F}}^{\dagger} }{\hat{ g}}$确定${\hat{ b}}$的初始值，再代入式(27)初始化${\hat{ W}}$，然后将${\hat{ W}}$代入式(23)更新${\hat{ b}}$，最后求解线性方程即可得到信号的仰角估计值$\left\{ {{{\hat \theta }_i}} \right\}_{i = 1}^K$。

同理，对方位角$\phi$进行估计，构造以${x_i} = {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\sin ({\theta _i}){\rm{cos}}({\phi _i})/\lambda }},\;\;\;{i = 1,2, ··· ,K}$为根的$K$次方程

只需要确定多项式系数${{c}} = {\left[ {{c_1},{c_2}, ··· ,{c_K}} \right]^{\rm{T}}}$，即可得到信号的方位角。构造$M \times \left( {M - K} \right)$维Toeplitz矩阵

易知${{{C}}^{\rm{H}}}{{{V}}_{\rm{S}}} = {{{0}}_{\left( {M - K} \right) \times K}}$，因此只需要将式(22)中的${{{U}}_{\rm{S}}}$替换成${{{V}}_{\rm{S}}}$，再进行相同的操作，即可确定系数${\hat{ c}}$，最后求解以${\hat{ c}}$为系数的方程得到$K$个根$\left\{ {{x_i}} \right\}_{i = 1}^K$。需要注意的是，不能通过${x_i}$直接估计方位角$\left\{ {{{\hat \phi }_i}} \right\}_{i = 1}^K$，因为${x_i} = {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\sin ({\theta _i}){\rm{cos}}({\phi _i})/\lambda }}$，与${\hat \theta _i}$之间存在配对问题，需要通过配对算法得到正确的配对次序。

3.3   角度配对

通过仰角估计值$\left\{ {{{\hat \theta }_i}} \right\}_{i = 1}^K$可以得到${{\hat{ A}}_{{z}}}$，通过$K$个根$\left\{ {{x_i}} \right\}_{i = 1}^K$可以得到${{\hat{ A}}_{{x}}}$。然而在上述算法中，仰角和方位角是分别估计得到的，不存在关联，即${{\hat{ A}}_{{z}}}$和${{\hat{ A}}_{{x}}}$之间存在顺序模糊。易知自协方差矩阵与互协方差矩阵中都含有信号协方差${{{R}}_{\rm{S}}}$，通过此共有信息可实现角度配对。将${{\hat{ A}}_{{z}}}$和${{\hat{ A}}_{{x}}}$分别表示为

${{{E}}_1}$和${{{E}}_2}$为列变换矩阵。易知，去噪后的自协方差矩阵${{{\bar R}}_Z}$,${{{\bar R}}_X}$以及互协方差矩阵${{R}}$中都含有信号协方差矩阵${{{R}}_{\rm{S}}}$，其中${{{\bar R}}_Z} = {{Z}}{{{Z}}^{\rm{H}}}/L - {\hat \sigma ^2}{{{I}}_M}$, ${{{\bar R}}_X} = {{X}}{{{X}}^{\rm{H}}}/L - {\hat \sigma ^2}{{{I}}_M}$。通过式(34)求出 ${{{R}}_{\rm{S}}}$

最小化如式(35)代价函数$p$，可以得到${{{E}}_1}$, ${{{E}}_2}$

不妨令${{H}} = {{{E}}_2}^{\rm{T}}{{{E}}_1}$，其中${{H}}$为列变换矩阵，在求解过程中需要加入本身的结构约束，即${{H}}$的每行只有一个元素为1，其余元素为0。式(35)可进一步简化为

由式(32)和式(33)可知，${{\hat{ A}}_{{z}}}$与${{\hat{ A}}_{{x}}}{{H}}$是完全匹配的。因此只需最小化代价函数$p$，求出变换矩阵${{H}}$即可完成角度配对。相比于文献[17]中的PSCM配对算法，本文所提的配对算法精度更高。

3.4   L-PUMA算法步骤

步骤 1　首先求出降噪后的互协方差矩阵${{R}}$；

步骤 2　将${{{U}}_{\rm{S}}}$按照式(19)和式(20)的形式构造矩阵${{F}}$和${{g}}$，得到${{b}}$的初始值${{\hat{ b}}_0} = {{\hat{ F}}^{\rm{ + }}}{\hat{ g}}$，再通过式(27)和式(26)的交替迭代优化${{b}}$，然后求解线性方程式(14)计算出$K$个根$\left\{ {{z_i}} \right\}_{i = 1}^K$，得到仰角的估计值${\hat \theta _i} = \arccos \left( { - \dfrac{\lambda }{{2{\rm{\pi }}d}}\angle {z_i}} \right)$；

步骤 3　将步骤2中的${{{U}}_{\rm{S}}}$替换为${{{V}}_{\rm{S}}}$，重复步骤2的迭代操作，得到${{c}}$，再求解式(30)计算出$K$个根$\left\{ {{x_i}} \right\}_{i = 1}^K$, ${x_i}$中同时包含了仰角和方位角的信息；

步骤 4　利用$\left\{ {{z_i}} \right\}_{i = 1}^K$和$\left\{ {{x_i}} \right\}_{i = 1}^K$得到估计的阵列流形${{\hat{ A}}_{{z}}}$和${{\hat{ A}}_{{x}}}$，再通过最小化式(36)求出变换矩阵${{H}}$，最后通过${{H}}$得到配对后的角度$\left\{ {{{\hat \theta }_i}} \right\}_{i = 1}^K$, $\left\{ {{{\hat \phi }_i}} \right\}_{i = 1}^K$。

4.   改进的L-PUMA(L-MPUMA)算法

L-PUMA算法的估计精度由线性方程${\hat{ F}}{{b}} \approx {\hat{ g}}$决定，${\hat{ F}}$的维数越大，方程的个数越多，则${{b}}$的计算精度越高。当信号完全相干时，主奇异矢量的个数减至1，大幅降低了线性方程的个数，算法性能下降严重，L-MPUMA通过构造反向主奇异矢量${{{U}}_{\rm{R}}} = {{{J}}_M}{{{U}}_{\rm{S}}}^*$来进一步扩大${\hat{ F}}$的维数，易知

式中$\phi = {\rm{diag}}\left( \left[ {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\cos ({\theta _1})/\lambda }},{{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\cos ({\theta _2})/\lambda }}, ··· , {{\rm{e}}^{ - {\rm{j}}2{\rm{\pi }}d\cos ({\theta _K})/\lambda }} \right] \right)$，令${{{U}}_{{\rm{SR}}}} = \left[ {{{{U}}_{\rm{S}}},{{{U}}_{\rm{R}}}} \right]$为增广主奇异矢量，可知${{B}}$与${{{U}}_{{\rm{SR}}}}$正交，因此式(22)可扩展为

式中${{\hat{ F}}_{{\rm{SR}}}} = {\left[ {{{{\hat{ F}}}^{\rm{T}}},{{\left( {{{{J}}_{M - K}}{{{\hat{ F}}}^*}} \right)}^{\rm{T}}}} \right]^{\rm{T}}}$, ${{\hat{ g}}_{{\rm{SR}}}} =\Bigr[ {{{\hat{ g}}}^{\rm{T}}}, {{\left( {{{{J}}_{M - K}}{{{\hat{ g}}}^*}} \right)}^{\rm{T}}}\Bigr]^{\rm{T}}$。式(38)的解为

式中${{\hat{ e}}_{{\rm{SR}}}} = \left( {{{{I}}_{2N}} \otimes {{{B}}^{\rm{H}}}} \right){{\hat{ u}}_{{\rm{SR}}}}$, ${{{W}}_{{\rm{SR}}}} = {\left( {{\rm{E}}\left\{ {{{{\hat{ e}}}_{{\rm{SR}}}}{{{\hat{ e}}}_{{\rm{SR}}}}^{\rm{H}}} \right\}} \right)^{ - 1}}$。${{\hat{ u}}_{{\rm{SR}}}} = {\rm{vec}}\left( {{{{\hat{ U}}}_{{\rm{SR}}}}} \right)$为矢量化后的增广主特征矢量， ${{{W}}_{{\rm{SR}}}}$为$2N$个矩阵组成的块对角阵。由第3节的推导可知式(39)可进一步写为

${{{W}}_{{\rm{SR}}}}$中的第$N + i$个块矩阵为

与第$i$个块矩阵完全相同。因此可以得到${{{W}}_{{\rm{SR}}}}$的近似解

式中

然后分别通过${{{U}}_{\rm{S}}}$和${{{V}}_{\rm{S}}}$构造增广主特征矢量${{{U}}_{{\rm{SR}}}}$和${{{V}}_{{\rm{SR}}}}$，最后通过迭代以及角度配对估计出仰角和方位角。

5.   多目标处理能力分析

为了方便描述，将文献[11]所提算法称为CCM算法，文献[12]所提的两种算法称为CVM和MCVM算法，文献[13]所提算法称为OMP算法。下面对L-PUMA算法和L-MPUMA算法的测向能力进行比较，在L-PUMA算法中，${{b}}$由式${\hat{ F}}{{b}} \approx {\hat{ g}}$确定，其中${\hat{ F}}$是$\left( {M - K} \right)N \times K$矩阵，相当于通过$\left( {M - K} \right)N$个方程确定$K$个未知数，当$\left( {M - K} \right)N < K$时，${{b}}$存在较大的求解误差，因此需要满足$\left( {M - K} \right)N \ge K$，即最大可分辨信号数${k_{\max }} = \dfrac{N}{{N + 1}}M$。同理，对于L-MPUMA算法，${{b}}$由式${{\hat{ F}}_{{\rm{SR}}}}{{b}} \approx {{\hat{ g}}_{{\rm{SR}}}}$确定，其中${{\hat{ F}}_{{\rm{SR}}}}$是$2\left( {M - K} \right)N \times K$维矩阵，因此最大可分辨信号数${k_{\max }} = \dfrac{{2N}}{{2N + 1}}M$，相比于L-PUMA, L-MPUMA算法提高了最大可分辨信号数，并且本文提出的两种算法较以往算法都提高了信号分辨数。

6.   仿真实验

将本文提出的L-PUMA, L-MPUMA算法与CCM, CVM, MCVM以及OMP算法进行比较。假设接收阵元噪声为高斯白噪声，噪声间相互独立，每次仿真均进行1000次蒙特卡罗实验。

实验1：配对算法性能比较

为了验证本文所提配对算法的有效性，将L-PUMA算法中的配对环节用文献[17]中的PSCM配对算法替代，简称为L-PUMA2算法。同样地，将L-MPUMA算法中的配对环节用PSCM替代，称为L-MPUMA2算法。假设有3个信号入射到17个阵元构成的L型阵列上，入射信号的仰角、方位角分别为$\left( {65^{\circ} ,120^{\circ} } \right)$, $\left( {80^{\circ} ,105^{\circ} } \right)$和$\left( {95^{\circ} ,90^{\circ} } \right)$。前两个信号相干，相干系数为$\left[ {1,0.4{{\rm{e}}^{{\rm{j\pi }}/3}}} \right]$，第3个信号独立，快拍数为200。由图2可知，当$- 5\;{\rm{dB}} < {\rm{SNR}} < 5\;{\rm{dB}}$时，L-PUMA算法的RMSE明显低于L-PUMA2, L-MPUMA算法的RMSE明显低于L-MPUMA2，可见改进的配对算法提升了算法性能。

图 2  配对算法RMSE随SNR的变化

下载: 全尺寸图片 幻灯片

再分析临近目标对配对算法性能的影响，假设${\rm{SNR}} = 10\;{\rm{dB}}$，快拍数为200，阵元数保持不变。有两个相干信号入射到接收阵列上，第1个信号的仰角、方位角为$\left( {50^{\circ} ,80^{\circ} } \right)$，第2个信号的仰角、方位角为$\left( {50^{\circ} + \Delta \alpha ,80^{\circ} - \Delta \alpha } \right)$, $\Delta \alpha$表示角度间隔，其变化范围为$1^{\circ} \sim 10^{\circ}$，步长为$0.5^{\circ}$。如图3所示，随着角度间隔的增加，入射信号的空间分辨率越来越大，配对成功率随之升高，其中经过配对算法改进的L-PUMA的成功率要高于L-PUMA2，同理，L-MPUMA的成功率也高于L-MPUMA2。

图 3  配对成功率随角度间隔的变化

下载: 全尺寸图片 幻灯片

实验2：信噪比变化对算法性能的影响

假设有3个信号入射到17个阵元构成的L型阵列上，入射信号的仰角、方位角分别为$\left( {65^{\circ} ,120^{\circ} } \right)$, $\left( {80^{\circ} ,105^{\circ} } \right)$和$\left( {95^{\circ} ,90^{\circ} } \right)$。前两个信号相干，相干系数为$\left[ {1,0.4{{\rm{e}}^{{\rm{j\pi /3}}}}} \right]$，第3个信号独立，快拍数为200。由图4可知，当${{\rm{SNR}} < - 5\;{\rm{dB}}}$时，CCM算法的RMSE最低，当${{\rm{SNR}} > - 5\;{\rm{dB}}}$时，L-MPUMA算法的RMSE最低，L-PUMA算法次之，而CVM和MCVM的RMSE均较大，OMP算法的性能在所有算法中处于中间水平。

图 4  RMSE随SNR的变化(部分相干)

下载: 全尺寸图片 幻灯片

仅改变入射信号之间的相干性，假设信号完全相干，相干系数为$\left[ {1,0.3{{\rm{e}}^{{\rm{j\pi }}/10}},0.8{{\rm{e}}^{ - {\rm{j\pi }}/6}}} \right]$。由图5可知，OMP算法在低信噪比条件下性能下降严重，L-PUMA算法的RMSE整体升高，性能完全被CCM算法超越，这是由于信号的完全相干导致了线性方程个数的减少，因此估计精度下降。而L-MPUMA算法仍能保持很好的性能，当${\rm{SNR}} > - 5\;{\rm{dB}}$时，其RMSE低于任何一种算法。

图 5  RMSE随SNR的变化(完全相干)

下载: 全尺寸图片 幻灯片

实验3：快拍数变化对算法性能的影响

设${{\rm{SNR}} = 10\;{\rm{dB}}}$，快拍数从100增加至1000，其余参数均与实验2相同。如图6所示，当信号部分相干时，可以看出算法受快拍数的影响不大，RMSE变化较为平缓，各算法的性能由差到好排序依次是：CVM, MCVM和CCM, L-PUMA, OMP, L-MPUMA。

图 6  RMSE随快拍数的变化(部分相干)

下载: 全尺寸图片 幻灯片

仅改变信号相干性，将部分相干的信号改成完全相干。如图7所示，当信号完全相干时，所有算法的RMSE均被抬高，其中L-PUMA的RMSE与CVM和MCVM相当，性能下滑严重；CCM和OMP算法的性能处于中间水平；而L-MPUMA的RMSE则最低，性能最好。可见信号的相干越大，解相干算法的性能越差，L-PUMA对于信号的相干程度十分敏感，而L-MPUMA受到的影响则较小，鲁棒性较好。

图 7  RMSE随快拍数的变化(完全相干)

下载: 全尺寸图片 幻灯片

实验4：多目标处理能力

假设${\rm{SNR}} = 10\;{\rm{dB}}$，快拍数为200，共有4个相干信号入射到11个阵元组成的L型阵列上，仰角、方位角分别为$\left( {65^{\circ} ,115^{\circ} } \right)$, $\left( {80^{\circ} ,100^{\circ} } \right)$, $\left( {95^{\circ} ,85^{\circ} } \right)$和$\left( {110^{\circ} ,70^{\circ} } \right)$。假设4个信号完全相干，都属于同一组相干信号，即$N = 1$，如图8所示，星号代表目标的真实位置，可以看出L-MPUMA算法估计出的角度与真实角度比较接近，而L-PUMA算法则出现估计错误，估计值严重偏离真实值。这是因为在$N = 1$, $M = 6$的条件下，L-PUMA最多可处理的信号数为${k_{\max }} = \dfrac{N}{{N + 1}}M = 3$, L-MPUMA可处理的信号数为${k_{\max }} = \dfrac{{2N}}{{2N + 1}}M = 4$。

图 8  目标的空间分布(1组相干信号)

下载: 全尺寸图片 幻灯片

假设4个信号部分相干，分为两组相干信号，即$N = 2$，且每组信号数为2个，仰角和方位角保持不变。从理论上分析，L-PUMA最多可处理的信号数为4个，L-MPUMA可处理的信号数为4.8个，向下取整为4个。如图9所示，两种算法估计出的目标均接近真实目标，其中L-MPUMA的偏差更小，L-PUMA的偏差稍大但也能做到准确估计。与图8中信号完全相干的情况相比，此时算法的分辨率有很大提升。

图 9  目标的空间分布(2组相干信号)

下载: 全尺寸图片 幻灯片

在其他条件不变的情况下，再增加一个独立信号，新增的仰角、方位角为$\left( {125^{\circ} ,55^{\circ} } \right)$，此时共有3组相干信号，即$N = 3$，信号总数为5个，阵元数仍保持11个。由分析可知L-PUMA最多可处理4个信号，L-MPUMA最多可处理5个信号。如图10所示，L-MPUMA算法仍可以得到和真实目标较为接近的结果，而L-PUMA算法则完全失效，无法显示出目标。

图 10  目标的空间分布(3组相干信号)

下载: 全尺寸图片 幻灯片

7.   结论

本文针对L型阵列相干信号DOA估计问题，提出了两种解相干算法：L-PUMA和L-MPUMA算法，较大程度提高了估计精度和最大可分辨信号数。虽然L-PUMA算法在信号完全相干条件下性能下降比较严重，但当信号部分相干时性能良好。L-MPUMA算法则对信号的相干程度不敏感，进一步提高了精度和鲁棒性。