CN111934711B - 一种时频混叠跳频信号的参数估计方法 - Google Patents

Abstract

本发明公开一种时频混叠跳频信号的参数估计方法，其步骤为：将单通道接收下的时频混叠信号，利用短时傅里叶变换转换到时频域，得到时频域信号；利用基于稀疏线性回归的矩阵优化算法对时频域信号进行去除背景噪声、信号畸变特征处理；利用基于二次包络优化的参数估计算法对时频域信号进行去除部分干扰信号特征和异常点处理，提取优化后时频域信号的平均时频脊线进行平滑处理；对平均时频脊线进行拐点检测，再对检测到的拐点进行时频域映射，完成跳频信号的跳周期估计；在跳频周期估计的基础上，基于Hough变换的优化方法完成跳频信号的频点估计。本发明针对各源信号在时、频域上发生混叠的场景，能够实现跳频信号的完整恢复与参数的高精度估计。

Description

一种时频混叠跳频信号的参数估计方法

技术领域

本发明属于信号处理技术领域，尤其是涉及一种时频混叠跳频信号的参数估计方法。

背景技术

跳频扩频信号以其截获概率低、抗干扰性能好、保密性好、安全可靠等优点被广泛运用于军事和商业通信领域，同时也给电子侦察带来了严峻的挑战。在电子侦察技术中，通常需要对信号进行检测以及参数估计，完成信号的侦察。正确的估计出跳频的信号参数是完成信号解跳的重要组成部分之一，因而，研究跳频信号的检测和参数估计具有重要的研究意义。

跳频信号在不同应用场景下的参数估计一直是通信领域的研究热点和难点。跳频信号参数估计是指在没有任何先验知识的情况下估计出跳频信号的跳周期、起跳时刻以及跳频频率。然而，在现代无线通信中，电磁环境日益复杂，各式各样的信号极易同时在时、频域上发生混叠，致使宽带接收下的信号无法直接用于后期的识别与信息提取，此时跳频信号的参数估计和追踪难度大大增加。

单通道接收场景下，跳频信号参数估计的参数化方法一般是将有效频率段建模为分段常数，并以更高的复杂度为代价实现更高的估计精度，但是整个问题的症结在于需要提前给定某些参数的先验信息，如跳速、频率范围等，然而在非协作接收场景下，传输信道和信号的先验信息往往是未知的，此时参数化的估计方法并不适用。非参数化的估计方法主要基于时频分析技术，其特点是不需要信号参数的相关先验信息，同时可以对非平稳信号的时、频两维信息进行展示，但是分辨率有限，信号特征容易发生模糊，在低信噪比下畸变较为严重，导致估计精度较差，大量的工作也对这一问题做了进一步的完善。目前，在时频分析的基础上，通过引入稀疏理论有效增强了参数估计方法在低信噪比下的鲁棒性，提高了估计精度。

目前，单通道跳频信号的参数估计方法主要集中在基于时频分析的研究，如表1所示，大致可以分为以下几类：

上述的参数估计方法都是基于时频稀疏性的优化方法，因此存在有一个明显的缺点：估计性能主要受限于信号在变换域上的稀疏性，且一般都假定输入信号是足够稀疏的；但是现实中存在的大部分接收信号在时、频域都表现为软稀疏性，即不同的源信号在时频域上发生互相混叠，形成时频混叠信号，此时部分时频点上有多个源信号共同存在，信号的稀疏性极差，而传统的参数估计方法无法满足假定的稀疏性要求，使得估计性能急剧下降。

发明内容

针对宽带接收中复杂电磁环境下单通道时频混叠跳频信号的参数估计问题，本发明的目的是提供一种基于稀疏线性回归和二次包络优化的时频混叠跳频信号的参数估计方法，其针对各源信号在时、频域上发生混叠的场景，对单通道跳频信号进行了有效的干扰特征去除和高精度的参数估计，有效降低对输入混叠信号的稀疏性限制，且在低信噪比下显著提高跳频信号的参数估计精度，在应用范围、参数估计精度和低信噪比下的鲁棒性方面都有较大的提升。

为实现上述发明目的，本发明采用如下技术方案：

一种时频混叠跳频信号的参数估计方法，其包括以下步骤：

S1、将单通道接收下的时频混叠信号f(t)，利用短时傅里叶变换STFT转换到时频域，得到的时频域信号F(t，f)；

S2、利用基于稀疏线性回归的矩阵优化算法对步骤S1中得到的时频域信号F(t，f)进行去除背景噪声、信号畸变特征且保留信号的主要特征处理，以改善低信噪比下信号在时频域上的稀疏性和时频分布准确性，处理后的时频域信号系数表示为B(t，f)；

S3、利用基于二次包络优化的参数估计算法，首先利用数据的连续性分析对步骤S2中得到的时频域信号B(t，f)进行去除部分干扰信号特征和异常点处理，然后提取优化后时频域信号B(t，f)的平均时频脊线，进行平滑处理，以减小剩余干扰对跳频信号的影响；

S4、对步骤S3中得到的平均时频脊线进行拐点检测，然后对检测到的拐点进行时频域映射，完成跳频信号的跳周期估计；

S5、最后在步骤S4中得到的跳频周期估计的基础上，基于Hough变换的优化方法完成跳频信号的频点估计。

进一步地，上述的步骤S1中，对单通道接收下的时频混叠信号f(t)进行短时傅里叶变换，过程为：

其中，f(t)表示观测信号，A∈R^1×n为一维的观测矩阵，s(t)表示源信号，v(t)表示加性噪声；

在观测时间[0，T]内，设跳频信号s_i(t)包含M个跳周期，即

