基于稀疏重构的混叠脉冲序列的周期估计

许成维，陶建武

空军航空大学 航空作战勤务学院，长春 130022

随着社会信息化的快速发展，周期估计广泛应用于电子侦察中脉冲重复间隔估计、节拍和节奏估计、神经学中峰值区间估计、呼吸率周期检测、蛋白质编码和脉冲星周期估计等各个军民应用领域。尤其在电子战领域的应用更为显著，多用于电子干扰与反干扰的脉冲重复间隔(PRI)估计[1-7]。PRI是现代电子战环境下雷达脉冲信号分选与识别的一个至关重要的参数。近年来，随着电子战技术的发展，电磁环境越发复杂、脉冲高度密集且相互交错，影响了现代电子侦察技术的精确识别，这对于现代电子侦察提出了新的挑战。当雷达发射脉冲信号时，为了防止对方轻易侦察到雷达发射的脉冲信号的重复间隔，往往会在发射脉冲时，加入人为干扰和抖动，导致了脉冲重复间隔的不确定性。由于接收设备自身原因和外界环境干扰因素，往往造成某些脉冲重复间隔上的脉冲缺失，导致测量脉冲串的不完整，使得测量脉冲串的周期性遭到破坏。尤其是在复杂电磁环境中，测量脉冲串通常是由多个具有不同周期的脉冲串混合叠加而成。这些因素都不利于信号周期的提取。

针对不完整脉冲串的周期估计，许多学者进行了大量的研究，如Fogel和Gavish[8]使用最大似然(ML)算法对带有缺失和抖动的数据进行周期估计。Sadler和Casey[9]基于修正的欧式算法对最大公约数进行计算。Sidiropoulos等[10]在此基础上，进一步发展了不完整脉冲串数据的周期估计方法。Clarkson[11]对最大似然估计做了改进，并验证了两种估计方法在周期估计上的性能。Ye等[12-13]指出了改进的Fogel周期图估计算法的性能是由空间网格的间隔决定，并从理论推导出合理的网格搜索间隔。McKilliam和Clarkson[14]通过使用Chirp z-变换和快速傅里叶变换，使其周期图估计算法的计算复杂性减少到O(NlgN)。虽然这些算法在周期估计的精确度和减小计算复杂性方面做出了很大的改进，同时对于低信噪比和高缺失率的信号具有较好的性能。但是对于由多个具有不同周期的脉冲串混合叠加的脉冲序列，这些方法不能得到真实的周期估计。因此，限制了这些方法在复杂电磁环境中的应用。面对这些挑战，需要我们提出全新的、可以应对各种复杂电磁环境的脉冲序列的周期估计方法。

近年来，随着压缩感知的提出与发展，其应用一直是关注热点。基于Farey序列理论，Vaidyanathan和Pal[15]构造了一种能够提取出信号隐藏周期的过完备字典。Tenneti和Vaidyana-than[16]将此字典用于心电图、蛋白质编码等真实信号的隐藏周期提取。但是，该字典是基于离散傅里叶变换(DFT)矩阵进行构造的，对于周期信号的表示不能形成一个完美的稀疏表示。基于Ramanujan结构字典(RD)理论，Vaidyanathan[17-18]提出了一种提取信号隐藏周期的方法。但是，在应用于带有人为抖动的、不完整的混叠脉冲序列的周期估计时，该方法不能达到较好的提取效果。因此，针对由多个具有不同周期的脉冲串混叠而成的脉冲序列，本文提出了一种基于稀疏重构的隐藏子周期估计新方法。从理论上探讨混叠脉冲序列的稀疏表示模型、稀疏字典的设计及隐藏子周期估计方法。通过仿真实验验证了此估计方法的可行性和有效性。

1 数据模型

假设周期脉冲信号的脉冲宽度远小于脉冲的周期。因此，可以将周期脉冲信号看成具有单位幅值的周期脉冲序列，从而形成周期脉冲信号的点过程模型。本文考虑如下周期脉冲信号的点过程模型：一个周期脉冲信号为T0k+θ0(k∈Z)是连续的整数，T0为周期，θ0为初相位。对其进行N次观测，其表达式为

