1.   引言

近年来，无人机的运用越来越广泛，在遥感测绘、抢险救援和打击犯罪等领域中发挥了重要的作用，然而，在许多有益应用前景的反面，无人机在机场净空区内非法升空等一系列威胁公共安全的事件时有发生，随着以无人机为代表的低慢小飞行器的应用快速增加，其带来的威胁和风险也逐渐显露。这类目标的飞行高度低、运动速度较慢、雷达反射截面积(Radar Cross Section, RCS)小，导致其雷达回波信噪比(Signal Noise Ratio, SNR)低、且多普勒频率较低，很难将目标和背景区分开来。无人机类目标的特点是同时存在整体径向运动和内部微动，微动是指物体的某些部分相对于物体质心的平动存在小幅振动、旋转和其他高阶运动的现象^[1–4]，在频域表现为微多普勒特征。微多普勒现象是Chen^[1]发现的一种多普勒谱调制现象，由于目标内部部件存在振动或旋转等运动，而其部件始终与目标主体连接，导致它们在雷达视角下呈现出速度的周期变化，所以微多普勒特征会表现出交叉、重叠等周期性现象^[5]。与径向运动带来的多普勒特征不同，微多普勒特征是对多普勒频率的高阶调制。随着雷达分辨率的提高，联合利用多普勒和微多普勒特征，提升低慢小目标的检测性能成为当前目标检测一种可能的技术方向^[6]。

现有对于微多普勒参数估计的研究主要分为两种方法^[5]，将回波信号经过时频变换后，一种方法借鉴光学图像处理中的提取曲线的经典思路，在时频图中搜索局部峰值点得到频率脊^[7–10]，使用一定方法将所有时刻的频率脊做排列，即可得到微动目标各部件的运动速度变化曲线，从而反演出微多普勒参数；另一种方法则使用逆拉登变换^[11]、Hough变换^[12,13]等方法将时频图转换到某一参数空间中，通过在参数空间检测峰值点，来估计微多普勒参数。使用微多普勒信息进行目标检测的研究还比较少。

将多普勒和微多普勒统一到1个数学框架下，是有效利用二者能量进行检测的基础。随机有限集(Random Finite Sets, RFS)滤波器是一种序贯贝叶斯的状态估计方法，由Mahler^[14]于21世纪初提出，Vo等人^[15–21]基于Mahler的工作，提出了一系列适用于多目标跟踪的滤波器，其观测模型主要分为两类^[22]，即检测器输出测量模型^[17]和强度测量模型^[23]，前者需要将回波先经过一个低门限检测器，这种方法会造成一些信息丢失，后者将回波直接输入滤波器，基于目标的运动模型和目标回波的观测模型，可以完整使用目标回波的所有信息。势均衡多目标多伯努利(Cardinality Balanced Multi-target Multi-Bernoulli, CBMeMBer)滤波器^[23]是随机有限集滤波器中的一种，可以在检测目标的同时获取目标状态的估计，其具备良好的多目标检测性能，同时计算量较小。

针对现有低慢小目标检测只使用目标径向运动产生的多普勒信息，丢失了一部分由于目标微动而分散的能量，导致目标检测性能下降的问题，本文将微多普勒参数估计问题与微动目标检测问题结合起来，对多普勒和微多普勒进行联合建模，提出一种基于CBMeMBer滤波器的无人机目标联合检测与估计方法。所提算法在传统点测量目标模型的基础上对多普勒和微多普勒进行联合建模，利用微多普勒带来的信息增量，提高了滤波器的检测灵敏度，同时可以获得目标位置以及微动状态的估计，具有适应多无人机目标检测与估计的潜力。仿真结果表明，本文所提算法充分利用两类信息蕴含的特征，有效实现了目标信号能量的积累，检测性能相较传统方法有所提升。

2.   目标多普勒与微多普勒联合建模

假设一无人机目标主体存在低速运动，内部微动部件存在旋转运动，场景如图1所示，可以对目标雷达回波进行建模。

图 1  实验场景

