莫禹涵,郭立民
哈爾濱工程大學 信息與通信工程學院 黑龍江 哈爾濱 150001
近年來隨著科技的發展,小型無人機逐漸出現在大眾的視野中。為大眾帶來方便的同時,也有些非合作無人機會給一些敏感地區的低空安全帶來很大的威脅。與其他無人機反制方案相比,雷達探測相比于其他手段具有全天候、作用距離大的優勢[1]。但低、慢、小目標的探測一直是雷達信號處理的難點,尤其是在無人機懸停的場景下,無人機回波在多普勒頻率上和地雜波混合在一起,為目標檢測帶來了很大的困難。不過無人機高速轉動的旋翼會給雷達回波帶來額外的頻率調制,使其多普勒特征不同于一般物體的雷達回波[2-3]。美國的著名學者Victor.C.Chen[4]在其著作中將其命名為“微多普勒效應”。目前國內和國際上對于無人機旋翼回波特征的研究大多在仿真階段,有實測數據的研究也都是對算法的佐證,實時性不強[5]。因此本文根據無人機回波的特點,結合相關文獻的研究結果,分析了無人機旋翼回波的頻域特性,根據無人機回波特點,選擇在FPGA 上構造一階遞歸型MTI 濾波器對背景雜波進行抑制,利用第三代雙倍速率(double data rate, DDR3)存儲器實現多周期雷達回波的MTD處理,在處理速度和吞吐率上進行了一定的優化,將FPGA 處理后的距離–多普勒(range-Doppler,R-D)二維譜和理論R-D 譜比較分析,證明FPGA處理結果的正確性。本文所述LFMCW 雷達系統的信號處理板是基于FPGA+DSP 架構的,FPGA負責信號的高速采集和處理,DSP 負責目標的參數估計以及系統主控。其中FPGA 為Xilinx V7系列的XC7V585T-2FFG1761I,開發平臺為Vivado 2019.2。
考慮到無人機的低、慢、小特性,本文選用的雷達體制為LFMCW 雷達。因為其無距離盲區,適合觀測近處無人機目標。而相比單頻連續波雷達,LFMCW 雷達不僅能夠觀測到目標的多普勒特性,還能得到目標的距離信息,更適用于監測無人機狀態的雷達應用場景。根據文獻[6-9]中的結論可知,無人機的時域回波特征是“閃爍”效應,理論上講提取時域特征對雷達硬件的要求更低。但在實際場景中,時域特征不如頻域特征穩定和易于分析。故本文僅分析無人機回波的頻域特征。無人機旋翼在雷達觀測下的幾何示意圖如圖1 所示。

圖1 無人機旋翼散射點示意
由圖1 可知,無人機距離雷達的方位角為α,俯仰角為β,模擬角速度為 Ω。假設無人機旋翼長度為l,旋翼上任意一與旋翼中心O′距離為r的散射點P在t時刻與雷達的距離為Rp(t),初相位為φ,旋翼中心距離雷達R0。由幾何關系可得到其表達式為

忽略初相,對P點處的回波與發射信號進行差頻處理并化簡后,得到差頻信號表達式:

式中:μ=B/T為調頻斜率,其中,B為調頻帶寬,T為調頻周期;c為光速。
對整個葉片的散射點做積分可得到整個旋翼的回波表達式:

對于多葉片的情況,回波表達式僅初相φ 的值不同。設一個旋翼上的葉片數目為N,則第k個葉的初相φk為

本文研究的無人機為大疆精靈3 標準版,葉片數為2。下面代入實際雷達參數對無人機旋翼回波進行仿真。仿真參數如表1 所示。

表1 仿真參數
需要說明的是,為減少LFMCW 雷達距離–速度耦合的影響,本文雷達選用三角調頻的方式,上下調頻周期各1 ms。后文給出的仿真和實測均為上掃頻結果。為保證差頻信號頻譜的純凈,需要將1 ms 內2 500 個時域采樣點頭尾數據刪去,用中間段2 048 個點做FFT。暫不考慮葉片的實際散射截面積(radar cross section,RCS)和噪聲,將仿真參數代入式(1),經過快時間距離維快速傅里葉變換(fast Fourier transform,FFT)和慢時間多普勒維FFT 后得到無人機旋翼回波的R-D 譜。仿真結果如圖2 所示。

