一種用于時(shí)差提取的卡爾曼-最優(yōu)階互相關(guān)算法

羅正華， 周方均， 雷 林， 李 霞， 劉一達(dá)

(1.成都大學(xué)信息科學(xué)與工程學(xué)院， 成都 610106； 2.電信科學(xué)技術(shù)第五研究所， 成都 610062)

時(shí)差提取是指對(duì)同一信號(hào)到達(dá)不同基站的時(shí)間差進(jìn)行提取的一種技術(shù)[1]。在時(shí)差提取技術(shù)中，基本互相關(guān)算法是一種最基本的算法，其算法原理簡(jiǎn)單，計(jì)算復(fù)雜程度低，計(jì)算量較小。另外，廣義互相關(guān)(generalized cross correlation)算法[2]通過(guò)對(duì)信號(hào)進(jìn)行加權(quán)濾波處理來(lái)提高信噪比，其主要的加權(quán)濾波器有ROTH處理器、平滑相干變換(smoothed coherence transform,SCOT)[3]以及相位變換(phase transform,PHAT)[4]。

文獻(xiàn)[5]說(shuō)明了最小均方(least mean square,LMS)算法在噪聲信號(hào)處理中的優(yōu)越性。然而，它的實(shí)現(xiàn)需要選擇一個(gè)收斂參數(shù)。由于選擇附加參數(shù)的必要性，使得LMS算法的實(shí)用性不如基本互相關(guān)算法。文獻(xiàn)[6]提出了一種基于分?jǐn)?shù)階傅里葉變換(Fourier transform,FFT)的廣義互相關(guān)(generialized cross-correlation,GCC)算法，但是廣義互相關(guān)算法需要合理地選擇加權(quán)函數(shù)才能取得更好的濾波效果和更高的時(shí)差提取精度。文獻(xiàn)[7]通過(guò)現(xiàn)場(chǎng)可編程邏輯門(mén)陣列(field programmable gate array，F(xiàn)PGA)實(shí)現(xiàn)了一種基于基本互相關(guān)算法的時(shí)域數(shù)字濾波算法，但是其計(jì)算互相關(guān)函數(shù)序列(n=N-99～N+99共199點(diǎn))采樣點(diǎn)數(shù)太少，不適用于數(shù)據(jù)吞吐量大的無(wú)人機(jī)定位。文獻(xiàn)[8]在FPGA上實(shí)現(xiàn)了正交頻分復(fù)用(orthogonal frequency division multiplexing，OFDM)的接收技術(shù)并提出了128點(diǎn)IFFT/FFT(inverse fast flourier transform/fast flourier transform)流水線結(jié)構(gòu)，但是其FFT變換的蝶形運(yùn)算階數(shù)固定，無(wú)法滿足階數(shù)的實(shí)時(shí)更新。更重要的是，在利用FPGA實(shí)現(xiàn)傅里葉變換功能時(shí)，變換階數(shù)的不合理選擇會(huì)造成大量的資源浪費(fèi)，同時(shí)，在兩路信號(hào)的時(shí)差值與傅里葉變換階數(shù)不耦合的情況下，將導(dǎo)致時(shí)差值提取出錯(cuò)。文獻(xiàn)[9]提出了一種利用卡爾曼濾波器對(duì)測(cè)量數(shù)據(jù)進(jìn)行預(yù)處理的定位技術(shù)，但其并沒(méi)有做到卡爾曼自適應(yīng)濾波。

針對(duì)上述問(wèn)題，提出卡爾曼-最優(yōu)階互相關(guān)算法，根據(jù)不同的實(shí)際時(shí)差值，選擇與之對(duì)應(yīng)的最優(yōu)的傅里葉變換階數(shù)，并進(jìn)行卡爾曼濾波，既減少了對(duì)FPGA資源的占用，也避免了因傅里葉變換階數(shù)導(dǎo)致的時(shí)差提取出錯(cuò)，同時(shí)，可以提高低信噪比條件下的時(shí)差提取精度。

1 時(shí)差提取算法

研究中，時(shí)差提取模塊需要處理來(lái)自4個(gè)基站接收到的無(wú)人機(jī)信號(hào)，提取出信號(hào)到達(dá)輔站和主站的時(shí)間差。信號(hào)接收系統(tǒng)以主站為節(jié)點(diǎn)，主站與輔站站距由近到遠(yuǎn)，分別為70、110、150 m的Y字型布設(shè)，如圖1所示。

1.1 廣義互相關(guān)算法