tn=T0sn+θ0+ωnn=1,2,…,N

(1)

式中:ω1,ω2,…,ωN为随机抖动(或称为脉冲滑动误差)；s1,s2,…,sN为不连续且未知的整数，由于外界环境干扰，会有部分脉冲未观测到。这造成了接收到的脉冲信号不是连续和完整的，即s1,s2,…,sN是不连续的整数。图1给出了从完整的周期脉冲信号到形成不完整且带有脉冲滑动误差的单脉冲序列的示意图。

假设有M个具有不同周期的周期脉冲信号同时到达接收装置，因此，整个脉冲序列是M个不完整且带有脉冲滑动误差的脉冲串混合而成，其表达式为

(2)

图2给出了从两个完整的脉冲周期信号到形成不完整且带有脉冲滑动误差的混叠脉冲序列的示意图。图2中，假设第1个脉冲序列的周期T1=15，第2个脉冲序列的周期T2=18。混合后，整个脉冲序列的周期为2个子脉冲序列周期的最小公倍数。本文所感兴趣的是从整个脉冲序列中估计子脉冲序列的周期Tm。

2 结构字典矩阵设计

2.1 Farey结构字典矩阵

给定一个整数P，考虑所有不可约的有理数k/q,其中：整数q∈{1,2,…,P}，k为小于等于q、且与q互质的整数，即gcd(k,q)=1(1≤k≤q)。gcd(·)=1表示为两个数的最大公约数为1，即它们是互质的。表1列出了在周期P=6时，有理数k/q的所有数值。将k/q的所有数值按从小到大的顺序进行排列，得到一个有理数序列FP，即Farey序列。

表1 周期P=6时，有理数k/q的数值Table 1 Period P=6, rational values k/q

(3)

Vq=[v(q,k1)v(q,k2) …v(q,kφ(q))]

(4)

这里gcd(q,ki)=1和1≤ki≤q，Wq=e-j2π/q。进一步，定义一个P×Φ(P)维Farey序列字典矩阵为

(5)

例如，当P为6时，6×12维Farey序列结构字典矩阵为

(6)

利用Farey序列结构的字典矩阵，能够将具有不同周期的混叠脉冲序列进行稀疏化表示，然后通过一定的稀疏恢复算法，获得混叠脉冲序列的不同子周期。但是，字典矩阵维数的大小是影响稀疏恢复算法复杂性的关键因素，随着P的增大，Φ(P)将急速增大，字典矩阵的维数增大。所以我们应设计一个合适的插值尺度来减少字典矩阵的维数。下面我们将Farey序列结构字典矩阵简称为Farey结构字典(FD)。

2.2 Ramanujan结构字典矩阵

Ramanujan 结构字典是由印度数学家Ramanujan提出的[19]，Vaidyanathan进一步发展了Ramanujan 结构矩阵， 其形式为[17-18]

(7)

其中:当0≤n≤q-1时，构成了周期序列cq(n)的主值区间。例如当q∈{1,2,…,6}时，则在主值区间，cq(n)的数值为

c1(n)={1},c2(n)={1,-1},c3(n)={2,-1,-1},c4(n)={2,0,-2,0},c5(n)={4,-1,-1,-1,-1},c6(n)={2,1,-1,-2,-1,1} 定义一个q×1矢量:

cq=[cq(0)cq(1) …cq(q-1)]T

它是由周期序列cq(n)主值区间的数值构成的。进一步，定义一个q×φ(q)矩阵：

(8)

如果序列的最大周期为Pmax，构造一个L×Φ(Pmax)维字典矩阵为

(9)

3 混叠脉冲序列的子周期估计

假设有M个具有不同周期的周期脉冲信号同时到达接收装置，对接收装置的测量值进行采样，得到的采样序列可表示为

(10)

式中：tk由式(2)给出。根据采样定理，采样间隔为

(11)

本文采用的均匀采样是一种简单而有效的方法，如果不考虑计算复杂性，尽可能使用较小的采样间隔。当所有tk为整数，且min(tk-tk-1)≥2 (k=2,…,K)时，取Δt=1,在这种情况下，n=1,2,…,tK。将采样序列y(n)写成矢量形式，即Y=[y(1)y(2) …y(tK)]T。利用采样序列和Ramanujan结构字典和Farey结构字典矩阵，采取稀疏重构方法估计混叠脉冲序列的子周期。

