禄晓飞, 靳硕静, 洪 灵, 戴奉周,*
(1. 中国人民解放军63620部队, 甘肃 兰州 732750; 2. 西安电子科技大学电子工程学院,陕西 西安 710071; 3. 陕西师范大学计算机科学学院, 陕西 西安 710061)
目标的微动是目标或其部分结构相对于雷达在径向上的非匀速运动,如振动、转动、加速运动等。微动在自然界中十分常见,如行人手脚的摆动,电机、直升机旋翼的旋转,桥梁的振动,弹道导弹目标的自旋、进动、章动等[1-4]。由微动引起的相位变化通常是非线性的,如目标自旋的微多普勒为正弦类曲线形式,目标进动的微多普勒为谐波函数形式,目标章动的微多普勒为更加复杂的时变函数形式。由雷达目标微动产生的微多普勒是时变的,这是由于在运动过程中，目标结构部件和其主体会产生相互作用[5-7]。由空间锥体目标微动产生的微多普勒特征具有独特性,能够反映出目标的精细特征[8-9],为雷达目标的检测与识别提供重要依据,具有重要的研究价值和意义。
时频分析方法作为微多普勒分析的一种有效手段,已被广泛应用。时频分析方法包括参数化时频分析和非参数化时频分析[10],在非参数化时频分析中,线性时频分析短时傅里叶变换(short-time fourier transform, STFT)方法的时频分辨率与窗长有关,而窗函数受Heisenberg不确定准则的限制,时间-频率分辨率相互制约,导致分辨率不高;双线性时频分布Wigner-Ville分布方法具有对称性、时移性、组合性等特征,可克服时频分辨率的Hisenbergur界,分辨率较高,但其对多分量信号或具有调制规律的信号会产生严重的交叉干扰,这是二次型时频分布的固有结果[11-12]。短时分数阶傅里叶变换[13]缓解了线性时频分析方法不确定性原理的局限性,计算简单,但在克服交叉项的干扰时,无法处理复杂信号。采用Hilbert-Huang变换获得了比传统Cohen类时频分析方法更好的性能,但其对噪声非常敏感。为了解决非参数化时频分析方法时频分辨率低、对噪声缺乏鲁棒性或交叉项严重等不足,有学者提出了基于参数化的时频分析方法。其中,块稀疏前后向时变自回归(block-sparse forward-backward time-varying autoregressive, BS-FBTVAR)模型利用模型参数中的块稀疏性对刚体目标进行微多普勒分析,利用了相对短时间窗内良好的频率分辨率,并将块稀疏贝叶斯学习(block-sparse Bayesian learning, BSBL)算法应用于模型的参数求解[14-16],不存在交叉项的影响。
然而，使用BSBL算法对BS-FBTVAR模型的求解并没有考虑到BS-FBTVAR模型中的块稀疏时不变系数与其邻域系数之间的相关性,与理想的时频分辨率有一定差距。在本文中,使用改进的BSBL(extended BSBL, EBSBL)算法[17]来求解BS-FBTVAR模型,通过对EBSBL模型中由人工构造的增广向量的积分形式进行分析,得到一个新的结构化稀疏先验，通过适当处理邻域的超参数来促进相邻稀疏系数之间的相关性。这种聚类先验采用了与(sparse Bayesian learning, SBL)理论中的模式耦合先验相似的公式[18-19]。刚体目标时变自回归(time-varying autoregre-ssive, TVAR)模型的时不变块稀疏系数的块边界是已知的,通过结合该先验信息,能够加快算法迭代的收敛速度。实验表明,所提算法的性能优于传统算法,不仅时频分辨率更高,抗噪声性能也较强。
本节介绍参数化时频分析方法中的BS-FBTVAR模型[16]。
前后向预测的TVAR模型用于分析非平稳随机信号的瞬时频率特性,弹道导弹目标的雷达回波数据为非平稳随机信号。设目标回波数据为x(n),可由p阶前后向TVAR模型表示。前后向预测的TVAR模型可分别表示为如下形式:
n=p,p+1,…,Nx-1
(1)
n=0,1,…,Nx-p-1
(2)
(3)
(4)
式中：“^ ”表示取估计值。
为了尽可能地利用观测数据,使前后向模型预测误差平均值最小,从而得到模型参数的解。根据式(1)和式(2)，可得前、后向预测误差分别为
n=p,p+1,…,Nx-1
(5)
n=0,1,…,Nx-p-1
(6)
进而,前、后向预测误差总和为
(7)
求解TVAR模型参数ai(n)(i=1,2,…,p)，需满足ξ最小。时变系数ai(n)不易直接求解,因此将前后向TVAR模型的时变系数用基函数与时不变常系数的线性组合进行表示,将TVAR模型中时变系数的估计转化为时不变块稀疏常系数向量的估计[20],进而使用压缩感知理论[21]对时不变常系数向量进行重构。将ai(n)由一组离散余弦基fj(n)(j=1,2,…,q)[22-25]展开为
(8)
式中:q为基函数的维数;aij(i=1,2,…,p；
j=1,2,…,q)均为常数。基函数的选择有多种方式,例如离散余弦(discrete cosine transform, DCT)基、离散傅里叶(discrete Fourier transform, DFT)基、Chebyshev基、多项式基等。由于在所使用的雷达目标回波数据中，TVAR系数在离散余弦基函数下展开的稀疏性最好,本文选用离散余弦基函数。将式(8)代入式(1)和式(2),可得
n=p,p+1,…,Nx-1
(9)
n=0,1,…,Nx-p-1
(10)
将式(9)和式(10)改写为矩阵形式:
Yf=-Zfβ+wf
(11)
Yb=-Zbβ*+wb
(12)
式中:Yf=[x(p),x(p+1),…,x(Nx-1)]T,Yb=[x(0),x(1),…,x(Nx-p-1)]T分别为TVAR模型的前向和后向观测数据。wf和wb分别为前、后向TVAR模型的高斯白噪声。Zf和Zb分别为TVAR模型的前向和后向观测矩阵
(13)
Zf=[A1,A2,…,Aq]
Zb=[B1,B2,…,Bq]
(14)
(15)
本节利用文献[16]中改进的EBSBL算法以及刚体目标块边界已知的先验条件来求解BS-FBTVAR模型。
2.1 β的聚类结构先验
EBSBL模型[17]用来解决未知块边界的块稀疏问题,该模型引入了一个扩展的隐藏块集zi,zi具有固定的块大小h,i=1∶g,块的数量g=Nβ-h+1,β的长度Nβ=pq，以及一个线性变换,将β表示为如下形式:
(16)
式中:Ei∈RNβ×h是一个除了第i行到第i+h-1行的部分被单位矩阵替代的零矩阵。分块向量z=[z1，z2，…，zg]是具有块大小h的等分割块稀疏信号。将式(16)代入式(11)和式(12),得到BS-FBTVAR模型的矩阵形式，如下所示:
(17)
式中:j∈{f,b1},f表示前向TVAR模型中用到的数据,b1表示后向TVAR模型中用到的数据。式(17)是一个具有已知块划分的块稀疏恢复问题,可由BSBL算法求解。由BSBL算法中假设的分块向量z服从p(z;{γi}i,B)=N(0,Σ0),Σ0=diag (γ1B,γ2B,…,γgB)来施加块稀疏性。每个块满足多元高斯分布:
p(zi;γi,B)=N(0,γi,B)
(18)
由于在本文的问题中,块稀疏系数向量β是块边界已知的,以p为块大小进行划分,因此h=p,则zi=[zi(0),zi(1),zi(p-1)]。根据式(16)对时不变系数向量β进行计算,则β的第i个元素为当前块zi和相邻的(p-1)个块中元素的和,考虑TVAR模型参数p的取值为奇、偶数的不同情况,总结得到
(19)
式(19)与模式耦合先验公式[18]相似。
2.2 求解BS-FBTVAR模型
利用β在式(18)中的聚类结构先验,本文的稀疏恢复算法[24]可以在贝叶斯框架下表述。考虑噪声为均值为0、方差为γ0的高斯白噪声,则式(11)和式(12)中信号模型的似然函数可以写为
(20)
(21)
(22)
超参数γ和γ0通过在系数β上边缘化，然后执行最大似然优化，从数据Yj中进行估计。边际分布可以由下式给出:
(23)
(24)
通过分析式(24)中的代价函数，并遵循文献[26]中的步骤,很容易证明式(24)的所有局部最小值都是作为稀疏解实现的。
为了获得超参数γ和γ0的估计,本文执行最大期望(expectation-maximization, EM)算法来等效地最大化联合分布p(Yj,β;γj,0,γj),可以将似然性和式(16)中的聚类结构先验组合为
p(Yj,β;γj,0,γj)=p(Yj|β;γj,0)p(β;γj)
(25)
在EM算法中,β可以被视为隐藏变量,可以最大化Eβ|Yj,γj,0,γj{p(Yj,β;γj,0,γj)}来找到γj,0和γj的估计值。其中，Eβ|Yj,γj,0,γj{·}表示关于β的后验分布的期望。
将式(24)代入Eβ|Yj,γj,0,γj{p(Yj,β;γj,0,γj)},本文可以通过最大化Eβ|Yj,γj{lnp(β;γj)}等效地估计γj:
Eβ|Yj,γj{ lnp(β;γj)}=
(26)
(27)
γj,0的更新形式如下:
(28)
由于前后向预测都是在恢复β,需要做如下约束:
(29)
(30)
(31)
将目标的时不变系数向量β估计出来后,通过离散余弦基的线性组合公式将β转换为时变系数ai(n),利用瞬时功率公式计算得到目标的瞬时功率谱,由式(3)给出。
2.3 模型参数的选取准则
BS-FBTVAR模型中的模型阶次p、基函数个数q与雷达回波数据的散射点个数和噪声环境有关,p和q的选取需要采用一定的准则,共包含p×q个未知量。本文使用最小描述长度(minimum description length, MDL)准则来实现p与q的选择[29]:
MDL (p,q)=
(32)
本节分别在电磁仿真工具软件产生雷达目标回波的情况下和使用目标实测回波数据的情况下对目标进行微多普勒分析。
3.1 仿真产生雷达回波数据的微多普勒分析结果
为了验证所提出的基于聚类先验的BS-FBTVAR模型求解算法的有效性,本文采用电磁仿真软件CST Studio Suite 2020产生空间锥体目标的雷达回波。目标使用锥体模型,锥顶为半球体,锥体材质为理想良导体。锥体参数设置如下:高度为1.5 m,底面半径为0.9 m,在距离锥体底面0.5 m处的锥体表面开出等间隔的半球体凹槽,半球体半径为0.05 m。仿真雷达目标回波所用的雷达参数设置如下:载频为10 GHz,带宽为2 GHz,脉冲重复频率(pulse repetition frequency, PRF)为800 Hz,驻留时间为1 s。目标的自旋频率为1 rad/s,进动频率为1 rad/s,雷达入射方向与锥体中心轴线垂直,模型图如图1所示。
图1 目标模型示意图Fig.1 Schematic diagram of the target model
图2所示为仿真雷达回波数据所对应的真实的时变功率谱。其中，两条正弦曲线对应目标的顶点A和顶点B,A为点散射中心,B为滑动散射中心。C点只有在雷达视线范围内才能被观测到,在图中是不可见的。中间的3条正弦曲线对应3个凹点D、E、F,即点散射中心。由于目标在旋转,因此空间锥体目标等间隔的3个凹槽会被雷达视线轮流观测,在时频图中呈现出不完整的曲线。
图2 真实的微多普勒时频图示意图Fig.2 Schematic diagram of real micro-Doppler time-frequency diagram
在实验中,将超参数设置为a0=1e-2,b0=1e-3,c0=1e-2,d0=1e-3,通过MDL准则得到p和q。在噪声较大时,p的取值会大于目标含有的散射点的个数。将收敛阈值设置为1e -3,对目标在无噪声、20 dB噪声、10 dB噪声环境中分别使用STFT时频分析方法、基于BSBL求解BS-FBTVAR模型的时频分析方法、本文算法进行微多普勒分析,得到目标回波的时变功率谱估计结果,设置功率谱幅度的动态范围一致,均为120 dB,如图3所示。
图3 仿真目标回波在无噪声、20 dB噪声以及10 dB噪声环境中使用3种时频分析算法的时变功率谱估计示意图Fig.3 Schematic diagram of time-varying power spectrum estimation of the use of three time-frequency analysis algorithms in noise-free and noise environment (20 dB and 10 dB)
从图2与图3可以看出,本文算法较为准确地估计出了目标回波的时变功率谱。STFT时频分析方法整体分辨率较低,对于目标模型中的D、E、F点散射中心而言,当处于即将被遮挡的状态时,回波信号能量较低,产生的微多普勒分辨率更低。对于由散射中心A、B微动导致的微多普勒特征模糊的问题,无法同时分辨较近的两条交叉曲线。基于BSBL求解BS-FBTVAR模型的时频分析方法的时频分辨率精度有所提高,能够将由目标模型滑动散射中心和点散射中心自旋和进动产生的微多普勒频率完整地分析出来,但是在对点散射中心D、E、F的微动进行分析时,在正弦曲线的各个周期的上顶点和下顶点处均出现了错误分量,时频图较为浑浊。本文算法能够将由空间锥体目标的滑动散射中心和点散射中心微动产生的微多普勒频率完整地提取出来,明显分辨出距离较近的两条交叉曲线,时频图更为清晰,时频分辨率更高。
由图3可以看出,在目标散射点处于即将被遮挡的状态或者低信噪比情况下,STFT时频分析方法时频能量较弱。而基于TVAR模型的时频分析方法时频能量较强,清晰可见,这是因为STFT方法对回波数据进行滑窗处理,只处理短时窗内的部分回波数据,而基于TVAR模型的时频分析方法对全部回波数据进行处理。
采用Renyi熵[30]对3种时频分析算法的时频分辨性能进行进一步量化分析。Renyi熵是描述系统混乱度的一种衡量指标,熵值越低，代表时频谱图的能量聚集性越好,其计算公式为
(33)
式中：P(t,f)为瞬时功率谱值;α为Renyi熵的次序。本文使用3次Renyi熵对时频聚集性进行评价,对噪声环境为0～20 dB的目标回波分别使用3种时频分析方法计算时频谱图的Renyi熵,取100次实验的平均结果如图4所示。由图4可以看出,基于BSBL求解BS-FBTVAR模型算法与本文算法的时频能量聚集性都明显高于STFT时频分析方法,本文算法的Renyi熵值在0～5 dB噪声环境下与基于BSBL的BS-FBTVAR的时频分析算法基本一致,在5～20 dB噪声情况下最低,时频聚集性更好,抗噪性能良好。
图4 目标回波的3种时频谱图的Renyi熵值Fig.4 Renyi entropy values of three time-frequency spectrograms of target echoes
3.2 利用实测数据对目标进行微多普勒分析
雷达参数设置如下:载频为15 GHz,带宽为2 GHz,PRF为200 Hz。目标使用光滑锥体,高为0.85 m,底面半径为0.225 m,中部和底部分别贴有2个半径为0.02 m的半球体,两个半球体对称分布,如图5所示。
图5 目标示意图Fig.5 Target schematic diagram
锥体目标自旋频率为1.5 Hz,雷达发射电磁波照射锥体目标并接收其回波,对其进行解线调、快时间傅里叶变换以及去除剩余视频项和包络斜置项处理,得到目标的高分辨距离像(high resolution range profile, HRRP)数据。取HRRP数据中锥顶处、中部贴有半球体处以及锥体底部距离单元内的数据,分别使用STFT、BSBL求解BS-FBTVAR模型算法和本文算法，进行微多普勒分析,得到目标回波的时变功率谱估计,设置功率谱幅度的动态范围一致,均为120 dB,如图6所示。图6(a)～图6(c)为对目标的锥顶距离单元进行微多普勒分析的结果。该距离单元含有1个点散射中心,在目标进行自旋时,此散射中心不会产生微多普勒频率,时频图为零频。对比3种算法可以看出，STFT算法和基于BSBL求解BS-FBTVAR模型算法都在零频附近的较大范围内出现了频谱能量,而本文算法在零频的能量聚集较为集中,时频分辨率最高。图6(b)～图6(d)为对目标的中部贴有半球体处的距离单元进行微多普勒分析的结果。该距离单元包含2个对称分布的点散射中心,自旋产生的微多普勒为两条正弦曲线。图6(g)～图6(i)为对目标的底部距离单元进行微多普勒分析的结果,包含1个滑动散射中心和2个点散射中心,2个点散射中心产生的微多普勒为两条正弦曲线,滑动散射中心不会产生微多普勒频率,在时频图中表现为零频[31]。
图6 目标实测回波数据的时变功率谱估计示意图Fig.6 Schematic diagram of time-varying power spectrum estimation of target measured echo data
由图6可以看出,在对目标的实测回波数据进行微多普勒分析时,本文算法能够准确估计出目标微动产生的微多普勒频率,并且较STFT算法和基于BSBL求解BS-FBTVAR模型的时频分析算法，时频图清晰度更高。
由Renyi熵值也可以表明本文算法的优越性,表1是对实测回波数据3种算法的熵值对比。由表1可以看出,本文算法Renyi熵值最低,时频聚集性最好。
表1 实测数据3种算法Renyi熵定量分析对比
本文提出了一种新的求解BS-FBTVAR模型的方法,进行目标的微多普勒分析,利用改进EBSBL算法,通过人工构造增广向量模型的潜在关系,得到一个结构化稀疏先验,然后考虑其超参数和其相邻超参数之间的相互作用,得到聚类结构先验,并施加刚体目标块边界已知的先验信息来求解BS-FBTVAR模型,从而对目标进行时变功率谱估计。实验表明所提算法性能优于传统的时频分析算法,给出的时频图时频分辨率更高,多普勒曲线的聚集性也更好,并且抗噪声性能良好。
猜你喜欢锥体时频时变搜集凹面锥体小猕猴智力画刊(2020年5期)2020-06-01锥体上滚实验的力学分析物理实验(2019年4期)2019-05-07基于时变Copula的股票市场相关性分析智富时代(2017年4期)2017-04-27基于时变Copula的股票市场相关性分析智富时代(2017年4期)2017-04-27烟气轮机复合故障时变退化特征提取广东石油化工学院学报(2016年6期)2016-05-17进动锥体目标平动补偿及微多普勒提取系统工程与电子技术(2016年2期)2016-04-16基于MEP法的在役桥梁时变可靠度研究中国铁道科学(2015年4期)2015-06-21基于时频分析的逆合成孔径雷达成像技术舰船科学技术(2015年8期)2015-02-27对采样数据序列进行时频分解法的改进电测与仪表(2014年17期)2014-04-04电针针刺锥体区即时镇痛发作期偏头痛218例中国中医药现代远程教育(2014年23期)2014-03-01