CN106778610A - 一种基于时频图像特征的脉内调制识别方法 - Google Patents

Abstract

本发明属于雷达辐射源信号识别技术领域，具体涉及一种基于时频图像特征的脉内调制识别方法。本发明包括：(1)对雷达信号进行Wigner‑Vill分布时频变换，得到反映信号能量随时间和频率的时频分布图；(2)对步骤(1)得到的时频图像进行图像预处理；(3)对预处理后的图像进行插值打散原轨迹；(4)将步骤(3)中被打散成单个点的轨迹类比成目标跟踪中的目标位置，并引入联合概率数据关联来识别出不同雷达信号的轨迹；(5)运用中心矩提取时频图像的形状特征等。本发明提出了一种将雷达信号通过时频分布变换转换为时频图像，通过数字图像处理将同一时刻到达并且频率互相交叠的情况下，分别识别出雷达信号调制方式的方法。

Description

一种基于时频图像特征的脉内调制识别方法

技术领域

本发明属于雷达辐射源信号识别技术领域，具体涉及一种基于时频图像特征的脉内调制识别方法。

背景技术

雷达辐射源信号识别是电子情报侦察、电子支援侦察和威胁告警系统中的关键处理过程，其识别水平直接关系到雷达对抗设备的技术先进程度。随着现代电子战的激烈对抗，新型复杂体制雷达不断投入使用并逐渐占据主导地位，电磁威胁环境的信号密度高达每秒120万个脉冲以上，雷达工作频率覆盖范围达到了0.1～20GHz，并正在向0.05～140GHz扩展，雷达信号波形在时频等多个域中同时变化，隐身和抗干扰能力大大增强。传统五参数(载频、脉冲到达时间、脉冲幅度、脉冲宽度和脉冲到达方向)识别方法难以适应这样密集、复杂和多变的信号环境，雷达辐射源信号识别面临着前所未有的挑战。近年来关于雷达辐射源信号识别的大量研究成果表明，脉内调制识别将有望在新型复杂体制雷达辐射源信号识别技术与装备方面获得重要突破。现有的雷达辐射源信号脉内调制识别方法主要有小波及小波包变换法、小波脊线法、相像系数法、经验模式分解法、相对无模糊相位重构法、熵特征法、围线积分双谱法、复杂度特征法、分数傅里叶-包络法、分形维数法、模糊函数主脊切面法、瞬时频率派生特征法和模糊函数主脊切面特征法等。现有的方法为使信号之间差异明显，都选择将采样信号做某种变换，对变换后的信号进行特征提取，进而完成对信号的调制方式的识别。这些方法都主要针对少数几种信号，对于其它信号的识别效果还有待于进一步研究，同时这些方法都是只考虑了同一时刻只有单一信号的情况，没有探讨同一时刻多信号频率交叠的情况。

发明内容

本发明的目的是提出一种既能够对不同调制类型的雷达信号、又能够利用时频图像区分同一时刻多信号频率交叠的基于时频图像特征的脉内调制识别方法。

本发明的目的是这样实现的：

(1)对雷达信号进行Wigner-Vill分布时频变换，得到反映信号能量随时间和频率的时频分布图；

对相关函数作滑窗处理，得到时变的局部相关函数公式，

当窗函数取时间冲击函数，不加限制，而在时域取瞬时值，

对时变局部相关函数作Fourier变换，即可得到WVD，

设有n个分量信号，得到多分量信号的WVD，

雷达信号经过Wigner-Vill分布时频变换后得到时频图像；

(2)对步骤(1)得到的时频图像进行图像预处理；

对时频图像灰度值进行归一化；设和σ²分别是所有像素点灰度值的均值和方差；

归一化后像素点的灰度值为：

采用自适应维纳滤波器初步去除时频图像的噪声点，对图像进行增强；采用形态学方法对灰度图像进行处理；形态学处理是应用具有一定形态的结构元素对灰度图像进行腐蚀和膨胀的操作，其中使用结构元素b对f的灰度膨胀记为

选择半径为3的菱型结构元素对时频图像进行闭运算；

通过检测信号时频分布的起止频率，将没有信号分布的图像区域剪切掉；

采用最近邻域插值法归一化时频图像的纵横比；

