基于卡尔曼滤波的高精度温度脉动仪设计

黎一宏，常帅，张恒，步宏坤

（长春理工大学 光电工程学院，长春 130022）

光学湍流效应是制约空间激光通信发展的重要因素，由于大气湍流造成折射率的随机起伏会影响空间激光通信的质量，因此研究并测量大气折射率结构常数是研究光波传输的重要环节［1-2］。国内外都对大气折射率结构常数的测量进行了大量的探索。 研究表明，大气湍流对光波在大气传输的影响主要是由于温度的起伏引起的。温度脉动仪通过对空间两点的微温脉动量实时测量，反演大气折射率结构常数，是一种较为常用的大气湍流强度测量仪器［3-4］。温度脉动仪探测的微温脉动量变化量小，频率高，现有的温度脉动仪设备对高频温度脉动信号采集时，噪音较大。 现采用卡尔曼滤波原理对采集到的温度脉动信号进行滤波处理，减小信号的噪音毛刺，提高设备的测量精度，并使用温度脉动仪对长春地区大气折射率结构常数进行测量。

1 实验原理与设备研制

1.1 温度脉动仪测量原理及设备研制

科尔莫戈罗夫（Kolmogorov）的“2/3”定律认为：在湍流的惯性区内，两点间的折射率结构常数只与两点间距离的2/3 次方有关，与两点的位置和相对方向无关。 在此定律基础上，可推导出折射率结构常数的计算公式如下：

由公式（1）和公式（2）可知，当气压P、温度T、距离r和温度差的系综平均值为已知量时，可计算得出折射率结构常数为。气压P和温度T可通过气压传感器和温度传感器获得，使用米尺和游标卡尺等工具可以精确地测得距离r，温度差的系综平均值需要对温度脉动仪测量的数据进行反演计算获得。由于温度差的系综平均值不能被直接测量，则需要一个能够反映其特征的可测量值来代替。温度脉动仪采用直径25 µm、电阻值为50 Ω 的铂丝作为传感器（两个铂丝探头的电阻差值不大于标准值的0.4%），感应大气湍流惯性区内空间两点温度的细微差异，将表示大气状态的物理量温度差转化为可测量的电信号电压差。铂丝的电阻值随温度变化的公式为：

由此得出：

式中，R0为T0下的电阻值，单位为Ω；α为电阻温度系数，单位为K-1；Tω为铂丝实际温度；T0为铂丝标定时的温度。

铂丝的温度系数为（3.908 3×10-3）（ppm/℃），因此根据公式（4），50 Ω 的铂丝电阻在发生1 ℃的温度变化时，其阻值变化为：

1 Ω 的电阻变化经过温度脉动仪测温电路的放大，最终引起的电压变化为5.2 V，因此可以推出当铂丝电阻探头发生1 ℃的温度变化时，最终引起电路上的电压变化为：

因此温度脉动仪的温度引起电压改变的系数为：

由此关系可以利用铂丝电阻的电压变化反演温度的脉动变化。

本文采用的温度脉动仪设计方案使用四线制测量电阻法，对铂丝电阻的阻值进行测量，最大限度消除导线电阻引入的误差，使用2.3 mA的恒流源对用于探测温度脉动的铂丝电阻通电，电阻流过的电流极低，减小了电流增温效应对测量的影响。采用24 位低噪音ADC 采集芯片采集电压值，ADC 采集芯片的采集频率最高可达到100 Hz，满足对高频温度脉动量的测量需求，经过单片机将数据输出，温度脉动仪的电子学框图如图1 所示。

图1 温度脉动仪电子学框图

铂丝探头受到温度脉动电阻值发生变化，两探头上的电阻结果测温电路差分放大后由ADC采集模块准确采集探头上的电压差值，传入信号处理模块。信号处理模块主要由两块单片机组成，使用一块STM32F103 单片机接收ADC 模块采集出的数据，对采集的数据做滤波算法处理后发送到一块STM32F407 单片机上，STM32F407单片机收集采集到的滤波后的电压差值数据与BME680 传感器传输来的气象数据，一并通过串口输出到上位机，同时控制OLED 显示屏实时显示数据，上位机通过软件对数据进行反演检测大气折射率结构常数。

1.2 基于卡尔曼滤波的温度脉动仪改进设计