廣義互相關(guān)算法是先進(jìn)行廣義加權(quán)，再將加權(quán)濾波器輸出的信號(hào)進(jìn)行互相關(guān)。研究中使用的廣義互相關(guān)加權(quán)[10]函數(shù)如表1所示。算法流程如圖2所示。

在表1中，ξ(f)為廣義互相關(guān)頻域的加權(quán)函數(shù)；Φ11(f)、Φ22(f)分別表示信號(hào)X1(t)和X2(t)的自相關(guān)函數(shù)。當(dāng)ξ(f)=1時(shí)，即為基本互相關(guān)算法。

表1 廣義互相關(guān)加權(quán)函數(shù)Table 1 Generalized cross-correlation weighting function

圖2 廣義互相關(guān)算法流程Fig.2 Generalized cross-correlation algorithm flow

1.2 卡爾曼-最優(yōu)階互相關(guān)算法

卡爾曼-最優(yōu)階互相關(guān)算法是將自適應(yīng)獲取最優(yōu)傅里葉變換階數(shù)[11]與卡爾曼濾波兩者進(jìn)行結(jié)合的算法，算法流程如圖3所示。

圖3 卡爾曼-最優(yōu)階互相關(guān)算法流程Fig.3 Kalman-optimal cross-correlation algorithm flow

遞推公式為

C1(t)=X1(t)/max{abs[X1(t)]}

(1)

C2(t)=X2(t)/max{abs[X2(t)]}

(2)

式中：abs(·)為取絕對(duì)值；max(·)為取信號(hào)中的最大值；C1(t)、C2(t)分別為對(duì)信號(hào)X1(t)和X2(t)進(jìn)行幅值歸一化。

R1(ω)=FFT[C1(t),2N]

(3)

R2(ω)=FFT[C2(t),2N]

(4)

R11(ω)=R1(ω)conj[R1(ω)]

(5)

R22(ω)=R2(ω)conj[R2(ω)]

(6)

式中：R11(ω)、R22(ω)為信號(hào)C1(t)、C2(t)的頻域自相關(guān)函數(shù)；N為傅里葉變換階數(shù)，conj(·)為取共軛復(fù)數(shù)。

G12(ω)=R1(ω)conj[R2(ω)]

(7)

K1(ω)=R11(ω)conj[R22(ω)]

(8)

Y12(t)=IFFT{G12(ω)conj[K1(ω)]}

(9)

(10)

(11)

式(11)中：δn為當(dāng)前時(shí)刻對(duì)應(yīng)階數(shù)的時(shí)差提取值，?n為[2(n-1)-1]e|mi-mi-1|由實(shí)驗(yàn)統(tǒng)計(jì)得到的修正參數(shù)，該值與此時(shí)刻的時(shí)差值和傅里葉變換階數(shù)有關(guān)。

令

?1=δn-δn-1

(12)

?2=δn-δn+1

(13)

當(dāng)?1≠0且?2=0時(shí)，n為最佳傅里葉變換階數(shù)。

(14)

(15)

式(15)中：將前4次最優(yōu)階數(shù)的時(shí)差估計(jì)值的平均值作為卡爾曼的初始時(shí)差提取值m1。

mnew.i=mi-1

(16)

式(16)將上一時(shí)刻的最優(yōu)估計(jì)值作為當(dāng)前時(shí)刻的時(shí)差預(yù)測(cè)值。

Pnew.i=Pi-1+Q

(17)

式(17)中：Pi-1為上一時(shí)刻最優(yōu)時(shí)差估計(jì)值的方差；Q表示連續(xù)兩個(gè)時(shí)刻最優(yōu)時(shí)差估計(jì)值的方差。

Ki=Pnew.i/(Pnew.i+R)

(18)

式(18)中：R表示兩次測(cè)量時(shí)差值間的方差；Ki為卡爾曼的計(jì)算增益。

mi=mnew.i+Ki(δn-mnew.i)

(19)

式(19)結(jié)合當(dāng)前時(shí)刻時(shí)差的提取值，對(duì)上一時(shí)刻的預(yù)測(cè)進(jìn)行校正，并輸出最優(yōu)時(shí)差提取值。

Pi=(1-Ki)Pnew.i

(20)

式(20)對(duì)最優(yōu)時(shí)差估計(jì)值的方差進(jìn)行更新。

2 算法仿真

實(shí)驗(yàn)采用MATLAB進(jìn)行仿真，仿真條件為：采集的大疆精靈3圖傳信號(hào)，采樣頻率Fs=200 MHz，采樣周期Ts=5 ns。在仿真中，兩段信號(hào)的時(shí)差為122Ts。信號(hào)波形如圖4所示。

