1.   引言

在许多无线遥感探测系统中，微动目标的检测和估计都有着广泛且重要的应用意义。例如，在地震救援、生物医疗等领域，微动目标的检测和估计可以应用于生命信号探测；在防空和监控领域，利用目标微振动引起的微多普勒效应，微动目标的检测和估计可以对无人机机翼进行检测和识别；同时，在机械性能分析方面，可以实现对机械轴承、车床流水线等器械的振动检测，从而判断机械是否正常工作^[1–5]。

关于微动目标参数估计，国内外许多研究机构对此开展了相应的研究。总的来说，傅里叶分析是实现微动目标微摄动参数估计的一种常用的方法^[6–8]。其中，文献[6]采用循环相关系数的方法，利用微动目标时频特征的周期性，实现了微动周期的估计。文献[7]采用时频分布-霍夫变换的方法，估计了目标微动周期和幅度等参数。文献[8]研究了在图像域的微动目标检测方法，并且利用超宽带雷达实现了对目标的定位。

值得指出的是，上述算法对微动目标的参数估计在高斯背景下能够呈现出较好的估计性能。然而，在非高斯背景下，基于傅里叶分析的相关研究方法的效果将发生明显的下降。在非高斯背景下，文献[9]采用了自相关熵的方法，从而估计了微动目标的微多普勒频率，然而，自相关熵的方法引入了非常明显的低频噪声。

为了进一步提高非高斯背景下微动参数的估计性能，本文采用单发多收(Single Input Multiple Output, SIMO)的雷达系统，提出了一种基于互相关熵的估计算法。具体来说，我们先采用了高斯核函数，对每个距离单元，计算每个接收通道的自相关熵矩阵；然后，对每两个通道间的自相关熵矩阵做互相关处理，得到多个互相关熵矩阵；紧接着通过算术融合的方式，对上一步骤得到的多组互相关熵矩阵做融合处理，得到融合后的互相关熵矩阵，该矩阵成像后，就可以得到最终的微动目标成像效果。最后，通过对多通道回波信号比相处理，本算法可以实现对目标的准确定位^[10]。仿真结果证明了该方法的有效性。此外，仿真结果也证实了互相关熵算法能对非高斯背景下的微动目标实现多微动参数估计的同时，提高雷达的微动目标成像质量和输出信噪比。

2.   信号模型

我们考虑一个单发多收的雷达系统，系统模型如图1所示，假设发射天线T位于坐标原点o(0, 0)，接收天线沿x轴等间隔分布，每个天线间的间距为半波长。雷达发射载频为f₀的线性调频连续波(Linear Frequency-Modulated Continuous Wave, LFMCW) s(t), s(t)的时频关系如图2所示，其第n周期发射信号的时域表达式为：

其中， $n = 1,2, ·\!·\!· ,N$, A₀为发射信号幅度， ${\rm j} = \sqrt { - 1}$, k=B/T为调频斜率，其中，B为信号带宽，T为信号时宽，u(t)是矩形函数，表示为：

图  1  SIMO雷达模型

Figure  1.  The SIMO system model

下载: 全尺寸图片 幻灯片

图  2  LFMCW时频关系

Figure  2.  Frequency-time relation of LFMCW

下载: 全尺寸图片 幻灯片

2.1   目标模型

假设一个远场微动目标初始时刻位于 ${P_0}\left( {{x_0},{y_0}} \right)$，由于目标的微动幅度通常来说非常小，因此目标的微动不会引起跨距离单元的问题。从而，本文假设由目标复合微动引起的微小距离变化对每个接收通道而言是一致的。对于两种振动模式的微动目标，微动引起的周期性的距离波动可以建模成双调谐信号

其中，A_B, f_B, ${\varphi _{\rm{B}}}$, A_H, f_H, ${\varphi _{\rm{H}}}$分别表示两种不同振动的幅度、频率和初相位。其中， ${A_{\rm{B}}},{A_{\rm{H}}} \ll$ $c/2B$, c/2B为雷达系统的距离分辨率。那么目标到第i个接收天线的距离可以表示为：

其中，r_i 是第i个接收天线到目标表面的平均距离。对于第i个接收天线，第k周期的目标回波可以表示为：