其中，K_m∈R和f_m∈[-f_max，f_max]分别为第m跳的幅度和频率，跳周期为T_hop＝t_m-t_m-1(t₀＝0)；

将其转换到时频域来改善信号的稀疏性，表示为：

其中，F(t，f)，S(t，f)和V(t，f)分别为f(t)，s(t)和v(t)在STFT下的时频系数。

进一步地，上述的步骤S2中，利用基于稀疏线性回归的矩阵优化算法对步骤S1中得到的时频域信号F(t，f)进行处理，过程为：

设优化后的时频矩阵表示为B(t，f)∈R^N×P，B(t，f)为二进制矩阵，即矩阵中的非零元素都为1；B(t，f)具有时频点稀疏性和双差分稀疏性两个特性，根据上述两个特性的限制，目标函数写为

其中，F⊙B表示矩阵的Hadamard积，即两个矩阵的同位置元素相乘，D∈R^(P-2)×P，且

在(4)式中，λ₁和λ₂分别控制了估计结果的稀疏性和平滑性，将(4)式中非凸的目标函数l₀范数转化为属于凸函数的l₁范数，同时加入可微的l₂范数，则有

其中，λ₁越大，B(t，f)的稀疏性越好，非零系数越少；λ₂越大，B(t，f)内同行系数的平滑性越好；

设F(t，f)＝[f₁，f₂，...，f_N]^T，f_i∈R^1×P，B(t，f)＝[b₁，b₂，...，b_N]^T，b_i∈R^1×P，则有

其中，w_i＝diag(f_i)，即w_i为对角矩阵，对角线上的元素为f_i；

用L2-ISTA算法对(7)式进行快速迭代求解，首先建立带约束条件的一般化模型：

min{F(b)≡f(b)+g(b)：b∈Rⁿ} (7)

其中，g(b)＝λ₁||b||₁，

对f(b)进行梯度下降来求解，则有

其中，L为下降步长，k代表迭代次数；因为f(b)属于凸函数，所以

有上界，即

||f(b_k)-f(b_k-1)||≤L(f)||b_k-b_k-1|| (9)

其中，||·||表示标准欧几里得范数，L(f)＞0；L＝min(L(f))；基于近端梯度法的思想，对f(b)在α处的泰勒展开为

且有

λ₀≥max[eig(w^Tw+2λ₂D^TD)] (11)

其中，max[eig(w^Tw+λ₂D^TD)]表示w^Tw+λ₂D^TD的最大特征值，φ(α)是与b无关的函数项，忽略不计；因此将(7)式改写为

其中，设λ₀＝max[eig(w^Tw+2λ₂D^TD)]；

由(12)式可知，λ₀等价为

的Lipschitz常数，即λ₀＝L；当α＝b_k-₁时，则有

b_k＝P(b_k-1) (13)

设

且

若H(b_k)取得最小值，则有

其中，sgn(·)为符号函数，基于软阈值算法解得

其中，

表示

的第i个元素，z_i表示z的第i个元素。

更进一步地，上述的(λ₁，λ₂)具有两个选择准则：

准则1、当λ₂＝0时，λ₁的有效取值范围为

且

当λ₂＝0时，(6)式则转化为l₁惩罚回归模型；当

时，则b_i→0，即b_i＝0为(6)式的求解结果；λ₁建议取值为

的5％-10％，此时(6)式能够获得良好的估计结果；

准则2、当λ₁＝0时，λ₂的有效取值范围为

且

当λ₁＝0时，则有

其中，

是非零元素全为1的下三角矩阵，M₀是M的第一列系数；当

时，则有b_i→c_i，其中c_i为常数，即常数向量为(6)式的求解结果；当λ₁≠0时，λ₂建议取值为

的5％-10％，此时(6)式能够获得良好的估计结果。

进一步地，上述的步骤S3中，利用基于二次包络优化的参数估计算法对步骤S2中得到的时频域信号B(t，f)进行处理，过程为：

根据跳频信号的时频分布特性，建立判别准则如下：

准则3、去除宽带信号和部分异常点：当c_w≤th₁，1≤w≤l时，使对应loc_w上的时频系数全为0，即

准则4、针对窄带信号中的定频干扰信号：当c_w≥th₂时，使对应loc_w上的时频系数全为0，即

通过上述的判别准去除时频域上大部分的干扰信号，包括宽带信号、窄带定频干扰信号以及信号因发生畸变而产生的边缘特征；B(t，f)经过准则3、准则4处理后，得到的结果表示为B₁(t，f)；

首先，对B₁(t，f)∈R^N×P每一列的非零值位置进行提取，记为

然后保留ind_R/2处的非零值；对每一列ind_R/2处的非零值进行重排，提取到的脊线称为平均时频脊线，记为mean_tfcurve，其时、频域上分别对应时间和暂态频率；

首先对异常点进行处理：检测mean_tfcurve中相邻点的差异分布，有

M_d＝DM(mean_tfcurve)＝max[Diff(mean_tfcurve)] (18)

其中，DM(·)用于提取目标数据的差分曲线最大值，max(·)表示取最大值操作，Diff(·)表示差分处理；

设平滑处理函数为Sm(·)，且有X＝[x₁，x₂，...，x_n]，则Sm(X)的处理结果为

其中，α为平滑处理阈值，为DM(X)的0.01，设SL＝Sm(mean_tfcurve)，α＝0.01M_d；

采用分段线性插值法对SL的上、下包络曲线U₁和L₁进行提取，过程为：

对U₁和L₁的拐点进行检测，设拐点检测函数为DI(·)，则有

DI(U₁)＝Diff(Diff(U₁))

DI(L₁)＝Diff(Diff(L₁)) (20)

