跳频通信因其优越的抗干扰、抗截获和多址组网等能力，在军民用信息系统中得到了广泛应用[1-2]．对截获的多跳频信号跳周期、频率及跳时参数进行准确盲估计是生成技术参数情报、网台分选及信息恢复等任务的核心环节[3-4]．目前针对该问题，一般基于多通道侦收，即使偶有涉及单通道算法，也大多面向异步或同步网台，极少能处理同步、异步网台信号共存情况．单通道对侦察载荷/平台要求低，是弹载、微/小无人机载等便携侦察系统的主体架构，故亟须开展基于单通道的多跳频信号参数盲估计算法研究．对截获信号在时-频二维域上进行联合处理能凸显其频率跳变特征，易于提取信号参数，是目前开展跳频信号参数估计的主流思路[5-6]．时频分析类方法最为典型，优点是需要的先验信息少，易实现盲估计，但大多适用于同步网台或单信号情况[7-9]．文献[10]提出了两个异步网台信号的盲估计算法，扩展了时频分析法的应用范围，开启了多跳频信号参数估计新思路，但无法用于同步网台信号．文献[11]通过谱熵结合脉冲重复间隔变换，实现了多异步网台跳周期估计，但也不能用于同步网台．文献[12]提出了一种短时傅里叶变换和平滑伪魏格纳分布相结合的估计算法，具有较好信噪比适应能力，但要求多通道接收，且须要信号分选．文献[13-14]将跳频信号分解为多个时频原子，再利用时频原子估计跳频参数．该类方法信噪比适应能力较强，但计算量较大，也只适用于单跳频信号．文献[15-18]将稀疏重构思路引入跳频信号时频图生成/修正，提出了基于稀疏线性回归[15]、稀疏贝叶斯等跳频参数盲估计方法[16-18]．此类方法理论上能处理多信号情况，可取得比时频分析更好的时频聚焦能力，但构建稀疏分解基、重构算法参数优化等难点问题还须进一步发展解决．本研究提出了一种基于单通道侦收，对同步、异步网台均具适应性的多跳频信号参数盲估计算法，且能实现异步网台频率集的自动分选．1 数学模型与问题描述设某跳频信号sn(t)跳周期为Rn，在观测时间Δt内共包含K个跳(hop)，第k(k=1,2,…,K)跳载频为ωnk，起始跳时长(跳时)为FT,n，sn(t)可表示为[3]sn(t)=vn(t)∑k=0K-1expj(ωnkt'+φnk)rect(t'/Rn)，(1)式中：vn(t)为sn(t)的基带复包络；t为时间；φnk为第k跳的初相；t'=t-(k-1)Rn-FT,n；rect(t)为单位矩形窗．假设在Δt内有N(≥2)个异步或同步跳频信号同时进入单通道接收机，并以采样频率fs对接收信号进行离散化，则截获样本z(l)可表示为z(l)=∑n=1Nsn(l)+v(l)  (1≤l≤M)，(2)式中：v(l)为零均值、方差为σ2的加性高斯白噪声；M为样本长度．各跳频信号第一个非完整跳为跳时参数．两个跳频信号的参数示意如图1所示，图中：f为频率；R1和R2分别为两信号跳周期；FT,1和FT,2分别为两信号跳时；αk (k=1,2,…,L)为所有起跳时刻．10.13245/j.hust.230624.F001图1两个跳频信号的参数示意图多跳频信号参数盲估计就是基于z(l)实现对各跳频信号频率、跳周期Ri(i=1,2,…,N)、跳时参数FT,i的准确估计．将图1所示所有αk 组成的集合称为起跳时刻序列，记为Th={α1,α2,…，αL}．相邻起跳时刻持续的时间段定义为频率驻留期，可见起跳时刻序列蕴含着跳频信号跳周期和跳时信息，其对应关系由不同情况下的起跳时刻分布特征决定，各频率驻留期内信号可建模为N个窄带定频叠加．2 跳频信号起跳时刻的分布特征可将截获的跳频样本归纳为三种场景，如图2所示，图中：Δt为观测时间；αk (k=1，2，…，L)为截获样本的L个起跳时刻．10.13245/j.hust.230624.F002图2不同场景下网台的跳频示意图场景1 截获样本为N个同步网台跳频信号同步网台(图2(a))同时跳频，理论上起跳时刻序列中各元素等间隔分布，Th为等差数列．α1等于跳时参数FT，Th差分结果对应跳周期Rh，即：FT=α1; (3)Rh=αl1+1-αl1, (4)式中l1=1,2,…,L-1．场景2 截获样本为N个异步网台跳频信号为不失一般性，假定N个异步网台中存在不等速网台，此时序列Th整体已无规律可循，但Th中各跳频信号相应起跳时刻元素构成的子序列仍为等差数列，其公差分别为各信号跳周期，即满足Th=∪i=1NThi; (5) Thi=αi|αi=αi,1+(ki-1)Ri  (i=1,2,…,N;ki=1,2,…,Li) ; (6)L=∑i=1NLi，(7)式中：Thi为Th中第i个跳频信号对应起跳时刻元素αi(子序列)所构成的子集；Ri为第i个跳频信号跳周期；Li为第i个跳频信号频率跳变次数，即子集Thi包含的元素数目．当N个异步网台等速跳频时，式(5)～(7)可简化，此时Th间隔N-1个元素的差分结果对应跳周期，前N个元素值分别对应各信号跳时，即FT,i=αi；(8)Rh=αli+N-αli，(9)式中：FT,i为各信号跳时参数；li=1,2,…,L-N-1．特别地，当N=2时，根据式(5)～(7)，不等速时(图1)，容易得到：Th1=αi|αi=αi,1+(k1-1)R1(k1=1,2,…,L1); (10)Th2=αj|αj=α2,1+(k2-1)R2(k2=1,2,…,L2) . (11)同理，根据式(8)和(9)，可得到等速时(图2(b))跳时参数分别为α1和α2，跳周期Rh可由下式估计Rh=αl1+2-αl1=αl2+2-αl2，(12)式中：l1(=1,3,5,…,L-3)为信号1起跳时刻；l2=(2,4,6,…,L-2)为信号2起跳时刻．场景3 截获样本为Ns(≥2)个同步网台与N-Ns个异步网台信号的组合基于对场景2分析可知：当场景3中同步与异步网台跳周期相同时，Th分布规律同式(5)、式(8)和(9)；反之只要存在不等速网台，Th分布规律同式(5)～(7)．通过以上分析可知：对于场景2和场景3，须对其起跳时刻序列进行分选，才能基于上述公式进行跳周期和跳时参数估计．3 参数盲估计与频率分选参数估计的基本思路和流程是：首先进行起跳时刻检测提取，然后分选各信号的起跳时刻子序列，最后利用上一节推导结论进行参数估计．3.1　基于改进谱熵的起跳时刻检测提取谱熵反映了信号频谱的分布状态，若信号在某时刻发生频率跳变，该时刻频谱被展宽，谱熵值出现正向突变[11,19]，可据此原理检测跳频信号起跳时刻．但直接通过文献[11]谱熵算法检测起跳时刻，信噪比适应能力较弱，且存在“频率选择性(随机出现谱熵结果不平坦，影响检测判决)”问题．本研究所提方法能解决上述问题，包括信号预处理、谱熵计算及起跳时刻检测提取三个环节．3.1.1　信号预处理对谱图矩阵SPEC(t,f)进行截断，得到STRU(t,f)再进行后续处理，以提高抗噪能力，计算公式为STRU(t,f)=SPEC(t,f)  (SPEC(t,f)≥ε);0              (SPEC(t,f)ε), (13)式中ε为截断门限，由SPEC(t,f)直方图得到．跳频信号分量在时频面上具有稀疏特性，各hop时频区域只占矩阵元素少部分．时频矩阵幅度的直方图峰值表征了噪声强度，截断门限ε的取值略大于该值即可．显然，直方图法计算截断门限可使门限值随着噪声强度自适应变化．3.1.2　新的谱熵计算方法设STRU(t,f)的离散形式为STRU(n,k)(n=1,2,…,(L-Nfft)/Lstep+1；k=1,2,…,Nfft)，式中：Nfft为短时傅里叶变换的短时窗长；Lstep为短时窗滑动步进．z(l)在第n窗的谱密度定义为P(n,k)=STRU(n,k)  /∑k=1NfftSTRU(n,k)．(14)设置频域累加因子Nsum(Nsum为整数且Nsum≥2)，将傅里叶变换长度Nfft平均划分成ΔN=Nfft/Nsum段，修正的谱密度为P^(n,k^)=∑(k^-1Nsum+1k^NsumP(n,k)，(15)式中k^=1,2,…,ΔN．容易证明：0≤P^(n,k^)≤1，∑k^=1ΔNP^(n,k^)=1．根据香农信息熵定义，有Dz(n)=-∑k^=1ΔNP^(n,k^)logP^(n,k^) (n=1,2,…,(L-Nfft)/Lstep+1) (16)为信号z(l)第n窗的谱熵值．不同窗时刻的谱熵值其集合构成跳频信号谱熵曲线．仿真生成跳速1 000 hop/s的两个跳频信号，信噪比为4 dB，短时傅里叶变换采取128点矩形窗，滑动间隔8个点，谱图修正结果及不同情况下的谱熵曲线见图3，其中：n为时间窗索引；Dz为各窗数据的谱熵值．10.13245/j.hust.230624.F003图3谱图修正结果及不同情况下的谱熵曲线由图3(a)可见：文献[11]谱熵曲线毛刺明显且在部分区域基底上移，难以进行起跳时刻检测判决．由图3(d)可见谱图修正结合新的谱熵方法得到的谱熵曲线非常平滑，在图3(b)各起跳时刻位置准确出现峰值脉冲，各频率驻留期内谱熵曲线平坦，原因为：谱图修正抑制了噪声，同时式(15)的频域累加处理将时频矩阵(若干)相邻行视为一个频率分量，解决了“频率选择性”问题(图3(c))．3.1.3　起跳时刻检测提取侦察属于非合作方，得到谱熵曲线后应先估计峰值个数(起跳时刻数目)及位置．本研究通过对谱熵曲线幅度进行直方图统计以生成峰值检测判决门限，如对于图3(d)，其直方图最大值为0.3，与图3(d)平坦段幅度对应正确．将直方图最大值记为Γ，谱熵判决门限β计算公式为β=Γ+ξ(max[Dz(n)]-Γ)，(17)式中ξ为门限调整因子，可根据信噪比条件灵活选择．由于谱熵曲线经门限β的判决结果为若干邻域(邻域数目等于跳频次数，邻域中心对应起跳时刻)，故本研究采取聚类思想结合自适应判决，对邻域数目和邻域中心进行估计．定义Dz(n)β的幅度集合为Λ={Λ1, Λ2,…,ΛQ}，显然QM，估计流程如下：步骤1 取出Λ中最大值元素max[Λ]及对应时刻kmax，即Dz(kmax)=max[Λ]；步骤2 将kmax向两侧逐点扩散，若同时满足Dz(kmax-l)≥β和Dz(kmax+l)≥β(l=1,2,…) ，认为kmax-l至kmax+l时刻处于同一邻域，反之认为该邻域结束，记录该邻域的位置、宽度及各元素幅度；步骤3 将Dz(n)中最大值邻域所有元素置零；步骤4 返回步骤2，直至Dz(n)中所有元素幅度均小于β，退出该流程．至此，得到了满足Dz(n)≥β的邻域数L和各邻域中心时刻{k1,k2,…,kL}，则起跳时刻序列为Th={k1,k2,…,kL}Lstep={α1,α2,…,αL}. (18)3.2　起跳时刻及频率分选本节所提方法能在分选场景2和场景3(基于序列Th判断信号场景的方法将在3.4节阐述)起跳时刻序列的同时，实现对各跳频信号频率估计与分选．首选根据Th将信号z(l)分成L+1段频率驻留期，得到L+1段样本zl(n) (l=1,2,…,L+1)；再利用魏格纳-威利分布结合快速傅里叶变换，估计各频率驻留期内的频率分量；最后通过检测相邻频率驻留期发生跳频的频率分量，利用跳频信号“频率跳变性”特性，实现对起跳时刻和异步网台频率的高效分选．3.2.1　频率驻留期内的频率分量估计利用快速傅里叶变换或者时间累加得到第l个驻留期的第i个频率的初估计值ρl,i，即ρl,i=argmax iZl(k)       (λl≥Nfft); argmaxi∑n=Th(l-1)Lstep +1n=Th(l)Lstep STRU(n,k)   (λlNfft), (19)式中：Zl(k)为zl(n)的快速傅里叶变换结果；λl=Th(l)Lstep-Th(l-1)Lstep为zl(n)长度．当λlNfft时，快速傅里叶变换可取得更好频率精度，故提取Zl(k)前N个峰值作为频率初估值，反之对谱图模值按时间累加后再进行峰值检测来实现频率初估计．对zl(n)进行魏格纳-威利分布，得到Wl(n,k)，则各频率驻留期内的N个频率估计值为f^l=argmaxk=1,2,…,N∑n=Th(l-1)Lstep+1n=Th(l)Lstep Wl(n,k) ．(20)Wl(n,k)存在交叉项，因此当对式(20)进行峰值提取时，只在ρl,i邻域内进行．不同频率驻留期的频率分量构成频率矩阵F，即F=[f^1,f^2,…,f^L+1]．(21)3.2.2　起跳时刻及频率分选图4为场景2的分选示意图．对图4的F矩阵各列元素从高至低排序，按照f^l,n(l=1,2,…，10.13245/j.hust.230624.F004图4场景2分选示意图L+1;n=1,2,…,N)标记．起跳时刻和频率分选的基本思路是：在相邻两个频率驻留期内发生跳频的频率属于同一跳频信号，同时该起跳时刻也属于该跳频信号；而在相邻两个频率驻留期内频率保持不变的属于对应信号(这些信号互斥于发生跳频的信号)的同一频率．因此，针对某一起跳时刻，可通过检测与其对应的相邻频率驻留期内各频率分量的跳变性或序贯连续性，来确定该起跳时刻及相应频率分量与各跳频信号的隶属关系(即实现分选)．将式(21)中f^1,1～f^1,N定义为信号1～N频率集的第1个频率分量，即f^1,n∈As,n(n=1,2,…,N)．如图4所示，对于α^1，第1频率驻留期内f1,1在第2频率驻留期跳变为f2,2，而f^1,2和f^1,3分别持续为f^2,1和f^2,3，因此α^1隶属于信号1的第1个起跳时刻，同时f1,1和f2,2均属于信号频率集As,1，并将f^1,2和f^2,1进行合并，判其隶属于信号2频率集As,2，f^1,3和f^2,3进行合并，判其隶属于信号3的频率集As,3；对于α^2及后续所有起跳时刻依次采取上述方法，最终实现对整个起跳时刻序列和频率矩阵F的分选．频率分量f^l,n跳变性和持续性的检测方法为{k,j}=arg mink,j{Ekl-Ejl+122}；(22)Cl={Ekl-Ejl+122}；(23)Hl={f^l}-{Ekl}；(24)Hl+1={f^l+1}-{Ejl+1}，(25)式中：Ekl为矩阵El的第k列矢量，El为由f^l所有N-1维子集构成的矩阵，维度为(N-1)N；Cl为在第l和第l+1个驻留期内不变的频率分量构成的集合；Hl和Hl+1分别为在第l和第l+1个驻留期内发生跳频的频率分量构成的集合．场景3的起跳时刻分选方法与场景2相同，并能得到同步网台频率集和异步网台跳图案．3.3　多网台信号参数盲估计3.3.1　场景1参数估计根据式(3)和(4)，跳周期和跳时估计值为：Rh =1L-1∑l=1kL-1k(l+1)-k(l)Lstep=kL-k1L-1Lstep；(26)FT=α1=k1Lstep．(27)频率通过式(20)估计即可．3.3.2　场景2和场景3参数估计根据3.2节方法得到时频矩阵F及分选得到各个跳频信号的起跳时刻子序列．根据式(8)和(9)可得第n个信号跳周期和跳时估计值分别为：Rn=∑l=0N-1kL-l-∑l=1Nkl/(L-N)Lstep; (28)FT,n=αn=knLstep  (n=1,2,…,N). (29)对于异步网台，根据起跳时刻对信号分段造成hop“断裂”，导致基于式(19)和(20)估计频率未将各hop数据完全利用．在起跳时刻分选后，各hop已可以在时频图上“定位”，此时可根据“定位”结果对样本z(l)进行对应截取，再通过式(20)计算各hop频率值，以取得更佳频率分辨性能．3.4　算法步骤假定截获的跳频信号数目已知或已由文献[7]提出的聚类方法得到，本算法步骤总结如下：步骤1 对信号z(l)求谱图，通过式(13)对谱图进行预处理；步骤2 根据式(14)～(16)求时频图各时刻谱熵值，得到谱熵曲线；步骤3 基于3.1.3节方法对谱熵曲线进行峰值提取，通过式(18)得到起跳时刻序列Th；步骤4 对起跳时刻序列Th求二阶差分，若均值和方差趋于0，则判定为场景1，转入步骤5，反之转入步骤6；步骤5 通过式(26)和(27)求跳周期、跳时参数，通过式(19)和(20)得到频率集；步骤6 利用3.2节方法分选起跳时刻序列和频率分量，基于分选结果截取各hop时域样点；步骤7 通过式(28)和(29)求跳周期、跳时参数，通过式(19)和(20)得到各跳频信号频率集/跳图案．3.5　算法复杂度分析忽略加法运算，主要考虑算法的乘法运算次数．谱图生成模块乘法次数为[log2Nfft(Nfft+2)(L-Lstep)]/(2Lstep)，谱图修正模块乘法次数为Nfft(L-Lstep)/Lstep，谱熵曲线生成模块运算量为[(Nfft+ΔN)(L-Lstep)]/Lstep次乘法与[Nfft(L-Lstep)]/Lstep次对数运算之和，跳频参数估计模块乘法次数可近似为(2L)log2L/(αL+1)，故总复杂度近似为[(4Nfft+2ΔN)+log2Nfft(Nfft+2)](L-Lstep)/(2Lstep)．由总复杂度表达式可知，减少快速傅里叶变换长度Nfft或增大滑动步进Lstep，均可显著降低复杂度．在本研究的仿真条件下，综合考虑算法复杂度和参数估计准确度，可取Nfft=128，Lstep=8．4 仿真与分析文献[10-11]是基于单通道和时频分析开展异步网台信号参数盲估计的典型代表，文献[11]提出了基于谱熵的异步网台信号跳周期估计算法，文献[7]提出了同步网台信号参数盲估计算法．在此构造三类侦察场景，对比所提算法与上述文献算法的性能．实验1 等速异步网台的估计性能设置两个1 000 hop/s跳频信号，两信号跳时因子分别为0.30和0.75；采样率为1.024 MHz，观测时长为6 500个样点；短时傅里叶变换(STFT)采用128点矩形窗，滑动间隔为8个点；信噪比RSN变化范围为-4～20 dB，各信噪比下进行1 000次蒙特卡罗试验的统计结果见图5，其中ERMS为均方根误差．图5表明：本算法优于文献[10-11]，原因是本研究采取的谱图修正结合新谱熵方法及起跳时刻检测分选具有更好的噪声抑制能力，同时频率估计采取了具有最佳分辨率的魏格纳-威利分布．10.13245/j.hust.230624.F005图5等速异步网台信号的参数估计结果由图5(a)可见：当RSN≥4 dB时，频率估计误差ERMS＜1×10-3，证明了算法对多信号频率的优良估计性能；当RSN进一步提高时，估计精度稳定于1×10-3～1×10-4，出现了“地板效应”(采样点有限，导致魏格纳-威利分布和快速傅里叶变换估精度受限)．图5(b)表明：跳时参数估计精度随RSN的提高不断提高，主要因为当RSN提高时，谱熵曲线各邻域呈对称形状，邻域中心时刻准确对应各跳时刻．由图5(c)可见：当RSN≥8 dB时，跳周期ERMS＜1×10-3；当RSN＜4 dB时，ERMS值急剧增大，存在门限效应，主要因为跳周期估计受跳时刻序列影响很大；当RSN较低时，起跳时刻序列各邻域中心与真实值存在一定偏差．通过增加采样数据的长度，能够降低RSN门限．实验2 对不等速异步网台的估计性能设置两个跳频信号，跳速分别为800 hop/s和1 000 hop/s，其余条件与实验1相同．算法估计结果如图6所示，可见：当不等速跳频时，本算法优于文献[10-11]算法，原因与实验1相同．由图6可见：总体来看，所提出算法信号2的性能略优于信号1，主要由于在有限采样时长内，信号2的hop数更多，当平均时取得了更好的结果．对比图6(a)和图5(a)可见：当不等速跳频时，频率估计精度略差于等速跳频，主要由于当不等速跳频时，跳时刻序列的分布更加随机，容易出现短数据段情况．文献[10]算法频率估计性能对异步网台跳速是否相同不敏感．对比图6(b)、图6(c)和图5(b)、图5(c)可见，本算法跳时和跳周期参数估计性能与等速跳频时相差无几．10.13245/j.hust.230624.F006图6不等速异步网台信号的参数估计结果实验3 对同步网台的估计性能设置两个同步跳频信号，跳时因子为0.30，其余条件与实验1相同．本算法与文献[7]算法的性能统计结果如图7所示，可见：本算法频率估计性能优于文献[7]算法的，原因在于本研究采用具有最佳时频聚焦能力的魏格纳-威利分布，且各hop样点全部参与了频率计算；本算法跳周期、跳时估计性能与文献[7]算法的相差无几，跳时性能甚至略逊于文献[7]算法的，原因在于本算法直接将起跳时刻序列第一元素作为跳时估计值，而文献[7]算法跳时估计分为粗估计和精估计两步处理．10.13245/j.hust.230624.F007图7同步网台信号的参数估计结果对比图5~7可见：本算法的同步网台估计性能优于异步，主要因为当同步组网时，无须信号分段和跳频时刻、频率分选，频率估计时各跳数据得到完全利用，起跳时刻序列估计结果也更加接近真实值．图7表明：当RSN≥4 dB时，频率误差ERMS＜1×10-3；当RSN≥6 dB时，跳周期误差ERMS＜1×10-3并持续减小．通过增加采样数据长度，能够进一步提高跳周期和跳频率估计性能．5 结语本研究在建立跳频信号单通道侦收模型基础上，将跳频信号侦察归结为三类典型场景，推导分析了各场景下起跳时刻序列理论分布和跳频信号时域参数估计数学表征；然后设计提出了起跳时刻序列检测分选、频率分量高精度估计方法，以及基于起跳时刻分选结果的跳频信号跳周期、跳时参数估计方法；最后实现了基于单通道接收的多跳频信号参数盲估计，并实现了异步网台频率分选．本研究所提出的算法门限设置少，计算效率高，对不同组网方式均具有适应性．下一步重点研究如何进一步提高本算法的信噪比适应能力和干扰剔除能力．