基于蒙特卡洛法的砝碼校準(zhǔn)不確定度評定

韓志國，李鎖印，馮亞南，許曉青

(中國電子科技集團(tuán)公司第十三研究所，河北 石家莊050051)

0 引言

測量是以確定量值為目的的一組操作。測量結(jié)果受人員、測量系統(tǒng)、測量程序、環(huán)境以及其它的因素影響而存在測量不確定度。測量不確定度是與測量結(jié)果相關(guān)聯(lián)的、表征合理的、賦予被測量值分散性的參數(shù)。ISO/IEC 發(fā)布的Guide98 －3 －2008 《測量不確定度表示指南》 (GUM)給出了測量不確定度評定與表示的通用方法。2008 年，計量指南聯(lián)合委員會(JCGM)在GUM 方法的基礎(chǔ)上推出了補(bǔ)充件——采用蒙特卡洛法評定測量不確定度，將蒙特卡洛法作為對GUM 方法的重要補(bǔ)充。我國新頒布的JJF1059.2 －2012《采用蒙特卡洛法評定測量不確定度》就是依據(jù)該補(bǔ)充件進(jìn)行編制的。該規(guī)范規(guī)定了用蒙特卡洛法評定與表示測量不確定度的方法，并就GUM 中未涉及的概率分布傳播的問題提供指導(dǎo)，擴(kuò)大了測量不確定度評定與表示的適用范圍，同時該規(guī)范提供了檢查GUM 法是否適用的驗證方法。

1 蒙特卡洛法介紹

在測量不確定度的評定中采用的蒙特卡洛法(MCM)是一種通過重復(fù)采樣實現(xiàn)分布傳播的數(shù)值方法。與GUM 法利用線性化模型傳播不確定度的解析方法不同，MCM 通過對輸入量Xi的概率密度函數(shù)(PDF)離散采樣，由測量模型傳播輸入量的分布，計算獲得輸出量Y 的PDF 的離散采樣值，進(jìn)而由輸出量的離散分布數(shù)值直接獲取輸出量的最佳估計值、標(biāo)準(zhǔn)不確定度和包含區(qū)間。該輸出量的最佳估計值、標(biāo)準(zhǔn)不確定度和包含區(qū)間等特性的計算質(zhì)量隨PDF 采樣數(shù)增加可得到改善。

圖1 描述的是，由輸入量Xi(i =1，…，N)的PDF，通過模型傳播，給出輸出量Y 的PDF 的一個過程示意。圖1 列出了分別為相互獨(dú)立的正態(tài)分布、三角分布和正態(tài)分布的3 個輸入量，而輸出量的PDF 顯示為分布不對稱的情形。

圖1 輸入量獨(dú)立時分布傳播的描述

當(dāng)測量模型是線性的，且輸出量概率分布是正態(tài)分布時，GUM 方法可以提供準(zhǔn)確的不確定度和包含區(qū)間。但是，如果輸出量的結(jié)果是圖1 所示的不對稱的情形，那么采用GUM 方法進(jìn)行不確定度評定時可能就會得出不切實際的包含區(qū)間，此時采用蒙特卡洛法就是有效的替代方法。由于使用蒙特卡洛法進(jìn)行不確定度評定時需要對概率分布進(jìn)行隨機(jī)采樣，因此進(jìn)行不確定度評定時，應(yīng)利用合適的軟件進(jìn)行計算 (例如MATLAB)。

2 用蒙特卡洛法評定砝碼校準(zhǔn)不確定度

雖然GUM 法在許多情況下被認(rèn)為是非常適用的，但是確定是否滿足其所有應(yīng)用條件并不總是一件容易的事。由于MCM 的適用范圍比GUM 法更廣泛，建議用MCM 及GUM 法兩種方法同時進(jìn)行不確定度評定，并對結(jié)果進(jìn)行比較，如果比較結(jié)果較好，則GUM 法適用于此場合及今后足夠類似的情形，否則，應(yīng)考慮采用MCM 或者其它合適的替代方法。

下面以本實驗室2012 年在砝碼能力驗證中校準(zhǔn)100g 砝碼為例，介紹一下如何使用MCM 方法進(jìn)行不確定度評定。本次校準(zhǔn)使用的標(biāo)準(zhǔn)砝碼為E2等級100g砝碼，天平為TG328A 型機(jī)械天平。

2.1 建立測量模型