其中，UI₁和LI₁分别表示U₁和L₁中检测到的有效拐点，拐点的检测阈值分别为β·max[DI(U₁)]和β·max[DI(L₁)]；

再进一步对U₁和L₁进行优化处理，过程为：对SL和U₁、L₁的重叠部分进行提取，有

W_u_loc＝Loc(SL∩U₁)

W_l_loc＝Loc(SL∩L₁) (22)

其中，W_u_loc和W_l_loc分别表示SL和U₁、SL和L₁的重叠部分的时域上的索引，Loc(·)为目标对象索引的提取函数，∩表示取交集；

分别对U₁中W_u_loc上的值和L₁中W_l_loc上的值进行平滑处理，则有

U₁_cs＝Sm[U₁(W_u_loc)]

L₁_cs＝Sm[L₁(W_l_loc)] (23)

其中，U₁_cs和L₁_cs分别表示对U₁(W_u_loc)和L₁(W_l_loc)进行平滑处理后的数据集合；对SL重新填充，得到SL_cs，即

对SL_cs进行平滑处理，完成最后的平均时频脊线优化，即

SL_cs＝Sm(SL_cs) (25)。

进一步地，上述的步骤S4中，对步骤S3中得到的平均时频脊线进行拐点检测完成跳频信号的跳周期估计，过程为：

再次对SL_cs进行包络曲线的提取，表示为U₂和L₂；对U₂和L₂的拐点进行检测，表示为UI₂和LI₂；通过提取UI₂和LI₂中平行分布的数据点，即成对平行于时间轴的拐点，并结合SL_cs完成对跳时刻

和跳周期

的高精度估计。

进一步地，上述的步骤S5中，利用基于Hough变换的优化方完成跳频信号的频点估计，过程为：对B₁(t,f)内单个跳周期

内的时频系数进行提取，表示为g_j(t,f)∈R^N×h，其中，h为

内所含列的数目；提取g_j(t,f)内非零值元素的位置，坐标位置表示为

并做极坐标转换：

当[x_z,y_z]与[x_z+i,y_z+i]分布在同一条直线l上时，对应的p_z与p_z+i会在θ＝θ_z处相交，其中θ_z为l和坐标轴的夹角；对

的交点位置和出现次数进行统计，并记出现次数最多的交点为

根据

筛选出相交于该点的曲线{p_z}，即可获得对应的{[x_z,y_z]}，保留g_j(t,f)中{[x_z,y_z]}上的非零值，其他系数全置为0，并重新记

内的时频系数为

提取

内所有非零值元素的纵轴值，并进行平均处理，即

其中，

为

内的频点估计值，f_i表示

内非零值元素的纵轴值，共S个；对B₁(t,f)内不同跳周期内的g_j(t,f)进行上述处理，即完成接收信号内所有跳频频点的估计。

由于采用如上所述的技术方案，本发明具有如下优越性：

该时频混叠跳频信号的参数估计方法，其将混叠信号转换到时频域上以改善信号的稀疏性，并提出一种基于稀疏线性回归的矩阵优化算法对时频矩阵进行优化，实现背景噪声和信号冗余、畸变特征的有效去除，同时保留信号的主要特征，有效改善低信噪比下信号在时频域上的稀疏性和分布准确性；提出二次包络曲线优化算法，首先利用数据的连续性分析对部分干扰源和异常点进行去除，然后提取最优的平均时频脊线进行处理，对跳频信号的信号特征进行了深度优化，并结合相关的特征提取工作完成了跳频时刻、跳周期和跳频频点的高精度估计，以减小时频域上剩余干扰对跳频信号的影响，提升参数的估计精度；在宽带接收中的复杂电磁环境下，在多信号时频混叠的情况下，实现跳频信号的完整恢复与参数的高精度估计，并在低信噪比下具有良好的鲁棒性。

附图说明

图1是b_i内连续分布的非零值段图；

图2是平均时频脊线mean_tfcurve示意图；

图3是SL的包络曲线分布图；

图4是UI₁和LI₁的分布图；

图5是UI₂和LI₂的分布图；

图6是g_j(t,f)的时频分布图；

图7是

的分布图；

图8是

的交点分布图；

图9是各信号源样本的时频分布图；

图10是输入混叠信号的时频分布和频谱图；

图11是B(t,f)的时频分布图；

图12是B₁(t,f)的时频分布图；

图13是优化后的平均时频脊线的拐点分布图；

图14是{t_m}和{f_m}的估计误差图；

图15是

拼接后的时频分布图；

图16是不同SNR下各方法的参数估计误差对比图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步详细说明。

如图1～16所示，一种时频混叠跳频信号的参数估计方法，其包括以下具体步骤：

S1、将单通道接收下的时频混叠信号f(t)，利用短时傅里叶变换STFT转换到时频域，得到的时频域信号系数表示为F(t，f)；过程为：

其中，f(t)表示观测信号，A∈R^1×n为一维的观测矩阵，s(t)表示源信号，v(t)表示加性噪声；这里设定在非协作异步接收场景下，各源信号之间在发射机之间是独立的；

在观测时间[0，T]内，设跳频信号s_i(t)包含M个跳周期，即

其中，K_m∈R和f_m∈[-f_max，f_max]分别为第m跳的幅度和频率，跳周期为T_hop＝t_m-t_m-1(t₀＝0)；

由于时频混叠信号在时频域上的低稀疏性，因此将其转换到时频域来改善信号的稀疏性，表示为：

其中，F(t，f)，S(t，f)和V(t，f)分别为f(t)，s(t)和v(t)在STFT下的时频系数；

