陀螺仪噪声对轨道高低不平顺检测的影响

程朝阳 赵延峰 陈仕明 韩志 张茂轩 魏世斌

1.中国铁道科学研究院集团有限公司基础设施检测研究所，北京 100081；2.中国铁道科学研究院研究生部，北京 100081

轨道检测是及时掌握轨道几何形位、指导工务养护、保证铁路安全运行的重要保障［1］。基于惯性基准法的轨道检测一般采用光纤陀螺提供惯性基准，因此有必要分析光纤陀螺噪声对轨道检测的影响［2］。

衡量光纤陀螺性能的基本指标是零偏稳定性。零偏稳定性是指输入角速率为零，采样周期一定，陀螺输出角速率的标准差。零偏稳定性是一个综合指标［3］，包含多项误差。Allan 方差法［4］是在时域上对频域特性进行分析的方法，可估计光纤陀螺角度随机游走、零偏不稳定性、角速率随机游走的功率谱密度［5］。三种噪声的功率谱密度分别与频率f的γ次方成反比，满足幂律谱模型，属于1∕f γ分形噪声，其功率谱密度的傅里叶逆变换可仿真三种噪声的时间样本序列。

用正弦波窄带信号模拟轨道不平顺空间连续波形，基于链式法则对时间和空间求导获取该空间波形对应的角速率时间序列，进行抗混滤波、空间采样、速度补偿、高通滤波等高低不平顺检测算法流程，对比叠加噪声的角速率通过相同算法流程的结果，可表征陀螺仪噪声对检测结果的影响。

1 用Allan方差估计陀螺仪噪声功率谱密度

1.1 用Allan方差和噪声功率谱密度的关系

令{x(t)，- ∞＜t＜∞}为各态历经的零均值平稳过程（t为时间），样本函数的时间平均可替代过程的总体平均，获取过程{x(t)}的一个样本函数xj，j=0、1、…、M，将其分成若干个连贯数组［7］，每组有连贯的n个数据点。相邻的第k、k+1 组的平均值分别表示为

式中：tk、tk-1分别为相邻两组时间起始值；τ为时间差，tk=tk-1+τ。

由式（1）、式（2）可知

因为x(t)为各态历经的平稳过程，自相关函数R(τ)与功率谱密度Sω为傅里叶变换对，且只与时间差τ有关，总平均功率与平均交流功率相同，自相关函数即为协方差，故协方差与功率谱密度为傅里叶变换对，表达式为

将式（5）代入式（4），可得

又因协方差为实函数，功率谱密度为偶函数，故Allan方差为

式（8）建立了Allan 方差和功率谱密度的关系，可有效分析陀螺各种随机误差［8］。

1.2 陀螺仪噪声功率谱对应的Allan方差

基于惯性基准法的高低不平顺检测要求陀螺测量输出的角速率在短时间内保持较高精度，故角度随机游走是影响检测的主要因素；同时轨道几何检测中陀螺仪工作时间通常较长，因此导致长期漂移误差的零偏不稳定性和角速率随机游走也成为影响检测精度的主要因素［9］。本节主要对光纤陀螺的角度随机游走、零偏不稳定性、角速率随机游走噪声的功率谱密度与Allan方差的关系进行说明。

1.2.1 角度随机游走

角度随机游走表现为陀螺输出的角速率白噪声。角度随机游走反映了角速率信号中宽带白噪声的特性，是宽带角速度白噪声积分的结果，即陀螺从0时刻起累积的总角增量误差表现为随机游走［10］，而每一时刻等效角速度误差表现为白噪声。这种噪声是服从高斯分布、具有零均值的广义平稳白噪声，其功率谱密度为

式中：Sω，ARW(f)为角度随机游走噪声功率谱密度；N为角度随机游走噪声功率谱密度系数。

因此，N由Allan方差表达为

1.2.2 零偏不稳定性

采样周期较短时，光纤陀螺的零偏稳定性随采样周期的增大而减小。采样周期增加到某一值时，零偏稳定性会限制在某一量级内，此时光纤陀螺输出中的低频漂移开始起作用，通过增加平均时间无法加以抑制。零偏不稳定性可用于描述陀螺角速率测量数据中的低频漂移。当光纤陀螺工作时间较长时，零偏不稳定性是角误差的主要来源。

