FPGA實(shí)時多普勒中心頻率估計(jì)

甘翼，李向陽

(中國電子科技集團(tuán)公司第十研究所，四川 成都 610036)

0 引言

合成孔徑雷達(dá)(synthetic aperture radar,SAR)受氣象條件的影響較小，具有全天候工作的能力，在測繪、軍事等領(lǐng)域具有廣范的應(yīng)用[1-2]。彈載SAR能夠?qū)崿F(xiàn)對目標(biāo)的快速高精度成像和持續(xù)跟蹤，是導(dǎo)彈實(shí)現(xiàn)遠(yuǎn)距離精確打擊的重要傳感器之一[3]。彈載SAR相對其他平臺體積小、質(zhì)量輕、功耗低，并且對處理實(shí)時性和大場景處理能力提出了更高的要求。因此，在彈載SAR實(shí)現(xiàn)方面，基于純FPGA的SAR成像實(shí)現(xiàn)方法，可以提高大場景處理能力、減小成像時間、提高系統(tǒng)集成度、降低系統(tǒng)功耗[4]，在并行處理能力和成像時間方面相對DSP處理器具備極大的優(yōu)勢[5]。

SAR成像算法包括距離向壓縮、距離走動矯正、方位向壓縮以及多普勒中心頻率估計(jì)、多普勒調(diào)頻率估計(jì)等模塊。其中距離向壓縮、距離走動矯正和方位向壓縮主要完成FFT變換、點(diǎn)乘和IFFT變換，具有單一的數(shù)據(jù)處理結(jié)構(gòu)，無迭代操作，因而適合采用FPGA(field-programmable gate array)進(jìn)行流水化處理。多普勒中心頻率(frequency of Doppler centroid,FDC)估計(jì)、多普勒調(diào)頻率(frequency of Doppler rate,FDR)估計(jì)[6]，因?yàn)樗惴ㄏ鄬^為復(fù)雜，使用FPGA實(shí)現(xiàn)具有較高的難度[7]，該部分通常采用DSP處理器完成[8]。這種折中方案會引入的數(shù)據(jù)交換延遲和FPGA處理等待延遲，大幅增加系統(tǒng)成像處理時間。以常用的DSP實(shí)現(xiàn)多普勒參數(shù)估計(jì)流程為例：首先，為了在DSP上實(shí)現(xiàn)參數(shù)估計(jì)，需要FPGA向DSP通過SRIO發(fā)送脈壓數(shù)據(jù)，而距離向4096、方位向2048的IQ單精度浮點(diǎn)數(shù)據(jù)高達(dá)64 MB，通過2Lane@5Gbps的SRIO單次傳輸需要64 ms；其次，在DSP進(jìn)行FDC估計(jì)時，因?yàn)楹罄m(xù)距離走動校正需要多普勒中心頻率估計(jì)結(jié)果，所以FPGA需要暫停當(dāng)前流程，等待DSP的FDC處理結(jié)果才能進(jìn)行后續(xù)操作，通常會產(chǎn)生高達(dá)80 ms的處理等待延遲；因此，基于FPGA+DSP的SAR處理流程在彈載條件硬件資源受限的情況下，無法滿足大場景實(shí)時成像的要求，只能通過減小成像區(qū)域，來滿足系統(tǒng)對實(shí)時性的要求。而通過將FDC估計(jì)采用FPGA實(shí)現(xiàn)即可大幅降低系統(tǒng)處理時間，提高系統(tǒng)的實(shí)時性。文獻(xiàn)[9]提出了一種基于Zynq系列FPGA和RD算法的SAR成像實(shí)現(xiàn)方法，其中運(yùn)動參數(shù)估計(jì)、多普勒參數(shù)估計(jì)、多普勒調(diào)頻率估計(jì)等放在ARM處理器上進(jìn)行處理，而距離向壓縮、方位向壓縮以及距離走動校正則是放在FPGA上進(jìn)行。文獻(xiàn)[10-12]提出了3種基于FPGA的SAR實(shí)現(xiàn)方法，其中文獻(xiàn)[10]沒有對參數(shù)估計(jì)部分進(jìn)行闡述，文獻(xiàn)[11-12]分別提出了基于FPGA和CS算法的SAR成像實(shí)現(xiàn)方法，并且采用FPGA實(shí)現(xiàn)多普勒參數(shù)估計(jì)。但是在多普勒參數(shù)估計(jì)部分均采用一次曲線進(jìn)行多普勒中心頻率擬合。而在彈載大斜視系統(tǒng)中，為了提高多普勒中心頻率估計(jì)的準(zhǔn)確性，通常采用2次或者3次曲線擬合進(jìn)行逼近[13]。