圖2 單個無人機旋翼的R-D 譜
從圖2 的仿真結果可以看出,無人機旋翼回波的譜峰在多普勒域上是成對出現的,并關于零頻對稱。這與文獻[2]中理論分析的結果一致。譜峰間隔為54.72 Hz,按照文獻[10]分析的結果,譜峰間隔Δf應為Δf=NBfrot。按此計算無人機旋翼轉速應為27.36 r/s,和仿真參數相近。
本文選擇在FPGA 上而非DSP 上實現MTI和MTD 處理,目的是提高算法的運算速度,為DSP進行更加復雜的算法爭取時間。下面分別敘述系統整體結構和數據處理流程,以及MTI、MTD的FPGA 實現方法。
LFMCW 雷達信號處理包括一維FFT、MTI雜波抑制和MTD 相參積累處理。一維FFT 的作用是獲取包含在回波頻域的目標距離信息,同時提高目標信號的信噪比增益,類似脈沖雷達的脈沖壓縮處理。MTI 則是通過構造多普勒域高通濾波器,達到抑制零頻雜波的目的。MTD 可以視為是多普勒濾波器組,將不同多普勒頻率的信號進行同相疊加輸出,得到目標的多普勒信息,區分開不同多普勒頻率的目標,也是獲取無人機微多普勒特征的重要手段。系統的整體結構如圖3所示。

圖3 系統整體結構
回波信號經天線進入到接收機,與發射機輸入到接收機的本振信號混頻、濾波,得到差頻信號。數模轉換器(analog to digital converter,ADC)將模擬差頻信號轉化為數字信號后將數字信號傳給FPGA,FPGA 對串行的AD 原始數據進行串并轉換、降采樣、數字下變頻,得到同相(I 路)和正交(Q 路)分量。參照FPGA 內部產生的上下掃頻指示信號,將上下2 個掃頻段的連續波信號分別進行加漢明窗抑制旁瓣、距離維FFT 運算、MTI雜波抑制、MTD 相參積累。最后通過串行高速總線(serial rapidIO,SRIO)接口將MTD 后的數據發送給DSP 進行后續信號處理。
本文雷達接收機為零中頻結構,直接采集實信號。若直接對實信號做FFT,則會產生鏡像頻率分量,即在頻域有多個峰值。因此需要對回波進行數字下變頻得到復信號,經過低通濾波器濾除鏡像頻率后再進行FFT 運算。數字下變頻過程中,差頻信號的頻譜變化如圖4 所示。

圖4 頻信號頻譜變化過程
圖中|S(f)|為差頻信號頻譜的幅度,f1為一個目標距離對應的差頻信號頻率。實線為正頻率譜峰,虛線為鏡像頻率的譜峰。系統中AD 采樣頻率為60 MHz,經過降采樣處理后數據速率下降到2.5 Hz,與圖中fs相對應。將降采樣后的2.5 Hz數據在2.5 MHz 時鐘下分別乘1、0、-1、0,即與cos(2πfst/4)相 乘,可以得到I 路輸出;分別乘0、-1、0、1,即與-sin(2πfst/4)相乘,可以得到Q 路輸出。在頻域上的頻譜搬移效果如圖4 所示,關于0 頻率對稱的譜峰被搬移到關于-fs/4 對稱。在經過截止頻率為fs/4 的低通濾波器后,鏡像頻率即被抑制,目標距離所對應的頻率由0~fs/2 變為-fs/4~fs/4。在經過FFT 后僅對這部分頻率的FFT 結果進行后續處理即可。
MTI 可由相鄰脈沖對消來實現。常用的MTI濾波器有兩脈沖對消、三脈沖對消等。對消次數越多,濾波器在零頻處的凹陷越深,凹口越大,能夠濾除的雜波范圍也更大。但脈沖對消的缺陷也很明顯。首先是通帶不夠平坦,而且對消次數越多越不平坦,無人機的旋翼信息很可能因此損失,恒虛警檢測時也容易將背景估計過高或過低,造成漏警或虛警;其次是對于無人機這種低、慢、小目標,其低速運動時的多普勒頻率很低,很容易被MTI 濾波器抑制,造成雷達跟蹤過程中航跡斷裂。若想解決這個問題,需要實現對MTI 濾波器的凹口大小和平坦度的實時調整。文獻[11]中提到了一種遞歸型MTI 濾波器的構造方法,在兩脈沖對消器基礎上引入了一條反饋支路,反饋系數為K。通過給系統函數增加新的極點z=K的方法實現了頻率響應的平坦化處理。遞歸MTI濾波器的結構如圖5 所示。

