程思备，骆骁，蒋博筹，王玉婷，吴寰宇

（重庆信息通信研究院，重庆 401336）

1 引言

在阵列信号处理中，更高的自由度（degree of freedoms，DOF）往往能够带来更高的输出信噪比、更低的旁瓣、更低的克拉美罗下界（Cramer-Rao lower bound，CRLB）等。与均匀线阵（uniformlinear array，ULA）相比，非均匀线阵（nonuniform linear array，NLA）因为其稀疏特性能够在使用相同阵元数的情况下，获得更高的 DOF[1]。常见的非均匀线阵有最小不冗余阵（minimum redundancy array，MRA）[2]、互质阵列（co-prime array）[3-4]以及嵌套阵列（nested array）[5]。

参考文献[2]结合 MRA结构提出了基于扩展矩阵的算法，该算法的确能达到增加DOF的目的，从而获得更优异的处理性能。但是在给定阵元数的情况下，MRA的布阵方式没有相应的解析表达式，只能通过迭代算法得到其阵元位置，因此该阵列通用性较差。2010年，嵌套阵列的概念和数学模型首次被提出，该阵列由两个或者多个阵元间隔不同的均匀线阵连接而成，当阵元数和阶数确定后其阵元位置表达式可由固定计算式直接给出[5]，近几年已经成为非均匀阵列领域研究的热点。参考文献[6]提出了超级嵌套阵列（super nested array）的布阵方式，该方式保留了嵌套阵列增加自由度优点的同时也能够减少阵元间的互耦。对于嵌套阵列的DOA估计，现有参考文献采用多重信号分类 （multiple signal classification，MUSIC）[7]、传播因子与 MUSIC 结合[8]、酉变换矩阵重构与MUSIC结合[9]等算法，上述算法需要在角度域进行峰值搜索，计算量较大。传统的旋转平移不变技术（estimating signal parameters via rotational invariance technique，ESPRIT）不需要进行峰值搜索，它利用均匀阵列的导向量矩阵为范德蒙德矩阵（Vandermonde matrix），其具有全局平移不变的特性，将阵列分为去顶和去底的子阵列，而通过求解子阵列协方差矩阵特征向量之间的旋转因子来完成 DOA估计[10]，但是全局平移不变特性不具有普遍性，只有均匀阵列具有该特性，大多数阵列只有局部平移不变特性，因此传统的 ESPRIT算法不适用于嵌套阵列的DOA估计。

参考文献[11-13]利用了非均匀线阵局部平移不变特性（将非均匀线阵分为若干个的子阵列，每个子阵列为均匀线阵，其导向量矩阵均为范德摩尔矩阵，具有平移不变性），将子阵列导向量矩阵构造转化为三阶张量的构造，最后利用典范分解求解三阶张量的水平切片，得到非均匀线阵的导向量矩阵，最后利用最小二范数法则完成DOA估计。ESPRIT算法与基于典范分解的DOA估计算法类似，都利用范德摩尔矩阵平移不变特性，不同的是ESPRIT算法利用全局特性而DOA估计算法利用的是局部特性，因此基于典范分解的DOA估计方法被认为是ESPRIT算法的延伸，被定义为“Multiple Invariance ESPRIT”[13]。但是在参考文献[13]中，快拍（Snapshot）信号模型不含噪声成分，作者指出该算法只适用于无噪环境的假设前提，而在实际环境中噪声普遍存在，通过仿真发现在有背景噪声的条件下，算法的鲁棒性急剧下降，甚至出现无解。如果将此算法应用于实际工程，必须对其进行改进优化。

本文结合含有噪声成分的二阶嵌套阵列信号模型，提出了改进的典范分解的嵌套阵列DOA估计算法，并对该算法进行了仿真，得到了相同信噪比和快拍数情况下，该算法和直接MUSIC、空间平滑算法，以及克拉美罗下界[14]的均方根误差对比图。仿真表明改进的典范分解算法适用于含噪条件下的二阶嵌套阵列DOA估计，且与参考文献[5，8]中提到的算法相比，具有更好的DOA估计性能。

2 二阶嵌套阵列信号模型

N阵元的二阶嵌套阵列结构如图1所示，由两个均匀线阵连接而成：第一阶被称为内部均匀线阵（inner ULA），含有N1个阵元且阵元间隔为d=λ/ 2 ，λ为接收信号波长；第二阶被称为外部均匀线阵（outer ULA）含有N2个阵元且阵元间隔为(N1+1)d。当N为偶数时，N1=N2=N/2；当N为奇数时N1=(N- 1 )/2，N2=(N+ 1 )/2。假设有R个窄带互不相关的信号从方向为θi, (i=1 ,2,… ,R)的远场空间入射该二阶嵌套阵列，其K次接收快拍矩阵X∈CN×K表示为：

图1 二阶嵌套阵列结构

