冷链温度在线计量的不确定度动态评定

周易明， 程银宝， 郑 重， 江文松， 李亚茹， 张佳璇

（1.中国计量大学 计量测试工程学院， 浙江 杭州 310018；2.湖州市生态环境局安吉分局 安吉县生态环境监测站， 浙江 安吉 313300）

0 引 言

在冷链运输全过程中，每年因药品疫苗、生鲜易腐食品的冷藏保温问题造成的经济损失巨大，故针对冷链运输过程中的温度准确计量是保证产品质量的关键因素。冷链温度的传统计量是将测温传感器采用送检的手段参与计量校准活动。在该过程中，冷链温度无法监测[1]；而冷链温度的在线计量是一种借助物联网技术实现对温度实时测量并保证量值准确、可靠的过程[2]。许多研究人员基于WiFi、NB-IoT、ZigBee 等技术设计了无线温度监测系统[3-5]，该系统具有精度高、性能稳定、能耗低等特点，对于在线计量技术的发展具有重要意义。

利用在线计量技术可实现工况下冷链温度的远程校准，但目前温度在线计量的精度无法估计，温度测量的不确定度动态评定是主要难点。目前，许多研究人员基于蒙特卡洛、灰色模型等方法建立关于温度测量的动态不确定度评定模型[6-8]，但由于测量信息利用不完整，评价结果存在偏颇性；而贝叶斯评定方法因具有融合历史信息和样本数据的特点，受到许多学者关注。

文献[9]中提出一种基于贝叶斯统计推断进行测量系统量值特性指标更新的方法，解决了因测量不确定度评定过后重复使用且无法反映测量系统中最新信息的问题；文献[10]中研究关于贝叶斯统计理论并且可运用于动态测量和专家判断的不确定度评定方法，在无信息先验和共轭先验的贝叶斯不确定度评定方法的基础上，提出一种基于最大熵原理的贝叶斯不确定度优化评定模型。由于温度参数具有时变性和自相关性的特点，上述的评估方法不能直接运用在温度领域。

马尔科夫链蒙特卡洛（Markov Chain Monte Carlo,MCMC）是一种可以构造动态稳定分布的方法，该方法考虑了随机变量的自相关，但目前与测量不确定度动态评定相关的研究较少。

文献[11]中研究关于MCMC 方法中的独立抽样和随机游走抽样的Metropolis-Hastings（M-H）算法，利用Matlab 程序来实现两种抽样算法并给出了详细的程序实施流程，分析了两种抽样方法的优缺点。

基于上述文献，本文以一种基于NB-IoT 技术的无线温度记录仪为实验仪器，测量冷库在工况下运行时的温度，针对A 类不确定度评定，提出一种结合马尔科夫链蒙特卡洛方法的数据融合贝叶斯估计方法，最终实现冷链温度在线计量的不确定度动态估计，兼顾了评定过程的实时性与稳定性。

1 冷链温度在线计量原理

基于NB-IoT 技术的冷链温度在线计量系统主要包括温度数据采集、信号调节与转换、数据传输与处理三部分，根据其实现的功能作用可以概括为感知层、调节层、传输层。温度在线计量系统结构如图1 所示。采集终端将温度信号转换成电信号后，通过调制放大以及模数转换输出数字信号，利用NB-IoT 通信技术将温度数据发送至云端服务器，数据监测终端在线获取温度测量结果，并对测量结果的不确定度进行评估。

图1 温度在线计量系统结构图

2 传统贝叶斯不确定度评定方法

传统贝叶斯的不确定度评定方法以贝叶斯统计推断原理为基础，将历史先验信息和测量样本信息利用贝叶斯公式融合，对后验分布进行参数估计以实现对测量不确定度的评定。由于后验分布中包含了总体、样本、经验等与测量对象有关的大量信息，对后验分布进行统计推断相较于只对测量样本进行分析来说，前者更为可靠。基于贝叶斯理论的测量不确定度评定模型[12]如下：

式中：θ为不确定度参数组成的向量；x为样本数据信息；π(θ)表示先验分布，是这些参数的联合概率密度函数；L(θ|x)为样本数据的似然函数；c为标准化常量，是描绘参数θ和数据x为某类模型时的可能性度量[13]；Θ为参数空间；h(θ|x)为后验分布的概率密度函数。

由式（1）得后验分布的期望值E(θ)和标准差S(θ)为：

3 结合MCMC 的改进贝叶斯方法

3.1 结合MCMC 的数据融合贝叶斯估计方法

在传统贝叶斯评定方法中，先验分布一经设定便不再改变，这一步骤将放大不确定度评定过程中主观因素的影响。针对这一问题，本文根据历史数据设定先验分布，将后验分布作为下一次评定过程的先验信息，不仅保证了不确定度评定过程的客观性，还实现了信息的融合与更新。假设在连续时间段内某一冷库内的温度测量结果为{ }Tn,n≥0 ，服从正态分布。本文方法对该测量结果进行不确定动态评定的步骤如下：

