基于协方差矩阵拟合的宽带无网格波达角估计

杨益新，张亚豪，杨 龙

（1.西北工业大学航海学院，陕西西安 710072；2.陕西省水下信息技术重点实验室，陕西西安 710072）

0 引言

宽带波达角（Direction of Arrival,DOA）估计是声呐系统阵列信号处理中一个非常重要的研究内容，目前比较流行的两类算法是非相干子空间法（Incoherent Signal-subspace Method,ISM）[1]和相干子空间法（Coherent Signal-subspace Method,CSM）[2]。ISM 将阵列接收信号通过离散傅里叶变换（Discrete Fourier Transform,DFT）划分至各个子带，在各子带内分别进行窄带DOA 估计；CSM 通过聚焦矩阵将各子带的信号转换至聚焦频率上，DOA 估计问题从而转换至在聚焦频率上的窄带DOA 估计。与ISM相比，CSM 在低信噪比（Signal to Noise Ratio,SNR）情况下性能更加稳健。

CSM 的核心是聚焦矩阵的设计，聚焦矩阵的性能直接影响DOA 估计的精度。文献[2-4]给出了聚焦矩阵不同的构造方法以降低聚焦误差。这些算法需预先获得DOA 信息作为聚焦角度来设计聚焦矩阵。因此，CSM 一般采用迭代的方式，以上一次迭代中的DOA 估计值作为下一次迭代的聚焦角度。然而，当聚焦角度与真实DOA 之间存在偏差时，将造成较大的聚焦误差。为提升聚焦矩阵对聚焦角度选取的稳健性，陈洪光等[5]将阵列流形中的角度变量和频率分离，仅对频率部分进行聚焦，避免了聚焦角度的选取，然而该方法产生的聚焦误差较大；Sellone[6]提出一种稳健聚焦矩阵，该算法可使聚焦角度周围扇区（即聚焦区域）的聚焦误差达到最小，因此无需得知DOA 的精确信息，仅需DOA 大致分布的先验信息即可，从而增强了对聚焦角度估计偏差的稳健性。

完成各子带的协方差矩阵聚焦后，现存方法一般采用子空间类算法如多重信号分类（Multiple Signal Classification,MUSIC）算法[7]和根 MUSIC（root MUSIC,rootMUSIC）[8]完成DOA 估计。然而，子空间类算法对SNR 要求较高。稀疏重构类算法[9-10]是近十几年内发展起来的高分辨DOA 估计算法，相较于子空间类算法，这类算法可用于SNR更低的情况[9]，因此受到广泛的关注。大部分稀疏重构类算法在DOA 估计时将空间划分为离散的网格，一般假设信号分布在该网格上。然而，当信号偏离网格时，DOA 估计精度会有所下降。为避免该问题，近年来针对直线阵阵形发展了许多无网格类稀疏重构算法，可分为原子范数类算法[11]和协方差矩阵拟合类算法[12]。其中，原子范数类算法利用原子范数作为稀疏惩罚项，可看作为l1范数类算法向连续空间的推广。然而，这类算法需要选取合适的正则参数控制数据拟合误差和解的稀疏性之间的平衡；协方差矩阵拟合类算法，如稀疏与参数方法（Sparse and Parameter Approach,SPA）[12]，在协方差矩阵拟合准则下利用采样协方差矩阵估计真实协方差矩阵，对所估计的协方差矩阵进行范德蒙分解（Vandermonde Decomposition,VD）[13]以获取阵列流形的估计值，并从中得到DOA 估计值。相较于原子范数类算法，协方差矩阵拟合类算法无需选取正则参数，因此更易用于实际的信号处理中。

本文将SPA 算法推广至宽带DOA 估计中，提出一种基于 CSM 的改进 SPA（Coherent Signal-subspace based Modified SPA,C-MSPA）算法，以在低SNR 情况下获得高DOA 估计精度。算法首先采用稳健聚焦矩阵[6]将各子带的信号投影至聚焦频率上，而在设计聚焦矩阵时则利用了信号空间分布的先验信息。为提升DOA 估计对聚焦矩阵输出的稳健性，本文基于频率选择的VD（Frequency-selective VD,FS-VD）原理[14]对协方差矩阵拟合准则进行改进，使得在新准则下重构的协方差矩阵中包含的DOA 信息严格限制在聚焦区域上，从而在DOA 估计时也有效利用了DOA 分布的先验信息。通过仿真和湖上实测数据分析可以看出，与基于ISM 的SPA（Incoherent Signal-subspace based SPA,I-SPA）算法和其他基于CSM 的无网格类算法相比，C-MSPA算法实现了高空间方位分辨能力且提高了DOA 估计精度。

