一种高精度低复杂度的改进Root-MUSIC算法

佘黎煌, 刘平凡, 张 石, 许方晗

(东北大学 计算机科学与工程学院, 辽宁 沈阳 110169)

阵列信号处理技术包括DOA估计技术及波束形成技术,广泛应用于空间信号探测、雷达及水声信号探测等诸多领域[1].其中DOA估计技术作为获得信号波达方向角度的重要技术,如何使DOA估计算法获得更高的估计精度和更低的时间复杂度一直是科研人员的主要研究方向.

在众多DOA估计算法中,Root-MUSIC算法以其较高的估计精度、较低的复杂度受到众多科研人员的重视[2],大量基于Root-MUSIC算法的改进算法被提出.传统Root-MUSIC算法存在的最大问题在于随着阵列规模的增大,算法的运算负担成倍增加.在较低阵列规模时,多项式求根使算法运行效率较高,而当阵列规模增大后,其求根多项式阶数快速升高,算法运行耗时甚至超过谱峰搜索类算法[3].目前众多基于Root-MUSIC算法的改进算法都着重针对求根多项式进行降维操作,通过较低的数据规模实现DOA估计.Ren等[4]提出了一种快速Root-MUSIC算法,通过将噪声子空间进行拆分,构造新的数据向量,使得求根多项式阶数降低至仅为信号源个数,大大降低了算法运行复杂度.由于噪声子空间信息利用率较低,算法精度损失严重,王新贺等[5]在Ren等[4]拆分噪声子空间矩阵的基础上重构数据向量,充分利用噪声子空间,算法精度有了大幅度提升,但运算复杂度也随之增大,同时算法在阵列规模较大时完全失效,不具备鲁棒性.Yan等[6]将均匀线阵拆分为两个相等子阵列,利用两个子阵列采集的数据构造3个无噪声互相关矩阵,用处理后的新数据矩阵进行DOA估计,并将此DOA估计数据作为参考值,重复以上方法，进行二次估计得到最终DOA估计值[7].算法通过拆分子阵构造新互相关矩阵的方法避免了特征分解运算,降低了求根运算的复杂度,但在低阵列规模时算法精度损失严重.探测阵列必须为偶数,可估计信源个数较少,鲁棒性较差.本文在Toeplitz矩阵重构及求根多项式降阶的基础上,提出了一种改进的Root-MUSIC算法,算法依据数据观测矩阵首行重构Toeplitz矩阵.特征值分解后得到噪声子空间,并将噪声子空间翻转拆分,重构新的求根多项式向量,进而通过求根方法得到DOA估计值.算法通过重构Toeplitz矩阵使得新自相关矩阵具备Hermitian性质,相比于传统自相关矩阵特征值分解算法,本文算法拥有更高的DOA估计精度.噪声子空间的翻转拆分重构了新的低阶求根多项式,有效降低了算法的时间复杂度.

1 阵列信号模型

设存在M个满足远场条件的各向同性阵元所组成的均匀线阵模型,同时有q个远场窄带目标信号源,其入射角为θk(k=1,2,…,q),并且满足M>q.选择首个阵元作为参考阵,阵元间距设为d,则此时阵列响应向量为

(1)

第m个阵元的输出为

(2)

式中:sk(t)为第k个由天线阵列接收到的信源信号;λ为信号的波长;q为信号源数目;nm(t)为阵元接收的高斯白噪声,方差为σ2.此时的阵列接收数据可以表示为向量形式:

X(t)=AS(t)+N(t) .

(3)

式中:

A=[a(θ1),a(θ2),…,a(θq)] ;

(4)

S(t)=[s1(t),s2(t),…,sq(t)]T;

(5)

N(t)=[n1(t),n2(t),…,nM(t)]T;

(6)

X(t)=[x1(t),x2(t),…,xM(t)]T.

(7)

通过式(3)可以得到理想观测数据的协方差矩阵:

R=E[X(t)XH(t)]=ARsAH+σ2IM.

(8)

式中:Rs=E[S(t)SH(t)]为信号的协方差矩阵;IM为M阶的单位矩阵.

2 本文方法描述

2.1 Toeplitz矩阵重构算法

由于硬件设备采样快拍数的限制,由数据观测矩阵获得的自相关矩阵并不能满足理想的Hermitian性质,仅能保证对角占优,因此存在估计误差.Toeplitz矩阵重构算法通过将数据观测矩阵的第一行即首个阵元获得的数据作为参考向量,构造新的自相关矩阵使得新矩阵具备Hermitian性质.相较于原始自相关矩阵更加逼近理想值,DOA估计精度更高,同时不会产生阵列孔径损失,估计信源个数不会受到影响[8].