1） 设定先验分布π(θ)～N(u^0,σ^0)，引入最大似然估计方法[14]，公式如下：

3） 使用MCMC 方法计算得到后验分布均值u2及标准差σ2，并更新先验分布π(θ)～N(u2,σ2)；

4） 待新测量结果出现后，循环步骤2）、3），实现温度测量过程的动态不确定度评定与更新。

温度动态测量过程具有时变性、随机性和自相关性的特点，其中自相关性具体表现为当前时刻的温度测量结果只与上一时刻的测量结果相关，与之前时刻的测量结果均不相关。这一动态测量过程等同于马尔科夫链过程，故本文结合MCMC 方法对后验分布进行参数估计。

3.2 马尔科夫链蒙特卡洛抽样方法

设{Xn,n≥0 }为描述测量过程中温度测量结果的随机过程，则该参数的马尔科夫链模型可表示为：

式中：P(Xt+1=it+1|xt=it)称为马尔科夫链的转移核，且对于随机过程{Xt,t≥0 }，将来的状态{Xt+1=it+1}只与现在的状态{Xt=it}有关，而与过去的状态{Xk=ik,k≤n}无关[15]。这种用条件概率分布来表征随机变量之间的自相关性的方法称为马尔科夫链过程。

MCMC 方法通过构造满足细致平衡条件的马尔科夫链，结合蒙特卡洛数值模拟获得温度测量结果X的样本，对样本结果进行分析得到统计推断结果。本文使用M-H 算 法 求 解MCMC 问 题。

3.3 Metropolis-Hastings 算法

Metropolis 抽样方法[16-17]只适用于目标分布为对称分布的情况，而M-H 算法在基于上述算法的前提下，对状态转移核进行了完善，使其也可以对非对称情况下的目标分布进行抽样。

具体实施步骤如下：

1） 构造建议分布q(·|θ(t))；

2） 给定初始值θ(0)；

3） 从建议分布q(·|θ(t))中产生扰动点θ(t+1)；

4） 计算接受转移概率：

6） 循环步骤3）～5），直至马尔科夫链收敛至平稳状态。

建议分布q(·|θ(t))为马尔科夫链中状态转移过程的一个转移规则，为保证经马尔科夫链收敛得到的后验分布是唯一的平稳分布，转移规则须满足不可约、非周期和正常返的特点。M-H 算法流程如图2 所示。

图2 Metropolis-Hastings 算法流程

4 实验分析

本文在设定温度为-20 ℃的冷库内进行实验，待温度稳定后，在14:00—16:30 的现实时间段中进行离散采样，使用无线温度记录仪1 min 采集一次数据，测量数据如表1 所示。

表1 冷库内温度测量结果

使用GUM 法、传统贝叶斯方法和本文方法分别对表1 测量数据评定其不确定度。将GUM 法引入贝塞尔公式，见式（11），对连续15 个时间节点的测量数据计算1 次标准不确定度；在传统贝叶斯方法中，按经验设定先验分布服从满足均值为-19.6 ℃、标准差为0.16 的正态分布，引入最大似然估计方法计算似然函数，由式（3）、式（4）计算后验分布均值和标准不确定度；在本文方法中，以前15 个时间节点的测量数据为先验信息，引入最大似然估计方法计算似然函数，以马尔科夫链蒙特卡洛的联合概率密度函数为状态转移函数构造马氏链并使用M-H 方法进行抽样，在收敛至稳态分布后，计算抽样结果的均值和标准不确定度。由三种方法分别计算得到的均值和标准不确定度，结果如表2 所示，其中均值计算结果如图3a）所示，动态不确定度评定结果如图3b）所示。

图3 GUM 法、传统贝叶斯、改进贝叶斯均值与动态不确定度对比折线图

由图3 可以看出，三种方法对于均值的估计具有一致性，但GUM 法评定结果波动较大，改进贝叶斯方法评定结果精准度更高。对于不确定度评定来说，GUM 法评定结果分散较大，且评定结果均大于其余两种方法；传统贝叶斯方法评定结果与GUM 法呈现出一致性，但稳定性优于GUM 法，评价结果波动性较小；改进贝叶斯方法评定结果呈下降趋势，经数据多次融合更新后，评定结果小于其余两种方法，且更加平稳。

5 结 语

本文提出一种结合马尔科夫链蒙特卡洛方法的数据融合贝叶斯估计方法，该方法可以针对冷链温度在线计量过程中的不确定度进行动态估计。为验证该方法的可行性，本文对冷库工况下的温度测量数据进行动态评定。实验结果表明，相比其他两种方法，本文方法在稳定的环境下评价结果精准度更高，稳定性更好，符合实际测量情况；改进后的贝叶斯测量不确定度评定方法较改进前最高改善了41.67%，提高了测量结果的价值。

注：本文通讯作者为李亚茹。