太赫茲頻段微動特征邊緣檢測及提取方法

田 坤，李 晉

(電子科技大學電子工程學院 成都 611731)

太赫茲頻段電磁波是指頻率在0.1～10 THz，即波長為30～3 mm范圍內的電磁輻射。太赫茲雷達大帶寬、高頻率等技術特點可以有效提高雷達系統對微動目標的探測精度[1]。此外，與紅外線和毫米波等相比，太赫茲波在散射、傳播、吸收、反射和穿透等方面都具有明顯的優勢。

從2000年起，歐美國家相關研究機構就相繼開展了關于微動目標雷達探測的研究工作。文獻[2]率先提出了微多普勒的概念，將目標回波內雷達頻率調制稱為微多普勒效應，他對直升機的微動特征進行估計得到旋翼的幾何結構和運動參數，從而實現了直升機的動態識別[3]。2002年，文獻[4]運用全相參X波段連續波雷達獲取人體行走微動數據，用于步態識別研究。2004年，文獻[3]研究了人體走動中胳膊和腿的擺動帶來的微多普勒調制，運用一臺小型X 波段雷達探測人體走動并運用時頻分析手段提取微多普勒特征。2007年，文獻[5-6]通過雷達測量行人獲取大量實測數據，結合小波分析與時頻分析算法提取人體目標微多普勒特征，并利用微多普勒特征精確估計出微動參數。文獻[7]針對傳統時頻分析方法對非平穩信號處理存在的缺陷，引入Hilbert-Huang變換分離噪聲和微動信號，有效地改善微動模式的識別效果。2014年，文獻[8]利用逆radon變換對時頻譜中的微動特征參數進行提取。文獻[9]利用微動特征提取對無人機旋翼目標展開特征識別。2015年，文獻[10]利用微動信息開展空間目標的高分辨率三維成像技術研究。文獻[11]利用微動特征開展了微動目標檢測技術研究，上述這些檢測和成像算法的性能好壞均取決于微動特征參數估計的準確與否。

然而，上述微動目標的探測與特征提取均僅建立在傳統微動信號處理層面，導致其在多散射中心情況下的微動特征信號分離及參數提取性能較差。本文擬在微動特征模型的基礎上，結合信號模型與圖像處理的方法，對微動特征進行提取，進一步改善多散射中心情況下的微動特征參數的提取性能。

1 典型微動特征模型及參數

振動模型是一種典型的微動模型，彈簧振子的振動、甚至發動機的機械振動都可以用振動模型來近似。振動模型具有周期性。圖1示出了雷達和振動質點的幾何關系。

圖1 雷達和振動質點的幾何關系

雷達位于坐標系(X,Y,Z)的原點，振動的目標可以看作一個質點，質點以點A為振動中心以固定的頻率ωv作振動，與雷達的距離為R0，A點在坐標系中方位角和俯仰角分別為α和β，再以質點振動中心A為原點建立坐標系(X′,Y′,Z′)，在該坐標系下，當質點振動至B點時，其位移為Ct，此時質點的方位角和俯仰角分別為αB和βB。此時，雷達的接收信號可以用下式表示：

式中，ρ為目標點的反射率；fc為信號載頻；φ(t)為相位函數：

對時間求導，得到多普勒頻率：

圖2為具有4個振動部分和1個平動部分的復雜目標回波信號仿真。從圖中可以看出，振動目標的回波信號呈典型的正弦特征。

此外，典型的雷達微動特征還有旋轉、翻滾等等，它們的特征模型均可由振動模型推導。

圖2 復雜運動目標回波信號仿真

2 基于邊緣檢測的微動參數提取

由于微動目標回波信號是典型的非平穩信號，通常情況下采用時頻分析的方法獲得信號的時頻譜，通過時頻譜上對頻譜曲線進行提取獲取目標微動特征的相關參數。常見的時頻曲線提取方法有脊線提取與質心曲線提取兩種。

對于時頻譜來說，在每個時間點的列上都會有一個最大值，脊線就是把每個時間點對應的列上的最大值作為該時刻的頻率值，然后將這些點相連得到一個以時間為變量的函數，以該函數近似即時頻率曲線。時頻譜脊線可以由下式得到：

如果對每個時刻都選取整個頻率范圍內求加權平均值，每個時刻計算的平均值所對應的頻率值即為該時刻的頻率，這就是時頻譜質心曲線提取方法。時頻譜質心曲線可由下式得到：