其中， $i = 1,2, ·\!·\!· ,M,n = 1,2, ·\!·\!· ,N$, K_i为反射后回波幅度衰减系数， ${\tau _i}\left( t \right)$表示第k周期发射信号的回波延时，其中

其中，c为光速， ${\tau _i}$是回波的平均延时， $\tilde \tau \left( t \right)$是由微动引起的随时间变化的回波延时。

2.2   回波模型

考虑噪声时，第i个接收天线的总回波可以表示为：

其中，z(t)为噪声。将回波信号与发射信号混频和低通滤波后，我们得到差拍信号为：

略去式(8)中的平方项，则

代入时延公式(6)，则

其中 $\varPsi = \displaystyle\frac{1}{2}{K_i}A_0^2\exp \left( {{\rm j}2{{π}} \left( {{f_0} + \mu t} \right){\tau _i}} \right)$，且 $\tilde \tau \left( t \right) \ll {\tau _i}$。代入 $\tilde \tau \left( t \right)$的具体展开式，式(10)可以表示为：

根据贝塞尔公式 ${{\rm e}^{{\rm j}x\sin \varphi }} = \sum\nolimits_{m = - \infty }^{ + \infty } {{J_m}} \left( x \right){{\rm e}^{{\rm j}m\varphi }}$，式(11)可以展开为：

由式(12)可知，微动目标的差拍信号中包含有直流分量、微动频率成分f_H和f_B、微动的谐波分量m₁f_H和m₂f_B，以及联合谐波分量m₁f_H+m₂f_B。由公式可看出f₀, A_H和A_B的大小决定了上述频率成分的信号强度。由于目标的微动引起的运动幅度A_H, A_B通常来说非常小，因此设置合适的f₀值可以使谐波分量的强度远弱于目标微动频率的信号成分，从而在对差拍信号去直流后可实现对多个微动频率的估计。

根据贝塞尔函数的特性，当载频f₀值选取合适时，多次谐波分量的强度远弱于目标微动频率的信号成分。我们可以将式(12)中的多次谐波项略去，仅保留直流分量和1次谐波，即m₁, m₂=0或1，此时，式(12)可以化简为：

由于f_B, f_H, ${f_{\rm{B}}} + {f_{\rm{H}}} \ll \mu {\tau _i}$，因此可以忽略微动引起的频率叠加。对式(13)做时域的傅里叶变换后，我们可以得到目标的距离谱，从而计算出第i个接收天线到微动目标的距离r_i，实现目标的测距^[11]。同时，对式(13)做频域的傅里叶变换，可以得到目标的多普勒谱，当f₀选取合适时， ${J_0} \gg {J_1}$，联合微动频率f_B+f_H的强度远远弱于微动频率f_B, f_H，因此经过互相关熵算法处理，除去大部分的噪声干扰成分后，微动频率分量可以从频率谱获得。

3.   互相关熵算法

3.1   互相关熵

互相关熵是由自相关熵的互相关处理得到的，自相关熵的原理如下。

自相关熵：对于一个平稳随机过程x(t)，定义自相关熵为^[12,13]：

其中，E[·]表示随机过程变量x(t)的数学期望，且

表示高斯核函数， $\sigma$是核的尺度参数值。将式(15)代入式(14)，自相关熵可变为：

自相关熵包含着随机变量 $\left\| {x(t) - x(t{\rm{ + }}\tau )} \right\|$所有的偶阶矩。如果要使相关熵具有时不变性质，即 $V\left( {t,t + \tau } \right) = V\left( \tau \right)$，那么输入随机过程的偶阶矩必须是时不变的，在本文中假设此条件是成立的。

对于离散时间平稳随机过程，自相关熵可表示为：

其中，x(n)是离散时间随机过程。自相关熵V[m]可以通过式(17)估计得到：

对V[m]做傅里叶变换可得到相关熵谱：

式(19)表示了随机过程的频域特性，被定义为自相关熵频谱密度(ACSD)。因此通过对自相关熵频谱密度的估计可以实现对微动目标的探测，并且在非高斯背景低信噪比情况下，可以获得比传统方法更高的输出信噪比^[9]。

