CN104635231A - 基于自适应卡尔曼滤波的机车雷达测速方法 - Google Patents

Abstract

本发明属于雷达近场测量技术领域，特别涉及基于自适应卡尔曼滤波的机车雷达测速方法，其具体步骤为：(1)计算回波的功率谱；(2)抑制低频干扰；(3)利用频谱重心法预估机车多普勒频率；(4)进行雷达位置自校准与速度解算；(5)利用自适应卡尔曼滤波器进行滤波；(6)机车速度输出。本发明可实现机车速度的实时测量，测量速度无滞后，测量精度高。

Description

基于自适应卡尔曼滤波的机车雷达测速方法

技术领域

本发明属于雷达近场测量技术领域，具体涉及机车速度的测量方法，特别涉及基于自适应卡尔曼滤波的机车雷达测速方法，可实现复杂电磁环境和不同的铁路路基环境下对机车速度稳定精确的测量。

背景技术

速度信息对于确保机车的运行安全具有重要意义，车载雷达测速传感器在高速重载机车上有着广泛的应用前景。在上个世纪七十年代，德国、法国、美国等已开始研发机车及汽车测速雷达。经过多年的发展和改进，在上个世纪末，测速雷达已开始在高速铁路上商业运营，在公路上的应用尤为广泛。其中德国在微型雷达前端的研发生产方面在世界上占据领先地位。目前，在我国的铁路系统中，机车速度的测量主要是依靠光电传感器，这种光电测速设备通过测量机车车轮或者一个刚性轴的转数来获得机车的速度值。光电传感器的原理简单，容易实现，但是这种测速方法在车轮打滑、空转、以及车轮直径因磨损而发生改变的情况下都会出现一定的系统偏差，不能满足高精度和高可靠性的要求。另外，由于磁悬浮列车没有车轮，因此，传统的光电测速设备无法测量磁悬浮机车的速度。正是由于这些问题的存在，雷达微波传感技术以非接触方式对机车速度测量成为了一种较好的选择。

雷达测速是利用多普勒效应来测量车辆对地的相对速度。多普勒效应是指当发射源和接收者之间有相对径向运动时，接收到的信号频率将发生变化。通过分析所获得的频率偏移，根据多普勒雷达公式即可计算出车辆的实际运动速度。所以雷达测出的速度与机车轮径无关，与车轮旋转的状态也无关，同时它是一种非接触式的传感器，无机械磨损，尤其适宜高速、重载机车和无旋转车轮的磁悬浮列车。如果把雷达测速结果与光电测速结果进行融合，就可以在现有国产监控设备不作任何改动的条件下，提供更为可靠的速度信号，从而实现机车的精确定位、相位防溜、防滑、防空转等功能。

虽然理论上多普勒测速雷达能够克服传统方式的缺点，稳定可靠，并能提供精确的测量精度，但实际上由于受机车电磁环境、地理环境、雷达天线波束宽度和雷达天线安装角度误差等因素的影响，现有的多普勒机车测速雷达的测量精度并不高，而且容易受到雷达本身安装误差的影响，同时也不能克服不同路基环境所造成的测速误差。传统的滑动平均方法将一段滑窗内的速度平均值作为最终的估计结果，在一定程度上减小了测速误差的波动，但是当机车加速或者减速时，容易出现估计速度滞后的现象。

发明内容

本发明的目的在于针对上述机车雷达测速方法的不足，提出了一种基于自适应卡尔曼滤波的测速方法，以避免现有的机车雷达测速精度容易受到机车电磁环境、不同路基环境、波束天线宽度和雷达安装角度误差等诸多因素的影响，避免出现估计速度滞后的现象，提高测速精度。

为实现上述技术目的，本发明采用如下技术方案予以实现。

基于自适应卡尔曼滤波的机车雷达测速方法包括以下步骤：

步骤1，在机车静止时，在机车底部安装测速雷达；测速雷达具有两个天线，包括位于测速雷达前部的前天线、以及位于测速雷达后部的后天线，前天线朝向机车的前方，后天线朝向机车的后方；

