基于分子群改进海象优化算法的矢量水听器宽带信号DOA估计

DOA Estimation of Wideband Signals from Vector Hydrophones Based on Molecular Group Improved Whale Optimization Algorithm

摘要

针对矢量水听器阵列宽带信号波达方向（DOA）估计中存在误差大， 求解精度不高等问题， 提出一种将分子群改进的海象优化算法（MGIWaOA）与极大似然（ML）估计算法相结合的宽带信号估计方法。首先利用信号子空间变换算法构建聚焦矩阵， 将宽带信号转换为窄带信号， 然后建立基于ML的DOA估计模型， 再通过MGIWaOA进行优化估计。仿真结果表明， 该算法比海象优化算法、 正弦余弦算法、 麻雀优化算法、 差分进化算法、 遗传算法、 粒子群算法的ML估计方法具有更快的收敛速度、 更低的均方根误差和更稳健的性能。

Abstract

To address the challenges such as significant estimation errors and low resolution in the estimation of the direction of arrival（DOA） of vector hydrophone array broadband signal， this paper proposes a novel approach that combines the molecular group improved whale optimization algorithm （MGIWaOA） with maximum likelihood（ML） estimation. The method begins by utilizing the signal subspace transformation algorithm to construct a focusing matrix， transforming broadband signals into narrowband signals. Subsequently， a ML-based DOA estimation model is established， and optimization is performed using the MGI-WaOA. Simulation results demonstrate that the proposed algorithm exhibits faster convergence， lower root mean square error， and more robust performance compared to ML estimation methods based on the whale optimization algorithm， sine cosine algorithm， sparrow optimization algorithm， differential evolution algorithm， genetic algorithm， and particle swarm optimization.

Graphical abstract

关键词

宽带信号 / 矢量水听器阵列 / 信号子空间变换 / 极大似然 / 混沌对立学习 / 海象优化算法

Key words

wideband signals / vector hydrophone array / signal subspace transformation / great likelihood / chaotic antithetical learning / walrus optimization

0 引　言

波达方向（Direction of Arrival， DOA）估计是矢量水听器阵列信号处理领域的重要研究方向， 宽带信号具有频带较宽、 携带信息丰富的特性， 使得其在多个领域得到广泛的应用^［1］， 尤其在DOA估计中展现出独特的优势。传统的窄带信号DOA估计方法不能直接应用于宽带信号， 需要对宽带信号进行相应的处理^［2‐4］。宽带信号的子空间类算法可以分为非相干信号子空间算法（Incoherent Signal Subspace Method， ISM）^［5‐6］和相干信号子空间算法（Coherent Signal Subspace Method， CSM）^［7‐9］。ISM主要是将宽带接收的数据通过频域变换方法分解到不重叠的频带上， 形成若干个窄带信号； 然后， 利用窄带信号子空间处理方法对各个窄带信号进行处理， 得到窄带子带空间谱； 最后， 将所有子带空间谱进行加权组合， 得到原宽带信号的DOA估计结果。CSM是一种经典的宽带阵列信号处理算法， 这类算法引入了聚焦的思想， 是将多频率点处的数据聚焦到参考频率上， 从而形成相关矩阵， 最后通过窄带信号处理方法实现DOA估计。CSM的关键在于聚焦矩阵的选择， 不同的选择对应于不同的算法， 本文使用CSM类算法中的信号子空间变换（Signal-Subspace Transformation， SST）^［10］算法。

目前， 窄带信号DOA估计的常用算法包括多重信号分类（Multiple Signal Classification， MUSIC）^{［11‐12］}算法、 旋转不变子空间算法（Estimated Signal Parameters via Rotational Invariance Technique， ESPRIT）^［13］、 极大似然（Maximum Likelihood， ML）^{［14‐15］}估计方法等。近年来， 群智能优化算法不断涌现， 将其应用于改善传统DOA估计方法的一些限制和不足引起了学者的关注， 例如麻雀搜索算法（Sparrow Search Algorithm， SSA）^［16］、 正弦余弦算法（Sine Cosine Algorithm， SCA）^［17］、 差分进化算法（Differential Evolution， DE）^［18］、 遗传算法（Genetic Algorithm， GA）^［19］、 粒子群算法（Particle Swarm Optimization， PSO）^［20］等经典算法。Trojovsky等^［21］提出的海象优化算法（Walrus Optimization Algorithm， WaOA）是一种新颖且具有较好性能的启发式优化算法， 它结合了进化算法和模拟退火算法的优点， 具有良好的全局寻优能力、 快速的收敛速度和更好的适应性， 但是在搜索中可能陷入局部极值， 该算法的计算精度和稳定性仍有待提高。

本文针对WaOA算法存在的缺点与不足对其进行分子群改进， 提出了一种新的优化算法——分子群改进的海象优化算法（Walrus Optimization Algorithm for Multiple Group Improvement， MGIWaOA）， 并将其成功应用于极大似然算法中进行宽带信号的DOA估计。通过一系列仿真实验， 将其与WaOA、 SCA、 SSA、 DE、 GA、 PSO算法的极大似然DOA估计性能进行比较。

1 宽带信号处理

1.1 宽带信号数学模型