下载: 全尺寸图片 幻灯片

2.1   目标多普勒模型

假设目标以恒定速度远离雷达运动，其回波会产生多普勒频移。雷达发射信号可表示为

其中， ${t_{\text{f}}}$表示脉冲内快时间， $\tau$表示脉冲宽度。

设目标为点目标，且在一个脉冲内没有走动，与雷达之间的距离为 $r({t_{\text{s}}})$，则目标一个脉冲的回波可表示为

其中， ${t_0}$表示目标在慢时间维的回波时延， ${t_{\text{s}}}$为脉冲间慢时间， $\lambda$为雷达载波波长， $v$为目标相对雷达的运动速度，远离雷达时为正。

2.2   目标微多普勒模型

假设目标内部存在一个双翼螺旋桨，桨叶关于旋转中心呈轴对称，其旋转轴方向与雷达视线方向垂直，桨叶初始方向与目标运动速度方向一致。旋转目标满足理想散射点模型^[3]，一片桨叶可由一个位于桨叶尖端的散射中心等效，将位于桨叶上的散射点投影到雷达平面如图2所示。

图 2  雷达-旋转散射点关系示意图

下载: 全尺寸图片 幻灯片

图2中， $O$表示近似的雷达波束中心点， ${O^\prime }$表示目标的旋转中心， $D$为 $O{O^\prime }$的长度，即雷达与目标之间的距离。目标的旋转部件散射点为 $P$，旋转半径为 $r$，旋转频率为 ${f_{\text{r}}}$，散射点运动初始相位为 ${\theta _0}$， $t$时刻的相位为 ${\theta _t} = 2\pi {f_{\text{r}}}t + {\theta _0}$，散射点 $P$到雷达的距离为

当目标与雷达的距离 $D$满足 ${r \mathord{\left/ {\vphantom {r D}} \right. } D} \to 0$时

假设 $D$为一常量，对式(4)求导，可以得到旋转部件运动速度

根据式(5)，可以得到目标旋转部件的多普勒频移为

2.3   目标多普勒与微多普勒联合模型

式(5)中的求导将式(4)中的 $D$作为常量忽略了，如果目标存在运动速度 $v$，那么式(3)可以写为

此时式(4)中的近似关系仍然满足，有

假设目标的一片桨叶可使用一个散射点近似，根据式(8)，单个桨叶的回波相位为

假设目标内部的螺旋桨存在 $B$个桨叶，则第 $b$个桨叶的微多普勒调制相位应为

则微动目标的雷达回波可表示为

3.   基于CBMeMBer滤波器的微弱目标联合检测估计

基于前述分析，对式(11)中回波慢时间维进行傅里叶变换

可以得到回波的距离-多普勒图，假设目标与前述假设一致，存在 $B = 2$片桨叶，则距离-多普勒图应存在3个谱峰，其频率如式(13)所示，其中， $t$表示该帧的起始时刻

3.1   目标运动模型

基于强度测量模型，假设目标可由式(14)中向量表示

其中，下标 $k$表示当前时刻， ${r_k}$表示目标与雷达的距离， ${v_k}$表示目标与雷达的相对速度(远离为正)， ${\rho _k}$表示目标的平均RCS系数， ${V_k}$表示最大相对速度差， ${\omega _k} = 2\pi {f_{{\text{r}},k}}$表示目标螺旋桨的旋转角速度， ${\theta _k}$表示微多普勒调制相位， ${B_k}$表示目标螺旋桨的桨叶个数。

向量中的 ${V_k}$表示距离-多普勒图中目标多普勒最大展宽所对应的最大速度变化值，而并非实际物理意义对应的旋转半径 ${r_k}$，这是因为目标在距离-多普勒图中的多普勒展宽并非由 ${r_k}$一个变量决定，其与旋转频率 ${f_{{\text{r}},k}}$也关联，而旋转频率又与目标距离维的运动周期相关联，选用 ${r_k}$作为滤波器状态输入容易造成滤波器无法收敛，故这里选择了 ${V_k}$和 ${\omega _k}$作为控制微多普勒频率调制的关键参数。

