文献阅读：Bayesian Nonparametric Hidden Markov Model for Agile Radar Pulse Sequences Streaming Analysis

用于敏捷雷达脉冲序列流分析的贝叶斯非参数隐马尔可夫模型

摘要：多功能雷达（MFR）是一类复杂传感器，具备实现复杂脉冲间敏捷调制和动态工作模式调度的能力。多功能雷达的发展对现代电子侦察系统或雷达告警接收机识别与推断其工作模式构成了巨大挑战。为解决这一问题，本文提出一种用于多功能雷达工作模式参数估计与变点检测的在线处理框架。首先，本文设计了一种具有特定先验分布的全共轭贝叶斯非参数隐马尔可夫模型（agile BNP-HMM），以表征多功能雷达的脉冲敏捷特性。随后，所提框架由两部分核心内容构成：第一部分是agile BNP-HMM模型，该模型可自动推断隐马尔可夫模型（HMM）的隐状态数量及对应隐状态的发射分布，并推导了估计性能的误差下界，实验表明所提算法相较于基准方法更接近该下界；第二部分结合流贝叶斯更新以提升计算效率，并基于加权序贯概率比检验（weighted SPRT）设计了在线工作模式变点检测框架。通过在多种仿真雷达信号数据和真实基准数据集上的实验验证，所提框架相较于基准方法始终具有更高的有效性和鲁棒性。源代码可在https://github.com/JiadiBao/Agile-BNP-HMM获取。

关键词：变点检测；概率图模型；脉冲间调制；多功能雷达；非参数贝叶斯模型；变分推断

1 引言

多功能雷达（MFR）可在雷达时间轴上调度多种同时工作模式以完成不同任务[1-3]，例如监视、目标跟踪和跟踪维持。随着敏捷性的提升，多功能雷达能够根据特定任务目标调整并优化其控制参数[4]，如脉冲重复间隔（PRI）、射频（RF）和脉冲宽度（PW）。多功能雷达日益增强的灵活性与适应性，对现代电子支援（ES）系统[5-7]构成了巨大挑战。为保护目标免受雷达威胁，电子支援系统需在可能缺乏先验信息的情况下识别动态复杂的工作模式，并尽快检测出多功能雷达工作模式的转换，从而及时调度有效的战术措施和干扰对抗手段。

传统上，多功能雷达工作模式识别基于有监督方法，研究者已设计了多种多功能雷达模型与识别方法[8-11]。然而，由于现代多功能雷达具备敏捷性和软件定义能力[12-13]，动态电磁环境中不断出现新的雷达及雷达信号，导致传统有监督识别方法所需的先验信息难以获取或失效。为降低对先验信息的依赖，研究者提出基于模型的时间序列聚类方法，用于多功能雷达脉冲序列的无监督识别与脉冲间调制参数估计[14]。但文献[14]中的方法在真实电磁环境下，面对杂波脉冲和缺失脉冲造成的非理想条件时，仍需进一步优化算法设计；此外，针对具有动态敏捷脉冲间调制的工作模式，该方法采用基于搜索的策略确定脉冲重复间隔（PRI）电平的真实数量[14]，这种搜索操作计算开销大，不适用于在线处理场景。因此，亟需研究一种能够在非理想条件下对不同雷达工作模式进行鲁棒识别与参数估计，并以在线方式检测连续工作模式间变点的方法。

快速变点检测[15]是解决上述多功能雷达相关问题的有效途径，其首要步骤是对不同工作模式的雷达脉冲序列进行有效建模。通常，雷达脉冲序列可通过参数化模型建模，该模型结合多雷达控制参数（如PRI、RF、PW）的脉冲间与脉冲内调制。例如，文献[14]采用四种参数化模型表征具有不同脉冲间调制的脉冲序列；文献[16-17]利用高斯混合模型（GMM）刻画雷达脉冲序列；文献[1-3]则采用隐马尔可夫模型（HMM）对著名的“Mercury”多功能雷达的脉冲序列进行建模①。隐马尔可夫模型能够将不同雷达脉冲序列的脉冲间调制类型与调制参数映射到统一的状态空间模型中：每个隐状态对应特定的脉冲参数值，且可通过引入额外隐状态来建模非理想条件的影响。在非协作环境下，隐马尔可夫模型推断存在一个问题——隐状态数量$K$未知且需预先设定；此外，由于雷达脉冲序列隐马尔可夫模型中的隐状态具有明确物理意义，$K$的不当设置必然会降低模型性能并削弱可解释性[18]。

