冯 成，龚晓峰，雒瑞森

(四川大学电气信息学院，四川 成都 610065)

1 引言

阵列信号处理[1]中的空间谱算法可实现波达角精确估计，一直以来是学者们研究的热点。测向技术主要有相关干涉仪测向[2]，TDOA时间差测向[3]，MUSIC类算法[4]和ESPRIT类算法[5]。其中，MUSIC 算法因分辨率高、稳定性良好，具有普遍的适用性，工程应用最为广泛。

MUSIC算法是利用噪声空间与信号空间的正交性实现对信号入射角估计。因其需要矩阵特征分解，谱函数计算和峰值搜索，涉及大量的复数运算，使得人们在应用MUSIC算法的过程中不得不增大扫描步进以减少算法计算量。

为降低MUSIC算法的运算复杂度，近来学者们提出了一些方法。文献[6]利用子空间投影快速构建信号子空间，文献[7]通过对一维均匀线阵的构造和重排得到等价噪声子空间。 这些方法避免了矩阵的特征分解，而谱曲线计算的运算复杂度仍然很高。为此，文献[8]通过变量代换将谱函数对应的复系数多项式转换为实系数进行求根，从而减少求根所需计算量。文献[9]利用部分噪声子空间进行变步进谱峰搜索。文献[10]提出先通过快速傅里叶变换(FFT)实现波束成形预估角度范围，进而在预测范围内用传统MUSIC方法进行步进扫描，从而避免了全空域的谱峰搜索。更进一步的，文献[11]在FFT变换预估角度的基础上，通过线性调频变换来计算谱函数值。然而，这些算法对计算量降低效果有限，复杂无线电环境下预估角度不准，导致测向失败率上升。并且，各种已有改进算法的核心思想依旧局限在步进扫描和峰值比较的思路中。

MUSIC算法谱峰搜索的关键在于寻找谱函数极值点对应的角度值，这促使作者思考从迭代求极值的角度出发，选取合适的目标函数，经多次迭代使自变量收敛至极值点，从而完成DOA估计。黄金分割法可有效求解任意函数的极值点，然而，它需要合适的搜索区间作为迭代初始值。文献[10-11]用FFT预处理信号向量阵得到角度预估范围的方法在信号夹角变小时会失效，本文在此基础上提出先对噪声阵列向量补0，再FFT变换、各列变换值累加的方法以提高角度预估分辨能力。本文提出的结合FFT和黄金分割的快速DOA估计算法，首先，可以灵活选择FFT变换点数调整角度预估能力，适应各种场景。其次，以构造的伪谱函数为目标函数做黄金分割，用利用迭代收敛快速精确的求解波达角，从而在保证精度的前提下避免了步进扫描繁冗的计算。而且，随着天线阵列数和信号数的增加，算法复杂度只会略微增加。

2 阵列信号模型和MUSIC算法原理