然而，這兩種傳統算法只能提取單一曲線，為了能同時提取出多條時頻譜曲線，本文提出了一種基于邊緣檢測的時頻譜曲線提取方法。其算法流程如圖3所示。

圖3 基于邊緣檢測的時頻譜曲線提取算法流程

1)對回波的時頻譜二值化處理。

其中x與時頻譜中的頻率值ω對應，y與時間t對應，Ty為閾值。

圖4示出了圖2中無噪聲情況下兩種回波時頻譜的直方圖，該直方圖顯示了場景與目標回波強度之間的對比關系，便于對閾值進行選取。

圖4 回波頻譜的直方圖及其二值化結果

2)對時頻譜二值化的結果g(x,y)進行邊緣檢測。

對于時頻譜來說，每一列有其特殊的意義，即每一列代表一個時刻該信號的頻譜，因此對時頻譜按列進行邊緣檢測才有意義，所以文本選用列向量[-1,1]′作為邊緣檢測算子對步驟1)中二值化的時頻譜進行邊緣檢測，其結果如圖5所示。

圖5 回波頻譜的邊緣

3)將每個頻譜條帶的邊緣合并至中心位置，得到時頻譜曲線。

由于使用的邊緣檢測算子是按列卷積，所以每列的點數均為偶數，且兩兩一組，因此大大簡化了對頻譜中心的提取，將每組的兩個邊緣點對應的頻率ω1,ω2找出，再求在該時刻兩個頻率的加權平均值即為即時頻率。為了簡化算法可以直接用ω1,ω2的平均值代替，得到時頻譜曲線點的矩陣h(x,y)。圖6顯示出了將邊緣合并為曲線的結果。

4)刪除部分曲線點，減小曲線點的結果中噪聲的影響。

從圖6可以看到，曲線點基本上都集中在時頻譜條帶的范圍內。但是在有噪聲的情況下，除了相對集中的時頻譜曲線點外，還有一些噪聲帶來的點會分散在其他頻率，因此需要去除這些相對分散的點，以減小噪聲的影響。可以用圖6的曲線點或者圖4時頻譜二值化后的結果進行統計。

圖6 回波頻譜曲線點

5)對步驟4)的結果h(x,y)進行Hough直線檢測。

Hough直線檢測得到的參數如圖7所示，該直線反應了主體運動的信息。由于微動是圍繞主體運動進行的，其時頻譜也圍繞主體運動，所以在進行下一步時，要先對h(x,y)進行校正，將其每列都以直線點對齊，以便對曲線參數的檢測。

6)對h(x,y)進行曲線檢測。

時頻譜曲線含有頻率、幅度和相位參數，所以設置頻率、幅度和初相的參數空間(f,A,θ)，參數空間中的每一個點(fi,Aj,θk)代入運動曲線模型：

再將h(x,y)中所有的點進行計算，若符合該參數空間點建立的上述模型，則對該參數空間點進行累加，得到整個參數空間的累加值再找出累加器中的局部最大點對應的坐標位置，其對應的3個參數，即反應了每條曲線的頻率、幅度和初相的參數。圖7給出了參數空間的累加情況及其對應的參數。

再根據直線和曲線的參數將曲線轉化為關于時間的函數，得到多條時頻譜曲線，如圖8所示。

上圖提取的即為圖2的運動回波時頻譜的曲線，可以看到，原時頻譜中的5個條帶對應的曲線都被很好地提取出來。每條曲線都對應了目標的一個微動散射中心，對這些曲線進行分析就能得到具體的微動參數。

圖7 直線和曲線參數檢測結果

圖8 運動回波時頻譜曲線

對于這種有多個散射中心的情況，由于其回波的時頻譜中每個中心都應該對應一條曲線，而傳統的脊線和質心曲線方法都只能從一個時頻譜圖像中提取出一條曲線。因此若使用傳統的曲線提取方法，則無法得各微動部分的時頻譜曲線，其脊線提取和質心曲線提取的結果如圖9所示。

圖9 目標回波時頻譜的脊線和質心曲線提取

可以看出，基于邊緣檢測提取時頻譜曲線的方法可以同時對多條曲線進行提取，克服了脊線提取和質心曲線提取這兩種方法的缺點。

3 實驗及結果分析

實驗使用的太赫茲雷達系統實物圖如圖10所示。

圖10 0.3 THz雷達系統實物圖

雷達系統發射信號的中心頻率在0.3 THz附近，信號帶寬10.08 GHz，調頻斜率3.36×1013Hz/s，發射脈沖信號的寬度300 μs，脈沖周期為1 ms。本文選取了單擺小球去近似人體模型進行實驗。實驗場景及示意圖如圖11所示。