互相关熵：假设第k个接收天线的自相关熵为V_k, $k = 1,2, ·\!·\!· ,M$，则k通道和l通道的互相关熵可以表示为：

对V_kl[m]做傅里叶变换可得到互相关熵谱：

根据互相关处理的原理，V_k和V_l中相关的目标信息被放大，而非相关的随机噪声和杂波成分被压制^[14]。因此，互相关熵可以在自相关熵的基础上进一步地去掉随机噪声和杂波成分，使微动目标实现更清晰的成像效果。同时，通过对多通道互相关熵的算术融合处理，可以更进一步地提高微动目标的成像效果和估计性能。

3.2   算法流程图

SIMO LFMCW雷达体制下的微动目标探测，在实现微动频率探测的同时还要完成目标距离信息的探测。互相关熵算法可以实现对上述微动参数的探测和估计。具体算法流程如下：

(1) 首先对雷达接收机M路通道所得的N×L维(N为每个通道的总的回波数，L为每个回波信号的采样点数)回波矩阵A₁, A₂, ···, A_M分别在快时间上(按行)进行加窗处理且做FFT(傅里叶变换)，使每个周期的回波信号压缩成sinc时域脉冲信号，sinc带有了目标的距离信息；

(2) 将N个扫频周期上进行的上述操作得到的N个sinc信号组合在一起，得到M路通道微动目标的距离-脉冲域矩阵Z₁, Z₂, ···, Z_M；对Z₁, Z₂, ···, Z_M分别在慢时间上(按列)作MTI滤波，得到去除了静止背景等零频杂波的距离-脉冲域矩阵D₁, D₂,···, D_M；

(3) 对Z₁, Z₂, ···, Z_M分别在慢时间上(按列)作MTI滤波，得到去除了静止背景等零频杂波的距离-脉冲域矩阵D₁, D₂, ···, D_M；

(4) 对距离-脉冲域矩阵D₁, D₂, ···, D_M分别在慢时间上(按列)计算自相关熵，得到自相关熵矩阵V₁, V₂, ···, V_M；

(5) 对自相关熵矩阵V₁, V₂, ···, V_M任选两个作为一组做互相关处理，得到 $C_M^2$组互相关熵 ${{V}_{il}}\, (i,l = 1,2,·\!·\!·M$， ${\text{且}}i \ne l)$；为了获得更好的成像效果以及更高的输出信噪比，我们对最后获得的多组互相关熵进行算术融合得到最后的互相关熵矩阵V，其中：

(6) 互相关熵矩阵V经过慢时间维上傅里叶变换处理得到互相关熵谱 ${P_{V}}\left( \omega \right)$。其对应的成像平面能清晰明确地显示出微动目标的距离、微动频率等参数；对 ${P_{V}}\left( \omega \right)$进行高通滤波，利用先验信息去除其它频段噪声干扰，得到多普勒-距离平面，可以进一步提高目标的成像质量，提升估计性能。

算法流程图如图3所示。

图  3  算法流程图

Figure  3.  The flowchart of the mutual correntropy algorithm

下载: 全尺寸图片 幻灯片

4.   多通道目标定位

单脉冲比相方法(Phase-Comparison Monopulse, PCM)是通过比较两天线接收信号的相位信息来确定目标角度的方法^[15]。比相雷达在每个坐标轴(方位和俯仰)方向采用了至少2个接收天线，如图4所示。

图  4  PCM测角示意图

Figure  4.  The schematic of PCM algorithm

下载: 全尺寸图片 幻灯片

由相关的公式推导^[15]，相邻天线间的相位差表示为：

其中， $\varphi$为目标所在方向，R₁, R₂分别为目标到相邻天线间的距离，d为两天线中心的间距且 $d \ll R$, $\lambda$为发射信号波长。如果 $\varphi {\rm{ = }}0$，则目标在主轴上。已知d和 $\lambda$，可以利用相位差 $\phi$来确定目标的方向 $\varphi$。