本文提出一種基于FPGA的FDC估計(jì)方法，通過調(diào)整FDC數(shù)據(jù)處理流程、實(shí)現(xiàn)了SAR信號處理機(jī)的“零等待”。在FDC估計(jì)算法實(shí)現(xiàn)過程中，提出了基于FPGA的流水線式相位解纏繞方法和基于最小二乘的參數(shù)曲線擬合算法，實(shí)現(xiàn)相位解纏繞的連續(xù)流水輸出以及中心頻率估計(jì)參數(shù)精確插值。

1 多普勒中心頻率估計(jì)算法

多普勒中心頻率的估計(jì)有2種途徑：①通過慣導(dǎo)獲得[14]；②通過算法從回波數(shù)據(jù)中提取。由于慣導(dǎo)本身存在測量誤差，為了獲得更好的聚焦性能，通常從回波數(shù)據(jù)中通過算法提取多普勒中心頻率參數(shù)15]。利用回波數(shù)據(jù)對多普勒中心頻率進(jìn)行估計(jì)有多種方法[16]。其中，能量均衡法、匹配相關(guān)法和最大似然法需要把距離向功率譜與某一參考函數(shù)做相關(guān)，把相關(guān)函數(shù)的零點(diǎn)作為多普勒中心[17]；相關(guān)函數(shù)法和數(shù)據(jù)序列符號自相關(guān)法則利用回波數(shù)據(jù)中的相位信息通過自相關(guān)提取多普勒中心，不需要額外的參考函數(shù)[18]?；贔PGA的可實(shí)現(xiàn)性，本文采用相位函數(shù)法實(shí)現(xiàn)多普勒中心估計(jì)[19]。由維納-辛欽定理可知，相關(guān)函數(shù)r(τ)可以表示為其功率譜函數(shù)E{S(f)}的IDFT變換。

r(τ)=F-1{E[S(f)]}.

(1)

對于雷達(dá)發(fā)射系統(tǒng)，假設(shè)距離向功率譜函數(shù)E{S(f)}為天線功率方向圖W(f)，則

r(τ)=F-1{W(f-fdc)}=r0(τ)exp{j2πfdcτ}.

(2)

從式(2)中可以看出自相關(guān)函數(shù)r(τ)是多普勒中心fdc的函數(shù)，當(dāng)采用與k個方位向上的數(shù)據(jù)進(jìn)行自相關(guān)時，也即τ=kT，其中T=1/PRF為重頻周期，fdc可以表示為

(3)

不同方位向上回波S(nT)的自相關(guān)函數(shù)可以表示為

r(kT)=E{s*(nT)s(nT+kT)}∝

(4)

從式(3)和式(4)通過不同方位向上相同距離門的回波函數(shù)相關(guān)，即可得到fdc。因?yàn)榉敲黠@場景中無法提取正確的多普勒信息，需根據(jù)回波能量對距離門進(jìn)行篩選，截取有效的場景數(shù)據(jù)進(jìn)行多普勒相關(guān)估計(jì)。同時為了降低計(jì)算量，在進(jìn)行方位向相關(guān)之前，會對距離門上的數(shù)據(jù)進(jìn)行抽取，最終只能得到部分多普勒頻率估計(jì)結(jié)果。為了得到完整的多普勒中心頻率，提高多普勒中心頻率估計(jì)的正確性，需要對多普勒頻率估計(jì)結(jié)果進(jìn)行二次曲線擬合和插值處理，通過求均值最終得到多普勒中心頻率，用于距離走動校正。