将原图定义为一个二维函数F(x,y)，(x,y)点的函数值对应像素点的灰度值，f(x,y)为归一化纵横比后的目标图；

Fw、Fh设为原图的宽度和高度；fw、fh设为目标图的宽度和高度；原图坐标(x,y)和目标图坐标(x′,y′)之间公式为：

x＝x′(Fw/fw)

y＝y′(Fh/fh)；

(3)对预处理后的图像进行插值打散原轨迹；

(4)将步骤(3)中被打散成单个点的轨迹类比成目标跟踪中的目标位置，并引入联合概率数据关联来识别出不同雷达信号的轨迹；

引入确认矩阵；确认矩阵定义为：

Ω＝(ω_jt),j＝1,2,...,m(k),t＝0,1,...,n

其中ω_jt表示量测j是否落入目标t的确认门之内；t＝0表示没有目标，此时对应的Ω的列的元素都是1；

设在杂波环境中已有T个目标，则他们的状态方程和测量方程分别表示为：

JPDA的状态更新方程为：

其中，表示k时刻目标t的状态矢量，表示目标t的状态预测矢量，V′(k)为联合新息：

其中，β_jt表示目标t与观测j的关联概率，V_j′(k)表示滤波新息；

关联概率表示为：

其中，表示在可行事件θ^j中，观测j是否源于目标t；如果源于目标则等于1，否则为0；在k时刻联合事件θ的条件概率表示为

其中，σ_t为目标检测指示器，如果目标在θ中与观测关联，则σ_t的值等于1，反之为零；Φ是虚假观测事件的个数，表示目标t的检测概率，V表示航迹有效门体积；

(5)运用中心矩提取时频图像的形状特征；

设大小为M*N的二值图像表示为：f(x,y)∈{1,0}，则图像的(p+q)阶原点矩定义为：

图像的(p+q)阶中心矩定义为：

其中表示水平方向上的质心；表示垂直方向的质心；

其中μ₀₂表示图像在垂直方向上的伸展度；μ₂₀表示图像在水平方向上的伸展度；μ₁₁表示图像的倾斜度；μ₀₃表示图像在垂直方向上的重心偏移度；μ₁₂表示图像垂直伸展的均衡程度；由于时频图的差异性主要体现在垂直方向的频域上，故舍弃μ₀₂，而构造出特征向量；

(6)选取支持向量机分类器的核函数，并根据训练样本的形态特征确定核函数的参数，以训练支持向量机分类器；

其中，支持向量机分类器的核函数为高斯核函数，高斯核函数为：

k(x₁,x₂)＝exp(-γ*|x₁-x₂|²)，

选择任意两类设计1个二分类SVM，总共设计n(n-1)/2个二分类SVM；然后将这个二分类SVM构成一个有向无循环图，该图有n(n-1)/2个非叶子节点和n个叶子节点，每个非叶子节点表示一个二分类SVM分类器，并与下一层的两个节点相连，每个叶子节点表示一个输出；

具体而言，选择上述的高斯核函数作为核函数，核函数参数的选择采用交叉验证的方法；其中，惩罚因子C从2^-5成倍增加到2⁵，核函数参数g从2^-5成倍增加到2⁵，对于(C,g)，选择5折交叉验证，即将训练样本分成5份，挑选1份作为测试样本，剩下4份作为训练样本，重复5次，计算平均正确率，以此衡量当前(C,g)的性能，选择性能相对最好的(C,g)作为最终的参数；最后，将步骤(5)得到的特征向量输入训练好的支持向量机。

本发明的有益效果在于：提出了一种将雷达信号通过时频分布变换转换为时频图像，通过数字图像处理将同一时刻到达并且频率互相交叠的情况下，分别识别出雷达信号调制方式的方法。

附图说明

图1是本发明雷达信号调制识别的流程图。

图2是本发明采用的支持向量机分类器示意图。

具体实施方式

下面结合附图对本发明做进一步描述。

首先将待识别的雷达信号进行Wigner-Vill分布时频变换，得到时频分布图。时频分析反映了信号能量随时间和频率的分布，在时频域精确地描述了信号。通过数字图像处理方法从信号的时频分布中提取出有效的识别特征，并将特征向量输入到分类器完成分类识别任务。

