一種基于雷達微多普勒效應(yīng)蒼蠅識別的方法

沈菲菲，李立鋼，賈配洋，宋振健

（1.中國科學(xué)院國家空間科學(xué)中心，中國科學(xué)院復(fù)雜航天系統(tǒng)電子信息技術(shù)重點實驗室，北京100190；2.中國科學(xué)院大學(xué) 北京100049）

一種基于雷達微多普勒效應(yīng)蒼蠅識別的方法

沈菲菲1，2，李立鋼1，賈配洋1，2，宋振健1，2

（1.中國科學(xué)院國家空間科學(xué)中心，中國科學(xué)院復(fù)雜航天系統(tǒng)電子信息技術(shù)重點實驗室，北京100190；2.中國科學(xué)院大學(xué) 北京100049）

為了有效地識別蒼蠅，提出了利用雷達的微多普勒效應(yīng)來識別的有效方法。首先，根據(jù)微多普勒效應(yīng)的基本概念，建立了蒼蠅飛行復(fù)合運動的雷達回波模型，推導(dǎo)了雷達回波信號和理論的多普勒頻率，提出了綜合運用時頻分析（偽Winger-Ville變換）、譜圖峰值估計、參數(shù)擬合、平動補償?shù)姆椒ǎ罱K提取出蒼蠅在飛行過程中的運動參數(shù)，包括：振翅頻率、振動幅度、運動速度、加速度。

微多普勒；復(fù)合運動；偽Winger-Ville變換；參數(shù)擬合；譜圖峰值估計；平動補償

隨著社會的發(fā)展，生活水平的提高，人們對食品安全越來越關(guān)注。而由媒介生物傳播的疾病史整個食品加工過程中非常值得關(guān)注的一個問題，在許多食品加工場所中，蒼蠅是最常見的一種害蟲，在食品的表面和加工操作表面停留，通過其直接接觸、嘔吐物以及排泄物傳播多種致病菌和寄生蟲，導(dǎo)致疾病的發(fā)生。所以食品加工和食品服務(wù)場所必須適時采取必要的放置措施[1]。

目前，防治蟲蠅的有效手段主要有以下幾種：物理防治技術(shù)：應(yīng)用各種物理因子如光譜、點、顏色、溫度、濕度、聲音等機械設(shè)備來防治[2]以及無公害粘蟲膠等物理手段。但是這些方法均不太適合于場所大的食品加工廠等場所。

雷達成像中，一般將目標(biāo)或者目標(biāo)部件除質(zhì)心平動之外的振動、旋轉(zhuǎn)和加速運動等微小運動統(tǒng)稱為微動[3]，而由目標(biāo)微動產(chǎn)生的對雷達回波信號的展寬頻率調(diào)制稱為微多普勒效應(yīng)（micro-Doppler effect）[4]。這種微動對雷達回波的調(diào)制稱為微多普勒現(xiàn)象[5-6]。微多普勒現(xiàn)象可以反映目標(biāo)的精細運動特征[3-7]，它在目標(biāo)的探測、分類、識別領(lǐng)域都具有很大的應(yīng)用潛力[8]。近年來，許多有效的微多普勒信號提取技術(shù)相繼被提出，例如時頻分析技術(shù)[9-13]、逆Radon變換[14]方法等。

文中提出了一種基于雷達微多普勒效應(yīng)識別目標(biāo)蒼蠅的方法。不僅無色無味沒有公害，而且成本低易于實現(xiàn)。

1 雷達識別蒼蠅的微動模型及數(shù)學(xué)原理

以雷達Q建立參考坐標(biāo)系Q（U，V，W），由于蒼蠅相對于雷達距離來說很小，可以假設(shè)為點目標(biāo)，其振翅中心記作O，以O(shè)為原點建立一個與Q（U，V，W）平行的坐標(biāo)系（X，Y，Z）。模型如圖1所示。

圖1 基于雷達的蒼蠅微動物理模型

