一种用于时差提取的卡尔曼-最优阶互相关算法

罗正华， 周方均， 雷 林， 李 霞， 刘一达

(1.成都大学信息科学与工程学院， 成都 610106； 2.电信科学技术第五研究所， 成都 610062)

时差提取是指对同一信号到达不同基站的时间差进行提取的一种技术[1]。在时差提取技术中，基本互相关算法是一种最基本的算法，其算法原理简单，计算复杂程度低，计算量较小。另外，广义互相关(generalized cross correlation)算法[2]通过对信号进行加权滤波处理来提高信噪比，其主要的加权滤波器有ROTH处理器、平滑相干变换(smoothed coherence transform,SCOT)[3]以及相位变换(phase transform,PHAT)[4]。

文献[5]说明了最小均方(least mean square,LMS)算法在噪声信号处理中的优越性。然而，它的实现需要选择一个收敛参数。由于选择附加参数的必要性，使得LMS算法的实用性不如基本互相关算法。文献[6]提出了一种基于分数阶傅里叶变换(Fourier transform,FFT)的广义互相关(generialized cross-correlation,GCC)算法，但是广义互相关算法需要合理地选择加权函数才能取得更好的滤波效果和更高的时差提取精度。文献[7]通过现场可编程逻辑门阵列(field programmable gate array，FPGA)实现了一种基于基本互相关算法的时域数字滤波算法，但是其计算互相关函数序列(n=N-99～N+99共199点)采样点数太少，不适用于数据吞吐量大的无人机定位。文献[8]在FPGA上实现了正交频分复用(orthogonal frequency division multiplexing，OFDM)的接收技术并提出了128点IFFT/FFT(inverse fast flourier transform/fast flourier transform)流水线结构，但是其FFT变换的蝶形运算阶数固定，无法满足阶数的实时更新。更重要的是，在利用FPGA实现傅里叶变换功能时，变换阶数的不合理选择会造成大量的资源浪费，同时，在两路信号的时差值与傅里叶变换阶数不耦合的情况下，将导致时差值提取出错。文献[9]提出了一种利用卡尔曼滤波器对测量数据进行预处理的定位技术，但其并没有做到卡尔曼自适应滤波。

针对上述问题，提出卡尔曼-最优阶互相关算法，根据不同的实际时差值，选择与之对应的最优的傅里叶变换阶数，并进行卡尔曼滤波，既减少了对FPGA资源的占用，也避免了因傅里叶变换阶数导致的时差提取出错，同时，可以提高低信噪比条件下的时差提取精度。

1 时差提取算法

研究中，时差提取模块需要处理来自4个基站接收到的无人机信号，提取出信号到达辅站和主站的时间差。信号接收系统以主站为节点，主站与辅站站距由近到远，分别为70、110、150 m的Y字型布设，如图1所示。

1.1 广义互相关算法

广义互相关算法是先进行广义加权，再将加权滤波器输出的信号进行互相关。研究中使用的广义互相关加权[10]函数如表1所示。算法流程如图2所示。

在表1中，ξ(f)为广义互相关频域的加权函数；Φ11(f)、Φ22(f)分别表示信号X1(t)和X2(t)的自相关函数。当ξ(f)=1时，即为基本互相关算法。

表1 广义互相关加权函数Table 1 Generalized cross-correlation weighting function

图2 广义互相关算法流程Fig.2 Generalized cross-correlation algorithm flow

1.2 卡尔曼-最优阶互相关算法

卡尔曼-最优阶互相关算法是将自适应获取最优傅里叶变换阶数[11]与卡尔曼滤波两者进行结合的算法，算法流程如图3所示。

图3 卡尔曼-最优阶互相关算法流程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)的频域自相关函数；N为傅里叶变换阶数，conj(·)为取共轭复数。

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

(7)

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

(8)

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

(9)

(10)

(11)

式(11)中：δn为当前时刻对应阶数的时差提取值，∂n为[2(n-1)-1]e|mi-mi-1|由实验统计得到的修正参数，该值与此时刻的时差值和傅里叶变换阶数有关。

令

ϑ1=δn-δn-1

(12)

ϑ2=δn-δn+1

(13)

当ϑ1≠0且ϑ2=0时，n为最佳傅里叶变换阶数。

(14)

(15)

式(15)中：将前4次最优阶数的时差估计值的平均值作为卡尔曼的初始时差提取值m1。

mnew.i=mi-1

(16)

式(16)将上一时刻的最优估计值作为当前时刻的时差预测值。

Pnew.i=Pi-1+Q

(17)

式(17)中：Pi-1为上一时刻最优时差估计值的方差；Q表示连续两个时刻最优时差估计值的方差。

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

(18)

式(18)中：R表示两次测量时差值间的方差；Ki为卡尔曼的计算增益。

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

(19)

式(19)结合当前时刻时差的提取值，对上一时刻的预测进行校正，并输出最优时差提取值。

Pi=(1-Ki)Pnew.i

(20)