步骤2，当机车开始运行后，测速雷达向外发射信号，在k时刻前天线和后天线分别获取对应的回波数据，k＝1,2,...；根据k时刻前天线和后天线分别获取对应的回波数据，得出k时刻前天线回波数据的功率谱和k时刻后天线回波数据的功率谱；

步骤3，根据k时刻前天线回波数据的功率谱和k时刻后天线回波数据的功率谱，得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k；得出k时刻前天线与水平面的夹角θ_1,k、以及k时刻后天线与水平面的夹角θ_2,k；得出k时刻前天线得到的机车速度观测值z_1,k、以及k时刻后天线得到的机车速度观测值z_2,k：

$z_{1,k}=\frac{f_{d1,k}\mathrm{\λ}}{2\cos {\mathrm{\θ}}_{1,k}},z_{2,k}=\frac{f_{d2,k}\mathrm{\λ}}{2{\cos \mathrm{\θ}}_{2,k}}$

其中，λ表示测速雷达发射信号的载波波长；得出k时刻机车速度观测值z_k，如果k≤3，则k时刻机车速度的估计值v_k＝z_k，机车雷达测速过程完毕，否则，则转至步骤4；

步骤4，当k≥4时，根据k-1时刻机车的运动状态矢量x_k-1，得出k时刻机车的运动状态的预测矢量k-1时刻机车的运动状态矢量x_k-1是3行的列向量，k-1时刻机车的运动状态矢量的元素为k-1时刻机车的速度、k-1时刻机车的加速度和k-1时刻机车的加加速度；设置k时刻增益矩阵K_k，K_k是大小为3×1的矩阵；得出k时刻机车的运动状态矢量x_k，H＝[1 0 0]；得出k时刻机车速度的估计值v_k，v_k为k时刻机车的运动状态矢量x_k的第1个元素。

本发明的特点和进一步改进在于：

在步骤2中，对k时刻前天线获取的回波数据进行N点傅里叶变换和求模平方运算，得到k时刻前天线回波数据的功率谱；对k时刻后天线获取的回波数据进行N点傅里叶变换和求模平方运算，得到k时刻后天线回波数据的功率谱。

在步骤3中，得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k的过程为：分别对k时刻前天线回波数据的功率谱和后天线回波数据的功率谱进行对数加权处理，k时刻前天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_f(i)_k为：E_f(i)_k＝P_f(i)_k·log₂(i)，i＝1,2,...,N，N表示对回波数据进行傅里叶变换的点数；P_f(i)_k表示k时刻前天线回波数据中第i个频点对应的功率谱幅度；

k时刻后天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_b(i)_k为：E_b(i)_k＝P_b(i)_k·log₂(i)，P_b(i)_k表示k时刻后天线回波数据中第i个频点对应的功率谱幅度；

利用以下公式得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k：

$f_{d1,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_f{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_f{\left(i\right)}_k},f_{d2,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_b{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_b{\left(i\right)}_k}$

其中，f(i)表示回波数据的功率谱的第i个频点。

在步骤3中，根据以下方程组，得出k时刻前天线与水平面的夹角θ_1,k、以及k时刻后天线与水平面的夹角θ_2,k：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1,k}+{\mathrm{\θ}}_{2,k}=\mathrm{\θ}\\ \frac{f_{d1,k}}{f_{d2,k}}=\frac{{\cos \mathrm{\θ}}_{1,k}}{\cos {\mathrm{\θ}}_{2,k}}\end{array}\right.,$

其中，θ表示前天线的法线方向与后天线的法线方向的夹角。

所述步骤4的具体子步骤为：

(4.1)当k＝4时，得出k-1时刻机车的运动状态矢量x_k-1，x_k-1＝[v₃ a₃ b₃]^T，其中，上标T表示矩阵或向量的转置，v₃＝z₃；T为每两个相邻时刻的时间间隔；得出k-1时刻状态协方差矩阵P_k-1，P_k-1＝I，I是大小为3×3的单位矩阵；

(4.2)得出k时刻机车的运动状态的预测矢量其中，x_k-1＝[v_k-1 a_k-1 b_k-1]^T，v_k-1表示k-1时刻机车的速度，a_k-1表示k-1时刻机车的加速度，b_k-1表示k-1时刻机车的加加速度；F_k|k-1表示大小为3×3的状态转移矩阵；

计算得出k时刻预测误差协方差矩阵其中，上标T表示矩阵或向量的转置，P_k-1表示k-1时刻状态误差协方差矩阵，Q_k-1表示设定的状态噪声协方差矩阵；

(4.3)得出k时刻观测噪声的协方差R_k，当k≤L+1时，R_k＝1，L为设定的自然数，当k>L+1时，k时刻观测噪声的协方差R_k为：

$R_k=\frac{1}{L}\underset{m=k-L-1}{\overset{k-1}{\mathrm{\Σ}}}{(v_m-z_m)}^2$

其中，v_m表示m时刻机车的速度，z_m表示m时刻机车速度观测值；

计算得出k时刻增益矩阵K_k，

(4.4)得到k时刻机车的运动状态矢量得出k时刻机车速度的估计值v_k，v_k为k时刻机车的运动状态矢量x_k的第1个元素；

得出k时刻状态协方差矩阵P_k，

(4.5)令k的值自增1，返回至子步骤(4.2)。

本发明的有益效果为：1)本发明采用对数加权的方法进行低频干扰抑制，可以更准确地估计机车速度，减少误报现象。2)本发明采用雷达位置自校准方法，可以实时地调整前后两天线的测速结果，防止由于机车颠簸震动时产生的测速误差。3)本发明采用自适应卡尔曼滤波方法，滤除了观测噪声，克服了传统滑动平均方法出现的速度滞后现象，可以更为准确得到机车速度的估计值。