2.1 廣義互相關(guān)算法仿真分析

研究中，廣義互相關(guān)算法仿真用于比較不同加權(quán)函數(shù)的時(shí)差提取精度。各種算法的均方根誤差如圖5所示。

圖4 未加噪聲的信號(hào)波形Fig.4 Signal waveform without noise

圖5 各種算法的均方根誤差Fig.5 Root mean sqare error of various algorithms

2.2 卡爾曼-最優(yōu)階互相關(guān)算法仿真分析

仿真時(shí)選取時(shí)差值為122個(gè)點(diǎn)的兩路信號(hào)。傅里葉變換的階數(shù)以21為公比，從21變換到220。實(shí)驗(yàn)結(jié)果如圖6所示。

圖6的結(jié)果表明，時(shí)差提取值隨傅里葉變換階數(shù)的增大而增大。對(duì)圖6的曲線進(jìn)行數(shù)據(jù)統(tǒng)計(jì)與分析，采用數(shù)學(xué)歸納法對(duì)曲線進(jìn)行修正，將階數(shù)與時(shí)差提取值的對(duì)應(yīng)關(guān)系統(tǒng)一起來(lái)，結(jié)果如圖7所示。

由圖7發(fā)現(xiàn)，在時(shí)差值為122的條件下，當(dāng)傅里葉變換的階數(shù)為29以后，時(shí)差提取值保持穩(wěn)定的輸出結(jié)果，故9階為最優(yōu)傅里葉變換階數(shù)。在此基礎(chǔ)上，連續(xù)改變兩路信號(hào)的時(shí)差值時(shí)，并進(jìn)行二階曲線擬合，最優(yōu)階數(shù)的變換如圖8所示。

在圖8中，最優(yōu)階數(shù)隨著時(shí)差值的增加而增大，說(shuō)明兩路信號(hào)的時(shí)差值與傅里葉變換階數(shù)存在關(guān)系。將兩值的關(guān)系統(tǒng)一起來(lái)，得到一個(gè)與時(shí)差變化相關(guān)的量?n，?n滿足任一時(shí)差與最優(yōu)階數(shù)的關(guān)系。將最優(yōu)階得到的時(shí)差進(jìn)行卡爾曼濾波，結(jié)果如圖9所示。

圖6 階數(shù)變換估計(jì)結(jié)果Fig.6 Order transformation estimation result

圖7 參數(shù)修正結(jié)果Fig.7 Parameter correction results

圖8 最優(yōu)階數(shù)隨時(shí)差值的變化結(jié)果Fig.8 The change result of the optimal order with time difference

在圖9中，在50 μs內(nèi)，卡爾曼-最優(yōu)階互相關(guān)算法具有較好的收斂值，可以更準(zhǔn)確地估計(jì)出實(shí)際的時(shí)差值，并提高算法在不同信噪比條件下的時(shí)差提取精度。誤差分析如圖10所示。

圖9 卡爾曼濾波結(jié)果Fig.9 The Kalman filter’s results

圖10 誤差分析Fig.10 Error analysis

2.3 5種算法仿真對(duì)比分析

對(duì)算法進(jìn)行曲線擬合[12]，比較5種算法的抗噪能力、時(shí)差提取精度和收斂速度。如圖11所示。

圖11 5種相關(guān)算法曲線擬合性能比較Fig.11 Comparison of curve fitting performance of five related algorithms

通過(guò)圖11對(duì)比分析可以發(fā)現(xiàn)，隨著信噪比的降低，相對(duì)于其他4種算法，卡爾曼-最優(yōu)階互相關(guān)算法無(wú)論是在抗躁能力和時(shí)差提取精度上，還是在收斂速度上，均有明顯的優(yōu)勢(shì)。其中，ROTH加權(quán)函數(shù)的誤差最大，說(shuō)明它對(duì)高斯白噪聲濾波效果并不理想?？梢?jiàn)，合適的加權(quán)函數(shù)是提高時(shí)差提取精度的重要因素。之后，將會(huì)加大卡爾曼-最優(yōu)階互相關(guān)算法在采樣誤差[13]、相關(guān)性信號(hào)[14]、多路徑信號(hào)[15]、內(nèi)插函數(shù)的分?jǐn)?shù)倍估計(jì)[16]領(lǐng)域內(nèi)的研究。

