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

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

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

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

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

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

微多普勒；復(fù)合運(yùn)動(dòng)；偽Winger-Ville變換；參數(shù)擬合；譜圖峰值估計(jì)；平動(dòng)補(bǔ)償

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

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

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

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

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

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

圖1 基于雷達(dá)的蒼蠅微動(dòng)物理模型

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

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

假設(shè)在雷達(dá)坐標(biāo)系中，目標(biāo)的振動(dòng)方向?yàn)椋悍轿唤铅羛，仰角βp，最大振幅Dv，振動(dòng)頻率fv，因此目標(biāo)的振動(dòng)形式可以表示為：

振動(dòng)方向的單位矢量是：

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

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

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

雷達(dá)回波：

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

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

t=0時(shí)刻，點(diǎn)目標(biāo)P在（X0，Y0，Z0），t0時(shí)刻移動(dòng)到：

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

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

從式（12）可看出第一項(xiàng)為平動(dòng)引起的多普勒頻率，第二項(xiàng)和第三項(xiàng)分別是由于加速運(yùn)動(dòng)和蒼蠅翅膀振動(dòng)引起的微多普勒頻率。如果補(bǔ)償?shù)羝絼?dòng)和加速運(yùn)動(dòng)，則翅膀振動(dòng)引起的微多普勒頻率是：

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

2 雷達(dá)微動(dòng)特征提取技術(shù)

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

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

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

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

對(duì)其進(jìn)行離散化處理，其離散化的表達(dá)式為：

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

圖2 偽winger-ville變換分布圖

2.2 譜圖峰值估計(jì)分析

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

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

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

2.3 曲線擬合工具Curve Fitting Tool

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

圖3 Curve Fitting Tool窗口

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

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

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

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

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

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

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

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

3 仿真實(shí)驗(yàn)

3.1 微動(dòng)參數(shù)估計(jì)

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

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

圖4 回波信號(hào)偽winger-ville變換圖

圖5 理想多普勒頻率和檢測(cè)出的瞬時(shí)頻率對(duì)比圖

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

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

補(bǔ)償?shù)羝絼?dòng)和加速運(yùn)動(dòng)后，目標(biāo)振翅帶來的微多普勒頻率為：

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

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

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

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

將振動(dòng)頻率代入上式即可得到：

振動(dòng)頻率Dv=0.066 516 m

圖6 平動(dòng)補(bǔ)償后的時(shí)頻曲線

3.2 性能分析

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

4 結(jié)束語

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

表1 微振動(dòng)參數(shù)理論值和測(cè)量值

[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)驗(yàn)?zāi)B(tài)分解的空中飛機(jī)目標(biāo)分類[J].電子與信息學(xué)報(bào)，2012（9）:2116-2121.

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

[7]張群，羅迎，何勁.雷達(dá)目標(biāo)微多普勒效應(yīng)研究概述[J].空軍工程大學(xué)學(xué)報(bào)：自然科學(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ù)化時(shí)頻分析的進(jìn)動(dòng)錐裙目標(biāo)瞬時(shí)微多普勒頻率提取方法[J].電子與信息學(xué)報(bào)，2016（1）:90-96.

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

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

[12]曹文杰，張磊，杜蘭，等.基于瞬時(shí)頻率估計(jì)的進(jìn)動(dòng)錐體目標(biāo)微多普勒頻率提取方法[J].電子與信息學(xué)報(bào)，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]梁必帥，張群，婁昊，等.基于微動(dòng)特征關(guān)聯(lián)的空間非對(duì)稱自旋目標(biāo)雷達(dá)三維成像方法[J].電子與信息學(xué)報(bào)，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稿件編號(hào)：201607009

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

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