目标状态转移函数为

其中

图3展示了目标在距离-多普勒图中的测量点与其可能的运动轨迹，其横轴与纵轴分别为距离轴与频率轴。 $t = 0$时框选的1个“ $\times$”与 $t = {1 \mathord{\left/ {\vphantom {1 {4{f_{\text{r}}}}}} \right. } {4{f_{\text{r}}}}}$时框选的3个“ $\times$”分别表示这两个时刻目标回波变换得到的谱峰。其中， $t = 0$时，目标的两个桨叶旋转到雷达视线所在直线上，此时，桨叶的速度方向与雷达视线垂直，其速度在雷达视角为 $0$，所以呈现出与目标主体相同的运动速度，假设该时刻微多普勒调制相位 $\theta = 0$。 $t = {1 \mathord{\left/ {\vphantom {1 {4{f_{\text{r}}}}}} \right. } {4{f_{\text{r}}}}}$时，目标的两个桨叶旋转到雷达视线的垂直线上，此时一个桨叶向远离雷达方向运动，一个桨叶向靠近雷达方向运动，所以其速度呈现出关于目标主体速度对称的特征，该时刻 $\theta = {\pi \mathord{\left/ {\vphantom {\pi 2}} \right. } 2}$。

图 3  微多普勒谱与状态参数关系示意图

下载: 全尺寸图片 幻灯片

各个桨叶虽然存在与主体不同的速度，但其仍然是目标上的一部分，所以其速度会在 $t = {1 \mathord{\left/ {\vphantom {1 {4{f_{\text{r}}}}}} \right. } {4{f_{\text{r}}}}}$之后渐渐又接近主体速度，并呈现周期变化，其最大频谱展宽为 $\Delta f = - {{4V} \mathord{\left/ {\vphantom {{4V} \lambda }} \right. } \lambda }$。假设目标运动速度不变，则目标螺旋桨旋转一个周期时目标相对雷达的运动距离为 $\Delta r = {v \mathord{\left/ {\vphantom {v {{f_{\text{r}}}}}} \right. } {{f_{\text{r}}}}}$。

3.2   目标观测模型

假设距离-多普勒图的分辨率为 ${N_r} \times {N_v}$，其距离维分辨率为 ${N_r}$，速度维分辨率为 ${N_v}$， $k$时刻测量数据可表示为 ${{\boldsymbol{Y}}_k} = {\{ {y_{1,k}}, \cdots ,{y_{m,k}}\} _{{N_r} \times {N_v}}}$，一个状态矢量 ${{\boldsymbol{x}}_k}$会影响集合 $T({{\boldsymbol{x}}_k})$中所包含的所有测量，即 ${y_{l,k}},\;l \in T({{\boldsymbol{x}}_k})$， ${y_{l,k}}$表示一个像素， $l$为该像素的绝对索引， $T$表示状态空间向观测空间点的映射函数。像素绝对索引与该像素所对应的距离和速度之间存在映射 $L(r,v) = l$，其逆映射为 ${L^{ - 1}}(l) = r,v$。

本文使用点测量假设，一个状态对应 $B + 1$个测量点，分为目标主体的一个测量点，以及目标螺旋桨产生的 $B$个测量点，各散射点的强度会因为雷达视线角不同而变化，数量上一般不会变化。

假设状态 ${{\boldsymbol{x}}_k}$共对应 ${B_k} + 1$个测量点，对于某个测量点集合 ${t_i} \subset T({{\boldsymbol{x}}_k}),i \in [1,{B_k} + 1]$，其距离为 ${r_k}$，速度为

测量点集合 ${t_i}$包含 ${r_k}$和 ${v_{k,i}}$周围一定范围内的全部像素，对于集合内的一个像素 ${y_{j,k}} \in {t_i}$，假设该像素绝对索引所对应的距离和速度为 $({r_{j,k}},{v_{j,k}})$，则

