朱立东，邓超升

(电子科技大学通信抗干扰技术国家级重点实验室，四川成都611731)

0 引言

在无线移动通信系统中，同频干扰(Co-Channel Interference，CCI)、多址干扰(Multiple Address Interference，MAI)是影响系统容量和性能的因素之一。智能天线技术为解决频率资源不足、提高移动通信系统容量和系统服务质量提供了一个有效的解决途径。在智能天线技术的研究中，入射信号的DOA(Direction Of Arrival)估计作为一个重要的研究内容受到广泛关注。根据入射信号带宽的不同，DOA估计方法可分为两大类:窄带信号DOA估计和宽带信号DOA估计。常用的宽带信号DOA估计方法大致分为信号子空间方法和最大似然估计方法。在信号子空间方法的几个分类中，非相干信号子空间方法(Incoherent Signal-Subspace Method，ISSM)[1]是最简单的宽带信号处理方法，该方法在高信噪比时简单有效，但在低信噪比下对某些频段信号的DOA估计效果非常不理想，降低了整体的估计性能，且该方法不能处理相干源。基于角度预估计的相干信号子空间方法(Coherent Signal-Subspace Method，CSSM)[2]通过构造聚焦矩阵，将预估角度处每个频点的阵列流形对齐到聚焦频率处，通过阵列流形替换实现DOA估计，该方法在克服ISSM方法缺陷的同时容易引入角度预估计误差，从而影响DOA估计性能。为进一步消除角度预估计误差带来的影响，近年来出现了很多无需角度预估计的宽带波达方向估计算法[3]。文献[4-6]提出基于波场模型的阵列流形内插方法(Array Manifold Interpolation，AMI)，将阵列流形分为仅与阵列结构和频率有关的采样矩阵和仅与波达方向有关的向量，并且对二维及三维阵列的波场模型分别进行建模;文献[7]针对稀疏均匀圆阵，给出基于流形分离技术(MST)的DOA估计方法;文献[8]将基于流形分离技术的波达方向估计方法应用到非均匀圆阵上;文献[9]在直线阵的基础上，给出一种采用傅里叶-勒让德级数展开构造变换矩阵的宽带聚焦估计方法;文献[10]研究了不同阵列结构对空间角度估计的影响，并用仿真分析给出几种常见阵列在低俯仰角下空间到达角的估计性能，验证均匀圆阵在信号空间到达角估计中的优越性。

在以上文献和方法的启发下，提出一种均匀圆阵下宽带相干信号的二维来波方向估计新方法。先使用多项式展开的方法实现阵列流形矢量的频率与角度的分离，再选取简易的与频率有关的采样矩阵构造聚焦矩阵，在保证聚焦前后无损失的基础上减小计算复杂度。此外，针对均匀圆阵在俯仰角较小或较大时对空间角估计精度较差的问题提出有效的解决方法。仿真分析表明，所提方法可以获得较优的2D-DOA估计效果，且该方法可推广至任意的阵列结构，并明显改善2D-DOA估计性能。

1 信号模型和算法描述

1.1 信号模型

对于任意一个含有M个阵元的几何平面阵列，设第m个阵元的位置为(xm，ym，0)，接收N(M≥N)个位于远场的宽带平稳随机信号，噪声为加性带通高斯白噪声，且噪声与源信号相互独立。则第m个阵元的接收信号在时域中的表达式为

其中，sn(t)为第n个信号源;nm(t)为第m个阵元上的加性高斯白噪声;τmn= (xmcosφn+ymsinφn)sinθn/c为第m个阵元接收源信号相对于参考点阵元接收信号的传播时间延迟，其中，c为光速，(θn，φn)为第n个源信号的俯仰角和方位角，m=1，2，…，M ，n=1，2，…，N 。

对阵列的接收信号进行离散傅里叶变换，并将整个信号带宽划分为J个子频带，各子频带的中心频率分别为 f1，f2…，fj，…，fJ，则接收阵列在频域内对应于频率fj的输出信号可以表示为

其中，j=1，2，…，J;X(fj)=[X1(fj)，X2(fj)，…，Xm(fj)，…，XM(fj)]T，Xm(fj)为第 m 个阵元对应于频率 fj的输出;A(fj)=[a1(fj)，a2(fj)，…，am(fj)，…，aM(fj)]T为对应于频率fj的阵列流形矩阵，am(fj)为导向矢量，其表达式为am(fj)= [ei2πfjτm1，ei2πfjτm2，…，ei2πfjτmn，…，ei2πfjτmN]，其中;S(fj)和 N(fj)分别为信号和噪声在第fj频段上的输出矩阵，S(fj)=[S1(fj)，S2(fj)，…，SM(fj)]T，N(fj)=[N1(fj)，N2(fj)，…，NM(fj)]T。

