電廠燃煤虛擬采樣及發(fā)熱量均值不確定度評定

劉福國 國欽光，殷炳毅，王守恩

(1.國網(wǎng)山東省電力公司電力科學(xué)研究院，山東 濟南 250003；2.山東電力研究院，山東 濟南 250003；3.華電濰坊發(fā)電有限公司，山東 濰坊 261000)

1 引 言

檢測發(fā)電機組效率時，鍋爐燃煤發(fā)熱量還無法連續(xù)測量，只能通過采集入爐煤樣品，在實驗室化驗分析得到[1]。用少量樣品煤的化驗數(shù)據(jù)表示檢測周期內(nèi)入爐煤發(fā)熱量，其不確定性來自樣品采集和實驗室分析等2個環(huán)節(jié)[2]，采樣環(huán)節(jié)引起的不確定度占95%以上[3]。文獻[4]研究了鍋爐采樣樣品煤發(fā)熱量的概率分布特性，對單個樣品發(fā)熱量的不確定度進行了評定。在實際生產(chǎn)中，機組效率常用供電煤耗表示，電廠日常供電煤耗檢測周期有時長達1個月、1個季度甚至1年。檢測周期內(nèi)的燃煤發(fā)熱量一般采用多個樣品煤發(fā)熱量的平均值[5]。因此，多次采樣樣品發(fā)熱量均值的不確定度評定對機組效率檢測具有重要意義。

μx=μ

(1)

(2)

當(dāng)總體N為正態(tài)分布或樣本容量n>30時，樣本均值的抽樣分布可作為正態(tài)分布，正態(tài)分布的參數(shù)μx、σx根據(jù)式(1)、式(2)計算，故可對n個樣品發(fā)熱量均值的不確定度進行評定。

但當(dāng)總體N為非正態(tài)分布且樣本容量n<30時，樣本均值的抽樣分布為非正態(tài)分布，此時，如何對n個樣品發(fā)熱量均值的不確定度進行評定，是本文要解決的問題。

2 虛擬采樣及發(fā)熱量均值不確定度

非正態(tài)總體N通常來自實際采樣實驗，為已知數(shù)據(jù)。要對n個樣品發(fā)熱量均值的不確定度進行評定，首先獲取樣本均值的抽樣樣本，傳統(tǒng)方法[9]是從總體的N個發(fā)熱量中隨機抽取一個容量為n的樣本，計算該樣本均值，采用同樣方式進行多次抽樣，得到樣本均值的抽樣樣本，利用該樣本對均值的不確定度進行評定。從總體的N個發(fā)熱量中抽取容量為n的樣本，重復(fù)抽樣總計可得到Nn個樣本，因此樣本均值抽樣樣本的最大容量為Nn。

獲得樣本均值的抽樣樣本的另一種方法是基于蒙特卡洛法[10，11]的隨機采樣實驗，與傳統(tǒng)抽樣方法相比，隨機采樣試驗獲得的虛擬樣本容量是任意的，采用隨機采樣試驗對n個樣品發(fā)熱量均值的不確定度進行評定的步驟為：

(1) 計算總體N的發(fā)熱量概率密度。

(2) 將總體N發(fā)熱量概率密度表示成多項式函數(shù)。

(3) 建立虛擬采樣系統(tǒng)。按照總體N的概率密度產(chǎn)生隨機數(shù)，作為虛擬采樣樣品的發(fā)熱量。這種利用給定概率密度函數(shù)進行隨機樣抽的系統(tǒng)，看作是樣品發(fā)熱量的虛擬采樣系統(tǒng)。

(4) 利用虛擬采樣系統(tǒng)進行n次采樣，得到n個虛擬樣品的發(fā)熱量數(shù)據(jù)，計算發(fā)熱量的平均值，得到第1個發(fā)熱量均值；采用同樣的方法得到m個發(fā)熱量均值，作為發(fā)熱量均值的抽樣樣本。

(5) 計算發(fā)熱量均值樣本的平均值、標準差及概率密度分布，在給定置信概率下求取包含系數(shù)，得到n個樣品發(fā)熱量均值擴展不確定度。

以某電廠鍋爐為例，對入爐煤發(fā)熱量均值不確定度進行評定。

