Deinterleaving of Mixtures of Renewal Processes

Deinterleaving of Mixtures of Renewal Processes

Jeremy Young
TSP

摘要

本文提出一种基于纯脉冲间隔的多输入单输出系统脉冲源分离（解交织）新方法。每个源的脉冲间隔建模为更新过程，为实现源分类，开发了基于脉冲间隔分布最大化脉冲分配似然度的算法。推导了性能的理论下界，且该算法被证明能接近此下界。将算法扩展至脉冲间隔分布参数未知的场景，并在多个应用实例中评估了性能。

索引词

解交织、点击序列、抖动脉冲重复间隔雷达、源分离、更新过程、海洋哺乳动物、脉冲源、多输入单输出

I. 引言

我们研究如下问题：源k∈{ 1,...,K}k \in \{1, ..., K\}k∈{ 1,...,K}在时刻${\tau }_{k,i}\in [0,T]$产生事件（事件可为数据包传输或脉冲信号发射），观测者获得所有脉冲集合{ τk,i}\{\tau_{k,i}\}{ τk,i​}但无任何标签，能否仅通过时序信息确定每个事件时刻所属的源？

该问题由多种应用场景驱动，包括：

• 计算机网络：数据包的时序能否揭示发送者身份？
• 海洋哺乳动物：能否仅基于时序分离海洋哺乳动物的回声定位点击序列[1], [2]？
• 雷达和声纳：能否仅利用脉冲到达时间被动分离（解交织[3]-[7]）多个脉冲雷达信号？
• 心跳：在混合心跳信号中（如胎儿心电图[8]），能否仅通过时序区分母亲和胎儿（多胎妊娠）的心跳？
• 神经系统：在神经峰电位序列记录中，能否分离不同源的峰电位序列[9]？

在许多应用中，除时序外还有脉冲幅度、形状等特征可辅助分离，但深入理解纯时序分离的基础问题仍具价值。例如，隐私/安全限制或传输带宽约束可能导致仅能获取时序信息；在次要特征未知的情况下，纯时序分离可作为初始分析手段；即使次要特征已知，本文方法遵循经典思路（如信息论[10]）：由实际问题构建简单理论/数学模型，深入研究后扩展至实际应用。未来可将本文纯时序分离算法扩展至融合次要特征，以适应真实场景。

脉冲序列分离并非新问题。例如，已有文献提出多种单传感器记录中海洋哺乳动物点击序列的分离方法：[11]利用点击互相关建立点击间关系并分离序列；[12]基于点击相关性、脉冲间隔细节和频率特征分离抹香鲸点击序列（最终依赖双水听器时延估计）；[13]对点击提取特征进行聚类分析；[14]逐次匹配连续点击的频谱；[15]基于频谱差异性分离序列。这些方法依赖点击序列内脉冲形状的慢变特性，仅有少数方法利用时序信息：[16]通过跟踪点击幅度和点击间隔（ICI）的慢变化分离两个点击序列；Baggenstoss的算法[17]-[19]基于包括ICI在内的点击对特征定义统计相似度度量，最大化关联点击的全局相似度；[20]与本文类似，仅基于时序分析点击序列，利用[21]的脉冲重复间隔（PRI）变换估计ICI并扩展至时变ICI，但未对每个脉冲进行分类。

仅基于脉冲到达时间的雷达信号解交织也已有研究（如[4]-[7], [22]-[25]）。雷达信号通常为周期性、交错式或抖动式，抖动PRI雷达中发射机为每个PRI添加随机时间，使其成为更新过程（见第二节），因此本文更新过程的结果可直接应用于抖动PRI雷达信号解交织。[6]考虑带噪声的到达时间，[7]采用PRI随机游走模型，[21]专门研究抖动PRI但仅用于估计PRI，[22]-[25]开发抖动PRI解交织的启发式方法。本文算法为该问题提供最优解，并在第VI-C1节与[22]-[25]进行比较。

本文专门开发基于纯时序的更新过程分离（解交织）算法，我们认为这是一个新问题。

---

II. 系统模型与问题描述

考虑包含$K$个源的系统模型，每个源产生事件或发射类脉冲信号。源$k$在时刻${\tau }_{k,i}$发射脉冲序列，假设脉冲响应的支撑集远小于脉冲间隔，因此无重叠（精确来说，重叠概率足够低可忽略）。不失一般性，假设第一个脉冲发生在$t=0$时刻。

