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

（1.中國計量大學 計量測試工程學院， 浙江 杭州 310018；2.湖州市生態(tài)環(huán)境局安吉分局 安吉縣生態(tài)環(huán)境監(jiān)測站， 浙江 安吉 313300）

0 引 言

在冷鏈運輸全過程中，每年因藥品疫苗、生鮮易腐食品的冷藏保溫問題造成的經(jīng)濟損失巨大，故針對冷鏈運輸過程中的溫度準確計量是保證產(chǎn)品質(zhì)量的關(guān)鍵因素。冷鏈溫度的傳統(tǒng)計量是將測溫傳感器采用送檢的手段參與計量校準活動。在該過程中，冷鏈溫度無法監(jiān)測[1]；而冷鏈溫度的在線計量是一種借助物聯(lián)網(wǎng)技術(shù)實現(xiàn)對溫度實時測量并保證量值準確、可靠的過程[2]。許多研究人員基于WiFi、NB-IoT、ZigBee 等技術(shù)設(shè)計了無線溫度監(jiān)測系統(tǒng)[3-5]，該系統(tǒng)具有精度高、性能穩(wěn)定、能耗低等特點，對于在線計量技術(shù)的發(fā)展具有重要意義。

利用在線計量技術(shù)可實現(xiàn)工況下冷鏈溫度的遠程校準，但目前溫度在線計量的精度無法估計，溫度測量的不確定度動態(tài)評定是主要難點。目前，許多研究人員基于蒙特卡洛、灰色模型等方法建立關(guān)于溫度測量的動態(tài)不確定度評定模型[6-8]，但由于測量信息利用不完整，評價結(jié)果存在偏頗性；而貝葉斯評定方法因具有融合歷史信息和樣本數(shù)據(jù)的特點，受到許多學者關(guān)注。

文獻[9]中提出一種基于貝葉斯統(tǒng)計推斷進行測量系統(tǒng)量值特性指標更新的方法，解決了因測量不確定度評定過后重復(fù)使用且無法反映測量系統(tǒng)中最新信息的問題；文獻[10]中研究關(guān)于貝葉斯統(tǒng)計理論并且可運用于動態(tài)測量和專家判斷的不確定度評定方法，在無信息先驗和共軛先驗的貝葉斯不確定度評定方法的基礎(chǔ)上，提出一種基于最大熵原理的貝葉斯不確定度優(yōu)化評定模型。由于溫度參數(shù)具有時變性和自相關(guān)性的特點，上述的評估方法不能直接運用在溫度領(lǐng)域。

馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo,MCMC）是一種可以構(gòu)造動態(tài)穩(wěn)定分布的方法，該方法考慮了隨機變量的自相關(guān)，但目前與測量不確定度動態(tài)評定相關(guān)的研究較少。

文獻[11]中研究關(guān)于MCMC 方法中的獨立抽樣和隨機游走抽樣的Metropolis-Hastings（M-H）算法，利用Matlab 程序來實現(xiàn)兩種抽樣算法并給出了詳細的程序?qū)嵤┝鞒?，分析了兩種抽樣方法的優(yōu)缺點。

基于上述文獻，本文以一種基于NB-IoT 技術(shù)的無線溫度記錄儀為實驗儀器，測量冷庫在工況下運行時的溫度，針對A 類不確定度評定，提出一種結(jié)合馬爾科夫鏈蒙特卡洛方法的數(shù)據(jù)融合貝葉斯估計方法，最終實現(xiàn)冷鏈溫度在線計量的不確定度動態(tài)估計，兼顧了評定過程的實時性與穩(wěn)定性。

1 冷鏈溫度在線計量原理

基于NB-IoT 技術(shù)的冷鏈溫度在線計量系統(tǒng)主要包括溫度數(shù)據(jù)采集、信號調(diào)節(jié)與轉(zhuǎn)換、數(shù)據(jù)傳輸與處理三部分，根據(jù)其實現(xiàn)的功能作用可以概括為感知層、調(diào)節(jié)層、傳輸層。溫度在線計量系統(tǒng)結(jié)構(gòu)如圖1 所示。采集終端將溫度信號轉(zhuǎn)換成電信號后，通過調(diào)制放大以及模數(shù)轉(zhuǎn)換輸出數(shù)字信號，利用NB-IoT 通信技術(shù)將溫度數(shù)據(jù)發(fā)送至云端服務(wù)器，數(shù)據(jù)監(jiān)測終端在線獲取溫度測量結(jié)果，并對測量結(jié)果的不確定度進行評估。