2.1 樣本N的獲得

表1為某電廠2017年5～7月入爐煤實際采樣試驗樣本N，該樣本包含276個樣品煤的發(fā)熱量[12]。這些發(fā)熱量數(shù)據(jù)既包含實驗室化驗分析引起的不確定性，又包含采制樣本環(huán)節(jié)引起的不確定性。研究表明，實驗室分析環(huán)節(jié)引起的不確定度約為0.05 MJ/kg[13～15]，而采制樣本環(huán)節(jié)引起的不確定度為0.5～2.5 MJ/kg[2，4，12]。因此，忽略實驗室分析環(huán)節(jié)的不確定度，可認為表1中發(fā)熱量數(shù)據(jù)的不確定性是由采制樣本引起的。利用表1的數(shù)據(jù)對n次采樣發(fā)熱量均值的不確定度進行評定。

表1給出的樣本平均值μN=19.686 MJ/kg，標準差σN=0.933 MJ/kg。

表1 總體樣本N的發(fā)熱量Tab.1 Calorific value of sample N MJ/kg

2.2 樣本N的概率密度函數(shù)

總樣本N中，發(fā)熱量最小值為16.894 MJ/kg，最大值為21.217 MJ/kg。將此區(qū)間等分成若干個子區(qū)間，統(tǒng)計每個子區(qū)間內(nèi)的樣品個數(shù)，計算落入?yún)^(qū)間的概率，概率密度是概率的變化率。利用Matlab統(tǒng)計工具箱中的hist函數(shù)可以統(tǒng)計落入每個子區(qū)間內(nèi)的樣品數(shù)，利用ksdensity函數(shù)可以對樣本N的概率密度進行核心平滑估計[16]。由圖1給出的通過區(qū)間統(tǒng)計得到的概率密度的直方圖以及概率密度平滑估計可以看出，樣本發(fā)熱量不服從正態(tài)分布。概率密度平滑估計結(jié)果可用8階多項式函數(shù)f表示為：

圖1 總體樣本N的概率密度Fig.1 Probability density of population sample N

(3)

對于總體樣本N，式(3)所示的多項式函數(shù)由分段函數(shù)f1和f2組成，在不同區(qū)間上，函數(shù)f1和f2的系數(shù)見表2。

表2 多項式函數(shù)f的系數(shù)Tab.2 Coefficient of polynomial function f

2.3 虛擬采樣

根據(jù)概率密度函數(shù)式(3)，采用挑選抽樣方法產(chǎn)生發(fā)熱量樣本，這種方法由馮·諾依曼在1951年提出[17]。如圖2所示，當(dāng)已知概率密度函數(shù)f(Q)時，利用挑選抽樣法產(chǎn)生發(fā)熱量隨機數(shù)的主要步驟為：

(1) 在發(fā)熱量取值區(qū)間[a，b]內(nèi)，產(chǎn)生均勻隨機數(shù)Q。

Q=(b-a)r1+a

(4)

式中：r1為區(qū)間[0，1]均勻隨機數(shù)，r1∈[0，1]。

(2) 在[0，fmax]區(qū)間內(nèi)，產(chǎn)生均勻隨機數(shù)f。

f=r2fmax

(5)

式中：r2為區(qū)間[0，1]均勻隨機數(shù)，r2∈[0，1]；fmax為函數(shù)f(Q)在區(qū)間[a，b]上的最大值。

(3) 當(dāng)f≤f(Q)時，接受Q為所選的隨機數(shù)，否則，返回第(1)步，重新抽取1對(Q，f)。

根據(jù)圖2，挑選樣抽法的幾何解釋為：隨機選取位于矩形abcd內(nèi)的點(Q，f)，選擇位于曲線f(Q)下面的點，它們對應(yīng)的發(fā)熱量Q服從概率密度為f(Q)的分布。

圖2 根據(jù)給定概率密度函數(shù)f(Q)產(chǎn)生隨機數(shù)QFig.2 Generating random numbers Q based on a given probability density function f(Q)

2.4 均值樣本的獲得及評定