式(20)对最优时差估计值的方差进行更新。

2 算法仿真

实验采用MATLAB进行仿真，仿真条件为：采集的大疆精灵3图传信号，采样频率Fs=200 MHz，采样周期Ts=5 ns。在仿真中，两段信号的时差为122Ts。信号波形如图4所示。

2.1 广义互相关算法仿真分析

研究中，广义互相关算法仿真用于比较不同加权函数的时差提取精度。各种算法的均方根误差如图5所示。

图4 未加噪声的信号波形Fig.4 Signal waveform without noise

图5 各种算法的均方根误差Fig.5 Root mean sqare error of various algorithms

2.2 卡尔曼-最优阶互相关算法仿真分析

仿真时选取时差值为122个点的两路信号。傅里叶变换的阶数以21为公比，从21变换到220。实验结果如图6所示。

图6的结果表明，时差提取值随傅里叶变换阶数的增大而增大。对图6的曲线进行数据统计与分析，采用数学归纳法对曲线进行修正，将阶数与时差提取值的对应关系统一起来，结果如图7所示。

由图7发现，在时差值为122的条件下，当傅里叶变换的阶数为29以后，时差提取值保持稳定的输出结果，故9阶为最优傅里叶变换阶数。在此基础上，连续改变两路信号的时差值时，并进行二阶曲线拟合，最优阶数的变换如图8所示。

在图8中，最优阶数随着时差值的增加而增大，说明两路信号的时差值与傅里叶变换阶数存在关系。将两值的关系统一起来，得到一个与时差变化相关的量∂n，∂n满足任一时差与最优阶数的关系。将最优阶得到的时差进行卡尔曼滤波，结果如图9所示。

图6 阶数变换估计结果Fig.6 Order transformation estimation result

图7 参数修正结果Fig.7 Parameter correction results

图8 最优阶数随时差值的变化结果Fig.8 The change result of the optimal order with time difference

在图9中，在50 μs内，卡尔曼-最优阶互相关算法具有较好的收敛值，可以更准确地估计出实际的时差值，并提高算法在不同信噪比条件下的时差提取精度。误差分析如图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对比分析可以发现，随着信噪比的降低，相对于其他4种算法，卡尔曼-最优阶互相关算法无论是在抗躁能力和时差提取精度上，还是在收敛速度上，均有明显的优势。其中，ROTH加权函数的误差最大，说明它对高斯白噪声滤波效果并不理想。可见，合适的加权函数是提高时差提取精度的重要因素。之后，将会加大卡尔曼-最优阶互相关算法在采样误差[13]、相关性信号[14]、多路径信号[15]、内插函数的分数倍估计[16]领域内的研究。

3 时差提取系统的设计

系统包含4单元天线阵列和射频接收芯片AD9361以及基于含有ARM Cortex-a9的Xilinx Zynq-7000全可编程片上系统的嵌入式处理器。

实验第一阶段为现场测试。天气状况：温度为16 ℃，2级东北风，伴有雨雾；环境：存在WIFI信号干扰；风级：<2级；主要设备：大疆PHANTOM 4型专业版、自研架构的全向高增益天线等。现场如图12所示。

图12 四川天空之眼现场测试Fig.12 Sichuan Sky Eye field test

系统时差提取部分硬件原理图如图13所示。

图13 系统部分硬件原理图Fig.13 Partial hardware schematic of system

在图13中，FPGA可编程模块读取AD9361中的数据，进行时差提取。ARM处理器通过GPIO(general-purpose input/output)读取时差值进行定位解算，并通过串口发送至上位机。上位机将解算坐标信息与真实坐标进行对比，可验证卡尔曼-最优阶互相关算法的精度。并且，在相同数据吞吐量条件下，可以对比不同处理平台的数据处理速度。上位机实验结果显示如图14所示。不同平台的对比结果如图15所示。

图14的对比实验结果表明，使用卡尔曼-最优阶互相关算法解算出的坐标(-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，可知在处理相同数据时，FPGA具有明显的速度优势。因此，设计的硬件系统与相关算法相结合后在实际工程应用上具有明显的优势。

图15 不同平台的对比Fig.15 Comparison of different platforms

4 结论

从算法原理和仿真角度，对基本互相关、广义互相关等传统时差提取算法进行了探讨，改进并提出了卡尔曼-最优阶互相关优化算法。卡尔曼-最优阶互相关算法的优点是通过对傅里叶变换阶数的扫描，可以使系统的硬件资源得以合理分配，同时，可以解决因傅里叶变换阶数导致的时差提取出错的问题。此外，卡尔曼滤波加快了最优时差提取的收敛速度，提升了时差提取的精度。对于无人机图传信号(OFDM信号)，卡尔曼-最优阶互相关算法展现出了更强的资源分配能力和抗噪能力，为进一步研究采样误差、相关性信号、多路径信号、内插函数的分数倍估计对时差提取精度的影响提供了参考，也为进一步的工程应用奠定了基础。