圖5 遞歸型MTI 濾波器的結構
根據遞歸型MTI 濾波器的結構可以得到其頻率響應,K取0.8 時,其頻率響應和脈沖對消器的對比如圖6 所示。兩脈沖對消器的極點為z=0,遞歸型MTI 濾波器的極點為z=0 和z=K。K的取值越接近1,2 個極點的距離越遠,系統的頻率響應就越平坦,零頻凹口就越小;反之系統頻率響應就越不平坦,零頻凹口就越大。在LFMCW 雷達系統工作時,可以根據雜波幅度大小和雜波多普勒頻率的范圍,靈活地更改K的值以獲取不同的雜波抑制效果,更適合本文所述的無人機探測場景。

圖6 不同MTI 濾波器的頻率響應對比
在FPGA 上實現遞歸型MTI 濾波器需要進行相鄰周期脈沖的相減、反饋支路數據的保存和與反饋系數相乘等工作。本文選用雙端口隨機存儲器(radom access memory,RAM)來緩存回波數據,經過前面預處理的回波經過波門選通邏輯寫入不同的RAM,反饋支路也選用雙端口RAM 保存。本文在從端口B 讀取出反饋支路RAM 數據的同時,將新計算的MTI 處理結果通過端口A 寫入RAM。無論是計算MTI 還是更新反饋支路數據,都采用流水線處理的方式,最大限度地提高處理速率。綜上,在FPGA 上實現MTI 處理的信號流向,如圖7 所示。

圖7 FPGA 實現MTI 信號流向
為驗證FPGA 程序設計的正確性,需要借助ModelSim 平臺對設計進行仿真。仿真時輸入的3 路回波數據為實測的室內開機下的雷達回波。FPGA 對回波進行數字下變頻、距離維FFT 后,不經后續處理直接截取每個回波的前256 個距離門數據上傳上位機。ModelSim 仿真時僅需要在測試文件中給出時鐘和輸入的一維FFT 數據,其余參數都由FPGA 程序設計決定,與表1 中參數相同。雜波數據如圖8 所示。

圖8 用于仿真的實測雜波數據
從圖8 中可以看出雜波的峰值出現在第3 個距離采樣點,且能量分布在1~17 個距離采樣點,遠遠高于背景噪聲,且不同周期的雜波幅度相近,符合固定目標回波特性,也間接證明了硬件系統的有效性和FPGA 前端處理的正確性。將雜波數據作為激勵輸入到testbench 中,得到的仿真結果如圖9 所示。

圖9 MTI 處理ModelSim 仿真
圖9 是以一次上掃頻的MTI 處理為例,當新的回波數據到來時,FPGA 先將數據寫入上掃頻RAM1。寫完后將上掃頻RAM1 和RAM2 的數據一同讀出輸入到減法器中,與此同時讀取上掃頻反饋支路RAM 中的數據乘以系數0xCCCC(即為16 位無符號定點小數0.8),再和減法器輸出的數據一同輸入到加法器中。最后得到的add_I_out和add_Q_out 的結果即為最后MTI 處理后的回波數據。可以看到在同一顯示尺度下,原始回波I、Q 路的高強度雜波已經得到了濾除。另外需要說明的是,反饋支路RAM 存在同時讀寫的情況,但新的MTI 結果一定是在讀取反饋支路RAM 的數據后才會輸出,所以并不存在對RAM 中同一個地址同時讀寫的情況,局部放大的結果如圖10所示。

圖10 反饋支路RAM 讀寫局部放大
圖10 中addra,wea,dina 分別對應雙端口RAM 的A 端口的地址、寫使能和數據線。同理另外3 組信號對應B 端口的地址、讀使能和數據。可以看到寫使能wea 是在讀使能enb 拉高后幾個時鐘才拉高。
對以上仿真結果進行時鐘計數,發現通過流水線并行處理,FPGA 能夠在10 μs 左右完成對雷達回波的遞歸MTI 處理,達到快速實現MTI 的目的。
對于傳統的信號處理機,MTD 處理多由DSP或嵌入式進階精簡指令集機器(advanced RISC machine,ARM)等處理器去完成。但串行邏輯在處理這種大量重復運算時的效率會很低,僅MTD處理就能占用大量的時間資源,導致后續留給參數估計的計算時間很有限。而一些在FPGA 上實現MTD 的方法只能做到固定點數,與串行處理器相比又少了靈活性。本文針對這些不足做出了一些改進,利用DDR3 存儲器在FPGA 上實現了點數可變的快速MTD 處理。具體時序操作由設計有限狀態機來完成。其信號流向如圖11 所示。