为获得更优的表征模型，假设隐状态数量为无穷大的贝叶斯非参数隐马尔可夫模型（BNP-HMM）可能是一种潜在选择。在贝叶斯非参数隐马尔可夫模型中，每个有效状态对应后验概率较大的隐状态。文献[19-20]首次尝试将无穷隐状态假设应用于隐马尔可夫模型，此后，基于该假设的研究广泛开展，例如在说话人分割[21]、雷达高分辨距离像（HRRP）识别[22-23]和人体运动预测[24]等领域。贝叶斯非参数理论为多功能雷达脉冲序列建模及少先验信息下的参数估计提供了理论基础。然而，据作者所知，目前尚未有将贝叶斯非参数理论应用于截获雷达信号分析的研究；同时，上述研究均基于隐状态转移的非敏捷或“粘性”假设，不适用于估计存在非理想条件干扰的敏捷动态多功能雷达工作模式隐状态转移。因此，对于多功能雷达应用场景，现有贝叶斯非参数隐马尔可夫模型方法需解决两个关键问题：1）需针对多功能雷达脉冲序列的动态敏捷特性与非理想条件进行特定设计；2）需设计具备高效模型学习与推断能力的在线处理框架。

多功能雷达工作模式快速变点检测的第二步是设计检测策略。由于以下两方面原因，隐马尔可夫模型（实际为状态空间模型）的快速变点检测具有挑战性：首先，若已知变前分布但未知变后分布，由于无法预先知晓可能发生的变化类型（包括结构变化和趋势变化），估计甚至建模变后分布往往难以实现[25-27]；其次，变点检测检验需计算KL散度（Kullback-Leibler散度），但在状态空间模型中计算KL散度通常十分困难[28-29]。对于隐马尔可夫模型，常用马尔可夫链蒙特卡洛（MCMC）仿真近似KL散度，并可推导其渐近性质[26]；文献[28]将隐马尔可夫模型视为一个大型马尔可夫链，基于极小极大准则构建递归累积和（CUSUM）；文献[29]还推导了对数似然，并提出一种非蒙特卡洛方法计算两状态隐马尔可夫模型的KL散度。然而，一方面，适用于任意有限状态隐马尔可夫模型KL散度的有效数值算法仍有待研究；另一方面，用于变前与变后分布建模的贝叶斯非参数隐马尔可夫模型具有非遍历性且隐状态空间为无穷基数，这导致经典隐马尔可夫模型的通用结论不再适用。

鉴于上述问题，本文提出一种新框架，用于多功能雷达工作模式的参数估计（PE任务）与变点检测（CPD任务），该框架称为Agile BNP-HMM-CUSUM，由多功能雷达脉冲序列建模、工作模式参数估计和变点检测三个连续步骤构成。首先，为实现更优的脉冲序列建模，提出贝叶斯非参数隐马尔可夫模型的一种变体——agile BNP-HMM，以更好地控制所推断的隐马尔可夫模型隐状态，这种控制对所研究的多功能雷达问题至关重要。其次，为agile BNP-HMM设计高效的变分推断方法：该方法可通过敏捷狄利克雷过程（agile DP）自动估计贝叶斯非参数隐马尔可夫模型的隐状态数量；敏捷狄利克雷过程是一种“测度上的测度”，由基础分布和正尺度参数表征；为实现敏捷狄利克雷过程，基于传统断棒构造[30]提出相应的敏捷断棒构造方法。此外，将agile BNP-HMM模型与流贝叶斯更新[31]处理相结合，流更新可避免参数估计任务陷入局部最优并提升计算效率。最后，为克服隐马尔可夫模型变点检测因实时处理需求导致的计算开销大问题并保持高性能，借鉴文献[29]的思路，将初始分布和转移矩阵视为冗余参数，从而将隐马尔可夫模型的变点检测问题转化为多元高斯序列的变点检测问题，并针对多功能雷达应用场景提出加权序贯概率比检验（weighted SPRT）②。

本文的主要贡献可总结如下：
1）提出agile BNP-HMM，用于对具有敏捷动态特性的多功能雷达工作模式进行建模；该模型为全共轭模型，支持高效的变分推断。
2）提出一种改进的断棒构造方法；与Sethuraman[30]提出的传统断棒构造方法相比，所提方法在面对敏捷脉冲间调制类型时可实现更准确的估计，并推导了非理想条件下估计性能的误差下界。
3）基于加权序贯概率比检验（weighted SPRT）[34]设计一种渐近最优的变点检测策略；该策略无需滑动窗口类方法所需的预设窗口大小信息，避免了窗口调整难题，且参数调优更简便。
4）将所提Agile BNP-HMM-CUSUM与流贝叶斯更新技术[31]相结合，实现在线处理，降低计算冗余。

本文其余部分结构如下：第2章阐述参数估计与变点检测任务的数学建模；第3章介绍用于参数估计与变点检测任务的agile BNP-HMM-CUSUM框架；第4章给出仿真设计、结果及分析；第5章介绍在真实基准数据集上的实验；第6章总结全文。本文主要符号如表1所示。

表1 本文所用主要符号

