A Track Initiation Method for FM-based Passive Radar Network Based on Multiple Elementary Hypotheses
-
摘要: 基于调频(FM)广播信号的外辐射源雷达有着检测概率低、虚警率高、量测精度差的特点,这给组网目标跟踪带来了极大挑战。一方面,较高的虚警率使计算量增加,组网算法的实时性受到考验;另一方面,检测概率低、方位角精度差造成冗余信息缺乏,量测关联与航迹起始变得困难。为解决这些问题,该文提出初级假设点和初级假设航迹的概念,以及基于此概念的FM广播外辐射源雷达网航迹起始算法。首先构造可能的低维关联假设,并解算出与其对应的初级假设点;随后关联不同时刻的初级假设点,形成多条可能的初级假设航迹;最后联合多场雷达网数据进行假设航迹判决,真实目标对应的初级假设航迹会得到确认,错误关联导致的虚假初级假设航迹会被剔除。相比于已有算法,所提算法有着更低的计算量,更快的航迹起始速度,仿真与实测结果均验证了所提算法的有效性。
-
关键词:
- FM广播外辐射源雷达 /
- 初级假设 /
- 航迹起始 /
- 雷达网信息融合
Abstract: Passive radars based on FM radio signals have low detection probability, high false alarm rates and poor accuracy, presenting considerable challenges to target tracking in radar networks. Moreover, a high false alarm rate increases the computational burden and puts forward high requirements for the real-time performance of networking algorithms. In addition, low detection probability and poor azimuth accuracy result in a lack of redundant information, making measurement association and track initiation challenging. To address these issues, this paper proposes an FM-based passive radar network based on the concepts of elementary hypothesis points and elementary hypothesis track, as well as a track initiation algorithm. First, we construct possible low-dimensional association hypotheses and solve for their corresponding elementary hypothesis points. Subsequently, we associate elementary hypothesis points from different frames to form multiple possible elementary hypothesis tracks. Finally, by combining multi-frame radar network data for hypothesis track judgment, we confirm the elementary hypothesis tracks corresponding to the real targets, and eliminate the false elementary hypothesis tracks caused by incorrect associations. Result reveal that the proposed algorithm has lower computational complexity and faster track initiation speed than existing algorithms. Moreover, we verified the effectiveness of the proposed algorithm using simulation and experimental results. -
1. 引言
近些年,频率步进(Stepped Frequency, SF)波形被广泛应用于雷达高精度探测领域[1−4]。它由一系列载频线性步进脉冲信号组成,通过相参合成这些信号获得大带宽效果,可以大幅提升雷达距离分辨率,同时避免因发射瞬时大宽带信号增加的硬件系统复杂度[5]。频率步进波形的一种改进方式是它的载频步进值随机改变,通常也被称之为随机步进频(Random Stepped Frequency, RSF)波形[6]。与常规SF波形相比,RSF波形具有两个显著的优点:首先,随机跳频的方式可以大幅降低信号被截获或“碰撞”概率,从而提升雷达抗干扰能力;其次,由于载频非线性变化,RSF波形解耦了高分辨距离和多普勒相位,产生图钉型模糊函数,从而获得了估计目标多普勒速度的能力[7]。然而,由于RSF波形对环境的感知形式在慢时间-载频(时-频)是稀疏的,因而雷达在时-频域对回波信号的采样是随机欠采样,造成相参信息缺失,传统匹配滤波(Matched Filter, MF)相参处理会演化为欠定估计,导致高分辨距离-多普勒谱中产生起伏高旁瓣[8,9]。此时,低信噪比目标可能被高信噪目标的旁瓣掩盖,导致漏检;或者起伏波动的旁瓣可能被检测为假目标,产生虚警。显然,上述问题限制了RSF波形的广泛应用。
为了解决该问题,文献[10−12]提出使用失配滤波器方法减轻高旁瓣对参数估计影响,但这类方法会降低目标信噪比或雷达距离分辨率。文献[13]提出使用RSF回波的包络和相位设计用于高分辨率距离像合成的旁瓣抑制滤波器。然而该方法只能消除单个多普勒频率下的起伏旁瓣;并且,当不能预先获得目标的多普勒信息时,其旁瓣抑制性能仍有待提高。从波形设计的角度出发,可以通过设计跳频序列来抑制旁瓣[14−16],然而这些方法并没有考虑到多目标应用场景,因此所设计的波形不能保证多目标的旁瓣在相互叠加后仍保持在低功率水平。
随着压缩感知(Compressed Sensing, CS)技术的兴起,该技术已被用于RSF波形中解决因随机欠采样造成的欠定估计问题[17−20]。其本质是利用回波信号数据中的稀疏性,通过设计特定的测量矩阵寻找信号在高分辨距离-多普勒域的稀疏表示,实现低旁瓣稀疏恢复。但CS恢复性能在很大程度上受设计的测量矩阵制约,当测量矩阵的稀疏基与场景中目标参数的分布不完全匹配时(又称“基不匹配”),恢复性能会受到显著影响[21]。文献[22]提出结合机器学习的压缩感知方法,通过网络自适应学习获得更适应场景特征的测量矩阵;然而,该方法的训练过程需要海量的样本数据,且CS方法在高维度空间中使用一维矢量运算,导致计算复杂度激增。
矩阵填充(Matrix Completion, MC)[23,24]是一种通过使用其部分有限观测数据在低秩约束下恢复矩阵中缺失数据的技术。文献[25−27]中已展示了MC在多输入多输出(Multiple Input Multiple Output, MIMO)雷达中的有效性和优势,其已被用于恢复存在随机稀疏轨迹的阵列方位或俯仰维的欠采样数据,保证多维稀疏成像质量。此外,MC技术也被证明在处理高维数据时具有效率高、恢复精度高等优势[28]。对于经过快时间脉冲压缩处理的RSF回波信号,由于其沿时-频域呈现二维随机欠采样特性,因而时-频维度中的回波数据也有所缺失,因此可以使用MC技术恢复未采样的数据。
基于上述分析,本文提出了一种基于Hankel重构矩阵填充的随机步进频雷达高分辨距离-多普勒谱估计方法,实现低旁瓣谱稀疏恢复。首先,推导并建立了随机步进频雷达单个粗分辨距离单元内回波时-频欠采样模型;并利用Hankel矩阵变换重构待恢复数据矩阵,证明了矩阵满足低秩先验特性,为基于矩阵填充的谱估计算法奠定了基础。随后,将低秩矩阵缺失数据恢复问题转化为以填充后矩阵截断核范数最小化为目标、以填充后矩阵原采样位置元素不变为约束的优化模型,进而采用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)求解模型,实现时-频数据矩阵填充,恢复缺失相参信息,解决匹配滤波欠定估计问题。所提方法直接以矩阵形式进行二维稀疏建模与求解,不仅有效避免了基不匹配问题,在保持良好重构性能的同时还进一步降低了运算量,极大地提升了运算速度。仿真和实测数据试验验证了方法的有效性和优越性。
2. 回波建模与信号特性分析
2.1 RSF回波建模及预处理
本文以调频RSF雷达为例,其发射信号由一系列脉内相同线性调频、脉间不同载频伪随机调制脉冲组成,信号预处理流程如图1所示。其中,第n个脉冲的载频为
fn=f0+UnΔf ,f0 为初始频率,Un 表示第n个脉冲(共发射N个脉冲)的载频编码,且Un∈{1,2,⋯,M} ,M为随机步进频率总点数,Δ f 为最小频率步进间隔,MΔf 为波形合成带宽。因此第n个脉冲的发射信号可以表示为st(t)=s(t)⋅exp[j2πfn(t−nT)]=rect(t−nTTp)⋅exp[j2πB2Tp(t−nT)2]⋅exp[j2πfn(t−nT)] (1) 其中,
s(t) 表示雷达基带信号,t表示脉内快时间,B 为脉内调频带宽,Tp 为发射调频信号的脉宽,T 为脉冲重复间隔。假设发射信号经理想点目标调制后幅度为A,回波信号可表示为
sr(t)=A⋅rect(t−nT−τTp)⋅exp[j2πB2Tp(t−nT−τ)2]⋅exp[j2πfn(t−nT−τ)] (2) 其中,
τ 表示回波信号的双程时延。基于“走-停”模型假设[29,30],时延τ 可以近似表示为τ≈2(r−vnT)c (3) 其中,r, v分别为目标的初始距离和相对雷达的径向速度,c为光速。
回波信号经模数转换、下变频、滤波处理后可表示为
sr(l,n)=A⋅exp{j2π[B2Tp(lts−nT−τ)2−fnτ]} (4) 其中,
l=0,1,⋯,L−1 为脉内快时间采样索引,L为一个脉冲内的快时间采样总点数,ts 代表快时间采样间隔。信号经快时间脉冲压缩处理后可表示为
y(l′,n)=sr(l,n)∗ˉs(l′−l,n)≈α⋅sinc(l′−2rBc)⋅exp(j4πUnΔfrc)⋅exp(j4πvnTf0c) (5) 其中,
∗ 表示卷积运算,ˉs(l′−l,n) 为式(1)中雷达基带信号s(t) 的共轭转置,α=A⋅exp(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)=K∑k=1αk⋅exp(j4πUnΔfrkc)⋅exp(j4πvknTf0c) (7) 其中,K表示一个CRRC内的目标数量,
αk ,rk 和vk 分别表示第k个目标的复振幅、初始距离和相对径向速度。2.2 回波信号矩阵形式变换及特性分析
考虑式(7)对目标
αk ,rk 和vk 的估计本质是对整个观测谱的重建,即估计每个高分辨距离-速度单元内的散射强度,再通过对应的栅格位置确定估计目标的幅度、距离、速度参数等。因此,完成MN个单元散射强度的重建即可获得目标参数rk 和vk 。显然,网格的格点数取决于重建观测谱的精细程度,精细程度越高,需要的网格越密,格点个数越多。考虑网格间距要小于或等于波形距离、速度分辨率,因此设置的距离和速度单元分别最少为M, N个。综上,待估计的高分辨距离和多普勒参数数量(M×N )远大于采样数据量(N),即RSF回波的两维随机欠采样问题。此时,若直接采用传统匹配滤波方法进行相参处理,参数估计会演化为欠定估计问题,导致谱中高旁瓣基底产生。针对上述问题,本文提出了一种基于矩阵填充的方法补全缺失的采样数据,避免欠定估计。首先,将式(7)中的回波信号重新排列为矩阵形式。流程如图2所示,根据载频编码顺序将采样数据重新排列至相应的时频采样位置,空白采样位置处设置为0。此时,在重新排列的欠采样矩阵Y中,第n行第m列元素可以表示为
(Y)n,m={y(n), m=Un0, m≠Un (8) 上述欠采样矩阵也可以视作是来自于完备矩阵的随机稀疏观测,即:
(Y)n,m={(Δ)n,m, (n,m)∈Ω0, (n,m)∉Ω (9) 其中,
Δ 为等效时频满采样矩阵,Ω 为载频编码模式Un 对应的时频数据集合。本文的目的即通过欠采样矩阵Y恢复满采样矩阵Δ 。基于式(7),由第k个目标构成的时频满采样矩阵中第n行第m列元素可以表示为
(Δk)n,m=αk⋅exp(j2π(nT2vkf0c+mΔf2rkc)) (10) 显然,式(10)可以改写为两个向量的乘积:
Δk=αkad(vk)ar(rk)H (11) 其中,两个向量分别为第k个目标关于速度和高分辨率距离的导向矢量,分别为
ad(vk)=[1, exp(j4πvkf0Tc), ⋯, exp(j4πvkf0(N−1)Tc)]H (12) ar(rk)=[1, exp(j4πrkΔfc), ⋯, exp(j4πrk(M−1)Δfc)]H (13) 根据矩阵乘积秩定理[31],两个矩阵相乘得到新矩阵,其秩不大于相乘前任意一个矩阵的秩,因此上述第k个目标的时频满采样矩阵秩为1。进一步地,由K个目标组成的等效时频满采样矩阵可以视为由单个目标组成的满采样矩阵的线性组合,因此可以得到:
Δ=K∑k=1Δk=K∑k=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(K∑k=1Δk)≤K∑k=1rank(Δk)=K (16) 2.3 Hankel矩阵变换及特性分析
矩阵填充方法要求待恢复矩阵满足低秩特性。为了保证该前提,本文引入了双重Hankel变换重新构造待恢复矩阵,从而增强回波信号等效满采样矩阵的低秩特性。具体地,首先定义一个维度为
N1×(N−N1+1) 的Hankel矩阵:ΔH=[Δ1Δ2⋯ΔN−N1+1Δ2Δ3⋯ΔN−N1+2⋮⋮Δn⋮ΔN1ΔN1+1⋯ΔN] (17) 其中,第n个块矩阵是一个维度为
M2×(M−M2+1) 的Hankel矩阵,该矩阵由矩阵Δ 的第n行所组成,即:Δn=[(Δ)n,1(Δ)n,2⋯(Δ)n,M−M2+1(Δ)n,2(Δ)n,3⋯(Δ)n,M−M2+2⋮⋮(Δ)n,m⋮(Δ)n,M2(Δ)n,M2+1⋯(Δ)n,M] (18) 其中,
N1 和M2 分别是构造一阶和二阶Hankel变换矩阵的维度参数。对具有双重Hankel形式的时频观测矩阵
ΔH 进行分解可以得到:ΔH=[Xr,LD1d⋮Xr,LDnd⋮Xr,LDN1d]Λ[D1dXHr,R⋯DN−n+1dXHr,R⋯DN−N1+1dXHr,R]=˜XdΛ˜Xr (19) 其中,
{Xr,L=(Xr)1:M2,:Dnd=diag{(Xd)n,:}DN−n+1d=diag{(Xd)N−n+1,:}Xr,R=(Xr)1:M−M2+1,: (20) 显然,由式(20)可知:
rank(ΔH)=rank(Δ)≤K (21) 而经Hankel构造变换后,矩阵的维度可达
N1M2×(N−N1)(M−M2+1) ,因此变换后的双重Hankel矩阵具有更优异的低秩特性。3. 基于ADMM的数据恢复算法
根据第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=(N−N1+1)(M−M2+1) 。由于低秩问题的求解是NP-hard的,通常通过核范数来近似表示矩阵的秩函数,即低秩优化问题可以建模为
minˆΔH‖ˆΔH‖∗s.t.(YH)p,q=(ˆΔH)p,q,(p,q)∈ΩH (24) 针对上述优化问题,现有算法通常将其定义为半正定规划问题进而采用内点法求解,然而由于双重Hankel变换后的矩阵维度较大,导致半正定规划求解十分困难;另一方面,核范数松弛是对所有奇异值的和求最小化处理,难以准确近似矩阵的秩[33]。对此,本文采用截断核范数代替传统核范数逼近矩阵的秩。截断核范数优化目标仅求解与矩阵秩不相关的较小奇异值的和,间接控制待恢复矩阵秩,从而确保了矩阵低秩估计;进一步地,本文通过交替迭代求解思想,将截断核范数优化问题转换为凸优化问题实现求解。
截断核范数优化模型如下:
minˆΔH‖ˆΔH‖rs.t.(YH)p,q=(ˆΔH)p,q,(p,q)∈ΩH (25) 其中,
‖ˆΔH‖r 表示矩阵ˆΔH 以第r个奇异值作为截断后的核范数,可以表示为‖ˆΔH‖r=∑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 )为截断奇异值个数。这里需要说明的是,矩阵A和BH在后续算法迭代求解过程中会随着
{\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) 由于A和B的秩为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]分析,整个算法处理过程中奇异值分解操作对算法运算复杂度影响最大。其中一个N行M列的矩阵奇异值分解计算复杂度为
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}} 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) 通常来说N1和M2设置为矩阵维度的一半,即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) 相比,为了重构未观测数据样本解决欠采样问题,矩阵填充方法付出了较高的计算负担。因此在未来我们还将研究在保持所提方法性能优势的同时进一步降低计算复杂度的改进算法。4. 仿真验证
本节设置了单/多点目标仿真试验以验证所提方法的有效性和优越性,波形仿真参数如表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窗 需要说明的是,本文研究的是如何实现对单个距离粗分辨单元内目标的高分辨距离-多普勒谱精确重建,因此试验中高分辨距离含义是以所在距离粗分辨单元对应距离为基准的相对值,不代表目标相对雷达的真实距离;速度也是在补偿了雷达自身载体速度后,目标相对雷达运动产生的径向速度。
4.1 单目标场景仿真
本节设置单目标信噪比为10 dB (快时间脉冲压缩后),高分辨距离和速度为0 m, 0 m/s。采用匹配滤波、所提方法估计得到的随机步进频雷达高分辨距离-多普勒谱如图3(b)、图3(c)所示;此外,时频等效满采样数据(信噪比同样也设置为10 dB)估计的高分辨距离-多普勒谱见图3(a)。可以看出,尽管匹配滤波处理可以有效地聚焦并分辨目标,但谱中起伏高旁瓣非常明显;相比之下,所提方法的估计谱与等效满采样数据估计结果非常接近,恢复效果一致,表明该方法可有效解决欠定估计问题。
进一步地,图4展示了不同方法在高分辨距离、速度以及对角线等剖面的点目标响应函数(Point Spread Functions, PSF)。可以看出,尽管所有方法估计的主瓣宽度一样,但与等效满采样数据恢复结果的PSF相比,匹配滤波处理的PSF旁瓣在距离、速度等维度仍增加了18 dB以上,而在对角线方向则增加了28 dB。与之相反,本文提出矩阵填充处理方法可以有效抑制旁瓣,特别是在对角线方向上,其旁瓣水平甚至比等效满采样恢复结果要低,这表明矩阵填充处理可以准确地恢复未观测的数据样本,且未引入额外噪声,间接提升了数据的输入信噪比。
为了量化评估提出方法性能,本文采用峰值旁瓣比(Peak Side Lobe Ratio, PSLR)和积分旁瓣比(Integral Side Lobe Ratio, ISLR)两个指标来评估分析提出方法的旁瓣抑制水平,定义如下。图5和图6展示了在不同输入信噪比下不同方法的PSLR和ISLR结果。其中每个结果是由300次蒙特卡洛仿真试验得到。
\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.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)则展示了通过匹配滤波和本文所提方法恢复数据后估计得到的高分辨距离-多普勒谱。结果表明,传统匹配滤波方法估计得到的频谱中存在严重的旁瓣基底,这直接导致了设置的小信噪比目标被大信噪比目标的随机起伏旁瓣所遮盖,无法有效检测;同时,在估计谱中出现若干个虚警目标。与之相比,经过矩阵填充处理后的估计谱旁瓣可以得到有效抑制。
图8显示了多目标在距离、速度剖面的PSF性能。可以看出,提出方法的聚焦效果良好,可以准确分辨多目标,且主瓣宽度与等效满采样恢复谱几乎一致,表明提出方法没有损失高分辨距离-多普勒的分辨率。尽管提出方法估计谱中仍存在一些随机波动的起伏旁瓣,但其旁瓣平均值也远低于匹配滤波结果,证明了所提方法可以在多目标环境下有效抑制旁瓣。需要注意的是,与时频等效满采样数据情况相比,提出方法的旁瓣抑制水平在目标个数增多时有所下降,这是因为随着目标数量的增加,待恢复矩阵秩也会随之增大,势必会影响算法恢复性能。此外,图8(c)还比较了多目标场景中对角线方向的PSF,结果表明传统匹配滤波方法恢复的弱目标主瓣(归一化后功率–8 dB左右)被淹没在强目标的高功率旁瓣(归一化后平均功率–13 dB)中,无法被正确地检测。相比之下,提出方法可以有效地抑制旁瓣,聚焦弱目标,取得较好的谱估计效果。
最后,本文计算了在输入信噪比10 dB情况下,不同方法的PSLR, ISLR与目标数量的关系,每次结果均进行了300次蒙特卡罗仿真试验,结果如图9和图10所示。可以看出所提方法的旁瓣抑制性能随着目标数量K的增大而下降。这是因为在多目标场景下,双重Hankel重构虽能一定程度上保持矩阵的低秩特性,但无法完全解决因为目标数增加导致矩阵秩快速变大,因此旁瓣抑制水平快速恶化。当然正如前文所说,快时间脉冲压缩预处理可以有效降低距离粗分辨单元的目标数,一定程度上使得待恢复矩阵满足低秩特性,进而保证提出方法的谱估计性能。
4.3 压缩感知方法性能比较仿真
进一步,本文还比较了所提方法与压缩感知方法估计高分辨距离-多普勒谱的性能。其中压缩感知方法优化模型如下,优化算法采用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),可以发现提出方法的聚焦效果更好;与之相反,压缩感知方法能量有所损失。
表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 此外,在处理高维数据时,对于CS方法,每个维度必须被重构为向量形式,处理后再形成高维矩阵,这导致相当大的内存使用,而本文方法将矩阵数据作为一个整体来处理,在减少内存使用方面具有巨大潜力。
5. 实测验证
为了进一步验证方法的有效性,对实测数据进行高分辨距离-速度估计,目标是一个运动的厢货卡车(如图13所示),回波信号通过一部93 GHz毫米波雷达样机采样获得,其中雷达参数如表3所示。
表 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}} 随机均匀分布 试验中首先通过Matlab生成遵循均匀分布的随机跳频序列,并将其存储到雷达样机中。回波信号经预处理(下混频、滤波、采样、脉冲压缩等)后存储到计算机中。采用匹配滤波、压缩感知以及本文所提方法估计的高分辨距离-多普勒谱如图14所示。
由于欠定估计,匹配滤波的谱中产生了较高的起伏旁瓣,导致恢复结果变差。提出方法的结果远优于匹配滤波方法,初步验证了矩阵填充方法在实际系统中的有效性。压缩感知方法虽能抑制旁瓣,但其依赖稀疏基的设计,需要目标的先验信息,对未知场景与目标成像效果差。
我们提取了v=–8.75 m/s的高分辨率距离像,以进一步分析CS和所提方法的估计性能,如图15所示。
显然,所提出的方法拥有相对更为优越的成像结果,具体表现在可以很好地在高分辨率距离像上呈现出目标的大体轮廓,从而提供更多关于该目标的特征信息;此外,所提出的方法结果与匹配滤波恢复结果的主瓣能量、散射点个数等匹配更好。相反,压缩感知恢复结果在目标分段区域没有被完全聚焦,估计的HRRP质量较差。因此,我们可以得出结论,与CS方法相比,本文方法可以获得更好的HRRP聚焦效果,对于目标参数可以实现高精度估计,更有利于后续的目标分类识别。
6. 结语
本文提出了一种基于矩阵填充的随机步进频雷达高分辨率距离-多普勒谱稀疏恢复方法。该方法利用矩阵填充技术恢复未观测的时频数据,解决了随机步进频波形固有的二维欠采样问题。与匹配滤波方法相比,该方法在保持能量、分辨率、模糊度等波形性能不变情况下,能够有效抑制抬高的旁瓣。此外,与压缩感知方法相比,所提方法不仅在非栅格点谱恢复方面取得了更好的性能,并且显著降低了计算负担。仿真和实测数据验证了该方法的有效性和优越性。
-
表 1 仿真目标信息
Table 1. Information of targets in simulation
目标 初始位置(km) 初始速度(m/s) 目标出现场序号 目标终止场序号 1 [150, 130] [–150, –90] 1 200 2 [80, –100] [20, –100] 1 200 3 [–150, 100] [180, –60] 1 200 4 [70, –120] [100, 100] 1 200 5 [55, 15] [20, –200] 1 200 6 [–100, 30] [80, 180] 1 200 表 2 各收发对不同时刻的检测概率
Table 2. Detection probability of radar stations at different times
收发对序号 场次 1~50 51~100 101~150 151~200 收发对1 0.7 0.7 0.5 0.5 收发对2 0.5 0.8 0.7 0.8 收发对3 0.8 0.5 0.8 0.7 -
[1] 万显荣, 易建新, 占伟杰, 等. 基于多照射源的被动雷达研究进展与发展趋势[J]. 雷达学报, 2020, 9(6): 939–958. doi: 10.12000/JR20143.WAN Xianrong, YI Jianxin, ZHAN Weijie, et al. Research progress and development trend of the multi-illuminator-based passive radar[J]. Journal of Radars, 2020, 9(6): 939–958. doi: 10.12000/JR20143. [2] 万显荣, 易建新, 程丰, 等. 单频网分布式外辐射源雷达技术[J]. 雷达学报, 2014, 3(6): 623–631. doi: 10.12000/JR14156.WAN Xianrong, YI Jianxin, CHENG Feng, et al. Single frequency network based distributed passive radar technology[J]. Journal of Radars, 2014, 3(6): 623–631. doi: 10.12000/JR14156. [3] KUSCHEL H. Approaching 80 years of passive radar[C]. 2013 International Conference on Radar, Adelaide, Australia, 2013: 213–217. doi: 10.1109/RADAR.2013.6651987. [4] COLONE F, BONGIOANNI C, and LOMBARDO P. Multifrequency integration in FM radio-based passive bistatic radar. Part I: Target detection[J]. IEEE Aerospace and Electronic Systems Magazine, 2013, 28(4): 28–39. doi: 10.1109/MAES.2013.6506827. [5] POULLIN D and FLECHEUX M. Recent progress in Passive Coherent Location (PCL) concepts and technique in France using DAB or FM broadcasters[C]. 2008 IEEE Radar Conference, Rome, Italy, 2008: 1–5. doi: 10.1109/RADAR.2008.4721009. [6] SLAVOV A, SANDENBERGH S, O’HAGAN D, et al. Multiple FM-based passive bistatic pairs for robust target detection with improved position accuracy[C]. 2022 23rd International Radar Symposium (IRS), Gdansk, Poland, 2022: 332–337. doi: 10.23919/IRS54158.2022.9905050. [7] PLŠEK R, STEJSKAL V, PELANT M, et al. FM based passive coherent radar: From detections to tracks[C]. 2011 Tyrrhenian International Workshop on Digital Communications - Enhanced Surveillance of Aircraft and Vehicles, Capri, Italy, 2011: 123–127. [8] ZHANG Chenqi, WU Yong, WANG Jun, et al. FM-based multi-frequency passive radar system[C]. 2016 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC), Hong Kong, China, 2016: 1–4. doi: 10.1109/ICSPCC.2016.7753668. [9] TUYSUZ B, URBINA J V, and MATHEWS J D. Effects of the equatorial electrojet on FM-based passive radar systems[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(7): 4082–4088. doi: 10.1109/TGRS.2017.2687830. [10] FU Yan, WAN Xianrong, ZHANG Xun, et al. Side peak interference mitigation in FM-based passive radar via detection identification[J]. IEEE Transactions on Aerospace and Electronic Systems, 2017, 53(2): 778–788. doi: 10.1109/TAES.2017.2665079. [11] O’HAGAN D W, KUSCHEL H, UMMENHOFER M, et al. A multi-frequency hybrid passive radar concept for medium range air surveillance[J]. IEEE Aerospace and Electronic Systems Magazine, 2012, 27(10): 6–15. doi: 10.1109/MAES.2012.6373907. [12] ZAIMBASHI A, DERAKHTIAN M, and SHEIKHI A. Invariant target detection in multiband FM-based passive Bistatic Radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2014, 50(1): 720–736. doi: 10.1109/TAES.2013.120248. [13] BONGIOANNI C, COLONE F, and LOMBARDO P. Performance analysis of a multi-frequency FM based passive Bistatic Radar[C]. 2008 IEEE Radar Conference, Rome, Italy, 2008: 1–6. doi: 10.1109/RADAR.2008.4720805. [14] BHARADWAJ L S, SWEETLIN P, and VISHNU O C. FM based passive radar tracking of targets under poor waveforms[C]. 2022 IEEE Microwaves, Antennas, and Propagation Conference (MAPCON), Bangalore, India, 2022: 1599–1603. doi: 10.1109/MAPCON56011.2022.10047271. [15] 何友, 修建娟, 关欣, 等. 雷达数据处理及应用[M]. 3版. 北京: 电子工业出版社, 2013: 108–125.HE You, XIU Jianjuan, GUAN Xin, et al. Radar Data Processing with Applications[M]. 3rd ed. Beijing: Publishing House of Electronics Industry, 2013: 108–125. [16] HELMICK R E and WATSON G A. Interacting multiple model integrated probabilistic data association filters (IMM-IPDAF) for track formation on maneuvering targets in cluttered environments[C]. Proceedings of SPIE 2235, Signal and Data Processing of Small Targets 1994, Orlando, USA, 1994: 460–471. doi: 10.1117/12.179097. [17] CHOI S, CROUSE D, WILLETT P, et al. Multistatic target tracking for passive radar in a DAB/DVB network: Initiation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(3): 2460–2469. doi: 10.1109/TAES.2015.130270. [18] BAEK J, LEE J, SHIM H, et al. Target tracking initiation for multi-static multi-frequency PCL system[J]. IEEE Transactions on Vehicular Technology, 2020, 69(10): 10558–10568. doi: 10.1109/TVT.2020.3012135. [19] DAUN M, NICKEL U, and KOCH W. Tracking in multistatic passive radar systems using DAB/DVB-T illumination[J]. Signal Processing, 2012, 92(6): 1365–1386. doi: 10.1016/j.sigpro.2011.09.005. [20] YI Jianxin, WAN Xianrong, CHENG Feng, et al. Deghosting for target tracking in single frequency network based passive radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(4): 2655–2668. doi: 10.1109/TAES.2015.130424. [21] SHU Kan, YI Jianxin, WAN Xianrong, et al. A hybrid tracking algorithm for multistatic passive radar[J]. IEEE Systems Journal, 2021, 15(2): 2024–2034. doi: 10.1109/JSYST.2020.2994009. [22] YI Jianxin, WAN Xianrong, LEUNG H, et al. MIMO passive radar tracking under a single frequency network[J]. IEEE Journal of Selected Topics in Signal Processing, 2015, 9(8): 1661–1671. doi: 10.1109/JSTSP.2015.2464188. [23] BOZDOĞAN A Ö and EFE M. Track initiation using multiple bistatic range and range rate measurements with multidimensional assignment algorithm[C]. 2010 IEEE 18th Signal Processing and Communications Applications Conference, Diyarbakir, Turkey, 2010: 435–438. doi: 10.1109/SIU.2010.5650348. [24] LIU Zongxiang, ZHU Xiaoping, and HUANG Bingjian. Track initiation technique and its application in PHD filter[C]. 2018 14th IEEE International Conference on Signal Processing (ICSP), Beijing, China, 2018: 822–826. doi: 10.1109/ICSP.2018.8652496. [25] BAR-SHALOM Y, LI X R, and KIRUBARAJAN T. Estimation with Applications to Tracking and Navigation: Theory, Algorithms and Software[M]. New York, USA: John Wiley & Sons, Inc., 2001: 267–295. doi: 10.1002/0471221279. [26] MALANOWSKI M and KULPA K. Two methods for target localization in multistatic passive radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2012, 48(1): 572–580. doi: 10.1109/TAES.2012.6129656. [27] SCHUHMACHER D, VO B T, and VO B N. A consistent metric for performance evaluation of multi-object filters[J]. IEEE Transactions on Signal Processing, 2008, 56(8): 3447–3457. doi: 10.1109/TSP.2008.920469. -