圖11 FPGA 實現MTD 的信號流向
做MTD 處理就需要先緩存一個相干處理間隔(coherent process interval,CPI)的回波數據。為實現大量數據的緩存以及突破速度瓶頸。本文選擇利用外掛的DDR3 存儲器來實現MTD 處理。其具有容量大、時鐘效率高等優勢。本文所述雷達系統的FPGA 外掛了容量為2 GB 的DDR3 內存條,最多能緩存近3 min 的回波數據,能夠實現MTD 處理的要求。
FPGA 對DDR3 存儲器的讀寫由Xilinx 的存儲器接口生成器(memory interface generator,MIG)IP 核來實現,由于DDR3 的讀寫方式有突發長度的限制,一次需要連續寫入或讀出8 個地址中的數據。所以對應到用戶接口時的數據位寬就是DDR3 實際位寬的8 倍。本文的設計中MTI 后的數據I 路和Q 路各32 位,總共64 位。對應到DDR3 的用戶接口的數據位寬就變成了512 位,位寬和時鐘域的轉換同樣可以通過雙端口RAM 來實現。位寬轉換后,只需要一個DDR3 用戶層時鐘就能將8 個距離采樣點的數據寫入存儲器。若用戶層時鐘為100 MHz,即相當于實現了數據以800 MHz 的速率進行存儲,遠超FPGA 的片內存儲器。本文即通過這種方式實現了數據的加速處理。
讀取數據和寫入時一致,一次讀出8 個有效數據。不同的是需要按照同一個距離維采樣點讀出數據。Xilinx 的FFT IP 核支持多通道FFT 且點數可實時配置。在FFT 處理前還需要加窗以減少旁瓣幅度。為充分利用FPGA 并行處理的優勢,本文將FFT 核配置為8 通道定點FFT。數據從DDR3中被讀出后從加窗到FFT 采用全流水并行操作的方式,獲得遠超于DSP 處理的速度和吞吐率。在實現點數可變的處理方面,FPGA 需要從DSP 得到后續信號處理所需的最佳MTD 點數,在下一個CPI 開始前配置FFT 的點數,修改只讀存儲器(read only memory,ROM)漢明窗表的尋址步長,進而實現窗函數和FFT 運算根據點數實時更新。定點數在經過FFT 后的數據位寬會變大,若截位處理又會損失一部分精度和動態范圍。所以本文在FFT 后將大位寬的定點數轉成32 位單精度浮點數,再通過高速SRIO 接口將結果給DSP 做后續處理。
在一些脈沖重復時間(pulse repetition time,PRT)較小的場景,FPGA 無法在一個PRT 內完成所有回波的MTD 處理。此時DDR3 會產生讀寫沖突。傳統的解決思路是用2 片DDR3 內存條實現乒乓緩存來獲得時間差。但這樣會造成存儲資源的浪費,1 片DDR3 完全能夠緩存2 個CPI 的回波數據。本文為了優化這個問題,利用了DDR3的塊、行、列存儲結構,將同一片DDR3 中的2 個塊作為乒乓緩存的2 個存儲區。對于讀寫的總線沖突問題,本文在設計中加入了寫請求中斷的機制。由于數據寫入是猝發的,本身不會持續很長時間,所以在讀數據的同時若接到寫請求,狀態機會在讀完這一次數據后跳到寫數據的狀態中,先將此次數據緩存到DDR3 后,再繼續執行讀上一個CPI 數據的操作。MTD 處理的有限狀態機狀態轉換如圖12 所示。

圖12 MTD 處理狀態轉換
將2.2 節MTI 處理后的數據輸入到MTD 模塊在ModelSim 進行仿真。為節省仿真時間,MTD點數取32,得到的結果如圖13 所示。

圖13 DDR3 讀寫時序仿真
從時序仿真中能夠看到,在狀態機剛剛進行完一次讀操作時,寫請求(ddr_flag)到達。狀態機則控制數據寫入到DDR3 的另一個塊內(2 個塊首地址分別為0x00 開頭和0x02 開頭,從圖11 中可以看出)。對MTD 仿真進行時鐘計數,得到FPGA 對三通道128 周期(上下掃頻各64 周期)的數據進行MTD 處理的時間為1.54 ms。將仿真后的結果保存為txt 文件,在Matlab 中畫出32 點上掃頻MTD 結果如圖14 所示。