卡尔曼滤波的原理包括两个主要过程：预估与校正。预估过程主要是利用时间更新方程建立对当前状态的先验估计，及时向前推算当前状态变量和误差协方差估计的值，以便为下一个时间状态构造先验估计值。校正过程负责反馈，利用测量更新方程在预估过程的先验估计值及当前测量变量的基础上建立起对当前状态改进的后验估计。 这个过程称之为预估-校正过程，对应的估计算法称为预估-校正算法［8-9］。以下给出卡尔曼滤波的时间更新方程和状态更新方程。

时间更新方程：

状态更新方程：

式中，Xk为系统状态变量；Uk是系统外部输入；A是系统状态转移矩阵；B是系统输入控制矩阵；是先验估计的预测协方差矩阵；k为算法的迭代次数；Q是控制过程噪音的协方差矩阵；Gk是卡尔曼增益；H是观测矩阵；R为观测噪音的协方差矩阵；Zk为观测值；I为单位矩阵。

大气湍流的扰动频率很高，通常可以达到上百赫兹，而温度脉动引起大气折射率结构常数变化是一个缓慢的过程，因此在采集温度脉动信号时可采用卡尔曼滤波的思路，利用上一个采样周期的系统状态定义为预估的当前状态，结合测温电路采集到的当前周期的数据进行加权组合得到系统的真实状态。在对温度脉动进行高频采集的同时降低由于硬件或环境异常扰动造成的噪音。由于是对采集到的单变量数据进行滤波，所以矩阵的维数都为1。在温度脉动仪信号处理模块中的STM32F103 单片机上建立卡尔曼滤波算法模型，对采集到的电压数据进行滤波。 首先定义卡尔曼滤波结构体包含：上次估算协方差Pk- 1，当前估算协方差Pk，卡尔曼滤波器输出Out，卡尔曼增益G，过程噪音协方差Q，观测噪音协方差R，并将估算协方差、过程噪音协方差、观测噪音协方差参数初始化。ADC 芯片采集的电压数据为Input，进入卡尔曼滤波循环，根据公式（5）～（9）的卡尔曼滤波原理，设计温度脉动仪的卡尔曼滤波模型：

预测协方差方程：

卡尔曼增益方程：

更新最优值方程：

更新预测协方差方程：

数据进入卡尔曼滤波循环后，将上一时刻输出值数据作为当前时刻的预测值，使用卡尔曼增益进行加权计算后得到这一时刻的状态变量最优值后，输出经过计算的数据得到滤波后的电压值，将电路噪音进行滤波，同时可以减小外界异常干扰造成电压值不稳定而形成的毛刺。

图2 为STM32F103 单片机的工作流程图。在系统开机后，首先其内部进行初始化工作，主要针对ADC 模块的工作模式进行配置，使其能够保障转换精度的同时，较快速地采集铂丝探头两端的电压信号。之后该单片机进入循环工作状态，读取ADC 模块输出的电压数据，在单片机内写入利用卡尔曼滤波算法对采集得到的信号进行迭代滤波，将上一次的输出值作为下一次循环的预测值进行数据迭代，对采集到的电压信号进行加权计算后得到当前时刻的最优值，降低电压信号内的噪声信息，提高采集精度，之后单片机将最优值数据输出，通过串口发送至STM32F407 内核单片机中，进行后续处理。

图2 温度脉动仪STM32F103 单片机工作流程图

2 实验结果与讨论

2.1 卡尔曼滤波器滤波效果

为验证卡尔曼滤波降低温度脉动仪噪音的效果，现将温度脉动仪中测温电路采集到的电压原始数据单独导出到上位机上的卡尔曼滤波程序中，需要手动调节卡尔曼滤波参数：过程噪声协方差Q，观测噪声协方差R，得到好的滤波效果后，将卡尔曼滤波器的初始值烧录到单片机中，使温度脉动仪在采集电压后立即滤波，简化后续的数据处理。原始电压数据经滤波后的对比如图3 所示。

图3 卡尔曼滤波器的电压滤波效果

由此可见经过滤波后的电压数据在变化趋势上符合原始数据的同时，极大减少了数据中的毛刺噪音，极大提高了数据的可靠程度，在后续通过电压反演温度脉动变化时，极大提升数据的精度。

