冷鏈溫度在線計量的不確定度動態評定

周易明， 程銀寶， 鄭 重， 江文松， 李亞茹， 張佳璇

（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%，提高了測量結果的價值。

注：本文通訊作者為李亞茹。