1.2 常规算法

文献[11]对傅里叶变换在平面、柱、球阵波束形成和DOA估计中的应用进行全面分析。对于平面阵列，其阵元坐标可通过相应的映射关系将其转换为球坐标，转换关系为

在球坐标下，设阵元坐标向量为 p=Rm[sinθmcosφm，sinθmsinφm，cosθm]T，同时设频率 fj对应的第n个信号Sn(fj)的波束向量为k=[kjsinθncosφn，kjsinθnsinφn，kjcosθn]T，其 中 kj=－2π/λj=－2πfj，则信号Sn(fj)在第m个阵元处的接收信号可以表示为

因此，设 N个信号的来波方向为 Ζ(N)={(θ1，φ1)，…，(θn，φn)，…，(θN，φN)}，第m个阵元对应于频率fj的导向矢量αm(fj)可以写为αm(fj)=[αm(fj，θ1，φ1)，…，αm(fj，θn，φn)，…，αm(fj，θN，φN)]。其中，αm(fj，θn，φn)= exp{ikjRm[sinθnsinθmcos(φn－φm)+cosθncosθm]}，αm(fj)中的元素 αm(fj，θn，φn)可用球形贝塞尔函数及球函数的级数展开表示。

其中jξ(·)为 ξ阶球形贝塞尔函数，Ylξ(·)为 ξ阶l度的球函数。Ylξ(θ，φ)可以表示为

其中，Plξ(·)为ξ阶l度的连带勒让德函数，表达式为

归一化的球函数可以表示为

针对球贝塞尔函数，其与贝塞尔函数(第一类)之间的转化关系如下

将式(8)和式(9)带入式(5)，则 αm(fj，θn，φn)可以表示为

其中，C=(2ξ+1)(ξ－l)!/((ξ+l)!)。

从式(11)可知，αm(fj，θn，φn)可以分为两项，一项仅与阵列的结构和频率有关，另外一项与信号的来波方向有关。因此，阵列流形矩阵A(fj)可以分为两个部分，一部分与阵列结构和频率有关，另一部分只与波达方向有关。设仅与阵列结构和频率有关的采样矩阵 为 G(fj)= [G1，…，Gm，…，GM]T，Gm=[g00，g1－1，g01，g11，…，]T，故导向矢量 αm(fj)又可以表示为αm(fj)=Gm·Ψ，Ψ为关于ψln的向量。Gm的表达式中对于矩阵G(fj)，根据贝塞尔函数的性质，可以选取适当的截断个数Nc(通常取Nc=2×(2πfjRm/c))，得到截断后的G(fj)为M×(Nc+1)2维的矩阵，之后可选取相应的方法构造聚焦矩阵。几种常见的方法有:基于波场模型的阵列流形内插方法[5]和旋转信号子空间的阵列流形内插方法[12]。使用聚焦矩阵可获得性质类似于窄带相关矩阵的宽带相关矩阵，进而获得宽带信号的DOA信息。

2 改进算法

考虑一个以圆心为参考点的均匀圆阵，设阵列的半径为宽带信号起始频率对应波长的p倍，则Rm为定值，即 Rm=pλ0=pc/f0。又 kj=－2π/λj=－2πfj/c，则 kjRm=2πfjp/f0。此时，式(10)中球形贝塞尔级数项的系数可化简得到

将(12)式带入(10)式，αm(fj，θn，φn)可以化简得到

通过带入归一化球函数的表达式，可以将αm(fj，θn，φn)进一步写为

其中，C1为常数项，表达式为C1=(2ξ+1)(ξl)!/(4π(ξ+l)!)。通过进一步化简与常数项合并，αm(θn，φn，fj)可以表示为

其中，C1为常数项，表达式为1)(ξ－l)!/((ξ+l)!)。

在此基础上，考虑信号带宽。对于宽带信号，其相对带宽定义为:ffoc=2(fH－fL)/(fH+fL)，其中fH和fL分别为宽带信号的最高频率和最低频率，fL=f0。若将整个频率带划分成J个子频率带，则第j(j=0，1，…，J－1)个子频率可以表示为

令ρ=fH/fL，显然ρ＞1，由此可以得到