由式(2)可以得到数据观测矩阵第一行即首个阵元的输出为

(9)

在实际情况下,由于采样快拍数有限,因此天线阵列所接收到的信源信号及高斯白噪声均表现为由若干采样值所组成的向量形式,向量维度即采样快拍数,这里将采样快拍数设为T,则有

x1=[x1(1),x1(2),…,x1(T)] .

(10)

同理,第k个阵元的接收数据向量表示为

xk=[xk(1),xk(2),…,xk(T)] .

(11)

将式(10)和式(11)作为参考向量,定义

(12)

式中,A(1)和A(k)分别代表矩阵A的第一行和第k行.随着k的变化,可以得到包含全部信源信息的相关矢量[rtp(0),rtp(1),…,rtp(M-1)],由M个相关函数构成的矩阵:

Rtp=

(13)

2.2 Root-MUSIC算法

谱峰搜索算法通过全范围角度的遍历来搜索DOA估计值,在阵元数较少时会导致算法运行耗时增大.而求根MUSIC算法(Root-MUSIC)将这一过程省去,替代为简化的多项式求根过程.

通过式(1)可以得到阵列响应向量,这里定义:

(14)

式(1)的阵列响应向量可以改写为

(15)

式(4)阵列流型矩阵可以表示为

A=[z1,z2,…,zq] .

(16)

根据式(8)可以得到自相关矩阵R,对R做特征分解,较大的q个特征值对应特征向量构成信号子空间Us,较小的M-q个特征值构成噪声子空间Un.依据MUSIC算法相关理论,信号子空间与噪声子空间相互正交,则噪声子空间也必与信号子空间的基即阵列流型矩阵A正交,满足

AHUn=Oq×(M-q).

(17)

式中，O为q×(M-q)阶零矩阵.根据式(17)定义噪声子空间Un:

Un=[un1,un2,…,un(M-q)] .

(18)

此时式(17)可以表示为

(19)

通过观察式(19),可将阵列响应方向向量定义为一般形式:

z=[1,z-1,…,z-(M-1)]T.

(20)

根据式(19)及式(20)构造多项式向量p(z):

p(z)=[p1(z),p2(z),…,pM-q(z)]T.

(21)

式中:

pi(z)=zHuni=β0,i+β1,iz+…+βM-1,izM-1,

1≤i≤M-q.

(22)

式中，uni=[β0,i,β1,i,…,βM-1,i]T为噪声子空间列向量.

由式(21)可知,M-q个多项式pi(z)的阶数均为M-1,且有q个公共根z1,z2,…,zq.

根据文献[7]构造多项式:

froot(z)=pH(z)p(z) .

(23)

求得式(23)的根即可获得有关信号源到达角的信息.因多项式存在z*项,这会使得求根过程变得复杂,因此对式(23)修正:

froot(z)=zM-1pT(z-1)p(z) .

(24)

在式(24)中,多项式阶数为2M-2,即有M-1对根,且每对根为相互共轭关系.在这些根中有q个根z1,z2,…,zq刚好分布在单位圆上.

在实际情况下因采样快拍数的限制无法获得理想的自相关矩阵,因此实际得到的近似协方差矩阵存在一定程度的误差,这时只需要对式(24)求近似的q个根即可(最靠近单位圆).在均匀直线阵下的DOA估计角度计算式为

(25)

2.3 本文提出的低复杂度求根运算方法

在Root-MUSIC算法中,时间复杂度最高的运算为多项式求根,求根运算的时间复杂度为O(2M-1)3,可知,随阵元规模的扩大,其时间复杂度迅速增加,实验也证实了这一点.当阵元规模达到32以上时,其耗时已经超过了谱峰搜索算法,因此本节将提出一种降低求根多项式阶数的低复杂度求根运算方法.

根据2.2节可知,设存在多项式p1(z),p2(z),…,pM-q(z)且同时拥有q个公共根z1,z2,…,zq,则它们之间存在线性组合Q(z):

(26)

接下来证明式(26)的合理性,若r为多项式p1(z),p2(z),…,pM-q(z)的公共根,则r为多项式p1(z),p2(z),…,pM-q(z)任意线性组合的根.

设存在c1,c2,…,cM-q使得Q(z)=∑cipi(x),则将公共根r代入Q(z)可得