圖1 溫度在線計量系統(tǒng)結(jié)構(gòu)圖

2 傳統(tǒng)貝葉斯不確定度評定方法

傳統(tǒng)貝葉斯的不確定度評定方法以貝葉斯統(tǒng)計推斷原理為基礎(chǔ)，將歷史先驗信息和測量樣本信息利用貝葉斯公式融合，對后驗分布進行參數(shù)估計以實現(xiàn)對測量不確定度的評定。由于后驗分布中包含了總體、樣本、經(jīng)驗等與測量對象有關(guān)的大量信息，對后驗分布進行統(tǒng)計推斷相較于只對測量樣本進行分析來說，前者更為可靠?；谪惾~斯理論的測量不確定度評定模型[12]如下：

式中：θ為不確定度參數(shù)組成的向量；x為樣本數(shù)據(jù)信息；π(θ)表示先驗分布，是這些參數(shù)的聯(lián)合概率密度函數(shù)；L(θ|x)為樣本數(shù)據(jù)的似然函數(shù)；c為標準化常量，是描繪參數(shù)θ和數(shù)據(jù)x為某類模型時的可能性度量[13]；Θ為參數(shù)空間；h(θ|x)為后驗分布的概率密度函數(shù)。

由式（1）得后驗分布的期望值E(θ)和標準差S(θ)為：

3 結(jié)合MCMC 的改進貝葉斯方法

3.1 結(jié)合MCMC 的數(shù)據(jù)融合貝葉斯估計方法

在傳統(tǒng)貝葉斯評定方法中，先驗分布一經(jīng)設(shè)定便不再改變，這一步驟將放大不確定度評定過程中主觀因素的影響。針對這一問題，本文根據(jù)歷史數(shù)據(jù)設(shè)定先驗分布，將后驗分布作為下一次評定過程的先驗信息，不僅保證了不確定度評定過程的客觀性，還實現(xiàn)了信息的融合與更新。假設(shè)在連續(xù)時間段內(nèi)某一冷庫內(nèi)的溫度測量結(jié)果為{ }Tn,n≥0 ，服從正態(tài)分布。本文方法對該測量結(jié)果進行不確定動態(tài)評定的步驟如下：

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

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

4） 待新測量結(jié)果出現(xiàn)后，循環(huán)步驟2）、3），實現(xiàn)溫度測量過程的動態(tài)不確定度評定與更新。

溫度動態(tài)測量過程具有時變性、隨機性和自相關(guān)性的特點，其中自相關(guān)性具體表現(xiàn)為當前時刻的溫度測量結(jié)果只與上一時刻的測量結(jié)果相關(guān)，與之前時刻的測量結(jié)果均不相關(guān)。這一動態(tài)測量過程等同于馬爾科夫鏈過程，故本文結(jié)合MCMC 方法對后驗分布進行參數(shù)估計。

3.2 馬爾科夫鏈蒙特卡洛抽樣方法

設(shè){Xn,n≥0 }為描述測量過程中溫度測量結(jié)果的隨機過程，則該參數(shù)的馬爾科夫鏈模型可表示為：

式中：P(Xt+1=it+1|xt=it)稱為馬爾科夫鏈的轉(zhuǎn)移核，且對于隨機過程{Xt,t≥0 }，將來的狀態(tài){Xt+1=it+1}只與現(xiàn)在的狀態(tài){Xt=it}有關(guān)，而與過去的狀態(tài){Xk=ik,k≤n}無關(guān)[15]。這種用條件概率分布來表征隨機變量之間的自相關(guān)性的方法稱為馬爾科夫鏈過程。

MCMC 方法通過構(gòu)造滿足細致平衡條件的馬爾科夫鏈，結(jié)合蒙特卡洛數(shù)值模擬獲得溫度測量結(jié)果X的樣本，對樣本結(jié)果進行分析得到統(tǒng)計推斷結(jié)果。本文使用M-H 算 法 求 解MCMC 問 題。

3.3 Metropolis-Hastings 算法

