Micro-Motion Feature Extraction Method Based on Edge Detection in Terahertz Band

振动模型是一种典型的微动模型，弹簧振子的振动、甚至发动机的机械振动都可以用振动模型来近似。振动模型具有周期性。图 1示出了雷达和振动质点的几何关系。

雷达位于坐标系$(X,Y,Z)$的原点，振动的目标可以看作一个质点，质点以点$A$为振动中心以固定的频率${{\omega }_{v}}$作振动，与雷达的距离为${{R}_{0}}$，$A$点在坐标系中方位角和俯仰角分别为$\alpha $和$\beta $，再以质点振动中心$A$为原点建立坐标系$({X}',{Y}',{Z}')$，在该坐标系下，当质点振动至$B$点时，其位移为${{C}_{t}}$，此时质点的方位角和俯仰角分别为${{\alpha }_{B}}$和${{\beta }_{B}}$。此时，雷达的接收信号可以用下式表示：

式中，$\rho $为目标点的反射率；${{f}_{c}}$为信号载频；$\phi (t)$为相位函数：

对时间求导，得到多普勒频率：

图 2为具有4个振动部分和1个平动部分的复杂目标回波信号仿真。从图中可以看出，振动目标的回波信号呈典型的正弦特征。

此外，典型的雷达微动特征还有旋转、翻滚等等，它们的特征模型均可由振动模型推导。

由于微动目标回波信号是典型的非平稳信号，通常情况下采用时频分析的方法获得信号的时频谱，通过时频谱上对频谱曲线进行提取获取目标微动特征的相关参数。常见的时频曲线提取方法有脊线提取与质心曲线提取两种。

对于时频谱来说，在每个时间点的列上都会有一个最大值，脊线就是把每个时间点对应的列上的最大值作为该时刻的频率值，然后将这些点相连得到一个以时间为变量的函数，以该函数近似即时频率曲线。时频谱脊线可以由下式得到：

如果对每个时刻都选取整个频率范围内求加权平均值，每个时刻计算的平均值所对应的频率值即为该时刻的频率，这就是时频谱质心曲线提取方法。时频谱质心曲线可由下式得到：

然而，这两种传统算法只能提取单一曲线，为了能同时提取出多条时频谱曲线，本文提出了一种基于边缘检测的时频谱曲线提取方法。其算法流程如图 3所示。

1) 对回波的时频谱二值化处理。

具体方法是对时频谱图像$\left| {\rm STFT}(t,\omega ) \right|$使用阈值分割的方法：

其中$x$与时频谱$\left| {\rm STFT}(t,\omega ) \right|$中的频率值$\omega $对应，$y$与时间$t$对应，${{T}_{y}}$为阈值。

图 4示出了图 2中无噪声情况下两种回波时频谱的直方图，该直方图显示了场景与目标回波强度之间的对比关系，便于对阈值进行选取。

2) 对时频谱二值化的结果$g(x,y)$进行边缘检测。

对于时频谱来说，每一列有其特殊的意义，即每一列代表一个时刻该信号的频谱，因此对时频谱按列进行边缘检测才有意义，所以文本选用列向量${{\left[ -1,1 \right]}^{\prime }}$作为边缘检测算子对步骤1) 中二值化的时频谱进行边缘检测，其结果如图 5所示。

3) 将每个频谱条带的边缘合并至中心位置，得到时频谱曲线。

由于使用的边缘检测算子是按列卷积，所以每列的点数均为偶数，且两两一组，因此大大简化了对频谱中心的提取，将每组的两个边缘点对应的频率${{\omega }_{1}},{{\omega }_{2}}$找出，再求$\left| {\rm STFT}(t,\omega ) \right|$在该时刻两个频率的加权平均值即为即时频率。为了简化算法可以直接用${{\omega }_{1}},{{\omega }_{2}}$的平均值代替，得到时频谱曲线点的矩阵h(x, y)。图 6显示出了将边缘合并为曲线的结果。

4) 删除部分曲线点，减小曲线点的结果中噪声的影响。

从图 6可以看到，曲线点基本上都集中在时频谱条带的范围内。但是在有噪声的情况下，除了相对集中的时频谱曲线点外，还有一些噪声带来的点会分散在其他频率，因此需要去除这些相对分散的点，以减小噪声的影响。可以用图 6的曲线点或者图 4时频谱二值化后的结果进行统计。