本发明的具体步骤是：

1)对雷达信号进行Wigner-Vill分布时频变换，得到反映信号能量随时间和频率的时频分布图。

2)对时频分布图进行预处理，抑制噪声，将图形分割为单独的点。各个信号时频图像灰度值的动态范围是不一样的，像素点灰度值大的数据对分类识别有着大的影响，为了减少数据间的不平衡性，首先对时频图像灰度值进行归一化，然后使用自适应维纳滤波器去除时频图像的噪声点，对图像进行增强。为了减少后期图像处理的计算量和存储空间，对灰度图像进行二值化，将灰度图转化为二值图，二值化过程中阈值的选取采用一维最大熵法。接着对二值图像进行闭运算(膨胀之后再进行腐蚀运算)，使得时频图像中信号分量的轮廓变得光滑，进一步减少噪声。从时频图中可以看出，并非所有区域都分布有信号，我们可以将没有信号分布的图像区域剪切掉，减小冗余信息，更有利于下一步信号特征的提取。最后采用最近邻插值法归一化时频图像的纵横比，使所有信号的时频图像大小都保持一致，并进一步减小数据量。

3)对步骤2预处理后的图像进行插值，将原来连成线的信号时频分布形状打散成单个一系列点迹。

4)通过联合概率数据关联算法对步骤3得到的图像进行轨迹识别，识别出几个轨迹代表可能包含几个不同的雷达信号，并将不同轨迹的信号分离。

5)采用图像处理的算法和工具对时频图像进行处理，提取时频图像的形状特征，本发明采用中心矩来提取形状特征。

6)采用支持向量机分类器识别雷达信号脉内调制方式。支持向量机是一种基于结构风险最小化的分类方法，通过核函数的非线性变换，将待分类样本从样本控件映射到高维特征控件，然后在特征控件中最大化分类间隔，确定最优分类超平面。从每种待识别的雷达信号中抽取一定的数量的样本，对其进行信号特征提取，并输入到对应的二分类SVM分类器中，通过计算，获得该SVM的参数值，直至最终确认整个分类器。将特征向量输入到支持向量机完成基于时频图像特征的脉内调制识别方法。

步骤1：对雷达信号进行Wigner-Vill分布时频变换，得到反映信号能量随时间和频率的时频分布图。

为了体现非平稳信号的局部时变特性，对相关函数作滑窗处理，得到时变的局部相关函数公式(1)，

当窗函数取时间冲击函数，对不加限制，而在时域取瞬时值，此时如公式(2)

对时变局部相关函数作Fourier变换，即可得到WVD，其表达式如公式(3)

WVD的时间带宽积达到Heisenberg不确定原理给出的下界，因此具有最理想的时频分辨率。设有n个分量信号，可以得到多分量信号的WVD，如公式(4)

雷达信号经过Wigner-Vill分布时频变换后得到时频图像。

步骤2：对步骤1得到的时频图像进行图像预处理。

各个信号的时频图像灰度值的动态范围是不一样的，像素点灰度值大的数据对分类识别有着大的影响，为了减少数据间的不平衡性，需要对时频图像灰度值进行归一化。设和σ²分别是所有像素点灰度值的均值和方差。

归一化后像素点的灰度值为：

接下来采用自适应维纳滤波器初步去除时频图像的噪声点，对图像进行增强。为了使得时频图像的信号分量轮廓变的光滑，降低噪声和交叉项时频点的影响，进一步采用形态学方法对灰度图像进行处理。形态学处理是应用具有一定形态的结构元素对灰度图像进行腐蚀和膨胀的操作，其中使用结构元素b对f的灰度膨胀记为

开运算可以去除比结构元素更小的时频独立点和细小的突出部分，而闭运算一般将狭窄的缺口连接起来，填充时频区域之间小的缝隙。考虑到二值化后信号项的时频区域会出现断裂的缝隙，本发明选择半径为3的菱型结构元素对时频图像进行闭运算。

时频图像剪切通过检测信号时频分布的起止频率，将没有信号分布的图像区域剪切掉，减小冗余信息，更有利于下一步信号特征的提取。