圖14 ModelSim 仿真的32 點MTD 結果
從仿真結果可以看出本應在多普勒零頻處的雜波已經被濾除干凈,且通帶較為平坦。證明了本文改進的MTI、MTD 實現方法的正確性和有效性。
為模擬強雜波場景,測試在室內進行。先將無人機放在靜止的平臺上,僅安裝一個旋翼。用遙控器控制其旋翼轉動,觀察旋翼速度對其多普勒特征的影響。根據文獻[3]的結果,對于無人機目標來說,1 個CPI 最好不超過0.2 s。因此在測試時將MTD 點數確定為64,即一個CPI 的時間為0.128 s。MTI 濾波器反饋支路的系數選為0.8。FPGA 計算得到的處理結果可從DSP 端通過仿真器讀取到PC 機上,圖15 為不同情況下的FPGA處理結果。

圖15 單旋翼不同轉速FPGA 的處理結果
作為對比,圖15(a)為無人機旋翼不旋轉時FPGA 的處理結果,即純雜波的處理結果。可以看到R-D 譜中無明顯峰值,背景幅度大約75 dB。圖15(b)、(c),(d)分別為無人機旋翼正常轉速、慢速、快速的處理結果。可以發現旋翼的多普勒特征與仿真類似,為關于零頻對稱的譜線。慢速旋轉時譜峰間隔變小,快速旋轉時譜峰間隔變大。但由于室內環境封閉,系統受多徑影響較大,造成能量峰值處的距離比實際距離遠。且系統PRT 較大,造成多普勒模糊,使得R-D 譜除了明顯的峰值以外還存在一些離散的譜線。為了分析不同旋翼數目和無人機不同運動狀態的影響,本文又對比了兩旋翼和四旋翼的情況以及無人機空中懸停和徑向飛行的情況,結果如圖16 所示。

圖16 無人機不同旋翼數目和飛行狀態的測試結果
兩旋翼與四旋翼的轉速與圖15(b)中單旋翼中速旋轉的轉速一致,可以看到旋翼數目增加后,離散的譜峰變得不再明顯。此時旋翼的多普勒頻率已經發生了很大的展寬,幾個旋翼的譜峰混合在了一起。造成這種現象的原因有2 個:一是由于旋翼的角度有差異,即便實際轉速相同,投影到雷達徑向的速度分量也會存在差異;二是雷達系統PRT 過大導致多普勒無模糊范圍較小,造成旋翼多普勒維信息的混疊。文獻[12]也對這種現象進行了分析。
對比圖16(c)和圖16(d),無人機懸停時本身不運動,但由于氣流擾動,實際上還是存在微小的位移。由于MTI 濾波器的凹口較小,導致一部分無人機主體的信息還是被保留了下來。旋翼信息方面,無人機懸停時的旋翼轉速要比靜止在平臺上時大,故距離零頻最近的峰值所對應的多普勒頻率也更大。而無人機在徑向運動時,其主體也有了速度,多普勒頻率位于MTI 濾波器的通帶范圍內。所以可以看到圖16(d)的目標峰值要比前面的旋翼以及懸停時的無人機主體能量要大得多。但運動的無人機旋翼轉速也要比靜止和懸停時大得多,故產生的多普勒模糊就更為嚴重。
1)從實際測試和仿真的結果能夠看到本文的方法利用FPGA 的并行處理優勢。通過流水線處理和乒乓操作以及DDR3 存儲器的一些獨特優勢,實現了對地雜波的抑制、對目標回波的相參積累和微多普勒信息顯示。一定程度上提高了系統的吞吐率,也緩解了DSP 的資源壓力。實測的旋翼多普勒特征也接近理論分析的結果,進一步證實了用FPGA 實時處理的可行性。
2)目前關于無人機檢測和識別的研究多在算法以及性能仿真上,少有立足于實際的雷達系統和硬件平臺上的研究。本文對傳統的MTI、MTD硬件實現方法做出了一些速度和吞吐率上的優化,目的是要在實際的系統上實現對無人機目標的回波特征,尤其是微多普勒特征的有效、快速和實時的提取,為雷達系統快速發現和識別無人機提供了一種思路。
3)從運動無人機的處理結果也可以看到,受硬件平臺參數所限,系統的多普勒無模糊測量范圍很小,無人機速度過快時會發生多普勒域的混疊。要解決這個問題可以適當減小PRT 以獲得更大的多普勒頻率測量范圍、更大的積累周期以及多普勒分辨力;也可以增大掃頻帶寬獲得更小的距離分辨力,進一步優化R-D 譜的精度,得到精確的旋翼參數估計結果。