文献阅读：Interleaved Hidden Markov Processes Inference for Deinterleaving Radar Pulse Sequences

用于雷达脉冲序列解交织的交织隐马尔可夫过程推断

摘要：隐马尔可夫过程（HMP）已被广泛用于建模雷达脉冲序列。在电子侦察系统的雷达信号解交织任务中，截获的雷达脉冲序列被假设为交织隐马尔可夫过程（IHMP）。在此背景下，本文提出一种生成模型来表示交织隐马尔可夫过程，并将解交织问题重新表述为后验推断任务。为计算后验概率，本文首先设计了一种精确推断算法。然而，由于隐状态表示的组合特性，精确推断在计算上具有不可解性。为解决这一局限性，本文进一步提出了一种基于采样的方法和两种基于变分的方法，为后验计算提供了可解方案。最后，基于似然比检验推导了错误概率的理论下界，且所提方法被证明能较好地接近该下界。在多种雷达脉冲信号数据集上的仿真验证表明，采用结构化近似的变分推断在解交织准确性和计算效率之间实现了更优的平衡，使其成为精确推断方法和基于搜索的方法的一种极具潜力的替代方案。

关键词：解交织；期望最大化；隐马尔可夫模型；马尔可夫链蒙特卡洛；变分推断

一、引言

隐马尔可夫过程[1]是一类离散时间有限状态齐次马尔可夫链，可通过离散时间无记忆不变信道观测得到。该信道由一组以马尔可夫链状态为索引的有限转移密度表征。隐马尔可夫过程通常用隐马尔可夫模型（HMM）表示。然而，在侦察系统[2,3]中，观测结果由多个隐马尔可夫过程的输出构成。本文研究如下问题：

设侦察系统的观测序列为p={pt}t=1Tp = \{p_t\}_{t=1}^Tp={pt​}t=1T​，其中下标$t$表示到达顺序（AO），上标$T$表示序列长度。该观测序列由多个隐马尔可夫模型的交织输出组成。交织机制通过马尔可夫切换过程$P_{sw}$建模。存在$M$个独立分量随机过程$P_1,\dots ,P_m,\dots ,P_M$，每个过程均用一个隐马尔可夫模型建模。设${\phi }^1,\dots ,{\phi }^m,\dots ,{\phi }^M$为一组非空有限集合，定义字母表$A={\phi }^1\cup \cdots \cup {\phi }^m\cup \cdots \cup {\phi }^M$。称$A$为字母表，每个${\phi }^m$为第$m$个子字母表。第$m$个隐马尔可夫模型定义在子字母表${\phi }^m$上，设φm={φ1m,…,φkm,…,φKmm}\varphi^m = \{\varphi_1^m, \dots, \varphi_k^m, \dots, \varphi_{K^m}^m\}φm={φ1m​,…,φkm​,…,φKmm​}，其中${\phi }_k^m$表示${\phi }^m$中的第$k$个字母，$K^m$表示${\phi }^m$的规模。每个字母被建模为高斯随机变量。设Π={φ1,…,φm,…,φM}\Pi = \{\varphi^1, \dots, \varphi^m, \dots, \varphi^M\}Π={φ1,…,φm,…,φM}表示字母表$A$的一个划分，且$P_{sw}$定义在字母表$\Pi$上。

基于上述定义，多个隐马尔可夫过程的交织过程$I(P_1,\dots ,P_M,P_{sw})$的生成方式如下：在每个时间步$t$，切换过程$P_{sw}$生成一个索引$m$；随后，分量过程$P_m$被激活[1]，并从字母${\phi }_k^m$中采样输出一个观测值，其中字母服从概率密度函数（pdf）$f_{{\phi }_k^m}(\cdot )$，且$k$对应于当前时刻分量过程$P_m$的隐状态索引。假设分量过程与切换过程相互独立。解交织问题即根据观测序列$p$重构由$P_{sw}$生成的索引序列。

图1 含两个分量过程的交织隐马尔可夫过程及解交织任务示意图

1.1 研究动机

上述问题源于电子侦察系统中的雷达脉冲解交织任务[4]。电子侦察接收机旨在截获并分析多个雷达源发射的雷达脉冲序列，需区分每个脉冲并将其归属到对应的辐射源。不同雷达的脉冲按时间顺序被观测，每个脉冲由脉冲描述字（PDW）表征，PDW包含到达时间（TOA）、脉冲宽度（PW）、射频（RF）和脉冲幅度（PA）等参数，其含义如图2所示。本文中，雷达发射的脉冲序列用隐马尔可夫过程表示，对应分量过程$P_m$；子字母表${\phi }^m$表示某一雷达脉冲序列的PDW值集合；切换过程$P_{sw}$反映多个雷达脉冲序列的交织机制。

图2 某雷达发射的两个脉冲示例，下标表示每个脉冲的到达顺序

1.2 相关工作

文献中已提出多种雷达脉冲序列解交织技术。传统方法通过累积量差直方图（CDIF）[5]、序贯差直方图（SDIF）[6]等技术识别脉冲重复间隔（PRI）模式。此外，文献[7]提出一种基于最优传输距离的凝聚层次聚类方法，利用PA和TOA信息关联PW与RF聚类。然而，这些方法的表示能力有限，无法有效捕捉信号序列中的潜在模式，在更复杂环境下性能可能下降。

近年来，基于模型的方法受到广泛关注。这类方法主要利用概率模型表示脉冲序列中的潜在模式。具体而言，文献[8]采用高斯模型表示单个雷达发射的脉冲序列，再通过高斯混合模型（GMM）对混合脉冲序列建模；文献[9]提出贝叶斯非参数高斯混合模型用于混合脉冲序列建模，并采用变分推断实现解交织，验证了变分推断在在线解交织中的有效性。但这些方法不适用于采用多射频或多脉冲宽度的捷变雷达。

针对雷达采用多射频或多脉冲宽度的复杂场景，相关研究已逐步展开。一种极具潜力的方法是利用马尔可夫链增强雷达脉冲序列的表示能力。文献[10]提出交织马尔可夫链的概念，随后Seroussi将其扩展为交织马尔可夫过程（IMP）[11,12,13]；文献[14]将交织马尔可夫过程用于雷达脉冲解交织，采用到达顺序（AO）而非到达时间（TOA）；为进一步利用可用信息，文献[15]采用更新过程表示TOA序列，文献[16]在此基础上利用马尔可夫更新过程同时融合PDW与TOA数据。然而，该研究方向存在以下局限性：第一，将字母视为确定性变量，难以适应电磁环境中常见的观测噪声；第二，不输出每个脉冲的标签，在非不交子字母表场景中的适用性受限；第三，依赖启发式优化技术，计算负担较重。

1.3 研究贡献

本文旨在构建具有强表示能力的模型并开发高效的解交织算法。首先，定义交织隐马尔可夫过程以表示混合雷达脉冲序列，并设计相应的生成模型；其次，将解交织问题转化为后验推断问题，设计基于期望最大化（EM）算法的精确推断方法，但发现由于组合特性，精确推断具有NP难特性；为解决该问题，设计了三种近似推断方法：基于马尔可夫链蒙特卡洛（MCMC）的算法和两种基于变分的算法；最后，基于似然比检验推导了解交织二状态隐马尔可夫模型的错误概率理论下界。本文的主要贡献总结如下：