由于FPGA內(nèi)部存儲資源的限制，無法保存所有的距離向數(shù)據(jù)，因而本文采用流水線的方位向相關(guān)累積方法，無需保存距離向數(shù)據(jù)，大幅降低對存儲資源的消耗。在基于最小二乘的二次曲線擬合中，需要進(jìn)行矩陣求逆操作，本文采用矩陣分解的形式，在FPGA上直接實(shí)現(xiàn)3×3矩陣的求逆操作。下面將詳細(xì)介紹基于FPGA的多普勒中心估計(jì)處理流程和方法。

2 “零等待”的多普勒中心估計(jì)處理流程

多普勒中心頻率估計(jì)包含距離向積累抽取、方位向相關(guān)積累、相位解纏繞和插值擬合等。FPGA的處理流程和DSP的處理流程[20]對比如圖1所示。在DSP中距離向脈壓塊數(shù)據(jù)的讀寫和多普勒中心頻率估計(jì)需順序進(jìn)行，相比FPGA實(shí)現(xiàn)，需等待一個孔徑或子孔徑的距離向脈壓數(shù)據(jù)全部存儲后，再完整讀取進(jìn)行多普勒參數(shù)估計(jì)，隨后還須增加一次DDR數(shù)據(jù)讀取時間，才能進(jìn)行RCMC處理。而RCMC依賴多普勒中心頻率的估計(jì)結(jié)果。而此時多普勒中心頻率估計(jì)尚未完成，因此，一個成像孔徑或子孔徑的距離向脈壓數(shù)據(jù)需要緩存。以孔徑大小為4 096×2 048個處理單元為例，距離向脈壓后的單精度IQ數(shù)據(jù)為64 MB，由于數(shù)據(jù)量較大，F(xiàn)PGA內(nèi)部BRAM資源無法滿足要求，需要轉(zhuǎn)儲到外部大容量DDR中，與前述DSP實(shí)現(xiàn)方法不同，F(xiàn)PGA可以發(fā)揮其并行處理的優(yōu)勢，實(shí)現(xiàn)邊存儲邊處理。如圖1所示，在對脈壓結(jié)果進(jìn)行存儲的過程中，即可同步進(jìn)行多普勒中心頻率的估計(jì)。距離向壓縮結(jié)果存儲完成時，即可同步完成多普勒中心頻率估計(jì)，隨即開展RCMC處理，實(shí)現(xiàn)系統(tǒng)數(shù)據(jù)處理的“零等待”，降低了系統(tǒng)的處理時間，同時也簡化了系統(tǒng)設(shè)計(jì)。

圖1 FPGA與DSP對比處理流程Fig.1 Comparison of DSP and FPGA processes

如圖2所示，基于FPGA的多普勒中心頻率估計(jì)可以采用流水線的方式進(jìn)行數(shù)據(jù)處理。傳統(tǒng)非流水線處理方式需要的時間為t=N(tddr+tfdc)，而流水線的處理方式其時間消耗為t=Ntddr+tfdc，其中tddr為DDR讀取一個距離向數(shù)據(jù)所需要的時間，tfdc為處理一個距離向所需要的時間，N為CPI周期內(nèi)方位向積累點(diǎn)數(shù)。在相同tddr和tfdc的情況下，通過流水線的處理方式，可以“隱藏”N-1次FDC處理時間，對于4 096個距離向脈壓數(shù)據(jù)，在數(shù)據(jù)接收完成后，只需額外的一個處理周期即可外出多普勒中心頻率的估計(jì)，亦即“隱藏”了4 095次FDC處理時間。

圖2 基于FPGA流水與非流水對比處理流程Fig.2 Comparison between FPGA-based pipelined and non-pipelined processes

以上是基于FPGA的多普勒中心估計(jì)采用的基本流程和優(yōu)化方法，下面將從距離向積累抽取、相位解纏繞、二次曲線擬合和插值3個方面詳細(xì)介紹實(shí)現(xiàn)過程。

2.1 距離向數(shù)據(jù)的積累抽取

