基于变换域APSO的任意阵列宽带DOA估计算法

刘学承 朱 敏 武岩波

（1.中国科学院声学研究所海洋声学技术中心，北京 100190；2.中国科学院大学电子电气与通信工程学院，北京 100049；3.北京市海洋声学装备工程技术研究中心，北京 100190；4.中国科学院声学研究所声场声信息国家重点实验室，北京 100190）

1 引言

波达方向（Direction-of-arrival，DOA）估计是水下目标定位的关键技术［1-3］。为了降低水声信道的频率选择性、混响等带来的不利影响，具有大时宽-带宽积的宽带信号被广泛应用于水声通信定位领域［4-6］，因此适用于宽带信号的DOA 估计成为了一个研究重点。目前宽带DOA 估计算法大都假设信号源发射的信号波形未知［7］，发射信号波形已知（如有源声呐中的水声定位信号或通信同步信号）条件下的宽带DOA估计研究相对较少。

针对发射信号波形已知情况，Li 和Compton 等首先将信号波形信息应用到确定性极大似然（De⁃terministic Maximum Likelihood，DML）算法中提高了DOA 估计精度，随后Li 等提出解耦最大似然估计（Decoupled ML，DEML）算法来降低算法运算量，提出白噪声环境下的解耦最大似然估计（White De⁃coupled ML，WDEML）算法来改善DOA 估计的空间分辨率［8-11］。在单信号情况下，Hu 等［12］提出了基于波形已知参考信号的增强干涉仪DOA 估计算法，该算法本质上仍属于DEML 算法，但计算量显著低于后者。这些波形已知DOA 估计算法不足之处在于：一是建立在窄带信号假设上，不能直接用于宽带信号，比如宽带线性调频（Linear Frequency Modu⁃lated，LFM）信号；二是需要大采样快拍和多维最优化搜索才能有效估计DOA，计算量仍然很大；三是受到均匀阵型的限制，如均匀线列阵（Uniform Lin⁃ear Array，ULA），当某一阵元出现故障时算法失效。因此，适用于任意阵列且具有低计算量的宽带DOA估计算法成了国内外学者的研究重点。

Sun 等［3］结合双阵元时延差估计方法，提出了一种任意阵型水声定位算法，该算法将阵列划分成多组双阵元子阵，每组子阵都对应一组时延差估计等式，联合所有等式即可估计宽带信号的DOA，但需要在高采样率和高信噪比下才能有效估计DOA。Liu 等［13］结合匹配滤波方法，提出了一种基于阻尼牛顿迭代的宽带DOA 估计方法，该方法适用于小快拍数、低信噪比和任意阵列情况，缺点是需要良好的DOA 初始值和对数似然函数的梯度信息才能有效估计DOA。近年来，群体智能优化算法被用来对基于DML 准则的DOA 估计空间谱函数进行最优化搜索［14-19］，避免了DOA 预估计和梯度求取，有效降低了计算复杂度。其中应用最广泛的是粒子群最优化（Particle Swarm Optimization，PSO）算法［20］，该算法简单易实现，收敛速度快，收敛精度高，并且适用于少快拍甚至单快拍数据的情况［18-19］。为了进一步提高PSO 算法的单目标（单峰）搜索性能，Yang［20］提出了加速粒子群最优化（Accelerated Particle Swarm Optimization，APSO）算法。APSO 算法无需计算粒子运动速度，计算更加简单，并具有更快的全局收敛速度和收敛精度。因此，在单信号DOA 估计中，使用APSO 算法搜索空间谱函数最大值具有潜在的计算优势，但目前尚未有研究者利用APSO 算法优势来提高DOA估计性能。

为了提高宽带信号DOA 估计精度并降低计算复杂度，本文结合论文［13］中的匹配滤波处理方法，提出一种适用于发射波形已知的基于变换域APSO 的宽带DOA 估计算法，该算法适用于任意阵列。首先根据波形已知的发射信号对阵列数据进行时域匹配滤波处理和傅里叶变换，得到宽带频域阵列数据模型，其次根据DML 准则构建宽带DOA估计的二维空间谱函数，然后采用变换域APSO 算法搜索得到该空间谱函数的最大值，也即DOA 估计结果。该算法的核心思想主要在于将变换域APSO算法用于宽带DOA 估计中的多维最优化搜索，有效提高了DOA 估计精度并降低了计算复杂度，同时避免了DOA 预估计和空间谱函数求导等复杂计算，简化了算法过程，易于工程实现。仿真结果验证了所提算法的有效性。

2 阵列数据模型及其DOA估计理论界

2.1 波形已知的宽带阵列数据模型

