一种基于相位噪声无意调制特征的雷达辐射源识别方法.pdf

一种基于相位噪声无意调制特征的雷达辐射源识别方法，涉及一种雷达辐射源的识别方法。本发明为了解决现有的基于相位噪声的辐射源识别方法的识别率不高的问题。本发明分析雷达发射机系统中锁相式频率合成器的结构，建立锁相式频率合成器产生的相位噪声模型，计算双谱对角切片特征和双谱反对角切片特征，然后将双谱对角切片特征矩阵A1和双谱反对角切片特征矩阵B1组成特征矩阵Y，通过PCA降维后，建立已知类型发射机向量机模型，然后利用已经建立的向量机模型完成对未知类型发射机的发射信号的识别；从而实现雷达辐射源的识别。本发明适用于雷达辐射源的识别。

1.  一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：包括以下步骤：
步骤1：根据雷达发射机系统中锁相式频率合成器的结构，建立锁相式频率合成器产生的相位噪声模型如下：
s(t)=sin(2πfct)+Σm′Δφm′2sin[2π(fc+fm′)t]-Σm′Δφm′2sin[2π(fc-fm′)t]---(7)]]>
其中，s(t)表示加入相位噪声的信号，f_c表示原始信号的频率，m′表示调制信号的序号，Δφ_m′示第m′个调制信号的调相系数，f_m′表示第m′个调制信号频率，t表示脉冲宽度；
步骤2：多次计算双谱对角切片特征，得到双谱对角切片样本组成矩阵A₁；多次计算双谱反对角切片特征，得到双谱反对角切片样本组成矩阵B₁；将双谱对角切片特征矩阵A₁和双谱反对角切片特征矩阵B₁组成特征矩阵Y=A1B1;]]>
步骤3、用PCA方法对特征矩阵Y进行降维，得到降低维数后的矩阵Z；
步骤4、选择矩阵Z的前50％列的样本作为训练样本，用支持向量机方法建立已知发射机类型的向量机模型，后50％列的样本作为测试样本，根据向量机模型得到能将其正确分类的概率，从而验证建立的向量机模型的分类识别性能；
计算得到某一未知类型发射机的发射信号的矩阵Z，并将其输入已知发射机类型的向量机模型；计算该未知类型发射信号与向量机模型对应的已知发射机类型的相似概率，并根据相似概率判断该发射信号对应的发射机类型。

2.  根据权利要求1所述的一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：
步骤1所述的建立锁相式频率合成器产生的相位噪声模型公式(7)具体步骤如下：
设原始信号的频率为f_c，加入单个相位噪声后信号表示为
s(t)＝sin(2πf_ct+Δφ(t))＝sin(2πf_ct+Δφ_msin(2πf_mt))   (1)
其中，Δφ(t)表示调制信号，Δφ_m为调相系数，t为脉冲宽度，f_m为调制信号频率；
将上式展开为
s(t)＝sin(2πf_ct)cos(Δφ_msin(2πf_mt))+cos(2πf_ct)sin(Δφ_msin(2πf_mt))   (2)
贝塞尔函数中(3)和(4)的表达式如下：
cos(Δφ_msin(2πf_mt))＝J₀(Δφ_m)+2[J₂(Δφ_m)cos(2πf_mt)+J₄(Δφ_m)cos(8πf_mt)+…]   (3)
sin(Δφ_msin(2πf_mt))＝2[J₁(Δφ_m)sin(2πf_mt)+J₃(Δφ_m)sin(6πf_mt)+…]   (4)
其中，上式中的J₀(Δφ_m)、J₁(Δφ_m)、J₂(Δφ_m)、…、J_n(Δφ_m)分别表示调相系数为Δφ_m的0、1、2、…、n阶贝塞尔函数；
贝塞尔函数的近似值为
Δφm<<1,J0(Δφm)&ap;1,J1(Δφm)&ap;Δφm2,Jn(Δφm)&ap;0,n&GreaterEqual;2---(5)]]>
将(3)-(5)带入(2)，得到
s(t)=sin(2πfct)+Δφm2sin[2π(fc+fm)t]-Δφm2sin[2π(fc-fm)t]---(6)]]>
加入多个相位噪声的信号可表示为
s(t)=sin(2πfct)+Σm′Δφm′2sin[2π(fc+fm′)t]-Σm′Δφm′2sin[2π(fc-fm′)t]---(7)]]>
其中，s(t)表示加入相位噪声的信号，f_c表示原始信号的频率，m′表示调制信号的序号，Δφ_m′示第m′个调制信号的调相系数，f_m′表示第m′个调制信号频率，t表示脉冲宽度。

