1.   引言

跳频(Frequency Hopping, FH)通信因其保密性好和抗截获、抗干扰能力强等特性在军事通信中得到广泛应用^[1–3]。FH信号的波达方向是辅助FH信号跟踪及网台分选的关键因素，因此对DOA进行精确估计具有重要意义。

可查文献中关于FH信号的DOA估计问题主要有基于空时频分析和多重信号分类(MUltiple SIgnal Classification, MUSIC)两种方法。文献[4]对信道化接收得到的数据分别进行时频分析从而得到全景时频图，提取有效跳之后建立其空时频矩阵然后完成每一跳的DOA估计，适用于欠定条件，但其估计性能依赖于时频分析的性能，且在多信号环境下提取有效跳存在困难；文献[5,6]将MUSIC算法与空时频分析相结合，先得到空时频矩阵，然后利用root-MUSIC算法对其处理以得到DOA的估计，但其性能在较少快拍数的情况下显著下降^[7]；文献[8]将极化信息引入DOA估计，利用空间极化时频分布矩阵中蕴含的极化-空域信息同时完成DOA和极化参数的估计，但是该算法对噪声敏感，且计算量大；文献[9]通过空时频分析得到时频域协方差矩阵，将共轭子空间的思想引入到MUSIC算法中，提出一种基于MUSIC对称压缩谱(MUSIC Symmetrical Compressed Spectrum, MSCS)的算法，降低了MUSIC算法的复杂度，但只能用于不相关的源信号。文献[4–6,8,9]均未考虑FH信号的频域稀疏性和DOA的空域稀疏性。文献[10,11]将稀疏重构方法用于DOA的估计，前者构造L1范数最小化模型来求解稀疏重构问题，并通过奇异值分解(Singular Value Decomposition, SVD)降低了运算量，后者在空间频率域建立过完备字典，先估计出信号的瞬时频率，然后得到DOA的估计，但是都没有考虑构建的字典与实际中连续的方位角之间的不匹配问题。稀疏贝叶斯学习理论最早被用于机器学习领域^[12]。文献[13]将SBL用于多测量向量(Multiple Measurement Vector, MMV)模型下的稀疏信号重构，获得了更高的估计精度。已有学者将SBL用于DOA估计^[14–17]。文献[14,15]利用SBL恢复联合行稀疏矩阵，相对于稀疏重构方法在保证估计性能的同时降低了运算的复杂度，文献[16]提出一种基于互质阵列协方差稀疏感知的DOA估计方法，利用线性变换从协方差向量中去除噪声方差以简化SBL，以上两种方法都只适用于真实DOA恰好落在空间角度网格上的情况。文献[17]把空间角度网格与实际DOA之间的偏差建模为超参数并通过SBL估计出来，同时利用奇异值分解减小了计算量，但对信号作了拉普拉斯先验假设，且只适用于连续信号。

综上所述，本文针对多FH信号DOA估计问题，利用DOA的空域稀疏性在空间角度域构造离散网格，建立包含网格点与真实DOA之间偏移量的阵列接收模型，在稀疏贝叶斯学习的基础上提出一种OGSBL(Off Grid-SBL)算法以得到信号与偏移量的联合估计。利用SBL算法得到包含超参数的输入信号后验分布，通过迭代得到收敛的超参数，从而完成DOA的估计。同时，为在整个空间角度域得到相同的估计精度，利用线性插值的方法建立偏移量的模型。仿真结果表明该算法在较少快拍数的情况下相比以往方法精确度和估计速度均有提升，在信噪比为0 dB时仍能取得较好的估计效果。

2.   数学模型

2.1   跳频信号阵列接收模型