考虑一个由M个相同阵元构成的任意阵列（如图1所示），第m个阵元的坐标为0，1，…，M-1，其中第0 个阵元（记为参考阵元）位于直角坐标系原点，即p0=0。假设一个中心频率为fc的远场宽带信号s(t)从二维角度ϑs入射到阵列上，ϑs=[φs，θs]T，其中φs为方位角，θs为俯仰角。信号来波方向ϑs在直角坐标系中对应的单位矢量为νs=[cosφscosθs，sinφscosθs，sinθs]T。

假设第m个阵元处接收到的信号（不考虑噪声）为

其中ρm、Tm分别为衰减幅度和信号到达时刻。式（1）可用参考阵元接收信号s0(t)=ρ0s(t-T0)表示为：

式中τm(ϑs)为信号到达第m个阵元时相对于参考阵元的时延差，即，c为信号传播速度，方便起见，后文将τm(ϑs)简写为τm。假设信号到达各阵元时衰减幅度相同，则有1。考虑噪声影响，第m个阵元输出为

nm(t)为xm(t)中的噪声分量。

由于阵列接收信号为带通信号，需要高采样率，从而使得数据量大大增加，若下降频到基带上进行处理，则可以降低采样率，有效减少数据量。假设发射信号s(t)的基带信号为同 理，和(t)分别为第m个阵元处接收信号sm(t)、阵元输出xm(t)和噪声分量nm(t)的基带信号。因此，式（3）对应的基带形式为：

由于(t)仍为宽带信号，无法直接在时域上将时延转化为相移以建立阵列数据模型，因此宽带阵列信号处理通常是在频域中进行的。

为了保证频域上每个频点处都有足够多的快拍，通常需要将阵列接收数据划分为多个不重叠的数据段，然后再对每段数据分别进行傅里叶变换。但分段傅里叶变换的方式只适用于整个观测时间内平稳的信道和频率与时间无关的信号，并不适用于时变的水声信道和宽带调频信号。针对该问题，假设在较短的观测时间内水声信道是平稳的，此时阵列采样数据数量较少，不做分段处理，直接对M个阵元输出进行L点快速傅里叶变换（Fast Fourier Transform，FFT），并写成矩阵形式：

2.2 波形已知的宽带DOA估计理论界

文献［21］给出了ULA 条件下发射波形已知的窄带DOA 估计的克拉美·劳界（Cramér-Rao Bound，CRB），将其扩展到宽带信号和任意阵列上，可得：

3 基于变换域APSO的宽带DOA估计算法

本节介绍了所提出的基于APSO 的宽带DOA估计算法，算法框图如图2 所示。首先介绍在单信号下基于DML 准则的宽带DOA 估计空间谱函数，然后介绍能够快速搜索该空间谱函数最优值的变换域APSO算法。

3.1 DML准则下的宽带DOA估计算法

根据论文［7］、［13］可知，式（9）在DML 准则下关于DOA的对数似然函数为：

式中ϑ=[φ，θ]T表示(φ，θ)空间（后文简称为ϑ空间）中的一个任意二维角度，φ∈[0，2π]、θ∈[0，]π/2，表示到a(ϑ，fl)所张成的列空间中的投影矩阵。L(ϑ)最大值对应的角度为信号方向ϑs的估计值，即将式（11）代入并化简，可得：

3.2 常规APSO算法和变换域APSO算法

假设APSO 算法中有N个粒子，迭代搜索K次，目的是搜索D维函数的最大值对应的解，即最大值在D维空间中的位置，则第k次迭代时粒子的位置更新公式为：

APSO 算法步骤主要分为两部分，首先是粒子位置初始化，将N个粒子均匀分布D维搜索空间中，获得然后迭代搜索，根据式（13）迭代更新全部粒子的位置和全局最优位置，第K次迭代后得到的即为函数F(g)最大值对应的解。在本文中空间谱函数F(ϑ)最大值对应的粒子位置即为DOA估计值。

使用APSO 算法在ϑ空间中搜索空间谱函数F(ϑ)的最大值时，τm为ϑ的二维非线性函数，且φ、θ具有不同取值范围，增加了APSO 算法的计算复杂度。而角度ϑ对应的单位矢量为ν=[νx，νy，νz]T，其中

表1 变换域APSO算法流程Tab.1 The flowchart of transformed-domain APSO algorithm

4 仿真实验与讨论