考虑下列稀疏优化问题：

(12)

本文采用了联合l2,0混合范数(JLZA)算法[20-22]求解稀疏优化问题，如式(12)所示。相比于文献[16]提出的求解方法，JLZA算法对周期抖动和数据缺失的鲁棒性更强。因为Ramanujan结构字典针对各个周期都有φ(q)列，因此，X的最优稀疏解X*是一个块矢量，即最优稀疏解中，对应于每个子周期的φ(q)个连续元素构成一个块矢量。利用式(13)计算最优稀疏解中的各个块矢量的总系数G(P)：

(13)

当混叠脉冲序列含有某个隐含子周期时，最优稀疏解中对应于这个子周期的总系数G(P)将出现较大数值。

4 仿真实验

4.1 不含有隐含子周期的单周期脉冲序列

假设一个不含有隐含子周期的单周期脉冲序列，如图1所示。常数T=15，采样间隔Δt=1，采样序列的长度为400，幅值为1的脉冲个数为26。在没有脉冲滑动误差和脉冲缺失的情况下，图3给出了3种方法的各个周期总系数G(P)随周期P的变化情况。

从图3(a)看出，在周期为15时，G(P)有明显峰值，这说明JLZA-RD方法能够正确估计出单周期脉冲序列的周期。从图3(b)看出，在周期为1、3和5时，G(P)有最大的峰值。但在周期为15时，其峰值相对较小。这说明很难从G(P)的大小来判断这个脉冲序列的周期。因此，2NCE-RD方法失效。从图3(c)看出，在周期为1、3、5和15时，G(P)有相等的峰值，因此，与2NCE-RD方法个脉冲序列的周期。在脉冲缺失个数一定时，图4给出了4种方法的成功概率随滑动误差变化曲线。在滑动误差一定时，图5给出了4种方法的成功概率随脉冲缺失个数变化曲线。估计成功的概率定义为运行100次Monte Carlo实验，正确估计出脉冲序列周期的次数除以100得到的数值。从图4看出，JLZA-RD方法的估计性能最好。当参数δ≤0.15时，JLZA-RD方法的成功概率为1。这说明JLZA-RD方法的抗噪能力强。在脉冲缺失个数较少时，JLZA-FD方法的估计性能很好。这说明JLZA-FD方法的抗噪能力较强。2NCE-RD方法的估计性能最差，其成功概率小于0.2。其原因是G(P)不仅在周期为15处出现峰值，而且在周期为15的除数处也出现虚假峰值，且在大多数情况下，其虚假峰值高度大于在周期为15处的峰值。对于PDMax方法，滑动误差较大时的成功概率大于滑动误差较小时的成功概率，其原因是滑动误差的增大，使得在周期为1、3和5处的峰值高度降低，其峰值高度小于在周期为15处的峰值高度。因此，其成功概率增加。从图5看出，JLZA-RD方法的估计性能最好。在参数δ=0.1时，其成功概率不随着脉冲缺失个数增加而减少，始终为1。这说明JLZA-RD方法的抗脉冲缺失个数能力较强。在脉冲缺失个数较少时，JLZA-FD方法的估计性能很好，但是，当脉冲缺失个数较多时，JLZA-FD方法的估计性能急剧下降,这说明JLZA-FD方法的抗脉冲缺失能力差。2NCE-RD方法的估计性能最差，其成功概率小于0.2。其原因如上所述。对于PDMax方法，其成功概率随着脉冲缺失个数增加而减少，但是，减少的程度不大。这说明PDMax方法有一定的抗脉冲缺失能力。

4.2 隐含子周期的混叠脉冲序列

假设一个具有隐含子周期的混叠脉冲序列，如图2所示。一个脉冲序列的周期T=15，初相θ0=1；另一个脉冲序列的周期T=18，初相θ0=1。采样间隔Δt=1，采样序列的长度为400，幅值为1的脉冲个数为48。在没有脉冲滑动误差和脉冲缺失的情况下，图6给出了3种方法的各个周期总系数G(P)随周期P的变化情况。