在均匀圆阵圆心的轴上增加一个全向阵元，设该阵元的编号为 M+1 ，为了保证 αM+1(fj，θn，φn)的可分离性，RM+1不能为0。设RM+1=qλ0，选取q为一个较小的数，可使包含阵元M+1的立体阵趋近于平面阵，此阵列布局易于工程实现。对于该阵元，其在球坐标系下的坐标可以表示为:RM+1=qλ0，θM+1=0°，φM+1可以为任意值。当阵元处于Z轴上时，φM+1取值已无意义，但考虑到构造(fj)的需要，因此对其进行赋值。在φM+1=0°的情况下，对于阵元M+1，αM+1(fj，θn，φn)=exp(ikjRM+1cosθn)，显然，该项可以用球贝塞尔函数及球函数的级数形式表示，并划分成两个部分:一部分只和该阵元的位置和频率有关，另一部分只和来波方向和对应该频率的衰减系数有关。

(19)式的拟合解为

其中，Uj、Vj分别为 G(f0)G(fj)H的左、右奇异矢量组成的矩阵。求得聚焦矩阵Tj后，进而可求得对应于每个频率的聚焦变换后的相关矩阵 ^R(fj)，其表达式为

对聚焦后的数据相关矩阵进行相加求平均，可得到相关矩阵R，其表达式为

由式(15)可知，R可视为一个窄带相关矩阵，且包含原宽带信号的所有信息。因此，可按照窄带信号DOA估计方法的步骤对该相关矩阵进行处理，从而得到宽带信号的DOA估计信息。此外，由于R的构成过程中对所划分的各个频率分量进行平滑处理，因而R具有解相干性能，能处理宽带相干信号。

根据以上理论分析，采用文中方法构造聚焦矩阵的宽带相干二维DOA估计方法的步骤可总结如下:

(1)将阵列接收信号的快拍数据分为K段，每段含有L个快拍，同时将频带划分为J个频段;

(2)分别对每段数据进行DFT变换，得到J个频段的阵列接收数据;

(3)分别对J个频段的接收数据构造频域采样空间相关矩阵:

(4)对G(f0)G(fj)H进行奇异值分解，求得Uj、Vj，再利用式(20)求出聚焦变换矩阵Tj;

(5)根据式(22)得到组合加权后的相关矩阵，并进行特征值分解，得到噪声子空间和信号子空间;

(6)利用二维MUSIC算法进行来波方向估计。

所述算法在计算聚焦变换矩阵Tj的过程中，尽可能剔除与频率无关的乘积项来化简采样矩阵(fj)，并且巧妙地结合均匀圆阵的特点，进一步化简(fj)，在提高聚焦效果的同时减小了计算量。此外，根据空间角估计精度随俯仰角的变化情况，可适当调整阵列结构，能够改善估计性能。

3 数值仿真分析

假设阵列的接收信号为远场宽带信号，中心频率为f0=100 Hz，信号相对带宽为40%，采样频率fs=256 Hz，噪声为加性零均值带通高斯白噪声，噪声与源信号相互独立，阵列接收数据的快拍数为4096个，把数据分为32段，每段含有128个快拍。信号带宽为40 Hz，可把整个带宽划分为41个子频带，选择中心频率f0为聚焦频率。根据贝塞尔函数的性质，选取截断点数Nc为9。为便于观察和分辨，将常规多项式展开方法(general polynomial expansion method)简记为GPEM;文中方法(improved polynomial expansion method)简记为IPEM;改进后的均匀圆阵(improved uniform circular array)简记为IUCA。

仿真1:验证文中算法与常规算法的均方根误差。假设进行n次独立重复实验，N个信号源的空间角的均方根误差定义为式中θd和φd为信号的真实俯仰角和方位角，^θdi和^φdi为估计值。设接收阵列的阵元个数为11，阵列构型分别为UCA和IUCA，仿真使用两个完全相干的远场时域平稳宽带信号，其均值为0，信号到达角分别为(40°，40°)、(20°，200°)，进行 100 次蒙特卡洛仿真，统计文中算法与常规算法在不同信噪比条件下的空间角的均方根误差，对比曲线如图1所示。

仿真2:统计文中算法与常规算法的成功分辨率。使用两个完全相干的信号进行仿真，其来波方向分别为(40°，50°)和(60°，75°)，其他仿真条件同仿真 1。当两个信号的俯仰角和方位角的估计值均小于0.2°时视为成功估计，进行200次蒙特卡洛仿真，成功分辨性能曲线如图2所示。

仿真3:验证文中方法在其他任意阵列上的空间角估计性能。仿真中使用的任意阵列的布局如图3所示，对任意阵的改进操作为:将任意阵中编号为M的阵元置于参考点上方合适位置，总体阵元个数不改变。两个相干信号的来波方向分别为(20°，50°)和(35°，55°)，其他仿真条件同仿真1。得到的均方根误差性能曲线如图4所示。