假设有K(M>K)个远场宽带信号s_{\mathrm{1}}(t),s_{\mathrm{2}}(t),\cdots ,s_K(t)入射到由M个矢量水听器按间距d组成的均匀线阵中， t时刻第k个入射信号的水平方位角为{\theta }_k， k=\mathrm{1,2},\cdots ,K， 则第k个入射信号的单位方向矢量{\mathit{h}}_k为

{\mathit{h}}_k=[\mathrm{1},\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_k,\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_k{]}^{\mathrm{T}}\mathrm{。}（1）

以阵列中第1个阵元为参考阵元， 则第m个阵元在t时刻的输出可表示为

x_m（t）={\displaystyle \sum _{k=\mathrm{1}}^K}{s}_k(t+{\tau }_m({\theta }_k))+n_m(t),（2）

式中： s_k(\cdot ) 为入射的宽带信号； n_m(t)为高斯白噪声； {\tau }_m({\theta }_k)为第k个信号到达第m个阵元时相对于参考阵元的时延。

对式（2）的两边做离散傅里叶变换（DFT）得到

x_m（f）={\displaystyle \sum _{k=\mathrm{1}}^K}{s}_k(f)e^{j\mathrm{2}\pi f{\tau }_m({\theta }_k)}+n_m(f)\mathrm{。}（3）

式（3）的矩阵形式为

\left[\begin{array}{l}{x}_{\mathrm{1}}(f)\\ x_{\mathrm{2}}(f)\\ \mathrm{︙}\\ x_M(f)\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{1}}({\theta }_{\mathrm{1}})}& {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{1}}({\theta }_{\mathrm{2}})}& \cdots & {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{1}}({\theta }_K)}\\ {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{2}}({\theta }_{\mathrm{1}})}& {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{2}}({\theta }_{\mathrm{2}})} & \cdots & {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_{\mathrm{2}}({\theta }_K)}\\ \mathrm{︙}& \mathrm{︙}& \mathrm{︙}& \mathrm{︙}\\ {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_M({\theta }_{\mathrm{1}})}& {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_M({\theta }_{\mathrm{2}})}& \cdots & {\mathrm{e}}^{j\mathrm{2}\pi f{\tau }_M({\theta }_K)}\end{array}\right]\left[\begin{array}{l}{s}_{\mathrm{1}}(f)\\ s_{\mathrm{2}}(f)\\ ⋮\\ s_K(f)\end{array}\right]+\left[\begin{array}{l}{n}_{\mathrm{1}}(f)\\ n_{\mathrm{2}}(f)\\ ⋮\\ n_M(f)\end{array}\right]\mathrm{。}（4）

将式（4）化为矢量形式

\mathit{X}(f)=\mathit{A}(f,\theta )\mathit{S}(f)+\mathit{N}(f),（5）

式中： \mathit{X}(f)\in {\mathbb{C}}^{\mathrm{3}M\times \mathrm{1}}为阵列接收信号矢量； \mathit{S}(f)\in {\mathbb{C}}^{K\times \mathrm{1}}为信号源矢量； \mathit{N}(f)为噪声矢量； 方向矩阵\mathit{A}(f,\theta )为

\mathit{A}(f,\theta )=[{\mathit{a}}_{\mathrm{1}}(f,{\theta }_{\mathrm{1}})\otimes h_{\mathrm{1}},{\mathit{a}}_{\mathrm{2}}(f,{\theta }_{\mathrm{2}})\otimes h_{\mathrm{2}},\cdots ,

{\mathit{a}}_K(f,{\theta }_K)\otimes h_K]\in {\mathbb{C}}^{\mathrm{3}M\times K},（6）

式中： {\mathit{a}}_k(f,{\theta }_k)为第k个信号声压阵对应的方向向量； \otimes为克罗内克积。

假设宽带信号的频率都在[f_l,f_h]之间， f_l和f_h分别为信号频率的上限和下限。将接收到的宽带信号分为不重叠的每段时间为T_{\mathrm{d}}的G段信号， 对分解后的每段信号进行J点离散傅里叶变换， 那么式（5）的宽带模型就变为

\begin{array}{c}{\mathit{X}}_g(f_j)=\mathit{A}(f_j,\theta ){\mathit{S}}_g(f_j)+{\mathit{N}}_g(f_j) ,\\ (g=\mathrm{1,2},\cdots ,G;j=\mathrm{1,2}，\cdots ,J),\end{array}（7）

式中： {\mathit{X}}_g(f_j)， {\mathit{S}}_g(f_j)， {\mathit{N}}_g(f_j)分别为第g段数据在频率f_j处的接收数据、 信号、 噪声的离散傅里叶变换； J为将带宽为B的宽带信号分割为子带的个数。式（7）在与每个子带对应的不同频率点f_{\mathrm{1}},f_{\mathrm{2}},\cdots ,f_J上均成立， 因此存在J个等式。

针对每个子带频率f_j， 计算频域采样数据的协方差矩阵

{\mathit{R}}_X(f_j)=\mathit{A}(f_j,\theta ){\mathit{R}}_S(f_j){\mathit{A}}^{\mathrm{H}}(f_j,\theta )+{\mathit{R}}_N(f_j),（8）