其中， $y_{j,k}^{\text{I}} = {\text{Re}}\left( {{y_{j,k}}} \right)$和 $y_{j,k}^{\text{Q}} = {\text{Im}}\left( {{y_{j,k}}} \right)$分别表示像素 ${y_{j,k}}$的实部和虚部， ${k_{\text{s}}}$为雷达系统增益， ${d \mathord{\left/ {\vphantom {d 2}} \right. } 2}$为信号方差， $\zeta _{\text{s}}^2$为噪声方差，目标回波信噪比为

一个状态所影响的集合 $T({{\boldsymbol{x}}_k})$为所有 ${B_k} + 1$个测量点子集的并集，即

对于状态影响范围外的像素 ${y_{j,k}} \notin T({{\boldsymbol{x}}_k})$，有

3.3   目标测量似然函数

对于某个状态 ${{\boldsymbol{x}}_k}$，观测似然比^[23]可表示为

其中， $d$需要根据式(18)计算得到。

3.4   目标联合检测估计滤波器构建

假设一个定义在状态空间 $\chi$上的伯努利RFS为空的概率为 $1 - r$，存在单个状态 ${\boldsymbol{x}}$的概率为 $r$， ${\boldsymbol{x}}$服从概率密度函数 $p$所定义的分布。

一个定义在 $\chi$上的多伯努利RFS是固定个数为 $M$的独立伯努利RFS ${{\boldsymbol{X}}^{(i)}}$的集合。其概率生成泛函(Probability Generating Functional, PGFl)定义为

其中， ${r^{(i)}}$和 ${p^{(i)}}$为 ${{\boldsymbol{X}}^{(i)}}$的单项存在概率和状态的概率密度函数。

则多伯努利RFS的概率密度可表示为

根据文献[23]的推导，CBMeMBer滤波可分为预测和更新两个过程，其序贯蒙特卡罗(Sequential Monte Carlo, SMC)实现可用算法1表示，其中，所有带有 $k$或 $k - 1$下标的变量分别表示 $k$时刻或者 $k - 1$时刻的该变量值。