1. 提出交织隐马尔可夫过程的生成模型，该模型能够表示观测噪声并适应非不交子字母表场景。
2. 提出一种精确推断方法和三种近似推断方法用于计算后验概率，这些方法均经过明确推导，显著提升了解交织效率。
3. 推导了二状态隐马尔可夫模型混合解交织的错误概率理论下界，仿真结果表明所提方法能较好地接近该下界。
4. 在雷达数据集上的仿真验证了所提方法的有效性，尤其证明了采用结构化近似的变分推断在准确性与计算效率之间实现了良好平衡。

1.4 论文结构

本文其余部分安排如下：第二章介绍脉冲间调制表示初步知识及交织隐马尔可夫过程的生成模型；第三章提出四种解交织方法；第四章分析错误概率；第五章给出在多种雷达数据集上的仿真结果；第六章总结全文。

[1] 当隐马尔可夫过程被激活时，其隐状态根据第二章A节所述的转移概率演化。

[2] 到达顺序（AO）是到达时间（TOA）的简化形式，脉冲按TOA值排序，排序后的索引即为AO。

二、预备知识与生成模型

2.1 雷达脉冲间调制表示初步

雷达通过排列有限数量的脉冲实现特定雷达任务，这些脉冲的动态模式称为脉冲间调制。常用的脉冲间调制类型有六种，包括参差、恒定、抖动、捷变、驻留切换（D&S）和滑动[17,18]。

受测量噪声影响，每个脉冲PDW值的观测函数采用参数模型建模。雷达脉冲序列用具有参数输出的隐马尔可夫模型（即高斯输出隐马尔可夫模型）表示，该模型可表示为三元组：
$$({\pi }^1,A^1,{\phi }^1)\text{\hspace{1em}(1)}$$
式中，上标“1”表示本节描述的单隐马尔可夫模型场景。在式（1）中，π1={πi1}i=1K1\pi^1 = \{\pi_i^1\}_{i=1}^{K^1}π1={πi1​}i=1K1​为规模为$K^1\times 1$的初始分布；A1={Aj1}j=1K1A^1 = \{A_j^1\}_{j=1}^{K^1}A1={Aj1​}j=1K1​为状态转移矩阵，其中Aj1={Aj,i1}i=1K1A_j^1 = \{A_{j,i}^1\}_{i=1}^{K^1}Aj1​={Aj,i1​}i=1K1​；隐状态$S_t^1$为规模为$K^1\times 1$的独热向量。隐状态根据初始分布和转移矩阵演化：
$$\begin{array}{rl}P(S_1^1)& ={\pi }^1,\\ P(S_t^1\mid S_{t-1}^1)& =S_{t-1}^{1\top }{A}^1{S}_t^1,\end{array}$$
式中，$v^{\top }$表示向量$v$的转置。时间$t$的观测PDW记为$p_t$，由高斯输出隐马尔可夫模型生成。该隐马尔可夫模型定义在字母表φ1={φ11,…,φk1,…,φK11}\varphi^1 = \{\varphi_1^1, \dots, \varphi_k^1, \dots, \varphi_{K^1}^1\}φ1={φ11​,…,φk1​,…,φK11​}上，其中${\phi }_k^1$为字母表中的第$k$个字母。从生成角度看：在每个时间$t$，隐状态$S_t^1$根据式（2）演化并输出索引$k$；随后从随机变量${\phi }_k^1$中采样得到观测值$p_t$，且${\phi }_k^1$服从高斯概率密度函数$f_{{\phi }_k^1}(\cdot )$：
$$f_{{\phi }_k^1}(p_t)\propto \exp \left(-\frac{1}{2}(p_t-{\mu }_k^1{)}^{\top }({\Sigma }_k^1{)}^{-1}(p_t-{\mu }_k^1)\right)\text{\hspace{1em}(3)}$$
式中，${\mu }_k^1$和${\Sigma }_k^1$分别为多元高斯分布的均值和协方差。该生成模型的概率表示如图3所示。

图3 高斯输出隐马尔可夫模型的图表示，由于本节场景中观测仅由一个隐马尔可夫模型生成，因此变量和参数的上标均为1

PDW序列遵循由转移矩阵$A$描述的潜在转移模式。例如，考虑某参差射频调制脉冲序列，其一个周期内的射频值为$(3,2,5)\text{MHz}$，则字母表定义为φ1={φ11,φ21,φ31}\varphi^1 = \{\varphi_1^1, \varphi_2^1, \varphi_3^1\}φ1={φ11​,φ21​,φ31​}，其中${\phi }_1^1=({\mu }_1=2,{\sigma }_1^2=1)$、${\phi }_2^1=({\mu }_2=3,{\sigma }_2^2=1)$、${\phi }_3^1=({\mu }_3=5,{\sigma }_3^2=1)$；同时，转移矩阵满足$A_{2,1}=1$、$A_{1,3}=1$、$A_{3,2}=1$，其余元素均为0。图4展示了三种不同类型的脉冲间调制。

图4 三种不同脉冲间调制类型示例：（a）参差调制；（b）驻留切换（D&S）调制；（c）滑动调制。左列展示三种调制类型的转移图，中列展示转移概率，右列展示脉冲间调制的序列示例。不同颜色代表雷达使用的不同PDW值（本文中为射频值）

此外，当检测概率为1时，在某些调制模式（如滑动或参差调制）下，观测脉冲序列可能呈现严格周期性，导致马尔可夫链具有周期性，从而不满足遍历性（遍历性需满足非周期性、不可约性和正常返性）[19]。在电磁环境中，检测概率通常小于1，会导致观测缺失；这些随机缺失会破坏PDW序列的周期性结构，使观测过程具备非周期性。因此，本文假设最终得到的马尔可夫链是遍历的。

2.2 交织隐马尔可夫过程的生成模型

多个雷达的观测序列构成交织隐马尔可夫过程。在实际场景中，准确描述交织机制具有复杂性。本文采用马尔可夫切换过程$P_{sw}$对交织机制建模，并采用一阶马尔可夫链作为该切换过程的生成模型。假设$P_{sw}$的初始相位服从均匀分布，所提模型的图表示如图5所示。

图5 交织隐马尔可夫过程的图模型，包含$M$个分量过程$P_1,\dots ,P_M$，每个分量过程的底层生成模型为隐马尔可夫模型。第$m$个分量过程在时间$t$的隐状态记为$S_t^m$；由马尔可夫链建模的切换过程$P_{sw}$控制每个时间步的激活分量，其在时间$t$的状态记为$Z_t$