式(13)中，由于目标微动幅度最大振动幅度 ${A_{\rm{B}}} + {A_{\rm{H}}} \ll \displaystyle\frac{c}{{2B}}$，可以假设微动对每个通道而言，引起的距离变化是一致的。因此，单脉冲比相定位时，微动引起的回波相位变化是一致的，不影响相位差 $\phi$的计算。为了简化定位过程的分析，令

则第n周期发射信号对应的第i通道的差拍信号表示为：

其中， ${\rm{ }}{f_{{b_i}}} = \displaystyle\frac{{2{r_i}}}{c}\mu$, ${\varphi _{{b_i}}} = \displaystyle\frac{{4{{π}} {r_i}}}{c}{f_0}$。因此，差频信号可近似为单频信号。采样后，差拍信号为：

其中， ${f_{{b_i}}}$为信号载频， ${\varphi _{{b_i}}}$为信号初相，T_s为采样间隔( ${T_{\rm{s}}} = 1/{f_{\rm{s}}}$), L为采样点数( $T = {T_{\rm{s}}}L$)。对 ${\tilde y_{i,n}}(l)$进行DFT变换得到差拍信号频谱，表示为：

由式(27)可以看出 ${{{\tilde{Y}}}_{i,n}}(w)$的幅度最大谱线位置 ${w_0} = [{f_{{b_i}}}T\;]$。此时 ${{{\tilde{Y}}}_{i,n}}(w)$的相位为：

设相邻的接收天线i与接收天线k的接收信号分别为 ${{\tilde{y}}_{i,n}}\left( l \right)$与 ${{\tilde{y}}_{k,n}}\left( l \right)$。分别对其DFT变换可得到 ${{\tilde{Y}}_{i,n}}\left( w \right)$与 ${{\tilde{Y}}_{k,n}}\left( w \right)$。由两天线间距 $\,d \ \ \ll \ R$，则 $[{f_b}_{_i}T\;] = [{f_b}_{_k}T\;]$。所以两天线接收信号的相位差 $\phi$为：

其中， ${{\psi }_{i}}({{w}_{0}})=\arctan \left( \frac{\operatorname{Im}\left( {{{\tilde{V}}}_{i,n}}\left( {{w}_{0}} \right) \right)}{\operatorname{Re}\left( {{{\tilde{V}}}_{i,n}}\left( {{w}_{0}} \right) \right)} \right),$ ${{\psi }_{k}}({{w}_{0}})=\arctan \left( \frac{\operatorname{Im}\left( {{{\tilde{V}}}_{k,n}}\left( {{w}_{0}} \right) \right)}{\operatorname{Re}\left( {{{\tilde{V}}}_{k,n}}\left( {{w}_{0}} \right) \right)} \right),$，则

根据式(23)，可计算目标的方位角：

结合测得的微动目标的距离R，可以实现对微动目标的定位，其坐标计算公式为：

5.   仿真及结果分析

为了验证互相关熵算法及定位算法的有效性，我们采用1发4收的天线阵列，并以发射天线坐标为原点，对远场复合微动目标进行参数估计和定位。假设微动目标以双调谐方式振动，振动频率分别为f₁和f₂，振动幅度分别为A₁和A₂，微动目标坐标为(0, 6) m，具体参数如表1所示。

表  1  仿真参数

Table  1.  Simulation parameters

参数	探测参数	参数值
	载波频率f0	1 GHz
	信号带宽B	500 M
系统参数	信号时宽T	3.31 ms
	LFM周期数N	512
	核尺度 $\sigma$	0.2
目标参数	距离R	6 m
	振动频率f1	4 Hz
	振动频率f2	8 Hz
	振动幅度A1	4 mm
	振动幅度A2	4 mm
	坐标(x, y)	(0, 6) m
	信噪比SNRin	–10 dB
场景噪声参数	噪声类型	韦布尔噪声
	韦布尔形状因子P	2

下载: 导出CSV 

| 显示表格

在本文中，为了验证互相关熵算法的有效性，非高斯噪声采用韦布尔噪声进行仿真。韦布尔噪声的概率密度函数为^[16]：

其中，x≥0, P>0, Q>0，且P为韦布尔噪声分布的形状因子，Q为尺度因子。韦布尔分布的数学期望 ${\mu _{\rm{w}}}$和方差 $\sigma _{\rm{w}}^2$分别为：

