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

羅正華， 周方均， 雷 林， 李 霞， 劉一達

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

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

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

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

1 時差提取算法

研究中，時差提取模塊需要處理來自4個基站接收到的無人機信號，提取出信號到達輔站和主站的時間差。信號接收系統(tǒng)以主站為節(jié)點，主站與輔站站距由近到遠，分別為70、110、150 m的Y字型布設，如圖1所示。

1.1 廣義互相關算法

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

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

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

圖2 廣義互相關算法流程Fig.2 Generalized cross-correlation algorithm flow

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

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

圖3 卡爾曼-最優(yōu)階互相關算法流程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(·)為取絕對值；max(·)為取信號中的最大值；C1(t)、C2(t)分別為對信號X1(t)和X2(t)進行幅值歸一化。

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

(3)

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

(4)

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

(5)

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

(6)

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

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

(7)

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

(8)

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

(9)

(10)

(11)

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

令

?1=δn-δn-1

(12)

?2=δn-δn+1

(13)

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

(14)

(15)

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

mnew.i=mi-1

(16)

式(16)將上一時刻的最優(yōu)估計值作為當前時刻的時差預測值。

Pnew.i=Pi-1+Q

(17)

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

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

(18)

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

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

(19)

式(19)結合當前時刻時差的提取值，對上一時刻的預測進行校正，并輸出最優(yōu)時差提取值。

Pi=(1-Ki)Pnew.i

(20)

式(20)對最優(yōu)時差估計值的方差進行更新。

2 算法仿真

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

2.1 廣義互相關算法仿真分析

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

圖4 未加噪聲的信號波形Fig.4 Signal waveform without noise

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

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

仿真時選取時差值為122個點的兩路信號。傅里葉變換的階數(shù)以21為公比，從21變換到220。實驗結果如圖6所示。

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

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

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

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

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

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

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

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

圖10 誤差分析Fig.10 Error analysis

2.3 5種算法仿真對比分析

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

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

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

3 時差提取系統(tǒng)的設計

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

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

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

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

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

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

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

圖15 不同平臺的對比Fig.15 Comparison of different platforms

4 結論

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