3.  根据权利要求1或2所述的一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：
步骤2所述的多次计算双谱对角切片特征得到双谱对角切片样本组成矩阵A₁的具体步骤如下：
步骤2.1.1：将接收到的信号s(t)简记为s，对其进行分段处理，每段取M个数据，共分为P段，对每一段中的数据进行去均值处理：
xi(n′)=si(n′)-1MΣn′=1Msi(n′)---(8)]]>
其中，信号s分段后第i段中的数据表示为sⁱ(n)，n′＝1，2，…M，i＝1，2，…P；xⁱ(n′)是第i段去均值后的数据；
步骤2.1.2：计算P段中各段中的三阶累积量函数τ₁,τ₂为函数的变量：
令τ₁＝τ₂＝τ时，三阶累积量为
cxi(τ,τ)=1MΣn′=p1p2xi(n′)xi(n′+τ)xi(n′+τ)---(9)]]>
其中，p₁＝max(1,1-τ)，p₂＝min(M,M-τ)，表示τ₁＝τ₂＝τ时各段的三阶累积量；
步骤2.1.3：.对各段的三阶累积量取均值，
cx(τ,τ)=1PΣi=1Pcxi(τ,τ)---(10)]]>
c_x(τ,τ)表示τ₁＝τ₂＝τ时各段的三阶累积量的均值；
步骤2.1.4：对c_x(τ,τ)进行傅里叶变换，得到双谱对角切片；
步骤2.1.5：有噪声时，仿真一次得到一个双谱对角切片样本a_j，仿真N次得到N个双谱对角切片，将N个双谱对角切片样本组成矩阵A₁＝[a₁^Ta₂^T……a_N^T]，将A₁记为双谱对角切片特征矩阵。

4.  根据权利要求3所述的一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：
步骤2所述的多次计算双谱反对角切片特征得到双谱反对角切片样本组成矩阵B₁的具体步骤如下：
步骤2.2.1：对s进行分段处理，每段取M个数据，共分为P段，对每一段中的数据进行去均值处理如下；
xi(n′)=si(n′)-1MΣn′=1Msi(n′)---(11)]]>
步骤2.2.2：计算三阶累积量令τ₁＝-τ₂＝τ，三阶累积量为
cxi(τ,-τ)=1MΣn′=p3p4xi(n′)xi(n′+τ)xi(n′-τ)---(12)]]>
其中，p₃＝max(1,1-τ,1+τ)，p₄＝min(M,M-τ,M+τ)，表示τ₁＝-τ₂＝τ时各段的三阶累积量；
步骤2.2.3：各段的三阶累积量取均值
cx(τ,-τ)=1PΣi=1Pcxi(τ,-τ)---(13)]]>
c_x(τ,-τ)表示τ₁＝-τ₂＝τ时各段的三阶累积量的均值；
步骤2.2.4：对c_x(τ,-τ)进行傅里叶变换，得到双谱反对角切片；
步骤2.2.5：有噪声时，仿真一次得到一个双谱反对角切片样本b_j，仿真N次得到N个双谱反对角切片样本，将N个双谱反对角切片样本组成矩阵B₁＝[b₁^Tb₂^T……b_N^T]，将B₁记为双谱反对角切片特征矩阵。