本节对所提算法进行仿真实验，仿真中假设深海垂直声信道下的多径可以忽略不计，且只有一个水下目标在ϑs=[44.5°，58.5°]T方向和水面母船进行通信。通信前导信号为对多普勒不敏感的宽带LFM信号，取其作为定位信号，该信号发射波形已知，对应的基带形式为，其中为调频斜率，带宽B=5 kHz，脉冲宽度T=51.2 ms，采样频率fs=10 kHz。接收阵列采用水声通信中常用的同心圆平面阵列，由19个相同阵元构成，相邻阵元间距为，其中发射信号中心频率fc=10 kHz，信号传播速度将阵列平面设置为直角坐标系xOy 平面，并将阵列几何中心放置在坐标系原点，如图3 所示。需要注意的是，这里为便于仿真，采用了具有规则几何特征的同心圆平面阵列，亦可换为其他任意阵列。

仿真中信噪比变化范围设置为-10 dB ∼20 dB，变化间隔为5 dB，每个信噪比下进行Q=500次独立的蒙特卡洛实验。定义二维DOA估计的均方根误差（Root Mean-Square Error，RMSE）为：

实验1为了说明APSO算法在搜索空间谱函数F(ν)最大值的优势，本实验对比了信噪比为0 dB时APSO 算法和PSO 算法在不同粒子数量下的全局收敛性能，包括收敛到空间谱函数最大值的收敛速度（即所需迭代次数，次数越小则收敛速度越大）、收敛精度（即最大值估计误差，误差越小则收敛精度越高），图4、图5分别为两种搜索算法的收敛速度、归一化收敛误差。由图4可知，在相同粒子数量下，APSO算法收敛速度明显大于PSO算法。由图5可知，在相同粒子数量下，APSO 算法的收敛精度也明显更高。因此，达到相同的DOA估计精度时APSO算法所需粒子数量和迭代次数均比PSO算法的要小，即APSO算法计算量要比PSO算法的小。并且图4、图5均表明粒子数量越大，APSO、PSO 算法的收敛速度越快，收敛精度越高。为了提高F(ν)最大值搜索精度并降低计算量，需要选择合适的粒子数量和迭代次数。经过大量仿真实验，发现APSO算法设置为30个粒子、40次迭代时具有合适的收敛精度和计算量。

实验2本实验主要考察所提算法的DOA 估计精度和计算复杂度。图6对比了已知发射波形条件下的宽带信号DOA 估计的CRB 和宽带DOA 估计算法在6种最优化搜索方法下的RMSE，表2对比了6 种算法的计算复杂度（不考虑6 种算法都包含有的匹配滤波处理和FFT 部分对应的计算复杂度）。这6种算法包括：

（1）基于单精度网格搜索的宽带DOA 估计算法，简记为Grid1 算法，其中空间网格间距取为1°，按其对应的弧度对ϑ空间划分出的空间网格数量为G1=20942；

（2）基于双精度网格组合搜索的宽带DOA 估计算法，简记为Grid2 算法，空间谱函数F(ϑ)最大值搜索时先在ϑ空间中进行全局粗略搜索，再进行局部精细搜索，全局搜索时空间网格间距取为5°，对应的空间网格数量为G2=887，局部搜索时空间网格间距取为0.1°，对应的空间网格数量为G3=5195；

（3）基于牛顿迭代搜索的宽带DOA 估计算法，其中DOA初始值通过单精度网格搜索获得［13］，将该算法简记为Grid-Newton 算法，为保证DOA 初始值靠近最大值，DOA 预估计时空间网格间距取为2°，对应的空间网格数量为G4=5313，此外牛顿迭代次数取为Knt=5；

（4）基于牛顿迭代搜索的宽带DOA 估计算法，其中DOA 初始值通过APSO 搜索获得，将该算法简记为APSO-Newton 算法，为保证DOA 初始值靠近最大值，DOA 预估计时APSO 粒子数量为N0=20、粒子位置迭代次数取为K0=15，此外牛顿迭代次数仍取为Knt=5；

（5）基于PSO 搜索的宽带DOA 估计算法，简记为PSO 算法，粒子数量取为N1=30，粒子位置迭代次数取为K1=40；

（6）基于APSO 搜索的宽带DOA 估计算法，简记为APSO 算法，粒子数量取为N1=30，粒子位置迭代次数取为K1=40。

需要注意的是，上述算法中的空间网格位置或粒子位置等同于二维空间谱中的谱点，且空间谱函数F(ϑ)或F(ν)的单个谱点值的计算复杂度为O(ML)。

