说明书一种基于稀疏频率的雷达通信波形设计方法
技术领域
本发明属于信号处理领域，特别涉及一种雷达嵌入通信的稀疏频率波形设计方法。
背景技术
雷达的主要功能是探测、目标定位和引导。通信系统和雷达的构成原理类似，均包含电磁波的发射和接收过程。在实现频段、天线等资源共享的前提下，同一个系统可以实现雷达探测和通信功能，这就是雷达通信综合化的概念。雷达通信综合化波形设计的三种途径：一、在雷达波形上调制通信数据；二、通信信号改造为雷达探测波形；三、雷达波形和通信信号独立产生，然后叠加合成。
一般雷达通信综合化研究中的通信系统，主要为两个或多个信号端(发射机/接收机)之间进行大带宽的数据通信。雷达通信综合化的另一种模式为脉内雷达嵌入通信。与一般的雷达通信综合化主要研究雷达波形与通信信号的共享能力相比，脉内雷达嵌入通信还重点考虑了通信信号的隐蔽性。这种通信的隐蔽性，对雷达识别和定位目标起到重要的作用。
雷达嵌入通信的基本原理为：在脉冲雷达探照范围内设置一个雷达响应标识(Tag)系统，接收雷达波形并进行再调制；再调制的结果是在回波中嵌入通信信号；Tag系统最后将混合波形发送至雷达接收机达到通信的目的。由于Tag系统中再调制过程的特殊性，通信过程可以达到隐蔽通信信息的目的。
然而，用线性调频(LFM)脉冲实现的雷达嵌入通信，所嵌入的通信信号只能利用较少的频带资源，从而导致通信样本数较少、通信效率低。
稀疏频率波形设计是雷达波形设计的一个重要研究方向。对于雷达系统，稀疏频率波形具有抑制干扰、改善检测性能等优点。对于设计一个特定的稀疏频率波形，可以用功率谱密度匹配的方法得到目标函数，用拟牛顿法可以对目标函数进行求解。
发明内容
有鉴于此，本发明针对用线性调频(LFM)脉冲实现的雷达嵌入通信，所嵌入的通信信号只能利用较少的频带资源，从而导致通信样本数较少、通信效率低的问题，提出了一种雷达嵌入通信的稀疏频率波形设计方法。该方法可以提高通信信号的频带利用率和通信样本数，同时保证低误码率、低拦截率和雷达的探测特性。
本发明解决其技术问题所采用的技术方案是：
一种基于稀疏频率的雷达通信波形设计方法，其实现步骤如下：(1)雷达标识(Tag)系统接收雷达散射回波，对雷达散射回波采用功率谱密度匹配方法设计稀疏频率雷达波形；(2)然后对该稀疏频率波形的序列进行特征值分解，得到特征向量组；用其中的一部分特征向量设计得出通信波形信号。(3)然后Tag系统在回波中嵌入该通信信号，然后将混合波形发送至雷达接收机。那么雷达接收机接收雷达散射回波信号与雷达嵌入通信波形，实现雷达探测与通信两种功能，其中通信过程可以达到隐蔽的目的。
其中，步骤(1)中采用功率谱密度匹配方法时，得到的目标函数为四次方无约束非凸优化问题，可以应用拟牛顿法求解，从而获得相应的稀疏频率雷达波形，实际上还可以采用交替投影法或遗传算法等经典优化算法求解；
其中，步骤(2)中通信波形的获得途径也可以采用如下步骤S21和S22实现，即S21：对稀疏频率雷达波形的序列进行特征分解，得到波形序列的主空间；S22：由正交伪随机向量在所述波形序列的主空间的垂直子空间进行投影，得到多个该子空间的基向量，这些基向量组成通信波形。在实际中，也可以通过对波形序列的特征值的其它变换获得，不同的变换方式会影响通信波形之间的正交性、通信波形与雷达波形的相关性。
与一般的雷达通信综合化主要研究雷达波形与通信信号的共享能力相比，这种混合信号重点考虑通信样本的隐蔽性，将Tag通信信号嵌入在雷达回波中，避免了被Tag所在目标区的敌方侦测设备发现。
所述步骤(1)中的基于雷达标识(Tag)系统应用，建立雷达嵌入通信模型。具体为：
设Tag系统输出信号yr(t)为：
yr(t)＝αkck(t)+ys(t)+n(t) (1)
其中，n(t)为系统噪声；ys(t)表示Tag系统接收到的雷达散射回波；αk表示路径损失等对通信信号的混合影响；ck(t)为所设计的通信信号；k＝1,2,...,K表示第k个通信样本，K表示通信样本总数。
对式(1)用向量运算表示为如下所示：
yr＝αkck+Sx+n (2)
其中，ck表示通信信号；n表示噪声向量；x为雷达散射的距离样本；Sx是雷达波形与散射响应卷积过程的离散表示。
所述步骤(1)通过功率谱密度匹配及拟牛顿法设计稀疏频率雷达波形的实施步骤如下：
步骤a.给定了功率谱密度u的分布后，稀疏频率雷达波形设计的目标函数用表示，其中x为所设计波形，其相位矢量A为离散傅里叶变换矩阵，A中的元素Amn＝exp(-i×2π×n×m/Nf)，Nf与雷达波形的采样数相等；表示理想稀疏频率波形的功率谱密度。在每一次迭代中计算目标函数的梯度，关于Θ的导数用g=∂f(Θ)∂Θ=2Σm=1Nf(x(Θ)AmxH(Θ)-um)2Im{diag{x(Θ)}AmHxH(Θ)}]]>表示；
步骤b.初始化Θ0和一个极小值ε，设k＝0，S0＝I，计算g0；I为单位阵；
步骤c.设dk＝-Skgk，计算f(Θk+αkdk)并找出令其最小的αk；
设δk＝αkdk，Θk+1＝Θk+δk；
步骤d.如果||δk||≤ε，那么迭代结束，输出结果Θoptim＝Θk+1，由x(Θ)得到稀疏频率波形；如果||δk||≥ε，那么进行步骤e；
步骤e.计算gk+1，设γk＝gk+1-gk，根据拟牛顿法迭代计算原理，采用下式计算Sk+1，
Sk+1=Sk+(1+γkHSkγkγkHδk)δkδkHγkHδk-δkγkHSk+SkγkδkHγkHδk]]>
然后令k自加1，接着回到步骤c。
所述步骤(2)对稀疏频率雷达波形的序列进行特征分解，得到波形序列的主空间，具体为：
由于Sx是雷达波形与散射响应卷积过程的离散表示，若波形矢量s＝[s0 s1 ... sN-1]T，则Toeplitz矩阵S可以表示为：
Sx=sN-1sN-2...s00...00sN-1sN-2...s0...0.....................00...sN-1sN-2...s0x]]>
其中，x是一个长度为2N-1的向量，表示雷达散射的距离样本。
对S进行特征值分解：
SSH＝VΛVH
其中，V＝[v0 v1 ... vN-1]表示N个特征向量；Λ是包含N个与V对应的相关特征值的对角矩阵。V的前L个特征向量属于非主空间，后N-L个特征向量属于主空间(L的具体值由稀疏频率波形的稀疏度决定)。
所述步骤(2)由正交伪随机向量在其垂直子空间的投影得出通信波形，具体为：
设稀疏频率雷达波形的序列的长度为N。首先，生成一组K个伪随机N-1维向量表示为 bk，k＝1,2,...,K。由正交伪随机向量bk在其垂直子空间的投影得出通信波形样本，第一个通信信号c1可以设计为：
c1=(I-V~D,0V~D,0H)b1]]>
其中：
V~D,0=vLvL+1...vN-1]]>是V的主空间组成的一个N-L维特征向量。
要得到第二个通信信号，将c1与矩阵S结合得到一个N×2N矩阵：
SP,1＝[S c1]
用它进行新的特征值分解：
SP,1SP,1H=VP,1ΛVP,1H]]>
这样，第二个通信信号可计算得：
c2=(I-V~D,1V~D,1H)b2]]>
其中V~D,1=vL-1vL...vN-1]]>是VP,1主空间组成的N-L+1维特征向量。如此类推，第k个通信信号ck可由N×(2N+k-1)维矩阵SP,k-1的特征值分解得到。
ck=(I-V~D,k-1V~D,k-1H)bk]]>
其中V~D,k-1=vL-k+1vL...vN-1]]>是VP,k-1主空间组成的N-L+k-1维特征向量。这样就得到了一组K个通信信号。这种顺序设计得到的K个通信信号是相互正交的，且与雷达散射回波有一定的相关性。
所述步骤(3)为根据雷达嵌入通信模型，将设计的通信波形样本嵌入至所得的稀疏频率雷达波形，并用伪随机序列虚拟实际噪声干扰，仿真验证所设计波形的特性。
本发明与现有技术相比优点在于：
(1)、本发明中，稀疏频率雷达波形的频谱有若干离散阻带，与拥有连续完整频带的线性调频波形的雷达嵌入通信相比，通信信号能量可以更多地分布于稀疏频率波形的阻带中，因此可以提高通信信号的频带占用空间，从而获得更高的通信速率；
(2)、本发明中，K个通信信号样本按一定的排列(发送方与接收方都知道的规律)嵌入到每个雷达脉冲波形中，这K个通信信号样本之间相互正交，且与所设计的稀疏频率雷达散射回波具有一定的相关性，可以保证通信信号的低误码率和低截获率；
(3)、实现简单，并且由于添加通信波形信号后对雷达波形的自相关特性影响很小，因此可同时保证雷达系统的探测特性。
附图说明
图1为本发明基于雷达嵌入通信的稀疏频率波形设计方法的流程图；
图2为本发明基于雷达嵌入通信的稀疏频率波形设计方法的实施流程图；
图3为雷达嵌入通信的原理图；
图4为稀疏频率波形序列的特征值曲线图；
图5为添加通信信号前后的稀疏频率波形功率谱密度对比图；
图6为解相关接收机与拦截接收机误码率对比图；
图7为所设计波形的拦截度随特征值数变化图；
图8为添加通信信号前后的雷达波形归一化自相关函数对比图。
具体实施方式
本发明提供了一种雷达嵌入通信的稀疏频率波形设计方法，其基本思想是：如图1所示，首先，通过功率谱密度匹配方法设计出一种稀疏频率波形；然后对该稀疏频率波形序列进行特征值分解，得到特征向量组；最后用其中的一部分特征向量设计得出通信信号。雷达接收机接收雷达散射回波信号与雷达嵌入通信波形，实现雷达探测与通信两种功能，其中通信过程可以达到隐蔽的目的。
图2为本发明一种较佳实施方式的流程图，如图2所示，具体包含以下5个步骤：
1、基于雷达标识(Tag)系统应用，建立雷达嵌入通信模型。
在脉冲雷达探照范围内设置一个雷达响应标识(Tag)系统，接收雷达波形并进行再调制；再调制的结果是在回波中嵌入通信信号；Tag系统最后将混合波形发送至雷达接收机达到通信的目的。由于Tag系统中再调制过程的特殊性，通信过程可以达到隐蔽通信信息的目的。雷达接收机接收雷达散射回波信号与雷达嵌入通信波形，实现雷达探测与通信两种功能。雷达嵌入通信的原理如图3所示。
设Tag系统输出信号yr(t)为：
yr(t)＝αkck(t)+ys(t)+n(t) (1)
其中，n(t)为系统噪声；ys(t)表示Tag系统接收到的雷达散射回波；αk表示路径损失等对通信信号的混合影响；ck(t)为所设计的通信信号；k＝1,2,...,K表示第k个通信样本，K表示通信样本总数。同一时刻K个通信波形样本中的1个嵌入稀疏频率雷达波形。
用向量运算表示如下所示：
yr＝αkck+Sx+n (2)
其中，ck表示通信波形信号的向量表示；n表示噪声向量；x为一个长度为2N-1的向量表示的雷达散射的距离样本；Sx表示连续雷达散射模型，也是稀疏频率雷达波形与散射响应卷积过程的离散表示。
由于要保证通信信号的隐蔽性，通信信号ck(t)的功率远小于回波信号ys(t)的功率。这样导致了当通信信号ck(t)与回波信号ys(t)具有相关性时，用标准相干滤波器难以从混合信号中分离出通信样本。针对这一点，有了以下关于解相关接收机的设计的分析。
由于通信接收功能的相对简单性(K个通信波形样本只嵌入其中1个)且并不知道信号的相对功率(由αk决定)，因此通信接收功能通过使用最大似然接收机就可以实现(也称为解相关接收机)。设序列长度为N，生成一个N×(2N+K-1)维矩阵C：
C＝[S c1 ... cK] (3)
假定噪声符合以z＝[xT α1 ... αK]为参数，接收信号yr的概率密度函数可以表示为:
p(yr;b)=1πσv2exp{-1σv2(yr-Cb)H(yr-Cb)}---(4)]]>
对b的最大似然估计可以通过最小化(yr-Cb)H(yr-Cb)获得，则b的最大似然估计为：
b^=(CC)-1CHyr---(5)]]>
因为只要找出在K个通信信号中最可能的估计值，因此αk的具体数值并不重要。第k个解相关滤波器通过计算：
wk＝(CCH)-1ck (6)
得到w1,w2,...wK，与接收信号yr运算得到k的极大似然估计为：
k^=arg{maxk{|wkHyr|}}---(7)]]>
式(6)所示的最大似然估计接收机是对线性模型(1)中参数b的最小方差无偏估计器，式(7)则是多重假设检验的最小距离判决器。式(6)所示的最大似然估计接收机是对线性模型(1)中参数b的最小方差无偏估计器。
2、通过功率谱密度匹配及拟牛顿法设计稀疏频率雷达波形。
通过功率谱匹配的方法，可以得到特定通带阻带的稀疏频率波形。给定了功率谱密度u的分布后，稀疏频率波形设计的目标函数可以表示为：
minΘ||(Ax(Θ))⊗(Ax(Θ))*-u||2---(8)]]>
其中x为所设计波形，其相位矢量Θ＝[φ1…φN]；A为离散傅里叶变换矩阵，A中的元素Amn＝exp(-i×2π×n×m/Nf)，i为虚数单位，Nf与N相等；表示理想稀疏频率波形的功率谱密度。
目标函数式(8)是一个四次方非凸最优化问题，较佳地，考虑用拟牛顿法求解局部最小值取代全局最小值。拟牛顿法利用目标函数值和一阶导数的信息，构造出目标函数的曲率近似，使其具有收敛速度快的优点。在每一次迭代中需要计算目标函数的梯度。
设am＝[Am1…AmN]T，则式(8)可以写为：
f(Θ)=||(Ax(Θ))⊗(Ax(Θ))*-u||2=a1Hx(Θ)xH(Θ)a1-u1...aNfHx(Θ)xH(Θ)aNf-uNf2=Σm=1Nf|xH(Θ)amamHxH(Θ)-um|2=Σm=1Nf|xH(Θ)AmxH(Θ)-um|2---(9)]]>
其中，然后关于Θ的导数可以表示为：
g=∂f(Θ)∂Θ=2Σm=1Nf(x(Θ)AmxH(Θ)-um)2Im{diag{x(Θ)}AmHxH(Θ)}---(10)]]>
求解该目标函数的迭代运算步骤如下所示。
(1)、初始化Θ0和一个极小值ε，设k＝0，S0＝I，用上式计算g0；
(2)、设dk＝-Skgk，计算f(Θk+αkdk)并找出令其最小的αk，设δk＝αkdk，Θk+1＝Θk+δk；
(3)、如果||δk||≤ε，那么迭代结束，输出结果Θoptim＝Θk+1，得出稀疏频率波形；如果||δk||≥ε，那么进行第(4)步；
(4)、计算gk+1，设γk＝gk+1-gk，计算
Sk+1=Sk+(1+γkHSkγkγkHδk)δkδkHγkHδk-δkγkHSk+SkγkδkHγkHδk---(11)]]>
(5)、设k＝k+1，然后回到步骤(2)。
设置好适当的ε值，迭代运算结束后可以得到与理想功率谱密度分布u接近的结果。在得 到稀疏频率雷达波形后，分析在雷达波形中嵌入通信信号设计。
3、对稀疏频率雷达波形序列进行特征分解，得到波形序列的主空间
首先对已设计好的稀疏频率雷达波形序列进行特征分解，得到波形序列的主空间，再由正交伪随机向量在其垂直子空间的投影得出K个通信波形样本。
设序列长度为N，则稀疏频率雷达波形可以表示为：
s＝[s0s1...sN-1]T (12)
连续雷达散射模型可以离散地表示为稀疏频率雷达波形s与集合周围散射的卷积。这个卷积操作可以表示为矩阵乘积形式，如式(13)所示：
Sx=sN-1sN-2...s00...00sN-1sN-2...s0...0.....................00...sN-1sN-2...s0x---(13)]]>
其中：
x是一个长度为2N-1的向量，表示雷达散射的距离样本。在设计通信信号的时候，并不需要知道x的特定值。
对S进行特征值分解：
SSH＝VΛVH (14)
其中
V＝[v0v1...vN-1]包含N个特征向量；
Λ是包含N个相关特征值的对角矩阵。对特征值从小到大排列并作图，可得到如图4所示曲线。
图4中对应通带的较大特征值，与对应阻带的较小特征值，有一个明显的分界。较大特征值对应的特征向量组属于主空间；较小特征值对应的特征向量组属于非主空间。计算通信信号ck(t)时，要充分利用雷达波形的阻带部分频谱。假设雷达嵌入通信波形为K个样本，则通信样本可以通过投影法得到。
4、由正交伪随机向量在其垂直子空间的投影得出通信波形样本
设稀疏频率雷达波形序列的长度为N。首先，生成一组K个伪随机N-1维向量表示为bk，k＝1,2,...,K(Tag系统和接收机通信双方都知道伪随机向量bk的具体数值)。V的前L个特征向量属于非主空间，而剩下的N-L个特征向量属于主空间，L的具体数值与稀疏频率波形的稀疏度有关。这样第一个通信信号可以设计为：
c1=(I-V~D,0V~D,0H)b1---(15)]]>
其中：
V~D,0=vLvL+1...vN-1]]>是V的主空间组成的一个N-L维特征向量。
要得到第二个通信信号，c1要与矩阵S结合得到一个N×2N矩阵：
SP,1＝[S c1] (16)
用它进行新的特征值分解：
SP,1SP,1H=VP,1ΛVP,1H---(17)]]>
这样，第二个通信信号可计算得：
c2=(I-V~D,1V~D,1H)b2---(18)]]>
其中V~D,1=vL-1vL...vN-1]]>是VP,1主空间组成的N-L+1维特征向量。如此类推，第k个通信信号ck可由N×(2N+k-1)维矩阵SP,k-1的特征值分解得到：
ck=(I-V~D,k-1V~D,k-1H)bk---(19)]]>
其中V~D,k-1=vL-k+1vL...vN-1]]>是VP,k-1主空间组成的N-L+k-1维特征向量。这样就得到了一组K个通信信号。这种顺序设计得到的K个通信信号是相互正交的，且与雷达散射回波有一定的相关性。
得到通信样本后，用一个变量衡量通信的拦截性(对敌方接收机为拦截性；对友方接收机则为隐蔽性)。由于脉内雷达嵌入通信的操作是通过在雷达脉冲频谱附近嵌入一个通信信号实现，而这个频谱已经被雷达散射回波占据，因此基于测量谱能量容量的拦截概念不再适用。
而如果允许接收机获得一些先验知识，就可以计算一个与拦截可能性相关的变量。假设拦截接收机知道Tag系统和雷达接收机所使用的雷达波形，执行特征分解步骤。然后利用特征向量组，可以得到一个基于接收信号在非主空间投影的变量。
N×j维矩阵由j个特征向量组成，这j个特征向量对应j个最大的特征值。这样，表示列为j的一个主空间。对每一个j∈[1,...,N]，投影矩阵可以表示为：
Pj=I-V~D,jV~D,jH---(20)]]>
第j个投影矩阵Pj用于接收信号yr从而得到第j个投影残留为
z~j=Pjyr---(21)]]>
其中，第k个通信信号的单位化相关度，表示为ηk,j，定义如下
ηk,j=|ckHz~j|(ckHck)(z~jHz~j)---(22)]]>
其中，k为通信信号样本序号；j为获知的特征值数。这个η就是用于衡量通信信号拦截 概率的变量，称为拦截度。η的取值为0-1。越高的η值表示检测的通信信号与实际嵌入信号有越高的相似度，拦截的可能性越高；η趋于1时，拦截概率趋于1，这时可以认为通信信号不具备隐蔽性能。上文中用投影法计算得出的通信信号具有较好的隐蔽性，计算η值可以作为衡量这一点的指标。
5、根据雷达嵌入通信模型，将设计的通信波形样本嵌入至所得的稀疏频率雷达波形，并用伪随机序列虚拟实际噪声干扰，仿真验证所设计波形的特性。
可见，采用设计出的K个通信波形信号中的一个，嵌入所设计出的稀疏频率雷达波形；这种混合波形对Tag通信具有良好的隐蔽性，通信信号样本之间相互正交保证了一定的误码率、低拦截率；并且稀疏频率雷达波形在嵌入了通信信号后，仍保持良好的自相关特性。这些特征均有利于在雷达响应标识Tag系统上应用。
下面通过仿真的方法对本发明进行验证。
用于仿真的稀疏频率波形序列长度设为N＝200；设置理想功率谱密度u具有4个通带，3个阻带；通带设置为5dB，阻带设置为-30dB；ε值设为0.001；根据上文所示迭代步骤，得到稀疏频率波形的功率谱密度如图5所示。在得到稀疏频率波形的基础上，按照投影法公式添加通信信号，设通信样本总数K＝16；计算得出c1,c2,...,c16；图5中表示嵌入了通信样本c1后的功率谱密度；嵌入通信样本c2-c16的情况与此类似。
得到雷达嵌入通信波形后，分析误码率随信噪比的变化情况。本设计中有K＝16种通信波形，仿真中以第一种通信波形为例(即：发送通信信号c1，接收端可以成功判断c1的概率)。改变噪声n的功率，使信噪比从-25dB至0dB变化；设非主导特征向量数L＝80，设最大特征值数j＝120，每个点用10000个样本分别计算解相关接收机和拦截接收机，记录错误的概率并画出如图6所示曲线。图6表示解相关接收机公式计算画出的误码率曲线与拦截度计算画出的曲线。由图6可以看出，在信噪比为-8dB时，使用解相关接收机可以达到10-3的误码率，远远低于拦截接收机的情况。由此得出通信信号具有良好隐蔽性的结论。在图6中只画出了嵌入通信样本c1的情况，仿真中嵌入通信样本c2-c16的曲线与此类似，在信噪比SNR为-5dB至-10dB时，解相关接收机仍能保持较低的误码率。
进一步对拦截性进行分析。令信噪比为-15dB；计算拦截度时将最大特征值数j从0至200变化(即从0％至100％变化，因为N＝200)；每个点用200个样本计算ηk,j并得出平均值，画出如图7所示曲线。图7中横轴为特征值数j的百分比变化；纵轴为ηk,j的归一化值，越接近1表示拦截可能性越高。根据拦截度公式，拦截度的计算需要知道其中一种通信信号。设定已知通信样本c1，计算c2,c3,...,c16的归一化相关度；而不失一般性对c1也进行计算。分析图 7中曲线，c2-c16对应ηk,j的最大值小于0.2，c1对应ηk,j的最大值约为0.3。可以看出，即使增加特征值百分比，用该方法设计的波形仍能保持良好的防拦截能力。图7仿真中仅画出了假设已知通信样本c1运算的情况，实际仿真中以通信样本c2-c16运算的情况所作曲线与此类似，结论不变。
嵌入通信信号后，用自相关函数来分析通信信号对雷达性能的影响。图8中表示雷达波形加入通信信号前与加入通信波形后的归一化自相关函数(ACF)，纵坐标取对数表示。图8中添加通信波形前自相关函数具体数值为：积分旁瓣电平ISL＝-10.19dB，峰值旁瓣电平PSL＝-15.62dB；添加通信波形后自相关函数具体数值为：积分旁瓣电平ISL＝-10.94dB，峰值旁瓣电平PSL＝-16.44dB。可以看出，加入了通信信号c1后，自相关性能接近，说明了这种嵌入通信信号的设计方法并不影响雷达系统的探测性能。图8中只画出了加入通信样本c1的情况；加入通信样本c2-c16后自相关函数曲线与此类似，结论不变。
本发明未详细阐述部分属于本领域技术人员的公知技术。