Joint Selection of Transmitters and Receivers in Distributed Multi-input Multi-output Radar Network for Multiple Targets Tracking

1.   引言

目标的振动、转动等微动产生的微多普勒效应包含了目标的结构和运动信息，常用于目标的分类和识别^[1–3]。目前，基于外辐射源雷达微多普勒效应目标分类和识别的研究还处于起步状态。外辐射源雷达是一种利用非合作照射源进行目标探测和分类识别的新体制雷达系统，其自身不辐射电磁能量，具有节约频谱资源，隐蔽性好，设备规模小，易于部署和组网等特点^[4–6]。在微多普勒效应目标分类和识别方面，外辐射源雷达表现出得天独厚的优势：(1)收发分置可实现空间分集，有效避免探测盲区。(2)第三方辐射源多为连续波，长时间相干积累可记录多个连续的回波闪烁，同时有利于提高对低雷达散射截面积(Radar Cross-Section, RCS)微动目标的探测与分类识别能力。(3)对微多普勒特征的提取不要求高距离分辨率，参数估计不受第三方辐射源带宽的限制^[7,8]。

针对微多普勒效应参数估计问题。文献[9,10]中依据微动目标正弦特征曲线，利用 Hough变换，在参数域中进行多维搜索提取出微动曲线进行参数估计。文献[11,12]通过正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法进行稀疏逼近实现了微动目标参量的估计。上述方法均具有较好的鲁棒性，但由于估计参量维数较高导致计算量巨大。文献[13]利用微动目标在时频域的周期性，采用循环相关系数方法，实现了目标微动周期的估计，但信号周期较长时计算量急剧增加。文献[14]计算了信号的高阶矩函数，通过检测在不同时延下，高阶虚函数部分傅里叶变换累计结果的峰值位置，快速获得目标的旋转速率，相比于图像处理方法和OMP分解方法，计算复杂度较小，但抗噪性能差。而外辐射源雷达所利用的第三方辐射源多为连续波信号，其发射波形不可控，信号能量主要覆盖地面，杂波环境复杂且对空中目标增益低，利用长时间相干积累来提高处理增益会带来数据量巨大的挑战。上述因素决定了外辐射源雷达参数估计方法需要有良好的抗噪性能且计算量要小。

直升机旋翼旋转时对雷达信号产生周期性调制，当叶片发生镜面反射时，旋翼回波出现峰值，即回波闪烁。闪烁信号在时频图像中表现为一定宽度的频率带，且闪烁时间、闪烁间隔与直升机旋翼微动参数密切相关。针对外辐射源雷达参数估计问题，本文结合上述时频域中闪烁信号的特点，通过时频分析和正交匹配追踪算法实现了直升机旋翼微动参数的估计。本文首先给出了外辐射源雷达直升机旋翼微动信号模型，其次介绍了如何在时频图中提取出闪烁信号参数及正交匹配追踪算法对直升机旋翼微动参数的估计，最后仿真和实测证明了本文方法的有效性。

2.   外辐射源雷达直升机旋翼回波模型的建立

直升机旋转叶片与外辐射源雷达的位置关系如图1所示。以直升机旋转叶片的中心点为原点 $o$ ，旋转叶片平面为 $xy$ 面， $x$ 轴平行于发射站与接收站所在直线，建立空间坐标系 $(x,y,z)$ 。直升机相对于发射站和接收站距离为 ${r_T}\;,\,$ ${r_R}$ ，方位角为 $\gamma ,\,$ $\alpha$ ，仰角为 ${\beta _T}\;,\,$ ${\beta _R}$ $(\cos {\beta _T} \approx \cos {\beta _R} = \cos \beta )$ 。叶片上某一散射点 $p$ 到原点 $o$ 距离为 ${l_P}$ ，方位角为 ${\varphi _t}$ 。

图  1  外辐射源雷达直升机旋翼回波模型

Figure  1.  Model of helicopter rotors echo for passive radar