Q(r)=∑cipi(r)=∑ci×0=0 .

(27)

通过式(27)证明了式(26)的合理性,构造多项式p1(z),p2(z),…,pM-q(z)的线性组合Q(z),对应系数设为α1,α2,…,αM-q,令αM-q=1,可以得到如下表达式:

Q(z)=α1p1(z)+…+αM-q-1pM-q-1(z)+

pM-q(z) .

(28)

由式(22)可知：

αipi(z)=αi(βM-1,izM-1+βM-2,izM-2+…+

β0,i).

因此,Q(z)可以被重写为

Q(z)=(α1βM-1,1+…+αM-q-1βM-1,M-q-1+

βM-1,M-q)zM-1+(α1βM-2,1+…+αM-q-1×

βM-2,M-q-1+βM-2,M-q)zM-2+(α1β1,1+…+

αM-q-1β1,M-q-1+β1,M-q)z+(α1β0,1+…+

αM-q-1β0,M-q-1+β0,M-q) .

(29)

由式(27)可知,将公共根r代入Q(z)的前M-q-1阶,即令式(29)中zM-1,zM-2,…,zq+1对应阶系数为零,得到如下方程组:

(30)

式(30)可以写成如下的矩阵形式:

(31)

将式(31)简化整理为

x=-H-1h.

(32)

通过式(32)求得x,得到了对应式(26)中的系数αi,1≤i≤M-q-1.要得到x,就必须得到矩阵H及向量h.

在2.2节的讨论中,对自相关矩阵R进行特征值分解得到噪声子空间Un,其列向量uni=[β0,i,β1,i,…,βM-1,i]T,则将噪声子空间列向量翻转可得：

(33)

(34)

(35)

依据式(32)及式(35)计算得到线性组合Q(z)的系数向量x,并将式(21)代入式(25)得到一个新的求根多项式:

(36)

式中:αM-q=1.通过式(36)所构造的新求根多项式来替代式(24)完成求根运算,寻找最靠近单位圆的q个根,最终通过式(25)求得DOA估计角度.

本文算法的整体实现过程如下:

1) 根据式(12)及式(13)构造近似数据观测矩阵下的Hermitian Toeplitz矩阵Rtp;

2) 对矩阵Rtp做特征值分解,得到噪声子空间矩阵Un;

3) 依据式(33)～式(35),对Un翻转拆分获得矩阵H及向量h;

4) 通过计算式(32)得到系数向量x,并根据式(36)构造新的求根多项式Q(z);

5) 完成求根运算,并寻找最靠近单位圆的q个根,通过式(25)求得DOA估计角度.

3 实验验证与比较

针对改进算法的DOA估计性能及算法运行耗时进行仿真测试.为验证本文算法的各项性能,仿真实验分别对比了文献[4]的快速Root-MUSIC算法、文献[6]的无需特征分解计算的两步Root-MUSIC算法和经典Root-MUSIC算法;同时,对几种算法的时间复杂度进行了分析.

均方根误差计算式为

(37)

3.1 DOA估计精度实验

实验1:两个独立入射信源角度分别为-10°和10°;阵列为8阵元均匀直线阵列,阵列间距为λ/2;采样快拍数为200,信噪比以2 dB间隔从-10 dB变化至10 dB，进行1 000次蒙特卡洛实验，计算均方根误差.不同信噪比下算法的DOA估计精度对比如图1所示.

图1 不同信噪比下算法的DOA估计精度对比

由图1可知,随着信噪比的增加,4种算法的均方根误差均有所下降,而本文算法的估计误差最接近经典Root-MUSIC算法,低于文献[4]及文献[6]所提算法;同时可以发现,本文算法及经典算法更加接近算法的无偏估计量下界(克拉美罗下界)CRLB,拥有更高的估计精度.

实验2: 独立入射信源数量为3个,分别为-10°,10°和20°,其余实验条件与实验1一致,同样使信噪比以2 dB间隔从-10 dB变化至10 dB进行1 000次蒙特卡洛实验,计算均方根误差.

不同信噪比下算法的DOA估计精度对比如图2所示.

由图2可知,当入射信源增加为3个时,由于阵列规模的限制,文献[4]和文献[6]的算法已经基本失效,而本文算法及经典Root-MUSIC算法在低信噪比环境下效果较差,随信噪比的增加,均方根误差都有较大幅度的下降.受限于阵列规模,相较于2个入射信源的情况本文算法及经典Root-MUSIC算法的DOA估计性能均有所下降.