仿真4:验证文中算法对两个来波方向相近的宽带相干信号的分辨性能。设两个相干信号来波方向分别为(10°，200°)和(15°，200°)，文中算法在11 个阵元组成的IUCA的阵列构型下对两个信号的俯仰角分辨性能如图5所示。

图1 不同方法在不同信噪比下的空间角估计性能

图2 不同方法的成功分辨概率

图1分别对比了改进的算法和常规算法分别在11阵元均匀圆阵、10阵元均匀圆阵与圆心轴上的阵元组成的11阵元阵列下的空间角估计性能。无论是使用均匀圆阵(UCA)或是改进的均匀圆阵(IUCA)，IPEM算法在空间角的均方根误差上的性能都优于GPEM。在UCA下，两种算法的均方根误差在13～20 dB内大致趋于稳定，上下波动在0.01°范围以内。其中，IPEM算法的均方根误差在20 dB时最小，大小为0.1227°;GPEM算法的均方根误差在20 dB时最小，大小为0.1468°。两种算法的均方根误差的差值在0.0220°和0.0707°之间波动。在IUCA下，使用两种算法得到的DOA估计的均方根误差逐渐减小，随着信噪比从0 dB到20 dB的变化，IPEM算法得到的均方根误差从0.1441°降低到0.054°，GPEM算法得到的均方根误差从0.2031°降低到0.079°。两种算法的均方根误差的差值在0.0250°和0.0589°之间波动。此外，GPEM算法在ICUA下的均方根误差相对于在UCA下降低的幅度大于0.04°。同时，IPEM算法在IUCA下的性能提升幅度同样大于 0.04°，并且最大达到0.0704°。

由图2可以看出，在UCA下，IPEM算法在信噪比为－1 dB时的成功分辨概率为81.5%，GPEM算法的成功分辨率为73%，IPEM算法在GPEM算法的基础上成功分辨概率提升了8.5%。在信噪比为－5～14 dB，改进后的 IPEM与常规的多项式展开方法GPEM相比成功分辨概率都有所提升，提升幅度范围为1%～8.5%。在阵列构型为IUCA的情况下，IPEM算法和GPEM算法在信噪比高于6 dB条件下的成功分辨概率都为100%。在信噪比为－5～5 dB，IPEM算法的成功分辨概率明显高于GPEM算法的成功分辨概率，提升幅度范围为0.25%～15%。其中，在信噪比为－5 dB时提升幅度为15%。总体来说，IPEM算法在相同阵列构型和相同信噪比的条件下，成功分辨概率高于GPEM算法。此外，同一算法在不同阵列构型下的成功分辨概率也有所不同。在信噪比为－5～15 dB，两种算法在ICUA阵列下的成功分辨率比在UCA阵列下有所提升，IPEM算法最高提升了12.75%，GPEM算法最高提升了11.25%。

分析出现图1和图2中统计结果的原因可知，两种算法的区别仅在于用于构造聚焦矩阵的采样矩阵G(fj)不同，在构造G(fj)的过程中，重点在于将不同频率上的阵列流形统一到聚焦频率上的阵列流形，IPEM算法尽可能剔除与频率无关的项，并巧妙结合均匀圆阵的特性，进一步将 ^G(fj)化简，提高了聚焦效果。在改进的阵列结构下，两种算法的性能在均匀圆阵的基础上提升明显。空间角的均方根误差不仅与算法有关，还与宽带信号的入射角度有关。仿真(1)中使用的俯仰角20°属于较低的俯仰角，如文中分析，该角度处sin(θ)的值较小，但sin(θ)的变化率较大，此时对方位角的估计精度不高，相当于有很小的阵列孔径，改进阵列在原阵列的基础上加上了cos(θ)项，该项的加入改善了低俯仰角下的方位角的精度，进而提升了空间角估计性能。在仿真(2)中使用的相干信号中，其中有一个信号的俯仰角为60°，属于较高俯仰角，该俯仰角的正弦值sin(θ)较大，但是sin(θ)的变化率很低，此时对俯仰角的估计精度不高，影响了算法的成功分辨性能。改进后的阵列有效地解决了均匀圆阵的不足，明显提高了阵列的成功分辨性能。

图3 某任意几何阵列示意图