下载: 全尺寸图片 幻灯片

假设直升机平动得到补偿。在 $t$ 时刻，从发射站经散射点 $p$ 到接收站的距离为：

其中 ${{{R}}_T}\;,\,$ ${{{R}}_R},\,$ ${{{R}}_P}(t)$ 分别为发射站、接收站、散射点 $p$ 在坐标系 $xyz$ 中的位置矢量。

参考文献[1]中单基地直升机建模，将叶片看作线模型，外辐射源雷达直升机旋翼回波可表示为：

其中，

${f_r}$ 为叶片转速， $L$ 为叶片长度， $N$ 为叶片数量，整数 $k$ $(0 < k \le N)$ 表示第 $k$ 个叶片， ${\varphi _0}$ 为叶片初相， $\lambda$ 为照射源信号波长。

由式(3)得第 $k$ 个叶片引起的瞬时多普勒频移为：

由式(2)可知时域信号幅值受 $\sin\!{\rm c}$ 函数调制，结合式(3)知当 ${\varphi _k}(t)$ 满足式(6)时， ${\phi _k}(t) = 0$ ，时域信号幅值最大，此刻即时域闪烁。

由式(2)知连续两个闪烁之间的时间间隔为：

3.   直升机旋翼参数估计

3.1   时频分析提取闪烁信号参数

直升机旋翼回波的微多普勒呈非线性变化，通过对目标回波信号进行时频分析能够揭示信号频率的时变特性。短时傅里叶变化(Short-Time Fourier Transform, STFT)计算简单，且不产生交叉项。对直升机旋翼回波信号 $s(t)$ 进行STFT到时频域

其中， $w(t)$ 为窗函数。旋翼微多普勒效应特征曲线为正弦曲线，对应时域闪烁出现的时刻出现垂直于时间横轴的频率带，即时频域“闪烁”^[15]。

图2为直升机旋翼回波的时频图。当直升机旋翼的叶片数为奇数时(图2(a))，时频域中正负多普勒“闪烁”交替出现；若旋翼叶片数为偶数(图2(b))，则是同时出现。

图  2  旋翼回波信号时频分析

Figure  2.  Time-frequency analysis of rotors echo

下载: 全尺寸图片 幻灯片

设时频域中正频率“闪烁”发生的时间为 ${t_0}$ ，由式(5)和式(6)知 ${t_0}$ 满足：

由式(4)和式(9)得第 $k$ 个叶片初相与叶片数量的关系：

其中

由于时频图像中闪烁信号频率带垂直于时间横轴，对正频率轴数据幅值进行累加计算，并判断累加后数据局部峰值点，可得到时频域中正频率“闪烁”发生的时间。同样，对负频率轴数据幅值进行累加计算得到时频域中负频率“闪烁”发生的时间。相应的也可得到闪烁间隔。

由式(7)知，闪烁间隔与旋翼转速、叶片数量密切相关。由式(10)知，闪烁发生的时间与叶片初相、叶片数量、整数 $k$ 密切相关。因此，可根据得到的闪烁间隔，用叶片数量表示出旋翼转速。根据得到的闪烁时间，用叶片数量、整数 $k$ 表示出第 $k$ 个叶片初相。

3.2   OMP估计直升机旋翼参数

由式(2)知时域回波信号可分解为：

其中， ${{{g}}_m}$ 为第 $m$ 个原子， ${{D}}$ 为以原子为列张成的字典矩阵 ${{D}} = [{{{g}}_1}\ {{{g}}_2}\ {{{g}}_3} ·\!·\!· {{{g}}_M}]$ $\in {{{C}}^{{N_t} \times M}},\,$ $M$ 为原子个数， ${N_t}$ 为时间 $t$ 离散后的取值个数， $\varLambda$ 为要估计的参量， ${c_m}$ 为原子系数， ${α} \in {{{C}}^M}$ 为系数矢量，是稀疏的。可转化最优 ${l_0}$ 范数问题进行稀疏向量求解。OMP常用于求解此类问题，通过构建字典矩阵，不断选定与信号最匹配的原子进行稀疏逼近^[16]。OMP将字典矩阵中原子正交化保证了迭代的最优性。