从图6(a)看出，G(P)有两个明显的峰值，分别位于周期为15和18处。这说明JLZA-RD方法能够正确估计出混叠脉冲序列的两个隐含子周期。从图6(b)看出，G(P)有多个峰值，其两个最大峰值分别位于周期为1和3处，而在周期为15和18处，其峰值较小。根据G(P)取得最大值处的周期是估计周期的判断原则，2NCE-RD方法不能正确估计出混叠脉冲序列的两个隐含子周期。因此，2NCE-RD方法失效。从图6(c)看出，PDMax方法与2NCE-RD方法相似，它也没有正确估计出混叠脉冲序列的两个隐含子周期。因此，PDMax方法失效。

在脉冲缺失个数一定时，图7给出了4种方法的成功概率随滑动误差变化曲线。在滑动误差一定时，图8给出了4种方法的成功概率随脉冲缺失个数变化曲线。

从图7看出，JLZA-RD方法的估计性能较好。在滑动误差较小时，JLZA-RD方法具有较高的成功概率。这说明JLZA-RD方法有一定的抗噪能力。在没有脉冲缺失时，JLZA-FD方法的估计性能很好。这说明JLZA-FD方法的抗噪能力较强。对于2NCE-RD和PDMax方法，其成功概率接近于0。这说明这两方法完全失效。其原因是G(P)具有虚假峰值，且虚假峰值的高度较大。从图8看出，JLZA-RD方法的估计性能较好。在滑动误差较小时，其成功概率ξ≤0.35时，始终为1。这说明JLZA-RD方法的抗脉冲缺失能力较强。JLZA-FD方法受脉冲缺失个数的影响较大。当脉冲缺失个数增加时，JLZA-FD方法的估计成功概率急剧下降。这说明JLZA-FD方法的抗脉冲缺失能力差。对于2NCE-RD方法和PDMax方法，其成功概率接近于0。这说明这两方法完全失效。

当改变混叠脉冲序列的初相时，即一个脉冲序列的周期T=15，初相θ0=5；另一个脉冲序列的周期T=18，初相θ0=6，而其他条件不变。在没有脉冲滑动误差和脉冲缺失的情况下，图9给出了3种方法的各个周期总系数G(P)随周期P变化情况。与图6比较，图9(a)没有变化，而对于图9(b)，峰值的位置和幅度都有较大变化，对于图9(c)，峰值的幅度有较大变化。这说明初相的改变对JLZA-RD方法没有影响，而2NCE-RD和PDMax方法对初相的改变是敏感的。因此，不同的初相，2NCE-RD和PDMax方法的估计性能可能会变化。

为了进一步验证初相变化对方法性能的影响。假设初相为均匀分布的随机变量，而其他条件与图9相同。在此条件下，对具有随机初相的混叠脉冲序列进行100次周期估计实验，图10给出成功估计出序列周期的概率曲线。从图10看出，在滑动误差和脉冲缺失都为零的情况下，JLZA-RD方法对初相的随机变化不敏感，其成功概率始终为1；而2NCE-RD和PDMax方法受初相改变的影响敏感，其成功概率始终为0。与图7和图8对比，在存在滑动误差和脉冲缺失的情况下，JLZA-RD方法的成功概率变化较小。这说明了本文方法在各种情况下，对初相变化都具有较强的鲁棒性。

通过上述实验，可以得到如下结果：

1) 基于结构字典的稀疏重构方法具有较好的性能。JLZA-RD方法具有较强抗噪能力和抗脉冲缺失能力，并且不受初相的改变影响。JLZA-FD方法具有较强抗噪能力，但是，抗脉冲缺失能力较差。其原因是稀疏重构方法使用了联合l2,0混合范数来进行优化，即确保未知矢量的非零元素的个数为最少。因此，通过稀疏优化迭代，不断减小具有较小能量的峰值幅度，使其趋于零；同时不断增大具有较大能量的峰值幅度，使其趋于最大。最终在脉冲序列的周期处，得到了幅度最大的峰值。