為了提高多普勒中心頻率估計(jì)的魯棒性，降低系統(tǒng)噪聲和場景突變對參數(shù)估計(jì)的影響，采用鄰域積累實(shí)現(xiàn)距離向數(shù)據(jù)的平滑濾波。雖然會降低多普勒中心頻率估計(jì)的精度，但是在彈載SAR面臨的復(fù)雜電磁環(huán)境中可以提高算法的魯棒性。同時為了降低方位向相關(guān)的數(shù)據(jù)處理量，對積累后的數(shù)據(jù)進(jìn)行1/4抽取。如圖3所示，在FPGA的實(shí)現(xiàn)過程中，采用定點(diǎn)累加的方式完成四鄰域數(shù)據(jù)的積累。通過2 bit計(jì)數(shù)器對數(shù)據(jù)進(jìn)行循環(huán)處理，當(dāng)循環(huán)計(jì)數(shù)值為‘0’時，使用當(dāng)前距離門數(shù)據(jù)對累加器進(jìn)行初始化，其他情況下進(jìn)行距離門數(shù)據(jù)累加，在累加完成4組數(shù)據(jù)后取走處理結(jié)果。雖然理論上應(yīng)該采用均值來表示累積的結(jié)果，但是每個距離門上的鄰域累積具有相同的長度，所以直接采用累加結(jié)果表示平滑濾波結(jié)果，從而可以減少邏輯資源的消耗。

圖3 距離向平滑抽取Fig.3 Smooth extraction of range dimension

根據(jù)FPGA數(shù)據(jù)處理的特點(diǎn)，以及距離向脈沖壓縮的實(shí)現(xiàn)過程，脈壓后的數(shù)據(jù)與回波數(shù)據(jù)一樣是按照固定重復(fù)周期采用流水線的方式依次接收和處理。在多普勒中心頻率估計(jì)中，距離向累積抽取可以直接在數(shù)據(jù)流上直接完成。而方位向相關(guān)累積數(shù)據(jù)則需要緩存。

數(shù)據(jù)緩存在FPGA中可以通過LE(logic elements)資源實(shí)現(xiàn)，也可以通過內(nèi)部BRAM實(shí)現(xiàn)。LE寄存器具有操作簡單、讀取延時低等特點(diǎn)。但是對于4 096點(diǎn)的64 bit復(fù)數(shù)據(jù)，其資源消耗是巨大的，會占用xc7v690t芯片20%的LE資源。通過BRAM進(jìn)行緩存可以使用FPGA內(nèi)部的專用存儲資源，存儲容量相對較多，但是BRAM會產(chǎn)生1～2個時鐘周期的讀寫延遲。為了實(shí)現(xiàn)數(shù)據(jù)的流水線處理，需要進(jìn)行預(yù)讀取。對于距離向脈壓數(shù)據(jù)，在4鄰域平滑累加的4個時鐘周期內(nèi)，需要完成一次讀取和一次寫入操作。對于采用的2時鐘延時的BRAM，數(shù)據(jù)讀取操作需要在4鄰域求和計(jì)算輸出完成2個時鐘周期之前進(jìn)行，寫入操作則在求和計(jì)算結(jié)果之后，為了避免讀寫操作的沖突，通常安排在求和結(jié)束后的第2個時鐘周期寫入。而對于方位向相關(guān)積累來說，則在5個時鐘周期內(nèi)需完成讀取、累加和寫入操作。其處理時間為3個時鐘周期，預(yù)讀取需要2個時鐘周期。

如圖4所示，方位向相關(guān)積累需要在CPI內(nèi)的所有PRT的每個距離門上進(jìn)行，抽取后的每個距離門需要獨(dú)立的寄存器用于緩存方位向相關(guān)積累結(jié)果，對于4 096的距離向長度，需要2組1 024點(diǎn)深度的寄存器組用于分別存放IQ數(shù)據(jù)的相關(guān)累加結(jié)果。

圖4 方位向相關(guān)積累Fig.4 Correlation accumulation of azimuth range

方位向相關(guān)累積完成后，對存放在BRAM中的各個距離門累加結(jié)果依次讀取，作為cordic模塊的輸入，進(jìn)行arctan操作。cordic生成的只是相關(guān)累積結(jié)果的瞬時相位信息，范圍在[-π,π)之間，因而瞬時相位信息和多普勒頻率相位之間需要進(jìn)行解纏繞操作。