圖11 實驗場景及示意圖

可以看出，目標有橫桿和小球兩個散射中心，場景中各參數在示意圖中示出，其中橫桿的半徑約為rg=1.2 cm，小球的半徑為rb= 0.8 cm，擺長為l=57 cm，擺角為θ= 9o，小球在最低點時與雷達的距離為R= 6 m，這里以小球在最低點時作為參考點，參考距離Rref= 6 m。

由單擺周期公式計算可得，單擺的周期：

橫桿高度為l+rg= 58.2 cm，與相對于參考點與雷達的距離差：

小球在最高點時，與雷達的水平距離Rx=R-lsinθ≈5.910 8 m，垂直方向距離為Ry=l(1-cosθ)≈0.007 0 m，所以在最高點小球與雷達距離為此時，小球相對于參考點與雷達的距離差為：

圖12示出了單個小球實驗的回波時頻譜。

得到其時頻譜曲線后，即可算出目標的即時距離及其他微動參數。結果如圖13所示。

直接提取出的參數與模型計算的參數對比于表1和表2中。

圖12 回波時頻譜及二值化

圖13 提取出的曲線

表1 目標微動提取參數對比

表2 計算獲得的參數對比

圖14 兩個小球實驗場景及曲線提取結果

表3 兩小球參數設置與提取結果對比

從上面的數據及分析可以看到，無論是直接提取的參數還是根據提取值進行進一步理論推導獲得的計算值，都與初始設置的模型參數基本一致。該實驗在驗證了該方法的可行性的同時也說明了該方法在太赫茲雷達下對目標微動參數提取具有較高的準確性。

4 結 束 語

本文采用了一種基于圖像處理的微動特征邊緣檢測技術，該方法可結合微動特征信號模型進行微動特征曲線及參數提取。實測數據表明，該方法具有較好的性能，能有效提取多散射中心微動目標的特征參數，為進一步開展微動信號處理及目標識別提供了有效的技術途徑。

[1]LI J, PI Y. Micro-doppler signature feature analysis in terahertz band[J]. Journal of Infrared, Millimeter, and Terahertz Waves, 2010, 31(3)： 319-328.

[2]CHEN V C. Analysis of radar micro-doppler signature with time-frequency transforms[C]//Proceedings of the 10th IEEE Workshop on Statistical Signal and Array Processing.Pocono Manor, USA： IEEE, 2000.

[3]CHEN V C, LI F Y. Micro-doppler effect in radar-phenomenon, model and simulation study[J]. IEEE Trans on Aerospace and Electronic Systems, 2006, 43(1)：2-21.

[4]GEISHEIMER J L, GRENEKER E, MARSHALL W S. A high-resolution Doppler model of human gait[C]//Proceedings of SPIE on Radar Technology. [S.l.]： SPIE,2002： 8-18.

[5]THAYAPARAN T, ABROL S, RISEBOROUGH E.Micro-Doppler radar signatures for intelligent target recognition[R]. Ottawa, Canada： Technical memorandum DRDC Ottawa, 2004.

[6]THAYAPARAN T, ABROL S, RISEBOROUGH E, et al.Analysis of radar micro-Doppler signatures from experimental helicopter and human data[J]. IET Radar Sonar Navigator, 2007, 1(4)： 289-299.

[7]LAI C P, RUAN Q, NARAYANAN R M. Hilbert-Huang transform (HHT)processing of through-wall noise radar data for human activity characterization[C]//IEEE Workshop on Signal Processing Applications for Public Security and Forensics. [S.l.]： IEEE, 2007： 1-6.

[8]QI Y, BIN D, HONGQIANG W, et al. Doppler aliasing free micro-motion parameter estimation algorithm based on the spliced time-frequency image and inverse Radon transform[C]//International Conference on Information and Communications Technologies. Nanjing： [s.n.], 2014.

[9]TAHMOUSH D. Detection of small UAV helicopters using micro-Doppler[J]. Proceedings of SPIE-The International Society for Optical Engineering, 2014, 9077(18)： 5450-5453.

[10]BAI X, ZHOU F, BAO Z. High-resolution threedimensional imaging of space targets in micro motion[J].IEEE J Sel Topics Appl Earth Observ Remote Sens, 2015,8(7)： 3428-3440.

[11]LI J, PI Y. Target detection for Terahertz radar networks based on micro-Doppler signatures[J]. International Journal of Sensor Networks, 2015, 17(2)： 115-121.