图5中包含$M$个分量过程$P_1,\dots ,P_M$，每个分量过程的底层生成模型为隐马尔可夫模型。第$m$个分量过程的隐状态记为$S_t^m$，规模为$K^m\times 1$，隐状态根据规模为$K^m\times 1$的初始分布${\pi }^m$和规模为$K^m\times K^m$的转移矩阵$A^m$演化。此外，存在一个由一阶马尔可夫链建模的切换过程$P_{sw}$，其在时间$t$的状态记为$Z_t$（规模为$M\times 1$），并根据规模为$M\times 1$的初始分布${\pi }^z$和规模为$M\times M$的转移矩阵$A^z$演化。状态变量$S_t^m$和$Z_t$均为独热向量。该生成模型的联合概率表示为：
P(Z,S,p,Az,A,πz,π,A)∝P(Z1∣πz)×∏t=2TP(Zt∣Zt−1,Az)×∏m=1M{P(S1m∣πm)∏t=1TP(pt∣Stm,Zt,m,φm)∏t=2TP(Stm∣St−1m,Zt,m,Am)},(4) \begin{aligned} &P(Z, S, p, A^z, A, \pi^z, \pi, \mathcal{A}) \propto \\ &P(Z_1 \mid \pi^z) \times \prod_{t=2}^T P(Z_t \mid Z_{t-1}, A^z) \times \prod_{m=1}^M \left\{ P(S_1^m \mid \pi^m) \right. \\ &\left. \prod_{t=1}^T P(p_t \mid S_t^m, Z_{t,m}, \varphi^m) \prod_{t=2}^T P(S_t^m \mid S_{t-1}^m, Z_{t,m}, A^m) \right\}, \end{aligned} \tag{4} ​P(Z,S,p,Az,A,πz,π,A)∝P(Z1​∣πz)×t=2∏T​P(Zt​∣Zt−1​,Az)×m=1∏M​{P(S1m​∣πm)t=1∏T​P(pt​∣Stm​,Zt,m​,φm)t=2∏T​P(Stm​∣St−1m​,Zt,m​,Am)},​(4)
式中，A={Am}m=1MA = \{A^m\}_{m=1}^MA={Am}m=1M​、Z={Zt}t=1TZ = \{Z_t\}_{t=1}^TZ={Zt​}t=1T​、p={pt}t=1Tp = \{p_t\}_{t=1}^Tp={pt​}t=1T​、π={πm}m=1M\pi = \{\pi^m\}_{m=1}^Mπ={πm}m=1M​，$Z_{t,m}$表示状态变量$Z_t$的第$m$个元素。具体而言，式（4）右侧第一、二项分别对应切换过程的初始概率和状态转移概率，满足以下马尔可夫性：
$$\begin{array}{rl}P(Z_1)& ={\pi }^z,\\ P(Z_t\mid Z_{t-1})& =Z_{t-1}^{\top }{A}^z{Z}_t.\end{array}\text{\hspace{1em}(5)}$$

式（4）右侧第三项对应分量过程的概率，因此有：
$$P(S_1^m)={\pi }^m,\text{\hspace{1em}(6)}$$
$$\begin{array}{c}P(S_t^m\mid S_{t-1}^m,Z_t)=S_{t-1}^{m\top }\left(Z_{t,m}{A}^m+{\overline{Z}}_{t,m}{I}^m\right)S_t^m,\\ P(p_t\mid S_t^m,Z_t,{\phi }^m)\propto \\ \exp \left(-\frac{1}{2}(p_t-{\mu }_t{)}^{\top }({\Sigma }_t{)}^{-1}(p_t-{\mu }_t)\right),\end{array}\text{\hspace{1em}(7)}$$
式中，${\overline{Z}}_{t,m}=1-Z_{t,m}$，$I^m$为规模为$K^m\times K^m$的单位矩阵，${\mu }_t$和${\Sigma }_t$分别为第$t$个观测字母${\phi }_t$的均值和协方差矩阵。字母${\phi }_t$定义为：
$${\phi }_t=\sum _{m=1}^M{Z}_{t,m}\sum _{k=1}^K{S}_{t,k}^m{\phi }_k^m.\text{\hspace{1em}(8)}$$

由式（5）和（6）可知，切换过程选择一个对应于某一分量过程的索引；当分量过程被激活时，其隐状态$S_t^m$根据状态转移矩阵$A^m$演化，否则隐状态保持不变（状态保持等效于根据单位矩阵$I^m$进行状态转移）。

2.3 解交织问题表述

基于交织隐马尔可夫过程的生成模型，根据贝叶斯定理，解交织问题可转化为后验推断问题：
$$P(Z,S,\Gamma \mid p)=\frac{P(Z,S,\Gamma )P(p\mid Z,S,\Gamma )}{\int P(p,Z,S,\Gamma )d(Z,S,\Gamma )},\text{\hspace{1em}(9)}$$
式中，模型参数Γ={A,Az,π,πz}\Gamma = \{A, A^z, \pi, \pi^z\}Γ={A,Az,π,πz}，$P(Z,S,\Gamma )$为先验分布。解交织结果由$Z^{\ast }$体现，通过最大化后验概率得到：
$$Z^{\ast }=\arg \underset{Z,S,\Gamma }{\max }P(Z,S,\Gamma \mid p).\text{\hspace{1em}(10)}$$

三、方法

在概率模型中，状态推断和参数学习是两个核心问题。状态推断旨在计算给定观测下隐状态的概率，参数学习则聚焦于根据观测估计模型参数。本文提出基于期望最大化（EM）的方法求解后验概率，这类EM类算法通过迭代执行以下两步：模型参数$\Gamma$估计（M步）和隐状态分配$Z,S$推断（E步）。

标准EM算法的目标是最大化隐变量与观测变量的期望联合对数似然[1]，该似然函数定义为：
$$\mathbb{Q}({\Gamma }^{\prime }\mid \Gamma )=\mathbb{E}\left(\log P(S,Z,p\mid {\Gamma }^{\prime })\mid \Gamma ,p\right),\text{\hspace{1em}(11)}$$
式中，$\mathbb{E}(\cdot )$表示求期望，${\Gamma }^{\prime }$为当前迭代中估计的新参数。E步计算式（11）定义的$\mathbb{Q}$函数，M步通过更新参数$\Gamma$最大化$\mathbb{Q}$函数。根据Jensen不等式，文献[20]已证明，经过每一轮E步和M步迭代，似然函数值会单调递增，直至收敛到局部最优。M步通过对$\mathbb{Q}$函数（式11）关于参数求导并令导数为零实现；初始分布$\pi$和转移矩阵$A$的M步更新采用经验估计器推导。关于$\mathbb{Q}$函数的详细描述及参数更新函数的推导过程见附录A。

以下小节将介绍四种用于执行E步的方法，并讨论它们各自的优缺点。为简化表示与说明，假设每个字母的协方差矩阵相同，即对所有$i\in [1,K^m]$和$m\in [1,M]$，有${\Sigma }_i^m=\Sigma$。

3.1 精确推断

精确推断方法的设计受文献[21]启发，但文献[21]假设切换过程无记忆性，因此表示能力有限。为实现交织隐马尔可夫过程生成模型的精确推断，本文将图5所示的图模型转化为等效的“大隐马尔可夫模型”。该等效隐马尔可夫模型的隐状态数量为$M{\prod }_{m=1}^M{K}^m$，因此可采用标准的前向-后向算法[20]。转化规则基于$M$个分量过程隐状态的所有可能组合。例如，若存在两个分量过程，每个过程均用二状态隐马尔可夫模型建模，则每个时间步存在8种可能的隐状态组合，因此等效隐马尔可夫模型将包含8个隐状态，其隐状态空间为：
$$\begin{array}{rl}& (\arg \max Z_t,\arg \max S_t^1,\arg \max S_t^2)\in \\ & \{(1,1,1),(1,1,2),(1,2,1),(1,2,2),\\ & (2,1,1),(2,1,2),(2,2,1),(2,2,2)\}.\end{array}$$

针对该等效隐马尔可夫模型的前向-后向算法，空间复杂度为$O\left(TM{\prod }_{m=1}^M{K}^m\right)$，时间复杂度为$O\left(T{M}^2{\prod }_{m=1}^M(K^m{)}^2\right)$。由于时间复杂度中包含$K^m$的乘积项，导致其呈现指数复杂度，使得精确推断在计算上不可解。尽管前向-后向算法能高效计算似然，但无法避免每个时间步分量过程隐状态所有可能组合的乘积运算。后文将该精确推断方法简称为“EM方法”。

为避免隐状态数量呈指数级增长，概率图模型中的近似推断受到广泛关注[22]。以下小节将介绍基于马尔可夫链蒙特卡洛（MCMC）和基于变分的近似推断算法。