5) 对步骤4) 的结果$h(x,y)$进行Hough直线检测。

Hough直线检测得到的参数如图 7所示，该直线反应了主体运动的信息。由于微动是围绕主体运动进行的，其时频谱也围绕主体运动，所以在进行下一步时，要先对$h(x,y)$进行校正，将其每列都以直线点对齐，以便对曲线参数的检测。

6) 对$h(x,y)$进行曲线检测。

时频谱曲线含有频率、幅度和相位参数，所以设置频率、幅度和初相的参数空间$(f,A,\ \theta )$，参数空间中的每一个点$({{f}_{i}},{{A}_{j}},\ {{\theta }_{k}})$代入运动曲线模型：

再将$h(x,y)$中所有的点进行计算，若符合该参数空间点建立的上述模型，则对该参数空间点进行累加，得到整个参数空间的累加值$\sum\limits_{i,j,k}{g({{f}_{i}},{{A}_{j}},\ {{\theta }_{k}})}$，再找出累加器中的局部最大点对应的坐标位置，其对应的3个参数，即反应了每条曲线的频率、幅度和初相的参数。图 7给出了参数空间的累加情况及其对应的参数。

再根据直线和曲线的参数将曲线转化为关于时间的函数，得到多条时频谱曲线，如图 8所示。

上图提取的即为图 2的运动回波时频谱的曲线，可以看到，原时频谱中的5个条带对应的曲线都被很好地提取出来。每条曲线都对应了目标的一个微动散射中心，对这些曲线进行分析就能得到具体的微动参数。

对于这种有多个散射中心的情况，由于其回波的时频谱中每个中心都应该对应一条曲线，而传统的脊线和质心曲线方法都只能从一个时频谱图像中提取出一条曲线。因此若使用传统的曲线提取方法，则无法得各微动部分的时频谱曲线，其脊线提取和质心曲线提取的结果如图 9所示。

可以看出，基于边缘检测提取时频谱曲线的方法可以同时对多条曲线进行提取，克服了脊线提取和质心曲线提取这两种方法的缺点。

实验使用的太赫兹雷达系统实物图如图 10所示。

雷达系统发射信号的中心频率在0.3 THz附近，信号带宽10.08 GHz，调频斜率3.36×10¹³ Hz/s，发射脉冲信号的宽度300 μs，脉冲周期为1 ms。本文选取了单摆小球去近似人体模型进行实验。实验场景及示意图如图 11所示。

可以看出，目标有横杆和小球两个散射中心，场景中各参数在示意图中示出，其中横杆的半径约为${{r}_{g}}$=1.2 cm，小球的半径为${{r}_{b}}$= 0.8 cm，摆长为$l$= 57 cm，摆角为$\theta $= 9°，小球在最低点时与雷达的距离为$R$= 6 m，这里以小球在最低点时作为参考点，参考距离${{R}_{{\rm ref}}}$= 6 m。

由单摆周期公式计算可得，单摆的周期：

横杆高度为$l+{{r}_{g}}$= 58.2 cm，与相对于参考点与雷达的距离差：

小球在最高点时，与雷达的水平距离${{R}_{x}}\text{=R}-l\sin \theta \approx $5.910 8 m，垂直方向距离为${{R}_{y}}\text{=}l(1-\cos \theta )\approx $0.007 0 m，所以在最高点小球与雷达距离为${{R}_{H}}\text{=}\sqrt{R_{x}^{2}+R_{y}^{2}}=$5.910 8 m，此时，小球相对于参考点与雷达的距离差为：

图 12示出了单个小球实验的回波时频谱。

得到其时频谱曲线后，即可算出目标的即时距离及其他微动参数。结果如图 13所示。

直接提取出的参数与模型计算的参数对比于表 1和表 2中。

从上面的数据及分析可以看到，无论是直接提取的参数还是根据提取值进行进一步理论推导获得的计算值，都与初始设置的模型参数基本一致。该实验在验证了该方法的可行性的同时也说明了该方法在太赫兹雷达下对目标微动参数提取具有较高的准确性。

本文采用了一种基于图像处理的微动特征边缘检测技术，该方法可结合微动特征信号模型进行微动特征曲线及参数提取。实测数据表明，该方法具有较好的性能，能有效提取多散射中心微动目标的特征参数，为进一步开展微动信号处理及目标识别提供了有效的技术途径。