图3中将原点作为参考点，各点的坐标表示各个阵元相对于参考点的位置，x轴和y轴的单位长度为宽带相干信号起始频率所对应的波长λ0。在图3所示的阵列构型下统计IPEM算法得到的两个宽带相干信号的俯仰角与方位角的均方根误差曲线如图4。由图4可以看出，在两个不同的阵列构型下，IPEM算法在俯仰角上的均方根误差变化不明显，尤其是当信噪比高于9 dB以后俯仰角的变化幅度都在0.004°范围内。但在方位角上的均方根误差变化明显，最大差值达到0.0516°，在改进阵列下，均方根误差平均减小了0.0417°。以上数据表明，方法可以推广至任意阵列，在改进阵列的基础上结合文中算法IPEM能够提高宽带相干信号的DOA估计性能。根据文中分析，对于来波方向分别为(20°，50°)和(35°，55°)的两个相干信号，其方位角的估计精度较低，如图4所示。若存在参考点上方的阵元，对具有较低俯仰角的接收信号而言，该阵列结构对俯仰角的估计精度影响不明显，但是可以有效地降低方位角的估计误差，从而提高总体的DOA估计性能。

图4 文中算法在任意阵列下的性能

图5 文中算法的分辨性能

由图5可以看出，在信噪比为0 dB、5 dB和10 dB的条件下，文中算法对仿真所设的两个相干宽带信号的俯仰角估计的平均误差分别为0.205°、0.135°和0.06°。通过仿真结果可以看出，采用文中方法进行宽带相干信号的二维DOA估计时，在信噪比为0 dB以上能够成功分辨出空间中相隔较近的相干源。以方位角相同、俯仰角相隔5°的两个宽带相干信号为例，随信噪比的增加，信号的空间谱峰的幅度不断增大。在信噪比从0 dB增长至10 dB的过程中，谱峰高度从70.52和102.8增大至1244和1754，且谱峰变得更加尖锐，谱峰所对应的俯仰角与真实俯仰角的平均误差从0.205°逐渐减小至0.06°，说明文中算法随着信噪比的不断增加分辨性能也逐渐提高。

4 结束语

均匀圆阵具有几乎无方向性、二维性和无模糊范围广等特点，应用均匀圆阵进行宽带测向能够得到良好的效果。以均匀圆阵为接收模型，首先使用多项式展开的方法对阵列流形矢量进行频率与角度的分离，通过分析阵列流形矢量的级数展开式，结合均匀圆阵的特点，构造较简易的与频率相关的采样矩阵构造聚焦矩阵，可以达到减少计算量的目的。在此基础上，对阵列构型进行一定调整可明显改善均匀圆阵或其他阵列在俯仰角较大时或较小时对空间角角估计精度较差的问题。数值仿真分析显示，文中方法在宽带相干信号的2D-DOA估计中，能够降低信号空间到达角度估计的均方根误差，提高相干信号的成功分辨概率。同时，该算法对空间中相隔较近的相干源具有良好的分辨性能，且该方法也可推广到任意的阵列结构。

[1] Wax M，Shan T J，Kailath T.Spatio-temporal spectral analysis by eigenstructure methods[J].IEEE Transactions on Acoustics Speech and Signal Processing，1984，32:817－827.

[2] M A Doron，A J WEISS.On focusing matrices for wide-band array processing[J].IEEE Trans on Signal Process-ing，1992，40(6):1295－130.

[3] Bo L，Xiong J，Luo L.A novel wideband DOA estimation method using direction-free focusing matrix[C]//Computer Science and Network Technology(ICCSNT)，2013 3rd International Conference on.IEEE，2013:1065－1069.

[4] M A Doron，E.Doron.Wavefield modeling and array processing.Ⅰ:Spatial sampling[J].IEEE Trans Signal Process，1994，42(10):2549－2559.

[5] M A Doron，E Doron Wavefield modeling and array processing.Ⅱ:Algorithms[J].IEEE Transaction Signal Process，1994，42(10):2560－2570.

[6] M A Doron，E Doron.Wavefield modeling and array processing.Ⅲ:Resolution capacity[J].IEEE Trans Signal Process，1994，42(10):2571－2580.

[7] 潘捷，周建江，汪飞.基于流形分离技术的稀疏均匀圆阵快速DOA估计方法[J].电子与信息学报，2010，32(4):963－966.

[8] 卢海杰，章新华，熊鑫.流形分离在非均匀圆阵上的应用[J].兵工学报，2011，32(9):1113－1117.

[9] 林晋美，侯卫民，涂英，等.傅立叶-勒让德级数展开的宽带聚焦[J].声学技术，2007，12，26(6):1235－1239.

[10] 顾建峰.空间角度估计的阵列结构及其算法的研究[D].成都:电子科技大学，2004.

[11] 陈丰.Fourier变换在平面、柱、球阵波束形成和DOA估计中的应用[D].杭州:浙江大学，2006.

[12] 李刚，赵春晖，宁海春.用于宽带测向的广义阵列流行内插算法[J].弹箭与制导学报，2008，28(1):277－279.