由式(2)知直升机旋翼回波信号由参数 $({f_r},L,{\varphi _0},N,k)$ 确定。利用叶片数量 $N$ 、整数 $k$ $(0 < k \le N)$ 与旋翼转速和叶片初相的关系式(7)和式(10)，时域回波可转化为参数 $(L,N,k)$ 来表示。设时间采样点数 ${N_t}$ ，目标回波为 ${N_t} \times 1$ 的矩阵。确定待估参数的取值范围并离散化，叶片长度取值： $L \in ({L_1}, ·\!·\!· ,{L_r}, ·\!·\!· ,{L_{{N_L}}})$ ，叶片数量 $N$ 的可能取值为： $N \in ({N_1}, ·\!·\!· ,{N_p}, ·\!·\!· ,{N_{{N_N}}})$ ，整数 $k$ 的取值为 $k \in (1, ·\!·\!· ,{k_q}, ·\!·\!· ,{k_{{N_q}}})$ $\;({k_{{N_q}}} \le {N_{{N_N}}})$ 。

由OMP算法原理可知，字典中的原子可按照待分解信号的内在特性来构造^[16]。根据微动目标的时域回波表达式(2)，第 $m$ 个原子可表示为：

其中

并对原子集里的每个原子进行能量归一化：

其中， ${\left\| \cdot \right\|_F}$ 表示矩阵的F范数。

将5参量 $({f_r},L,{\varphi _0},N,k)$ 的估计转换为3参量 $(L,N,k)$ 估计， ${N_L},\,$ ${N_N},\,$ ${N_k}$ 分别为 $L,\,$ $N,\,$ $k$ 的取值个数，由于常见直升机主旋翼叶片数量为：3片、5片、7片(奇数)，2片、4片、8片(偶数)， ${N_N},\,$ ${N_k}$ 较小，降低字典维数为： ${N_L} \times {N_k} \times {N_N}$ ，可达到降低计算量的目的。

3.3   直升机旋翼参数估计流程

直升机旋翼参数估计具体步骤如下：

步骤1 对直升机旋翼信号进行短时傅里叶变换，得到时频图像 ${\rm{TF}}(t,f)$ 。

步骤2 对时频图中正频率轴数据幅值进行累加计算，并判断累加后数据局部峰值点，对应时频域正频率“闪烁”发生的时间。同样，对负频率轴数据幅值进行累加计算得到时频域中负频率“闪烁”发生的时间。

步骤3 根据步骤2中正负频率“闪烁”发生的时间，判别时频域中正负多普勒“闪烁”是否交替出现。若是，则旋翼叶片数为奇数，否则，旋翼叶片数为偶数。

步骤4 读取某一正频率闪烁发生的时间 ${t_0}$ 及闪烁间隔 $\Delta t$ 。依据式(7)用叶片数量 $N$ 表示出旋翼转速，依据式(10)和式(11)用叶片数量 $N$ 及整数 $k$ 表示出第 $k$ 个叶片初相。

步骤5 确定 $(L,N,k)$ 的取值范围并离散化： $L \in ({L_1}, ·\!·\!· ,{L_r}, ·\!·\!· ,{L_{{N_L}}})$ , $N \in ({N_1}, ·\!·\!· ,{N_p},·\!·\!· ,{N_{{N_N}}})$ , $k \in (1, ·\!·\!· ,{k_q}, ·\!·\!· ,{k_{{N_q}}})$ $\;({k_{{N_q}}} \le {N_{{N_N}}})$ 。利用步骤4中表示出的旋翼转速及初相，依据式(13)和式(15)构建字典矩阵。