S2、利用基于稀疏线性回归的矩阵优化算法对步骤S1中得到的时频域信号F(t，f)进行去除背景噪声、信号畸变特征且保留信号的主要特征处理，以改善低信噪比下信号在时频域上的稀疏性和时频分布准确性，处理后的时频域信号系数表示为B(t，f)；过程为：

设优化后的时频矩阵表示为B(t，f)∈RN^×P，为了减小计算复杂度，这里设B(t，f)为二进制矩阵，即矩阵中的非零元素都为1；为了实现求解过程中的变量选择，这里考虑B(t，f)两个特性：

1)、时频点稀疏性：B(t，f)中的大部分元素为0，仅存在信号的时频点b(t，f)∈B(t，f)对应的时频系数非零，以保证B(t，f)在时频域上是稀疏的；

2)、双差分稀疏性：根据跳频信号的时频分布特性，由于跳频信号的每一跳具有一定的持续时间，所以如果不发生频率跳变，则B(t，f)中相邻列的时频系数相同，同时考虑到相邻列的时频系数的平滑性，则有2b(t，f)-b(t-1，f)-b(t+1，f)＝0。

根据上述两个特性的限制，目标函数写为

其中，F⊙B表示矩阵的Hadamard积，即两个矩阵的同位置元素相乘，D∈R^(P-2)×P，且

在(4)式中，第一项考虑了信号优化前后的误差，而λ₁和λ₂则分别控制了估计结果的稀疏性和平滑性，然而(4)式中的目标函数是非凸的，属于NP-hard问题，因此，为了使(4)式更容易直接最小化，这里将l₀范数转化为属于凸函数的l₁范数，同时加入可微的l₂范数以简化迭代求解过程，则有

其中，λ₁越大，B(t，f)的稀疏性越好，非零系数越少；λ₂越大，B(t，f)内同行系数的平滑性越好，可以看到(5)式是严格的凸函数，因此具有唯一解。

为了简化分析，将B(t，f)的优化问题分解为逐行的单独求解，这里设F(t，f)＝[f₁，f₂，...，f_N]^T，f_i∈R^1×P，B(t，f)＝[b₁，b₂，...，b_N]^T，b_i∈R^1×P，则有

其中，w_i＝diag(f_i)，即w_i为对角矩阵，对角线上的元素为f_i；

利用L2-ISTA算法对(7)式进行快速迭代求解，首先建立带约束条件的一般化模型：

min{F(b)≡f(b)+g(b)：b∈Rⁿ} (7)

其中，g(b)＝λ₁||b||₁，

对f(b)进行梯度下降来求解，则有

其中，L为下降步长，k代表迭代次数；因为f(b)属于凸函数，所以

有上界，即

||f(b_k)-f(b_k-1)||≤L(f)||b_k-b_k-1|| (9)

其中，||·||表示标准欧几里得范数，L(f)＞0；为了防止梯度下降过程中参数更新变化太大，降低梯度爆炸的发生概率，所以L＝min(L(f))；基于近端梯度法的思想，对f(b)在α处的泰勒展开为

且有

λ₀≥max[eig(w^Tw+2λ₂D^TD)] (11)

其中，max[eig(w^Tw+λ₂D^TD)]表示w^Tw+λ₂D^TD的最大特征值，φ(α)是与b无关的函数项，忽略不计；因此将(7)式改写为

其中，设λ₀＝max[eig(w^Tw+2λ₂D^TD)]；

观察(12)式可知，λ₀等价为▽f的Lipschitz常数，即λ₀＝L；当α＝b_k-1时，则有

b_k＝P(b_k-1) (13)

设

且

若H(b_k)取得最小值，则有

其中sgn(·)为符号函数，基于软阈值算法解得

其中，

表示

的第i个元素，z_i表示z的第i个元素；

接下来给出关于(λ₁,λ₂)的选择准则：

准则1、当λ₂＝0时，λ₁的有效取值范围为

且

当λ₂＝0时，(6)式则转化为l₁惩罚回归(Lasso)模型；当

时，则b_i→0，即b_i＝0为(6)式的求解结果；λ₁建议取值为

的5％-10％，此时(6)式获得良好的估计结果；

准则2、当λ₁＝0时，λ₂的有效取值范围为

且

当λ₁＝0时，则有

其中，

是非零元素全为1的下三角矩阵，M₀是M的第一列系数；当

时，则有b_i→c_i，其中c_i为常数，即常数向量为(6)式的求解结果；当λ₁≠0时，λ₂的建议取值为

的5％-10％，此时(6)式能够获得良好的估计结果；

S3、利用基于二次包络优化的参数估计算法，首先利用数据的连续性分析对步骤S2中得到的时频域信号B(t,f)进行去除部分干扰信号特征和异常点处理，然后提取优化后时频域信号B(t,f)的平均时频脊线，进行平滑处理，以减小剩余干扰对跳频信号的影响；过程为：

通过分析B(t,f)中每行连续分布的非零系数对干扰源进行针对性去除，已知B(t,f)＝[b₁,b₂,...,b_N]^T，b_i∈R^1×P，对b_i＝[b_i1,b_i2,...,b_iP]内连续的非零值(非零值全为1)进行数目和位置的统计，分别表示为{c₁,c₂,...,c_l}和{loc₁,loc₂,...,loc_l}，其中c_j表示第j段连续非零值的数目，loc_j表示其位置，如图1所示；

根据跳频信号的时频分布特性，建立判别准则如下：

准则3、为去除宽带信号和部分异常点，当c_w≤th₁，1≤w≤l时，使对应loc_w上的时频系数全为0，即

一般th₁＝η₁P，η₁≤1/50；