3.2 马尔可夫链蒙特卡洛推断

马尔可夫链蒙特卡洛算法可解决计算不可解模型的推断问题[23]。现有多种采样方法，本文选择采用带祖先采样的粒子吉布斯（PGAS）方法[24]，该方法属于粒子马尔可夫链蒙特卡洛（PMCMC）方法家族[25]。

在PGAS方法中，通过运行序贯蒙特卡洛（SMC）采样器构建马尔可夫核，其中一个粒子被设为确定性参考输入粒子$r_t$，该参考粒子对应前一次迭代的输出。新的参考轨迹从粒子轨迹中根据其重要性权重采样得到。具体而言，每个时间步设置$N$个粒子，代表所提生成模型的隐状态；每个粒子用规模为$(M+1)\times 1$的向量表示，每个元素对应$M+1$个马尔可夫链的隐状态索引；第$i$个粒子记为$x_t(i)$，第$N$个粒子被确定性地设为参考粒子$r_t$。同时，PGAS算法利用祖先索引at(i)∈{1,…,N}a_t(i) \in \{1, \dots, N\}at​(i)∈{1,…,N}实现粒子传播，粒子轨迹可递归定义为$x_{1:t}(i)=\text{concat}(x_{1:t-1}(a_t(i)),x_t(i))$（$t>1$），其中$\text{concat}$表示拼接操作。

采样过程描述如下：在每个时间步$t$，前$N-1$个粒子的祖先索引从分类分布中采样，即$a_t(i)\sim \text{Categorical}(w_{t-1}(1),\dots ,w_{t-1}(i),\dots ,w_{t-1}(N))$，其中重要性权重$w_t(i)$定义为：
$$w_t(i)=P(p_t\mid x_t(i)),\text{\hspace{1em}(12)}$$
式中，$P(p_t\mid x_t(i))$为给定隐状态下$p_t$的条件概率。根据采样得到的祖先索引，粒子按照转移规则$P(x_t\mid x_{t-1})$传播。同时，第$N$个粒子设为参考粒子，其祖先索引$a_t(N)$从分类分布$a_t(N)\sim \text{Categorical}({\tilde{w}}_{t|T}(1),\dots ,{\tilde{w}}_{t|T}(i),\dots ,{\tilde{w}}_{t|T}(N))$中采样，其中：
$${\tilde{w}}_{t|T}(i)\propto w_{t-1}(i)P(x_t(N)\mid x_{t-1}(i)),\text{\hspace{1em}(13)}$$
式中，$P(x_t(N)\mid x_{t-1}(i))$为转移概率的乘积。算法1总结了PGAS方法的流程。

算法1 带祖先采样的粒子吉布斯（PGAS）
输入：时间$t=1,\dots ,T$的参考轨迹$r_t$，以及参数$\Gamma$
输出：采样轨迹$r_{out}$

1. 根据式（6）采样$x_1(i)$（$i=1,\dots ,N$）
2. 设置$x_1(N)=r_1$
3. 根据式（12）计算$w_1(i)$（$i=1,\dots ,N$）
4. 对于$t=2,\dots ,T$：
5. 采样$a_t(i)\sim \text{Categorical}(w_{t-1}(1),\dots ,w_{t-1}(N))$（$i=1,\dots ,N-1$）
6. 根据式（13）计算${\tilde{w}}_{t-1|T}(i)$（$i=1,\dots ,N$）
7. 采样$a_t(N)\sim \text{Categorical}({\tilde{w}}_{t-1|T}(1),\dots ,{\tilde{w}}_{t-1|T}(N))$
8. 根据式（6）采样$x_t(i)$（$i=1,\dots ,N-1$）
9. 设置$x_t(N)=r_t$
10. 对$i=1,\dots ,N$，设置$x_{1:t}(i)=\text{concat}(x_{1:t-1}(a_t(i)),x_t(i))$
11. 根据式（12）计算$w_t(i)$（$i=1,\dots ,N$）
12. 结束循环
13. 采样$k\sim \text{Categorical}(w_T(1),\dots ,w_T(N))$
14. 返回$r_{out}=x_{k,1:T}$

PGAS算法单次迭代的时间复杂度为$O(NTM)$，空间复杂度为$O((M+2)NT)$。然而，基于马尔可夫链蒙特卡洛的方法耗时较长，因为只有当马尔可夫链经过预热期[26]后，采样过程才有效，即需经过一定迭代后，采样轨迹才能作为近似推断结果。此外，还需计算每个粒子的权重，进一步增加了计算负担。因此，基于马尔可夫链蒙特卡洛的方法不适用于电子侦察场景。

3.3 变分推断

变分方法可解决计算不可解模型的推断问题，并为后验计算提供高效方案。变分推断的核心思想是引入变分分布$Q$，最小化真实后验分布$P$与变分分布$Q$之间的距离。根据文献[27]，观测序列的对数似然函数可表示为：
$$

\log P(p \mid Z, S, \Gamma) =
\mathcal{L}(Q(Z, S, \Gamma)) + \text{KL}(Q(Z, S, \Gamma) \parallel P(Z, S, \Gamma \mid p)), \tag{14}

KaTeX parse error: Can't use function '$' in math mode at position 5: 式中，$̲\mathcal{L}$表示证…
\mathcal{L} = \mathbb{E}_Q\left{ \log P(\Gamma, S, Z, p) \right} - \mathbb{E}_Q\left{ \log Q(S, Z, \Gamma) \right}, \tag{15}
$$
式中，${\mathbb{E}}_Q$表示关于变分分布$Q$的期望。根据Blei[28]的研究，变分推断的复杂度由条件独立关系决定。因此，需构建可计算的变分结构，通过最大化证据下界（式15）优化变分参数，以获得最紧的下界。本节将明确推导两种基于不同变量耦合程度的变分推断方法：基于平均场假设的变分推断和基于结构化假设的变分推断，两种近似的图表示如图6所示。

图6 当$M=2$时，平均场近似与结构化近似的图表示：（a）平均场近似；（b）结构化近似

3.3.1 平均场变分推断（MFVI）

一种简单的选择是假设所有隐变量相互独立，如图6（a）所示。此时变分分布定义为：
$$Q(S,Z,\Gamma )=Q(\Gamma )\prod _{t=1}^T\prod _{m=1}^{M}Q(S_t^m\mid {\theta }_t^m)\prod _{t=1}^{T}Q(Z_t\mid {\varphi }_t),\text{\hspace{1em}(16)}$$
式中，${\theta }_t^m$和${\varphi }_t$分别为隐状态$S_t^m$和$Z_t$的期望。基于上述符号，可明确表示变分分布为：
$$Q(S_t^m\mid {\theta }_t^m)=\prod _{k=1}^{K^m}({\theta }_{t,k}^m{)}^{S_{t,k}^m},Q(Z_t\mid {\varphi }_t)=\prod _{m=1}^M({\varphi }_{t,m}{)}^{Z_{t,m}}.\text{\hspace{1em}(17)}$$

为简化计算，设$Q(\Gamma )=1$。令${\overline{A}}_t^m={\overline{Z}}_{t,m}\log I^m+Z_{t,m}\log A^m$，则隐状态的更新函数如式（18）和（19）所示。式中，${\Delta }_m^n={\mu }^m{\Sigma }^{-1}{\mu }^{n\top }$；$\text{idiag}(\cdot )$表示对输入方阵取对角线元素构成向量的算子；$\text{diag}(\cdot )$表示将输入向量作为对角线元素构成方阵的算子；$\text{tr}(\cdot )$表示返回输入矩阵迹的算子。证据下界的表达式及式（18）、（19）的推导过程见附录B。