1 信号模型建立

假设K个非相关信号从θ=[θ1θ2…θK]T入射至M元均匀线列阵（Uniform Linear array,ULA），信号的中心频率为fc，带宽为H。将阵列接收信号划分为N段，在每一段内进行L点的DFT，从而将时域信号转换至频域。在第l个子带上，指向角度θ的阵列流形为

式中：fl和λl分别表示第l个子带的中心频率和对应的波长；d为阵元间距。第l个子带上阵列接收信号可建模为

式中：xl(n) ∈ℂM×1，sl(n) ∈ℂK×1以及el(n) ∈ℂM×1分别表示第n段阵列接收信号、目标信号以及环境噪声在第l个子带上的DFT 系数；Al(θ)=[a(f l,θ1) …a(f l,θK)]∈ℂM×K表示阵列流形矩阵；N为频域快拍数。

假设噪声为均匀高斯白噪声，且与信号不相关，则第l个子带上阵列接收信号的协方差矩阵Rl∈ℂM×M为

式中：plk和σl分别为第l个子带上第k个信号功率和噪声功率；IM表示维度为M×M的单位矩阵；diag(x)表示以向量x为对角线元素的对角矩阵；E[⋅]表示期望算子。在实际信号处理中，真实的协方差矩阵一般难以获取。因此，采用采样协方差矩阵来近似。采样协方差矩阵的计算公式为

2 算法原理

2.1 稳健聚焦矩阵设计

聚焦矩阵的设计原理是将频率f l（l=1,…,L）上的信号子空间转变为聚焦频率f0上的信号子空间，则第l个子带上聚焦矩阵的设计问题表示为[3]

2.2 修正稀疏与参数方法

SPA 算法的原理是在协方差矩阵拟合准则下基于采样协方差矩阵重构真实协方差矩阵。聚焦后的协方差矩阵见式（8）。第i次迭代中聚焦后的采样协方差矩阵表示为

基于此，SPA 算法表示为[12]

从A0(θ(i))中可很容易获得第i次迭代中的DOA 估计值θ(i)。

3 仿真分析

下面进行仿真实验。仿真中考虑阵元数为5 的ULA，阵元间距为0.3 m，两个中心频率为2.5 kHz、带宽为1 kHz 的宽带远场目标信号入射至阵列。目标1 的DOA 为80.3°，目标2 的DOA 为95.2°。阵列的采样率为104 kHz。将接收信号均匀划分为10段，即N=10，每一段包含10 400 个快拍。每一段内进行10 400 点的DFT，将时域信号转换至频域，即频率分辨率为10 Hz。分别采用C-MSPA 算法、基于 CSM 的 SPA 算法（Coherent Signal-subspace based SPA,C-SPA）、基于 CSM 的 rootMUSIC 算法（Coherent Signal-subspace based rootMUSIC,C-rootMUSIC）以及I-SPA 算法进行DOA 估计。为了保证聚焦矩阵对聚焦角度的稳健性，各基于CSM 的算法中聚焦矩阵均由式（12）计算。迭代初始的聚焦角度由CBF 算法给出，聚焦区域的初始宽度B设定为 CBF 的 3 dB 波束宽度。聚焦频率设定为2.5 kHz，迭代终止条件设定为ε=0.01。

分别从DOA 估计精度以及分辨概率两方面对各算法性能进行比较。其中，DOA 估计精度由均方根误差（Root Mean Square Error,RMSE）来衡量，计算公式为

则称两个目标被成功分辨[17]，式中的θ1和θ2分别表示目标1 和目标2 的真实DOA。

图1 和图2 分别给出各算法的RMSE 和分辨概率随SNR 的变化曲线。由图1、2 可以看出，C-MSPA 算法在SNR 小于−10 dB 时的DOA 估计精度和分辨概率都优于其他算法。主要原因是C-MSPA 算法在进行DOA 估计时利用了信号分布的先验信息，增强了其对噪声的稳健性。相较于传统的高分辨算法，稀疏重构类算法对SNR 的要求较低，C-MSPA 和C-SPA 的DOA 估计精度在SNR小于0 dB 时优于C-rootMUSIC，同时分辨概率在 SNR 小于−7.5 dB 时也高于C-rootMUSIC。I-SPA算法在各个子带内独立进行DOA 估计，并将各子带的DOA 估计结果进行平均。若某一子带的估计结果偏差较大，对该算法的估计精度会有很大影响。因此，与基于CSM 的算法相比，I-SPA 算法对SNR 的要求较高。当SNR 大于0 dB 时，I-SPA 算法的性能才慢慢接近其他算法。