5.  根据权利要求4所述的一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：
步骤3的具体步骤如下：
Y=A1B1=y11,y12...y1Ny21,y22...y2N...yr1,yr2...yrN]]>
Y的每一列为一个观察样本y＝(y_1k,y_2k,…,y_rk)^T为r维随机矢量，k∈[1,N]，每一行代表一维数据；
步骤3.1：计算样本矩阵Y每一行的均值，分别记为E[y₁]、…、E[y_r]，并对矩阵中每个观察样本y做如下处理：
y＝y-E[y]＝(y_1k-E[y₁],…,y_rk-E[y_r])^T   (14)
去均值后的样本矩阵记为Y′；
步骤3.2：计算去均值后的样本矩阵Y′的协方差矩阵：

其中，cov(·)表示求协方差；
步骤3.3：计算协方差矩阵C_y的特征值λ_u及相应特征向量q_u，其中u＝1,2,…,r；
步骤3.4：将特征值λ_u按降序排列，计算前g个较大特征值对整体信息保存的完整程度，即累积贡献率η(g)；当累积贡献率η(g)大于预设率值时，选择此时对应的g，执行步骤3.5；
η(g)=Σk=1gλk/Σk=1rλk---(16)]]>
步骤3.5：利用取得的前g个较大特征值所对应的特征向量构成变换矩阵Q′：
Q′＝[q₁,q₂,…,q_g]，g<r   (17)
步骤3.6：通过(18)计算前g个主成分，得到降低维数的矩阵Z；
Z＝(Q′)^TY   (18)。

6.  根据权利要求5所述的一种基于相位噪声无意调制特征的雷达辐射源识别方法，其特征在于：
步骤3.4所述的预设率值为85％，即累积贡献率η大于85％。

一种基于相位噪声无意调制特征的雷达辐射源识别方法
技术领域
本发明涉及一种雷达辐射源的识别方法。
背景技术
在早期的辐射源识别中，主要是分析其常规特征参数和辐射源信号的脉内调制特征，脉内调制特征分析还主要集中在信号脉内的有意调制特征上。后来，人们才意识到基于脉内有意调制特征提取也不能适应当今愈加复杂的战场环境和配置愈加先进的各类新式雷达，因为新式雷达大多数具有复杂的多种工作方式和调制方式，其脉内有意调制可以任意改变，非合作方很难检测到其雷达工作体制。在这种情况下，无意调制特征的研究渐渐进入了研究人员的视线。
当发射机激活时，在经过一段时间后，从本机显示发出的信号，这个信号具有独特的特点，可以用于明确地标识单个发射机，这就是信号的“指纹”，也称为无意调制特征。无意调制特征包括幅度变化、频率漂移和相位噪声，基于无意调制特征的辐射源识别方法的研究较少，相位噪声对信号的影响较大且特征较稳定，但目前的研究中，基于相位噪声的辐射源识别方法的识别率不高。
发明内容
本发明为了解决现有的基于相位噪声的辐射源识别方法的识别率不高的问题。
一种基于相位噪声无意调制特征的雷达辐射源识别方法，包括以下步骤：
步骤1：根据雷达发射机系统中锁相式频率合成器的结构，建立锁相式频率合成器产生的相位噪声模型如下：
      s(t)=sin(2πfct)+Σm′Δφm′2sin[2π(fc+fm′)t]-Σm′Δφm′2sin[2π(fc-fm′)t]---(7)]]>      