尽管采用平均场假设近似隐变量的后验分布，但仍保留了时间序列依赖性，式（18）和（19）的最后两项体现了这种时间依赖性。每个隐状态向量根据式（18）和（19）更新，单次迭代的时间复杂度为$O\left(T{\sum }_{m=1}^M(K^m{)}^2+T{M}^2\right)$，空间复杂度为$O\left(T{\sum }_{m=1}^M{K}^m+MT\right)$。

3.3.2 结构化变分推断（SVI）

平均场变分推断假设所有状态变量统计独立，这一假设过于严格。另一种思路是在保证计算可解性的同时保留模型结构。为此，本文保留横向耦合（如图6（b）所示），可采用前向-后向算法计算数据似然。相应的变分分布定义为：
$$\begin{array}{rl}Q(S,Z,\Gamma )& =Q(\Gamma )\prod _{m=1}^M\left\{Q(S_1^m)\prod _{t=2}^{T}Q(S_t^m\mid S_{t-1}^m)\right\}\\ & \times Q(Z_1)\prod _{t=2}^{T}Q(Z_t\mid Z_{t-1}).\end{array}\text{\hspace{1em}(20)}$$

与前文一致，设$Q(\Gamma )=1$。初始分布和转移分布可表示为：
$$\begin{array}{rl}& Q(S_1^m)=\prod _{k=1}^{K^m}(h_{1,k}^m{\pi }_k^m{)}^{S_{1,k}^m};Q(Z_1)=\prod _{k=1}^M(g_{1,k}{\pi }_{1,k}^z{)}^{Z_{1,k}},\\ & Q(S_t^m\mid S_{t-1}^m)=\prod _{j=1}^{K^m}{\left(\prod _{i=1}^{K^m}(h_{t,i}^m{A}_{j,i}^m{)}^{S_{t,i}^m}\right)}^{S_{t-1,j}^m},\\ & Q(Z_t\mid Z_{t-1})=\prod _{j=1}^M{\left(\prod _{i=1}^M(g_{t,i}{A}_{j,i}^z{)}^{Z_{t,i}}\right)}^{Z_{t-1,j}},\end{array}\text{\hspace{1em}(21)}$$
式中，$h_t^m=[h_{t,1}^m,\dots ,h_{t,K^m}^m{]}^{\top }$为规模$K^m\times 1$的向量，$g_t=[g_{t,1},\dots ,g_{t,M}{]}^{\top }$为规模$M\times 1$的向量，均表示观测概率。与平均场变分推断类似，定义${\theta }_t^m$和${\varphi }_t$分别为隐状态$S_t^m$和$Z_t$的期望；观测概率$h_t^m$和$g_t$可通过对${\theta }_t^m$和${\varphi }_t$应用前向-后向算法得到。基于上述表述，通过调整$h_t^m$和$g_t$可最大化证据下界，其更新函数如式（22）和（23）所示（推导过程见附录C）。直观上，结构化变分推断将$M+1$个马尔可夫链解耦，为每个状态变量分配独立的观测。结构化变分推断单次迭代的时间复杂度为$O\left(T{\sum }_{m=1}^M(K^m{)}^2+T{M}^2\right)$，空间复杂度为$O\left({\sum }_{m=1}^M{K}^{m}T+MT\right)$。

θtm=exp⁡{μmΣ−1pt⊤ϕt,m−12∑n=1n≠mMΔmn(ϕt,nθtn)ϕt,m−12idiag{Δmmϕt,m}+EQ(A‾tm)θt−1m+EQ(A‾tm)⊤θt+1m},(18) \theta_t^m = \exp\left\{ \mu^m \Sigma^{-1} p_t^\top \phi_{t,m} - \frac{1}{2} \sum_{\substack{n=1 \\ n \neq m}}^M \Delta_m^n (\phi_{t,n} \theta_t^n) \phi_{t,m} - \frac{1}{2} \text{idiag}\{\Delta_m^m \phi_{t,m}\} + \mathbb{E}_Q(\overline{A}_t^m) \theta_{t-1}^m + \mathbb{E}_Q(\overline{A}_t^m)^\top \theta_{t+1}^m \right\}, \tag{18} θtm​=exp⎩⎨⎧​μmΣ−1pt⊤​ϕt,m​−21​n=1n=m​∑M​Δmn​(ϕt,n​θtn​)ϕt,m​−21​idiag{Δmm​ϕt,m​}+EQ​(Atm​)θt−1m​+EQ​(Atm​)⊤θt+1m​⎭⎬⎫​,(18)

$${\varphi }_{t,m}=\exp \left\{p_t{\Sigma }^{-1}{\mu }^{m\top }{\theta }_t^m-\frac{1}{2}\text{tr}\left({\Delta }_m^m\text{diag}({\theta }_t^m)\right)-\frac{1}{2}\text{tr}\left\{{\Delta }_m^n{\varphi }_{t,n}{\theta }_t^n{\theta }_t^{m\top }\right\}+2\log A_{m,m}^z+\sum _{\begin{array}{c}n=1\\ n\ne m\end{array}}^M{\varphi }_{t,n}\log A_{m,n}^z\right\}.\text{\hspace{1em}(19)}$$

$$

h_t^m = \exp\left{ \mu^m \Sigma^{-1} p_t^\top \phi_{t,m} - \frac{1}{2} \text{diag}{\Delta_m^m \phi_{t,m}} + (1 - \phi_{t,m}) \left[ (\log I^m - \log A^m) \theta_{t-1}^m + (\log I^m - \log A^m)\top \theta_{t+1}^m \right] \right. \
\left. - \frac{1}{2} \sum_{\substack{n=1 \ n \neq m}}^M \Delta_m^n (\phi_{t,n} \theta_t^n) \phi_{t,m} \right}, \tag{22}

$$

gt,m=exp⁡{ptΣ−1μm⊤θtm−12tr{Δmmdiag(θtm)}−12∑n=1n≠mMtr{Δmnϕt,nθtnθtm⊤}−θtm⊤(log⁡Im−log⁡Am)θt−1m}.(23) g_{t,m} = \exp\left\{ p_t \Sigma^{-1} \mu^{m\top} \theta_t^m - \frac{1}{2} \text{tr}\left\{ \Delta_m^m \text{diag}(\theta_t^m) \right\} - \frac{1}{2} \sum_{\substack{n=1 \\ n \neq m}}^M \text{tr}\left\{ \Delta_m^n \phi_{t,n} \theta_t^n \theta_t^{m\top} \right\} - \theta_t^{m\top} (\log I^m - \log A^m) \theta_{t-1}^m \right\}. \tag{23} gt,m​=exp⎩⎨⎧​pt​Σ−1μm⊤θtm​−21​tr{Δmm​diag(θtm​)}−21​n=1n=m​∑M​tr{Δmn​ϕt,n​θtn​θtm⊤​}−θtm⊤​(logIm−logAm)θt−1m​⎭⎬⎫​.(23)

四、错误分析

本节基于最大似然准则，推导二状态隐马尔可夫模型混合解交织的错误概率下界。该下界用于评估所提算法的性能，且可推广到多隐马尔可夫模型混合场景及非二状态场景。为推导错误概率下界，首先引入以下两个假设：

假设1：每个马尔可夫链均为遍历链，在两个状态上具有唯一平稳分布，平稳分布为${\Xi }^m=[{\xi }_1^m,{\xi }_2^m]$（$m\in [1,2]$）和${\Xi }^z=[{\xi }_1^z,{\xi }_2^z]$；子字母表互不相交；每个马尔可夫链的初始状态从平稳分布中采样。