步骤6 利用OMP算法寻找叶片数量，叶片长度的最优值，代入式(7)计算出旋翼转速，代入式(10)和式(11)计算出叶片初相。

4.   仿真验证

4.1   仿真结果分析

结合上述模型对直升机旋翼回波信号进行仿真，仿真参数设置如表1所示。

表  1  外辐射源雷达直升机旋翼回波模型仿真参数

Table  1.  Simulation parameters of helicopter rotor echo model for passive radar

信号载频	叶片数	叶片长度	旋转速率	发射站方位角	接收站方位角	发射站仰角	接收站仰角	SNR
658 MHz	3	5 m	200 rpm	33°	76°	23°	23°	–5 dB

下载: 导出CSV 

| 显示表格

图3(a)显示了信号的联合时频域特征，可看出闪烁信号及噪声严重影响直升机旋翼微多普勒特征曲线的检测，使微多普勒特征曲线提取困难。

分别对时频图像中正负频率轴数据幅值进行累加计算，得到时频域中正负多普勒“闪烁”时间，如图3(b)所示，图中正负多普勒“闪烁”等间隔交替出现，则旋翼叶片数为奇数。读取闪烁信号时间间隔为0.05 s，根据式(7)表示出旋翼转速为：

读取某一正频率闪烁信号对应时刻为0.066 s(此处选择了图3(b)中的第1个正频率闪烁信号)，根据式(10)和式(11)表示出第 $k$ 个叶片初相为：

图3(c)为利用OMP方法对 $(L,N,k)$ 的估计结果，得到叶片数为3片，图中给出了其对应的切面图，3叶片长度分别4.99 m, 5.00 m, 4.98 m，均值4.99 m，与理论基本一致，代入式(17)得3叶片初相分别为1.14 rad, 3.24 rad, 5.34 rad，代入式(16)得旋翼转速为200 rpm，与理论值一致，本文方法准确实现了直升机旋翼参数估计。

图  3  本文方法参数估计结果

Figure  3.  Parameter estimation by this article method

下载: 全尺寸图片 幻灯片

图4为利用常规Hough变换，通过微多普勒曲线 $f = {f_{\max }}\sin(2{{π}} {f\!_r}t + {\varphi _0})$ 检测对参数 $({f\!_r},{\varphi _0},L)$ 的估计结果。其中 ${f_{\max }}$ 为最大频移。

图4中给出了参数空间中局部峰值点中心位置。可得到直升机旋翼转速为200 rpm。3叶片最大频移分别为385.6 Hz, 393.9 Hz, 389.8 Hz，平均值为390.0 Hz，由式(18)计算得叶片长度为4.96 m，与理论值基本一致。3叶片初相分别为0.91 rad, 3.16 rad, 5.24 rad，利用式(4)对初相进行修正，得到3叶片初相位为1.86 rad, 4.11 rad, 6.19 rad，存在较大的误差，是由于STFT受不确定原理的限制，时频图像中时频分辨率受限使参数空间中的局部峰值点扩展范围较大，只能大致估计局部峰值点的位置，估计结果精度较低。

图  4  常规Hough变换参数估计结果

Figure  4.  Parameter estimation by traditional Hough transform

下载: 全尺寸图片 幻灯片

4.2   计算复杂度分析

设待处理的时频图像大小为 ${N_t} \times {N_f}$ 像素， ${N_f} \approx {N_t}$ ，利用常规的Hough变换对微多普勒曲线 $f = {f_{\max }}\sin(2{{π}} {f\!_r}t{\rm{ + }}{\varphi _0})$ 进行检测，参数 $({f\!_r},{\varphi _0},L)$ 分别量化为 ${N_{f_{r}}}$ , ${N_{{\varphi _0}}}$ , ${N_L}$ 份。乘法次数可近似表示为： $2{N_{f_{r}}}{N_{{\varphi _0}}}{N_L}{N_t}^{\!2}$ 。