其中，s(t)表示加入相位噪声的信号，f_c表示原始信号的频率，m′表示调制信号的序号，Δφ_m′示第m′个调制信号的调相系数，f_m′表示第m′个调制信号频率，t表示脉冲宽度；
步骤2：多次计算双谱对角切片特征，得到双谱对角切片样本组成矩阵A₁；多次计算双谱反对角切片特征，得到双谱反对角切片样本组成矩阵B₁；将双谱对角切片特征矩阵A₁和双谱反对角切片特征矩阵B₁组成特征矩阵Y=A1B1;]]>
步骤3、用PCA方法对特征矩阵Y进行降维，得到降低维数后的矩阵Z；
步骤4、选择矩阵Z的前50％列的样本作为训练样本，用支持向量机方法建立已知发射机类型的向量机模型，后50％列的样本作为测试样本，根据向量机模型得到能将其正确分类的概率，从而验证建立的向量机模型的分类识别性能；
计算得到某一未知类型发射机的发射信号的矩阵Z，并将其输入已知发射机类型的向量机模型；计算该未知类型发射信号与向量机模型对应的已知发射机类型的相似概率，并根据相似概率判断该发射信号对应的发射机类型。
本发明具有以下优点：
与仅仅利用双谱对角切片或者反对角切片相比，本发明计算双谱对角切片和双谱反对角切片联合特征的计算量并没有增加很多，但却更能完整地反应信号的特征，从而建立的识别模型更加准确，在信噪比为-8dB到15dB时，本发明通过将正反对角切片特征联合起来，再用PCA降维处理，建立模型的准确识别率达到91％以上；由于建立的支持向量机的识别模型更加准确，从而对雷达辐射源识别判断上更加准确，取得了更好的雷达辐射源识别分类效果。
附图说明
图1为双谱对角切片反对角切片计算流程图；
图2为PCA计算流程图；
图3(a)为加入相位噪声的信号源1的仿真图；
图3(b)为信号源1的双谱图；
图3(c)为信号源1的双谱对角切片图；
图3(d)为信号源1的双谱反对角切片图；
图4为不同信噪比下的不同方法建立支持向量机模型的识别率效果图。
具体实施方式
具体实施方式一：一种基于相位噪声无意调制特征的雷达辐射源识别方法，包括以下步骤：
步骤1：根据雷达发射机系统中锁相式频率合成器的结构，建立锁相式频率合成器产生的相位噪声模型如下：
      s(t)=sin(2πfct)+Σm′Δφm′2sin[2π(fc+fm′)t]-Σm′Δφm′2sin[2π(fc-fm′)t]---(7)]]>      
其中，s(t)表示加入相位噪声的信号，f_c表示原始信号的频率，m′表示调制信号的序 号，Δφ_m′示第m′个调制信号的调相系数，f_m′表示第m′个调制信号频率，t表示脉冲宽度；
步骤2：如图1，多次计算双谱对角切片特征，得到双谱对角切片样本组成矩阵A₁；多次计算双谱反对角切片特征，得到双谱反对角切片样本组成矩阵B₁；将双谱对角切片特征矩阵A₁和双谱反对角切片特征矩阵B₁组成特征矩阵Y=A1B1;]]>
步骤3、用PCA方法对特征矩阵Y进行降维，得到降低维数后的矩阵Z；
步骤4、选择矩阵Z的前50％列的样本作为训练样本，用支持向量机方法建立已知发射机类型的向量机模型，后50％列的样本作为测试样本，根据向量机模型得到能将其正确分类的概率，从而验证建立的向量机模型的分类识别性能；
计算得到某一未知类型发射机的发射信号的矩阵Z，并将其输入已知发射机类型的向量机模型；计算该未知类型发射信号与向量机模型对应的已知发射机类型的相似概率，并根据相似概率判断该发射信号对应的发射机类型。
具体实施方式二：本实施方式所述的
步骤1所述的建立锁相式频率合成器产生的相位噪声模型公式(7)具体步骤如下：
相位噪声可看作是对原始信号的相位进行频率为f_m的正弦波调制，以单频信号为例，设原始信号的频率为f_c，加入单个相位噪声后信号表示为
s(t)＝sin(2πf_ct+Δφ(t))＝sin(2πf_ct+Δφ_msin(2πf_mt))  (1)
其中，Δφ(t)表示调制信号，Δφ_m为调相系数，t为脉冲宽度，f_m为调制信号频率；
将上式展开为
s(t)＝sin(2πf_ct)cos(Δφ_msin(2πf_mt))+cos(2πf_ct)sin(Δφ_msin(2πf_mt))  (2)
贝塞尔函数中(3)和(4)的表达式如下：
cos(Δφ_msin(2πf_mt))＝J₀(Δφ_m)+2[J₂(Δφ_m)cos(2πf_mt)+J₄(Δφ_m)cos(8πf_mt)+…]  (3)
sin(Δφ_msin(2πf_mt))＝2[J₁(Δφ_m)sin(2πf_mt)+J₃(Δφ_m)sin(6πf_mt)+…]  (4)
其中，上式中的J₀(Δφ_m)、J₁(Δφ_m)、J₂(Δφ_m)、…、J_n(Δφ_m)分别表示调相系数为Δφ_m的0、1、2、…、n阶贝塞尔函数；
贝塞尔函数的近似值为
      Δφm<<1,J0(Δφm)&ap;1,J1(Δφm)&ap;Δφm2,Jn(Δφm)&ap;0,n&GreaterEqual;2---(5)]]>            
