Deinterleaving of Mixtures of Renewal Processes

Jeremy Young
TSP

摘要

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

索引词

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

I. 引言

我们研究如下问题：源k∈{1,…,K}k in {1, …, K}k∈{1,…,K}在时刻τk,i∈[0,T] au_{k,i} in [0, T]τk,i∈[0,T]产生事件（事件可为数据包传输或脉冲信号发射），观测者获得所有脉冲集合{τk,i}{ au_{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. 系统模型与问题描述

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

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

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

脉冲间隔τk,i−τk,i−1>0 au_{k,i} – au_{k,i-1} > 0τk,i−τk,i−1>0服从参数化先验分布Tk(t)T_k(t)Tk(t)，其参数可能已知或未知；对于k≠lk
eq lk=l或i≠ji
eq ji=j，脉冲间隔τk,i−τk,i−1 au_{k,i} – au_{k,i-1}τk,i−τk,i−1与τl,j−τl,j−1 au_{l,j} – au_{l,j-1}τl,j−τl,j−1相互独立。

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

在许多应用中，先验分布Tk(t)T_k(t)Tk(t)的类型未知，若无进一步信息，截断高斯分布是合理假设。从中心极限定理或最大熵角度可论证截断高斯分布的合理性：给定均值μmuμ和方差σ2sigma^2σ2，具有最大熵的正分布正是截断高斯分布[28]。为满足脉冲间隔非负约束，需在0处截断分布。定义均值为μmuμ、方差为σ2sigma^2σ2的0截断高斯分布为N0(μ,σ2)N_0(mu, sigma^2)N0(μ,σ2)，其概率密度函数为：

通过错误概率评估算法性能，即脉冲被错误分类的概率。设Di∈{1,…,K}D_i in {1, …, K}Di∈{1,…,K}表示时刻τi au_iτi脉冲的真实源，D^ihat{D}_iD^i为估计值，定义源kkk的错误概率为：

III. 已知参数的时序分离

基础方法中，假设源数量和先验分布完全已知，问题转化为利用已知先验分布Tk(t)T_k(t)Tk(t)，仅基于脉冲时刻{τi}{ au_i}{τi}对脉冲分类。考虑最大似然（ML）解：

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

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

若D101=1D_{101} = 1D101=1：log⁡T1(τ101−τ98)log mathcal{T}_1( au_{101} – au_{98})logT1(τ101−τ98)若D101=2D_{101} = 2D101=2：log⁡T2(τ101−τ100)log mathcal{T}_2( au_{101} – au_{100})logT2(τ101−τ100)若D101=3D_{101} = 3D101=3：log⁡T3(τ101−τ94)log mathcal{T}_3( au_{101} – au_{94})logT3(τ101−τ94)

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

若D101=1D_{101} = 1D101=1：log⁡T1(τ101−τ98)log mathcal{T}_1( au_{101} – au_{98})logT1(τ101−τ98)若D101=2D_{101} = 2D101=2：log⁡T2(τ101−τ100)log mathcal{T}_2( au_{101} – au_{100})logT2(τ101−τ100)若D101=3D_{101} = 3D101=3：log⁡T3(τ101−τ93)log mathcal{T}_3( au_{101} – au_{93})logT3(τ101−τ93)

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

对于大MMM，复杂度仍过高，因此考虑简化算法。首先注意，许多最优解路径的最终成为最优解的可能性极低。对路径PiP_iPi，除dk(Pi)d_k(P_i)dk(Pi)外，定义rk(Pi)r_k(P_i)rk(Pi)为距上一个分配给源kkk的脉冲的实际时间。设TkT_kTk满足：

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

输入：KKK个源、TkT_kTk、观测τi au_iτi（i=1,…,Mi=1,…,Mi=1,…,M）

初始化τ1 au_1τ1的分配：前KKK条路径中τ1 au_1τ1分别标记为1,…,K1,…,K1,…,K，设i=1i=1i=1。对每条路径，考虑τi+1 au_{i+1}τi+1分配给源k=1,…,Kk=1,…,Kk=1,…,K的扩展：
a. 对每个扩展路径Pi+1P_{i+1}Pi+1，更新(d1(Pi),…,dK(Pi))(d_1(P_i), …, d_K(P_i))(d1(Pi),…,dK(Pi))：若扩展为D^i+1=khat{D}_{i+1} = kD^i+1=k，则

初始化τ1 au_1τ1的分配：前KKK条路径中τ1 au_1τ1分别标记为1,…,K1,…,K1,…,K，剩余a−Ka-Ka−K条路径分配任意，设i=1i=1i=1。对每条路径，考虑τi+1 au_{i+1}τi+1分配给源k=1,…,Kk=1,…,Kk=1,…,K的扩展：
a. 对每个扩展路径Pi+1P_{i+1}Pi+1，更新(r1(Pi),…,rK(Pi))(r_1(P_i), …, r_K(P_i))(r1(Pi),…,rK(Pi))：若扩展为D^i+1=khat{D}_{i+1} = kD^i+1=k，则

进一步简化的算法是保留固定数量（或最大数量）的路径（记为aaa），每步仅维持似然度最大的路径，同时剔除rk(Pi)>Tkr_k(P_i) > T_krk(Pi)>Tk的路径（避免源超时错误），该算法如算法3（记为A3）所示。或者，不剔除路径，而是当最可能路径出现源超时时重启算法（记为算法4，A4）。A4通常比A3更快，因无需检查每条路径的超时或管理路径删除。实验表明，a=K3a=K^3a=K3对两种版本均适用（见第V-A节），复杂度为O(MK3)O(MK^3)O(MK3)。

算法4：（A4）已知分布下带固定状态数和STO重启的简化时序分离算法

输入：KKK个源、TkT_kTk、观测τi au_iτi（i=1,…,Mi=1,…,Mi=1,…,M）、内存路径数aaa

下界基于粗略分析，聚焦一种特定错误事件：考虑源1的第iii个脉冲，j(i)j(i)j(i)为源2的最近脉冲，AiA_iAi表示脉冲iii和j(i)j(i)j(i)为最近邻且分类交换（即iii分配给2，j(i)j(i)j(i)分配给1）的事件。若AiA_iAi发生，源1产生一个错误；由于iii和j(i)j(i)j(i)为最近邻，该对与其他最近邻对不重叠。忽略其他错误类型（因推导下界），源1的错误概率下界为：

P(Ai)P(A_i)P(Ai)仍难以计算，但可通过“精灵”辅助得到下界：精灵提供除iii和j(i)j(i)j(i)外所有脉冲的正确决策，并告知iii和j(i)j(i)j(i)来自不同源。参考图1的符号说明：设δ=τj(i)−τidelta = au_{j(i)} – au_iδ=τj(i)−τi，tkpt_{kp}tkp为iii和j(i)j(i)j(i)中点到源kkk最近前向脉冲的时间差，tkst_{ks}tks为到最近后向脉冲的时间差；t=(t1p,t2p,t1s,t2s)t = (t_{1p}, t_{2p}, t_{1s}, t_{2s})t=(t1p,t2p,t1s,t2s)，d∈{(12),(21)}d in {(12), (21)}d∈{(12),(21)}（d=(12)d=(12)d=(12)表示第一个脉冲来自源1，第二个来自源2）。最优决策为：

若时序过程为泊松过程，T1T_1T1和T2T_2T2为指数分布（即Ti(t)=λie−λitT_i(t) = lambda_i e^{-lambda_i t}Ti(t)=λie−λit），则T1(t1p−12δ)T1(t1s+12δ)=λ12e−λ1(t1p+t1s)=T1(t1p+12δ)T1(t1s−12δ)mathcal{T}_1(t_{1p}-frac{1}{2}delta)mathcal{T}_1(t_{1s}+frac{1}{2}delta) = lambda_1^2 e^{-lambda_1(t_{1p}+t_{1s})} = mathcal{T}_1(t_{1p}+frac{1}{2}delta)mathcal{T}_1(t_{1s}-frac{1}{2}delta)T1(t1p−21δ)T1(t1s+21δ)=λ12e−λ1(t1p+t1s)=T1(t1p+21δ)T1(t1s−21δ)，T2T_2T2同理。因此两种概率相等，最优决策为随机选择假设，得到如下命题：

命题1：泊松过程无法仅基于时序分离。

对于截断高斯分布，可计算错误概率（错误即错误决策的f(t∣d)f(t|d)f(t∣d)大于正确决策的f(t∣d)f(t|d)f(t∣d)），得到如下定理：

定理2：设Ti∼N0(μi,σi2)T_i sim N_0(mu_i, sigma_i^2)Ti∼N0(μi,σi2)，则源1的错误概率下界为：

B. 性能比较

本节比较第三节的算法及定理2的下界。生成两个源的大量脉冲（T1∼N0(1,σ2)T_1 sim N_0(1, sigma^2)T1∼N0(1,σ2)，T2∼N0(π,σ2)T_2 sim N_0(pi, sigma^2)T2∼N0(π,σ2)），混合后取前500个脉冲作为接收信号，应用A1、A2、A3、A4，计算有限时间的源错误概率：

所有提出的算法性能相近：下界与算法性能相差约5倍常数因子，但斜率相同。不确定差距源于下界较松还是精确ML算法A1非直接最优解。随着σsigmaσ增大，简化算法（A3、A4）比“最优”算法性能略有下降，因此对于最优算法复杂度过高的大问题，可用简化算法替代。测试中两种简化算法无显著差异。

图2：两个合成源（红色μ=1mu=1μ=1，蓝色μ=πmu=piμ=π，标准差相同（x轴））的错误概率，A1（圆圈）、A2（方块）、A3（×）、A4（圆点）的结果为500次试验的平均值，虚线为理论下界（6）。

IV. 未知参数

前文假设分布TkT_kTk已知，现扩展至源数量KKK已知但TkT_kTk由未知参数向量θk heta_kθk（如θk=[μk,σk2] heta_k = [mu_k, sigma_k^2]θk=[μk,σk2]，μkmu_kμk和σk2sigma_k^2σk2为分布的均值和方差）定义的场景，需估计Θ=[θ1,…,θK]Theta = [ heta_1, …, heta_K]Θ=[θ1,…,θK]。不寻求单一“最优”ΘThetaΘ，而是生成多个ΘThetaΘ估计，对每个估计运行脉冲分离算法，选择使式(2)最大的分配作为最终输出。

可利用现有算法生成ΘThetaΘ估计（如[3], [4], [32], [33]），本文以delta-T直方图[3]为例：对交织脉冲序列的脉冲时刻，计算所有可能成对的脉冲间隔，构建直方图，峰值对应每个源的实际脉冲间隔（及其倍数），从峰值提取ΘThetaΘ估计。

delta-T直方图等价于脉冲时间函数自相关在直方图区间的积分[3]。在数字系统中，对采样脉冲时间函数计算自相关，结果与先计算自相关再积分一致。设脉冲时刻集合{τi}i=1M{ au_i}_{i=1}^M{τi}i=1M，脉冲时间函数为：

f(t)f(t)f(t)的采样版本定义为：

通常r[s]r[s]r[s]需计算至s=1,2,…,⌊T/γ⌋s=1,2,…,lfloor T/gamma
floors=1,2,…,⌊T/γ⌋；若已知最大脉冲间隔TmaxT_{max}Tmax，则计算至对应索引即可（已知最小间隔时同理）。γgammaγ的选择至关重要，后文详细说明。

直方图峰值表示脉冲间隔的重复次数多，但次谐波（实际脉冲间隔的整数倍）也会形成峰值，因此不能仅选择最高峰值作为μkmu_kμk，而应选择接近预期直方图高度的峰值。直方图高度为总样本数乘以样本落入区间的概率：若区间sss有峰值，假设μk=sγmu_k = sgammaμk=sγ（中心），TkT_kTk为单峰（如高斯）且不重叠，则峰值的理论高度为：

选择超过Γ(s)Gamma(s)Γ(s)的峰值，再挑选最接近Γ(s)Gamma(s)Γ(s)的KKK个峰值（排除次谐波）。次谐波的高度远大于Γ(s)Gamma(s)Γ(s)，因此基于差值（r[s]−Γ[s]r[s] – Gamma[s]r[s]−Γ[s]）选择可排除多数次谐波。若μjmu_jμj为μkmu_kμk的整数倍，μjmu_jμj会增强μkmu_kμk次谐波的峰值，此时需额外选择与Γ(s)Gamma(s)Γ(s)差值最大的KKK个峰值，共2K2K2K个候选峰值，提取μkmu_kμk和σksigma_kσk（位置和宽度），考虑所有无放回组合作为潜在ΘThetaΘ；若超过Γ(s)Gamma(s)Γ(s)的峰值数<K<K<K，则考虑有放回组合。完整过程如算法5（记为A5）所示。

γgammaγ过大时，所有数据落入同一区间，峰值数为0或1；γgammaγ过小时，同一分布的数据分散到多个区间，形成多个小峰值（如图4）。文献建议将γgammaγ设为近似抖动（即σksigma_kσk）[3], [20]，实际中σksigma_kσk未知，需根据观测调整：可尝试多个γgammaγ，对每个生成的ΘThetaΘ应用时序分离算法，选择使式(2)最大的结果（计算成本高）；或判断候选直方图，不满意时调整γgammaγ（如小γgammaγ下有多个小峰值，增大γgammaγ至超过Γ(s)Gamma(s)Γ(s)的峰值数接近∑k=1K⌊Tmaxμ^k⌋sum_{k=1}^K lfloor frac{T_{max}}{hat{mu}_k}
floor∑k=1K⌊μ^kTmax⌋，μ^khat{mu}_kμ^k为升序排列的峰值位置）。

delta-T直方图受次谐波峰值影响，但本文算法可拒绝次谐波，因此性能不受影响。其他算法（如CDIF、SDIF直方图、PRI变换）更直接处理次谐波，但专为周期性脉冲序列设计，对更新过程的抖动鲁棒性差；而delta-T直方图对抖动鲁棒（抖动仅加宽峰值，用于估计σksigma_kσk）。

算法5：时序参数初始化算法

输入：KKK个源的观测脉冲时刻{τi}{ au_i}{τi}、TmaxT_{max}Tmax、精度γgammaγ

构建脉冲时间序列f[n]f[n]f[n]：脉冲时刻为1，其余为0
a. 时间索引为⌊τ1/γ⌋lfloor au_1/gamma
floor⌊τ1/γ⌋到⌊τM/γ⌋lfloor au_M/gamma
floor⌊τM/γ⌋，每个点间隔γgammaγ。计算f[n]f[n]f[n]在0<s<⌊Tmaxγ⌋0 < s < lfloor frac{T_{max}}{gamma}
floor0<s<⌊γTmax⌋的自相关r[s]r[s]r[s]。找到r[s]r[s]r[s]的峰值及其宽度。对每个峰值，以其延迟作为μkmu_kμk、宽度作为σksigma_kσk，按式(7)计算Γ(s)Gamma(s)Γ(s)；若峰值被选中，将(μk,σk2)(mu_k, sigma_k^2)(μk,σk2)纳入ΘThetaΘ。比较峰值高度与Γ(s)Gamma(s)Γ(s)，考虑超过Γ(s)Gamma(s)Γ(s)的子集：
a. 若无峰值超过，调整γgammaγ重新开始：
i. 若有多个未超过的小峰值，增大γgammaγ；
ii. 若峰值极少（或无），减小γgammaγ。
b. 若峰值数<K<K<K，输出这些峰值的所有有放回组合。
c. 否则，计算r[s]−Γ[s]r[s] – Gamma[s]r[s]−Γ[s]，选择差值最小和最大的各KKK个峰值，输出这2K2K2K个峰值的所有无放回组合。

V. 仿真性能评估

通过软件生成符合指定截断高斯分布{Tk}k=1K{T_k}_{k=1}^K{Tk}k=1K的脉冲序列，混合后取MMM个连续脉冲作为待分离混合信号，用式(1)的有限MMM版本PeP_ePe作为性能指标，源编号任意，取所有标签排列的最小PeP_ePe。

A. 已知参数的算法

首先考虑两个合成源（均值μk∈[1,100]mu_k in [1,100]μk∈[1,100]，50个等间隔值；标准差σ=0.1sigma=0.1σ=0.1），生成100个脉冲，应用时序分离算法（使用真实参数），结果来自A2（已知参数下所有方法性能相近）。每个均值组合重复500次，计算平均PeP_ePe，结果如图5所示。

总体而言，除源均值小（μk<20mu_k < 20μk<20）且相近（μi≈μjmu_i approx mu_jμi≈μj）外，平均错误概率较低。原因：μi≈μjmu_i approx mu_jμi≈μj时，源频繁近乎同时发射脉冲，难以分离；μi=μjmu_i = mu_jμi=μj且σ≪μsigma ll muσ≪μ时，近乎同时发射会持续发生；小μkmu_kμk时，∣μk−σk∣vert mu_k – sigma_k vert∣μk−σk∣小，脉冲时刻易重叠；大∣μk−σk∣vert mu_k – sigma_k vert∣μk−σk∣可降低重叠概率。

实验表明，a=K3a=K^3a=K3对近似算法（A3、A4）足够。对a=K,K2,K3,K4,K5a=K, K^2, K^3, K^4, K^5a=K,K2,K3,K4,K5重复仿真，性能指标Pˉear{P}_ePˉe为(μ1,μ2)∈[1,100]2(mu_1, mu_2) in [1,100]^2(μ1,μ2)∈[1,100]2的平均值，结果如表I所示。路径数越多性能越好，但改进随KKK增大而减弱，K3K^3K3到K4K^4K4的改进<0.001%<0.001\%<0.001%，因此选择a=K3a=K^3a=K3。

表I 不同aaa值的错误率

aaa	KKK	K2K^2K2	K3K^3K3	K4K^4K4	K5K^5K5
Pe(×10−2)P_e( imes 10^{-2})Pe(×10−2)	0.83198	0.25730	0.25368	0.25283	0.25225
改进率(%)	–	0.5747	0.0036	0.0008	0.0006

B. 未知分布参数的算法

使用A5（γ=0.5gamma=0.5γ=0.5，Tmax=100T_{max}=100Tmax=100）生成候选ΘThetaΘ，两个合成源参数同上，错误率为100次试验的平均值，脉冲分离使用A2，比较参数估计与真实参数的错误率，结果如图6所示。

图6(a)为参数估计算法的错误率细节，图6(b)将(μ1,μ2)↦μ1−μ2(mu_1, mu_2) mapsto mu_1 – mu_2(μ1,μ2)↦μ1−μ2，沿相同索引平均以叠加比较。真实参数的A2平均错误率Pˉe=0.0025086ar{P}_e=0.0025086Pˉe=0.0025086，参数估计算法的Pˉe=0.011957ar{P}_e=0.011957Pˉe=0.011957，性能略有下降但仍较低。

分析图6(a)与图5（真实参数）的差异：μ1=μ2mu_1 = mu_2μ1=μ2时错误率较高（已知参数仅μk<20mu_k < 20μk<20时差，估计参数μk>20mu_k > 20μk>20时也差）；μ1<12μ2mu_1 < frac{1}{2}mu_2μ1<21μ2（μ1<40mu_1 < 40μ1<40）时错误率升高（源于次优ΘThetaΘ候选，μi≈2μjmu_i approx 2mu_jμi≈2μj时μimu_iμi峰值与μjmu_jμj次谐波叠加，候选选择错误）。

为验证算法适用于多源，考虑K∈{2,3,4,5}K in {2,3,4,5}K∈{2,3,4,5}（μkmu_kμk随机取自[1,100][1,100][1,100]，σ=0.1sigma=0.1σ=0.1），生成1000个混合脉冲，应用A4，比较ΘThetaΘ已知和估计的情况，每个KKK重复100次（每次生成新均值），平均PeP_ePe如表II所示。错误率随KKK增大而升高（错误类型增多），估计参数的错误率约为已知参数的10倍，但K=2K=2K=2时仍很低（0.18×10−20.18 imes 10^{-2}0.18×10−2，低于M=100M=100M=100的0.25×10−20.25 imes 10^{-2}0.25×10−2，更多数据降低错误）。

表II K∈[2,5]K in [2,5]K∈[2,5]、μk∈[1,100]mu_k in [1,100]μk∈[1,100]、σ=0.1sigma=0.1σ=0.1时A4的平均错误率(×10−2)( imes 10^{-2})(×10−2)

KKK	2	3	4	5
已知	0.182	0.457	0.758	3.098
估计	1.725	6.115	12.111	24.05

结论：简单参数估计方法在多数情况下性能良好，未来可改进参数估计算法（如第IV节末尾提及的方法），但本文重点为时序分离算法，不再深入参数估计。

VI. 应用

A. 神经元

提取单个神经元的活动通常需侵入式或专用设备，不切实际；而细胞外方法（如脑电图）可通过单个电极获取局部所有活动，识别单个峰电位（峰电位分选[9]）。假设峰电位间隔可用于分选，将两个细胞内记录的峰电位时刻混合（源1：θ=(0.17,0.022) heta=(0.17, 0.02^2)θ=(0.17,0.022)；源2：θ=(0.45,0.022) heta=(0.45, 0.02^2)θ=(0.45,0.022)），生成2秒16个脉冲的混合信号（图7上），应用A2（已知参数）得到完全正确的分配（图7下）。

图7：（上）基于记录动作电位生成的仿真混合信号；（下）A2的分配结果完全准确。

B. 抹香鲸

齿鲸通过点击序列（点击间隔ICI）回声定位和定向，ICI为微秒至2秒级，物种、功能、场景不同而变化，且序列内ICI慢变[1], [2]，适合时序分离。

为获得有真实标签的数据集，将单个抹香鲸记录的三个片段的时刻混合（参数：θ={(1.3940,0.12852),(0.8814,0.07082),(0.9212,0.06242)} heta={(1.3940, 0.1285^2), (0.8814, 0.0708^2), (0.9212, 0.0624^2)}θ={(1.3940,0.12852),(0.8814,0.07082),(0.9212,0.06242)}），生成40秒116个点击的混合信号（图8上），应用A2（已知参数），错误率8.6%（116个点击中10个错误，图8下黑色“×”标记），错误源于点击近乎同时到达。

图8：（上）基于记录抹香鲸点击生成的仿真混合信号；（下）A2的分配结果，10个错误用黑色“×”标记。

C. 雷达

雷达脉冲和时序为人为生成，仿真与真实数据等效。参考[7]，PRI为1μs级，抖动为PRI的8%，生成3个源（Tk∼N0(μk,σk2)T_k sim N_0(mu_k, sigma_k^2)Tk∼N0(μk,σk2)，Θ={(50,0.802),(72,0.402),(84,0.042)}Theta={(50, 0.80^2), (72, 0.40^2), (84, 0.04^2)}Θ={(50,0.802),(72,0.402),(84,0.042)}，单位μs和μs²），混合生成1000个脉冲，应用A3和A5（γ=1gamma=1γ=1，Tmax=100T_{max}=100Tmax=100）估计参数，错误率仅1.9%，估计参数为：

1) 比较

时序解交织的现有算法：[22]用聚类处理PRI抖动和脉冲缺失，[23][24]分别用方波正弦变换和傅里叶变换，[25]用模糊自适应谐振理论，核心均为“获取候选PRI→序贯搜索（SS）分配并移除序列”。SS为恒PRI设计，本文算法专为抖动PRI优化，比较已知参数下的SS与A2。

SS实现：脉冲到达时间TTT在目标时间x±εx pm varepsilonx±ε内视为匹配，ε=3σvarepsilon=3sigmaε=3σ（适应抖动），接受长度>3的序列。重复第III-B节仿真，σsigmaσ扩展至10010^0100，平均100次试验，结果如图9所示。A2全程优于SS，两者在高抖动下均失效。

考虑脉冲缺失情况，重复图9的实验（σk=0.1sigma_k=0.1σk=0.1，缺失率ςvarsigmaς从0到0.5，每次试验生成100个脉冲，随机移除ς×100varsigma imes 100ς×100个），结果如图10所示。即使未修改处理缺失脉冲，本文算法仍优于或等于SS（内置缺失处理）。

图9：SS算法与A2的错误率（μ1=1mu_1=1μ1=1，μ2=πmu_2=piμ2=π）随σksigma_kσk变化的100次试验平均值。

图10：不同缺失脉冲率下，SS算法与A3的错误率（μ1=1mu_1=1μ1=1，μ2=πmu_2=piμ2=π，σk=0.1sigma_k=0.1σk=0.1）的100次试验平均值。

VII. 结论

本文证明仅基于时序信息可实现脉冲源分离，开发了已知时序分布和源数量的脉冲分类算法，并扩展至未知参数场景。

提出多种算法，实验结论：已知参数时，可变路径数的精确ML算法（A1、A2）优于固定路径数的近似算法（A3、A4），但差异小，近似算法计算成本更低；未知参数时，改进现有参数估计算法可生成足够好的ΘThetaΘ，使分离算法收敛到良好解，识别了参数估计的失效场景，未来将改进。

纯时序分离在部分场景可行，实际数据中应结合脉冲形状（通常可用）提升性能，未来将研究这一方向。

附录定理2的证明

假设d=(12)d=(12)d=(12)，错误发生当（见式(5)）：

若Ti∼N0(μi,σi2)T_i sim N_0(mu_i, sigma_i^2)Ti∼N0(μi,σi2)，错误发生当：

等价于：

用脉冲间隔r1p=t1p−12δr_{1p}=t_{1p}-frac{1}{2}deltar1p=t1p−21δ、r1s=t1s+12δr_{1s}=t_{1s}+frac{1}{2}deltar1s=t1s+21δ等重写，得：

对于独立随机变量X,YX,YX,Y，X>a,Y>b⇒X+Y>a+bX>a, Y>b Rightarrow X+Y>a+bX>a,Y>b⇒X+Y>a+b，故P(X+Y>a+b)≥P(X>a)P(Y>b)P(X+Y>a+b) geq P(X>a)P(Y>b)P(X+Y>a+b)≥P(X>a)P(Y>b)，因此错误概率下界为：

由设定r1p>δr_{1p}>deltar1p>δ，随机变量1σ1(r1p−μ1)frac{1}{sigma_1}(r_{1p}-mu_1)σ11(r1p−μ1)的概率密度（t>−μ1+δσ1t>frac{-mu_1+delta}{sigma_1}t>σ1−μ1+δ）为：

r1s>2δr_{1s}>2deltar1s>2δ，随机变量1σ1(r1s−μ1)frac{1}{sigma_1}(r_{1s}-mu_1)σ11(r1s−μ1)的概率密度（t>−μ1+2δσ1t>frac{-mu_1+2delta}{sigma_1}t>σ1−μ1+2δ）为：

随机变量X1=−1σ1(r1p−μ1)+1σ1(r1s−μ1)X_1 = -frac{1}{sigma_1}(r_{1p}-mu_1) + frac{1}{sigma_1}(r_{1s}-mu_1)X1=−σ11(r1p−μ1)+σ11(r1s−μ1)的概率密度通过卷积积分计算，得：

源2的相关量同理。设gXi(x,t)g_{X_i}(x,t)gXi(x,t)为对应函数，D=δσ2D=frac{delta}{sigma_2}D=σ2δ，式(8)的第二个因子为∫D∞gX2(D,t)dtint_D^infty g_{X_2}(D,t)dt∫D∞gX2(D,t)dt，因此：

由更新理论[26, 定理10.3.5]，D=δσ2D=frac{delta}{sigma_2}D=σ2δ的分布为：

代入得错误概率下界为式(6)。