其中，矩阵A∈CN×R为入射信号的导向量矩阵A=[a(θ1),a(θ2) , … ,a(θR)]，向量a(θi)表示入射角度为θi的信号的导向量，向量si表示其K次采样，矩阵S=[s1,s2,… ,sR]T∈CR×K，上标符号T表示矩阵转秩，矩阵Wn∈CN×K表示所有阵元K次快拍接收的噪声分量，噪声为高斯白噪声满足正态分布，均值为0，方差为。

3 基于改进典范分解的嵌套阵列 DOA估计算法

3.1 典范分解算法基本原理

根据参考文献[13]，假设K次接收快拍不含噪声成分，对其做奇异值分解（singular value decomposition，SVD）得到上标符号H表示矩阵共轭转秩，是有奇异值从大到小排列而成的对角矩阵，分别为与各个奇异值对应的左特征向量矩阵和右特征向量矩阵。一定存在 矩 阵使 得对 矩 阵做双线性谱图函数映射得到矩阵，双线性谱图函数映射目的是通过线性映射改变矩阵元素使其满足典范分解的约束条件[12]，求解矩阵的核得到矩阵即求解矩阵广义特征值，即上标符号-1表示矩阵求逆，故导向量矩阵A=，利用求解得到的导向量矩阵估计入射信号角度。

典范分解算法直接对接收快拍做SVD和双线性谱图函数映射，实质上是对接收数据做一阶矩处理，在接收快拍不含噪声成分时，矩阵Q不含有噪声成分，为非满秩矩阵，它的核存在非零解，但是当接收快拍中含有噪声成分时，即接收快拍为X=AS+Wn而不再是=AS，一阶距处理无法将信号和噪声分离，矩阵含有噪声成分，利用直接求核的方法得到的矩阵只含0元素，无法估计导向量矩阵和入射角度。

3.2 改进典范分解与嵌套阵列

得到二阶嵌套阵列的K次接收快拍X后对其做SVD得到：

其中，∑∈CK×K是有奇异值从大到小排列而成的对角矩阵，矩阵U∈CN×R，矩阵V∈CK×K分别为与各个奇异值对应的左特征向量矩阵和右特征向量矩阵。UR=U(:,1:R)，Σs∈CR×R为信号空间对应的奇异值对角矩阵，且ΣR=Σ( 1:R,1:R)。参考文献[11，14]，一定存在矩阵F∈CR×R，使得URΣR=AFT。定义矩阵Y∈CN×R且Y=URΣR。

Λ(θ)=表示R×R的对角矩阵，被称为阶间旋转因子矩阵，其维度由入射信号个数决定，对角元素由入射信号角度和阶间阵元距离差共同决定，根据二阶嵌套阵列结构特性可得Y(2)=Y(1)×Λ(θ)。定义矩阵表示去掉矩阵第一行得到的子矩阵，=Y(1)(2:N1,:) , 矩阵表示去掉矩阵最后一行得到的子矩阵，同理可得矩阵表示向量J的第n个元素J=[N1- 1N2- 1 ]，定义矩阵

其中，i,j∈ [1,… ,Jn]，k,l∈ [1 ,2]。

如第 3.1节所述，矩阵Q含有噪声成分，如果按照典范分解的方法采用直接求核的方法得到的矩阵只含0元素，无法估计导向量矩阵和入射角度。参考文献[16]利用信号和噪声在二阶统计量空间正交的特性，采用奇异值分解将矩阵Q的协防差矩阵分解为正交的信号子空间和噪声子空间，再利用噪声子空间得到矩阵Mr，过程如下：对QHQ做奇异值分解得到中R个最小的奇异值构成递增排列集合σ=[σ1σ2...σR],即σ1≤σ2≤ … ≤σR，该集合中每个元素作为奇异值时所对应的右特征向量被定义为,将m(r)的元素下角标定义为

重新排列各元素得到矩阵：

得到矩阵Mr后，将求解矩阵F∈CR×R的问题转化问联合对角化的问题，即求解矩阵G-1=F：

当R=2时，即入射信号个数为两个，联合对角化问题退化为广义特征值分解（generalized eigenvalue decomposition，GEVD）问题，当R≥3时，需用迭代算法求解数值解，可参考文献[16-18]。求得矩阵F后，利用URΣR=AFT，求解得到导向量矩阵A：

定义导向量矩阵A的第i列列向量为ci A，a(θ)(1)=i可得，定义以下函数则入射角iθ的计算式为：

算法步骤归纳如下：

步骤 1输入接收快拍X=AS+Wn,做奇异值分解得到Y=URΣR；

步骤 2对Y进行分组重排得到Yred(n,r),利用函数φr,s(n),得到矩阵QHQ；

步骤3对QHQ进行奇异值分解，对奇异值和右奇异值向量进行重排得到矩阵Mr；

步骤4对Mr进行联合对角化，求得矩阵F；

步骤5利用URΣR=AFT，求得导向量矩阵A；

步骤6利用式（13）求得入射角度θi。