将(3)-(5)带入(2)，得到
      s(t)=sin(2πfct)+Δφm2sin[2π(fc+fm)t]-Δφm2sin[2π(fc-fm)t]---(6)]]>      
相位噪声可看作是受调制的单频信号组合而成，所以，考虑到多个不同频率分量的调制信号对雷达信号相位的影响，加入多个相位噪声的信号可表示为
      s(t)=sin(2πfct)+Σm′Δφm′2sin[2π(fc+fm′)t]-Σm′Δφm′2sin[2π(fc-fm′)t]---(7)]]>      
其中，s(t)表示加入相位噪声的信号，f_c表示原始信号的频率，m′表示调制信号的序号，Δφ_m′示第m′个调制信号的调相系数，f_m′表示第m′个调制信号频率，t表示脉冲宽度。
其它步骤与具体实施方式一相同。
具体实施方式三：本实施方式所述的
步骤2所述的多次计算双谱对角切片特征得到双谱对角切片样本组成矩阵A₁的具体步骤如下：
步骤2.1.1：将接收到的信号s(t)简记为s，对其进行分段处理，每段取M个数据，共分为P段，对每一段中的数据进行去均值处理：
      xi(n′)=si(n′)-1MΣn′=1Msi(n′)---(8)]]>      
其中，信号s分段后第i段中的数据表示为sⁱ(n)，n′＝1，2，…M，i＝1，2，…P；xⁱ(n′)是第i段去均值后的数据；
步骤2.1.2：计算P段中各段中的三阶累积量函数τ₁,τ₂为函数的变量：
令τ₁＝τ₂＝τ时，三阶累积量为
      cxi(τ,τ)=1MΣn′=p1p2xi(n′)xi(n′+τ)xi(n′+τ)---(9)]]>      
其中，p₁＝max(1,1-τ)，p₂＝min(M,M-τ)，表示τ₁＝τ₂＝τ时各段的三阶累积量；
步骤2.1.3：.对各段的三阶累积量取均值，
      cx(τ,τ)=1PΣi=1Pcxi(τ,τ)---(10)]]>            
c_x(τ,τ)表示τ₁＝τ₂＝τ时各段的三阶累积量的均值；
步骤2.1.4：对c_x(τ,τ)进行傅里叶变换，得到双谱对角切片；
步骤2.1.5：有噪声时，仿真一次得到一个双谱对角切片样本a_j，仿真N次得到N个双谱对角切片，将N个双谱对角切片样本组成矩阵A₁＝[a₁^Ta₂^T……a_N^T]，将A₁记为双谱对角切片特征矩阵。
其它步骤与具体实施方式二相同。
具体实施方式四：本实施方式所述的
步骤2所述的多次计算双谱反对角切片特征得到双谱反对角切片样本组成矩阵B₁的具体步骤如下：
步骤2.2.1：对s进行分段处理，每段取M个数据，共分为P段，对每一段中的数据进行去均值处理如下；
      xi(n′)=si(n′)-1MΣn′=1Msi(n′)---(11)]]>      