Metropolis 抽樣方法[16-17]只適用于目標分布為對稱分布的情況，而M-H 算法在基于上述算法的前提下，對狀態(tài)轉(zhuǎn)移核進行了完善，使其也可以對非對稱情況下的目標分布進行抽樣。

具體實施步驟如下：

1） 構(gòu)造建議分布q(·|θ(t))；

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

3） 從建議分布q(·|θ(t))中產(chǎn)生擾動點θ(t+1)；

4） 計算接受轉(zhuǎn)移概率：

6） 循環(huán)步驟3）～5），直至馬爾科夫鏈收斂至平穩(wěn)狀態(tài)。

建議分布q(·|θ(t))為馬爾科夫鏈中狀態(tài)轉(zhuǎn)移過程的一個轉(zhuǎn)移規(guī)則，為保證經(jīng)馬爾科夫鏈收斂得到的后驗分布是唯一的平穩(wěn)分布，轉(zhuǎn)移規(guī)則須滿足不可約、非周期和正常返的特點。M-H 算法流程如圖2 所示。

圖2 Metropolis-Hastings 算法流程

4 實驗分析

本文在設(shè)定溫度為-20 ℃的冷庫內(nèi)進行實驗，待溫度穩(wěn)定后，在14:00—16:30 的現(xiàn)實時間段中進行離散采樣，使用無線溫度記錄儀1 min 采集一次數(shù)據(jù)，測量數(shù)據(jù)如表1 所示。

表1 冷庫內(nèi)溫度測量結(jié)果

使用GUM 法、傳統(tǒng)貝葉斯方法和本文方法分別對表1 測量數(shù)據(jù)評定其不確定度。將GUM 法引入貝塞爾公式，見式（11），對連續(xù)15 個時間節(jié)點的測量數(shù)據(jù)計算1 次標準不確定度；在傳統(tǒng)貝葉斯方法中，按經(jīng)驗設(shè)定先驗分布服從滿足均值為-19.6 ℃、標準差為0.16 的正態(tài)分布，引入最大似然估計方法計算似然函數(shù)，由式（3）、式（4）計算后驗分布均值和標準不確定度；在本文方法中，以前15 個時間節(jié)點的測量數(shù)據(jù)為先驗信息，引入最大似然估計方法計算似然函數(shù)，以馬爾科夫鏈蒙特卡洛的聯(lián)合概率密度函數(shù)為狀態(tài)轉(zhuǎn)移函數(shù)構(gòu)造馬氏鏈并使用M-H 方法進行抽樣，在收斂至穩(wěn)態(tài)分布后，計算抽樣結(jié)果的均值和標準不確定度。由三種方法分別計算得到的均值和標準不確定度，結(jié)果如表2 所示，其中均值計算結(jié)果如圖3a）所示，動態(tài)不確定度評定結(jié)果如圖3b）所示。

圖3 GUM 法、傳統(tǒng)貝葉斯、改進貝葉斯均值與動態(tài)不確定度對比折線圖

由圖3 可以看出，三種方法對于均值的估計具有一致性，但GUM 法評定結(jié)果波動較大，改進貝葉斯方法評定結(jié)果精準度更高。對于不確定度評定來說，GUM 法評定結(jié)果分散較大，且評定結(jié)果均大于其余兩種方法；傳統(tǒng)貝葉斯方法評定結(jié)果與GUM 法呈現(xiàn)出一致性，但穩(wěn)定性優(yōu)于GUM 法，評價結(jié)果波動性較小；改進貝葉斯方法評定結(jié)果呈下降趨勢，經(jīng)數(shù)據(jù)多次融合更新后，評定結(jié)果小于其余兩種方法，且更加平穩(wěn)。

5 結(jié) 語

本文提出一種結(jié)合馬爾科夫鏈蒙特卡洛方法的數(shù)據(jù)融合貝葉斯估計方法，該方法可以針對冷鏈溫度在線計量過程中的不確定度進行動態(tài)估計。為驗證該方法的可行性，本文對冷庫工況下的溫度測量數(shù)據(jù)進行動態(tài)評定。實驗結(jié)果表明，相比其他兩種方法，本文方法在穩(wěn)定的環(huán)境下評價結(jié)果精準度更高，穩(wěn)定性更好，符合實際測量情況；改進后的貝葉斯測量不確定度評定方法較改進前最高改善了41.67%，提高了測量結(jié)果的價值。

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