符号	说明
$P_t$	索引为$t$的数据批次
$S$	隐马尔可夫模型的底层状态序列
$\pi$、$A$	隐马尔可夫模型的初始分布与转移分布
${\Theta }_t$	时刻$t$的高斯分布集合
${\rm Y}$	贝叶斯非参数隐马尔可夫模型的超参数集合
$Z_t$	索引为$t$的数据批次的潜变量
${\varphi }_k$	第$k$个隐状态的参数化模型
$T$	数据批次长度
$L$	截断水平
$K$	隐状态数量
$N$	检测到的变点（停止时间）
$\mathcal{N}$	高斯分布
$\nu$	真实变点
${\theta }_0$、${\theta }_1$	变前、变后分布的参数
${\mathbb{E}}_q$	关于分布$q$的期望

2 问题建模

为不失一般性，本文以脉冲重复间隔（PRI）参数定义多功能雷达工作模式的参数估计与变点检测任务，但所提方法可扩展至其他脉冲描述字（PDW）及多参数场景。首先，建立不同PRI调制的概率生成模型，随后给出参数估计与变点检测任务的数学表述。

2.1 基于概率图模型的PRI调制表征

相关文献中常用的PRI调制类型主要有6种[5]。根据隐状态自转移倾向，可将常见PRI类型分为两类：敏捷调制类型与非敏捷调制类型。其中，敏捷调制类型包括参差PRI、滑动PRI、敏捷PRI；非敏捷调制类型包括抖动PRI、驻留切换（D&S）PRI。由于PRI序列是时序数据，需通过到达时间（TOA）序列的一阶差分提取，因此杂波脉冲和缺失脉冲造成的非理想条件会影响其前序脉冲的一阶差值（为与真实PRI值区分，该差值记为$\Delta \text{TOA}$）。脉冲序列中单个缺失脉冲与单个杂波脉冲的影响如图1所示：杂波脉冲会使原始PRI值分裂为两个相邻的较小$\Delta \text{TOA}$值；缺失脉冲则会导致一个较大的$\Delta \text{TOA}$值，该值等于缺失脉冲与其前序脉冲的PRI值之和。为增强非理想条件下的建模能力，本文采用概率图模型表征各类调制类型。

（注：杂波脉冲会使原始PRI值分裂为两个相邻的较小$\Delta \text{TOA}$值；缺失脉冲会导致一个较大的$\Delta \text{TOA}$值，该值等于原始PRI值之和。）

多功能雷达通过排列有限数量的有序脉冲实现特定雷达工作模式。受测量噪声影响，每个脉冲PRI值的观测函数可通过概率密度函数或概率质量函数建模。本文将工作模式建模为具有高斯发射分布的隐马尔可夫模型（G-HMM），该模型可表示为如下三元组：
$$\lambda =(\pi ,A,\Theta )(1)$$
其中，$A$为$K$个状态的状态转移矩阵，$\pi$为初始状态分布，$\Theta$为高斯分布集合。下文定义将明确G-HMM与不同PRI定义的工作模式之间的映射关系。

定义：雷达工作模式的PRI序列可由高斯分布集合$\Theta =({\phi }^1,\dots ,{\phi }^k,\dots ,{\phi }^K)$表征。本文中，${\phi }^k$表示对应第$k$个隐状态的参数化模型；在高斯假设下，${\phi }^k\sim \mathcal{N}({\mu }_k,{\sigma }_k^2)$，其概率密度函数（PDF）定义为：
$$f_{{\phi }^k}(p_t)=\frac{1}{\sqrt{2\pi {\sigma }_k^2}}\exp \left(-\frac{(p_t-{\mu }_k{)}^2}{2{\sigma }_k^2}\right),1\le k\le K(2)$$
其中，$p_t$为时刻$t$接收到的脉冲PRI值；PRI序列遵循底层转移模式$A=(a_j{)}_{j=1}^K$，且$a_j=(a_{ji}{)}_{i=1}^K$。例如，考虑某工作模式采用参差PRI调制，其单个周期内的参差序列为(2, 3, 5)，则在测量噪声（假设噪声方差为0.1）下，接收到的参差序列可由$\Theta =({\phi }^1,{\phi }^2,{\phi }^3)$表征，其中${\phi }^1=({\mu }_1=2,{\sigma }_1^2=0.1)$、${\phi }^2=({\mu }_2=3,{\sigma }_2^2=0.1)$、${\phi }^3=({\mu }_3=5,{\sigma }_3^2=0.1)$；同时，转移概率满足$a_{12}=1$、$a_{23}=1$、$a_{31}=1$，且当ji∉{ 12,23,31}ji \notin \{12, 23, 31\}ji∈/{ 12,23,31}时，$a_{ji}=0$。图2给出了不同调制PRI序列及其对应的马尔可夫链与状态转移矩阵示例。