其中， $\Gamma \left[ \cdot \right]$为伽马函数，方差 $\sigma _{\rm{w}}^{\rm{2}}$即为韦布尔噪声的噪声功率。

定义输入信噪比(input SNR)为：

当形状因子P和输入信噪比SNR_in给定时，由于发射信号的幅度A₀已知，根据上述理论分析，逆向推导，即可计算韦布尔分布的尺度因子Q，从而生成韦布尔噪声。仿真所采用的韦布尔噪声参数如表1所示。

图5是传统傅里叶分析下微动目标参数估计结果。图5(a)是目标的距离像，图中表明目标距离为6.024 m，与实际目标距离相差不大。图5(b)是目标的2维成像，从图中可以看出目标的微动频率1，但是目标的成像背景强干扰成分很多，微动频率2被淹没，大大降低了参数估计性能。

图  5  传统傅里叶分析结果

Figure  5.  Conventional Fourier analysis results

下载: 全尺寸图片 幻灯片

为了对成像效果对比分析，图6给出了基于自相关熵融合前后的分析结果。图6(a)是自相关熵分析方法融合前的目标成像，目标参数如图中标明，与实际目标参数相差不大，相比传统傅里叶分析效果，自相关熵的成像效果中，大部分的随机噪声被消除，但低频分量依旧存在。图6(b)是融合后的目标成像，从图中可以看出算术融合后，低频分量被轻微抑制。

图  6  自相关熵分析结果

Figure  6.  The auto-correntropy processing results

下载: 全尺寸图片 幻灯片

图7给出了互相关熵算法融合前后的对比分析结果。图7(a)是互相关熵分析方法融合前的目标成像，目标参数如图中标明，与实际目标参数相差不大，相比图5和图6所呈现的成像效果，图7(a)中，更多的干扰成分被消除。图7(b)是互相关熵算法融合后的成像效果，算术融合后，从图中可以看出成像背景的干扰成分基本被消除，提高了成像效果，提升了参数估计性能，进一步证明了该算法的有效性。

图  7  互相关熵分析结果

Figure  7.  The mutual correntropy processing results

下载: 全尺寸图片 幻灯片

为了证明非高斯背景下互相关熵算法的有效性，我们计算不同方法成像时对应的输出信噪比(output SNR)从而量化不同算法的成像性能。计算方法为：首先我们计算相应成像矩阵中所有元素的幅值之和从而获得输出功率P_S，然后减去目标所在位置处的部分幅值之和P_T，最后得到成像时的噪声功率P_N，即

因此，输出信噪比可以定义为

图8是非高斯背景下不同方法对应的输入输出信噪比结果。从图中可以看出，基于互相关熵算法的输出性能明显超过自相关熵算法和传统傅里叶算法的输出性能。再者，比较融合前后的输出信噪比曲线，可以看出，算术融合可以进一步提高目标的成像性能。

图  8  非高斯背景下输入输出信噪比关系

Figure  8.  The output SNR versus different input SNR

下载: 全尺寸图片 幻灯片

图9是微动目标的定位结果，蓝点代表目标的实际位置，坐标为(0, 6) m，红点是检测到的目标位置，坐标为(0.05, 6.02) m，可以看出，定位出来的坐标误差非常小，证明了定位算法的有效性。

图  9  微动目标定位结果

Figure  9.  The location result of PCM

下载: 全尺寸图片 幻灯片

6.   结论

本文研究了非高斯背景下微动目标的参数估计和目标定位的问题。为了得到更好的估计性能，本文利用单发多收的雷达系统，提出了基于互相关熵的估计算法估计目标的微摄动参数。总的来说，互相关算法去除了传统傅里叶分析方法下的随机噪声干扰，在一定程度上减轻了自相关熵分析时所引入的低频分量干扰。该方法能对非高斯背景下的微动目标实现多微动参数估计，在提高雷达成像质量和输出信噪比的同时，不会造成微动参数的丢失。仿真结果证明了该方法的有效性。同时，通过对多通道回波信号比相单脉冲定位处理，本文实现了对微动目标的准确定位，仿真结果证明了该结论。