要對n次采樣發(fā)熱量均值的不確定度進行評定，首先利用第2.3節(jié)的虛擬采樣技術(shù)獲得均值樣本。按式(3)給定的概率密度函數(shù)，產(chǎn)生n個發(fā)熱量隨機數(shù)，計算它們的平均值，作為發(fā)熱量均值的第1個采樣數(shù)據(jù)。用相同方式得到m個發(fā)熱量均值的采樣數(shù)據(jù)，作為發(fā)熱量均值樣本。

取n=4和n=50，采用虛擬采樣方法分別得到發(fā)熱量均值樣本，樣本容量m取為m=10 000。對樣本進行統(tǒng)計計算，n=4時，樣本平均值μ4=19.689 MJ/kg，標準差σ4=0.480 MJ/kg；n=50時，樣本平均值μ50=19.689 MJ/kg，標準差σ50=0.135 MJ/kg。利用ksdensity函數(shù)對2個樣本的概率密度進行核心平滑估計，n=4時，結(jié)果見圖3，其概率密度可擬合成式(3)所示的多項式函數(shù)f(Q)，限于篇幅，這里不再給出多項式系數(shù)；n=50時，區(qū)間統(tǒng)計得到的概率密度直方圖以及函數(shù)ksdensity的平滑估計結(jié)果見圖4，發(fā)熱量均值樣本非常接近正態(tài)分布，圖4中還給出了正態(tài)分布N(19.686，0.1352)的概率密度，可以看出，它和平滑估計結(jié)果吻合良好。

圖3 n=4時發(fā)熱量平均值的概率密度分布Fig.3 Probabilistic density distribution of the mean calorific value at n=4

圖4 n=50時發(fā)熱量平均值的概率密度分布與正態(tài)分布的對比Fig.4 Comparison of probability density of average calorific value with normal distribution at n=50

如圖3，對n=4時發(fā)熱量均值的不確定度進行評定，就是求給定的置信概率ε下的包含因子k，根據(jù)置信概率的定義[18]得到關(guān)于包含因子k的方程如式(6)。

求解式(6)，可得到包含因子k。求解時取置信概率ε=0.95，式(6)中μ4=19.689，σ4=0.480。采用迭代法求解式(6)，得到k=1.934。因此，n=4時發(fā)熱量均值的擴展不確定度為kσ4=0.928 MJ/kg，此時發(fā)熱量均值結(jié)果可表示為19.689±0.928 MJ/kg。

(6)

根據(jù)虛擬采樣方法，同樣可對n=50時發(fā)熱量均值的不確定度進行評定。利用Matlab的jbtest函數(shù)對n=50時發(fā)熱量均值樣本進行正態(tài)性檢驗[19]。結(jié)果表明，均值樣本服從正態(tài)分布，置信概率ε=0.95對應(yīng)的包含因子k=1.96[20,21]，發(fā)熱量均值的擴展不確定度為kσ50=0.265 MJ/kg。因此，n=50時發(fā)熱量均值的采樣化驗結(jié)果可表示為19.689±0.265 MJ/kg。

3 結(jié) 論

(1) 在發(fā)電機組效率檢測時，鍋爐燃煤多次采樣樣品發(fā)熱量均值的不確定度評定對于機組效率檢測具有重要意義。

(2) 已知燃煤發(fā)熱量樣本，要對n個樣品發(fā)熱量均值的不確定度進行評定，當(dāng)發(fā)熱量樣本正態(tài)分布或n>30時，則n個樣品發(fā)熱量均值服從正態(tài)分布，該正態(tài)分布的參數(shù)、根據(jù)式(1)和式(2)計算，據(jù)此可對n個樣品發(fā)熱量均值的不確定度進行評定。

(3) 當(dāng)發(fā)熱量樣本不服從正態(tài)分布且n<30時，要對n個樣品發(fā)熱量均值的不確定度進行評定，可采用挑選樣抽法進行虛擬采樣，得到n個樣品發(fā)熱量的均值樣本，然后根據(jù)該均值樣本的概率密度分布對發(fā)熱量均值的不確定度評定。

(4) 采用虛擬采樣均值樣本對n個樣品發(fā)熱量均值的不確定度評定結(jié)果與根據(jù)中心極限定理得到的結(jié)論一致。