附图说明

图1为本发明的基于自适应卡尔曼滤波的机车雷达测速方法的流程框图；

图2为雷达测速装置与光电测速装置进行融合后的机车测速的结构示意图；

图3为本发明的测速雷达的自校准模块示意图；

图4为自适应卡尔曼滤波模块的处理过程示意图；

图5为对实测数据采用传统机车雷达测速方法得出的机车速度与时间的关系示意图；

图6为对实测数据采用本发明得出的机车速度与时间的关系示意图；

图7为对仿真数据分别采用本发明和滑动平均方法得出的机车速度与时间的关系示意图；

图8为图7的第一局部放大图；

图9为图7的第二局部放大图。

具体实施方式

下面结合附图对本发明作进一步说明：

本发明提出了一种基于自适应卡尔曼滤波的机车雷达测速方法。参照图1，为本发明的基于自适应卡尔曼滤波的机车雷达测速方法的流程框图。该机车雷达测速方法包括计算功率谱模块、抑制低频干扰模块、频谱重心法预估速度模块、雷达位置自校准与速度解算模块以及自适应卡尔曼滤波模块。计算功率谱模块，是将前天线和后天线AD采样后的两路数据首先进行FFT运算，然后进行求模平方运算，得到回波的功率谱，从而将回波数据由时域转到频域，利于后续对机车多普勒频率的提取。

参照图2，为雷达测速装置与光电测速装置进行融合后的机车测速的结构示意图。如果雷达安装在在机车底部，则雷达受到低频干扰的影响非常严重。在本发明中，抑制低频干扰模块，是对前天线和后天线的回波功率谱分别乘以一个对数函数，从而减弱低频部分的强度，提高高频部分的强度，达到抑制低频干扰的目的。

频谱重心法预估速度模块，用于预估计雷达的前后天线各自的回波功率谱重心位置对应的机车速度，并将前后天线对应的预估计多普勒频率送入速度解算模块。机车测速雷达的回波功率谱曲线是一个受到路基反射面变化、机车震动和噪声等因素影响的随机起伏且不对称的包络，直接采用最大值作为机车的多普勒估计值会有很大的误差。为了减小测速误差，本模块采用功率谱的重心位置进行预估。

雷达位置自校准与速度解算模块，是通过对雷达前后两天线相对位置予以固定，保证二者与水平面的夹角和一定，通过解方程实时校准天线与水平线的夹角，然后再由多普勒雷达公式解算机车运行速度，避免由于机车颠簸震动引起的测速误差。

在自适应卡尔曼滤波模块中，首先用前一时刻的机车状态预测机车当前时刻的运动状态，并利用滑窗方法实时计算观测噪声的方差，以适应不同路面情况下雷达的观测误差波动，然后再利用观测值(进过自校准和解算后得到的前后天线速度平均值)对机车速度的预测值进行更新，得到最终机车运行速度的估计值。