本次測量使用的是雙次替代法，用密度為ρr的標(biāo)準(zhǔn)砝碼R 對密度為ρt的砝碼進(jìn)行校準(zhǔn)，兩個砝碼標(biāo)稱值相同，需要在密度為ρa(bǔ)的空氣中進(jìn)行平衡配重。通常ρt和ρr不同，因此需考慮浮力的影響，應(yīng)用阿基米德定律，得到模型如下

式中:ms是密度為ρr靈敏度小砝碼質(zhì)量，用來測量機(jī)械天平的分度值。ΔI 為被測砝碼與標(biāo)準(zhǔn)砝碼標(biāo)尺讀數(shù)偏差;ΔIs為放上靈敏度小砝碼后，標(biāo)尺偏離的讀數(shù)。以折算質(zhì)量mct，mcr，mcs為變量，模型(1)轉(zhuǎn)變?yōu)?/p>

將公式(2)整理得

設(shè)Δmc= mct－ mnom為mct與標(biāo)稱質(zhì)量mnom=100 g的偏差。因此本例中使用的模型為

2.2 確定輸入量的概率密度分布函數(shù)(PDF)

可根據(jù)一系列測量值的分析，或根據(jù)某些歷史數(shù)據(jù)、校準(zhǔn)數(shù)據(jù)和專家判斷之類的信息所得到的科學(xué)判斷，為各輸入量設(shè)定PDF。各輸入量服從的分布、最佳估計值以及標(biāo)準(zhǔn)不確定度均列入表1。其中ρa(bǔ)0=1.2 kg/m3，是約定空氣密度，不包含不確定度。

表1 關(guān)于質(zhì)量校準(zhǔn)模型的輸入量及其概率密度函數(shù)

表1 中，標(biāo)準(zhǔn)砝碼的相關(guān)信息是通過檢定證書和JJG99 －2006 《砝碼檢定規(guī)程》得到。ΔI 與ΔIs的相關(guān)信息是通過使用雙次替代法重復(fù)測量10 次得到。砝碼的密度信息是通過JJG99 －2006 《砝碼檢定規(guī)程》中給出的典型材料密度表中查表得到?？諝饷芏圈補(bǔ)的信息是由JJG99 －2006 《砝碼檢定規(guī)程》中得到。

2.3 蒙特卡洛試驗次數(shù)

在規(guī)定的數(shù)值容差下MCM 所提供的結(jié)果所需的試驗次數(shù)跟輸出量的PDF “形狀”及包含概率有關(guān)，應(yīng)合理選擇蒙特卡洛試驗次數(shù)及樣本量的大小M，一般情況下取M=106，這樣通常會為輸出量提供95%包含區(qū)間，該包含區(qū)間長度被修約為1 或2 位有效十進(jìn)制數(shù)字。

2.4 程序

本程序是使用MATLAB 進(jìn)行編制的，具體如下:

clc;clear; %清除MATLAB 中的變量。

M=1000000; %設(shè)定蒙特卡洛試驗次數(shù)。

m0 =100; %砝碼標(biāo)稱質(zhì)量(單位:g)。

dr=0.04; %標(biāo)準(zhǔn)砝碼折算質(zhì)量修正值(單位:mg)。

ur=0.033; %標(biāo)準(zhǔn)砝碼不確定度(單位:mg)。

mr=normrnd ((m0* 10^3 +dr)，ur，1，M); %標(biāo)準(zhǔn)砝碼采樣值(單位:mg)。

ms=normrnd (10，0.002，1，M); %靈敏度小砝碼采樣值(單位:mg)。

dI= normrnd (4.2，0.2，1，M); % 讀數(shù)重復(fù)性 (單位:格)。

dIs = normrnd (99.6，0.3，1，M); % 靈敏度讀數(shù)重復(fù)性(單位:格)。

rou0 =normrnd (8，0.07，1，M); %標(biāo)準(zhǔn)砝碼材料密度(g/cm3)。

rou1 =normrnd (8，0.07，1，M); %被檢砝碼材料密度(g/cm3)。

ra=0.24* rand (1，M) +1.08; %空氣密度(mg/cm3)。

r0 =1.2; %標(biāo)準(zhǔn)空氣密度(mg/cm3)。

dm= (mr+dI.* ms. /dIs). * (rou0* 10^3 －ra). * (rou1*10^3 －r0). / ((rou1* 10^3 －ra).* (rou0* 10^3 －r0)) －m0

* 10^3;%雙次交換衡量法數(shù)學(xué)模型。

s=std (dm); %輸出量標(biāo)準(zhǔn)不確定度。

avg=mean (dm); %輸出量估計值。

y=sort (dm); %將試驗結(jié)果從小到大進(jìn)行排列。