算法1 无人机目标CBMeMBer滤波器
输入：当前时刻测量 ${{\boldsymbol{Y}}_k} = {\{ {y_{1,k}}, \cdots ,{y_{m,k}}\} _{{N_r} \times {N_v}}}$、上一时刻滤波状态 $r_{k - 1}^{(i)},\{ w_{k - 1}^{(i,j)},{\boldsymbol{x}}_{k - 1}^{(i,j)}\} _{j = 1}^{L_{k - 1}^{(i)}},\quad i \in [1,{M_{k - 1}}]$
输出：本时刻滤波状态 $r_k^{(i)},\{ w_k^{(i,j)},{\boldsymbol{x}}_k^{(i,j)}\} _{j = 1}^{L_k^{(i)}},\quad i \in [1,{M_k}]$
(1) 根据建议分布采样新生粒子状态 ${\boldsymbol{x}}_{\varGamma ,k}^{(i,j)}{\text{～}}b_k^{(i)}( \cdot |{{\boldsymbol{Y}}_k}),j = 1,2, \cdots ,L_{\varGamma ,k}^{(i)}$
(2) 计算新生粒子权重 $w_{\varGamma ,k}^{(i,j)} = \dfrac{{{p_{\varGamma ,k}}({\boldsymbol{x}}_{\varGamma ,k}^{(i,j)})}}{{b_k^{(i)}({\boldsymbol{x}}_{\varGamma ,k}^{(i,j)}|{{\boldsymbol{Y}}_k})}}$ 　(3) 新生粒子权重归一化 $\tilde w_{\varGamma ,k}^{(i,j)} = \dfrac{{w_{\varGamma ,k}^{(i,j)}}}{{\displaystyle\sum\limits_{j = 1}^{L_{\varGamma ,k}^{(i)}} {w_{\varGamma ,k}^{(i,j)}} }}$
(4) 根据建议分布采样预测粒子状态 ${\boldsymbol{x} }_{P,k|k - 1}^{(i,j)}{\text{～} }q_k^{(i)}( \cdot |{\boldsymbol{x} }_{k - 1}^{(i,j)},{ {\boldsymbol{Y} }_k}),j = 1,2, \cdots ,L_{k - 1}^{(i)}$
一般设 $q_k^{(i)}( \cdot |{\boldsymbol{x}}_{k - 1}^{(i,j)},{{\boldsymbol{Y}}_k}) = {f_{k|k - 1}}( \cdot |{\boldsymbol{x}}_{k - 1}^{(i,j)})$， ${f_{k|k - 1}}( \cdot |{\boldsymbol{x}}_{k - 1}^{(i,j)})$由式(15)给出
(5) 根据式(16)为 ${\boldsymbol{x}}_{P,k|k - 1}^{(i,j)}$添加噪声
(6) 更新预测粒子权重 $w_{P,k|k - 1}^{(i,j)} = \dfrac{{w_{k - 1}^{(i,j)}{f_{k|k - 1}}({\boldsymbol{x}}_{P,k|k - 1}^{(i,j)}|{\boldsymbol{x}}_{k - 1}^{(i,j)}){p_{S,k}}({\boldsymbol{x}}_{k - 1}^{(i,j)})}}{{q_k^{(i)}({\boldsymbol{x}}_{P,k|k - 1}^{(i,j)}|{\boldsymbol{x}}_{k - 1}^{(i,j)},{{\boldsymbol{Y}}_k})}}$ 　(7) 预测粒子权重归一化 $\tilde w_{P,k|k - 1}^{(i,j)} = \dfrac{{w_{P,k|k - 1}^{(i,j)}}}{{\displaystyle\sum\limits_{j = 1}^{L_{k - 1}^{(i)}} {w_{P,k|k - 1}^{(i,j)}} }}$
(8) 根据上一时刻滤波状态计算本时刻单个伯努利项存在概率 $r_{P,k|k - 1}^{(i)} = r_{k - 1}^{(i)}\displaystyle\sum\limits_{j = 1}^{L_{k - 1}^{(i)}} {w_{k - 1}^{(i,j)}} {p_{S,k}}({\boldsymbol{x}}_{k - 1}^{(i,j)})$
(9) 根据新生模型计算新生多伯努利密度存在概率 $r_{\varGamma ,k}^{(i)}$
(10) 将新生粒子合并到预测粒子中 $\{ w_{k|k - 1}^{(i,j)},{\boldsymbol{x} }_{k|k - 1}^{(i,j)}\} _{j = 1}^{L_{k - 1}^{(i)} } = \{ w_{\varGamma ,k|k - 1}^{(i,j)},{\boldsymbol{x} }_{\varGamma ,k|k - 1}^{(i,j)}\} \cup \{ w_{P,k|k - 1}^{(i,j)},{\boldsymbol{x} }_{P,k|k - 1}^{(i,j)}\}$
(11) 使用式(22)计算后验粒子权重 $\hat w_k^{(i,j)} = w_{k|k - 1}^{(i,j)}{g_{{y_k}}}({\boldsymbol{x}}_{k|k - 1}^{(i,j)})$ 　(12) 计算粒子权重总和 $\varrho _k^{(i)} = \displaystyle\sum\limits_{j = 1}^{L_{k|k - 1}^{(i)}} {\hat w_k^{(i,j)}}$ 　(13) 计算后验多伯努利密度存在概率 $r_k^{(i)} = \dfrac{{r_{k|k - 1}^{(i)}\varrho _k^{(i)}}}{{1 - r_{k|k - 1}^{(i)} + r_{k|k - 1}^{(i)}\varrho _k^{(i)}}}$
(14) 后验粒子权重归一化 $w_k^{(i,j)} = \dfrac{{\hat w_k^{(i,j)}}}{{\varrho _k^{(i)}}}$
(15) 后验粒子状态为上一时刻预测状态 ${\boldsymbol{x}}_k^{(i,j)} = {\boldsymbol{x}}_{k|k - 1}^{(i,j)}$
(16) 根据粒子权重分布对 $\{ w_{k - 1}^{(i,j)},{\boldsymbol{x}}_{k - 1}^{(i,j)}\} _{j = 1}^{L_{k - 1}^{(i)}}$重采样，然后将权重设置为相同值并归一化