本发明的基于自适应卡尔曼滤波的机车雷达测速方法包括以下步骤：

步骤1，在机车静止时，在机车底部安装测速雷达；参照图3，为本发明测速雷达的自校准模块示意图。测速雷达包括位于测速雷达前部的前天线、以及位于测速雷达后部的后天线，分别用于接收回波数据。前天线的上端(即前天线从下端到上端的方向，也称前天线的辐射方向)朝向机车的前方，后天线的上端(后天线从下端到上端的方向，也称后天线的辐射方向)朝向机车的后方，也就是说，前天线与后天线的辐射方向相反。前天线与路面的夹角(为锐角)表示为θ₁，后天线与路面的夹角(为锐角)表示为θ₂，前天线的法线方向(即与前天线垂直的方向)与后天线的法线方向(与后天线垂直的方向)的夹角为θ(为已知量)，由几何结构可知，θ₁+θ₂＝θ。在初始时刻(安装上测速雷达的时刻)θ₁＝θ₂，当机车开始运行时，由于机车的颠簸等因素，测速雷达与路面的相对位置可能发生变化，此时θ₁不一定与θ₂相等，但是测速雷达的前天线与后天线的相对位置会保持不变，此时，依然满足θ₁+θ₂＝θ。

步骤2，计算功率谱模块的处理过程。

当机车开始运行后，测速雷达向外发射信号，在k时刻前天线通过AD采样获取对应的回波数据，在k时刻后天线通过AD采样获取对应的回波数据，k＝1,2,...；对k时刻前天线获取的回波数据进行FFT运算(N点傅里叶变换)和求模平方运算，得到k时刻前天线回波数据的功率谱；对k时刻后天线获取的回波数据进行FFT运算(N点傅里叶变换)和求模平方运算，得到k时刻后天线回波数据的功率谱。

步骤3，得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k；得出k时刻前天线与路面的夹角θ_1,k、以及k时刻后天线与路面的夹角θ_2,k；由以下公式计算得出k时刻前天线得到的机车速度观测值z_1,k、以及k时刻后天线得到的机车速度观测值z_2,k：

$z_{1,k}=\frac{f_{d1,k}\mathrm{\λ}}{2\cos {\mathrm{\θ}}_{1,k}},z_{2,k}=\frac{f_{d2,k}\mathrm{\λ}}{2{\cos \mathrm{\θ}}_{2,k}}$

其中，λ表示测速雷达发射信号的载波波长；得出k时刻机车速度观测值z_k，如果k≤3，则k时刻机车速度的估计值v_k＝z_k，机车雷达测速过程完毕，否则，则转至步骤4。

其具体子步骤为：

(3.1)抑制低频干扰模块的处理过程。对k时刻前天线回波数据的功率谱进行对数加权处理，对k时刻后天线回波数据的功率谱进行对数加权处理，以抑制低频部分的干扰。k时刻前天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_f(i)_k为：E_f(i)_k＝P_f(i)_k·log₂(i)，i＝1,2,...,N，N表示对回波数据进行傅里叶变换的点数；P_f(i)_k表示k时刻前天线回波数据中第i个频点对应的功率谱幅度。

k时刻后天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_b(i)_k为：E_b(i)_k＝P_b(i)_k·log₂(i)，P_b(i)_k表示k时刻后天线回波数据中第i个频点对应的功率谱幅度。

(3.2)频谱重心法预估速度模块的处理过程。

利用以下公式得出k时刻前天线回波数据的功率谱进行对数加权处理之后的回波功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱进行对数加权处理之后的回波功率谱包络区域的频点重心f_d2,k：