为了使所有信号的时频图像大小都保持一致，减小数据量，采用最近邻域插值法归一化时频图像的纵横比。将原图定义为一个二维函数F(x,y)，(x,y)点的函数值对应像素点的灰度值，f(x,y)为归一化纵横比后的目标图。

Fw、Fh设为原图的宽度和高度。fw、fh设为目标图的宽度和高度。那么原图坐标(x,y)和目标图坐标(x′,y′)之间公式为：

坐标转化过程中，目标图的坐标值可能会出现小数，对此在原图像中选择距离该坐标点最近点灰度值作为目标图中的灰度值。

步骤3：对预处理后的图像进行插值以便打散原轨迹。

图像插值一般都是利用已知邻近像素点的灰度值来产生未知像素点的灰度值，以便由原始图像再生出具有更高分辨率的图像。在此处使用插值是直接将空白的像素点插入图像，以便获得被打散的轨迹。

步骤4：将步骤3中被打散成单个点的轨迹类比成目标跟踪中的目标位置，并引入联合概率数据关联来识别出不同雷达信号的轨迹。

为表示有效回波和各目标跟踪门的复杂关系,引入确认矩阵的概念。确认矩阵定义为：

Ω＝(ω_jt),j＝1,2,...,m(k),t＝0,1,...,n (9)

其中ω_jt用以表示量测j是否落入目标t的确认门之内。t＝0表示没有目标，此时对应的Ω的列的元素都是1，这是因为任意测量都可能源于杂波或虚警。量测落入跟踪门相交区域的情形，对应某些量测可能源于多个目标，联合概率数据关联的目的就是计算每个量测与其可能的各种源目标相关联的概率。

假设在杂波环境中已有T个目标，则他们的状态方程和测量方程分别表示为：

JPDA的状态更新方程为：

其中，表示k时刻目标t的状态矢量，表示目标t的状态预测矢量，V′(k)为联合新息：

其中，β_jt表示目标t与观测j的关联概率，V_j′(k)表示滤波新息。

在JPDAF中，由于它考虑了所有目标和航迹的可行联合事件，关联概率可以表示为：

其中，表示在可行事件θ^j中，观测j是否源于目标t。如果源于目标则等于1，否则为0。另外观测还必须服从两个假设：第一，每一个观测有唯一源；第二，对于一个给定的目标，最多只有一个观测以其为源。在k时刻联合事件θ的条件概率可表示为

其中，σ_t为目标检测指示器，如果目标在θ中与观测关联，则σ_t的值等于1，反之为零。Φ是虚假观测事件的个数，表示目标t的检测概率，V表示航迹有效门体积。

步骤5：运用中心矩提取时频图像的形状特征。

设大小为M*N的二值图像可以表示为：f(x,y)∈{1,0}，则图像的(p+q)阶原点矩可以定义为：

图像的(p+q)阶中心矩可以定义为：

其中表示水平方向上的质心。表示垂直方向的质心。

图像的不同阶数中心矩表征不同的物理意义，其中μ₀₂表示图像在垂直方向上的伸展度；μ₂₀表示图像在水平方向上的伸展度；μ₁₁表示图像的倾斜度；μ₀₃表示图像在垂直方向上的重心偏移度；μ₁₂表示图像垂直伸展的均衡程度。由于时频图的差异性主要体现在垂直方向的频域上，故舍弃μ₀₂，μ₀₂这些描述水平方向时域的特征值，对此可以选用μ₀₂，μ₀₂，μ₀₂，μ₀₂，μ₀₂作为时频图像特征。而构造出特征向量。

步骤6：选取支持向量机分类器的核函数，并根据训练样本的形态特征确定核函数的参数，以训练支持向量机分类器。

其中，支持向量机分类器的核函数为高斯核函数，高斯核函数为：

k(x₁,x₂)＝exp(-γ*|x₁-x₂|²)，其中，采用交叉验证的方法确定所述高斯核函数的参数。

对于n类雷达辐射源信号脉内调制方式识别问题，选择任意两类设计1个二分类SVM，这样，总共可设计n(n-1)/2个二分类SVM；然后将这n(n-1)/2个二分类SVM构成一个有向无循环图，如图2所示，该图有n(n-1)/2个非叶子节点和n个叶子节点，每个非叶子节点表示一个二分类SVM分类器，并与下一层的两个节点相连，每个叶子节点表示一个输出。