式中： {\mathit{R}}_X(f_j)和{\mathit{R}}_S(f_j)分别为频率f_j处满足窄带信号条件的阵元接收数据协方差矩阵和信号协方差矩阵； {(\cdot )}^{\mathrm{H}}为共轭转置； {\mathit{R}}_N(f_j)={\sigma }^{\mathrm{2}}\mathit{I}为窄带噪声协方差矩阵， \mathit{I}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathrm{1},\frac{\mathrm{1}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}},\cdots ,\mathrm{1},\frac{\mathrm{1}}{\mathrm{2}},\frac{\mathrm{1}}{\mathrm{2}}\right)。

{\mathit{R}}_X(f_j)实际是通过G次傅里叶变换（即G段， 每段包含有限数量的频点， 每个频点对应一个值）得到的， 如式（9）所示。

{\mathit{R}}_X(f_j)=\frac{\mathrm{1}}{G}{\displaystyle \sum _{g=\mathrm{1}}^G}{\mathit{X}}_g(f_j) {\mathit{X}}_g^{\mathrm{H}}(f_j)\mathrm{。}（9）

1.2 信号子空间变换算法

信号子空间变换（SST）算法利用不同频率点的信号子空间与参考频率点信号子空间之间的关系生成聚焦变换矩阵。

定理^［22］　聚焦矩阵\mathit{T}(f_j)应满足

\mathit{T}(f_j)=\mathit{U}(f_{\mathrm{0}}){\mathit{U}}^{\mathrm{H}}(f_j),（10）

式中： f_j为信号带宽内的任意频率； f_{\mathrm{0}}为参考频率； \mathit{U}(f_{\mathrm{0}})和\mathit{U}(f_j)为M\times M方阵， 前K列为Span\left\{\mathit{A}(f_{\mathrm{0}},\theta )\right\}和Span\left\{\mathit{A}(f_j,\theta )\right\}的正交基。

式（10）定义的矩阵\mathit{T}(f_j)就是SST算法的聚焦矩阵， 描述了不同频率点的导向矢量子空间与参考频率点的导向矢量子空间之间的关系。

SST算法可以概括为以下步骤：

1） 利用常规波束形成（Conventional BeamForming， CBF）算法^［23］得到波达方向的预估值 {\widehat{\theta }}_{\mathrm{1}},{\widehat{\theta }}_{\mathrm{2}},\cdots ,{\widehat{\theta }}_K， 并选定参考频率点f_{\mathrm{0}}。

2） 使用初始值构造各频率点上的阵列流形。

3） 对每个阵列流形矩阵进行奇异值分解， 即对方向矢量\mathit{A}(f_{\mathrm{0}},\theta )和\mathit{A}(f_j,\theta )进行奇异值分解。

4） 根据式（10）计算每个相应频率点处的聚焦矩阵。

5） 利用一系列聚焦矩阵对阵列接收数据进行聚焦变换， 由式（11）获得单一频率点的数据协方差矩阵\mathit{P}。

\mathit{P}=\frac{\mathrm{1}}{J}{\displaystyle \sum _{j=\mathrm{1}}^J}\mathit{T}(f_j) \mathit{P}(f_j) {\mathit{T}}^{\mathrm{H}}(f_j),（11）

式中： \mathit{P}(f_j)为阵列输出数据去噪后的协方差矩阵。

6） 利用式（16）对\mathit{P}进行处理， 即最大似然估计方法， 就可以得到宽带信号的DOA估计。

SST算法可以有效地将带宽内的信号能量聚焦到参考频率的子空间上， 1.3节将具体介绍如何利用最大似然估计算法对数据协方差矩阵\mathit{P}进行处理。

1.3 最大似然估计算法

最大似然估计是在已知存在白噪声的情况下进行的贝叶斯估计方法， 该方法利用多维数据模型和接收数据间的最佳拟合来实现DOA估计， 表现出鲁棒性高、 估计精度优越等优点， 尤其在小样本量和低信噪比（Signal Noise Ratio， SNR）的情况下。

假设信源个数已知， 则可得到L次快拍联合概率密度函数

\begin{array}{c} f_{\mathrm{D}\mathrm{M}\mathrm{L}}\left\{x_{\mathrm{1}},\cdots ,x_L\right\}=\\ {\displaystyle \prod _{l=\mathrm{1}}^L}\frac{\mathrm{1}}{\mathrm{d}\mathrm{e}\mathrm{t}\left\{\mathrm{\pi }{\sigma }^{\mathrm{2}}\mathit{I}\right\}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{1}}{{\sigma }^{\mathrm{2}}}{‖\mathit{X}-\mathit{A}\mathit{S}‖}^{\mathrm{2}}\right),\end{array}（12）

式中： \mathrm{d}\mathrm{e}\mathrm{t}\left\{·\right\}表示该矩阵的行列式； ‖·‖表示求矩阵的Frobenius范数。

对式（12）两边同时取负对数， 有