图1 各算法DOA 估计的RMSE 随SNR 的变化 Fig.1 RMSE of each method versus SNR

图2 各算法的分辨概率随SNR 的变化 Fig.2 Resolution probability of each method versus SNR

为了对比算法的空间方位分辨能力，固定目标1，将目标2 从85.2°，以2°为步长移动至99.2°，即两个目标之间的角度间隔从5°变化至19°。SNR设定为−5 dB，其他仿真条件不变。

图3、4 分别给出了RMSE 和分辨概率随角度间隔的变化曲线。由图3、4 中可知，C-MSPA 算法在目标角度间隔小于9°时的RMSE 和分辨概率均优于其他算法，说明当DOA 估计时考虑信号分布的先验信息有效提升了算法的空间方位分辨能力。由于稀疏重构类算法利用信号在空间分布的稀疏性，获得了更高的空间方位分辨能力，因此C-MSPA 和C-SPA 算法在目标小角度间隔情况下的RMSE 和分辨概率均优于C-rootMUSIC。由图1、2 可知，I-SPA 算法对SNR 的要求较高。由于在该仿真条件下设置的SNR 较低，导致I-SPA 算法在所有考虑的角度间隔情况下RMSE 均大于其他算法，同时其分辨概率也小于1。由此也可说明将各子带相干处理相较于非相干处理可有效提升算法性能。

图3 各算法的RMSE 随角度间隔的变化 Fig.3 RMSE of each method versus angle separation

图4 各算法的分辨概率随角度间隔的变化 Fig.4 Resolution probability of each method versus angle separation

4 实验分析

为了验证C-MSPA 算法的性能，对实验数据进行处理和分析。实验数据来源于2016 年在新安江水库试验站进行的湖上实验，湖深约为60 m。实验设备布放如图5 所示。

图5 实验设备布放图 Fig.5 The layout of experimental equipment

实验中，接收阵为一条自制的5 阵元ULA，阵元间距为0.3 m，水平放置在水下。实验采用两个声源，其中声源1 为低频宽带声源UW350，声源2为溢流环。根据几何关系可以算出，声源1 的方位角大致为75°，声源2 的方位角大致为57°。两个声源同时发射中心频率为2.5 kHz、带宽为1 kHz 的宽带信号，持续时间为17 s。接收阵列以104 kHz 的采样率采集信号。将接收到的信号划分为17 组，每一组持续时间为1 s。将每组信号划分为10 段，即N=10，在每一段中进行10 400 点的DFT，即频率分辨率为10 Hz。

采用第3 节仿真中的算法对每一组数据进行DOA 估计，其余的参数设置与仿真中一致，DOA估计结果如图6 所示。考虑到CBF 算法的稳健性，图6 同时给出了CBF 的估计结果作为参考。图6中背景色即为CBF 的估计结果，三角符号则为各无网格算法的估计结果。

从图6 中可以看出，由于CBF 算法的分辨能力较低，无法分辨出两个目标，但从结果中可以看出，在区间[50，80]°存在目标。C-MSPA 算法在DOA 估计时利用了信号大致分布的先验信息，该算法可以很好地估计出两个目标。C-SPA 算法在DOA 估计时则没有考虑信号的分布情况，其在某些时刻的DOA 估计值存在较大偏差。与C-MSPA和C-SPA 算法将各子带进行相干处理不同，I-SPA算法是在各个频点上进行非相干DOA 估计。由于部分子带内的估计结果存在较大偏差，导致该算法无法成功分辨出两个目标。相较于稀疏重构类算法，C-rootMUSIC 算法对SNR 要求较高，同时空间方位分辨能力弱于稀疏重构类算法，该算法也无法很好地估计出两个目标信号。通过对实验数据的处理，验证了所提出的算法在宽带DOA 估计中的优异性能。

图6 各算法的DOA 估计结果 Fig.6 DOA estimation results of each method

5 结论

本文将协方差矩阵拟合类无网格DOA 估计算法推广至宽带情况。算法首先利用聚焦矩阵将各子带的数据转移至聚焦频率上，以进行相干DOA 估计。不同于已有协方差矩阵拟合类算法，本文对协方差矩阵拟合准则进行改进，以利用DOA 分布的先验信息，使在该准则下重构的协方差矩阵中所包含的DOA 信息严格限制在聚焦区域内。仿真和实验数据处理结果表明，与基于ISM 的无网格算法相比，本文所提出的算法通过聚焦矩阵对各个子带进行相干处理；相较于其他基于CSM 的无网格类算法，所提出的算法在DOA 估计时利用了DOA 分布的先验信息，实现了较高的空间方位分辨能力且提高了DOA 估计的精度。