具体而言，选择上述的高斯核函数作为核函数，核函数参数的选择采用交叉验证的方法。其中，C(惩罚因子)从2^-5～2⁵(成倍增加)，g(核函数参数γ)从2^-5～2⁵(成倍增加)，对于某个(C,g)，选择5折交叉验证，即将训练样本分成5份，挑选1份作为测试样本，剩下4份作为训练样本，重复5次，计算平均正确率，以此衡量当前(C,g)的性能，选择性能相对最好的(C,g)作为最终的参数。

最后，将步骤5得到的特征向量输入训练好的支持向量机，得到同一时刻多信号频率交叠的雷达信号脉内调制方式。

Claims (1)

1.一种基于时频图像特征的脉内调制识别方法，其特征在于，包括如下步骤：

(1)对雷达信号进行Wigner-Vill分布时频变换，得到反映信号能量随时间和频率的时频分布图；

对相关函数作滑窗处理，得到时变的局部相关函数公式，

$R(t,\mathrm{\τ})={\∫}_{\mathrm{\∞}}^{\mathrm{\∞}}s(u+\frac{\mathrm{\τ}}{2})s^{*}(u-\frac{\mathrm{\τ}}{2})\mathrm{\φ}(u-t,\mathrm{\τ})du$

当窗函数取时间冲击函数，不加限制，而在时域取瞬时值，

$R(t,\mathrm{\τ})={\∫}_{\mathrm{\∞}}^{\mathrm{\∞}}s(u+\frac{\mathrm{\τ}}{2})s^{*}(u-\frac{\mathrm{\τ}}{2})\mathrm{\σ}(u-t,\mathrm{\τ})du=s(t+\frac{\mathrm{\τ}}{2})s^{*}(t-\frac{\mathrm{\τ}}{2})$

对时变局部相关函数作Fourier变换，即可得到WVD，

$MVD(t,f)=\∫s(t+\mathrm{\τ}/2)s^{*}(u-\mathrm{\τ})du=s(t+\frac{\mathrm{\τ}}{2})s^{*}(t-\frac{\mathrm{\τ}}{2})$

设有n个分量信号，得到多分量信号的WVD，

${\mathrm{MVD}}_x(t,f)=\underset{k}{\Σ}{\mathrm{MVD}}_{x_k}(t,f)+\underset{k}{\Σ}\underset{l\≠k}{\Σ}2\mathrm{Re}\[{\mathrm{MVD}}_{x_k{x}_l}(t,f)\]$

雷达信号经过Wigner-Vill分布时频变换后得到时频图像；

(2)对步骤(1)得到的时频图像进行图像预处理；

对时频图像灰度值进行归一化；设和σ²分别是所有像素点灰度值的均值和方差；

$\stackrel{\‾}{x}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{x}_i$

${\mathrm{\σ}}^2=\frac{1}{N-1}\underset{i-1}{\overset{N}{\Σ}}{(x_i-\stackrel{\‾}{x})}^2$

归一化后像素点的灰度值为：

采用自适应维纳滤波器初步去除时频图像的噪声点，对图像进行增强；采用形态学方法对灰度图像进行处理；形态学处理是应用具有一定形态的结构元素对灰度图像进行腐蚀和膨胀的操作，其中使用结构元素b对f的灰度膨胀记为

选择半径为3的菱型结构元素对时频图像进行闭运算；

通过检测信号时频分布的起止频率，将没有信号分布的图像区域剪切掉；

采用最近邻域插值法归一化时频图像的纵横比；

将原图定义为一个二维函数F(x,y)，(x,y)点的函数值对应像素点的灰度值，f(x,y)为归一化纵横比后的目标图；

Fw、Fh设为原图的宽度和高度；fw、fh设为目标图的宽度和高度；原图坐标(x,y)和目标图坐标(x′,y′)之间公式为：

$\begin{array}{c}x=x^{\′}(Fw/fw)\\ y=y^{\′}\left(Fh/fh\right)\end{array};$

(3)对预处理后的图像进行插值打散原轨迹；

(4)将步骤(3)中被打散成单个点的轨迹类比成目标跟踪中的目标位置，并引入联合概率数据关联来识别出不同雷达信号的轨迹；