步骤2.2.2：计算三阶累积量令τ₁＝-τ₂＝τ，三阶累积量为
      cxi(τ,-τ)=1MΣn′=p3p4xi(n′)xi(n′+τ)xi(n′-τ)---(12)]]>      
其中，p₃＝max(1,1-τ,1+τ)，p₄＝min(M,M-τ,M+τ)，表示τ₁＝-τ₂＝τ时各段的三阶累积量；
步骤2.2.3：各段的三阶累积量取均值
      cx(τ,-τ)=1PΣi=1Pcxi(τ,-τ)---(13)]]>      
c_x(τ,-τ)表示τ₁＝-τ₂＝τ时各段的三阶累积量的均值；
步骤2.2.4：对c_x(τ,-τ)进行傅里叶变换，得到双谱反对角切片；
步骤2.2.5：有噪声时，仿真一次得到一个双谱反对角切片样本b_j，仿真N次得到N个双谱反对角切片样本，将N个双谱反对角切片样本组成矩阵B₁＝[b₁^Tb₂^T……b_N^T]，将B₁记为双谱反对角切片特征矩阵。
其它步骤与具体实施方式三相同。
具体实施方式五：结合图2说明本实施方式，本实施方式所述的步骤3的具体步骤如下：
      Y=A1B1=y11,y12...y1Ny21,y22...y2N............yr1,yr2...yrN]]>      
Y的每一列为一个观察样本y＝(y_1k,y_2k,…,y_rk)^T为r维随机矢量，k∈[1,N]，每一行代表一维数据；
步骤3.1：计算样本矩阵Y每一行的均值，分别记为E[y₁]、…、E[y_r]，并对矩阵中每个观察样本y做如下处理：
y＝y-E[y]＝(y_1k-E[y₁],…,y_rk-E[y_r])^T  (14)
去均值后的样本矩阵记为Y′；
步骤3.2：计算去均值后的样本矩阵Y′的协方差矩阵：
            
其中，cov(·)表示求协方差；
步骤3.3：计算协方差矩阵C_y的特征值λ_u及相应特征向量q_u，其中u＝1,2,…,r；
步骤3.4：将特征值λ_u按降序排列，计算前g个较大特征值对整体信息保存的完整程度，即累积贡献率η(g)；当累积贡献率η(g)大于预设率值时，选择此时对应的g，执行步骤3.5；
      η(g)=Σk=1gλk/Σk=1rλk---(16)]]>      
步骤3.5：利用取得的前g个较大特征值所对应的特征向量构成变换矩阵Q′：
Q′＝[q₁,q₂,…,q_g]，g<r  (17)
步骤3.6：通过(18)计算前g个主成分，得到降低维数的矩阵Z；
Z＝(Q′)^TY  (18)。
其它步骤与具体实施方式四相同。
具体实施方式六：本实施方式所述的
步骤3.4所述的预设率值为85％，即累积贡献率η大于85％。
其它步骤与具体实施方式五相同。
具体实施例
利用公式(7)建立加入相位噪声的信号源模型，定义原始信号f_c＝4MHz，t的取值范围为0～0.01ms，采样频率为40MHz；以发射机个体1的信号源为例，加入相位噪声后信号的双谱及双谱切片如图3(a)—3(d)所示；
假设五个雷达发射机个体发射的理想单频信号是相同的，但由于各个发射机不理想，混入了不同的相位噪声，调相系数和调制信号频率分别定义如表1所示；
表1相位噪声类型
            
分别对原始信号加入五种不同的相位噪声，完成信号源建模，然后分别对五种信号源通过步骤2、3计算其特征，最后根据步骤4进行分类识别；对比只用双谱对角切片特征和双谱反对角切片特征的分类效果，三种方法的分类识别率效果对比如图4所示。