如图1所示，假设$K$个远场FH信号${{s}}(t) =$${[{s_1}(t), {s_2}(t), ·\!·\!· , {s_K}(t)]^{\rm{T}}}$入射到一$M$元均匀线阵上，阵元间距为$d$，入射方向为${{θ }} = {[{\theta _1}, {\theta _2}, ·\!·\!· , {\theta _K}]^{\rm{T}}}$，第$k$个跳频信号${s_k}(t)$的幅度为${a_k}$，跳周期为${T_k}$，在观测时间内共有$B$跳，其表达式为${s_k}(t) ={a_k}$$\cdot\displaystyle\sum\nolimits_{b = 1}^{\rm{B}}\exp{({\rm j}2{{π}} {f_{k, b}}{t'})} {\rm rect}({{{t'}}}/{{{T_k}}})$，式中${t'} = t \!-\!(b \!-\! 1)$$\cdot{T_k}$, ${f_{k, b}}$为${s_k}(t)$第$b$跳的载频^[18]，则阵列接收信号${{y}}(t)$可表示为

图 1  均匀线阵接收模型

下载: 全尺寸图片 幻灯片

其中，${{y}}(t) = {[{y_1}(t), {y_2}(t), ·\!·\!· , {y_M}(t)]^{\rm{T}}}$, ${{v}}(t)=[{v_1}(t),$${v_2}(t), ·\!·\!· ,{v_M}(t){]^{\rm{T}}}$, ${y_m}(t)$和${v_m}(t)$, $m = 1,$$2, ·\!·\!· , M$分别代表第$m$个阵元上的输出和噪声。用相移来表示信号入射到不同阵元上的时间差，则有阵列流型矩阵${{A}}({{θ }}) = [{{a}}({\theta _1}), {{a}}({\theta _2}), ·\!·\!· , {{a}}({\theta _K})],{{a}}({\theta _k}) = {[{e^{{\rm j}2{{π}} {f_k}{\tau _{k, 1}}}}},$${ {e^{{\rm j}2{{π}} {f_k}{\tau _{k, 2}}}}, ·\!·\!· , {e^{{\rm j}2{{π}} {f_k}{\tau _{k, M}}}}]^{\rm{T}}}$为第$k$个FH信号的导向矢量，其中${\tau _{k, m}} = (m - 1)d\cos ({\theta _k})/c$, $c$代表真空中的光速。

FH信号是宽带信号，但是在同一频率驻留时间内可当作窄带信号进行处理。在已得到FH信号跳变时刻的基础上，将FH信号以跳变时刻为界进行分割，得到多个分段信号，然后对其分别进行处理。将入射角空间$[0, {{π}} ]$进行离散划分以得到一个入射方向的有限网格${{ϑ}} = {[{\vartheta _1}, {\vartheta _2}, ·\!·\!· , {\vartheta _N}]^{\rm{T}}}$, $N$代表网格的个数且$N \gg M > K$。当$K$个FH信号的入射角恰好都落在网格点上，即$\theta \in {{ϑ}}$时，有式(2)所示的单测量向量(Single Measurement Vector, SMV)模型：

其中，${{A}}({{ϑ}} ) = [{{a}}({\vartheta _1}), {{a}}({\vartheta _2}), ·\!·\!· , {{a}}({\vartheta _N})]$, $\overline {{s}} (t) =$$\left[ {{{\bar s}_1}(t),} \right.{\left. {{{\bar s}_2}(t), ·\!·\!· ,{{\bar s}_N}(t)} \right]^{\text{T}}}$是将${{s}}(t)$按行进行扩展后得到的新信号矩阵，$\overline {{s}} (t)$的非零行对应的网格点就是信号实际的入射方向，即

当接收到$L$个快拍的数据时，式(2)所表示的SMV模型变为如式(4)MMV模型：

其中，${{Y}} = [{{y}}({t_1}), {{y}}({t_2}), ·\!·\!· , {{y}}({t_L})]$为阵列输出矩阵，$\overline {{S}} = [\overline {{s}} ({t_1}), \overline {{s}} ({t_2}), ·\!·\!· , \overline {{s}} ({t_L})]$为扩展后的信号矩阵，${{V}} = [{{v}}({t_1}), {{v}}({t_2}), ·\!·\!· , {{v}}({t_L})]$为阵列的噪声矩阵。已被证明^[19,20]，与SMV模型相比，MMV模型下重构的成功率更高，参数的估计性能也更好。

2.2   离散网格偏移模型

无论网格${{ϑ}}$划分得多密，连续的DOA与离散的网格之间总是存在不匹配，这严重限制了DOA的估计精度，同时随着网格密度的增大，计算量也会增加。为了避免网格不匹配问题带来的性能恶化，引入一组描述偏移量的参数${{ρ }} = {[{\rho _1}, {\rho _2}, ·\!·\!· , {\rho _{N - 1}}]^{\rm{T}}}$来对导向矢量${a}({\vartheta _n})$进行修正：

其中，${\underline \vartheta _n}$和${\overline \vartheta _n}$分别表示真实DOA${\vartheta _n}$左侧和右侧相邻的网格点，于是流型矩阵变为

${{A}}(u:v)$表示取矩阵${{A}}$的第$u$列到第$v$列。为计算方便，将0补在${{ρ }}$的最后使之成为一个$N \times 1$的向量，定义${Δ} = {\rm diag}({{ρ }})$，${{{A}}_b} = {{A}}(2:N, {\bf{0}})$，则新的流型矩阵可表示为

于是式(4)所表示的MMV模型变为

在${{ρ }} = {\bf{0}}$时式(8)与式(4)相同。我们的目标就是通过阵列输出${{Y}} = [{{y}}({t_1}), {{y}}({t_2}), ·\!·\!· , {{y}}({t_L})]$和一组对应关系${{ϑ}} \to {{A}}({{ϑ}} )$对${{θ }} = {[{\theta _1}, {\theta _2}, ·\!·\!· , {\theta _K}]^{\rm{T}}}$和${{ρ }} = {[{\rho _1}, {\rho _2}, ·\!·\!· , }$${{\rho _{N - 1}}, 0]^{\rm{T}}}$进行估计。

3.   算法原理

3.1   OG-SBL算法

假设信号矩阵$\overline {{S}}$的每一列之间相互独立，且服从均值为0的高斯分布：

其中，${{γ }} = {[{\gamma _1}, {\gamma _2}, ·\!·\!· , {\gamma _N}]^{\rm{T}}}$, ${Γ} = {\rm diag}({{γ }})$, ${\gamma _i}$是控制$\overline {{S}}$行稀疏性的非负超参数，当${\gamma _i} = 0$时$\overline {{S}}$的第$i$行为零，${Γ}$为矩阵$\overline {{S}}$的第$l$列的协方差矩阵。于是信号矩阵$\overline {{S}}$的概率密度函数为

假设每个阵元上的噪声是相互独立的，即噪声矩阵${{V}}$的每一行之间相互独立且服从均值为$0$方差为$\lambda$的复高斯分布，有$p({{{v}}_n}({t_l}))\sim\cal {N}(0, \lambda )$，故${{Y}}$的似然函数为$p({{Y}}|\overline {{S}} ;\lambda , {{ρ }})\sim{{\cal {N}}_{{{Y}}|\overline {{S}} }}(\overline {{A}} \overline {{S}} , \lambda {{I}})$, $\overline {{S}}$的先验概率由$p(\overline {{S}} ;{Γ})\sim{{\cal {N}}_{\overline {{S}} }}({\bf{0}}, {Γ} )$给出，然后利用贝叶斯准则可得到$\overline {{S}}$的后验概率密度函数：

其均值和协方差矩阵分别为

式中${Σ}\nolimits_{y}^{{-1}}= \lambda {{I}} + \overline {{A}}{Γ} {\overline {{A}} ^{\rm{H}}}$。利用EM算法得到未知超参数的估计后，就能得到$\overline {{S}}$的最大后验估计，当${\gamma _i} = 0$时令$\overline {{S}}$的相应行等于零向量，然后根据$\overline {{S}}$中非零行对应的网格点及偏移量${{ρ }}$即可得到DOA的估计。

3.2   超参数的估计

为得到超参数集$\varTheta = \{ \lambda , {{ρ }}, \varGamma \}$的估计，利用EM算法使观测矩阵${{Y}}$关于超参数集$\varTheta$的概率密度函数$p({{Y}}, \overline {{S}} ;\varTheta )$最大，等价于使$- \lg p({{Y}}, \overline {{S}} ;\varTheta )$最小。在迭代过程中将$\overline {{S}}$视为隐藏变量，则有：

${\varTheta ^{({\rm old})}}$代表上一次迭代后得到的超参数。式(14)中的第1项与${Γ}$无关，故Q函数$Q({Γ} )$可简化为

${\mu _{\overline {{S}} }}$和${{Σ} _{\overline {{S}} }}$为根据${\varTheta ^{({\rm old})}}$得到的均值和协方差矩阵。令式(15)关于${\gamma _i}(i = 1, 2, ·\!·\!· , N)$的导数为零，得到超参数${\gamma _i}$的学习规则：

其中，${({{Σ} _{\overline {{S}} }})_{(i, i)}}$表示位于${{Σ} _{\overline {{S}} }}$第$i$行第$i$列的元素，${{{({{{μ}} _{\overline {{S}} }})}_{i \cdot }}}$表示$\ {{{μ }}_{\overline {{S}} }}$的第$i$行。

同样地，式(14)中的第2项与$\lambda$和${{ρ }}$无关，故Q函数$Q(\lambda , {{ρ }})$可简化为

令式(17)关于$\lambda$的导数为零，得到超参数$\lambda$的学习规则：

将$\overline {{A}} = {{A}} + ({{{A}}_b} - {{A}}){Δ}$代入式(17)，得到超参数${{ρ }}$的学习规则：

${{A}} \odot {{B}}$表示两个矩阵的Hadamard积，此处利用了文献[17]中的推导。至此经过1次迭代后的超参数集$\Theta$得到了更新，在算法收敛后，根据${{γ }} = {[{\gamma _1}, {\gamma _2}, ·\!·\!· , {\gamma _N}]^{\rm{T}}}$得到$\overline {{S}}$非零行的位置，根据其对应的网格点及偏移量即可完成DOA的估计：

$\delta$为网格${{ϑ}}$的划分间隔。

4.   算法流程与复杂度分析

4.1   算法流程总结

步骤 1  输入阵列输出矩阵${{Y}}$；

步骤 2  由划分好的离散网格${{ϑ}}$得到阵列流型矩阵$\overline {{A}}$；

步骤 3  初始化超参数集$\varTheta$，由式(12)，式(13)计算${{{μ }}_{\overline {{S}} }}$和${{Σ} _{\overline {{S}} }}$；

步骤 4  利用EM算法对超参数进行迭代，直到收敛：

(1) 依据式(16)更新${{γ }}$；

(2) 依据式(18)更新$\lambda$；

(3) 依据式(19)更新${{ρ }}$；

(4) 由超参数$\Theta$更新${{{μ }}_{\overline {{S}} }}$和${{Σ} _{\overline {{S}} }}$；

步骤 5  由式(20)得到DOA的估计。

4.2   算法复杂度分析

以复数乘法的次数作为运算复杂度的标准，忽略标量乘法等运算计算量，为便于分析，假设$N \gg L > M$。超参数每次迭代之前都需要先计算出更新后的${{{μ }}_{\overline {{S}} }}$和${{Σ} _{\overline {{S}} }}$，计算${{{μ }}_{\overline {{S}} }}$和${{Σ} _y}$都需要$O(N{M^2})$次复数乘法，同时对${{Σ} _y}$求逆需要$O({M^3})$次复数乘法，然后计算${{Σ} _{\overline {{S}} }}$的复杂度为$M{N^2}$。$\lambda$和${{ρ }}$的更新复杂度分别为$O(M{N^2})$和$O({N^3})$，所以总的算法复杂度约为$O(N{M^2}) + O({M^3}) + M{N^2} + \zeta [O(M{N^2}) + O({N^3})]$, $\zeta$表示迭代的次数。可见减少阵元数$M$和离散网格的网格点数$N$均可降低复杂度，但$M$与$N$对于算法的精度是有影响的，将在下一节通过仿真实验具体说明。

5.   仿真实验与分析

均匀线阵结构如图1所示，阵元间距设为$1.5\;{\rm{m}}$, 4个远场FH信号记为FH1～FH4，跳周期均为${\rm{2\; {μ} s}}$，采样率为$100\;{\rm{MHz}}$, DOA分别为${\theta _1} = - {37.5^ \circ }$, ${\theta _2} \!= - {48.3^ \circ }$, ${\theta _3} \!=\! {24.1^ \circ }$和${\theta _4} \!=\! {39.7^ \circ }$，在观测时间内的跳频图案分别为$[35.4, 47.5, 39.6]\;{\rm{MHz}}$, $\left[ {{\text{48}}{\text{.5, }}}{\text{42.8,}} \right.$35.6] MHz, [43.8, 34.8, 49.6] MHz和$\left[ {39.7, 30.6, } \right.$$\left. {42.8} \right]\;{\rm{MHz}}$，每次取20个快拍组成的分段信号进行处理。离散网格的划分范围从$- {90^ \circ }$~${90^ \circ }$，信噪比从–10 dB以2 dB为步长递增至20 dB，本文实验在每个信噪比下均进行100次蒙特卡罗实验，以均方根误差(Root Mean Square Error, RMSE)作为算法性能的评价标准，本文中RMSE定义为

其中，T为蒙特卡罗实验的次数。

5.1   实验1

为验证阵元数$M$对算法精度和复杂度的影响，假设空间内存在FH1和FH2两个跳频信号，设置$M = 5, 7, 9$，离散网格的间隔$\delta = {2^ \circ }$, RMSE与算法运行时间随信噪比的变化情况如图2所示(RMSE采用对数坐标)。

图 2  不同阵元数下算法性能与运行时间的比较

下载: 全尺寸图片 幻灯片

可以看出，随着阵元数$M$的增加，算法在同一信噪比下的RMSE有所降低，但其运行时间也会增长。随着阵元数的增加，阵列流型矩阵$\overline {{A}}$和接收信号矩阵${{Y}}$的行数增多，因此DOA估计精度会提高但算法时间复杂度也会变高。考虑到在信噪比大于0 dB时$M = 7$与$M = 9$下RMSE的差别不大，而运行时间有较大差距，综合考虑在后面实验中取$M = 7$。

5.2   实验2

为验证网格点数对算法精度和复杂度的影响，假设空间内存在FH1和FH2两个跳频信号，设置离散网格的间隔$\delta = {1^ \circ }, {2^ \circ }, {3^ \circ }$, RMSE与算法运行时间随信噪比的变化情况如图3所示(RMSE采用对数坐标)。

图 3  不同网格间隔下算法性能与运行时间的比较

下载: 全尺寸图片 幻灯片

可以看出，网格的划分间隔越大，算法在同一信噪比下的RMSE也就越大，同时运行时间越短。这是因为$\delta$增大会使网格点数变少，阵列流型矩阵各列之间的相关性就会降低，网格变得更疏，在超参数迭代时的复杂度大大降低，所以运行时间更短，而DOA估计精度会相应降低。综合考虑估计精度和时间复杂性，在后面实验中取$\delta = {2^ \circ }$。

5.3   实验3

空间内4个FH信号同时存在，为比较本文算法与文献[9]MSCS算法及文献[11]稀疏重构算法的空间谱差异，选取快拍数$L$分别等于20和80，各算法在信噪比为10 dB下的空间谱如图4。

图 4  各算法的空间谱比较

下载: 全尺寸图片 幻灯片

图4(a)中放大的是第2个谱峰，图4(b)中放大的是第3个谱峰。由图4可见相对于其他两种算法，OG-SBL算法在$L = 20$时空间谱的谱峰更加锋利，且更为接近真实的DOA，在$L = 80$时OG-SBL算法与稀疏重构算法的性能相当，都优于MSCS算法。故本文算法在只观测到较少快拍数据的情况下性能更好。

5.4   实验4

实验条件同实验3，为比较本文算法与文献[11]稀疏重构算法及文献[17]离格稀疏贝叶斯学习(Off Grid Sparse Bayesian Inference, OGSBI)算法的性能差异，选取快拍数$L$分别等于20和80, 3种算法的性能曲线如图5所示。

图 5  不同快拍数下算法性能的比较

下载: 全尺寸图片 幻灯片

可见随着快拍数的增加，3种算法的性能均有所提升，在$L = 20$时OG-SBL算法相较于其他两种算法具有更低的RMSE，在$L = 80$时，信噪比低于4 dB时OG-SBL算法更优，信噪比高于4 dB时稀疏重构算法更优，这是因为快拍数的增加意味着时域采样点数的增加，采样点数越多稀疏重构的结果也就越精确，稀疏重构算法的这种优点在高信噪比时得到体现。可见本文算法在较少采样点的情况下精度更高。

5.5   实验5

实验条件同实验3，为比较本文算法与文献[11]稀疏重构算法及文献[17]OGSBI算法复杂度的差异，设置在快拍数分别为20和80下每种算法进行100次蒙特卡罗实验，3种算法在信噪比为10 dB下所用时间如表1所示。

表 1  不同快拍数下算法运行时间的比较(s)

快拍数	20	80
本文算法所用时间	0.3804	0.5915
稀疏重构算法所用时间	0.4066	0.3935
OGSBI算法所用时间	0.6454	0.8428

下载: 导出CSV 

| 显示表格

可见本文算法在不同快拍数下所用时间均少于OGSBI算法，在$L = 20$时本文算法相比稀疏重构算法复杂度更低，而在$L = 80$时稀疏重构算法复杂度更低。

6.   结束语

跳频信号源的DOA是辅助网台分选、跟踪和引导干扰的重要信息，实际电子对抗等场合对参数估计实时性提出了更高的要求。本文基于SBL提出一种适用于多跳频信号的DOA估计方法，考虑到真实DOA与离散网格之间的偏差，适用于较少快拍的情况，能更好地满足实时性的要求。理论分析和仿真实验表明本文方法能够解决多跳频信号的DOA精确估计问题，如何将本文方法用于相关信号源以及提升在低信噪比下的性能有待进一步研究。