假设2：每个隐马尔可夫模型定义在某一字母表上，字母表中的每个字母为一元高斯随机变量，服从$\mathcal{N}({\mu }_k^m,{\sigma }_k^m)$分布；且每个字母的标准差（SD）相同，即对m∈{1,2}m \in \{1,2\}m∈{1,2}和k∈{1,2}k \in \{1,2\}k∈{1,2}，有${\sigma }_k^m=\sigma$。

在遍历的参数输出隐马尔可夫模型中，平稳输出为该参数分布族的混合分布[29]。因此，每个遍历隐马尔可夫模型等效于高斯混合模型：
$$p_t^m\sim \sum _{k=1}^2{\xi }_k^m\mathcal{N}(p_t^m;{\mu }_k^m,\sigma ),$$
式中，$p_t^m$为第$m$个隐马尔可夫模型在时间$t$的输出。

4.1 理论错误下界

一般而言，两个二状态隐马尔可夫模型生成的观测可能以多种方式混淆。由于本文旨在推导下界，因此聚焦于某一类特定的错误事件。考虑错误事件${\mathcal{E}}_i^{xy}$：第$i$个观测$p_i$由第$y$个隐马尔可夫模型生成，但被错误地识别为第$x$个隐马尔可夫模型的输出，且前一个观测由第$x$个隐马尔可夫模型正确生成并正确识别。由于推导的是下界，因此忽略所有其他潜在错误事件。错误概率可表示为：
$$\begin{array}{rl}{P}_e& \ge \underset{n\to \infty }{\lim }\mathbb{E}\left[\frac{1}{n}\sum _{i=1}^n\sum _{x=1}^2\sum _{y=1}^2{\xi }_x^z{A}_{x,y}^z\mathbb{I}({\mathcal{E}}_i^{xy})\right]\\ & =\sum _{x=1}^2\sum _{y=1}^2{\xi }_x^z{A}_{x,y}^{z}P({\mathcal{E}}_i^{xy}),\end{array}$$
式中，$\mathbb{I}(\cdot )$为指示函数，$A_{x,y}^z$表示转移矩阵$A^z$第$x$行第$y$列的元素。

定理1：假设有两个二状态隐马尔可夫模型，满足${\mu }_1^i<{\mu }_2^i$、${\sigma }_1^i={\sigma }_2^i=\sigma$，则错误概率满足：
$$\begin{array}{rl}{P}_e\ge & \sum _{x=1}^2\sum _{y=1}^2\{{\xi }_x^z{A}_{x,y}^z\times \\ & \sum _{k=1}^2\sum _{l=1}^2{\xi }_k^y{\xi }_l^{\backslash y}\mathcal{Q}\left(\frac{(-1{)}^{u({\mu }_k^y-{\mu }_l^{\backslash y})}({\gamma }_{kl}^{xy}-{\mu }_l^{\backslash y})}{\sigma }\right)\},\end{array}\text{\hspace{1em}(26)}$$
式中，$\mathcal{Q}(\cdot )$为标准高斯分布的右尾（生存）函数，$u(\cdot )$为单位阶跃函数，${\gamma }_{kl}^{xy}$定义为：
$${\gamma }_{kl}^{xy}=\frac{({\mu }_k^y{)}^2-({\mu }_l^{\backslash y}{)}^2}{2({\mu }_k^y-{\mu }_l^{\backslash y})}+2{\sigma }^2\frac{\log \frac{{\xi }_k^y{A}_{x,y}^z}{{\xi }_k^{\backslash y}{A}_{x,y}^z}+{\Xi }^{y\top }\log A_{\cdot ,k}^y-{\Xi }^{\backslash y\top }\log A_{\cdot ,l}^{\backslash y}}{2({\mu }_k^y{\xi }_k^y-{\mu }_l^{\backslash y}{\xi }_l^{\backslash y})},\text{\hspace{1em}(27)}$$
其中，$\backslash y=2$（若$y=1$），$\backslash y=1$（若$y=2$），$A_{\cdot ,l}^m$表示转移矩阵$A^m$的第$l$列。

证明：证明过程见附录D。

4.2 性能对比

本节对比四种所提算法的错误概率与定理1定义的下界。具体而言，考虑两个隐马尔可夫模型，其转移矩阵均为$A^1=A^2=\left[\begin{array}{cc}0.1& 0.9\\ 0.9& 0.1\end{array}\right]$；每个字母${\phi }_i^m$的发射函数服从$\mathcal{N}({\mu }_i^m,{\sigma }^2)$分布。设置三种不同场景，分别采用不同的均值${\mu }_i^m$：$[{\mu }_1^1,{\mu }_2^1,{\mu }_1^2,{\mu }_2^2]\in \{[1,2,3,4],[1,3,2,4],[1,3,4,2]\}$[3]；切换过程$P_{sw}$的转移矩阵设为$A^z=\left[\begin{array}{cc}0.1& 0.9\\ 0.9& 0.1\end{array}\right]$；取前$T=900$个交织数据作为混合隐马尔可夫模型的观测数据。解交织错误概率定义为：
$$P_e=\frac{{\sum }_{i=1}^T\mathbb{I}({\hat{d}}_i\ne d_i)}{T},\text{\hspace{1em}(28)}$$
式中，${\hat{d}}_i$为第$i$个观测的估计索引，$d_i$为第$i$个观测的真实索引。将每个字母的标准差从0调整至0.5（步长0.1），计算100次实验的平均错误概率，并与下界对比，仿真结果如图7所示。

图7 三种场景下的错误概率，结果与理论下界（式26，红色线）对比：场景1、场景2、场景3分别对应$[{\mu }_1^1,{\mu }_2^1,{\mu }_1^2,{\mu }_2^2]\in \{[1,2,3,4],[1,3,2,4],[1,3,4,2]\}$

如图7所示，四种所提算法均表现出有效的性能。需注意的是，随着标准差趋近于0，PGAS算法的错误概率可能无法趋近于0，且其曲线呈现非单调趋势。PGAS性能较差的主要原因是，基于采样的方法得到的后验轨迹会将不确定性引入采样轨迹中。相比之下，精确EM方法和两种变分方法的错误概率均随标准差减小而趋近于0。从以下两个角度分析这一现象：

首先，分析精确方法与变分方法的关系。随着标准差趋近于0，精确方法与变分方法的性能逐渐等效，这一特性在文献[30]中被称为“渐近有效性”。该特性表明，尽管变分证据下界与对数似然函数之间可能存在差距，但所提变分方法在渐近意义下会收敛到精确解。本文通过仿真验证了变分方法的渐近有效性，其理论证明超出本文范围，将在未来研究中展开。目前关于变分方法渐近有效性和渐近一致性的研究较少，但已有部分相关成果：文献[30]讨论了线性状态空间系统变分方法的渐近有效性和渐近一致性缺失问题（其中“一致性”指当样本量趋近于无穷大时，精确EM方法与变分方法等效）；文献[31]针对特定马尔可夫模型，证明了通过最大化变分证据下界得到的参数估计具有渐近一致性。

