GUM 法和蒙特卡洛法在飼料中粗蛋白含量不確定度評定的比較

◆作者：余舒寧 高志存 崔惠娟＊ 石露莎 高蕓

◆單位：紅河州獸藥飼料監(jiān)察所

測量不確定度是表征賦予被測量值分散性的非負參數(shù)，是判定測量結(jié)果可信度的參數(shù)。測量不確定度的評定方法通常采用GUM 法；蒙特卡洛法（MCM）是一種新的評定不確定度方法，提供了比較和驗證GUM 法的一種途徑。本文以飼料中粗蛋白含量測定為例，采用GUM 法和MCM 法對其含量測定結(jié)果的不確定度進行評定，對兩種方法的評定結(jié)果進行比較驗證，得出相應結(jié)論。

1 材料和方法

1.1 儀器

AG-204 電子天平，瑞士梅特勒公司；Kjeltec 8200 凱氏定氮儀，F(xiàn)OSS 福斯分析儀器公司；Titrette 25mL 數(shù)字瓶口滴定器，德國BRAND 公司。

1.2 試劑

硫酸、無水硫酸鈉、硫酸銅均為分析純，0.1018mol/L 鹽酸滴定液，400g/L，氫氧化鈉，20g/L 硼酸，甲基紅—溴甲酚綠混合指示劑。

1.3 樣品及測定方法

樣品為仔豬配合飼料，按GB/T 6432-2018 飼料中粗蛋白的測定凱氏定氮法測定，用全量法，數(shù)據(jù)為兩平行樣測定結(jié)果，測量相關(guān)參數(shù)見表1。

2 GUM 法評定不確定度

2.1 測量模型與不確定度來源分析

式中：VHCl為試樣消耗鹽酸滴定液的體積，mL；V0為空白消耗鹽酸滴定液的體積，mL；CHCl為鹽酸滴定液的濃度，moL/L；為試樣質(zhì)量，g；0.014 為氮的毫摩爾質(zhì)量，g/mmol；6.25 為氮換算粗蛋白的平均系數(shù)；rep 為方法重復性因子。

由測量模型得知，該測量過程中不確定度有以下幾方面：重復測定的不確定度；試樣稱量的不確定度；鹽酸滴定液濃度的不確定度；鹽酸滴定液消耗體積的不確定度。

2.2 GUM 測量不確定度的評定

2.2.1 重復測定的不確定度u（rep）

根據(jù)表中的平行測定結(jié)果，用極差法來評定兩次測定的不確定度：

2.2.2 試樣稱量的不確定度u (m)，實驗天平最大允許誤差為±0.0002g，按均勻分布，則換算成標準不確定度：=0.00012g，稱量是兩次操作，u相對不確定度=0.0004。

表1 飼料中粗蛋白含量測定的數(shù)據(jù)

2.2.3 鹽酸滴定液濃度的不確定度u(c)，已評定出，直接引用，（0.1018±0.0002）mol/L，k=2。按正態(tài)分布，標準不確定度為：u相對不確定度為:0.0010。

2.2.4 鹽酸滴定液消耗體積的不確定度u(v-v0)，按扣除空白體積后計算。數(shù)字瓶口滴定器的最大允許誤差為±0.040mL，按三角分布，標準不確定度：；溫度對體積影響引入的不確定度，實驗溫度控制在（20±3）℃，水的膨脹系數(shù)為2.1×10-4℃-1，由此引起體積的變化為mL，按正態(tài)分布，u（V-V0）t=

兩項合成u（V-V0）=體積的相對不確定度為

2.2.5 飼料中粗蛋白含量的相對合成標準不確定度：

2.2.6 飼料中粗蛋白含量的合成標準不確定度：

uw=0.0037×19.6%=0.0725%

2.2.7 取k=2 時，則擴展不確定度：

U=k×uw=2×0.0725%=0.1%

3 MCM 法評定不確定度

3.1 MCM 法評定不確定度原理