2.2 解纏繞

多普勒中心頻率反映了脈沖發(fā)射單元與目標(biāo)之間的徑向速度差。不同PRT之間的相關(guān)性反映了不同距離門上的多普勒變化率，這是一種相位累積效應(yīng)，這種累積結(jié)果通過反正切變換為角度信息時，映射在[-π,π)的范圍內(nèi)。解纏繞操作可以恢復(fù)相位對應(yīng)于角度上的累加信息。當(dāng)角度信息在-π和-π附近變化時，通過判斷鄰域數(shù)據(jù)的變化率，增加或者減去2π，來達(dá)到相位連續(xù)累加的目的。在得到多普勒相位值后，根據(jù)其對應(yīng)圓周角度的比例，通過PRT換算即可得到多普勒頻率。由于PRT的離散數(shù)據(jù)采集造成的多普勒模糊[21]，則通過慣導(dǎo)信息來提取多普勒模糊數(shù)進(jìn)行修正，從而可以得到實(shí)際的多普勒頻率值。解纏繞的FPGA實(shí)現(xiàn)偽代碼如表1所示。

表1 解纏繞FPGA實(shí)現(xiàn)偽代碼Table 1 Unwrapping pseudo code

2.3 基于最小二乘的曲線擬合

如前所述，為了降低計(jì)算量，多普勒中心頻率估計(jì)是1/4抽取后進(jìn)行處理的，但是在方位向壓縮中，需要計(jì)算每個距離門上的多普勒頻率信息，為了得到完整的多普勒頻率估計(jì)值，需要對多普勒頻率的估計(jì)值進(jìn)行插值。本文采用了二次曲線擬合的方法來實(shí)現(xiàn)缺失部分多普勒中心頻率的插值?；谧钚《朔ǖ那€擬合方法提供基于均方誤差最小的最優(yōu)解，可以有效地降低多普勒中心頻率估計(jì)誤差，提高估計(jì)精度。也能在一定程度上抑制由于平臺抖動等帶來的多普勒中心頻率劇烈變化，對多普勒中心頻率估計(jì)結(jié)果具有濾波的作用。多普勒中心頻率在標(biāo)準(zhǔn)掃描SAR場景中可以用二次曲線對多普勒中心頻率進(jìn)行近似描述。

下面詳細(xì)介紹基于最小二乘的二次曲線擬合方法中矩陣求逆的FPGA簡化實(shí)現(xiàn)方法[22]。對于線性函數(shù)y=Ax，若A,y為已知量，x為待解量，則x可以表示為：x=inv(ATA)ATy。其中，由于矩陣A可能不是方陣，通過左乘AT得到方陣ATA。

對于二次曲線：

(5)

其矩陣A：

(6)

x=(a,b,c).

對于3×3矩陣，矩陣求逆操作可以通過代數(shù)余子式來實(shí)現(xiàn)[23]：

(7)

定義矩陣M=ATA，亦即

(8)

(9)

|A|=a11|A11|+a12|A12|+a13|A13|，

(10)

(11)

式中：

(12)

在FPGA實(shí)現(xiàn)過程中，矩陣M的行列式可以通過其代數(shù)余子式進(jìn)行簡化實(shí)現(xiàn)：

|M|=a11·detA11+a12·detA12+a13·detA13=

fx4·detA11+fx3·detA12+fx2·detA13.

(13)

在得到矩陣M的代數(shù)余子式和行列式后，即可以通過cordic模塊實(shí)現(xiàn)除法操作，實(shí)現(xiàn)矩陣的求逆，從而可以求解方陣ATA的逆矩陣，從而可以實(shí)現(xiàn)二次曲線系數(shù)(a，b，c)。再得到二次曲線的系數(shù)后，可以通過動態(tài)生成x及其相關(guān)x2，從而完成二次曲線的插值。

綜上所述，在FPGA實(shí)現(xiàn)基于最小二乘的二次曲線擬合中起主要關(guān)鍵點(diǎn)在于：①動態(tài)生成采用流水線方式生成x通過乘法器的流水線操作從x中得到x2，隨后通過流水線方式直接輸出fx,fx2,fx3,fx4。②通過代數(shù)余子式的形式求解方陣ATA的逆矩陣inv(ATA)，從而為使用FPGA實(shí)現(xiàn)最小二乘算法提供有力保證。