（a）驻留切换（D&S）PRI；（b）参差PRI；（c）滑动PRI；（d）抖动PRI）

2.2 多功能雷达工作模式参数估计与变点检测

本文将脉冲间调制类型或调制参数的差异定义为相邻工作模式间的变点。假设{ pt}t≥1\{p_t\}_{t \geq 1}{ pt​}t≥1​为参数化序列，其参数为{ λt}t≥1\{\lambda_t\}_{t \geq 1}{ λt​}t≥1​；在变点$\nu$之前，参数保持为${\lambda }_{<\nu }$，在变点之后，参数变为${\lambda }_{\ge \nu }$，具体如式(3)所示：
$$\mathcal{H}(p_t,{\lambda }_t)=\{\begin{array}{ll}{\lambda }_{<\nu }& \text{若 }t<\nu \\ {\lambda }_{\ge \nu }& \text{若 }t\ge \nu \end{array}$$

在传统变点检测研究中，变前与变后分布可能已知（或至少其结构已知）；但在雷达截获应用中，先验信息往往不足，因此本文假设变前与变后分布的结构和参数均未知。变前与变后分布的结构通过所提agile BNP-HMM建模，且变前分布可根据贝叶斯定理从数据中估计得到：
$$p({\lambda }_t|P^t)=\frac{p({\lambda }_t)p(P^t|{\lambda }_t)}{\int p(P^t,{\lambda }_t)d{\lambda }_t}$$
其中，$p({\lambda }_t)$为先验分布，Pt={ pi}i=1tP^t=\{p_i\}_{i=1}^tPt={ pi​}i=1t​为数据批次，$\int p(P^t,{\lambda }_t)d{\lambda }_t$为证据积分。本文通过隐马尔可夫模型表征雷达PRI序列的脉冲间调制特性，具体而言，参数估计任务的目标包括：一个周期内PRI电平的数量$K$（即隐马尔可夫模型的隐状态数量）、每个PRI电平的高斯参数$({\phi }^1,\dots ,{\phi }^K)$，以及时序调制特性$A$（即隐马尔可夫模型的状态转移矩阵）。在推断出这些隐变量后，基于估计的参数序列执行变点检测任务。

为明确检测准则，定义检测器的停止时间$N$为变点告警时间。从在线处理角度，告警时间应始终晚于实际变点时间。本文采用Lorden提出的“最坏情况”平均检测延迟（MDD）[35]：
$${\overline{\tau }}^{\ast }=\underset{t\ge 1}{\sup }\text{ess sup }{\mathbb{E}}_{{\lambda }_{\ge \nu }}(N-\nu |N\ge \nu ,P)$$

在平均虚警时间（MT2FA）约束下，平均检测延迟应尽可能小；平均虚警时间定义为：
$$\overline{T}={\mathbb{E}}_{{\lambda }_{<\nu }}(N)(6)$$
其中，$\text{ess sup}$表示本质上确界。参数估计任务需最小化平均检测延迟${\overline{\tau }}^{\ast }$，同时最大化平均虚警时间$\overline{T}$，即：
$$\text{minimize }\text{MDD}(N)\text{subject to }\text{MT2FA}(N)\ge \gamma (7)$$
其中，$\gamma$为平均虚警时间约束；上述优化问题为极小极大问题。

3 方法

本文提出一种用于多功能雷达工作模式参数估计与变点检测的在线处理框架。本节首先介绍全共轭生成模型（agile BNP-HMM）与新的先验分布（敏捷先验分布），随后阐述在线处理流程及相应实现细节。

3.1 敏捷贝叶斯非参数隐马尔可夫模型构建

在多功能雷达侦察的非协作任务中，工作模式的隐状态数量及对应参数通常未知，需通过推断获取。本文设计一种新颖的敏捷先验分布，用于控制G-HMM隐状态的自转移倾向，以建模具有敏捷脉冲间调制类型的多功能雷达工作模式转移过程。本节首先给出多功能雷达工作模式的生成模型，随后介绍用于实现敏捷调制类型先验分布的改进断棒构造方法，最后基于所提图模型推导变分贝叶斯迭代的闭式解。

3.1.1 敏捷贝叶斯非参数隐马尔可夫模型

agile BNP-HMM的图模型如图3所示，其中包含截获脉冲P={ pt}t≥1P=\{p_t\}_{t \geq 1}P={ pt​}t≥1​、隐状态分配S={ st}t≥1S=\{s_t\}_{t \geq 1}S={ st​}t≥1​、充分统计量$\overline{\Theta }=({\phi }^k{)}_{k=1}^{\infty }$、初始分布$\pi =({\pi }_k{)}_{k=1}^{\infty }$、转移矩阵$A=(a_j{)}_{j=1}^{\infty }$（且$a_j=(a$