$f_{d1,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_f{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_f{\left(i\right)}_k},f_{d2,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_b{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_b{\left(i\right)}_k}$

其中，f(i)表示回波数据的功率谱的第i个频点(对于前天线和后天线来说，f(i)是一致的。每个时刻回波数据的功率谱的第i个频点均相同)。

(3.3)雷达位置自校准与速度解算模块的处理过程。

根据以下方程组，得出k时刻前天线与路面的夹角θ_1,k、以及k时刻后天线与路面的夹角θ_2,k：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1,k}+{\mathrm{\θ}}_{2,k}=\mathrm{\θ}\\ \frac{f_{d1,k}}{f_{d2,k}}=\frac{{\cos \mathrm{\θ}}_{1,k}}{\cos {\mathrm{\θ}}_{2,k}}\end{array}\right.,$

其中，θ表示前天线的法线方向(即与前天线垂直的方向)与后天线的法线方向(与后天线垂直的方向)的夹角，其为已知量。f_d1,k与f_d2,k通过子步骤(3.2)得到。

然后由以下公式计算得出k时刻前天线得到的机车速度观测值z_1,k、以及k时刻后天线得到的机车速度观测值z_2,k：

$z_{1,k}=\frac{f_{d1,k}\mathrm{\λ}}{2\cos {\mathrm{\θ}}_{1,k}},z_{2,k}=\frac{f_{d2,k}\mathrm{\λ}}{2{\cos \mathrm{\θ}}_{2,k}}$

其中，λ表示测速雷达发射信号的载波波长。

最后得出k时刻机车速度观测值z_k

$z_k=\frac{z_{1,k}+z_{2,k}}{2}$

如果k≤3，则k时刻机车速度的估计值v_k＝z_k，本发明的机车雷达测速过程完毕，否则，如果k>3，则转至步骤4。

步骤4，自适应卡尔曼滤波模块的处理过程。

参照图4，为自适应卡尔曼滤波模块的处理过程示意图。步骤4的具体子步骤为：

(4.1)当k＝4时，初始化机车的运动状态矢量，得出k-1时刻机车的运动状态矢量x_k-1，x_k-1＝[v₃ a₃ b₃]^T，其中，上标T表示矩阵或向量的转置，v₃表示k-1时刻机车的速度，v₃＝z₃；a₃表示k-1时刻机车的加速度，b₃表示k-1时刻机车的加加速度，T为观测间隔，即每两个相邻时刻的时间间隔。

初始化状态误差协方差矩阵，得出k-1时刻状态协方差矩阵P_k-1，P_k-1＝I，其中，I是大小为3×3的单位矩阵。

(4.2)由k-1时刻机车的运动状态矢量x_k-1，得出k时刻机车的运动状态的预测矢量其中，x_k-1＝[v_k-1 a_k-1 b_k-1]^T，v_k-1表示k-1时刻机车的速度，a_k-1表示k-1时刻机车的加速度，b_k-1表示k-1时刻机车的加加速度。F_k|k-1表示大小为3×3的状态转移矩阵，F_k|k-1为：

$F_{k|k-1}=\left[\begin{array}{ccc}1& T& T^2/2\\ 0& 1& T\\ 0& 0& 1\end{array}\right];$

计算得出k时刻预测误差协方差矩阵其中，上标T表示矩阵或向量的转置，P_k-1表示k-1时刻状态误差协方差矩阵，Q_k-1表示设定的状态噪声协方差矩阵，Q_k-1＝q·I，q为状态噪声强度，I是大小为3×3的单位矩阵。

(4.3)得出k时刻观测噪声的协方差R_k，当k≤L+1时，R_k＝1，L为设定的自然数(滑窗的长度)，当k>L+1时，k时刻观测噪声的协方差R_k为：

$R_k=\frac{1}{L}\underset{m=k-L-1}{\overset{k-1}{\mathrm{\Σ}}}{(v_m-z_m)}^2$

其中，v_m表示m时刻机车的速度，z_m表示m时刻机车速度观测值。

计算得出k时刻增益矩阵K_k，其中，H表示观测矩阵，H＝[1 0 0]；上标T表示矩阵或向量的转置。

(4.4)利用观测值更新当前时刻的预测值，得到k时刻机车的运动状态矢量x_k， $x_k={\stackrel{\‾}{x}}_k+K_k[z_k-H{\stackrel{\‾}{x}}_k].$

得出k时刻机车速度的估计值v_k，v_k为k时刻机车的运动状态矢量x_k的第1个元素。

得出k时刻状态误差协方差矩阵P_k，

(4.5)令k的值自增1，返回至子步骤(4.2)。

本发明的效果可以通过以下实测数据和仿真数据进一步说明：

参照图5，为对实测数据采用传统机车雷达测速方法(未进行自适应卡尔曼滤波处理)得出的机车速度与时间的关系示意图。参照图6，为对实测数据采用本发明得出的机车速度与时间的关系示意图。图5和图6中，横轴表示时间采样点，纵轴表示机车速度，单位为km/h。比较图5和图6可以发现，经过自适应卡尔曼滤波之后，机车速度的波动明显减小，估计误差的方差也减小。因此本发明对机车速度的估计稳定，而且精度较高。

参照图7，为对仿真数据分别采用本发明和滑动平均方法得出的机车速度与时间的关系示意图；参照图8，为图7的第一局部放大图；参照图9，为图7的第二局部放大图。图7至图9中，横轴表示时间采样点，纵轴表示机车速度，单位为km/h。图7至图9中，真实数据代表仿真数据的理论上的车速与时间的关系曲线，滑动平均结果代表利用滑动平均方法得出的机车速度与时间的关系曲线，卡尔曼滤波结果代表利用本发明得出的机车速度与时间的关系曲线。通过图7、图8和图9可以发现，采用滑动平均方法，估计误差较大，而且存在速度滞后的现象，即当机车加速运行时，测得的速度值比真实值偏小，而当机车减速运行时，测得的速度值比真实值偏大。本发明采用自适应卡尔曼滤波方法，得到的估计结果误差较小，较为稳定，而且不存在速度滞后的现象。

综上，本发明充分考虑到机车测速雷达的实际应用问题，根据雷达近场回波功率谱曲线的特征，首先通过对数加权抑制低频干扰，然后利用频谱重心法对机车多普勒频率进行预估，再通过固定前后两天线的相对位置实现机车颠簸震动时雷达位置的自校准，并平均前后天线的估计值，最后进行自适应卡尔曼滤波，最终得到机车运行速度的估计值。本发明估计误差小，稳定性好，实际运行情况良好。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.基于自适应卡尔曼滤波的机车雷达测速方法，其特征在于，包括以下步骤：

步骤1，在机车静止时，在机车底部安装测速雷达；测速雷达具有前后两个天线，前天线朝向机车的前方，后天线朝向机车的后方；

步骤2，当机车开始运行后，测速雷达向外发射信号，在k时刻前天线和后天线分别获取对应的回波数据，k＝1,2,...；根据k时刻前天线和后天线分别获取对应的回波数据，得出k时刻前天线回波数据的功率谱和k时刻后天线回波数据的功率谱；

步骤3，根据k时刻前天线回波数据的功率谱和k时刻后天线回波数据的功率谱，得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k；得出k时刻前天线与水平面的夹角θ_1,k、以及k时刻后天线与水平面的夹角θ_2,k；得出k时刻前天线得到的机车速度观测值z_1,k、以及k时刻后天线得到的机车速度观测值z_2,k：

$z_{1,k}=\frac{f_{d1,k}\mathrm{\λ}}{{2\cos \mathrm{\θ}}_{1,k}},z_{2,k}=\frac{f_{d2,k}\mathrm{\λ}}{{2\cos \mathrm{\θ}}_{2,k}}$

其中，λ表示测速雷达发射信号的载波波长；得出k时刻机车速度观测值z_k，如果k≤3，则k时刻机车速度的估计值v_k＝z_k，机车雷达测速过程完毕，否则，则转至步骤4；

步骤4，当k≥4时，根据k-1时刻机车的运动状态矢量x_k-1，得出k时刻机车的运动状态的预测矢量k-1时刻机车的运动状态矢量x_k-1是3行的列向量，k-1时刻机车的运动状态矢量的元素为k-1时刻机车的速度、k-1时刻机车的加速度和k-1时刻机车的加加速度；设置k时刻增益矩阵K_k，K_k是大小为3×1的矩阵；得出k时刻机车的运动状态矢量x_k，H＝[1 0 0]；得出k时刻机车速度的估计值v_k，v_k为k时刻机车的运动状态矢量x_k的第1个元素。

2.如权利要求1所述的基于自适应卡尔曼滤波的机车雷达测速方法，其特征在于，在步骤2中，对k时刻前天线获取的回波数据进行N点傅里叶变换和求模平方运算，得到k时刻前天线回波数据的功率谱；对k时刻后天线获取的回波数据进行N点傅里叶变换和求模平方运算，得到k时刻后天线回波数据的功率谱。

3.如权利要求2所述的基于自适应卡尔曼滤波的机车雷达测速方法，其特征在于，在步骤3中，得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k的过程为：分别对k时刻前天线回波数据的功率谱和后天线回波数据的功率谱进行对数加权处理，k时刻前天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_f(i)_k为：E_f(i)_k＝P_f(i)_k·log₂(i)，i＝1,2,...,N，N表示对回波数据进行傅里叶变换的点数；P_f(i)_k表示k时刻前天线回波数据中第i个频点对应的功率谱幅度；

k时刻后天线回波数据中第i个频点对应的功率谱幅度的对数加权处理后结果E_b(i)_k为：E_b(i)_k＝P_b(i)_k·log₂(i)，P_b(i)_k表示k时刻后天线回波数据中第i个频点对应的功率谱幅度；

利用以下公式得出k时刻前天线回波数据的功率谱包络区域的频点重心f_d1,k、以及k时刻后天线回波数据的功率谱包络区域的频点重心f_d2,k：

$f_{d1,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_f{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_f{\left(i\right)}_k},f_{d2,k}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}f\left(i\right)E_b{\left(i\right)}_k}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_b{\left(i\right)}_k}$