2.2 温度脉动仪噪音测试

为进一步研究卡尔曼滤波对温度脉动仪精度的提升，对温度脉动仪进行噪音测试。 现将两个精密的50 Ω 定值电阻替换掉原来的铂丝传感器探头安装到温度脉动仪的探测天线两端，记录下原始的电压数据，然后将通过卡尔曼滤波器后的数据与原始数据对比，结果如图4 所示。

图4 温度脉动仪噪音对比

使用定值电阻替代传感器探头，两端的电压差值应该恒定无变化，电压差值应该呈现一条平稳的直线，此情况下可以对温度脉动仪的噪音进行测量。

实验前经过精密电阻计测得两精密定值电阻的电阻值之差为0.011 Ω，此阻值差值经过温度脉动仪测温电路差分放大后带来的电压差值为0.059 V，因此噪音测试时电压差值始终有0.059 V 左右的输出。结果对比可知，在添加滤波算法前温度脉动仪的原始数据噪音在0.002 V范围内，噪音总体的毛刺还是比较突出；经过卡尔曼滤波后，噪音控制在0.001 V 以内，噪音降低50%，消除数据中的毛刺噪音，滤波效果极好。按照本文中温度脉动仪测量电压和探测温度的转换关系，改进后温度脉动仪的噪音对应的温度脉动在0.001 K 范围内。根据公式（1）和公式（2）的大气折射率结构常数计算公式，0.001 K 的温度脉动噪音计算为大气折射率结构常数为10-18m-2/3量级，处于温度脉动仪的测量下限，远小于正常测量时大气湍流数据，对测量结果影响极低，提高了设备的测量准确性。

2.3 温度脉动仪大气折射率结构常数测量实验

使用改进后的温度脉动仪在长春地区搭建大气折射率结构常数测量平台，地理位置为长春市南关区生态大街6666 号创业服务中心西附楼。测量点在四楼天台，海拔212 m，经度125.39，纬度43.78。温度脉动仪搭载在距下垫面1 m 的水泥台座上，下垫面环境为人工塑料草皮。 实验现场如图5 所示，温度脉动仪放置在天台边缘，四周开阔无遮挡，无高热辐射源或强电磁场辐射源，由于长春冬季天气寒冷，实验场地有积雪覆盖。

图5 温度脉动仪实验现场

实验结果如图6 所示，温度脉动仪在观测点成功绘制大气折射率结构常数在一整天的变化曲线。12 月23 日长春地区的日出时间为7 点10分，日中时间为11 点37 分，日落时间为16 点05分。 在日出之前温度脉动量较低，大气湍流较弱，大气折射率结构常数在10-14m-2/3附近波动，7点钟日出以后，温度脉动量开始增加，大气折射率结构常数开始升高，在日中前的11 点附近达到一个峰值，最高达到5 × 10-13m-2/3，然后在日中时温度脉动量较为平稳，大气折射率结构常数有所回落，到10-13m-2/3附近波动，在日中过后13点30 分时又迎来一个回升，达到5 × 10-13m-2/3。 之后便逐渐降低，在日落以后大气折射率结构常数回到当日低点，在10-14m-2/3附近波动。

图6 大气折射率结构常数测量结果

3 结论

本文提出了一种使用卡尔曼滤波算法，提高温度脉动仪测量准确性的方案。经过卡尔曼滤波后温度脉动仪采集的电压值更加平稳准确，设备的测量噪音能够控制在0.001 K 以内，噪音对大气折射率结构常数的影响量级为10-18m-2/3，远低于正常测量的数据，对大气折射率结构常数的测量影响极低，提高了温度脉动仪的测量准确度。利用研制的温度脉动仪对长春地区冬季的大气湍流进行全天监测，成功绘制大气折射率结构常数特性曲线，长春冬季大气折射率结构常数在10-12m-2/3～10-15m-2/3之间波动，在日落后日出前，大气湍流较弱，折射率结构常数在10-14m-2/3附近波动，日出后受太阳辐射影响，温度脉动量变大，光学湍流较强，折射率结构常数在10-13m-2/3附近波动，最高达到5 × 10-13m-2/3。全天变化趋势与先前的研究结果［10-11］具有很好的一致性。