零偏不稳定性具有低频特性，该噪声表现为偏置不稳定性，最简单的表达形式为随机常数。考虑到陀螺的动态时变特性，选择用1∕f来模拟偏置不稳定性噪声功率谱，其表达式为

式中：Sω，BI(f)为零偏不稳定性噪声功率谱密度；B为零偏不稳定性噪声功率谱密度系数。

因此，B由Allan方差表达为

1.2.3 角速率随机游走

陀螺角加速率误差表现为宽带白噪声，角速率随机游走（Rate Random Walk，RRW）是角加速率经过积分后叠加在陀螺角速率观测量上的一部分误差，这类误差的方差随采样周期增大而增大。积分后的角速率误差表现为随机游走，可通过对白噪声的积分来近似得到，其随机统计特性服从布朗运动特征。角速率功率谱密度表达式为

式中：Sω，RRW(f)为角速率随机游走功率谱密度；K为角速率随机游走功率谱密度系数。

K对应的角速率误差S∫ωt(f)为

角速率随机游走对应的Allan方差σ2RRW(τ)为

因此，K由Allan方差表达为

1.3 用Allan方差计算陀螺仪噪声功率谱密度系数

按照上述分析，陀螺仪输出角速率构成的一个样本空间采用Allan 方差法进行处理，可以辨识出角度随机游走、零偏不稳定性、角速率随机游走各项误差的噪声系数。计算Allan方差的具体步骤如下。

①采集陀螺仪角速率静态输出数据，以采样频率f采集光纤陀螺的输出角速率，获取长度为L的样本空间。

②将长度为L的样本空间每m个数据分成一组（m＜L∕2），得到L∕m个独立的数组。

上述步骤的Allan方差σ2(τ)计算公式为

按照上述方法采集陀螺仪静态输出角速率，计算Allan方差曲线，结果见图1。

图1 陀螺仪Allan方差曲线

用最小二乘法拟合Allan 方差曲线，获取噪声功率谱密度系数，结果见图2。可知，角度随机游走噪声功率谱密度系数N= 12.608 2°∕h1∕2，零偏不稳定性噪声功率谱密度系数B= 1.030 4°∕h，角速率随机游走功率谱密度系数K= 0.060 0°∕h3∕2。

图2 陀螺仪噪声功率谱密度系数

2 基于FFT生成陀螺仪噪声

2.1 陀螺仪噪声的仿真

仿真生成角度随机游走、零偏不稳定性、角速率随机游走等1∕f γ分形噪声一般有两种方法：①利用小波［11］方法通过Karhunen-Loeve 展开式来合成；②利用FFT（Fast Fourier Transform）方法将白噪声［12］通过一个线性滤波器，利用测定的N、B、K即可建立所有功率谱密度频率的比例因子，采用IFFT（Inverse Fast Fourier Transform）得到分形噪声的时间序列样本。

本文选用第二种方法，如图3 所示，随机信号x(t)经线性滤波器后得到分形信号y(t)，其功率谱密度为|H(f)|2Sx(f)。

图3 幂律噪声仿真算法

如图4所示，基于同一个白噪声信号，仿真了一组功率谱密度分别为N∕f0、B∕f1、K∕f2（自上而下）的分形噪声，分别对应陀螺仪输出角速率信号中的角度随机游走d、零偏不稳定性b和角速率随机游走ωg，幅值归一化，单位为°∕h。

图4 分形噪声仿真结果

2.2 仿真陀螺仪噪声的验证

理论上分形噪声的功率谱系数-2～0是噪声功率谱密度对应的傅里叶频率f的幂，故功率谱密度线性拟合的斜率可以验证仿真噪声的真实性。

如图5 所示，自上而下依次是角度随机游走d、零偏不稳定性b和角速率随机游走ωg的功率谱密度，仿真生成的噪声功率谱密度为蓝色曲线，其一次线性拟合的结果为红色直线，直线的斜率即为噪声功率谱密度的系数。

图5 仿真幂律噪声拟合结果

陀螺仪噪声功率谱密度系数见表1。可知，FFT方法生成的分形噪声的实际仿真值与理论值的最大相对误差不超过1%，且随着取样点数的增加，仿真值将更接近理论值。本文所提方法的实际仿真值与理论值的最大相对误差不超过2%。

表1 陀螺仪噪声功率谱密度系数

将三种噪声叠加生成陀螺噪声信号（图6），即为影响陀螺仪高低不平顺检测精度的主要噪声信号。