直接使用OMP进行参数 $({f\!_r},{\varphi _0},L)$ 估计时，设迭代次数为K，乘法次数近似表示为： $K{N_{f\!_{r}}}{N_{{\varphi _0}}}$ ${N_L}{N_t}^{\!2}$ 。

本文方法计算量集中在OMP阶段，根据提取的时频域中的闪烁时间，依据式(7)和式(10)，最终转化为对参数 $(L,N,k)$ 的估计，常见直升机的叶片数只有若干个取值，且由3.1节方法可判断出叶片数量的奇偶性， ${N_N}{N_k}$ 远小于 ${N_{f_{r}}}{N_{{\varphi _0}}}$ 。本文方法乘法次数近似表示为： $K{N_N}{N_k}{N_L}{N_t}^{\!2}$ 。

在对直升机旋翼微动参数估计时，一般 ${N_N}{N_k}$ 取值量级为10¹～10²，迭代次数K的取值量级为 ${10^0} {\text ～} {10^1}$ ，当初相 ${\varphi _0}$ 的估计精度为7°时， $2{N_{{\varphi _0}}}$ 取值量级为10²，当转速 ${f_r}$ 的估计精度为10 rpm时， ${N_{f_{r}}}$ 的取值量级为 ${10^1}$ ，在乘法次数上，常规 Hough变换参数估计方法为本文方法 ${10^0}{\text ～} {10^2}$ 倍，当进一步提高 ${f_r},$ ${\varphi _0}$ 的估计精度时，算法之间的计算量差距将进一步变大。本文在相同的配置环境下，利用Matlab仿真平台，常规Hough变换方法运行时长13006 s，而本文方法运行总时长只有145 s。

5.   实验结果

5.1   实验场景

武汉大学电波传播实验室对EC_120B直升机进行了微多普勒效应探究外场实验，EC_120B直升机主旋翼3叶片，叶片长度5 m，额定转速406 rpm，实验中以武汉龟山电视塔数字电视信号为照射源，信号中心频率为658 MHz，带宽8 MHz，接收站位于武汉大学电波传播实验室楼顶，距离发射站7.56 km，实验场景如图5所示。本组实测数据相干积累时间0.8 s，可近似认为目标在这段时间位置不变，直升机旋翼转速为常量。

图  5  实验场景图

Figure  5.  Experimental scene map

下载: 全尺寸图片 幻灯片

5.2   实验结果分析

图6(a)为去除目标主体影响后，对直升机旋翼回波信号进行短时傅里叶变换后的时频图像。可以观察到闪烁信号，但微多普勒特征曲线已观察不到。分别对时频图像中正负频率轴数据幅值进行累加计算，得到时频域中正负多普勒“闪烁”时间，如图6(b)所示，图中正负多普勒“闪烁”等间隔交替出现，则旋翼叶片数为奇数。

图  6  本文方法参数估计结果

Figure  6.  Parameter estimation by this article method

下载: 全尺寸图片 幻灯片

读取闪烁信号时间间隔为26.2 ms，根据式(7)用 $N$ 表示出旋翼转速。读取某一正频率闪烁时间对应时间为0.25 s，利用式(10)和式(11)表示出第 $k$ 个叶片初相，图6(c)为利用OMP方法对 $(L,N,k)$ 的估计结果，得到叶片数量为3片，3叶片长度分别为4.93 m, 5.00 m, 4.66 m，均值4.86 m，存在较小的误差，与仰角，方位角估计不精确有关，由式(7)知旋翼转速均为382 rpm，符合实际情况。

6.   结束语

本文根据外辐射源雷达直升机旋翼微动信号模型，充分利用时频域中闪烁信号特征和微动信号内在特性进行了参数估计。通过时频分析和正交匹配追踪算法，估计出了旋翼转速、叶片长度、叶片数量和初相。同时开展了外场实验。仿真数据和实测数据处理都表明本文方法对外辐射源雷达直升机旋翼参数估计的可行性。