准则4、针对窄带信号中的定频干扰信号，根据其能量分布集中、连续的特点，当c_w≥th₂时，使对应loc_w上的时频系数全为0，即

一般th₂＝η₂P，η₂≥1/2；

通过上述的判别准去除时频域上大部分的干扰信号，包括宽带信号、窄带定频干扰信号以及信号因发生畸变而产生的边缘特征等；B(t,f)经过准则3、准则4处理后，得到的结果表示为B₁(t,f)；

当存在频率调制类型的窄带信号时，因其与跳频信号的时频分布具有较大的相似性而无法被有效限制；此时，当该类窄带信号与跳频信号在部分频点上存在混叠时，由于窄带信号的能量分布较为集中，使得这些频点上的时频分布很难被辨别，从而极大概率会恶化后续的时频脊线提取性能，使得参数估计发生严重错误；因此，考虑到这种极端的混叠情况，提出基于平均时频脊线的包络优化方法以突破传统方法下信号的稀疏性限制。

B₁(t,f)内主要由跳频信号特征和频率调制类信号的部分特征构成，其内部的非零值全为1；提取跳频信号的时频特征，简单地对每一列的第一个或者最后一个非零值进行选取可能会直接使得时频分布在同一个频点上存在巨大差异或者不同频点之间存在微小差异，进而破坏脊线的周期性。也可以通过获取整段信号在时频域上分布的先验知识来改善特征提取的准确性，但工作量非常大。因此，为了避免非零值提取过程中的随机性，这里提取全局的平均时频脊线为基础对参数进行估计。

首先，对B₁(t,f)∈R^N×P每一列的非零值位置进行提取，记为

然后保留ind_R/2处的非零值；由于每一列只取一个非零值的位置，此时基本上可以保证相邻跳的能量分布不在同一个频率点上；对每一列ind_R/2处的非零值进行重排，提取到的脊线称为平均时频脊线，记为mean_tfcurve，其时、频域上分别对应时间和暂态频率；根据跳频信号的时频分布特性可知，平均时频脊线在一定程度上减小了同一跳内存在的频点差异，同时扩大了不同跳之间的频点差异。

此时mean_tfcurve存在较为明显的跳周期，但仍然存在大量的跳变数据点和异常点，跳时刻检测的精度较低。虽然跳变数据点具有较强的聚集性，但在mean_tfcurve中是非连续的，这里将其表示为W＝[W₁,W₂,...,W_p]，其中W_i表示单个聚集的跳变数据点集合，如图2所示。

首先对异常点进行处理：检测mean_tfcurve中相邻点的差异分布，有

M_d＝DM(mean_tfcurve)＝max[Diff(mean_tfcurve)] (18)

其中，DM(·)用于提取目标数据的差分曲线最大值，max(·)表示取最大值操作，Diff(·)表示差分处理；

设平滑处理函数为Sm(·)，且有X＝[x₁,x₂,...,x_n]，则Sm(X)的处理结果为

其中，α为平滑处理阈值，一般为DM(X)的0.01；设SL＝Sm(mean_tfcurve)，α＝0.01M_d；由于SL中依然存在多个跳变数据点集合W_i和异常点，使得脊线的连续性较差，因此需要对SL采取进一步的优化处理。

根据SL的曲线特性，同时也便于后续的估计工作，这里采用分段线性插值法对SL的上、下包络曲线U₁和L₁进行提取，如图3所示。

从图3中可以看出，U₁和L₁对SL的分布进行了良好的描述，且将多个W_i分别进行了连接，改善了SL的完整性和连续性；为验证包络曲线的提取效果和W_i的连接效果，对U₁和L₁的拐点进行检测；设拐点检测函数为DI(·)，则有

DI(U₁)＝Diff(Diff(U₁))

DI(L₁)＝Diff(Diff(L₁)) (20)

其中，UI₁和LI₁分别表示U₁和L₁中检测到的有效拐点，拐点的检测阈值分别为β·max[DI(U₁)]和β·max[DI(L₁)]，一般β≤0.01；UI₁和LI₁分布如图4所示。

为了获得更准确的拐点分布，进一步对U₁和L₁进行优化处理。首先提取W_i的位置，当SL中的数据点数目较多时，需要自动定位W_i的分布，然后对相应的取值进行平滑处理。根据图3可以发现，U₁和L₁有效提取了SL中的平行部分，而这些平行分布的数据段中包含了W＝[W₁,W₂,...,W_p]，因此对SL和U₁、L₁的重叠部分进行提取，有

W_u_loc＝Loc(SL∩U₁)

W_l_loc＝Loc(SL∩L₁) (22)

其中，W_u_loc和W_l_loc分别表示SL和U₁、SL和L₁的重叠部分的时域上的索引，Loc(·)为目标对象索引的提取函数，∩表示取交集。

由于W_i在U₁和L₁中是连续的，所以接下来分别对U₁中W_u_loc上的值和L₁中W_l_loc上的值进行平滑处理，则有

U₁_cs＝Sm[U₁(W_u_loc)]

L₁_cs＝Sm[L₁(W_l_loc)] (23)

其中，U₁_cs和L₁_cs分别表示对U₁(W_u_loc)和L₁(W_l_loc)进行平滑处理后的数据集合，参数估计依然要回归到SL的优化问题，因此对SL重新填充，得到SL_cs，即

由于重新填充后SL在W_u_loc和W_l_loc的周边可能会出现跳变值，因此依然对SL_cs进行平滑处理，完成最后的平均时频脊线优化，即

SL_cs＝Sm(SL_cs) (25)；

S4、对步骤S3中得到的平均时频脊线进行拐点检测，然后对检测到的拐点进行时频域映射，完成跳频信号的跳周期估计；过程为：