y_ low=prctile (y，2.5); %包含概率為95%時，包含區(qū)間的下限值。

y_ high=prctile (y，97.5); %包含概率為95%時，包含區(qū)間的上限值。

2.5 報告結(jié)果

MATLAB 程序運(yùn)行得到如下結(jié)果:

被校砝碼折算質(zhì)量偏差為0.46 mg;標(biāo)準(zhǔn)不確定度為0.04 mg;包含概率為95%;包含區(qū)間為［0.38 mg，0.54 mg］。圖2 給出了上述程序得到的輸出量的近似概率分布。圖中，柱狀圖是輸出量的近似概率分布，實線是正態(tài)分布曲線，可見輸出量近似服從正態(tài)分布。

圖2 輸出量近似概率分布

3 蒙特卡洛法與其它不確定度評定方法的比較

3.1 GUM 法評定砝碼校準(zhǔn)不確定度

由于使用偏導(dǎo)的解析計算對公式(4)求偏導(dǎo)過于復(fù)雜，因此將公式(4)作近似處理，得

對上述模型求偏導(dǎo)，得到各輸入量靈敏系數(shù)如下:

則合成標(biāo)準(zhǔn)不確定度為

將表1 中輸入量的相關(guān)信息帶入上式，計算得uc=0.04mg，取包含概率為95%，得到如下結(jié)果:

被校砝碼折算質(zhì)量偏差為0.46 mg;標(biāo)準(zhǔn)不確定度為0.04 mg;包含概率為95%;包含區(qū)間為［0.38 mg，0.54 mg］。

3.2 JJG99 －2006 《砝碼檢定規(guī)程》中的方法評定砝碼校準(zhǔn)不確定度

規(guī)程中認(rèn)為砝碼的校準(zhǔn)不確定度由以下幾個方面引入:衡量過程，標(biāo)準(zhǔn)砝碼，空氣浮力修正，天平。各不確定度分量列入表2。

表2 測量不確定度分量一覽表

取包含概率為95%，得到如下結(jié)果:

被校砝碼折算質(zhì)量偏差為0.46 mg;標(biāo)準(zhǔn)不確定度為0.05 mg;包含概率為95%;包含區(qū)間為［0.36 mg，0.56 mg］。

三種方法得到的結(jié)果見表3。

表3 砝碼校準(zhǔn)計算結(jié)果

由三種方法得到的比較結(jié)果可知，MCM 與GUM得到的結(jié)果一致，而規(guī)程給出的不確定度略顯保守。這主要是由于某些不確定度分量被重復(fù)計入和計量人員對不確定度分量的不同理解程度而造成的。按照J(rèn)JF1059.2 －2012 《用蒙特卡洛法評定測量不確定度》中使用MCM 驗證GUM 的要求，GUM 法通過了驗證，在以后進(jìn)行砝碼校準(zhǔn)不確定度評定時，可以使用GUM法進(jìn)行評定。

4 結(jié)論

MCM 尤其適用于以下兩種情況:不宜對測量模型進(jìn)行線性化等近似的場合，在這種情況下，按GUM 確定輸出量的估計值和標(biāo)準(zhǔn)不確定度可能會變得不可靠;輸出量的概率密度函數(shù)較大程度地偏離正態(tài)分布或t 分布，例如分布明顯不對稱的場合，在這種情況下，可能會導(dǎo)致對包含區(qū)間或擴(kuò)展不確定度的估計不切實際。另外，如果GUM 法通過了MCM 的驗證，即GUM 法若明顯適用，則GUM 法依然是不確定度評定的主要方法。在實際工作中可以根據(jù)實際情況選用合適的方法進(jìn)行測量不確定度的評定。

［1］國家質(zhì)量監(jiān)督檢驗檢疫總局.JJF1059.2 －2012 用蒙特卡洛法評定測量不確定度［S］. 北京:中國計量出版社，2012.

［2］Evaluation of measurement data — Supplement 1 to the“Guide to the expression of uncertainty in measurement”— Propagation of distributions using aMonte Carlo method ［S］. ISO/IEC，2008.

［3］國家質(zhì)量監(jiān)督檢驗檢疫總局.JJG99 －2006 砝碼檢定規(guī)程［S］. 北京:中國計量出版社，2006.

［4］陳懷艷，曹蕓，韓潔. 基于蒙特卡羅法的測量不確定度評定［J］. 電子測量與儀器學(xué)報，2007，25.

［5］倪育才. 實用測量不確定度評定［M］. 北京:中國計量出版社，2010.