其中，f(i)表示回波数据的功率谱的第i个频点。

4.如权利要求1所述的基于自适应卡尔曼滤波的机车雷达测速方法，其特征在于，在步骤3中，根据以下方程组，得出k时刻前天线与路面的夹角θ_1,k、以及k时刻后天线与路面的夹角θ_2,k：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{1,k}+{\mathrm{\θ}}_{2,k}=\mathrm{\θ}\\ \frac{f_{d1,k}}{f_{d2,k}}=\frac{{\cos \mathrm{\θ}}_{1,k}}{{\cos \mathrm{\θ}}_{2,k}}\end{array}\right.,$

其中，θ表示前天线的法线方向与后天线的法线方向的夹角。

5.如权利要求1所述的基于自适应卡尔曼滤波的机车雷达测速方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)当k＝4时，得出k-1时刻机车的运动状态矢量x_k-1，x_k-1＝[v₃ a₃ b₃]^T，其中，上标T表示矩阵或向量的转置，v₃＝z₃,T为每两个相邻时刻的时间间隔；设定k-1时刻状态协方差矩阵P_k-1，P_k-1＝I，I是大小为3×3的单位矩阵；

(4.2)得出k时刻机车的运动状态的预测矢量 其中，x_k-1＝[v_k-1 a_k-1 b_k-1]^T，v_k-1表示k-1时刻机车的速度，a_k-1表示k-1时刻机车的加速度，b_k-1表示k-1时刻机车的加加速度；F_k|k-1表示大小为3×3的状态转移矩阵；

计算得出k时刻预测误差协方差矩阵 其中，上标T表示矩阵或向量的转置，P_k-1表示k-1时刻状态误差协方差矩阵，Q_k-1表示设定的状态噪声协方差矩阵；

(4.3)得出k时刻观测噪声的协方差R_k，当k≤L+1时，R_k＝1，L为设定的自然数，当k>L+1时，k时刻观测噪声的协方差R_k为：

$R_k=\frac{1}{L}\underset{m=k-L-1}{\overset{k-1}{\mathrm{\Σ}}}{(v_m-z_m)}^2$

其中，v_m表示m时刻机车的速度，z_m表示m时刻机车速度观测值；

计算得出k时刻增益矩阵K_k， $K_k={\stackrel{\‾}{P}}_k{H}^T{(H{\stackrel{\‾}{P}}_k{H}^T+R_k)}^{-1};$

(4.4)得到k时刻机车的运动状态矢量x_k，得出k时刻机车速度的估计值v_k，v_k为k时刻机车的运动状态矢量x_k的第1个元素；

得出k时刻状态协方差矩阵P_k，

(4.5)令k的值自增1，返回至子步骤(4.2)。