2) 2NCE-RD方法具有较差的性能，尤其是在具有隐含子周期的混叠脉冲序列的情况。其原因是2NCE-RD方法使用了l2范数来进行优化，即确保未知矢量的每个元素平方之和为最少。因此，相对于联合l2,0混合范数优化，l2范数优化会使非零元素的个数增加。最终在多个周期处，出现幅度不等的峰值。这会给正确识别脉冲序列的周期带来很大困难。另外，2NCE-RD方法对混和脉冲序列初相的改变是敏感的。

3) PDMax方法具有较差的性能，尤其是在具有隐含子周期的混叠脉冲序列的情况。其原因是PDMax方法使用了脉冲序列的功率谱来估计序列的周期。因此，与l2范数类似，会在多个周期处，出现幅度不等的功率谱峰值。这会给正确识别脉冲序列的周期带来很大困难。

4.3 数据长度对脉冲序列估计性能的影响

在这个实验中，对于单周期脉冲序列，我们给出了能够正确估计出脉冲序列周期所需要满足的最小序列长度。这里，我们将脉冲序列周期分解成素数的乘积。图11给出了当脉冲序列周期P为素数(周期从2到60)时，序列长度随周期P的变化。在图11中，由星号标记的数据表示没有脉冲滑动误差的情况下，估计出脉冲序列周期所需要满足的最小序列长度；而由正方形标记的数据表示有脉冲滑动误差δ=0.1的情况下，估计出脉冲序列周期所需要满足的最小序列长度。两条曲线是由指数函数拟合序列长度数据所得到的结果。从图11看出，周期的大小是影响所需最小序列长度的关键因素。周期增大，所需最小序列长度近似按指数曲线增大。

当脉冲序列周期P为两个以上素数(素数1除外)乘积时，在没有脉冲滑动误差的情况下，序列长度随周期P(周期从4到40)的变化如表2～表4所示。从表2看出，所需最小序列长度不仅与周期的大小有关，还与分解的两个素数中的最大素数有关。最大素数的数值增大，所需最小序列长度相应增加。例如，周期为14(素数分解为2×7)和15(素数分解为3×5)时，虽然14比15小，但是由14分解的两个素数中的最大素数是7，而由15分解的两个素数中的最大素数是5，因而，估计周期14所需最小序列长度大于估计周期15所需最小序列长度。

从表3和4看出，所需最小序列长度还与分解的素数多少有关。分解素数的个数增多，所需最小序列长度相应增加。例如，在周期为36(素数分解为2×2×3×3)时，分解素数的个数是4，因而，估计周期36所需最小序列长度大于估计周期37、38(素数分解为2×19)、39(素数分解为3×13)所需最小序列长度，其中，周期37分解素数的个数是1，而周期38和39分解素数的个数是2。总之，对于单周期脉冲序列，能够正确估计出脉冲序列周期所需要满足的最小序列长度不但与脉冲序列周期的大小有关、还与分解素数的个数和其中的最大素数有关。

表2 序列长度随周期P的变化(P分解为2个素数乘积)

表3 序列长度随周期P的变化(P分解为3个素数乘积)Table 3 Variation of sequence length with period P (P is decomposed into prouct of three prime numbers)

表4序列长度随周期P的变化(P分解为4个以上素数乘积)

Table4VariationofsequencelengthwithperiodP(Pcanbedecomposedintoproductoffourormoreprimenumbers)

复合数24364032素数分解2×2×2×32×2×3×32×2×2×52×2×2×2×2序列长度260498550338

5 结 论

1) 仿真实验表明，基于结构字典的稀疏重构方法具有较好的性能。JLZA-RD方法具有较强抗噪能力和抗脉冲缺失能力，并且不受初相的改变影响。JLZA-FD方法具有较强抗噪能力，但是，抗脉冲缺失能力较差。2NCE-RD和PDMax方法具有较差的性能，尤其是在具有隐含子周期的混叠脉冲序列的情况。

2) 对于单周期脉冲序列，能够正确估计出脉冲序列周期所需要满足的最小序列长度不但与脉冲序列周期的大小有关、还与分解素数的个数和其中的最大素数有关。