O在雷達坐標(biāo)系中，其方位角和仰角分別為α和β，O與雷達的距離為R0。因此振動中心O點的坐標(biāo)可以表示為：

雷達發(fā)射方向的單位向量為：

假設(shè)在雷達坐標(biāo)系中，目標(biāo)的振動方向為：方位角αp，仰角βp，最大振幅Dv，振動頻率fv，因此目標(biāo)的振動形式可以表示為：

振動方向的單位矢量是：

任意時刻，目標(biāo)與雷達的距離為：

所以t時刻，目標(biāo)與雷達的距離：

通常情況下Rt＞＞Dt，則（6）式可以近似為

雷達回波：

其中，ρ是點目標(biāo)的反射率，f0是發(fā)射信號的載頻，表示相位調(diào)制函數(shù)。對Sr（t）利用K階Bessel函數(shù)，

展開可得微多普勒頻譜是以中心頻率f為中心的對稱分布，相鄰譜線的間隔為

t=0時刻，點目標(biāo)P在（X0，Y0，Z0），t0時刻移動到：

對應(yīng)的速度為：

根據(jù)式，對應(yīng)的微多普勒頻率為：

從式（12）可看出第一項為平動引起的多普勒頻率，第二項和第三項分別是由于加速運動和蒼蠅翅膀振動引起的微多普勒頻率。如果補償?shù)羝絼雍图铀龠\動，則翅膀振動引起的微多普勒頻率是：

式（12）表明：由平動引起的多普勒頻率是一個常量，由于加速運動引起的微多普勒頻率是隨時間成正比變化。而由于振動引起的微多普勒頻率是隨時間成周期性變化的，變化周期與振動周期相同，變化幅度與振動幅度和振動頻率成正比。由此可見先補償?shù)羝絼雍图铀龠\動帶來的多普勒頻率，然后就可以通過振動引起的微多普勒來提取各個微動的各個參數(shù)。

2 雷達微動特征提取技術(shù)

雷達目標(biāo)微動特征提取主要是從雷達回波信號中獲取反映目標(biāo)運動、結(jié)構(gòu)等信息的特征量，基于特征量實現(xiàn)對目標(biāo)微動參數(shù)估計，為目標(biāo)的識別提供基礎(chǔ)。

2.1 基于偽Winger-Ville變換的曲線參數(shù)提取

從公式（9）得到的雷達回波信號的形式可以看出，信號頻率是隨著時間而不斷變化的，即微多普勒信號是一個非平穩(wěn)時間信號。因此傳統(tǒng)的傅里葉變換方法（Fourier Transform）求出的頻率將不再能反應(yīng)信號頻率隨時間變換的情況。因此有效的提取微多普勒需要使用高分辨的聯(lián)合時間-頻率分析方法。時頻分析方法是提取目標(biāo)微動特征的有效方法之一，也是應(yīng)用最廣泛的一類方法，它能夠觀察到信號的頻率隨時間變換的規(guī)律。常見的時頻分析方法主要有兩種：一種是線性變換，短時傅里葉變換（STFT），小波變換 （Gabor）；另一種是雙線性變換，Winger-Ville分布，偽 Winger-Ville分布（PWVD），平滑偽Winger-Ville分布（SPWV）[5]。

Winger-Ville分布是時頻分析方法中最早問世的，它是一種基于自相關(guān)系數(shù)的非線性時頻分布，但WVD受交叉項干擾嚴重影響，故本文引入時域加窗的偽魏格納分布（Pseduo，PWVD，PWVD），其表達式如下：

對其進行離散化處理，其離散化的表達式為：

典型的偽Winger-Ville變換分布如圖3所示。本論文，從實際雷達布局出發(fā)，對公式（9）中的回波信號進行偽Winger-Ville變換，和圖2中的理想多普勒頻率圖對比可知，偽Winger-Ville變換可以得到振動目標(biāo)微多普勒效應(yīng)的輪廓。