图2 不同信噪比下算法的DOA估计精度对比

实验3: 2个独立入射信源角度为-10°和10°,阵列仍采用8阵元均匀直线阵列,阵列间距为λ/2,信噪比固定为5 dB.采样快拍数以50为间隔从100变化至600进行1 000次蒙特卡洛实验,计算均方根误差.不同快拍数下算法的DOA估计精度对比如图3所示.

图3 不同快拍数下算法的DOA估计精度对比

由图3可知,采样快拍数对算法的DOA估计精度有显著影响,随采样快拍数量的增加,4种算法的均方根误差都呈现出减小的趋势.但从仿真图中可以清晰看出,本文算法的整体均方根误差曲线在文献[4]和文献[6]之下,其均方根误差曲线更加接近经典Root-MUSIC算法及无偏估计量下界CRLB,拥有更高的DOA估计精度及更强的鲁棒性.

实验4:独立入射信源数量变为3个,分别为-10°,10°和20°,其余实验条件与实验3一致.采样快拍数同样以50为间隔从100变化至600进行1 000次蒙特卡洛实验,计算均方根误差.不同快拍数下算法的DOA估计精度对比如图4所示.

图4 不同快拍数下算法的DOA估计精度对比

由图4可知,算法的均方根误差均随着快拍数的增加而下降.但在入射信源个数增加至3个时,由于阵列规模的限制,文献[6]由于算法自身需要拆分阵元进行互相关矩阵重构,其可利用阵元数量严重下降,DOA估计精度损失严重;文献[4]则由于算法本身利用噪声子空间信息较少,DOA估计精度同样不足;本文算法及经典Root-MUSIC算法的均方根误差相较于2个入射信源时同样有所上升,但总体更加趋近于无偏估计量下界CRLB,算法的DOA估计精度更高.

3.2 时间复杂度分析

本文算法的目的旨在解决目前多数低复杂度Root-MUSIC算法的DOA估计精度不足问题,那么所提算法自身的时间复杂度也不能有明显提升,对3.1节中的4种算法的时间复杂度进行分析(其中T为采样快拍数).算法各主要运算部分时间复杂度对比如表1所示.

表1 算法各主要运算部分时间复杂度对比

由表1可知,文献[4]算法的求根运算时间复杂度最低;本文算法在求根多项式运算部分的时间复杂度上高于文献[4],但低于文献[6]及经典Root-MUSIC算法;本文算法的协方差矩阵运算时间复杂度最低;在阵元规模较小时,本文算法及文献[4]算法的子空间提取时间复杂度略高于文献[6]算法,但随着阵元规模的增大,本文及文献[4]算法的子空间提取时间复杂度将低于文献[6]算法.整体看,本文算法的综合时间复杂度略高于文献[4]算法,但低于文献[6]及经典Root-MUSIC算法.下面通过实验5来验证这一结果.

实验5:两个独立入射信源角度分别为-10°和10°,信噪比固定为5 dB,采样快拍数为200,阵元间距为λ/2.阵元数量以2为间隔从6变化至40,每种算法运行1 000次计算耗时.算法运行1 000次的耗时对比如图5所示.

图5 算法运行1 000次的时间耗时对比

由图5可知,随着阵元数量的增加,算法的运行耗时均有所增加.其中文献[4]的算法耗时始终最低,这也说明了多项式求根运算是整个算法运行过程中耗时占比最高的部分;本文算法耗时高于文献[4]算法,但始终低于文献[6]算法及经典Root-MUSIC算法,从表1的时间复杂度分析也可以验证这一点.多项式求根运算的时间复杂度决定了算法整体的时间复杂度,而本文算法在不提高其他运算部分时间复杂度的同时,降低了多项式求根运算的时间复杂度,使得算法整体运行耗时降低.

4 结 论

针对多数低复杂度Root-MUSIC算法的DOA估计精度损失问题,本文提出了一种高精度低复杂度的改进Root-MUSIC算法.算法通过近似数据观测矩阵首行重构具有 Hermitian 性质的Toeplitz矩阵,并将其作为新的自相关矩阵做特征值分解运算得到噪声子空间;对噪声子空间作翻转拆分处理,重构新的求根多项式;最终通过求根运算得到DOA估计值.仿真实验结果表明,相较于部分前人算法,本文算法拥有更高的DOA估计精度,其均方根误差更接近于经典Root-MUSIC算法及无偏估计量下界,同时拥有更高的鲁棒性及稳定性，其时间复杂度则基本不高于前人算法.