我们需确定每个脉冲所属的源，目标仅基于未分类的脉冲时刻{ τi}\{\tau_i\}{ τi​}进行分类，假设脉冲已通过低误差概率的检测算法提取。

为仅基于时序信息分类{ τi}\{\tau_i\}{ τi​}，做出如下假设：

• 脉冲间隔${\tau }_{k,i}-{\tau }_{k,i-1}>0$服从参数化先验分布$T_k(t)$，其参数可能已知或未知；
• 对于$k\ne l$或$i\ne j$，脉冲间隔${\tau }_{k,i}-{\tau }_{k,i-1}$与${\tau }_{l,j}-{\tau }_{l,j-1}$相互独立。

基于这些假设，单个源的事件时刻集合构成更新过程或计数过程[26, 第8.3节]，最著名的例子是泊松过程。更新过程广泛用于建模多种过程[26, 第8.3节]，尤其适用于生物信号建模。从推测角度，生物系统可能没有内部节拍器控制事件（如心跳、点击）产生，而是通过计时器倒计时至下一个事件；为生成近周期性信号，时序分布可在均值附近高度集中。时序分布可能非独立，但一阶近似通常假设独立，从而形成更新过程。虽然点击序列的最优模型未知，但心跳常以此方式建模[27]。

在许多应用中，先验分布$T_k(t)$的类型未知，若无进一步信息，截断高斯分布是合理假设。从中心极限定理或最大熵角度可论证截断高斯分布的合理性：给定均值$\mu$和方差${\sigma }^2$，具有最大熵的正分布正是截断高斯分布[28]。为满足脉冲间隔非负约束，需在0处截断分布。定义均值为$\mu$、方差为${\sigma }^2$的0截断高斯分布为$N_0(\mu ,{\sigma }^2)$，其概率密度函数为：
$$f(x)=\frac{1}{\Phi (\frac{\mu }{\sigma })\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right)(x>0)$$
其中$\Phi$为标准正态累积分布函数[26]。在多数情况下，$\mu \gg 0$且$\mu \gg \sigma$，因此截断修正项$\Phi (\frac{\mu }{\sigma })$接近1，在许多计算中可忽略。

通过错误概率评估算法性能，即脉冲被错误分类的概率。设Di∈{ 1,...,K}D_i \in \{1, ..., K\}Di​∈{ 1,...,K}表示时刻${\tau }_i$脉冲的真实源，${\hat{D}}_i$为估计值，定义源$k$的错误概率为：
$$P_{ϵ,k}=\underset{T\to \infty }{\lim }P({\hat{D}}_j\ne D_j|D_j=k)$$
总错误概率为：
$$P_e=\underset{T\to \infty }{\lim }P({\hat{D}}_j\ne D_j)\text{\hspace{1em}(1)}$$

---

III. 已知参数的时序分离

基础方法中，假设源数量和先验分布完全已知，问题转化为利用已知先验分布$T_k(t)$，仅基于脉冲时刻{ τi}\{\tau_i\}{ τi​}对脉冲分类。考虑最大似然（ML）解：
$$\underset{D}{\max }L;L=\sum _{k=1}^K\sum _i\log {\mathcal{T}}_k\left({\tau }_{S_k(i)}-{\tau }_{S_k(i-1)}\right)$$
其中$S_k(i)$表示分配给第$k$个源的第$i$个脉冲，是$\hat{D}=(D_1,D_2,...)$的函数。例如，若两个源的$\hat{D}=(1,1,2,1,2,2)$，则$S_1=(1,2,4)$，$S_2=(3,5,6)$。

暴力最大化式(2)的复杂度为$O(K^M)$（$M$为接收信号中的总脉冲数）；借鉴维特比算法[30], [31]的思想，可得到更低复杂度的算法。为开发这些算法，首先以适应当前场景的方式描述维特比算法：假设需在时刻$1,2,3,...$从$K$个选项中决策，且过程为马尔可夫过程。时刻$i$的路径为决策集合Pi={ D^j;j=1,...,i}P_i = \{\hat{D}_j; j=1,...,i\}Pi​={ D^j​;j=1,...,i}，$P_{i,k}$表示指向${\hat{D}}_i=k$的路径集合。由于马尔可夫性，时刻$i+1,i+2,...$的最优决策与导致${\hat{D}}_i=k$的路径无关，因此只需记忆每个$k$对应的最优路径。每个阶段有$K$条最优路径，对每条指向${\hat{D}}_i=k$的最优路径，考虑扩展至时刻$i+1$的$K$种可能，再为每个${\hat{D}}_{i+1}=k$选择$K$条最优路径，复杂度为$O(MK)$。