其次，讨论理论下界与性能曲线的关系。随着标准差增大，算法的错误概率和理论下界均随之增大，且理论下界与错误概率曲线之间的差距扩大。产生这一差距的主要原因是，多种错误事件都可能导致解交织错误，但本文分析仅聚焦于特定错误事件${\mathcal{E}}_i^{xy}$（即第$i$个脉冲识别错误而第$i-1$个脉冲识别正确）。随着标准差增大，多个连续脉冲被错误识别的概率升高，进一步扩大了差距。未来研究可从两个方向缩小下界与仿真曲线的差距：一是分析更广泛的解交织错误事件；二是融合到达时间等额外信息，改进算法性能分析，助力缩小差距。

尽管如此，所提变分方法仍表现出较低的错误概率，因此有理由认为，在实际应用中，变分方法可作为精确方法的有效替代方案。

五、仿真

通过对仿真混合脉冲序列开展综合实验，验证所提算法的有效性与优越性。5.1节介绍数据集、评估指标及当前主流的基准方法；5.2节呈现并讨论实验结果。

5.1 实验设置

5.1.1 数据集描述

本研究生成由3部雷达发射的混合脉冲序列。尽管雷达信号解交织可利用脉冲宽度（PW）、脉冲幅度（PA）等多个PDW参数，但本文遵循文献[32]的设置，仅聚焦于射频（RF）信息。在此设定下，字母表中的每个字母均建模为一元高斯随机变量。生成5类数据集以评估所提算法的有效性，数据集详情如表1所示。

表1 5类数据集的描述

雷达索引	射频调制类型	射频值（MHz）	PRI调制类型	PRI值（μs）
实验类别I（不交子字母表）	-	-	-	-
1	滑动	1025, 1050, 1075, 1100, 1125	参差	70, 110
2	参差	900, 1000, 960, 980, 1120	参差	60, 100
3	驻留切换（D&S）	910, 920, 930, 940, 950	参差	50, 120
实验类别II（非不交子字母表）	-	-	-	-
1	滑动	1025, 1050, 1075, 1100, 1125	参差	70, 110
2	参差	910, 1000, 960, 980, 1125	参差	60, 100
3	驻留切换（D&S）	910, 920, 930, 940, 950	参差	50, 120
实验类别III（时间性能）	-	-	-	-
1	滑动	[1000, 1100]，共$K$个点	参差	70, 110
2	参差	[1100, 1200]，共$K$个点	参差	60, 100
3	驻留切换（D&S）	[1200, 1300]，共$K$个点	参差	50, 120
实验类别IV（相同雷达）	-	-	-	-
1	参差	1025, 1050, 1075, 1100, 1125	抖动	50
2	参差	1025, 1050, 1075, 1100, 1125	抖动	50
实验类别V（缺失脉冲）	-	对实验类别I和II的脉冲按均匀分布随机丢弃，构建缺失脉冲场景	-	-

• 实验类别I：聚焦于单部雷达采用多频率时的雷达脉冲序列解交织，该场景下雷达无共享频率（即子字母表不交），脉冲缺失概率设为0.05。
• 实验类别II：聚焦于雷达采用多频率且允许频率共享的脉冲序列解交织，此时交织隐马尔可夫过程的子字母表非不交。具体而言，第2部雷达与第1部雷达共享1个频率、与第3部雷达共享1个频率（共享频率以粗体标注），脉冲缺失概率同样设为0.05。
• 实验类别III：评估算法处理时间。频率从预定义范围中均匀采样，每部雷达选择的频率数量对应$K^m$（表1中方括号定义频率范围）；为简化计算，假设$K^1=K^2=K^3=K$，脉冲缺失概率设为0.05。
• 实验类别IV：针对两部“相同雷达”的脉冲解交织。“相同”指两部雷达的射频调制类型及数值、PRI调制类型及数值完全一致，仅初始射频值与初始PRI值不同；实验中，射频序列的初始相位从均匀分布中采样，脉冲缺失概率设为0.05。
• 实验类别V：评估算法对缺失脉冲的鲁棒性。脉冲缺失概率从0增至0.5（步长0.1）；需说明的是，本节中驻留切换（D&S）调制的驻留点数量设为20。

5.1.2 评估指标

采用解交织准确率评估算法性能，计算公式为：
$$P_a=\frac{{\sum }_{i=1}^T\mathbb{I}({\hat{d}}_i=d_i)}{T}\text{\hspace{1em}(29)}$$
其中，${\hat{d}}_i$为第$i$个观测的估计索引，$d_i$为第$i$个观测的真实索引；$\mathbb{I}(\cdot )$为指示函数。采用Munkres算法[33]将估计索引序列的随机索引映射到与真实索引序列重叠最大的索引集合，以消除索引随机性对准确率计算的影响。

5.1.3 基准方法

重新实现3种当前主流的解交织算法用于性能对比，具体如下：

1. IMP（交织马尔可夫过程）：文献[13]提出的解交织方案，将数据生成过程建模为交织马尔可夫过程$P=I_{\Pi }^I(P_1,\dots ,P_M;P_{sw})$，融合到达顺序（AO）与射频（RF）信息。由于该方法未估计单个脉冲的标签，无法直接适配本文实验，需对其改进（改进细节见附录E）。
2. TEDS（马尔可夫更新链）：文献[16]提出的解交织方案，将每个分量过程建模为马尔可夫更新链$(p[{\phi }^m],\tau [{\phi }^m])$（其中τ={τi}i=1T\tau=\{\tau_i\}_{i=1}^Tτ={τi​}i=1T​为字母的TOA序列，$\tau [{\phi }^m]$为子字母表${\phi }^m$对应字母的TOA序列），融合时序信息与PDW信息。同理，该方法需改进以适配实验（改进细节见附录E）。
3. OT（最优传输距离）：文献[7]提出的解交织算法，流程为：先对观测序列聚类形成字母表，再计算所有字母间的最优传输距离矩阵，最后基于该矩阵通过层次聚类得到解交织结果。实验中假设聚类可形成正确字母表，且最优传输距离基于每个字母的平均PRI值计算；由于OT仅适用于不交子字母表场景，因此仅在实验类别I中参与对比。

为保证对比公平性，IMP与TEDS算法通过遗传算法（GA）[34]根据改进的代价函数优化。借助改进的代价函数，IMP与TEDS可适配非不交子字母表场景。

5.2 性能验证

由于基于EM的算法仅能保证收敛到局部最优，因此良好的初始化至关重要。实验中，采用狄利克雷过程混合模型（DPMM）[35]提取字母并形成字母表，再根据IMP[13]搜索字母表划分结果。基于该划分，验证所提方法与基准方法的性能。各方法的参数设置如表2所示；在相同雷达解交织仿真（实验类别IV）中，假设字母表划分已知。对每个数据集计算解交织准确率，并在100次蒙特卡洛仿真中取平均值。

表2 各方法的参数设置

方法	参数	描述	数值
EM	$\Pi$	字母表划分	由IMP初始化
	$ϵ$	收敛阈值	$10^{-5}$
PGAS	$\Pi$	字母表划分	由IMP初始化
	$N$	粒子数量	1000
	$ϵ$	收敛阈值	$10^{-5}$
SVI	$\Pi$	字母表划分	由IMP初始化
	$ϵ$	收敛阈值	$10^{-5}$
MFVI	$\Pi$	字母表划分	由IMP初始化
	$ϵ$	收敛阈值	$10^{-5}$
IMP	$\Pi$	字母表划分	由IMP初始化
	$\beta$	惩罚参数	1
	$pop$	遗传算法种群数量	3000
	$cr$	遗传算法交叉率	0.99
	$sr$	遗传算法选择率	0.05
	$vr$	遗传算法变异率	0.99
	$ϵ$	收敛阈值	$10^{-5}$