3 時(shí)差提取系統(tǒng)的設(shè)計(jì)

系統(tǒng)包含4單元天線陣列和射頻接收芯片AD9361以及基于含有ARM Cortex-a9的Xilinx Zynq-7000全可編程片上系統(tǒng)的嵌入式處理器。

實(shí)驗(yàn)第一階段為現(xiàn)場(chǎng)測(cè)試。天氣狀況：溫度為16 ℃，2級(jí)東北風(fēng)，伴有雨霧；環(huán)境：存在WIFI信號(hào)干擾；風(fēng)級(jí)：<2級(jí)；主要設(shè)備：大疆PHANTOM 4型專業(yè)版、自研架構(gòu)的全向高增益天線等?，F(xiàn)場(chǎng)如圖12所示。

圖12 四川天空之眼現(xiàn)場(chǎng)測(cè)試Fig.12 Sichuan Sky Eye field test

系統(tǒng)時(shí)差提取部分硬件原理圖如圖13所示。

圖13 系統(tǒng)部分硬件原理圖Fig.13 Partial hardware schematic of system

在圖13中，F(xiàn)PGA可編程模塊讀取AD9361中的數(shù)據(jù)，進(jìn)行時(shí)差提取。ARM處理器通過(guò)GPIO(general-purpose input/output)讀取時(shí)差值進(jìn)行定位解算，并通過(guò)串口發(fā)送至上位機(jī)。上位機(jī)將解算坐標(biāo)信息與真實(shí)坐標(biāo)進(jìn)行對(duì)比，可驗(yàn)證卡爾曼-最優(yōu)階互相關(guān)算法的精度。并且，在相同數(shù)據(jù)吞吐量條件下，可以對(duì)比不同處理平臺(tái)的數(shù)據(jù)處理速度。上位機(jī)實(shí)驗(yàn)結(jié)果顯示如圖14所示。不同平臺(tái)的對(duì)比結(jié)果如圖15所示。

圖14的對(duì)比實(shí)驗(yàn)結(jié)果表明，使用卡爾曼-最優(yōu)階互相關(guān)算法解算出的坐標(biāo)(-12.654 5 m,88.006 5 m, 249.003 9 m)與真實(shí)坐標(biāo)(-12.5 m,87.5 m,249.3 m)的誤差相對(duì)較小，此時(shí)無(wú)人機(jī)處于主站，東偏北81.817 76°，距離為264.401 m，高度為249.003 9 m。由圖15的實(shí)驗(yàn)結(jié)果可以分析在不同平臺(tái)的處理速度，其中，MATLAB的處理時(shí)間為254 000 μs，C#的處理時(shí)間為40 836.2 μs，而FPGA的處理時(shí)間為579 μs，可知在處理相同數(shù)據(jù)時(shí)，F(xiàn)PGA具有明顯的速度優(yōu)勢(shì)。因此，設(shè)計(jì)的硬件系統(tǒng)與相關(guān)算法相結(jié)合后在實(shí)際工程應(yīng)用上具有明顯的優(yōu)勢(shì)。

圖15 不同平臺(tái)的對(duì)比Fig.15 Comparison of different platforms

4 結(jié)論

從算法原理和仿真角度，對(duì)基本互相關(guān)、廣義互相關(guān)等傳統(tǒng)時(shí)差提取算法進(jìn)行了探討，改進(jìn)并提出了卡爾曼-最優(yōu)階互相關(guān)優(yōu)化算法?？柭?最優(yōu)階互相關(guān)算法的優(yōu)點(diǎn)是通過(guò)對(duì)傅里葉變換階數(shù)的掃描，可以使系統(tǒng)的硬件資源得以合理分配，同時(shí)，可以解決因傅里葉變換階數(shù)導(dǎo)致的時(shí)差提取出錯(cuò)的問(wèn)題。此外，卡爾曼濾波加快了最優(yōu)時(shí)差提取的收斂速度，提升了時(shí)差提取的精度。對(duì)于無(wú)人機(jī)圖傳信號(hào)(OFDM信號(hào))，卡爾曼-最優(yōu)階互相關(guān)算法展現(xiàn)出了更強(qiáng)的資源分配能力和抗噪能力，為進(jìn)一步研究采樣誤差、相關(guān)性信號(hào)、多路徑信號(hào)、內(nèi)插函數(shù)的分?jǐn)?shù)倍估計(jì)對(duì)時(shí)差提取精度的影響提供了參考，也為進(jìn)一步的工程應(yīng)用奠定了基礎(chǔ)。