MCM 是實現(xiàn)概率分布傳播的一種數(shù)值方法，它通過對輸入量Xi 的概率密度函數(shù)（PDF）離散抽樣，由測量模型傳播輸入量的分布，計算獲得輸出量Y 的PDF 的離散抽樣值，進而由輸出量的離散分布數(shù)值直接獲得輸出量的最佳估計值、標準不確定度和包含區(qū)間。MCM 中關(guān)鍵是對所有輸入量的PDF 進行隨機抽樣，用數(shù)學方法產(chǎn)生隨機數(shù)，Matlab 是一款功能強大的數(shù)值計算軟件，完全能勝任產(chǎn)生隨機數(shù)的功能。本文采用MathWorks公司的Matlab 軟件（2012b）編寫MCM 程序。

3.2 測量模型及經(jīng)驗參數(shù)

式中：α 為水的膨脹系數(shù)，2.1×10-4℃-1；t 為實測溫度，℃；t0為校準溫度。其他同GUM 模型。

據(jù)GUM 評定引用相應參數(shù)計算后列于表2。

3.3 MCM 實施的Matlab 程序

據(jù)周桃庚、劉存成（2014）的MCM 的Matlab 程序，略作改動及注釋后如下：

clear;

randn('state',0);%正態(tài)分布隨機數(shù)發(fā)生器初始狀態(tài);

rand('state',0);%均勻分布隨機數(shù)發(fā)生器初始狀態(tài);

mu=19.6;u=0.0725;k=1.96;%GUM 參數(shù);

dt=0.5*10^-4;%數(shù)值容差;

m=10^5;%試驗次數(shù);

format short;

hold off;

for h=1:100%循環(huán)次數(shù)上限

rep=normrnd (1.0000,0.0029,1,m);%重復性因子分布;

m_s=unifrnd(0.3838-0.0002,0.3838+0.0002,1,m);%樣品質(zhì)量分布g;

c_HCl=normrnd(0.1018,0.0001,1,m);%滴定用HCl的濃度分布mol/L;

a=8.44-0.04;

b=8.44+0.04;

v_HCl=a+(b-a)/2*(rand(1,m)+rand(1,m));%滴定消耗HCl 的體積分布mL；

表2 輸入量的大小，不確定度和設(shè)定的分布

t_t0=normrnd(0.0,1.5,1,m); %測定時的溫度分布℃；

alfa=2.1*10^-4; %水的膨脹系數(shù)；

y(h,:)=c_HCl.*v_HCl.*(1-alfa*t_t0).*0.014*6.25*100.*rep./m_s;%測量數(shù)學模型;

y_mean(h,:)=mean(y(h,:));%計算平均數(shù);

y_std(h,:)=std(y(h,:));%計算標準差（不確定度）;

y1=y (:);%將h-1 行m 列輸出量矩陣y 轉(zhuǎn)置為行向量y;

Y1=sort(y1);%對輸出量排序

for i=1:m*5/100;%輸出樣本的5%，約定概率為95%;

W1(:,i)=Y1(m*95/100+i)-Y1(i);%計算輸出量約定95%概率的包含區(qū)間，共m*95/100 個包含區(qū)間;

end;

n1=find(W1==min(W1));%找出包含區(qū)間最小值;

y_low(h,:)=Y1(n1);%包含區(qū)間左端點，2.5%;

y_high (h,:)=Y1(m*95/100 +n1);%包含區(qū)間右端點，97.5%

if h＞1;

y_mean_std=2*std (y_mean)/sqrt(h);%平均數(shù)的2 倍標準差

y_std_std=2*std (y_std)/sqrt(h);%標準差的2 倍標準差；

y_low_std=2*std (y_low)/sqrt(h);%包含區(qū)間左端點的2 倍標準差；

y_high_std=2*std(y_high)/sqrt(h);%包含區(qū)間右端點的2 倍標準差；

if (y_mean_std＜=dt/5&y_std_std＜=dt/5&y_low_std＜=dt/5&y_hi gh_std＜=dt/5)==1%進行穩(wěn)定性判斷；