由于W_i在SL_cs中仍然是非连续的，且数量较多，因此再次对SL_cs进行包络曲线的提取以连接W_i，这里表示为U₂和L₂；对U₂和L₂的拐点进行检测，表示为UI₂和LI₂，如图5所示。

对比图4可以发现，多个W_i内的数据点得到了有效的平滑处理，具体变现为多余的拐点被去除，使得基本上有且只有一对拐点与W_i对应；因此可以通过提取UI₂和LI₂中平行分布的数据点，即成对平行于时间轴的拐点，并结合SL_cs完成对跳时刻

和跳周期

的高精度估计；

S5、由于平均时频脊线中的脊线部分对应的为暂态频率，而非真实的跳频频点，因此需要额外设计跳频频点的估计方法。

对B₁(t,f)内单个跳周期

内的时频系数进行提取，表示为g_j(t,f)∈R^N×h，其中h为

内所含列的数目。由于B₁(t,f)内仍有残余的频率调制信号的信号特征，因此g_j(t,f)内也可能会含有除跳频信号以外的时频分布。如图6所示，其中，上部分为频率调制信号的信号特征，下部分为跳频信号的信号特征，两者都表现为不同的线性特征。

由于信号的时频分布经过了之前的迭代和准则处理，因此相比于干扰信号，跳频信号会产生聚集度更高的线性特征。为了实现对应频点的自动估计，提出基于Hough变换的优化方法来检测跳频信号的信号特征，具体过程为：提取g_j(t,f)内非零值元素的位置，坐标位置表示为

并做极坐标转换：

当[x_z,y_z]与[x_z+i,y_z+i]分布在同一条直线l上时，对应的p_z与p_z+i会在θ＝θ_z处相交，其中θ_z为l和坐标轴的夹角；对

进行计算，其分布如图7所示；对

的交点位置和出现次数进行统计，并记出现次数最多的交点为

如图8所示；根据

筛选出相交于该点的曲线{p_z}，即可获得对应的{[x_z,y_z]}；保留g_j(t,f)中{[x_z,y_z]}上的非零值，其他系数全置为0，并重新记

内的时频系数为

如图9所示，图9中，(a)、LFM源信号,(b)、EQFM源信号,(c)、BPSK源信号,(d)、4FSK源信号,(e)、FH源信号,(f)、混合信号；提取

内所有非零值元素的纵轴值，并进行平均处理，即

其中

为

内的频点估计值，f_i表示

内非零值元素的纵轴值，共S个；对B₁(t,f)内不同跳周期内的g_j(t,f)进行上述处理，即可完成接收信号内所有跳频频点的估计。

对SNR＝0dB下的混叠信号f(t)进行处理，首先给出其时频分布F(t,f)和频谱，如图10所示；利用基于稀疏线性回归的矩阵优化算法对F(t,f)优化后为B(t,f)，其中

如图11所示；可以发现，F(t,f)中的噪声基本被去除，且主要信号特征被有效提取，时频分布精度得到了大幅改善。

对B(t,f)进行二次包络优化处理，然后提取其平均时频脊线，其中η₁＝1/100，η₂＝1/5，得到处理后的B₁(t,f)，如图12所示；最终优化后的平均时频脊线的拐点分布UI₂和LI₂，如图13所示。

舍弃数据段末尾部分的冗余信号特征，根据B₁(t,f)与平均时频脊线上的拐点分布，可以得到跳时刻的估计值集合

和跳周期的估计值集合

接下来通过

对B₁(t,f)进行分段，得到{g_j(t,f)}；利用基于Hough变换的优化方法对g_j(t,f)进行处理得到

并自动得到跳频频点的估计值集合

图14分别给出了{t_m}和{f_m}的估计误差，这里同时给出

拼接后的时频分布，如图15所示；对比图10可以看出，干扰信号的信号特征被有效去除，

中只保留了跳频信号的信号特征；可以得到跳时刻的估计误差

dB，跳周期的估计误差

跳频频点的估计误差

综上所述，在低信噪比下，本发明的方法对混叠信号中跳频信号的参数估计具有良好的性能和鲁棒性。

为对比现有技术中的其他主流传统算法，选取了两种主流的单通道跳频信号参数估计方法与本发明时频混叠跳频信号的参数估计方法进行性能对比，分别是基于SPWVD的滤波方法和基于STFT的二次迭代稀疏重构方法。

基于SPWVD的方法侧重于提高全局的时频分辨率，而基于二次迭代稀疏重构方法则更好地去除了噪声和信号的冗余、畸变特征。在此基础上，两者继续对信号进行形态学滤波，以尽可能地去除干扰信号对目标信号的影响，然后结合不同的优化算法对参数做出估计。其中形态学滤波过程复杂繁琐，自适应性较差。

分别在不同信噪比下对算法进行100次蒙特卡洛实验，图16给出各方法的参数估计误差NMSE，图16中，(a)、跳时刻的估计误差

曲线对比，(b)、跳频频点的估计误差

曲线对比。

实验结果表明，本发明时频混叠跳频信号的参数估计方法，其在应用范围和参数估计性能方面都要优于传统方法，且在低信噪比下具有良好的鲁棒性。

以上所述仅为本发明的较佳实施例，而非对本发明的限制，在不脱离本发明的精神和范围的情况下，凡依本发明申请专利范围所作的均等变化与修饰，皆应属本发明的专利保护范围之内。

Claims (7)

1.一种时频混叠跳频信号的参数估计方法，其特征是：其包括以下步骤：