3 實(shí)驗(yàn)結(jié)果

基于FPGA的多普勒中心頻率估計(jì)為彈載SAR純FPGA實(shí)現(xiàn)的一部分，已經(jīng)在某項(xiàng)目上實(shí)現(xiàn)，并完成了各項(xiàng)測試，所有指標(biāo)均達(dá)到或優(yōu)于相關(guān)要求。SAR成像模塊采用一片Xilinx公司的Virtex 7系列xc7vx690t芯片外擴(kuò)2片DDR3芯片，完成多普勒中心頻率估計(jì)、RCMC、多普勒調(diào)頻率估計(jì)以及方位向壓縮，生成4 096×2 048大小的圖像。系統(tǒng)信號處理模塊主頻為200 MHz，系統(tǒng)LE資源占用約占芯片總資源的58%。信號處理板數(shù)量由原理樣機(jī)的4塊減少為2塊，處理功耗下降85 W以上。FDC模塊資源使用情況如表2所示。其中arcTan函數(shù)為采用cordic完成相關(guān)數(shù)據(jù)的反正切計(jì)算；CorrComplex函數(shù)計(jì)算復(fù)數(shù)據(jù)相關(guān)；corrMemRe函數(shù)為使用LUT資源例化相關(guān)累積。lsqFDC函數(shù)使用最小二乘法實(shí)現(xiàn)二次曲線擬合系數(shù)估計(jì)。genFDC函數(shù)使用lsqFDC計(jì)算的系數(shù)，實(shí)時生成2 048點(diǎn)有效場景的FDC值；thetaWarp函數(shù)完成相位角解纏繞。FDC估計(jì)模塊在RC數(shù)據(jù)接收過程中直接流水完成，不消耗額外的運(yùn)行時間，接收完數(shù)據(jù)可以直接得到相關(guān)累積結(jié)果。為了使參數(shù)擬合和插值模塊不產(chǎn)生等待，F(xiàn)DC相關(guān)累積模塊拋棄最后一個方位向累積，直接轉(zhuǎn)入多普勒參數(shù)擬合狀態(tài)，因基于最小二乘的二次曲線擬合和插值模塊以及后續(xù)的多普勒模糊數(shù)疊加模塊總處理時間為4 124個時鐘周期，在200 MHz處理時鐘下，為2.1 μs，遠(yuǎn)小于重頻周期。因而可以在RC數(shù)據(jù)接收完最后一包距離向脈壓數(shù)據(jù)前直接輸出FDC估計(jì)結(jié)果，直接轉(zhuǎn)入RCMC處理，實(shí)現(xiàn)系統(tǒng)“0等待”，因此FDC估計(jì)模塊系統(tǒng)消耗時間為0。

表2 FDC模塊資源使用統(tǒng)計(jì)Table 2 FDC module resource usage statistics

4 結(jié)束語

本文根據(jù)彈載SAR對信號處理軟硬件實(shí)現(xiàn)的需求，提出了一種基于FPGA的快速多普勒中心頻率估計(jì)方法，采用流水線并行處理的模式，實(shí)現(xiàn)FDC估計(jì)時間的零消耗。為了實(shí)現(xiàn)基于FPGA的FDC估計(jì)，創(chuàng)新性地提出了基于FPGA的流水線式相位解纏繞方法和基于最小二乘的參數(shù)曲線擬合算法，實(shí)現(xiàn)相位解纏繞的連續(xù)流水輸出以及中心頻率估計(jì)的二次取下擬合和插值，并在實(shí)際工程中進(jìn)行了有效驗(yàn)證。測試結(jié)果說明，該方法充分發(fā)揮了FPGA并行處理的優(yōu)勢，大大降低了SAR信號處理對硬件資源的需求，在彈載SAR和無人機(jī)載SAR等對結(jié)構(gòu)、質(zhì)量和功耗限制較嚴(yán)格的條件下，有著非常廣闊的應(yīng)用前景。