图6 陀螺仪噪声仿真信号

3 高低不平顺的陀螺仪信号

3.1 高低不平顺陀螺仪实测信号的仿真

轨道不平顺的空间曲线幅值y（单位：mm）和轨道不平顺对应的陀螺仪角速率ω（前后两空间采样点之间的平均角速率，单位：rad∕s）的表达式分别为

式中：x为空间里程，m；λ为轨道不平顺波长，m；A为轨道不平顺幅值，mm；v为速度，m∕s。

3.2 陀螺仪实测信号融合噪声

陀螺仪信号融合噪声算法见图7。其中，ωn为二阶环节的自然振荡频率；s为复频率；ωout为陀螺的输出频率；在忽略标度因数误差、安装误差及频带宽度的情况下，陀螺仪信号融合噪声后输出为ωm，其表达式为

图7 陀螺仪信号融合噪声算法

4 高低不平顺算法仿真

4.1 高低不平顺算法（图8）

图8 高低不平顺算法

仿真后的陀螺仪信号融合噪声的信号ωm在通过抗混滤波、空间采样、速度补偿、高通滤波等算法流程后即可得出高低不平顺。

4.2 抗混滤波及速度补偿

抗混滤波H1(s)计算公式为

式中：Ω1为滤波器固有频率，rad∕s。

时间空间的频率对应关系为

式中：Ω为时间角频率，rad∕s；l为空间波长，m；φ为空间频率，1/m。

由s= jΩ及式（24），式（23）转化为

式中：H1(φ)关于空间频率和车速v的函数。

对于同样的时间截止角频率，车速不同，则对应的空间截止频率不同。为消除速度对抗混滤波器在空间截止波长的影响，采用去移变滤波器H1(z)消除H1(s)的移变特性，计算式为

式中：a1、a2为常系数；z为拉普拉斯变换因子，z= esLi∕vi，Li为空间采样距离，vi为采样距离对应的车速。

则H1(z)的拉普拉斯变换表达式为

由H1(s)H1z(s) ≡1，取其物理意义为空间波长为零时满足条件，此时Z平面的零点与S平面的极点相互抵消。H1(z)可写成

式（28）中，H1(z)的系数中含有Ti，与列车运行速度有关，从而实现滤波特性随车速的变化而变化。

4.3 高通滤波器

高通滤波器P(z)采用矩形窗、三角窗并联的形式，其表达式为

式中：b1、b2、b3、p、q、r、u均为常系数，随截止频率变化而变化。

按照式（29）可将滤波器系统函数利用转化为差分方程形式，或单个滤波器时域补零后卷积，二者均可快速实现该高通滤波器，25 ～ 120 m 不同截止波长的对应的幅频特性曲线如图9所示。

图9 高通滤波器幅频特性

4.4 仿真结果

按照图8 算法流程，对于模拟的幅值为12.5 mm的正弦轨道高低不平顺波形，陀螺仪叠加仿真噪声与非叠加仿真噪声的差值见图10。可以看出角度随机游走、零偏不稳定性、角度随机游走、三种噪声同时叠加对轨道高低不平顺检测精度的影响。可知，三种噪声中，角度随机游走对陀螺仪高低不平顺检测结果影响较小，而角速率随机游走和零偏不稳定性噪声在算法中并未受抑制，是检测中惯性基准误差的主要来源。

图10 高低不平顺仿真结果对比

5 结论

1）在陀螺仪轨道高低不平顺检测中，零偏不稳定性、角速率随机游走是影响系统检测精度的主要因素。

2）利用Allan 方差拟合出陀螺仪噪声的功率谱密度系数，并基于该系数使用FFT 逆变换仿真符合该功率谱密度的噪声时域信号，能够有效仿真生成符合陀螺仪特性的噪声。

3）陀螺仪噪声对检测结果的重复性具有一定影响，选取合适精度的陀螺仪，尤其是零偏不稳定性、角速率随机游走噪声较小的陀螺仪可有效提高检测结果的重复性。

4）本文分析结果为轨道几何高低不平顺检测系统的惯性组件选型提供了参考，也为提升高低不平顺检测的重复性提供了方向。惯性器件是轨道几何检测中的重要传感器，也应用于超高、水平检测项，可以基于相同思路对超高、水平检测算法流程的惯性噪声误差进行分析，研究影响测量精度的主要因素。