S1、将单通道接收下的时频混叠信号f(t)，利用短时傅里叶变换STFT转换到时频域，得到的时频域信号F(t,f)；其中，t表示时间，f表示频率；

S2、利用基于稀疏线性回归的矩阵优化算法对步骤S1中得到的时频域信号F(t,f)进行去除背景噪声、信号畸变特征且保留信号的主要特征处理，以改善低信噪比下信号在时频域上的稀疏性和时频分布准确性，处理后的时频域信号表示为B(t,f)；

S3、利用基于二次包络优化的参数估计算法，首先利用数据的连续性分析对步骤S2中得到的时频域信号B(t,f)进行去除部分干扰信号特征和异常点处理，然后提取优化后时频域信号B(t,f)的平均时频脊线，进行平滑处理，以减小剩余干扰对跳频信号的影响；

S4、对步骤S3中得到的平均时频脊线进行拐点检测，然后对检测到的拐点进行时频域映射，完成跳频信号的跳频周期估计；

S5、最后在步骤S4中得到的跳频周期估计的基础上，基于Hough变换的优化方法完成跳频信号的频点估计。

2.根据权利要求1所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S1中，对单通道接收下的时频混叠信号f(t)进行短时傅里叶变换，过程为：

其中，f(t)表示观测信号，A∈R^1×n为一维的观测矩阵，其中，R是实数，n是正整数，a_n表示第n个源信号的能量系数，s_n(t)表示第n个源信号，T是信号f(t)的观测时长，s(t)表示源信号，v(t)表示加性噪声；

在观测时间[0,T]内，设跳频信号s_i(t)包含M个跳频周期，即

其中，m表示跳频信号的序号，K_m∈R和f_m∈[-f_max,f_max]分别为第m跳的幅度和频率，跳频周期为T_hop＝t_m-t_m-1，t₀＝0；

将其转换到时频域来改善信号的稀疏性，表示为：

其中，F(t,f),S(t,f)和V(t,f)分别为f(t),s(t)和v(t)在短时傅里叶变换STFT下的时频系数。

3.根据权利要求1所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S2中，利用基于稀疏线性回归的矩阵优化算法对步骤S1中得到的时频域信号F(t,f)进行处理，过程为：

设优化后的时频矩阵表示为R(t,f)∈R^N×P，其中，N是矩阵的行数目，P是矩阵的列数目，B(t,f)为二进制矩阵，即矩阵中的非零元素都为1；B(t,f)具有时频点稀疏性和双差分稀疏性两个特性，根据两个特性的限制，目标函数写为

其中，F⊙B表示矩阵的Hadamard积，即两个矩阵的同位置元素相乘，D∈R^(P-2)×P，且

在(4)式中，λ₁和λ₂为可变的、大于0的实数，对应的取值分别控制了估计结果的稀疏性和平滑性，将(4)式中非凸的目标函数

范数转化为属于凸函数的

范数，同时加入可微的

范数，则有

其中，λ₁越大，B(t,f)的稀疏性越好，非零系数越少；λ₂越大，B(t,f)内同行系数的平滑性越好；

设F(t,f)＝[f₁,f₂,...,f_N]^T,f_i∈R^1×P，f_N表示矩阵的行数目；B(t,f)＝[b₁,b₂,...,b_N]^T,b_N表示第N行的矩阵信号,b_i∈R^1×P，P是矩阵的列数目，则有

其中，w_i＝diag(f_i)，即w_i为对角矩阵，对角线上的元素为f_i；

用L2范数-迭代阈值收缩优化算法对(7)式进行快速迭代求解，首先建立带约束条件的一般化模型：

min{F(b)≡f(b)+g(b):b∈Rⁿ} (7)

其中，g(b)＝λ₁||b||₁，

对f(b)进行梯度下降来求解，则有

其中，L为下降步长，k代表迭代次数；因为f(b)属于凸函数，所以

有上界，即

||f(b_k)-f(b_k-1)||≤L(f)||b_k-b_k-1|| (9)

其中，||·||表示标准欧几里得范数，L(f)>0；L＝min(L(f))；基于近端梯度法的思想，对f(b)在α处的泰勒展开为

且有

λ₀≥max[eig(w^Tw+2λ₂D^TD)] (11)

其中，max[eig(w^Tw+λ₂D^TD)]表示w^Tw+λ₂D^TD的最大特征值，φ(α)是与b无关的函数项，忽略不计；因此将(7)式改写为

其中，设λ₀＝max[eig(w^Tw+2λ₂D^TD)]；

由(12)式可知，λ₀等价为

的Lipschitz常数，即λ₀＝L；当α＝b_k-1时，则有

b_k＝P(b_k-1) (13)

设

且

若H(b_k)取得最小值，则有

其中，sgn(·)为符号函数，基于软阈值算法解得

其中，

表示

的第i个元素，z_i表示z的第i个元素。

4.根据权利要求3所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S2中的(λ₁,λ₂)具有两个选择准则：

准则1、当λ₂＝0时，λ₁的有效取值范围为

且

当λ₂＝0时，(6)式则转化为

惩罚回归模型；当

时，则b_i→0，即b_i＝0为(6)式的求解结果；λ₁取值为

的5％-10％，此时(6)式能够获得良好的估计结果；

准则2、当λ₁＝0时，λ₂的有效取值范围为

且

当λ₁＝0时，则有

其中，

是非零元素全为1的下三角矩阵，M₀是M的第一列系数；当

时，则有b_i→c_i，其中c_i为常数，即常数向量为(6)式的求解结果；当λ₁≠0时，λ₂取值为

的5％-10％，此时(6)式能够获得良好的估计结果。

5.根据权利要求3所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S3中，利用基于二次包络优化的参数估计算法对步骤S2中得到的时频域信号B(t,f)进行处理，过程为：

