Loading [MathJax]/jax/element/mml/optable/BasicLatin.js

基于矩阵填充的随机步进频雷达高分辨距离-多普勒谱稀疏恢复方法

胡雪瑶 梁灿 卢珊珊 王在洋 郑乐 李阳

赵博, 陈基, 黄磊. 单比特多模态SAR干扰方法研究[J]. 雷达学报, 2022, 11(6): 1119–1130. doi: 10.12000/JR22176
引用本文: 胡雪瑶, 梁灿, 卢珊珊, 等. 基于矩阵填充的随机步进频雷达高分辨距离-多普勒谱稀疏恢复方法[J]. 雷达学报(中英文), 2024, 13(1): 200–214. doi: 10.12000/JR23176
ZHAO Bo, CHEN Ji, and HUANG Lei. One-bit multi-modality jamming method against SAR[J]. Journal of Radars, 2022, 11(6): 1119–1130. doi: 10.12000/JR22176
Citation: HU Xueyao, LIANG Can, LU Shanshan, et al. Matrix completion-based range-Doppler spectrum estimation for random stepped-frequency radars[J]. Journal of Radars, 2024, 13(1): 200–214. doi: 10.12000/JR23176

基于矩阵填充的随机步进频雷达高分辨距离-多普勒谱稀疏恢复方法

DOI: 10.12000/JR23176
基金项目: 国家自然科学基金(62388102),国家重点研发计划(2018YFE0202101, 2018YFE0202103)
详细信息
    作者简介:

    胡雪瑶,博士,副研究员,硕士生导师,主要研究方向为毫米波雷达系统设计与雷达信号处理技术等

    梁 灿,博士生,主要研究方向为毫米波雷达系统设计与雷达信号处理技术等

    卢珊珊,硕士,工程师,主要研究方向为雷达信息研究等

    王在洋,硕士生,主要研究方向为毫米波雷达信号处理技术等

    郑 乐,博士,教授,博士生导师,主要研究方向为毫米波雷达信号处理与数据处理技术等

    李 阳,博士,研究员,博士生导师,主要研究方向为宽带雷达系统、雷达信号处理与雷达系统设计等

    通讯作者:

    郑乐 le.zheng@bit.edu.cn

  • 责任主编:兰岚 Corresponding Editor: LAN Lan
  • 中图分类号: TN958.6

Matrix Completion-based Range-Doppler Spectrum Estimation for Random Stepped-frequency Radars