波长为λ，互不相干的M(M

X(k)=A(θ)S(k)+N(k)

(1)

其中A(θ)=[a(θ1)，a(θ2)，…，a(θN)]是阵列导向矢量矩阵[12]，S(k)是入射信号，是噪声信号。均匀直线阵的方向响应向量a(θi)如式(2)

(2)

其中λ是入射信号中心波长，d是阵元间距，θi是信号与法线的夹角。

阵列接受数据的协方差矩阵记为Rxx，经Jacobi旋转[13]或QR算法[14]特征分解得到信号子空间Us和噪声子空间Un，如式(3)所示

Rxx=E[X(k)X(k)H]=UΣUH

(3)

其中Σs是信号特征值构成的对角阵，Σn是噪声特征值构成的对角阵。MUSIC算法指出阵列方向响应向量和噪声子空间正交，既

a(θi)Un=0;i=1，2，…，M

(4)

实际数据受噪声影响，上式不精确为0。故取矩阵乘法后的向量模值的倒数，再取对数为谱函数值，如式(5-6)所示。谱曲线峰值所对应的角度便是来波方向的估计值。

(5)

Pmusic=-10log10(c2+d2)

(6)

若用经典MUSIC算法完成谱曲线计算，以9元直线阵为例，[-90° 90°]范围内以1°为扫描步进，式(5-6)合计需要180*3次三角函数计算，180*91次复数乘法，180*80次复数加法，180次对数运算。至于圆阵，因需要在[0° 360°]范围扫描，计算量翻倍。当阵元数N增加或者需要更高精度时，计算量也会成倍增加。可见降低运算复杂度对MUSIC算法十分重要。当然，若增大扫描步进，计算量将降低相应倍数，但同时牺牲了测向精度，这也是工程应用中常采用的折中办法。

3 基于FFT和黄金分割的新型DOA

本文所提算法旨在提高测向精度的同时降低MUSIC算法的运算复杂度。基于FFT和黄金分割的MUSIC算法分为两步：首先，用FFT变换预处理噪声向量矩阵，得到角度预估范围；然后，以预估的角度范围作为迭代初始值，用黄金分割法求预估范围内的伪谱函数极小值点，即为波达角的精确估计值。

3.1 FFT粗估计DOA

空间谱估计的核心思想在于阵列导向矢量a(θ)与噪声向量Un的正交性，两者向量点积的模值应当接近0。构造伪谱函数P(θ)如式(7)所示

(7)

其中R为噪声向量阵的总列数。均匀直线阵的导向矢量a(θ)具有范德蒙结构，将式(2)带入化简式(7)可得

(8)

其中Ui(n)为噪声向量阵的第i列，第n行对应的元素。令f=-dsinθ/λ，代入式(8)有

(9)

易发现式(9)中含有Ui的离散傅里叶变换(DFT)的表达形式，使得可以考虑用快速傅里叶变换(FFT)计算式(9)在一系列指定f值对应的伪谱值。为此，需要将相位均分，令f=l/N，代入式(9)有

(10)

用FFT变换得到的N个伪谱值点，相当于将[-90° 90°]的空域范围分成了N个波束，其中第l个波束对应的中心角度θl为

(11)

这里需要特别指明：当l的增加使反三角函数的自变量超出定义值时，需要先根据傅里叶变换的周期性对相位进行折算，再代入反三角函数。此外，N的取值不一定非得选择天线阵元个数，根据FFT的数学定义可知，可以在噪声列向量后添加任意个0，再取FFT变换以提高此预估方法的分辨率。因添加的是0值，故这只需付出极少的计算代价。

虽然相邻两波束的相移量是等间隔的，但每个波束的中心指向却是不等间隔。表1给出了做32点FFT变换，d=λ/2时，用本文所提方法预处理噪声向量阵时各波束对应的子空域角度值。

表1 各波束对应空域角度

比较FFT得到的伪谱值的大小，若某个谱值同时小于其两侧的谱值，且相邻谱值差的绝对值大于2(筛掉假极小值点)，则认为其对应角度θl是波达角的粗略估计值，选取其相邻两点对应的角度值作为下一步精确估计DOA的搜索范围。

以一个信号源，入射角47°为例，图1是对噪声向量阵分别做16点和32点FFT变换所得的一系列伪谱值对比。理论分析和仿真都表明不同点数的FFT变换相当于对空域做不同比例的均分，可以看做是伪谱函数p(θ)在不同采用率下得到的一系列采样点。对于32点FFT，20号点是符合条件的极小值点，查表格1知应选取[43.43° 54.34°]作为下一步黄金分割的搜索范围。

图1 FFT预处理所得的伪谱值曲线

3.2 黄金分割法精确估计DOA

黄金分割法是建立在区间消去法基础上的试探方法，对函数f(x)在搜索区间[ab]内适当插入两点x0和x1，如式(10-11)所示。计算对应函数值f(x0)和f(x1)，如图2所示。比较两者大小，按式(12-13)更新区间。

x0=0.382*(b-a)

(12)

x1=0.618*(b-a)

(13)

iff(x0)>=f(x1)；a=x0；b=b；

(14)

iff(x0)

(15)

完成一次计算后，以新区间[ab]进行下一次黄金分割点插值，如此多次迭代，直到满足条件b-a<ε(ε是一个预设的较小值)时，退出迭代，并取x=(a+b)/2作为极小值点的估计值。

图2 黄金分割法示意图

本文将伪谱函数p(θ)做为黄金分割的目标函数，则波达角是使p(θ)取接近0值的极小值点，FFT预处理噪声阵的预估角度范围是迭代初始值，多次迭代后，自变量θ便会收敛至波达角的精确估计值。选取p(θ)做黄金分割的目标函数是因为其极小值是接近0的值，可以辅助验证极值点求解的正确性。另外，相比选择传统的谱函数值，如式(6)的Pmusic(θ)，可以避免一些不必要的对数运算。

由2.1节中的示例信号绘制的p(θ)曲线如图3所示，上一节预估的角度范围[43.43° 54.34°]是迭代初始值，取ε=0.001，MATLAB仿真有：黄金分割法经15次迭代后，θ收敛至47.0092°。这说明了本文算法的正确性。

图3 伪谱函数曲线图

4 算法复杂度分析

本节重点从理论上对比分析经典MUSIC和本文所提MUSIC算法的计算量。当天线阵元数取9，FFT变换点数取32时，两种方法的谱峰搜索模块对应计算量汇总于表格2。

表2 计算复杂度对比

其中K表示经典MUSIC需要计算的谱函数点数，取决于扫描精度。I表示黄金分割的迭代次数，后文仿真结果可知，黄金分割MUSIC求单个极值点的迭代次数总是迭代15次左右后收敛，远远小于传统谱峰搜索方法需要计算的谱值点数K。80次复数乘法和160次复数加法是32点FFT变换所需计算量。

从表格1可直观的对比出本文的新型DOA算法与经典MUSIC算法的运算复杂度比值近似为(I+1)/K。I的大小依赖于信号源个数m，因为m个信源意味着m次黄金分割，可以认为I=15*m。而用经典MUSIC的方法以1°为步进扫描直线阵对应的空域[-90° 90°]，需要扫描180个点，即取K=180。有m=1时，新算法复杂度只有经典MUSIC的8.89%，m=2时，计算复杂度是原方法的17.22%。

综上，改进算法先通过FFT变换快速计算一系列特定角度点对应的伪谱值，预估出角度范围，再利用黄金分割迭代逼近DOA。避免了全空域的谱函数值的计算，而且完全规避了峰值搜索过程，计算量大幅度减少的同时还可以提高测向精度。

5 仿真验证

本节通过MATLAB仿真验证所提算法的准确性和收敛性。仿真参数设置如下：均匀直线阵阵元数N=9，阵元间距d=0.5λ，快拍数n=1024，噪声是加性高斯白噪声，FFT变换点数选取32，黄金分割法退出迭代的预设阀值ε=10-3。

实验1：选定两个信号以±27.5°入射，信噪比snr=10dB，图4是本文所提的FFT变换预处理噪声向量阵而得的一系列伪谱值。通过比较此32个伪谱值大小，可得7号点和25号点是符合条件的极小值点，查表格1知，应选取[-30°-22.02°]和[22.02° 30°]作为角度预估范围。

图4 实验1 FFT预处理

图5是根据式(7)绘制的伪谱值函数P(θ)与角度θ的函数曲线图，即黄金分割法的目标函数。

图5 实验1伪谱函数曲线

图6 实验2 FFT预处理

实验2：增加信号源数量，减小入射信号夹角，降低信噪比。选定3个信号源以(-47.5°，42°，63.5°)入射，信噪比snr=0dB。图6为FFT预处理得到的曲线图，图7是伪谱值P(θ)与角度θ的函数曲线图。用相同的方法选取极小值点和查表1选取角度搜索范围。

图7 实验2伪谱函数曲线

MUSIC算法改进前后，两次实验的估计角度θ，仿真时间T(仅指从噪声子空间到DOA估计的计算时间)汇总于表格3。实验显示黄金分割法迭代15次左右收敛，验证了之前的复杂度分析，相应仿真时间也降低至预计值。

表3 MUSIC改进前后DOA精度仿真时间对比

实验3：为了比较算法改进前后的DOA精度做均方根(RMSE)误差实验[15]，均方根误差计算如式(13)所示

(16)

其中θi为波达角预设值，i为算法估计值，num为独立仿真次数。实验选取两个信号源，分别在[-70°-20°]和[20° 70°]范围内取随机角度入射9元直线阵，信噪比范围为[-10dB 20dB]，步进为2dB。不同信噪比下独立仿真100次计算均方根误差。

图8给出了经典MUSIC算法和黄金分割MUSIC算法的均方根误差随信噪比变化的曲线。由图可得，在信噪比低端两者DOA估计精度略有差距，随着信噪比的提高，新型DOA估计算法的RMSE明显小于以1°为步径的传统MUSIC算法，此时迭代达到了更高精度扫描的效果。

图8 两种方法估计精度对比

从理论上分析，由于噪声子空间的构造方式一致，所以经典MUSIC算法与基于FFT和黄金分割的MUSIC算法对应谱函数的极值点也应该一致，差异在于本文所提算法是用黄金分割法迭代逼近极值点，从而使得在计算复杂度降低的情况下，测向精度反而提升。

6 结束语

本文以降低MUSIC算法谱峰搜索模块运算复杂度为目的，提出一种结合FFT和黄金分割的快速DOA估计算法，解决工程实践在测向精度和实时性两者间做取舍的困境。改进DOA算法主要优点有：1：谱函数模块的运算复杂度降低至经典方法的17%左右；2：改进算法测向精度比1°步径扫描的经典MUSIC算法更高；3：在多个信号源同时存在，低信噪比下也能成功预估角度并实时测向；4：改进算法将迭代求极值的思想引入以规避传统方法的峰值搜索，为无线电信号处理提供了一些新的数学思路。黄金分割法并不依赖于特定的天线阵列形状，而FFT预处理噪声阵的方法却只适用于均匀直线阵，寻求复杂无线电环境下任意阵列形状的角度预估方法和更加精简的迭代算法是迭代求解波达角的下一步研究方向。