下载: 导出CSV 

| 显示表格

在使用算法1对测量进行滤波后，需要计算当前时刻的估计，可以使用

4.   仿真实验

为了验证本文所提算法的有效性，仿真使用高重频脉冲多普勒(Pulse Doppler, PD)雷达对远距离微动目标进行检测估计。为了简化计算量，仿真以PD雷达距离-多普勒图作为算法输入，其参数设置如表1所示。总共40个时刻，目标自0时刻出现，一直存在至最后一个时刻，其部分参数如表2所示。

表 1  仿真雷达参数

参数	值
帧间隔	ε=1 s
最大不模糊距离	200 m
最大不模糊速度	20 m/s
帧分辨率	Nr×Nv=256×512
距离测量误差	$\zeta_r=3 \;{\rm{m}}$
速度测量误差	$\zeta_v=5 \;{\rm{m/s}}$
强度测量误差	$\zeta_{\rm{s}}=0.3$
雷达系统增益	ks=1×1019

下载: 导出CSV 

| 显示表格

表 2  仿真目标参数

参数	值
初始距离	r0=20 010 m
运动速度	v=1 m/s
旋转角速度	$\omega$=0.06π rad/s
初始旋转相位	θ0=0
桨叶个数	B=2
测量强度	ρ=1

下载: 导出CSV 

| 显示表格

仿真场景如图4所示，雷达静止不动，目标主体存在径向运动，同时目标内部存在旋转运动。

图 4  仿真场景示意图

下载: 全尺寸图片 幻灯片

滤波器使用SMC实现，最大粒子数1 500，最小粒子数450，轨迹存活概率恒定为0.9，初始距离分布为 ${r_{{\text{init}}}}{\text{～}}\mathcal{N}\left( {20\;010,{2^2}} \right){\text{m}}$，初始速度分布为 ${v_{{\text{init}}}}{\text{～}}\mathcal{N}\left( {1,{{0.4}^2}} \right){\text{m}}/{\text{s}}$，初始目标强度分布为 ${\rho _{{\text{init}}}}{\text{～}} \mathcal{N}\left( {1,{{0.1}^2}} \right)$，其中， $\mathcal{N}$表示正态分布。为了在相同条件下对比本算法与不使用微多普勒信息时滤波器的性能，假设本算法滤波器已知目标微多普勒信息。

根据仿真设置，可以生成各时刻目标的距离-多普勒图，第8帧与第16帧如图5所示。从图5可以清晰地看到由目标所产生的3个多普勒谱峰，这表明在同一距离上出现了3个不同速度的散射点。对比图5(a)和图5(b)两幅子图可以观察到，3个峰值点中，有1个峰值点的速度始终不变，两个峰值点的速度变化呈现出周期性调制的特性。由此可以得到与式(13)一致的结论，即目标螺旋桨的回波呈现出周期性的多普勒频率调制，且当其旋转轴与雷达视线垂直时，该调制函数为正弦函数。

图 5  第8帧与第16帧距离-多普勒图

下载: 全尺寸图片 幻灯片

采用本文的方法对无人机目标进行检测估计，滤波得到的目标状态估计结果如图6所示，从中可以看出，滤波器可以准确估计出各时刻目标的距离、速度信息，并且能够同时准确估计目标的微多普勒信息。

图 6  2dB信噪比下滤波器估计结果与真值

下载: 全尺寸图片 幻灯片

为了验证本文算法在不同信噪比下的检测性能，可以调整 ${k_{\text{s}}}$使目标在相同距离下的信噪比变化，在不同信噪比下进行100次蒙特卡罗实验，统计目标存在情况，可以得到目标检测概率随信噪比变化曲线，如图7中圈号实线所示。可以看到，随着信噪比增加，滤波器对目标的检测概率逐渐趋近于1。