Funds: The National Natural Science Foundation of China (62388102), The National Key R&D Program of China (2018YFE0202101, 2018YFE0202103)
More Information
  • 摘要: 随机步进频雷达通过合成大带宽,能在较低硬件复杂度下获得距离高分辨效果,同时由于其每个脉冲的载频随机捷变,因而具有强的抗干扰、电磁兼容能力,在复杂电磁环境高精度探测领域具有重要的应用价值。然而,由于其波形在时频域稀疏的感知形式,造成回波相参信息有所缺失,因而传统匹配滤波方法在估计高分辨距离-多普勒时会演化为欠定估计,导致估计谱中产生起伏高旁瓣,严重影响探测性能。为此,该文提出一种基于Hankel重构矩阵填充的随机步进频雷达高分辨距离-多普勒谱低旁瓣稀疏恢复方法。该方法采用低秩矩阵填充思想补全波形在时频域稀疏感知时造成的缺失采样,恢复目标连续相参信息,可以有效解决欠定估计问题。文章首先构建了随机步进频雷达的慢时间-载频(时-频)回波欠采样数据矩阵;然后,重构待恢复数据矩阵为双重Hankel型,并分析证明了矩阵满足低秩先验特性;最后,利用ADMM算法补全未采样时频数据,恢复相参信息,保证了高分辨距离-多普勒谱低旁瓣稀疏恢复。仿真和实测试验证明了该文所提方法的有效性和优越性。

     

  • 近些年,频率步进(Stepped Frequency, SF)波形被广泛应用于雷达高精度探测领域[14]。它由一系列载频线性步进脉冲信号组成,通过相参合成这些信号获得大带宽效果,可以大幅提升雷达距离分辨率,同时避免因发射瞬时大宽带信号增加的硬件系统复杂度[5]。频率步进波形的一种改进方式是它的载频步进值随机改变,通常也被称之为随机步进频(Random Stepped Frequency, RSF)波形[6]。与常规SF波形相比,RSF波形具有两个显著的优点:首先,随机跳频的方式可以大幅降低信号被截获或“碰撞”概率,从而提升雷达抗干扰能力;其次,由于载频非线性变化,RSF波形解耦了高分辨距离和多普勒相位,产生图钉型模糊函数,从而获得了估计目标多普勒速度的能力[7]。然而,由于RSF波形对环境的感知形式在慢时间-载频(时-频)是稀疏的,因而雷达在时-频域对回波信号的采样是随机欠采样,造成相参信息缺失,传统匹配滤波(Matched Filter, MF)相参处理会演化为欠定估计,导致高分辨距离-多普勒谱中产生起伏高旁瓣[8,9]。此时,低信噪比目标可能被高信噪目标的旁瓣掩盖,导致漏检;或者起伏波动的旁瓣可能被检测为假目标,产生虚警。显然,上述问题限制了RSF波形的广泛应用。

    为了解决该问题,文献[1012]提出使用失配滤波器方法减轻高旁瓣对参数估计影响,但这类方法会降低目标信噪比或雷达距离分辨率。文献[13]提出使用RSF回波的包络和相位设计用于高分辨率距离像合成的旁瓣抑制滤波器。然而该方法只能消除单个多普勒频率下的起伏旁瓣;并且,当不能预先获得目标的多普勒信息时,其旁瓣抑制性能仍有待提高。从波形设计的角度出发,可以通过设计跳频序列来抑制旁瓣[1416],然而这些方法并没有考虑到多目标应用场景,因此所设计的波形不能保证多目标的旁瓣在相互叠加后仍保持在低功率水平。

    随着压缩感知(Compressed Sensing, CS)技术的兴起,该技术已被用于RSF波形中解决因随机欠采样造成的欠定估计问题[1720]。其本质是利用回波信号数据中的稀疏性,通过设计特定的测量矩阵寻找信号在高分辨距离-多普勒域的稀疏表示,实现低旁瓣稀疏恢复。但CS恢复性能在很大程度上受设计的测量矩阵制约,当测量矩阵的稀疏基与场景中目标参数的分布不完全匹配时(又称“基不匹配”),恢复性能会受到显著影响[21]。文献[22]提出结合机器学习的压缩感知方法,通过网络自适应学习获得更适应场景特征的测量矩阵;然而,该方法的训练过程需要海量的样本数据,且CS方法在高维度空间中使用一维矢量运算,导致计算复杂度激增。

    矩阵填充(Matrix Completion, MC)[23,24]是一种通过使用其部分有限观测数据在低秩约束下恢复矩阵中缺失数据的技术。文献[2527]中已展示了MC在多输入多输出(Multiple Input Multiple Output, MIMO)雷达中的有效性和优势,其已被用于恢复存在随机稀疏轨迹的阵列方位或俯仰维的欠采样数据,保证多维稀疏成像质量。此外,MC技术也被证明在处理高维数据时具有效率高、恢复精度高等优势[28]。对于经过快时间脉冲压缩处理的RSF回波信号,由于其沿时-频域呈现二维随机欠采样特性,因而时-频维度中的回波数据也有所缺失,因此可以使用MC技术恢复未采样的数据。

    基于上述分析,本文提出了一种基于Hankel重构矩阵填充的随机步进频雷达高分辨距离-多普勒谱估计方法,实现低旁瓣谱稀疏恢复。首先,推导并建立了随机步进频雷达单个粗分辨距离单元内回波时-频欠采样模型;并利用Hankel矩阵变换重构待恢复数据矩阵,证明了矩阵满足低秩先验特性,为基于矩阵填充的谱估计算法奠定了基础。随后,将低秩矩阵缺失数据恢复问题转化为以填充后矩阵截断核范数最小化为目标、以填充后矩阵原采样位置元素不变为约束的优化模型,进而采用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)求解模型,实现时-频数据矩阵填充,恢复缺失相参信息,解决匹配滤波欠定估计问题。所提方法直接以矩阵形式进行二维稀疏建模与求解,不仅有效避免了基不匹配问题,在保持良好重构性能的同时还进一步降低了运算量,极大地提升了运算速度。仿真和实测数据试验验证了方法的有效性和优越性。

    本文以调频RSF雷达为例,其发射信号由一系列脉内相同线性调频、脉间不同载频伪随机调制脉冲组成,信号预处理流程如图1所示。其中,第n个脉冲的载频为fn=f0+UnΔff0为初始频率,Un表示第n个脉冲(共发射N个脉冲)的载频编码,且Un{1,2,,M}M为随机步进频率总点数,Δf 为最小频率步进间隔,MΔf为波形合成带宽。因此第n个脉冲的发射信号可以表示为

    图  1  调频RSF信号预处理流程图
    Figure  1.  Diagram of the RSF chirp signal pre-processing flow
    st(t)=s(t)exp[j2πfn(tnT)]=rect(tnTTp)exp[j2πB2Tp(tnT)2]exp[j2πfn(tnT)] (1)

    其中,s(t)表示雷达基带信号,t表示脉内快时间,B为脉内调频带宽,Tp为发射调频信号的脉宽,T为脉冲重复间隔。

    假设发射信号经理想点目标调制后幅度为A,回波信号可表示为

    sr(t)=Arect(tnTτTp)exp[j2πB2Tp(tnTτ)2]exp[j2πfn(tnTτ)] (2)

    其中,τ表示回波信号的双程时延。基于“走-停”模型假设[29,30],时延τ可以近似表示为

    τ2(rvnT)c (3)

    其中,r, v分别为目标的初始距离和相对雷达的径向速度,c为光速。

    回波信号经模数转换、下变频、滤波处理后可表示为

    sr(l,n)=Aexp{j2π[B2Tp(ltsnTτ)2fnτ]} (4)

    其中,l=0,1,,L1为脉内快时间采样索引,L为一个脉冲内的快时间采样总点数,ts代表快时间采样间隔。

    信号经快时间脉冲压缩处理后可表示为

    y(l,n)=sr(l,n)ˉs(ll,n)αsinc(l2rBc)exp(j4πUnΔfrc)exp(j4πvnTf0c) (5)

    其中,表示卷积运算,ˉs(ll,n)为式(1)中雷达基带信号s(t)的共轭转置,α=Aexp(jϕC)表示经脉冲压缩后目标的复振幅,其中ϕC为脉冲压缩后回波信号中的相位常数项,sinc(·)表示辛格函数,为脉冲压缩后回波信号,l则表示经脉冲压缩后的粗分辨距离单元(Coarse-Resolution Range Cell, CRRC)索引。需要注意的是,由于f0Δf且本文关注的目标速度远小于光速,因此式(5)中忽略了相位交叉项。

    提取不同脉冲目标所在CRRC(如图1所示)信号,回波时-频采样模型可以表示为

    y(n)=αexp(j4πUnΔfrc)exp(j4πvnTf0c) (6)

    将上述点目标回波模型推广至多目标情景,可以得到RSF雷达时-频回波信号一般表达式为

    y(n)=Kk=1αkexp(j4πUnΔfrkc)exp(j4πvknTf0c) (7)

    其中,K表示一个CRRC内的目标数量,αk, rkvk分别表示第k个目标的复振幅、初始距离和相对径向速度。

    考虑式(7)对目标αk, rkvk的估计本质是对整个观测谱的重建,即估计每个高分辨距离-速度单元内的散射强度,再通过对应的栅格位置确定估计目标的幅度、距离、速度参数等。因此,完成MN个单元散射强度的重建即可获得目标参数rkvk。显然,网格的格点数取决于重建观测谱的精细程度,精细程度越高,需要的网格越密,格点个数越多。考虑网格间距要小于或等于波形距离、速度分辨率,因此设置的距离和速度单元分别最少为M, N个。综上,待估计的高分辨距离和多普勒参数数量(M×N)远大于采样数据量(N),即RSF回波的两维随机欠采样问题。此时,若直接采用传统匹配滤波方法进行相参处理,参数估计会演化为欠定估计问题,导致谱中高旁瓣基底产生。针对上述问题,本文提出了一种基于矩阵填充的方法补全缺失的采样数据,避免欠定估计。

    首先,将式(7)中的回波信号重新排列为矩阵形式。流程如图2所示,根据载频编码顺序将采样数据重新排列至相应的时频采样位置,空白采样位置处设置为0。此时,在重新排列的欠采样矩阵Y中,第n行第m列元素可以表示为

    图  2  回波信号矩阵填充处理流程图
    Figure  2.  Matrix completion flow
    (Y)n,m={y(n), m=Un0, mUn (8)

    上述欠采样矩阵也可以视作是来自于完备矩阵的随机稀疏观测,即:

    (Y)n,m={(Δ)n,m, (n,m)Ω0, (n,m)Ω (9)

    其中,Δ为等效时频满采样矩阵,Ω为载频编码模式Un对应的时频数据集合。本文的目的即通过欠采样矩阵Y恢复满采样矩阵Δ

    基于式(7),由第k个目标构成的时频满采样矩阵中第n行第m列元素可以表示为

    (Δk)n,m=αkexp(j2π(nT2vkf0c+mΔf2rkc)) (10)

    显然,式(10)可以改写为两个向量的乘积:

    Δk=αkad(vk)ar(rk)H (11)

    其中,两个向量分别为第k个目标关于速度和高分辨率距离的导向矢量,分别为

    ad(vk)=[1, exp(j4πvkf0Tc), , exp(j4πvkf0(N1)Tc)]H (12)
    ar(rk)=[1, exp(j4πrkΔfc), , exp(j4πrk(M1)Δfc)]H (13)

    根据矩阵乘积秩定理[31],两个矩阵相乘得到新矩阵,其秩不大于相乘前任意一个矩阵的秩,因此上述第k个目标的时频满采样矩阵秩为1。进一步地,由K个目标组成的等效时频满采样矩阵可以视为由单个目标组成的满采样矩阵的线性组合,因此可以得到:

    Δ=Kk=1Δk=Kk=1αkad(vk)ar(rk)H= XdΛXHr (14)

    其中,

    Xd=[ad(v1) ad(v2)  ad(vK)]Xr=[ar(r1) ar(r2)  ar(rK)]Λ=diag{[α1 α2  αK]} (15)

    根据矩阵秩的不等式性质[32]可知,由K个目标组成的时频满采样矩阵的秩满足:

    rank(Δ)=rank(Kk=1Δk)Kk=1rank(Δk)=K (16)

    矩阵填充方法要求待恢复矩阵满足低秩特性。为了保证该前提,本文引入了双重Hankel变换重新构造待恢复矩阵,从而增强回波信号等效满采样矩阵的低秩特性。具体地,首先定义一个维度为N1×(NN1+1)的Hankel矩阵:

    ΔH=[Δ1Δ2ΔNN1+1Δ2Δ3ΔNN1+2ΔnΔN1ΔN1+1ΔN] (17)

    其中,第n个块矩阵是一个维度为M2×(MM2+1)的Hankel矩阵,该矩阵由矩阵Δ的第n行所组成,即:

    Δn=[(Δ)n,1(Δ)n,2(Δ)n,MM2+1(Δ)n,2(Δ)n,3(Δ)n,MM2+2(Δ)n,m(Δ)n,M2(Δ)n,M2+1(Δ)n,M] (18)

    其中,N1M2分别是构造一阶和二阶Hankel变换矩阵的维度参数。

    对具有双重Hankel形式的时频观测矩阵ΔH进行分解可以得到:

    ΔH=[Xr,LD1dXr,LDndXr,LDN1d]Λ[D1dXHr,RDNn+1dXHr,RDNN1+1dXHr,R]=˜XdΛ˜Xr (19)

    其中,

    {Xr,L=(Xr)1:M2,:Dnd=diag{(Xd)n,:}DNn+1d=diag{(Xd)Nn+1,:}Xr,R=(Xr)1:MM2+1,: (20)

    显然,由式(20)可知:

    rank(ΔH)=rank(Δ)K (21)

    而经Hankel构造变换后,矩阵的维度可达N1M2×(NN1)(MM2+1),因此变换后的双重Hankel矩阵具有更优异的低秩特性。

    根据第2节的分析,单个粗分辨距离单元目标信号具有稀疏性质,且经Hankel重构后的时频回波矩阵还具有了低秩特性。

    因此结合矩阵稀疏-低秩特性,利用矩阵填充算法可以恢复存在数据缺失的观测矩阵YH(由原始观测矩阵经双重Hankel构造变换得到),进而得到时频满采样矩阵估计值ˆΔH,建立下述优化模型:

    minˆΔHrank(ˆΔH)s.t.YH=PΩ(ˆΔH) (22)

    其中,PΩ是观测投影操作,可以表示为

    (YH)p,q={(ˆΔH)p,q, (p,q)ΩH0, (p,q)ΩH  (23)

    其中,ΩH表示矩阵经双重Hankel变换后数据所在位置的集合,(·)p,q表示矩阵第p行、第q列元素,其中,P=N1M2, Q=(NN1+1)(MM2+1)

    由于低秩问题的求解是NP-hard的,通常通过核范数来近似表示矩阵的秩函数,即低秩优化问题可以建模为

    minˆΔHˆΔHs.t.(YH)p,q=(ˆΔH)p,q,(p,q)ΩH (24)

    针对上述优化问题,现有算法通常将其定义为半正定规划问题进而采用内点法求解,然而由于双重Hankel变换后的矩阵维度较大,导致半正定规划求解十分困难;另一方面,核范数松弛是对所有奇异值的和求最小化处理,难以准确近似矩阵的秩[33]。对此,本文采用截断核范数代替传统核范数逼近矩阵的秩。截断核范数优化目标仅求解与矩阵秩不相关的较小奇异值的和,间接控制待恢复矩阵秩,从而确保了矩阵低秩估计;进一步地,本文通过交替迭代求解思想,将截断核范数优化问题转换为凸优化问题实现求解。

    截断核范数优化模型如下:

    minˆΔHˆΔHrs.t.(YH)p,q=(ˆΔH)p,q,(p,q)ΩH (25)

    其中,ˆΔHr表示矩阵ˆΔH以第r个奇异值作为截断后的核范数,可以表示为

    ˆΔHr=i=r+1min (26)

    其中, {\sigma _i}( \cdot ) 表示取矩阵的第i个奇异值。

    由于矩阵的截断核范数不遵循凸性质,因此我们难以直接求解式(25)。对此,我们采用如下变换将原问题转化为凸问题进行求解。

    首先,对待估计的双重Hankel时频满采样矩阵 {\hat{\boldsymbol{\varDelta}} _{\mathrm{H}}} 进行奇异值分解可得到:

    \left[ {{\boldsymbol{U}},{\boldsymbol{\Sigma}} ,{\boldsymbol{V}}} \right] = {\mathrm{svd}}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right) (27)

    其中,

    {\boldsymbol{U}} = \left( {{{\boldsymbol{u}}_1},{{\boldsymbol{u}}_2}, \cdots ,{{\boldsymbol{u}}_P}} \right) \in {\mathbb{R}^{P \times P}} (28)
    {\boldsymbol{V}} = \left( {{{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2}, \cdots ,{{\boldsymbol{v}}_Q}} \right) \in {\mathbb{R}^{Q \times Q}} (29)

    对左奇异向量U和右奇异向量V取前r行以实现核范数截断的效果:

    {\boldsymbol{A}} = \left( {{{\boldsymbol{u}}_1},{{\boldsymbol{u}}_2}, \cdots ,{{\boldsymbol{u}}_r}} \right) \in {\mathbb{R}^{r \times P}} (30)
    {{\boldsymbol{B}}^{\mathrm{H}}} = \left( {{{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2}, \cdots ,{{\boldsymbol{v}}_r}} \right) \in {\mathbb{R}^{Q \times r}} (31)

    其中,(·)H表示矩阵共轭转置操作,r ≤ min( P,Q )为截断奇异值个数。这里需要说明的是,矩阵ABH在后续算法迭代求解过程中会随着 {\hat {\boldsymbol{\varDelta}} _{\mathrm{H}}} 的更新而更新。

    进一步,根据迹不等式可知:

    \begin{split} {\mathrm{Tr}}\left( {{\boldsymbol{A}}{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}{{\boldsymbol{B}}^{\mathrm{H}}}} \right) \, &= {\mathrm{Tr}}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}{{\boldsymbol{B}}^{\mathrm{H}}}{\boldsymbol{A}}} \right) \\ & \le \sum\limits_{i = 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\Bigr( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \Bigr){\sigma _i}\Bigr( {{{\boldsymbol{B}}^{\mathrm{H}}}{\boldsymbol{A}}} \Bigr)} \end{split} (32)

    由于AB的秩为r,因此矩阵BHA的秩≤ r。由此可以得到:

    \begin{split} & \sum\limits_{i = 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right){\sigma _i}\Bigr( {{{\boldsymbol{B}}^{\mathrm{H}}}{\boldsymbol{A}}} \Bigr)} \\ & \quad= \sum\limits_{i = 1}^s {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right){\sigma _i}\Bigr( {{{\boldsymbol{B}}^{\mathrm{H}}}{\boldsymbol{A}}} \Bigr)} \\ & \qquad + \sum\limits_{i = s + 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right){\sigma _i}\Bigr( {{{\boldsymbol{B}}^{\rm H}}{\boldsymbol{A}}} \Bigr)} \\ & \quad= \sum\limits_{i = 1}^s {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right) \cdot 1} + \sum\limits_{i = s + 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right) \cdot 0} \\ & \quad = \sum\limits_{i = 1}^s {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} \le \sum\limits_{i = 1}^r {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} \end{split} (33)

    综上,可以得到:

    \begin{split} {\mathrm{Tr}}\left( {{\boldsymbol{A}}{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}{{\boldsymbol{B}}^{\rm H}}} \right)\, & \le \sum\limits_{i = 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right){\sigma _i}\Bigr( {{{\boldsymbol{B}}^{\mathrm{H}}}{\boldsymbol{A}}} \Bigr)} \\ & \le \sum\limits_{i = 1}^r {{\sigma _i}\left( {{{\hat {\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} \end{split} (34)

    基于以上分析并参考文献[34,35],我们通过引入中间变量W= {\hat {\boldsymbol{\varDelta}} _{\mathrm{H}}} ,可将式(25)中的目标函数转化为如下形式:

    \begin{split} & \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} {\text{ }}{\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right\|_r} \\ &\Leftrightarrow \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} {\text{ }}\sum\limits_{i = r + 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} \\ & \Leftrightarrow \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} {\text{ }}\sum\limits_{i = 1}^{\min \left( {P,Q} \right)} {{\sigma _i}\left( {{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} - \sum\limits_{i = 1}^r {{\sigma _i}\left( {{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right)} \\ & \Leftrightarrow \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} {\text{ }}{\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right\|_*} - \mathop {\max }\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} {\mathrm{Tr}}\left({\boldsymbol{A}}{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}{{\boldsymbol{B}}^{\mathrm{H}}}\right) \\ &\Leftrightarrow \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} {\text{ }}{\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\mathrm{H}}}} \right\|_*} - {\mathrm{Tr}}\Bigr({\boldsymbol{AW}}{{\boldsymbol{B}}^{\mathrm{H}}}\Bigr) \end{split} (35)

    至此,我们优化模型可表示为

    \begin{split} & \mathop {\min}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} {\text{ }}{\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right\|_*} - {\mathrm{Tr}}\Bigr({\boldsymbol{AW}}{{\boldsymbol{B}}^{\mathrm{H}}}\Bigr) \\ & {\mathrm{s.t}}. {{\hat{\boldsymbol{\varDelta}} }_{\rm H}} = {\boldsymbol{W}},{P_\varOmega }\left( {\boldsymbol{W}} \right) = {P_\varOmega }\left( {{{\boldsymbol{Y}}_{\rm H}}} \right) \end{split} (36)

    其中,矩阵的核范数与 {\mathrm{Tr}}({\boldsymbol{AW}}{{\boldsymbol{B}}^{\rm H}}) 均遵循凸性质。

    对于上述凸优化问题,我们采用ADMM算法进行求解[36,37]。我们先将其转化为无约束优化问题,对应的增广拉格朗日函数可以写为

    \begin{split} L \left({{\hat{\boldsymbol{\varDelta}} }_{\rm H}},{\boldsymbol{W}},{\boldsymbol{Z}}\right) =\,& {\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right\|_*} - {\mathrm{Tr}}\left({\boldsymbol{AW}}{{\boldsymbol{B}}^{\rm H}}\right) \\ &+ \frac{\beta }{2}\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}} - {\boldsymbol{W}}} \right\|_{\mathrm{F}}^2 \\ & + {\mathrm{Tr}}\left({{\boldsymbol{Z}}^{\rm H}}\left({{\hat{\boldsymbol{\varDelta}} }_{\rm H}} - {\boldsymbol{W}}\right)\right) \end{split} (37)

    其中, \beta 是惩罚因子,Z为拉格朗日乘子。

    初始化 {( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} )_1} = {{\boldsymbol{Y}}_{\rm H}} , {{\boldsymbol{W}}_1} = {( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} )_1}

    随后,关于 {\hat{\boldsymbol{\varDelta}} _{\rm H}} ,第k+1次迭代的结果可根据式(38)优化得到:

    {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} = \mathop {{\mathrm{argmin}}}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} L \left({\hat{\boldsymbol{\varDelta}} _{\rm H}},{{\boldsymbol{W}}_k},{{\boldsymbol{Z}}_k}\right) (38)

    进一步,式(38)可以表示为

    \begin{split} & {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} \\ & = \mathop {{\mathrm{argmin}}}\limits_{{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \left\{ {{{\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right\|}_*} + \frac{\beta }{2}\left\| {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}} - \left({{\boldsymbol{W}}_k} - \frac{1}{\beta }{{\boldsymbol{Z}}_k}\right)} \right\|_{\mathrm{F}}^2} \right\} \end{split} (39)

    式(39)的解可表示为

    {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} = {D_{\textstyle\frac{1}{\beta }}}\left({{\boldsymbol{W}}_k} - \frac{1}{\beta }{{\boldsymbol{Z}}_k}\right) (40)

    其中, {D_\tau }\left( \cdot \right) 为奇异值收缩算子,可以表示为

    {D_\tau }\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right) = {\boldsymbol{U}}{\mathrm{diag}}\left( {\max\left\{ {{\sigma _i}\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right) - \tau ,0} \right\}} \right){{\boldsymbol{V}}^{\rm H}} (41)

    其中, \tau 表示奇异值收缩算子门限。

    然后,第k+1次迭代的中间变量 {\boldsymbol{W}}_{{k}\text{+1}} 可根据式(42)得到:

    {{\boldsymbol{W}}_{k + 1}} = \mathop {{\mathrm{argmin}}}\limits_{{P_\varOmega }\left( {\boldsymbol{W}} \right) = {{\boldsymbol{Y}}_{\rm H}}} L \left({\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}},{\boldsymbol{W}},{{\boldsymbol{Z}}_k}\right) (42)

    式(42)关于中间变量W的子问题解为

    {{\boldsymbol{W}}_{k + 1}} = {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} + \frac{1}{\beta }\Bigr({{\boldsymbol{A}}^{\rm H}}{\boldsymbol{B}} + {{\boldsymbol{Z}}_k}\Bigr) (43)

    根据 {P_\varOmega }({\boldsymbol{W}}) = {{\boldsymbol{Y}}_{\text{H}}} 约束项,我们对中间变量 {\boldsymbol{W}}_{\text{k}\text{+1}} 更新:

    {{\boldsymbol{W}}_{k + 1}} = {P_{\varOmega \notin {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{W}}_{k + 1}}} \right) + {P_{\varOmega \in {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{Y}}_{\rm H}}} \right) (44)

    最后,根据上述更新结果,我们可以得到第k+1次迭代中关于拉格朗日乘子Z的子问题解为

    {{\boldsymbol{Z}}_{k + 1}} = {{\boldsymbol{Z}}_k} + \beta \left( {{{\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)}_{k + 1}} - {{\boldsymbol{W}}_{k + 1}}} \right) (45)

    至此,按照上述流程通过反复迭代即可恢复 {\hat{\boldsymbol{\varDelta}} _{\rm H}} 中缺失数据,然后对 {\hat{\boldsymbol{\varDelta}} _{\rm H}} 进行Hankel解码可以获得完备的时频满采样矩阵 \hat{\boldsymbol{\varDelta}} 。最后,本文利用二维快速傅里叶变换(Fast Fourier Transform, FFT)对目标进行匹配滤波以实现高分辨距离-多普勒观测谱重建。二维快速傅里叶变换的过程可以表示如下:

    {\boldsymbol{S}}{ = }{\boldsymbol{F}_{{{\mathrm{range}}}}}\hat{\boldsymbol{\varDelta}} {\boldsymbol{F}_{{{\mathrm{doppler}}}}} (46)

    其中, {\boldsymbol{F}_{\rm{range}}} {\boldsymbol{F}_{\rm{doppler}}} 表示傅里叶正交基矩阵,S表示重建后的的距离-多普勒观测谱。

    综上,基于ADMM的Hankel变换矩阵填充算法流程如算法1所示。根据文献[38]分析,整个算法处理过程中奇异值分解操作对算法运算复杂度影响最大。其中一个NM列的矩阵奇异值分解计算复杂度为

    1  基于ADMM的Hankel变换矩阵填充算法
    1.  Matrix completion algorithm combined with Hankel transformation based on ADMM
     输入:时频欠采样矩阵Y,载频编码{{U}_{n}}
     输出:时频满采样恢复矩阵 \hat{\boldsymbol{\varDelta}}
     1. 对时频欠采样矩阵Hankel变换得到 {{\boldsymbol{Y}}_{\rm H}}
     2. ADMM矩阵填充:
      (a) 初始化: {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_1} = {{\boldsymbol{W}}_1} = {{\boldsymbol{Y}}_{\rm H}} , Z = 0,迭代次数K,惩罚
       因子 \beta ,步长因子 \rho
      (b) for k = 1 to K do
        i. 更新 {\hat{\boldsymbol{\varDelta}} _{\rm H}} {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} = {D_{\textstyle\frac{1}{\beta }}}\left({{\boldsymbol{W}}_k} - \dfrac{1}{\beta }{{\boldsymbol{Z}}_k}\right)
        ii. 更新W {{\boldsymbol{W}}_{k + 1}} = {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} + \dfrac{1}{\beta }({{\boldsymbol{A}}^{\rm H}}{\boldsymbol{B}} + {{\boldsymbol{Z}}_k})
       保持原有观测采样数据不变,即:
        {{\boldsymbol{W}}_{k + 1}} = {P_{\varOmega \notin {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{W}}_{k + 1}}} \right) + {P_{\varOmega \in {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{Y}}_{\rm H}}} \right)
        iii. 更新Z {{\boldsymbol{Z}}_{k + 1}} = {{\boldsymbol{Z}}_k} + \beta \left( {{{\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)}_{k + 1}} - {{\boldsymbol{W}}_{k + 1}}} \right)
        iv. 更新 \beta {\beta _{k + 1}} = \rho {\beta _k}
      (c) end
     3. 对 {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_K} 进行Hankel逆变换得到 \hat{\boldsymbol{\varDelta}}
    下载: 导出CSV 
    | 显示表格
    O\left( {NM \cdot \min \left( {N,M} \right)} \right) (47)

    自然地,经双重Hankel矩阵构造变换后,算法计算复杂度为

    O\left( \begin{gathered} {N_1}{M_2} \\ \cdot \left( {N - {N_1} + 1} \right)\left( {M - {M_2} + 1} \right) \\ \cdot \min \left( {{N_1}{M_2},\left( {N - {N_1} + 1} \right)\left( {M - {M_2} + 1} \right)} \right) \\ \end{gathered} \right) (48)

    通常来说N1M2设置为矩阵维度的一半,即N1=N/2, M2=M/2,此参数已在文献[39]中被证明是利用双重Hankel变换提升矩阵填充性能的维度最佳选择。因此,式(48)可简化表示为 O\left( {{N^3}{M^3}} \right)

    与经典匹配滤波方法的理论计算复杂度 O\left( NM \cdot \log \left( {NM} \right) \right) 相比,为了重构未观测数据样本解决欠采样问题,矩阵填充方法付出了较高的计算负担。因此在未来我们还将研究在保持所提方法性能优势的同时进一步降低计算复杂度的改进算法。

    本节设置了单/多点目标仿真试验以验证所提方法的有效性和优越性,波形仿真参数如表1所示。此外,本文仿真中设置ADMM算法的最大迭代次数K为300,步长因子ρ为1.1,初始化惩罚因子 \beta 为0.05 ,截断奇异值数r等于目标个数。

    表  1  仿真参数
    Table  1.  Simulation parameters
    参数 数值
    中心频率{{f}_0} 93 GHz
    最小步进带宽\Delta f 18.75 MHz
    脉冲重复时间T 25 μs
    脉冲个数N 64
    捷变频点数M 64
    载频编码 {{U}}_{{n}} 随机均匀分布
    二维FFT处理时窗函数 Hamming窗
    下载: 导出CSV 
    | 显示表格

    需要说明的是,本文研究的是如何实现对单个距离粗分辨单元内目标的高分辨距离-多普勒谱精确重建,因此试验中高分辨距离含义是以所在距离粗分辨单元对应距离为基准的相对值,不代表目标相对雷达的真实距离;速度也是在补偿了雷达自身载体速度后,目标相对雷达运动产生的径向速度。

    本节设置单目标信噪比为10 dB (快时间脉冲压缩后),高分辨距离和速度为0 m, 0 m/s。采用匹配滤波、所提方法估计得到的随机步进频雷达高分辨距离-多普勒谱如图3(b)图3(c)所示;此外,时频等效满采样数据(信噪比同样也设置为10 dB)估计的高分辨距离-多普勒谱见图3(a)。可以看出,尽管匹配滤波处理可以有效地聚焦并分辨目标,但谱中起伏高旁瓣非常明显;相比之下,所提方法的估计谱与等效满采样数据估计结果非常接近,恢复效果一致,表明该方法可有效解决欠定估计问题。

    图  3  单点目标高分辨距离-速度谱估计结果
    Figure  3.  Range-velocity spectrum for single point target

    进一步地,图4展示了不同方法在高分辨距离、速度以及对角线等剖面的点目标响应函数(Point Spread Functions, PSF)。可以看出,尽管所有方法估计的主瓣宽度一样,但与等效满采样数据恢复结果的PSF相比,匹配滤波处理的PSF旁瓣在距离、速度等维度仍增加了18 dB以上,而在对角线方向则增加了28 dB。与之相反,本文提出矩阵填充处理方法可以有效抑制旁瓣,特别是在对角线方向上,其旁瓣水平甚至比等效满采样恢复结果要低,这表明矩阵填充处理可以准确地恢复未观测的数据样本,且未引入额外噪声,间接提升了数据的输入信噪比。

    图  4  点目标响应函数
    Figure  4.  Point spread functions response

    为了量化评估提出方法性能,本文采用峰值旁瓣比(Peak Side Lobe Ratio, PSLR)和积分旁瓣比(Integral Side Lobe Ratio, ISLR)两个指标来评估分析提出方法的旁瓣抑制水平,定义如下。图5图6展示了在不同输入信噪比下不同方法的PSLR和ISLR结果。其中每个结果是由300次蒙特卡洛仿真试验得到。

    图  5  PSLR随输入信噪比变化关系
    Figure  5.  PSLR versus input SNR
    图  6  ISLR随输入信噪比变化关系
    Figure  6.  ISLR versus input SNR
    \text{PLSR}=10\mathrm{lg}\frac{最大峰值旁瓣功率}{主瓣功率} (49)
    \text{ILSR}=10\mathrm{lg}\frac{平均峰值旁瓣功率}{主瓣功率}\;\; (50)

    图5显示了PSLR与输入信噪比(Signal-to-Noise Ratio, SNR)的关系。可以看出,所提方法无论在何种输入SNR下,均优于匹配滤波方法,特别是在输入SNR≥–5 dB后,性能优势更加明显;此外,当输入SNR>0 dB后,提出方法的PSLR缓慢下降,几乎达到与等效满采样数据恢复结果一样的性能。需要注意的是,当输入SNR>10 dB后,该方法的PSLR不再随输入SNR的增加而变化,表明该方法的峰值旁瓣抑制性能已达到最优。

    图6显示了ISLR与输入SNR的关系。结果表明,相比匹配滤波方法,提出方法获得了更高的平均旁瓣抑制性能。在输入SNR>−2 dB时,由于提出方法恢复的缺失采样数据不包含噪声,其数据恢复后的平均输入SNR高于等效满采样数据的平均输入SNR,因此等效满采样数据经匹配滤波处理后其平均旁瓣功率(加窗处理后实际为平均噪声功率)高于本文方法,但其随输入信噪比的增加线性降低。此外,与PSLR结果类似,当−2 dB<输入SNR<2 dB,所提方法ISLR急速下降,表明矩阵填充方法有效的前提是相对较高的输入SNR。

    综上,在本文使用波形参数和处理算法参数边界下,输入信噪比0 dB是提出方法可有效抑制旁瓣的关键边界。一般情况下,随机步进频雷达都会预先进行脉内脉冲压缩处理,在快时间维聚焦以提高目标信噪比并进行距离粗分辨,因此在大多数情况下进行高分辨距离-多普勒谱估计时目标输入信噪比基本能满足0 dB要求。

    在4.2节中,我们开展了多目标场景下不同方法估计性能的对比,仿真中共设置了4个理想点目标,其中强目标3个,输入信噪比为10 dB,以高分辨距离0 m、速度0 m/s为中心等间隔分布;弱目标1个,输入信噪比为2 dB,与强目标相对距离1.8 m、速度18 m/s。图7(a)为等效时频满采样回波数据处理结果;图7(b)图7(c)则展示了通过匹配滤波和本文所提方法恢复数据后估计得到的高分辨距离-多普勒谱。结果表明,传统匹配滤波方法估计得到的频谱中存在严重的旁瓣基底,这直接导致了设置的小信噪比目标被大信噪比目标的随机起伏旁瓣所遮盖,无法有效检测;同时,在估计谱中出现若干个虚警目标。与之相比,经过矩阵填充处理后的估计谱旁瓣可以得到有效抑制。

    图  7  多目标高分辨距离-速度谱估计结果
    Figure  7.  Range-velocity spectrum for multi-target

    图8显示了多目标在距离、速度剖面的PSF性能。可以看出,提出方法的聚焦效果良好,可以准确分辨多目标,且主瓣宽度与等效满采样恢复谱几乎一致,表明提出方法没有损失高分辨距离-多普勒的分辨率。尽管提出方法估计谱中仍存在一些随机波动的起伏旁瓣,但其旁瓣平均值也远低于匹配滤波结果,证明了所提方法可以在多目标环境下有效抑制旁瓣。需要注意的是,与时频等效满采样数据情况相比,提出方法的旁瓣抑制水平在目标个数增多时有所下降,这是因为随着目标数量的增加,待恢复矩阵秩也会随之增大,势必会影响算法恢复性能。此外,图8(c)还比较了多目标场景中对角线方向的PSF,结果表明传统匹配滤波方法恢复的弱目标主瓣(归一化后功率–8 dB左右)被淹没在强目标的高功率旁瓣(归一化后平均功率–13 dB)中,无法被正确地检测。相比之下,提出方法可以有效地抑制旁瓣,聚焦弱目标,取得较好的谱估计效果。

    图  8  多目标响应函数
    Figure  8.  Point spread functions response for multi-target

    最后,本文计算了在输入信噪比10 dB情况下,不同方法的PSLR, ISLR与目标数量的关系,每次结果均进行了300次蒙特卡罗仿真试验,结果如图9图10所示。可以看出所提方法的旁瓣抑制性能随着目标数量K的增大而下降。这是因为在多目标场景下,双重Hankel重构虽能一定程度上保持矩阵的低秩特性,但无法完全解决因为目标数增加导致矩阵秩快速变大,因此旁瓣抑制水平快速恶化。当然正如前文所说,快时间脉冲压缩预处理可以有效降低距离粗分辨单元的目标数,一定程度上使得待恢复矩阵满足低秩特性,进而保证提出方法的谱估计性能。

    图  9  多目标场景PSLR随目标数量变化关系
    Figure  9.  PSLR versus the number of targets
    图  10  多目标场景ISLR随目标数量变化关系
    Figure  10.  ISLR versus the number of targets

    进一步,本文还比较了所提方法与压缩感知方法估计高分辨距离-多普勒谱的性能。其中压缩感知方法优化模型如下,优化算法采用CVX工具包[40]

    \begin{split} & \mathop {\min}\limits_{\boldsymbol{s}} {\left\| {\boldsymbol{s}} \right\|_1} \\ & {\mathrm{s.t.}} {\left\| {{\boldsymbol{y}} - {\boldsymbol{Cs}}} \right\|_2} \le \varepsilon \end{split} (51)

    其中,s为矢量化后的待恢复高分辨距离-速度谱,y为时频欠采样信号,其第n个元素由式(7)表示,C是根据{{U}_{n}}构造的感知矩阵, \varepsilon 为算法迭代中止门限(与本文方法门限保持一致,均为1E–6)。

    仿真中考虑两个点目标场景,并且目标参数均没有位于预设的栅格点之上。所提方法与CS方法估计的距离速度谱如图11所示。结果表明,尽管压缩感知方法可以有效抑制旁瓣,然而谱中也出现了很多虚假目标。比较两者在对角线的PSF (见图12),可以发现提出方法的聚焦效果更好;与之相反,压缩感知方法能量有所损失。

    图  11  高分辨距离-速度谱估计结果
    Figure  11.  Range-velocity spectrum for multi-target
    图  12  距离-速度对角线点目标响应函数
    Figure  12.  Point spread function response in diagonal line

    表2显示了在配备AMD Ryzen9 4900HS, 3 GHz处理器和16 GB内存的计算机上测试所提方法和CS方法的MATLAB运行时间。可以发现,CS方法的计算时间随着网格数量的增加而显著提升,而本文所提方法的运行时间(约5 s)总是远小于CS方法,这表明提出的方法比CS方法需要更低的计算成本。

    表  2  运行时间对比结果(s)
    Table  2.  Runtimes comparison (s)
    栅格点数 所提方法 压缩感知 栅格点数 所提方法 压缩感知
    64×64 5.6 23.6 256×256 5.7 412.6
    128×128 5.6 90.5 320×320 5.7 810.9
    192×192 5.7 187.1
    下载: 导出CSV 
    | 显示表格

    此外,在处理高维数据时,对于CS方法,每个维度必须被重构为向量形式,处理后再形成高维矩阵,这导致相当大的内存使用,而本文方法将矩阵数据作为一个整体来处理,在减少内存使用方面具有巨大潜力。

    为了进一步验证方法的有效性,对实测数据进行高分辨距离-速度估计,目标是一个运动的厢货卡车(如图13所示),回波信号通过一部93 GHz毫米波雷达样机采样获得,其中雷达参数如表3所示。

    图  13  厢货卡车实测数据采集场景
    Figure  13.  Radar prototype and test scenario
    表  3  雷达波形参数
    Table  3.  Waveform parameters of the test radar
    参数 数值
    中心频率{{f}_0} 93.31 GHz
    最小步进带宽\Delta f 18.75 MHz
    脉冲重复时间T 25.8 μs
    捷变频点数M 64
    脉冲个数N 64
    调频带宽B 25 MHz
    脉冲宽度Tp 2 μs
    载频编码 {{U}}_{{n}} 随机均匀分布
    下载: 导出CSV 
    | 显示表格

    试验中首先通过Matlab生成遵循均匀分布的随机跳频序列,并将其存储到雷达样机中。回波信号经预处理(下混频、滤波、采样、脉冲压缩等)后存储到计算机中。采用匹配滤波、压缩感知以及本文所提方法估计的高分辨距离-多普勒谱如图14所示。

    图  14  距离-速度谱估计结果(红框内为目标强散射点)
    Figure  14.  Range-velocity spectrum of MF, CS and MC (the target strong scattering points are in the red box)

    由于欠定估计,匹配滤波的谱中产生了较高的起伏旁瓣,导致恢复结果变差。提出方法的结果远优于匹配滤波方法,初步验证了矩阵填充方法在实际系统中的有效性。压缩感知方法虽能抑制旁瓣,但其依赖稀疏基的设计,需要目标的先验信息,对未知场景与目标成像效果差。

    我们提取了v=–8.75 m/s的高分辨率距离像,以进一步分析CS和所提方法的估计性能,如图15所示。

    图  15  目标高分辨距离像
    Figure  15.  High resolution range profiles of target

    显然,所提出的方法拥有相对更为优越的成像结果,具体表现在可以很好地在高分辨率距离像上呈现出目标的大体轮廓,从而提供更多关于该目标的特征信息;此外,所提出的方法结果与匹配滤波恢复结果的主瓣能量、散射点个数等匹配更好。相反,压缩感知恢复结果在目标分段区域没有被完全聚焦,估计的HRRP质量较差。因此,我们可以得出结论,与CS方法相比,本文方法可以获得更好的HRRP聚焦效果,对于目标参数可以实现高精度估计,更有利于后续的目标分类识别。

    本文提出了一种基于矩阵填充的随机步进频雷达高分辨率距离-多普勒谱稀疏恢复方法。该方法利用矩阵填充技术恢复未观测的时频数据,解决了随机步进频波形固有的二维欠采样问题。与匹配滤波方法相比,该方法在保持能量、分辨率、模糊度等波形性能不变情况下,能够有效抑制抬高的旁瓣。此外,与压缩感知方法相比,所提方法不仅在非栅格点谱恢复方面取得了更好的性能,并且显著降低了计算负担。仿真和实测数据验证了该方法的有效性和优越性。

  • 图  1  调频RSF信号预处理流程图

    Figure  1.  Diagram of the RSF chirp signal pre-processing flow

    图  2  回波信号矩阵填充处理流程图

    Figure  2.  Matrix completion flow

    图  3  单点目标高分辨距离-速度谱估计结果

    Figure  3.  Range-velocity spectrum for single point target

    图  4  点目标响应函数

    Figure  4.  Point spread functions response

    图  5  PSLR随输入信噪比变化关系

    Figure  5.  PSLR versus input SNR

    图  6  ISLR随输入信噪比变化关系

    Figure  6.  ISLR versus input SNR

    图  7  多目标高分辨距离-速度谱估计结果

    Figure  7.  Range-velocity spectrum for multi-target

    图  8  多目标响应函数

    Figure  8.  Point spread functions response for multi-target

    图  9  多目标场景PSLR随目标数量变化关系

    Figure  9.  PSLR versus the number of targets

    图  10  多目标场景ISLR随目标数量变化关系

    Figure  10.  ISLR versus the number of targets

    图  11  高分辨距离-速度谱估计结果

    Figure  11.  Range-velocity spectrum for multi-target

    图  12  距离-速度对角线点目标响应函数

    Figure  12.  Point spread function response in diagonal line

    图  13  厢货卡车实测数据采集场景

    Figure  13.  Radar prototype and test scenario

    图  14  距离-速度谱估计结果(红框内为目标强散射点)

    Figure  14.  Range-velocity spectrum of MF, CS and MC (the target strong scattering points are in the red box)

    图  15  目标高分辨距离像

    Figure  15.  High resolution range profiles of target

    1  基于ADMM的Hankel变换矩阵填充算法

    1.   Matrix completion algorithm combined with Hankel transformation based on ADMM

     输入:时频欠采样矩阵Y,载频编码{{U}_{n}}
     输出:时频满采样恢复矩阵 \hat{\boldsymbol{\varDelta}}
     1. 对时频欠采样矩阵Hankel变换得到 {{\boldsymbol{Y}}_{\rm H}}
     2. ADMM矩阵填充:
      (a) 初始化: {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_1} = {{\boldsymbol{W}}_1} = {{\boldsymbol{Y}}_{\rm H}} , Z = 0,迭代次数K,惩罚
       因子 \beta ,步长因子 \rho
      (b) for k = 1 to K do
        i. 更新 {\hat{\boldsymbol{\varDelta}} _{\rm H}} {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} = {D_{\textstyle\frac{1}{\beta }}}\left({{\boldsymbol{W}}_k} - \dfrac{1}{\beta }{{\boldsymbol{Z}}_k}\right)
        ii. 更新W {{\boldsymbol{W}}_{k + 1}} = {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_{k + 1}} + \dfrac{1}{\beta }({{\boldsymbol{A}}^{\rm H}}{\boldsymbol{B}} + {{\boldsymbol{Z}}_k})
       保持原有观测采样数据不变,即:
        {{\boldsymbol{W}}_{k + 1}} = {P_{\varOmega \notin {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{W}}_{k + 1}}} \right) + {P_{\varOmega \in {\varOmega _{\rm H}}}}\left( {{{\boldsymbol{Y}}_{\rm H}}} \right)
        iii. 更新Z {{\boldsymbol{Z}}_{k + 1}} = {{\boldsymbol{Z}}_k} + \beta \left( {{{\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)}_{k + 1}} - {{\boldsymbol{W}}_{k + 1}}} \right)
        iv. 更新 \beta {\beta _{k + 1}} = \rho {\beta _k}
      (c) end
     3. 对 {\left( {{{\hat{\boldsymbol{\varDelta}} }_{\rm H}}} \right)_K} 进行Hankel逆变换得到 \hat{\boldsymbol{\varDelta}}
    下载: 导出CSV

    表  1  仿真参数

    Table  1.   Simulation parameters

    参数 数值
    中心频率{{f}_0} 93 GHz
    最小步进带宽\Delta f 18.75 MHz
    脉冲重复时间T 25 μs
    脉冲个数N 64
    捷变频点数M 64
    载频编码 {{U}}_{{n}} 随机均匀分布
    二维FFT处理时窗函数 Hamming窗
    下载: 导出CSV

    表  2  运行时间对比结果(s)

    Table  2.   Runtimes comparison (s)

    栅格点数 所提方法 压缩感知 栅格点数 所提方法 压缩感知
    64×64 5.6 23.6 256×256 5.7 412.6
    128×128 5.6 90.5 320×320 5.7 810.9
    192×192 5.7 187.1
    下载: 导出CSV

    表  3  雷达波形参数

    Table  3.   Waveform parameters of the test radar

    参数 数值
    中心频率{{f}_0} 93.31 GHz
    最小步进带宽\Delta f 18.75 MHz
    脉冲重复时间T 25.8 μs
    捷变频点数M 64
    脉冲个数N 64
    调频带宽B 25 MHz
    脉冲宽度Tp 2 μs
    载频编码 {{U}}_{{n}} 随机均匀分布
    下载: 导出CSV
  • [1] PANDA S S S, PANIGRAHI T, PARNE S R, et al. Recent advances and future directions of microwave photonic radars: A review[J]. IEEE Sensors Journal, 2021, 21(19): 21144–21158. doi: 10.1109/JSEN.2021.3099533.
    [2] 向寅, 张凯, 胡程. 基于NUFFT的调频步进频高分辨成像与目标识别算法[J]. 雷达学报, 2015, 4(6): 639–647. doi: 10.12000/JR15083.

    XIANG Yin, ZHANG Kai, and HU Cheng. A NUFFT based step-frequency Chirp signal high resolution imaging algorithm and target recognition algorithm[J]. Journal of Radars, 2015, 4(6): 639–647. doi: 10.12000/JR15083.
    [3] AL-HOURANI A, EVANS R J, MORAN B, et al. Efficient range-Doppler processing for random stepped frequency radar in automotive applications[C]. IEEE 85th Vehicular Technology Conference, Sydney, Australia, 2017: 1–7. doi: 10.1109/VTCSpring.2017.8108414.
    [4] SAPONARA S, GRECO M S, and GINI F. Radar-on-chip/in-package in autonomous driving vehicles and intelligent transport systems: Opportunities and challenges[J]. IEEE Signal Processing Magazine, 2019, 36(5): 71–84. doi: 10.1109/MSP.2019.2909074.
    [5] JIANG Yuan, WANG Yanhua, LI Yang, et al. Eigenvalue-based ground target detection in high-resolution range profiles[J]. IET Radar, Sonar & Navigation, 2020, 14(11): 1747–1756. doi: 10.1049/iet-rsn.2020.0002.
    [6] AXELSSON S R J. Analysis of random step frequency radar and comparison with experiments[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(4): 890–904. doi: 10.1109/TGRS.2006.888865.
    [7] HUANG Tianyao, LIU Yimin, LI Gang, et al. Randomized stepped frequency ISAR imaging[C]. IEEE International Radar Conference, Atlanta, USA, 2012: 553–557. doi: 10.1109/RADAR.2012.6212202.
    [8] WEHNER D R. High Resolution Radar[M]. Norwood, MA, Artech House, 1987.
    [9] LIU Yimin, MENG Huadong, LI Gang, et al. Range-velocity estimation of multiple targets in randomised stepped-frequency radar[J]. Electronics Letters, 2008, 44(17): 1032–1034. doi: 10.1049/el:20081608.
    [10] ACKROYD M H and GHANI F. Optimum mismatched filters for sidelobe suppression[J]. IEEE Transactions on Aerospace and Electronic Systems, 1973, AES-9(2): 214–218. doi: 10.1109/TAES.1973.309769.
    [11] DAVIS R M, FANTE R L, and PERRY R P. Phase-coded waveforms for radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(1): 401–408. doi: 10.1109/TAES.2007.357142.
    [12] KAJENSKI P J. Mismatch filter design via convex optimization[J]. IEEE Transactions on Aerospace and Electronic Systems, 2016, 52(4): 1587–1591. doi: 10.1109/TAES.2016.140556.
    [13] LIU Shuai, CAO Yunhe, YEO T S, et al. Range sidelobe suppression for randomized stepped-frequency chirp radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2021, 57(6): 3874–3885. doi: 10.1109/TAES.2021.3082670.
    [14] LIAO Zhikun, HU Jiemin, LU Dawei, et al. Motion analysis and compensation method for random stepped frequency radar using the pseudorandom code[J]. IEEE Access, 2018, 6: 57643–57654. doi: 10.1109/ACCESS.2018.2873784.
    [15] LIAO Zhikun, LU Dawei, HU Jiemin, et al. Waveform design for random stepped frequency radar to estimate object velocity[J]. Electronics Letters, 2018, 54(14): 894–896. doi: 10.1049/el.2018.1033.
    [16] QUAN Yinghui, LI Yachao, HU Wen, et al. FM sequence optimisation of chaotic-based random stepped frequency signal in through-the-wall radar[J]. IET Signal Processing, 2017, 11(7): 830–837. doi: 10.1049/iet-spr.2015.0565.
    [17] HUANG Tianyao, LIU Yimin, XU Xingyu, et al. Analysis of frequency agile radar via compressed sensing[J]. IEEE Transactions on Signal Processing, 2018, 66(23): 6228–6240. doi: 10.1109/TSP.2018.2876301.
    [18] YOON Y S, HONG Yunseog, and KIM S. Simple strategies to build random compressive sensing matrices in step-frequency radars[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(9): 1357–1361. doi: 10.1109/LGRS.2018.2841189.
    [19] HUANG Tianyao, LIU Yimin, MENG Huadong, et al. Cognitive random stepped frequency radar with sparse recovery[J]. IEEE Transactions on Aerospace and Electronic Systems, 2014, 50(2): 858–870. doi: 10.1109/TAES.2013.120443.
    [20] QUAN Yinghui, LI Yachao, WU Yaojun, et al. Moving target detection for frequency agility radar by sparse reconstruction[J]. Review of Scientific Instruments, 2016, 87(9): 094703. doi: 10.1063/1.4962700.
    [21] CHI Yuejie, SCHARF L L, PEZESHKI A, et al. Sensitivity to basis mismatch in compressed sensing[J]. IEEE Transactions on Signal Processing, 2011, 59(5): 2182–2195. doi: 10.1109/TSP.2011.2112650.
    [22] WANG Mou, WEI Shunjun, SHI Jun, et al. CSR-Net: A novel complex-valued network for fast and precise 3-D microwave sparse reconstruction[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 4476–4492. doi: 10.1109/JSTARS.2020.3014696.
    [23] CANDES E J and TAO T. The power of convex relaxation: Near-optimal matrix completion[J]. IEEE Transactions on Information Theory, 2010, 56(5): 2053–2080. doi: 10.1109/TIT.2010.2044061.
    [24] 马宇欣, 海宇, 李中余, 等. 稀疏轨迹毫米波雷达三维高分辨成像算法[J]. 雷达学报, 2023, 12(5): 1000–1013. doi: 10.12000/JR23001.

    MA Yuxin, HAI Yu, LI Zhongyu, et al. 3D high-resolution imaging algorithm with sparse trajectory for millimeter-wave radar[J]. Journal of Radars, 2023, 12(5): 1000–1013. doi: 10.12000/JR23001.
    [25] SUN Shunqiao and ZHANG Y D. 4D automotive radar sensing for autonomous vehicles: A sparsity-oriented approach[J]. IEEE Journal of Selected Topics in Signal Processing, 2021, 15(4): 879–891. doi: 10.1109/JSTSP.2021.3079626.
    [26] SUN Shunqiao, BAJWA W U, and PETROPULU A P. MIMO-MC radar: A MIMO radar approach based on matrix completion[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(3): 1839–1852. doi: 10.1109/TAES.2015.140452.
    [27] HU Xiaowei, TONG Ningning, WANG Jianye, et al. Matrix completion-based MIMO radar imaging with sparse planar array[J]. Signal Processing, 2017, 131: 49–57. doi: 10.1016/j.sigpro.2016.07.034.
    [28] ZHANG Yilong, LI Yuehua, CHEN Jianfei, et al. Sparse millimeter-wave InSAR imaging approach based on MC[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(5): 714–718. doi: 10.1109/LGRS.2018.2810234.
    [29] HU Xueyao, LI Yang, LU Man, et al. A multi-carrier-frequency random-transmission Chirp sequence for TDM MIMO automotive radar[J]. IEEE Transactions on Vehicular Technology, 2019, 68(4): 3672–3685. doi: 10.1109/TVT.2019.2900357.
    [30] LIANG Can, LI Yang, HU Xueyao, et al. Coherent-on-Receive synthesis using dominant scatterer in millimeter-wave distributed coherent aperture radar[J]. Remote Sensing, 2023, 15(6): 1505. doi: 10.3390/rs15061505.
    [31] STRANG G. Introduction to Linear Algebra[M]. 6th ed. Wellesley, USA: Wellesley-Cambridge Press, 2022.
    [32] HORN R A and JOHNSON C R. Matrix Analysis[M]. Cambridge, UK: Cambridge University Press, 1985. doi: 10.1017/CBO9780511810817.
    [33] YE Hailiang, LI Hong, CAO Feilong, et al. A hybrid truncated norm regularization method for matrix completion[J]. IEEE Transactions on Image Processing, 2019, 28(10): 5171–5186. doi: 10.1109/TIP.2019.2918733.
    [34] HU Yao, ZHANG Debing, YE Jieping, et al. Fast and accurate matrix completion via truncated nuclear norm regularization[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2013, 35(9): 2117–2130. doi: 10.1109/TPAMI.2012.271.
    [35] CAI Jianfeng, CANDÈS E J, and SHEN Zuowei. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization, 2010, 20(4): 1956–1982. doi: 10.1137/080738970.
    [36] 邓理康, 张双辉, 张弛, 等. 一种基于多维交替方向乘子法的多输入多输出逆合成孔径雷达成像方法[J]. 雷达学报, 2021, 10(3): 416–431. doi: 10.12000/JR20132.

    DENG Likang, ZHANG Shuanghui, ZHANG Chi, et al. A multiple-input multiple-output inverse synthetic aperture radar imaging method based on multidimensional alternating direction method of multipliers[J]. Journal of Radars, 2021, 10(3): 416–431. doi: 10.12000/JR20132.
    [37] 范文, 蔚保国, 陈镜, 等. 基于波形优化和天线位置选择的MIMO雷达波束扫描算法研究[J]. 雷达学报, 2022, 11(4): 530–542. doi: 10.12000/JR22135.

    FAN Wen, YU Baoguo, CHEN Jing, et al. Joint waveform optimization and antenna position selection for MIMO radar beam scanning[J]. Journal of Radars, 2022, 11(4): 530–542. doi: 10.12000/JR22135.
    [38] OH T H, MATSUSHITA Y, TAI Y W, et al. Fast randomized singular value thresholding for low-rank optimization[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2018, 40(2): 376–391. doi: 10.1109/TPAMI.2017.2677440.
    [39] CHEN Yuxin and CHI Yuejie. Robust spectral compressed sensing via structured matrix completion[J]. IEEE Transactions on Information Theory, 2014, 60(10): 6576–6601. doi: 10.1109/TIT.2014.2343623.
    [40] GRANT M, BOYD S, and YE Yinyu. MATLAB software for disciplined convex programming[EB/OL]. https://web.stanford.edu/~boyd/papers/disc_cvx_prog.html, 2006.
  • 加载中
图(15) / 表(4)
计量
  • 文章访问数: 697
  • HTML全文浏览量: 141
  • PDF下载量: 207
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-10-03
  • 修回日期:  2023-12-28
  • 网络出版日期:  2024-01-09
  • 刊出日期:  2024-02-28

目录

/

返回文章
返回