图6 表明，基于网格搜索的Grid1、Grid2 算法受限于空间网格划分精度，RMSE 相对较大，而基于牛顿迭代的Grid-Newton 算法、APSO-Newton 算法与基于粒子群随机搜索的PSO 算法、APSO 算法都能够有效克服网格搜索的缺点，具有高DOA 估计精度，即RMSE 更趋近于CRB。其中APSO算法和Grid-Newton 算法、APSO-Newton 算法在不同信噪比下的DOA 估计误差均相同，意味着APSO 算法具有和牛顿类算法相当的DOA 估计性能。低信噪比下PSO和APSO 算法DOA 估计误差相同，信噪比为0dB 时RMSE 约为0.22°，高信噪比下APSO算法具有比PSO 算法的更小DOA 估计误差，信噪比为20 dB 时APSO算法RMSE约为0.02°，因此APSO算法DOA估计性能优于PSO算法。综上所述，所提的基于APSO搜索的宽带DOA估计算法具有高DOA估计精度。

通过比较单个谱点值的计算复杂度和表2 中6种算法计算复杂度可知，6 种算法可以通过比较等效谱点数量来衡量计算复杂度高低。表2 表明，包含网格搜索的算法（即Grid1、Grid2、Grid-Newton 算法）计算复杂度较高，而包含APSO 或PSO 的算法（即APSO-Newton、PSO、APSO 算法）计算复杂度较小，因此采用APSO 或PSO 搜索的算法具有明显的低计算量优势。APSO-Newton 算法计算复杂度和PSO 算法、APSO 算法的大致相当，但Newton 类算法依赖于空间谱函数的一阶、二阶导数和良好的DOA预估值，步骤复杂，而APSO 算法可以有效克服上述问题。综上所述，所提的基于APSO 搜索的宽带DOA 估计算法具有低计算复杂度，且无需DOA 预估计，步骤简单，易于工程实现。

表2 6种宽带DOA估计算法计算复杂度对比Tab.2 Comparison for computing complexity of 6 wideband DOA estimation algorithms

实验3本实验主要考察所提算法在不同阵列结构下的性能。由于DOA 估计性能与阵元数量、阵元间距有关，为了有效对比所提算法在不同阵列结构下的性能，需要保持不同阵列的阵元数量、相同相邻阵元间距的一致性。因此，仿真所选取的阵列均由19个阵元组成，且相邻阵元间距为除同心圆阵列（图3）外，其他阵列阵型分别如图7～图10所示，包括三角形阵列（19 个阵元按等边三角形分布）、方形随机阵列（5×5均匀矩形阵列，从其中随机选择19 个阵元）、伪立方阵列（3×3×3 立方体阵列，顶层只保留中心阵元）和圆柱阵列（均匀圆柱阵列，在底层圆心处添加一个阵元）。其中，三角形阵列、方形随机阵列以及同心圆阵列为二维平面阵列，伪立方阵列和圆柱阵列为三维立体阵列。需要注意的是，在满足阵元数量相同、相邻阵元间距相同的条件下，各个二维平面阵列的平面孔径（在xOy平面上的面积）大致相当，各个三维立体阵列对应的平面孔径也大致相当，但比二维阵列的要小。

图11 仿真了不同阵列下DOA 估计的CRB（如图中虚线所示）和所提算法在不同阵列下的DOA 估计的RMSE（如图中实线所示）。由图11 可知，不同阵列下DOA 估计的CRB 基本一致，所提算法在不同阵列下的RMSE也基本一致。并且当阵列平面孔径大小大致相当时，不同阵列下所提算法的RMSE近似相等。仿真中所提算法在两个三维立体阵列下的RMSE 要比在其他二维平面阵列下的RMSE 稍大，这是因为三维立体阵列对应的平面孔径小于二维平面阵列对应的平面孔径。由仿真结果还可以看到，在多种拓扑形式的阵列结构下，所提算法均具有高DOA 估计精度。由此可知，所提算法能够适用于任意阵列，且当阵列平面孔径大致相当时DOA估计精度基本相等。

5 结论

本文提出了一种基于变换域APSO 的宽带DOA估计算法，可用于对发射波形已知的宽带信号的俯仰角和方位角进行联合估计，并适用于任意阵列。所提算法利用已知的发射波形对阵列接收数据进行了匹配滤波处理，接着根据DML 准则构建宽带DOA 估计对应的空间谱函数，然后利用变换域APSO 搜索空间谱峰值位置，搜索结果即为宽带DOA 估计结果。仿真结果表明，所提算法具有高DOA 估计精度，在低信噪比0 dB 时RMSE 约为0.22°，高信噪比20 dB时RMSE约为0.02°。相比于牛顿迭代算法，所提算法无需DOA 预估值以及空间谱函数的导数信息，步骤简单易实现，计算复杂度低，有利于工程中的实时应用。