TEDS	$\Pi$	字母表划分	由IMP初始化
	$\beta$	惩罚参数	1
	$pop$	遗传算法种群数量	3000
	$cr$	遗传算法交叉率	0.99
	$sr$	遗传算法选择率	0.05
	$vr$	遗传算法变异率	0.99
	$ϵ$	收敛阈值	$10^{-5}$
OT	$\vartheta$	阈值	0.01

5.2.1 解交织准确率

实验类别I（不交子字母表）：评估不同射频标准差下的解交织准确率，标准差从0增至4（步长1），结果如图8（a）所示。可观察到：当标准差趋近于0时，除PGAS外，所有算法的准确率均趋近于1；PGAS性能较差的原因是其基于采样的后验轨迹会将不确定性引入结果。随着标准差增大，所有算法的准确率均下降；其中TEDS性能最佳，其次是IMP。

实验类别II（非不交子字母表）：评估不同射频标准差下的解交织准确率，结果如图8（b）所示。与实验类别I一致，TEDS与IMP仍表现最优，PGAS与MFVI性能相对较差；但值得注意的是，SVI性能超过EM算法——这是因为EM算法需在大状态空间执行E步，而在非不交子字母表场景中，过多隐状态输出相同字母会导致参数估计误差增大（参数输出隐马尔可夫模型中，估计误差与隐状态数量成正比[29]），最终使EM性能下降。

图8 不交与非不交子字母表下的解交织准确率：（a）实验类别I；（b）实验类别II

进一步分析所提4种算法的性能差异：

• EM算法在不交子字母表中准确率最高，但在非不交场景中受大状态空间影响性能下降；
• SVI在两种场景中均保持较高准确率，且在非不交场景中优于EM，体现出对复杂状态空间的适配性；
• MFVI因“隐变量独立”假设过于严格，准确率始终低于SVI；
• PGAS因采样引入不确定性，准确率在所有场景中均最低。

基准方法中，OT性能最差——其采用PRI值计算传输距离，而参差PRI调制会破坏OT对“字母相似性”的假设（原始OT研究[7]采用脉冲幅度与TOA信息计算距离，更适配字母相似场景），导致性能下降，因此后续实验不再考虑OT。

5.2.2 时间性能对比

为评估各方法的计算效率，在实验类别III中改变$K$（每个雷达的频率数量），射频标准差设为3 MHz，在搭载Intel Core i7-13700K CPU的设备上通过MATLAB实现实验，记录处理时间与解交织准确率，结果如表3所示。

表3 所提方法与基准方法的时间性能及解交织准确率（最佳值用粗体标注）

方法	$K=5$（时间/s，准确率）	$K=10$（时间/s，准确率）	$K=15$（时间/s，准确率）
EM	1.7327，0.9723	2.7853，0.9899	6.6417，0.9969
PGAS	5.4185，0.9754	38.8000，0.9909	365.7490，0.9964
SVI	0.0446，0.9969	0.0617，0.9969	0.0680，0.9969
MFVI	0.0396，0.9754	0.0634，0.9959	0.1170，0.9989
IMP	7.9785，0.9909	17.3492，0.9949	18.0041，0.9989
TEDS	8.0974，1.0000	18.8247，1.0000	253.5000，0.6886

由表3可得出以下结论：

1. SVI效率最优：单次迭代时间复杂度为$O\left(T{\sum }_{m=1}^M(K^m{)}^2+T{M}^2\right)$，且不受$K$增大的显著影响，始终保持最快处理速度；
2. EM效率随$K$增长显著下降：由于EM需处理$M{\prod }_{m=1}^M{K}^m$个隐状态，$K$增大导致状态数量呈指数增长，处理时间急剧增加；
3. PGAS、IMP、TEDS效率较低：PGAS需计算每个粒子的重要性权重，IMP与TEDS需通过遗传算法计算种群适应度，均导致计算负担加重，不适用于实时处理场景；
4. MFVI效率接近SVI，但准确率较低：其“隐变量独立”假设简化了计算，但牺牲了模型表示能力，导致准确率低于SVI。

综合准确率与时间性能，SVI实现了最优平衡：相比IMP、TEDS，SVI处理速度提升两个数量级，且准确率仅略有下降；相比EM、PGAS、MFVI，SVI在准确率与效率上均占优。

5.2.3 对相同雷达的鲁棒性

相同雷达的解交织是公认的难题（两部雷达的PDW参数与调制模式完全一致，仅初始相位不同）。实验类别IV中，分别测试“初始相位差（IPD）变化”与“PRI标准差变化”对解交织准确率的影响，结果如图9所示。

图9 两部相同雷达的解交织准确率：（a）改变初始相位差；（b）固定初始相位差（20 μs），改变PRI标准差

• 初始相位差的影响（图9（a））：当IPD完全重叠（IPD=0,50,100,…）时，无法通过PDW信息区分脉冲，所有算法准确率均较低；当IPD不完全重叠时，SVI与IMP准确率迅速提升至近100%，而TEDS性能逐步提升——这是因为TEDS依赖TOA信息，而TOA对相位差的敏感度低于AO（SVI与IMP采用AO信息）。
• PRI标准差的影响（图9（b））：当PRI标准差小于8时，SVI与IMP性能相近；当标准差超过8时，SVI性能略有下降，但仍优于TEDS；TEDS受PRI偏差影响最大，因为其基于马尔可夫更新链建模TOA序列，PRI偏差会直接破坏时序相关性。

该实验表明：采用AO信息的方法（SVI、IMP）对相同雷达的鲁棒性优于依赖TOA的方法（TEDS），尤其在PRI存在偏差的场景中优势更显著。

5.2.4 对缺失脉冲的鲁棒性

实验类别V中，测试脉冲缺失概率（0~0.5，步长0.1）对解交织准确率的影响，射频标准差设为1 MHz，结果如图10所示。

图10 不交与非不交子字母表下缺失脉冲对解交织准确率的影响：（a）实验类别I；（b）实验类别II

由图10可见：

• IMP与TEDS性能随缺失概率增大显著下降：两者均基于“最小化熵值”的准则，脉冲缺失会破坏观测序列的完整性，导致熵值计算偏差增大；
• SVI性能保持稳健：其通过变分分布近似后验，能够有效处理观测不确定性，对脉冲缺失的容忍度更高；即使缺失概率达到0.5，SVI的准确率仍保持在0.85以上，显著优于其他方法。

六、结论

本文提出一种生成模型用于表示交织隐马尔可夫过程（IHMP），并将隐马尔可夫过程（HMP）混合解交织问题转化为后验推断问题。为计算后验概率，设计了1种精确推断算法（EM）与3种近似推断算法（PGAS、MFVI、SVI）；同时，基于似然比检验推导了解交织错误概率的理论下界，用于评估算法性能。

在多类雷达数据集上的仿真验证了以下核心结论：

1. SVI性能最优：在解交织准确率与计算效率之间实现了最佳平衡，可作为精确推断（EM）的有效替代方案；
2. AO信息更鲁棒：依赖AO信息的方法（SVI、IMP）对PRI偏差与相同雷达场景的鲁棒性优于依赖TOA的方法（TEDS）；
3. SVI抗缺失能力强：相比IMP、TEDS，SVI对脉冲缺失的容忍度更高，在缺失概率0.5时仍保持高准确率。

未来研究将聚焦两个方向：

1. 开发和积算法（SPA）[40,41]以降低精确推断的计算负担，该算法在多分量关联任务中已被证明有效；
2. 探索模型基数自动确定方法[42]，实现雷达数量与字母表划分的自适应学习，进一步提升算法的实用性。