已知B(t,f)＝[b₁,b₂,...,b_N]^T，b_i∈R^1×P，对b_i＝[b_i1,b_i2,...,b_iP]内连续的非零值进行数目和位置的统计，分别表示为{c₁,c₂,...,c_l}和{loc₁,loc₂,...,loc_l}，其中，c_j表示第j段连续非零值的数目，loc_j表示其位置；

根据跳频信号的时频分布特性，建立判别准则如下：

准则3、去除宽带信号和部分异常点：当c_w≤th₁，1≤w≤l时，使对应loc_w上的时频系数全为0，即时频系数

c_w表示第w段连续非零值的数目，loc_w表示其位置，th₁＝η₁P，η₁≤1/50；

准则4、针对窄带信号中的定频干扰信号：当c_w≥th₂时，使对应loc_w上的时频系数全为0，即时频系数

c_w表示第w段连续非零值的数目，loc_w表示其位置，th₂＝η₂P，η₂≥1/2；

通过准则3、准则4去除时频域上大部分的干扰信号，包括宽带信号、窄带定频干扰信号以及信号因发生畸变而产生的边缘特征；B(t,f)经过准则3、准则4处理后，得到的结果表示为B₁(t,f)；

首先，对B₁(t,f)∈R^N×P每一列的非零值位置进行提取，此位置记为

然后保留ind_R/2处的非零值，即B₁(t,f)∈R^N×P每一列的非零值中间位置处的非零值；对每一列ind_R/2处的非零值进行重排，提取到的脊线称为平均时频脊线，记为mean_tfcurve，其时、频域上分别对应时间和暂态频率；

首先对异常点进行处理：检测mean_tfcurve中相邻点的差异分布，有

M_d＝DM(mean_tfcurve)＝max[Diff(mean_tfcurve)] (18)

其中，DM(·)用于提取目标数据的差分曲线最大值，max(·)表示取最大值操作，Diff(·)表示差分处理；

设平滑处理函数为Sm(·)，且有X＝[2₁,x₂,...,x_n]，则Sm(X)的处理结果为

其中，α为平滑处理阈值，为DM(X)的0.01，设SL＝Sm(mean_tfcurve)，α＝0.01M_d；

采用分段线性插值法对SL的上、下包络曲线U₁和L₁进行提取，过程为：

对U₁和L₁的拐点进行检测，设拐点检测函数为DI(·)，则有

DI(U₁)＝Diff(Diff(U₁))

DI(L₁)＝Diff(Diff(L₁)) (20)

其中，UI₁和LI₁分别表示U₁和L₁中检测到的有效拐点，拐点的检测阈值分别为β·max[DI(U₁)]和β·max[DI(L₁)]，β≤0.01；

再进一步对U₁和L₁进行优化处理，过程为：U₁和L₁有效提取了SL中的平行部分，这些平行分布的数据段中包含了W＝[W₁,W₂,...,W_p]，因此对SL和U₁、L₁的重叠部分进行提取，有

W_u_loc＝Loc(SL∩U₁)

W_l_loc＝Loc(SL∩L₁) (22)

其中，W_u_loc和W_l_loc分别表示SL和U₁、SL和L₁的重叠部分的时域上的索引，Loc(·)为目标对象索引的提取函数，∩表示取交集；

分别对U₁中W_u_loc上的值和L₁中W_l_loc上的值进行平滑处理，则有

U₁_cs＝Sm[U₁(W_u_loc)]

L₁_cs＝Sm[L₁(W_l_loc)] (23)

其中，U₁_cs和L₁_cs分别表示对U₁(W_u_loc)和L₁(W_l_loc)进行平滑处理后的数据集合；对SL重新填充，得到SL_cs，即

对SL_cs进行平滑处理，完成最后的平均时频脊线优化，即

SL_cs＝Sm(SL_cs) (25)。

6.根据权利要求5所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S4中，对步骤S3中得到的平均时频脊线进行拐点检测完成跳频信号的跳频周期估计，过程为：

再次对SL_cs进行包络曲线的提取，表示为U₂和L₂；对U₂和L₂的拐点进行检测，表示为UI₂和LI₂；通过提取UI₂和LI₂中平行分布的数据点，即成对平行于时间轴的拐点，并结合SL_cs完成对跳频时刻

和跳频周期

的高精度估计。

7.根据权利要求5所述的时频混叠跳频信号的参数估计方法，其特征是：其步骤S5中，利用基于Hough变换的优化方法完成跳频信号的频点估计，过程为：对B₁(t,f)内单个跳频周期

内的时频系数进行提取，表示为g_j(t,f)∈R^N×h，其中，h为

内所含列的数目；提取g_j(t,f)内非零值元素的位置，坐标位置表示为

Z表示非零值元素的数目，并做极坐标转换：

当[x_z,y_z]与[x_z+i,y_z+i]分布在同一条直线l上时，对应的p_z与p_z+i会在θ＝θ_z处相交，其中θ_z为l和坐标轴的夹角；对

的交点位置和出现次数进行统计，并记出现次数最多的交点为

根据

筛选出相交于该点的曲线{p_z}，即可获得对应的{[x_z,y_z]}，保留g_j(t,f)中{[x_z,y_z]}上的非零值，其他系数全置为0，并重新记

内的时频系数为

提取

内所有非零值元素的纵轴值，并进行平均处理，即

其中，

为

内的频点估计值，f_i表示

内非零值元素的纵轴值，共S个；对B₁(t,f)内不同跳频周期内的g_j(t,f)进行处理，即完成接收信号内所有跳频频点的估计。

Images