图 7  两种算法的检测概率随信噪比变化曲线

下载: 全尺寸图片 幻灯片

为了对比本文方法的性能，同时给出仅使用多普勒特征进行目标检测估计的检测性能，如图7中加号虚线所示。通过对比可以看到，由于微多普勒特征的加入，算法可以使用的信息增加了，在回波信噪比为1 dB以下时的检测概率相较仅使用多普勒信息的情况下更高，在检测概率为0.92时，检测灵敏度提高了2 dB左右。

为了分析各信噪比下滤波器的估计性能，可以计算滤波器估计结果与真值之间的最优子模式分配(Optimal Sub Pattern Assignment, OSPA)距离^[14]，对所有时刻的OSPA距离取平均，可以得到如图8所示平均OSPA距离随信噪比变化曲线。可以看到，随着信噪比增加，平均OSPA距离不断减小，在0.5 dB开始趋于稳定，这说明算法在该信噪比下趋于收敛。

图 8  本文算法的平均OSPA误差随信噪比变化曲线

下载: 全尺寸图片 幻灯片

通过改变假设的桨叶个数 $B$，可以得到如图9所示不同桨叶个数下检测概率随信噪比变化曲线，可以观察到，随着桨叶数量的增加，在低信噪比下的目标检测概率整体上呈现增加的趋势，这说明桨叶个数的增加会带来额外的信息增量。

图 9  不同桨叶数 $B$下检测概率随信噪比变化曲线

下载: 全尺寸图片 幻灯片

图10所示为不同旋转速度 $\omega$下检测概率随信噪比变化曲线，可以看到，虽然旋转速度 $\omega$有所变化，但是各信噪比下的检测概率整体上呈现出近似的变化。这说明旋转速度在一定范围内变化时，不会对检测概率带来非常大的影响。

图 10  不同旋转速度 $\omega$下检测概率随信噪比变化曲线

下载: 全尺寸图片 幻灯片

为了对比本文算法的性能，尝试将粒子滤波器中粒子数分别设置为700, 1 000, 1 500，运行100帧，计算每一帧平均消耗时间，结果如表3所示。在不同粒子数设置情况下，对比不同信噪比下算法检测概率，如图11所示。计算机使用Windows 10操作系统，CPU为AMD R4750U，内存为48 GB，运行MATLAB R2022b实验软件。

表 3  两种算法运行 100 次平均每帧使用时间

算法	最大粒子数	时间(s)
仅使用多普勒特征	1 500	0.068
本文算法	700	0.066
本文算法	1 000	0.101
本文算法	1 500	0.143

下载: 导出CSV 

| 显示表格

图 11  不同的最大粒子数下检测概率随信噪比变化曲线

下载: 全尺寸图片 幻灯片

结合表3与图11，可以看到，本文算法在3种粒子数设置下的检测概率差别不大，在检测概率0.9时，检测灵敏度提高了2 dB左右。在粒子数设置为700时，计算性能优于仅使用多普勒信息的情况。

5.   结束语

本文提出联合多普勒和微多普勒的目标检测前跟踪方法，通过对微动目标所产生的多普勒和微多普勒现象进行分别分析，建立了两者的联合信号模型，通过选择合适的参数，计算了滤波器所需要的运动方程，并分析了由目标状态得到目标测量的观测模型；对本文所提滤波器进行了100次蒙特卡罗实验，给出了检测概率性能曲线，并对滤波结果做了误差分析，对比了在不同桨叶数、旋转速度下滤波器的检测性能，验证了滤波器的估计准确性。

与仅使用多普勒信息的传统算法相比，本文算法得益于充分利用了多普勒和微多普勒两个特征，在检测概率为0.9时的检测灵敏度提高了2 dB左右，算法在信噪比为0.5 dB时能够稳定实现目标检测，同时OSPA误差趋近于0.1，估计精度较高。后续工作可以根据实测数据验证本文算法的有效性，针对实际可能遇到的回波信号强闪烁问题，修改观测模型。