圖2 偽winger-ville變換分布圖

2.2 譜圖峰值估計分析

從圖2知，我們不能夠精確的得到每一時刻對應(yīng)的微多普勒頻率，即瞬時頻率無法得到。同樣振動幅度也無法精確得知。而譜圖峰值估計算法估計的瞬時頻率是其峰值在時頻平面上的投影，是對目標(biāo)瞬時頻率的提取，進一步表征信號的特征參數(shù)，從而達到對目標(biāo)的識別。假設(shè)信號的時頻變換頻譜為：W（m，k），其中1≤m≤M，1≤k≤N，譜圖峰值估計算法原理為：

式中，arg{·}為計算在任意時間t，幅度達到最大值對應(yīng)的頻率值。

本論文對得到的偽Winger-Ville變換圖采用譜圖峰值估計法來提取目標(biāo)的瞬時頻率，從而確定信號譜圖峰值的具體位置。因此可以得到瞬時頻率，以及每一時刻對應(yīng)的信號幅度。

2.3 曲線擬合工具Curve Fitting Tool

Matlab有一個功能強大的曲線擬合工具箱cftool，能夠?qū)崿F(xiàn)多種類型的線性、非線性曲線擬合。我們要得到微動特征的各個參數(shù)，可以使用cftool方法進行數(shù)據(jù)擬合。cftool曲線擬合窗口如圖3所示。

圖3 Curve Fitting Tool窗口

式（12）表明，可以先用cftool中的直線擬合法f（x）=p1*x+p2來得到蒼蠅飛行的平動速度和加速度。綜合p1，p2和式（12），可估計出：

從而補償?shù)羝絼訋淼亩嗥绽疹l率和加速運動帶來的微多普勒頻率。之后再利用cftool中的正弦擬合法 f（x）=a1*sin（b1*x+c1），對應(yīng)的求解到振翅頻率fv和振動幅度Dv

綜上分析，綜合利用時頻分析法、參數(shù)擬合法結(jié)合擬合工具來快速提取識別蒼蠅的具體步驟如下：

步驟1構(gòu)建基于雷達的蒼蠅微動物理模型，理論上求得微動帶來的雷達回波信號和微多普勒頻率值，為下一步信號處理以及驗證理論正確性提供數(shù)據(jù)支持；

步驟2對雷達回波信號進行時頻分析，得到目標(biāo)微多普勒信號的偽Winger-Ville時頻分布骨架S（t，f）；

步驟3對S（t，f）采用譜圖峰值估計法，求得微多普勒效應(yīng)的瞬時頻率；

步驟4對經(jīng)過譜圖峰值估計法得到的瞬時頻率進行matlab線性參數(shù)擬合，利用式（16）求得蒼蠅運動的速度和加速度；

步驟5利用得到的速度v0，加速度a補償?shù)艟€性的多普勒頻率，在利用matlab正弦參數(shù)擬合得到振翅頻率fv和振動幅度Dv；

3 仿真實驗

3.1 微動參數(shù)估計

設(shè)雷達工作在X波段，載頻發(fā)射信號載頻f0=10 GHz。蒼蠅的振翅頻率約三百五十二次，飛行速度是6～8 km/h，所以為了驗證理論的正確性，本文假設(shè)振動頻率fv=350 Hz，振動幅度Dv=0.06 m，平動速度v0=2 m/s，加速度a=10 m/s2，蒼蠅沿著雷達視線方向（LOS）飛行和加速。對雷達接收到的回波信號進行偽Winger-Ville變換。信號采樣率是fs=50 000 Hz，畫圖時候的時間間隔選取0.05 s。PWVD仿真結(jié)果如下圖4所示，其中選擇的窗函數(shù)寬度為47。從圖中大致可以看出頻率隨時間變換的輪廓，但不能得到精確值。所以，仍需進一步處理。