\begin{array}{l}-\mathrm{I}\mathrm{n}{f}_{\mathrm{D}\mathrm{M}\mathrm{L}}=L\mathrm{I}\mathrm{n}\mathrm{\pi }+ML\mathrm{I}\mathrm{n}{\sigma }^{\mathrm{2}}+\\ \frac{\mathrm{1}}{{\sigma }^{\mathrm{2}}}{\displaystyle \sum _{l=\mathrm{1}}^L}{‖\mathit{X}-\mathit{A}\mathit{S}‖}^{\mathrm{2}}。\end{array}（13）

要得到参数的极大似然估计， 就要求得到该对数似然函数在参数空间的极大值。式（13）中的f_{\mathrm{D}\mathrm{M}L}为关于未知参量\theta、 {\sigma }^{\mathrm{2}}和S的函数。

{\sigma }_{\mathrm{D}\mathrm{M}\mathrm{L}}^{\mathrm{2}}=\frac{\mathrm{1}}{M}\mathrm{t}\mathrm{r}\left({\mathit{P}}_A^{\perp }\widehat{\mathit{R}}\right),（14）

{\mathit{S}}_{\mathrm{D}\mathrm{M}\mathrm{L}}={\mathit{A}}^{+}\mathit{X}\mathrm{。}（15）

将式（14）和式（15）代入式（13）中可以得到关于未知参量\theta的最大似然估计，

\begin{array}{c}{\theta }_{\mathrm{D}\mathrm{M}\mathrm{L}}=\underset{\theta }{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}\left\{\mathrm{t}\mathrm{r}\left\{{\mathit{P}}_A\widehat{\mathit{R}}\right\}\right\}=\\ \underset{\theta }{ \mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\mathrm{t}\mathrm{r}\left\{{\mathit{P}}_A^{\perp }\widehat{\mathit{R}}\right\}\right\}，\end{array}（16）

式中： {\mathit{A}}^{+}为\mathit{A}的伪逆矩阵； {\mathit{P}}_A^{\perp } =\mathrm{1}-{\mathit{P}}_A为{\mathit{A}}^{\mathrm{H}}零空间上的正交投影； {\mathit{P}}_A={\mathit{A}}^{+}\mathit{A}； \widehat{\mathit{R}}为样本协方差矩阵， 即

\begin{array}{c}\widehat{\mathit{R}}=\frac{\mathrm{1}}{J}{\displaystyle \sum _{j=\mathrm{1}}^J}\mathit{T}(f_j) \mathit{X}(f_i) {\mathit{X}}^{\mathrm{H}}(f_j){\mathit{T}}^{\mathrm{H}}(f_j)=\\ \mathit{P}+\frac{\mathrm{1}}{J}{\displaystyle \sum _{j=\mathrm{1}}^J}\mathit{T}(f_j) {\mathit{R}}_N(f_i) {\mathit{T}}^{\mathrm{H}}(f_j),\end{array}（17）

式中： \mathit{P}为1.2节宽带阵列信号经过SST算法处理得到的数据协方差矩阵； {\mathit{R}}_N(f_j)为噪声协方差矩阵。

2 分子群改进的海象优化算法的宽带信号DOA估计

2.1 WaOA算法

WaOA算法是基于种群的元启发式算法， 每只海象在搜索空间中的位置决定了问题变量的候选解。算法中， 海象个体的位置在搜索空间内随机初始化， 通过计算每只海象位置的适应度值来挑选群体中最强壮的海象以引导群体中其他海象寻找食物。基于海象的自然行为， 对种群数量为N的海象个体位置进行3个不同阶段的建模。

阶段1： 喂养

根据式（18）生成海象的新位置， 如果新位置提高了目标函数的值， 则新位置取代之前的位置， 如式（19）。

q_{i,j}^{P_{\mathrm{1}}}=q_{i,j}+ran{d}_{i,j}(S{W}_j-I_{i,j}\cdot q_{i,j}),（18）

（19）

式中： Q_i^{P_{\mathrm{1}}}为第i个海象在阶段1生成的新位置； q_{i,j}^{P_{\mathrm{1}}}为在阶段1第i个海象在第j维的位置分量； F_i^{P_{\mathrm{1}}}为位置Q_i^{P_{\mathrm{1}}}的目标函数值； ran{d}_{i,j}为［0， 1］中的随机数； SW为目标函数值最佳的个体； I_{i,j}为1~2之间随机选取的数。

阶段2： 迁徙

该建模假设每只海象迁移到搜索空间另一个区域的另一个海象（随机选择）位置。首先根据式（20）生成提议的新位置， 如果这个新位置改善了目标函数的值， 则替换海象之前的位置， 如式（21）。

（20）

（21）

式中： Q_i^{P_{\mathrm{2}}}为第i个海象在阶段1生成的新位置； q_{i,j}^{P_{\mathrm{2}}}为在阶段2第i个海象在第j维的位置分量； F_i^{P_{\mathrm{2}}}为位置Q_i^{P_{\mathrm{2}}}的目标函数值； Q_k为所选择的第k只海象向第i只海象的位置迁徙， k\in \left\{\mathrm{1,2},\cdots ,N\right\}且k\ne i； q_{k,j}和F_k分别为第j维的位置向量和目标函数值。

阶段3： 逃避和对抗捕食者