改进的典范分解算法与传统的典范算法相比，主要改进体现在步骤1和步骤3上。步骤1中，传统算法的接收快拍不含噪声，而实际应用中噪声是普遍存在的，所以改进算法的接收快拍中加入噪声。步骤2对接收快拍进行了分组和线性映射，传统算法得到的不含噪声，故在步骤 3中对其直接求核可以得到含有可行解的，而改进算法得到的Q中含有噪声，故在步骤（3）中对其协方差矩阵进行SVD，利用正交特性得到含有角度信息的Mr，剩下步骤两种算法相同。

假设接收快拍不含噪声，在步骤3中对的协方差矩阵做SVD，得到的噪声子空间奇异值均为0，然后求解其对应的右特征向量矩阵，问题等效为=0与直接求核结果一样。因此在无噪条件下改进算法与传统算法性能一致。

综上所述，改进算法与传统算法计算复杂度基本相当，但是传统算法只能应用于无噪条件，而改进算法在无噪和有噪条件下均适用。

4 仿真实验

仿真实验采用8阵元二阶嵌套阵列，即内部均匀线阵含有N1=4个阵元，各阵元间隔为d，外部均匀线阵含有N2=4，各阵元间隔为5d，阵元分布集合为[d, 2d, 3d, 4d, 5d, 10d, 15d, 20d]。仿真在不同信噪比和快拍次数的情况下，分别进行200次蒙特卡洛实验，对比MUSIC、空间平滑和本文所提的改进典范分解算法 DOA估计结果的均方根误差，并以克拉美罗下界值作为参考，定量分析算法估计性能。

实验中信噪比（SNR）、均方根误差（ERMS）和克拉美罗下界值（ C RBL）计算式分别为:

其中，表示二阶范数，R表示信号个数，MC表示蒙特卡洛实验次数，本文中MC=2 00，表示第mc次实验所得估计角集合。

4.1 仿真实验1

3个功率相同的信号分别从30°、45°、60°入射8阵元二阶的嵌套阵列，信源数R=3，入射集合快拍次数K=50，信噪比从0 dB变化到20 dB，3种算法的均方根误差如图2所示。

图2中，8阵元二阶嵌套阵列比8阵元均匀线阵列拥有更高的自由度，因此在相同信噪比和快拍次数的情况下，8阵元二阶嵌套阵列的克拉美罗界明显低于8阵元均匀线阵列，体现了相同条件下，二阶嵌套阵列的性能优越性。对于直接MUSIC、空间平滑和改进典范分解算法的均方根误差，信噪比小于2 dB 时，3种算法均方根误差接近；而信噪比大于2 dB时，本文所提的改进典范分解算法的均方根误差明显低于其余两种算法；当信噪比大于10 dB时，改进典范分解算法对8阵元二阶嵌套阵列的均方根误差比直接 MUSIC算法小一个数量级，且低于8阵元均匀线阵的克拉美罗界。图2验证了改进典范分解算法对嵌套阵列的可行性和体现出性能优越性。

图2 均方根误差随信噪比变化曲线

4.2 仿真实验2

采用仿真1中的入射角度分别为30°、45°、60°的3个信号信噪比为10 dB，快拍次数从5变化到500，均方根误差如图3所示。

图3 均方根误差随快拍数变化曲线

图3中，当快拍次数小于50时，3种算法的均方根误差受快拍次数变化影响较大，当快拍次数大于50时，变化趋势趋于平缓，而实际情况中快拍数常常高于50，在这一区间内3种算法的均方根误差受快拍数影响不大。仿真结果表明当信噪比为10 dB时，在相同快拍次数的情况下，改进典范分解算法均方根误差小于直接 MUSIC算法和空间平滑算法，表明所提算法DOA估计精度高于另外两种算法。

4.3 仿真实验3

该仿真参考了文献[19]所采用的在相同信噪比和快拍数的前提下，通过比较运行时间来对比运算复杂度的方式，使用同一计算机和同一软件（实验环境：CPU 2.5 GHz，8 GB RAM，MATLAB R2014a，Windows 8 x64）对本文所提算法、直接MUSIC算法和空间平滑算法，在快拍数为200和信噪比为10 dB（参考图2和图3，角度域搜索步长取0.05°）的情况下，分别做100次、200次和500次蒙特卡洛分析得到3种算法的运行时间，结果见表1。从表1可以看出，本文所提算法不需要峰值搜索，复杂度明显低于直接MUSIC算法和空间平滑算法。

表1 算法运行时间比较

5 结束语

本文结合二阶嵌套阵列含噪信号模型，首次提出了基于改进典范分解的嵌套阵列 DOA估计算法。该算法无需进行峰值搜索，利用局部平移不变特性，通过张量构造和典范分解求解得到导向量矩阵并求得入射角度。该算法利用SVD解将矩阵Q的协防差矩阵分解为正交的信号子空间和噪声子空间,再利用噪声子空间得到矩阵r M，解决了含噪条件下典范分解算法对嵌套阵列 DOA估计的问题，且在相同信噪比、快拍数情况下比直接MUSIC算法、空间平滑算法有更好的估计性能和更少的运算时间。