引入确认矩阵；确认矩阵定义为：

Ω＝(ω_jt),j＝1,2,...,m(k),t＝0,1,...,n

其中ω_jt表示量测j是否落入目标t的确认门之内；t＝0表示没有目标，此时对应的Ω的列的元素都是1；

设在杂波环境中已有T个目标，则他们的状态方程和测量方程分别表示为：

$\left\{\begin{array}{c}{X}^{\′}(k+1)=F^{\′}\left(k\right)X^{\′}\left(k\right)+W^{\′}\left(k\right),k=0,1,2,...;t=1,2,...,T\\ Z\left(k\right)=H\left(k\right)X\left(k\right)+V\left(k\right),k=0,1,2,...;\end{array}\right.$

JPDA的状态更新方程为：

其中，表示k时刻目标t的状态矢量，表示目标t的状态预测矢量，V′(k)为联合新息：

$V^{\′}\left(k\right)=\underset{j}{\Σ}{\mathrm{\β}}_{jt}{V}_j^{\′}\left(k\right)$

其中，β_jt表示目标t与观测j的关联概率，V′_j(k)表示滤波新息；

关联概率表示为：

其中，表示在可行事件θ^j中，观测j是否源于目标t；如果源于目标则等于1，否则为0；在k时刻联合事件θ的条件概率表示为

$P\{\mathrm{\θ}/Z^k\}=\frac{1}{c}\frac{\mathrm{\Φ}!}{v^{\mathrm{\Φ}}}\underset{j=1}{\overset{m_k}{\Π}}{\[N_{tj}\left(Z_j\right(k\left)\right)\]}^{{\mathrm{\τ}}_j}\underset{t=1}{\overset{T}{\Π}}{\left(P_D^t\right)}^{\mathrm{\σ}t}{(1-P_D^t)}^{1-\mathrm{\σ}t}$

其中，σ_t为目标检测指示器，如果目标在θ中与观测关联，则σ_t的值等于1，反之为零；Φ是虚假观测事件的个数，表示目标t的检测概率，V表示航迹有效门体积；

(5)运用中心矩提取时频图像的形状特征；

设大小为M*N的二值图像表示为：f(x,y)∈{1,0}，则图像的(p+q)阶原点矩定义为：

$m_{pq}=\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}{x}^p{y}^{q}f(x,y)p,q\∈N$

图像的(p+q)阶中心矩定义为：

${\mathrm{\μ}}_{pq}=\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}{(x-\stackrel{\‾}{x})}^p{(y-\stackrel{\‾}{y})}^{q}f(x,y)$

其中表示水平方向上的质心；表示垂直方向的质心；

其中μ₀₂表示图像在垂直方向上的伸展度；μ₂₀表示图像在水平方向上的伸展度；μ₁₁表示图像的倾斜度；μ₀₃表示图像在垂直方向上的重心偏移度；μ₁₂表示图像垂直伸展的均衡程度；由于时频图的差异性主要体现在垂直方向的频域上，故舍弃μ₀₂，而构造出特征向量；

(6)选取支持向量机分类器的核函数，并根据训练样本的形态特征确定核函数的参数，以训练支持向量机分类器；

其中，支持向量机分类器的核函数为高斯核函数，高斯核函数为：

k(x₁,x₂)＝exp(-γ*|x₁-x₂|²)，

选择任意两类设计1个二分类SVM，总共设计n(n-1)/2个二分类SVM；然后将这个二分类SVM构成一个有向无循环图，该图有n(n-1)/2个非叶子节点和n个叶子节点，每个非叶子节点表示一个二分类SVM分类器，并与下一层的两个节点相连，每个叶子节点表示一个输出；

具体而言，选择上述的高斯核函数作为核函数，核函数参数的选择采用交叉验证的方法；其中，惩罚因子C从2^-5成倍增加到2⁵，核函数参数g从2^-5成倍增加到2⁵，对于(C,g)，选择5折交叉验证，即将训练样本分成5份，挑选1份作为测试样本，剩下4份作为训练样本，重复5次，计算平均正确率，以此衡量当前(C,g)的性能，选择性能相对最好的(C,g)作为最终的参数；最后，将步骤(5)得到的特征向量输入训练好的支持向量机。