逃离和对抗捕食者的策略导致海象在它们所在位置附近的位置发生变化。假设每只海象周围有一个邻域， 首先使用式（22）和式（23）在这个邻域内随机生成一个新的位置， 如果目标函数的值得到改进， 则根据式（24）替换之前的位置。

q_{i,j}^{P_{\mathrm{3}}}=q_{i,j}+(l{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t+(u{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t-rand\cdot l{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t)),（22）

Local bounds:\left\{\begin{array}{l}l{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t=\frac{l{b}_j}{t},\\ u{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t=\frac{u{b}_j}{t},\end{array}\right.（23）

（24）

式中： Q_i^{P_{\mathrm{3}}}为第i个海象在阶段3生成的新位置； q_{i,j}^{P_{\mathrm{3}}}为在阶段3第i个海象在第j维的位置分量； F_i^{P_{\mathrm{3}}}为位置Q_i^{P_{\mathrm{3}}}的目标函数值； t为当前的迭代次数； l{b}_j和u{b}_j分别为第j维变量的下界和上界； l{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t和u{b}_{\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{a}\mathrm{l},j}^t为迭代次数为t时第j维变量所允许的下界和上界， 模拟在候选解领域内的局部搜索。

2.2 分子群改进的海象优化算法

2.2.1 基于混沌对立学习策略的种群初始化

初始种群的设定直接影响着算法的收敛速度和最终结果的质量。传统元启发式算法在初始阶段通常采用随机化的方式来生成种群， 经过多次迭代过程， 探索并趋近于全局最优解。随机生成的种群个体容易出现聚集的现象， 为增强种群多样性， 可以选择优秀的初始化种群， 使得个体能够更均匀地遍历整个搜索空间， 为算法的全局搜索奠定基础。

Tent映射运算速度快， 混沌序列在［0， 1］内均匀分布， 其表达式为

{\lambda }_{t+\mathrm{1}}=\left\{\begin{array}{l}{\lambda }_t/\alpha ,{\lambda }_t\in \left.\left[\mathrm{0},\alpha \right.\right),\\ (\mathrm{1}-{\lambda }_t)/(\mathrm{1}-\alpha ),{\lambda }_t\in \left.\left[\alpha ,\mathrm{1}\right.\right),\end{array}\right.（25）

式中： {\lambda }_t为第t次迭代产生的混沌数， t=\mathrm{0,1},\cdots ,T， T为最大迭代次数； \alpha为［0， 1］中的常数， 本文选取\alpha=0.7。

对立学习策略（Opposition-Based Learning， OBL）^［24］在生成初始种群时， 主要依赖于对当前种群中个体对立点（对立解）的计算与利用。通过引入对立元素来增加算法的多样性和全局搜索能力， 从而更全面地探索潜在的解空间， 有效提高种群的质量。

假设m维空间内的一个候选解为

Q_i=(q_{i,\mathrm{1}},q_{i,\mathrm{2}},\cdots ,q_{i,m}),（26）

其下界和上界分别为l{b}_j， u{b}_j， 则对立点位置为

{\overline{Q}}_i=({\overline{q}}_{i,\mathrm{1}},{\overline{q}}_{i,\mathrm{2}},\cdots ,{\overline{q}}_{i,m}),（27）

{\overline{q}}_j=l{b}_j+u{b}_j-q_j\mathrm{。}（28）

混沌对立（Tent and Opposition-based Learning， TOBL）机制^［25］将混沌理论和对立学习策略结合起来以产生更有利于全局搜索的初始种群。TOBL机制利用目标函数的上下界之和作为中心点， 通过Tent映射均匀变化的特性， 动态地调整和优化初始种群的分布范围， 进而提高种群的质量。具体实施步骤如下： 首先， 对随机生成的n个位置应用Tent映射调整和优化； 然后， 根据OBL策略分别计算n个新位置的对立解； 最后， 将这\mathrm{2}n个位置按照适应度值进行排序， 选取适应度最佳的前n个位置作为初始化种群， 混沌对立机制的数学模型如式（29）所示。

{\overline{q}}_{i,j}=l{b}_j+u{b}_j-{\lambda }_i\otimes q_{i,j},i=\mathrm{1,2},\cdots ,N,（29）

式中： {\overline{q}}_{i,j}为第i只海象混沌对立位置的第j维分量。

2.2.2 加入分子群改进策略

为了解决WaOA算法全局搜索能力不足， 容易陷入局部最优的问题， 海象种群由最强壮的成年海象引导捕食和迁徙， 其他成年海象也具有一定的空间搜索能力， 而幼年海象行动和跟随能力较弱， 或多或少会影响种群的整体寻优能力。对此， 本文在种群位置进行迭代后， 将获得的种群适应度值从小到大排序， 接着分为3个子种群进行不同的更新操作。取排名\mathrm{2}N/\mathrm{10}+\mathrm{1}到\mathrm{8}N/\mathrm{10}分为一个子群， 命名为成年子群； 取排名\mathrm{8}N/\mathrm{10}+\mathrm{1}到N分为一个子群， 命名为幼年子群； 排名在前\mathrm{2}N/\mathrm{10}的子群命名为领导子群。

领导子群具有优异的寻优能力， 其位置保持不变； 成年子群根据式（30）进行随机差分变异更新； 幼年子群根据式（31）更新。

\begin{array}{c}{Q_i}^{\mathrm{\text{'}}}=ran{d}_{\mathrm{1}}\times (Q_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-Q_i)-ran{d}_{\mathrm{2}}\times (Q_k-Q_i) , \\ \mathrm{2}N/\mathrm{10}+\mathrm{1}\le i\le \mathrm{8}N/\mathrm{10},\end{array}（30）

{Q_i}^{\mathrm{\text{'}}}=Q_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}-Q_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}-ran{d}_{\mathrm{3}}\times (l{b}_i+ran{d}_{\mathrm{4}}\times (u{b}_i-l{b}_i)),

\mathrm{8}N/\mathrm{10}+\mathrm{1}\le i\le N,（31）

式中： Q_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}为所有个体位置矢量的平均值； l{b}_i和u{b}_i分别为动态边界的下界和上界， 通过将全局最优个体参与到幼年个体的更新过程中， 使幼年个体具备了更好的寻优能力； ran{d}_{\mathrm{1}}， ran{d}_{\mathrm{2}}， ran{d}_{\mathrm{3}}， ran{d}_{\mathrm{4}}为服从[\mathrm{0,1}]均匀分布的随机数。

种群中个体由当前最优个体引导向最优解靠拢， 随着搜索的进行， 个体均聚集在局部最优解附近， 导致种群失去多样性。式（30）通过对当前个体位置Q_i、 当前最优个体位置Q_{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}以及群体中随机选择的个体位置Q_k进行随机变异操作， 生成新的个体位置{Q_i}^{\mathrm{\text{'}}}， 有效减少了群体陷入局部最优的可能性。

2.2.3 引入自适应t分布

引入自适应t分布变异算子对每次更新完的种群进行扰动， 能够提高算法在不同阶段搜索的灵活性和鲁棒性。自适应t分布在迭代初期近似于柯西分布， 使算法具有良好的全局探索能力； 在迭代后期近似于标准正态分布， t分布算子大概率取较小值， 使得算法具备了良好的局部开发能力， 如下式

{Q_i}^{\mathrm{\text{'}}}=Q_i+Q_i\cdot t(iter)\mathrm{。}（32）

2.3 基于MGIWaOA算法的极大似然DOA估计

在SST-ML-MGIWaOA算法中， 海象的初始种群Q=\{Q_{\mathrm{1}},Q_{\mathrm{2}},\cdots ,Q_N\}代表待估计方向角\theta =\left\{{\theta }_{\mathrm{1}},{\theta }_{\mathrm{2}},\cdots ,{\theta }_N\right\}， 式（16）看作MGIWaOA算法的目标函数， 每次迭代中最强海象的位置对应待估计方向角的一个可行解， 要求解的DOA估计角度对应全局最优解。算法流程图如图 1 所示。

3 仿真实验

实验中假设远场信号分别为\theta =\left\{\mathrm{30}°,\mathrm{60}°\right\}与\theta =\left\{\mathrm{30}°,\mathrm{60}°,\mathrm{80}°\right\}入射至由8个阵元组成的均匀线性阵列中， 阵元间隔为中心频率f_c=(f_l+f_h)/\mathrm{2}对应波长的一半， 子带个数N=\mathrm{10}， 水下声速c=1 500 m/s。

在统计性能比较中， 进行了100次蒙特卡罗实验， 将均方根误差（Root Mean Square Error， RMSE）作为衡量本文仿真实验方位估计的重要指标， 计算公式为

RMSE=\frac{\mathrm{1}}{K}{\displaystyle \sum _{k=\mathrm{1}}^K}\sqrt[]{\frac{\mathrm{1}}{Z}{\displaystyle \sum _{z=\mathrm{1}}^Z}\sqrt[]{({\overline{\theta }}_{kz}-{\theta }_k{)}^{\mathrm{2}}}},（33）

式中： K和Z分别为信号数和实验次数； {\theta }_k为第k个信号的真实入射角； {\overline{\theta }}_{kz}为在第z次实验中对第k个入射信号的DOA估计值。

表 1 给出了6种优化算法的参数值。

实验 1： 仿真比较各算法的收敛性能

收敛性能可以很好地反映算法的估计效率， 在SNR分别为-10， 0， 10 dB条件下， 图 2 和图 3 分别给出了入射信号\theta =\left\{\mathrm{30}°,\mathrm{60}°\right\}和\theta =\left\{\mathrm{30}°,\mathrm{60}°,\mathrm{80}°\right\}时， MGIWaOA、 WaOA、 SCA、 SSA、 DE、 PSO 6种算法的适应度值收敛曲线。

通过观察图 2 和图 3 可以发现， 在不同信噪比条件下， MGIWaOA算法均保持着最低的适应度值和最快的收敛速度。随着迭代次数的增加， MGIWaOA算法可以准确找到适应度函数的最小值。

在图 2 中， 当SNR为0 dB和10 dB时， MGIWaOA算法在迭代次数不足70时， 适应度值已经收敛到接近于0的水平。特别是在SNR为10 dB时， MGIWaOA算法表现出最快的收敛速度， 其次是WaOA算法。从整体分析， 在SNR为-10 dB时， SCA和PSO算法陷入局部最优， DE算法的适应度值也相对较高。实验证明了MGIWaOA算法在宽带信号DOA估计方面能够更快更准确地到达收敛。

实验 2： 仿真比较各算法RMSE值和成功估计概率

图 4 为SNR在-20~20 dB范围内变化时， MGIWaOA等6种算法在DOA估计中的RMSE值以及对应的成功估计概率。值得注意的是， 无论信源数为多少， MGIWaOA算法均表现出更好、 更稳定的性能， 其次是WaOA算法。在信源数为3时， SCA、 SSA、 DE、 PSO算法变得不稳定， 而MGIWaOA算法虽然出现了一些动荡， 但它的RMSE值和估计成功概率仍然优于其他算法。实验结果说明MGIWaOA算法在不同信噪比和信源数条件下都有着较好的DOA估计性能。

实验 3： 仿真比较种群大小对RMSE值的影响

种群多样性是群智能优化算法最基本的特征， 为种群的生存和进化提供了基本前提。在解决极大似然DOA估计问题中， 种群的规模是一个关键参数， 会直接影响算法的全局寻优能力。当种群规模太小时， 种群多样性不足， 可能导致种群陷入局部最优， 影响算法的全局寻优能力； 当种群规模太大时， 可能会增加计算量， 降低算法的计算效率。

图 5 为固定SNR是10 dB， 分别在信号源个数为2和3时， 种群个数由0变化到100得到的不同算法下DOA估计的RMSE曲线。可以看到， 无论信源数为2还是3， 在种群数为20时， MGIWaOA算法相较于其他5种算法都可以保持较低的RMSE值和较高的估计精度。在信源数为3时， SCA、 SSA、 PSO算法的RMSE值受到明显影响。实验证明， 与其他算法相比， MGIWaOA算法仅需要很小的种群规模就能够实现较高的估计精度， 即MGIWaOA算法具有在保证较高估计精度的同时可以兼顾较小计算量的优点。

实验 4： 仿真比较各算法的分辨率

在本实验中， 假设入射信号分别为\theta=\left\{\mathrm{30}°\right.，（30+Δs）\left.°\right\}与\theta=\left\{\mathrm{30}°\right.，(\mathrm{30}+\mathrm{\Delta }s)°，30+2*Δs）\left.°\right\}， 角度差\mathrm{\Delta }s由0变化到10°， 固定SNR为10 dB， 得到不同角度差下DOA估计的RMSE值， RMSE值越低， 意味着估计方法的分辨率越高。

从图 6 中可以看出， 当信源个数为2时， MGIWaOA， WaOA， SSA， DE算法的RMSE值在角度差达到1°时都降低到了5以下， 并且MGIWaOA算法的RMSE值最低。随着角度差的增大， WaOA， SSA， DE算法的RMSE值不断下降， 在角度差为5°时逼近于最小值， 而MGIWaOA算法在角度差为1°时， RMSE值曲线就已经接近于最小值， 说明MGIWaOA算法在角度差为1°时已具有较高的分辨率， WaOA， SSA， DE算法次之， PSO和SCA算法的分辨率最差； 当信源个数为3°时， PSO， SSA， SCA、 DE算法受到较大的影响， MGIWaOA算法的RMSE值曲线在角度差为2°时趋于平缓且接近于最小值， 说明MGIWaOA算法在此时已具有较高的分辨率， WaOA算法次之。

4 结　论

本文将MGIWaOA算法成功应用于宽带信号DOA估计中， MGIWaOA算法的引入有效地克服了搜索容易陷入局部最优的问题， 显著提升了算法的全局搜索能力。仿真实验的结果显示， 在不同信源个数、 信噪比、 种群规模以及较小的角度差等多种情况下， MGIWaOA算法相较于WaOA、 SCA、 SSA、 DE、 GA、 PSO算法都展现出更优越的收敛效果、 更低的均方根误差以及较高的分辨率， 充分验证了基于MGIWaOA算法的宽带信号DOA估计方法的有效性。

参考文献

[1] 张进锋. 宽带信号阵列测向技术及应用研究［D］. 郑州： 战略支援部队信息工程大学， 2021.

[2] 王永良， 陈辉， 彭应宁， 等. 空间谱估计理论与算法［M］.北京： 清华大学出版社， 2004.

[3] 刘付刚. 宽带信号波达方向估计方法研究［D］. 哈尔滨： 哈尔滨工程大学， 2013.

[4] 姚鹏， 周章旺. 对宽带数字信号占用带宽测量方法的探讨［J］. 中国无线电， 2023（9）： 53-54.

[5] YAO Peng， ZHOU Zhangwang. Discussion on the measurement method of broadband digital signal occupied bandwidth［J］. China Radio， 2023（9）： 53-54. （in Chinese）

[6] AHMAD Z， SONG Y， DU Q. Wideband DOA estimation based on incoherent signal subspace method［J］. COMPEL-The International Journal for Computation and Mathematics in Electrical and Electronic Engineering， 2018， 37（3）： 1271-1289.

[7] 韩金金. 基于矢量水听器的声源定位算法研究［D］. 南京： 南京信息工程大学， 2018.

[8] TANG Y， DENG W， LI J， et al. Direction of arrival estimation of coherent wideband sources using nested array［J］. Sensors， 2023， 23（15）： 6984.

[9] HAMZA A， AMMAR M， NABIL E K M， et al. Coherent wideband signals direction finding using subspace-based methods［J］. Circuits， Systems， and Signal Processing， 2023， 42（3）： 1663-1684.

[10] 王建伟. 低截获概率信号测向技术［D］. 南京： 南京理工大学， 2008.

[11] MA F， ZHANG X. Wideband DOA estimation based on focusing signal subspace［J］. Signal， Image and Video Processing， 2019， 13（4）： 675-682.

[12] PENG S， CHEN B， YANG M. Joint sparse recovery for direction of arrival based on the generalized MUSIC criterion［J］. Digital Signal Processing， 2022， 122： 103382.

[13] 张远驰， 胡进. 一种基于MUSIC算法的宽带信号DOA估计［J］. 电声技术， 2023， 47（10）： 97-99.

[14] ZHANG Yuanchi， HU Jin. A wideband signal DOA estimation based on MUSIC algorithm［J］. Audio Engineering， 2023， 47（10）： 97-99. （in Chinese）

[15] 马娟. 雷达通信融合系统中基于数字波束形成的DOA估计研究［D］. 长沙： 湖南大学， 2016.

[16] WANG P， HUANG J， HE W， et al. Maximum likelihood DOA estimation based on improved invasive weed optimization algorithm and application of MEMS vectorhydrophone array［J］. AIMS Mathematics， 2022， 7（7）： 12342-12363.

[17] 张靖奇， 王鹏， 梅银珍. 基于野马优化的宽带信号波达方向估计［J］. 测试技术学报， 2023， 37（3）： 208-213.

[18] ZHANG Jingqi， WANG Peng， MEI Yinzhen. Wideband signal DOA estimation based on wild horse optimization［J］. Journal of Test and Measurement Technology， 2023， 37（3）： 208-213. （in Chinese）

[19] BAIRATHI D， GOPALANI D. Salp swarm algorithm （SSA） for training feed-forward neural networks［C］//Bansal J， Das K， Nagar A， et al. Soft Computing for Problem Solving. Singapore： Springer， 2019： 521-534.

[20] MIRJALILI S. SCA： a Sine Cosine Algorithm for solving optimization problems［J］. Knowledge-Based Systems， 2016， 96： 120-133.

[21] CIVICIOGLU P， BESDOK E， GUNEN M A， et al. Weighted differential evolution algorithm for numerical function optimization： a comparative study with cuckoo search， artificial bee colony， adaptive differential evolution， and backtracking search optimization algorithms［J］. Neural Computing and Applications， 2020， 32（8）： 3923-3937.

[22] EL-SHORBAGY M A， EL-REFAEY A M. Hybridization of grasshopper optimization algorithm with genetic algorithm for solving system of non-linear equations［J］. IEEE Access， 2020， 8： 220944-220961.

[23] GAD A G. Particle swarm optimization algorithm and its applications： a systematic review［J］. Archives of Computational Methods in Engineering， 2022， 29（5）： 2531-2561.

[24] TROJOVSKÝ P， DEHGHANI M. A new bio-inspired metaheuristic algorithm for solving optimization problems based on walruses behavior［J］. Scientific Reports， 2023， 13（1）： 8775.

[25] 翟丽杰， 李志臣. 基于改进SST算法的宽带信号DOA估计［J］. 科学技术与工程， 2012， 12（23）： 5771-5774.

[26] ZHAI Lijie， LI Zhichen. Based on improved SST algorithm for estimation of angles of arrival of multiple wide-band sources［J］. Science Technology and Engineering， 2012， 12（23）： 5771-5774. （in Chinese）

[27] NIE W， ZHANG X， XU J， et al. Adaptive direction-of-arrival estimation using deep neural network in marine acoustic environment［J］. IEEE Sensors Journal， 2023， 23（13）： 15093-15105.

[28] TIZHOOSH H R. Opposition-based learning： a new scheme for machine intelligence［C］//International Conference on Computational Intelligence for Modelling， Control and Automation and International Conference on Intelligent Agents， Web Technologies and Internet Commerce （CIMCA-IAWTIC'06）， 2005： 695-701.

[29] 马驰， 曾国辉， 黄勃， 等. 融合混沌对立和分组学习的海洋捕食者算法［J］. 计算机工程与应用， 2022， 58（22）： 271-283.

[30] MA Chi， ZENG Guohui， HUANG Bo， et al. Marine predator algorithm based on chaotic opposition learning and group learning［J］. Computer Engineering and Applications， 2022， 58（22）： 271-283. （in Chinese）

基金资助

国家自然科学基金资助项目(61774137)

山西省基础研究计划资助项目(202103021224195)

山西省基础研究计划资助项目(202103021224212)

山西省基础研究计划资助项目(202103021223189)

山西省基础研究计划资助项目(20210302123019)

山西省回国留学人员科研资助项目(2020-104)

山西省回国留学人员科研资助项目(2021-108)

山西省回国留学人员科研资助项目(2022-149)