對經(jīng)過偽Winger-Ville的信號進行譜圖峰值估計，譜圖峰值檢測之后得到單值的時頻曲線，峰值檢測后的瞬時頻率與式（12）對應(yīng)求得的理論多普勒頻率如圖5所示。

圖4 回波信號偽winger-ville變換圖

圖5 理想多普勒頻率和檢測出的瞬時頻率對比圖

由圖5可知，經(jīng)過譜圖峰值估計檢測到的瞬時頻率和理想的微多普勒頻率在幅度和頻率上相差很小，可以初步確定方法的正確性。但精確的數(shù)值從波形上無法用肉眼得到。對經(jīng)過峰值檢測后瞬時頻率波形進行采樣，采樣率仍然是fs1=50 000 Hz，為了能用肉眼看得更清晰一些，本文時間只選取前0.02s時間內(nèi)可得到10 000組數(shù)據(jù)，對采樣后得到的10 000組數(shù)據(jù)進行matlab曲線線性擬合f（x）=p1*x+p2。

擬合出的斜率和截距值分別是 p1=637.8；p2=137.1，所以根據(jù)式（16）可求得

補償?shù)羝絼雍图铀龠\動后，目標(biāo)振翅帶來的微多普勒頻率為：

其中，f（m）是經(jīng)過譜圖峰值估計檢測后的頻率隨時間變換的數(shù)據(jù)組。

經(jīng)過上式修正后的蒼蠅振翅帶來的微多普勒頻率如圖6所示。

平動補償后再利用matlab正弦擬合f（x）=a*sin（b*x+c），擬合出的數(shù)據(jù)值為：a=5648，b=2199，c=1.571

由參數(shù)b和式（14）可以得到

將振動頻率代入上式即可得到：

振動頻率Dv=0.066 516 m

圖6 平動補償后的時頻曲線

3.2 性能分析

微振動參數(shù)的求解質(zhì)和實際真實的理論值對比如表1，可知誤差控制在0.0052%-10%之間，其中振動頻率誤差只有0.0193 Hz，由此可見能夠利用微多普勒效應(yīng)來準(zhǔn)確的識別到振翅的蒼蠅，這對于食品行業(yè)生產(chǎn)加工儲藏過程中捕捉蒼蠅價值重大。同時，能夠檢測到蒼蠅飛行的速度、加速度，為進一步捕捉飛行中的蒼蠅提供數(shù)據(jù)支持。從表格可以看出，振動幅度DV誤差較大，原因可能是:1）振動幅度是通過先計算出振動頻率 fv，代入（19）中的fv本身就有一定的誤差，可能因此引入二次誤差；2）在推導(dǎo)數(shù)學(xué)原理公式過程中做了一些近似。

4 結(jié)束語

微動目標(biāo)的雷達回波相位中包含著目標(biāo)微動的運動信息，而特征提取和微多普勒現(xiàn)象是一個新的研究方向。目前，國內(nèi)外對于復(fù)合微運動帶來的微多普勒效應(yīng)研究較少，尤其是對于像蒼蠅振翅頻率約300多Hz頻率的微運動研究更少。本文提出了一種利用微多普勒效應(yīng)來識別蒼蠅的新方法，并為捕捉蒼蠅提供了數(shù)據(jù)支持。在提取微動參數(shù)過程中，綜合利用了時頻分析、譜圖峰值估計、平動補償、參數(shù)擬合的方法得到了微動所有參數(shù)確切的值，即振動頻率、震動幅度、平動速度、加速度，誤差控制在0.0052%～10%之間。

表1 微振動參數(shù)理論值和測量值

[1]JamesAustin，DawnCalibeo-Hayes.有效防治害蟲，捍衛(wèi)食品安全[J].食品安全導(dǎo)刊，2009（5）:71-72.

[2]周建波.無公害粘蟲膠在害蟲綜合治理中的應(yīng)用技術(shù)研究[D].保定：河北農(nóng)業(yè)大學(xué)，2008.