break%達到數(shù)值容差1/5 要求，直接退出循環(huán)體；

end；

end；

end；

y=y(:);%將h-1 行m 列輸出量矩陣y 轉(zhuǎn)置為行向量y；

Y=sort(y);%對輸出量排序；

for r=1:m*(h-1)*5/100%建立循環(huán)體，輸出樣本5%包含區(qū)間；

W2(:,r)=Y(m*(h-1)*95/100+r)-Y (r);%輸出樣本95%包含區(qū)間；

end；

n=find(W2==min(W2))%找出包含區(qū)間最小值，共m*(h-1)*95/100 個包含區(qū)間；

d_low =abs(mu-k*u-Y(n));%計算左包含區(qū)間絕對差值；

d_high=abs (mu+k*u-Y(m*(h-1)*95/100+n));%計算右包含區(qū)間絕對差值；

圖1 飼料粗蛋白含量的概率密度分布

Y_mean=mean (y),Y_std=std(y),Y_low=Y(n),Y_high=Y(m*(h-1)*95/100+n);

d_low=d_low,d_high=d_high%顯示計算參數(shù)值;

lower=min(y);upper=max(y);%以下為繪制概率分布直方圖設(shè)置；

xc=lower: (upper-lower)/49:upper;

n=hist(y,xc);

bar (xc,n./(((upper-lower)/49)*m),1);

hh=findobj(gca,'Type','patch');

set (hh,'Facecolor','w','Edgecolor','b');

title('評估的概率密度分布');

xlabel('粗蛋白的含量（%）');

ylabel('概率密度');

hold on;

z=normrnd (mu,u,1,10^6);%GUM 正態(tài)分布

[Z,z]=ksdensity(z);%估計樣本的概率密度；

plot (z,Z*10^2,'g','linewidth',2.0);

legend('直方圖MCM 分布','正態(tài)曲線GUM 分布')；

hold off；

運行程序結(jié)果為：

n=299177，Y_mean=19.5881，Y_std =0.0715，Y_low =19.4540，Y_high =19.7230，d_low=0.0039，d_high=0.0191；

MCM 得到的平均值19.6%，標準不確定度為0.07%，95%概率的包含區(qū)間為[19.5%,19.7%]，擴展不確定度U95=0.07%×2=0.1%。

表3 GUM 法與MCM 法不確定度評定結(jié)果比較

4 GUM 法與MCM 法不確定度評定結(jié)果的比較驗證

對GUM 法和MCM 法評定結(jié)果進行比較，當dlow=｜y-u(y)-ylow｜及dhigh=｜y+u(y)-yhigh｜均小于δ 時，GUM 法可通過驗證（JJG 1036-2008，JJG 196-2006）。GUM 法與MCM 法不確定度評定結(jié)果如表3 所示。

由表3可知，MCM 法對GUM 法不確定度框架通過了驗證，且MCM 與GUM 兩者評定的結(jié)果一致，達到預期。

5 分析與結(jié)論

通過上述評定過程我們可以看出，飼料中粗蛋白含量測定結(jié)果的不確定度主要來源于重復測定的不確定度和鹽酸滴定液消耗體積的不確定度。以平行樣數(shù)據(jù)（n=2）評定不確定度，結(jié)果準確性會差一些。要準確評定不確定度，應按標準要求，試驗平行樣不少于10 次（n≥10），用貝塞爾法計算標準差，評定不確定度。但現(xiàn)有檢測技術(shù)要求，化學實驗中檢測樣品，基本都是要求檢測平行樣（n=2），這種用平行樣評定不確定度的方法，可粗略評估出不確定度，作為初步的判斷，在工作中有一定的實用價值。

MCM 評定不確定度雖然是純數(shù)學方法，但還是以GUM 法為基礎(chǔ)，所有的經(jīng)驗參數(shù)值，都是以實驗數(shù)據(jù)為依據(jù)，進行仿真數(shù)學模擬，結(jié)果可信度是有保障的。