对于时序分离，每个源的时序过程假设为马尔可夫过程（即更新过程），但混合过程通常不满足马尔可夫性，时刻$i+1,i+2,...$的决策并非仅依赖时刻$i$的决策。然而，仍可丢弃大量路径，具体如下：对路径$P_i$，定义$d_k(P_i)$为距上一个分配给源$k$的脉冲的脉冲数。例如，若$P_6=\hat{D}=(1,1,2,1,2,2)$，则$d_1(P_6)=2$，$d_2(P_6)=0$，$d_3(P_6)=\infty$。关键观察：与维特比算法类似，具有相同$K$元组$(d_1(P_i),d_2(P_i),...,d_K(P_i))$的路径对未来决策等价。通过示例可清晰理解：假设$P_{100}=(...,3,1,1,2,1,2,2)$，当脉冲101分配给三个源时，似然度分别添加如下项：

• 若$D_{101}=1$：$\log {\mathcal{T}}_1({\tau }_{101}-{\tau }_{98})$
• 若$D_{101}=2$：$\log {\mathcal{T}}_2({\tau }_{101}-{\tau }_{100})$
• 若$D_{101}=3$：$\log {\mathcal{T}}_3({\tau }_{101}-{\tau }_{94})$

显然，这些项与路径中“…”表示的3之前的部分无关；若$P_{100}=(...,3,2,1,1,2,1,2,2)$，则修正项为：

• 若$D_{101}=1$：$\log {\mathcal{T}}_1({\tau }_{101}-{\tau }_{98})$
• 若$D_{101}=2$：$\log {\mathcal{T}}_2({\tau }_{101}-{\tau }_{100})$
• 若$D_{101}=3$：$\log {\mathcal{T}}_3({\tau }_{101}-{\tau }_{93})$

因此，时刻101决策的似然度仅依赖路径的尾部，且由$K$元组$(d_1(P_i),...,d_K(P_i))$完全表征。因此，对于具有相同$K$元组的路径，只需记忆最优路径。根据定义，$d_k(P_i)$可取$i+1$个值（即$0,1,...,i-1,\infty$），且不同$k$的$d_k(P_i)$不同。因此，时刻$i$的独特路径数为$(i+1)i(i-1)\cdots (i+1-(K-1))\sim i^K$，总体复杂度（需记忆的路径数）为$O(M^{K+1})$，对于中等$M$和$K$可行。该ML算法如算法1（记为A1）所示。

对于大$M$，复杂度仍过高，因此考虑简化算法。首先注意，许多最优解路径的最终成为最优解的可能性极低。对路径$P_i$，除$d_k(P_i)$外，定义$r_k(P_i)$为距上一个分配给源$k$的脉冲的实际时间。设$T_k$满足：
$${\int }_{T_k}^{\infty }{\mathcal{T}}_k(\tau )d\tau <ϵ$$
其中$ϵ$为小值。若$r_k(P_i)>T_k$，则将路径$P_i$扩展为$j>i$时${\hat{D}}_j=k$，似然度至少添加$\log ϵ$项。若$ϵ$很小，该项为较大负值，使得路径$P$的扩展很难成为最大似然解，称为“源超时（STO）”。因此，可从保留路径集中剔除此类路径，与最优解相比损失极小（若$ϵ$选择合理）。该STO算法如算法2（记为A2）所示。

算法1：（A1）已知分布下仅基于时序分离的精确最大似然解

输入：$K$个源、$T_k$、观测${\tau }_i$（$i=1,...,M$）

1. 初始化${\tau }_1$的分配：前$K$条路径中${\tau }_1$分别标记为$1,...,K$，设$i=1$。
2. 对每条路径，考虑${\tau }_{i+1}$分配给源$k=1,...,K$的扩展：
   a. 对每个扩展路径$P_{i+1}$，更新$(d_1(P_i$