[3]Chen VC，LI F，HO S S，et al.Micro-Doppler effect in radar:Phenomenon，model and simulation study [J].IEEE Transactions on Aerospace and Electronic Systems，2006，42（1）:2.

[4]Chen VC.Micro-Doppler effect of micro-motion dynamics:a review[C]//Processing of SPIE on Independent Component Analysis，Wavelets，and Neural Networks[C].Orlando，USA:SPIE Press，2003.

[5]王寶帥，杜蘭，劉宏偉，等.基于經(jīng)驗?zāi)B(tài)分解的空中飛機目標(biāo)分類[J].電子與信息學(xué)報，2012（9）:2116-2121.

[6]關(guān)永勝，左群聲，劉宏偉，等.空間進動目標(biāo)微動參數(shù)估計方法[J].電子與信息學(xué)報，2011（10）:2427-2432.

[7]張群，羅迎，何勁.雷達目標(biāo)微多普勒效應(yīng)研究概述[J].空軍工程大學(xué)學(xué)報：自然科學(xué)版，2011，12（2）:22-26.

[8]LuoYing，ZhangQun，Qiu Cheng-wei，etal.Micro-Doppler effect analysis and feature extraction in ISAR imaging with stepped-frequency chirp signal[J]. IEEE transactions on Geoscience and Remote Sensing，2010，48（4）:2087-2098.

[9]肖金國，杜蘭，韓勛，等.基于參數(shù)化時頻分析的進動錐裙目標(biāo)瞬時微多普勒頻率提取方法[J].電子與信息學(xué)報，2016（1）:90-96.

[10]羅迎，張群，朱仁飛，等.多載頻MIMO雷達中目標(biāo)旋轉(zhuǎn)部件三維微動特征提取方法 [J].電子學(xué)報，2011（9）:1975-1981.

[11]何勁，張群，羅迎，等.逆合成孔徑成像激光雷達微多普勒效應(yīng)分析及特征提取 [J].電子學(xué)報，2011（9）:2052-2059.

[12]曹文杰，張磊，杜蘭，等.基于瞬時頻率估計的進動錐體目標(biāo)微多普勒頻率提取方法[J].電子與信息學(xué)報，2015（5）:1091-1096.

[13]Gao Hongwei，Xie Lianggui，Wen Shuliang，et al. Micro-Doppler signature extraction from ballistic target with micro-motions[J].IEEE Transactions on Aerospace and Electronic Systems，2010，46（4）: 1969-1982.

[14]梁必帥，張群，婁昊，等.基于微動特征關(guān)聯(lián)的空間非對稱自旋目標(biāo)雷達三維成像方法[J].電子與信息學(xué)報，2014（6）:1381-1388.

A method of recognizing flies based on radar micro-doppler effect

SHEN Fei-fei1，2，LI Li-gang1，JIA Pei-yang1，2，SONG Zhen-jian1，2
（1.National Space Science Center/Chinese Academy of Sciences，Beijing 100190，China；2.University of Chinese Academy of Sciences，Beijing 100049，China）

In order to effectivelyidentify flies，this paper proposes an effective method of using radar micro-Doppler effect to recognize flies.Firstly，according to the basic concept of the micro-Doppler effect，the radar echo model of composite motion of a fly is established，then the radar echo signal and theory Doppler frequency is derived.A method of comprehensive application of the frequency analysis（pseudo Winger-Ville transform），the peak spectrum estimation，parameterfitting，motion compensation is proposed.Finally，the motion parameters of the fly during the flight are extracted. Including:flapping frequency vibration amplitude，velocity，and the acceleration.

MicroDoppler；compositemotion；pseudo Winger-Ville transform；parameter fitting；peak spectrum estimation；motioncompensation

TN95

A

1674-6236（2017）09-0022-05

2016-07-02稿件編號：201607009

中科院創(chuàng)新基金（CXJJ-15S146）

沈菲菲（1992—），女，河南濮陽人，碩士研究生。研究方向：計算機仿真以及信號處理。