基于擬蒙特卡羅方法的虛擬儀器不確定度評定

石艷軍　翁煒霖

摘要：針對蒙特卡羅方法在虛擬儀器不確定度評定時存在收斂速度慢、計算結(jié)果不穩(wěn)定、評定效果受仿真次數(shù)影響的缺陷，該文提出使用擬蒙特卡羅方法替代蒙特卡羅方法進行虛擬儀器的不確定評定。其核心是產(chǎn)生低偏差數(shù)列，通過轉(zhuǎn)換成與虛擬儀器不確定度源相同的概率分布，進行仿真模擬。文中以數(shù)字稱重傳感系統(tǒng)為例，介紹了具體的評定步驟與流程。通過與GUM方法、蒙特卡羅方法的比較，驗證了該方法的可取得較好的評定效果，可廣泛應(yīng)用了虛擬儀器的不確定度評定。

關(guān)鍵詞： 虛擬儀器；不確定度評定；蒙特卡羅方法；擬蒙特卡羅方法；低偏差數(shù)列

中圖分類號：TP311 文獻標識碼：A 文章編號：1009-3044（2015）12-0230-04

Evaluation of Virtual Instrument Measurement Uncertainty Based on Quasi Monte-Carlo Method

SHI Yan-jun，WEI Wei-lin

（ Guangzhou Meteorological Ssatellite Ground Station， Guangzhou 510640， China）

Abstract： Because of the limitations of low convergence， unstable results and sample quantities affecting evaluation quality of Monte-Carlo method， quasi Monte-Carlo method was used to assess virtual instrument uncertainty instead of Monte-Carlo method. The key point was how to produce low discrepancy numbers. These numbers were then transformed into the desirable distributed random numbers which conformed to uncertainty sources of virtual instrument. The paper take a digital weighing system for example， introducing the evaluation steps in detail. The result demonstrated that quasi Monte-Carlo is a better method on uncertainty evaluation， through compared with GUM and Monte-Carlo method. This method can widespread use on evaluation of virtual instrument measurement uncertainty.

Key words： virtual instrument； evaluation of uncertainty； Meonte-Carlo method； quasi Monte-Carlo method； low discrepancy numbers

20世紀80年代在美國興起的虛擬儀器以其高度的集成性、配置的靈活性和經(jīng)濟性等諸多優(yōu)點迅速成為研究的熱點，并廣泛應(yīng)用于測量控制領(lǐng)域[1]。測量不確定度是用以表征測量結(jié)果不能確定的程度，是衡量測量質(zhì)量的重要指標，ISO 17025[2]標準明確規(guī)定：校準實驗室或進行自校準的檢測實驗室，對所有的校準和各種校準類型都應(yīng)具有并應(yīng)用評定測量不確定度的程序。

《測量不確定度評定指南》（GUM）提出了A類和B類兩種測量不確定度的評定，但由于虛擬儀器的復(fù)雜性和多樣性，其數(shù)學(xué)表達式不具有GUM所設(shè)定的顯示解析、可導(dǎo)、近線性的適用條件，收到相關(guān)性的限制[3]。針對GUM的不足，國內(nèi)外學(xué)者對虛擬儀器的不確定度進行了廣泛的研究，出現(xiàn)了很多新方法。Salvatore Nuccio等[4]首先提出用數(shù)值仿真方法對虛擬儀器的AD轉(zhuǎn)換模塊進行不確定度評定，并驗證了其準確性。Ghiani E，Locci N等[5-6]提出用蒙特卡羅的數(shù)值方法取代GUM對虛擬儀器采樣輸入信號的數(shù)字處理，并開發(fā)了一套虛擬儀器不確定度的自動評估軟件。王中宇等[7]應(yīng)用徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)構(gòu)建數(shù)學(xué)模型，使用差分方程計算誤差傳播系數(shù)，應(yīng)用灰色系統(tǒng)理論解決了小樣本虛擬儀器的測量。袁敏等[8]對蒙特卡羅方法進行了改進，評定了計算機舍入誤差的影響，并應(yīng)用于數(shù)字稱重傳感器。詹惠琴[3]、王偉[9]等也對蒙特卡羅在虛擬儀器測量不確定度上進行了相關(guān)研究。

然而蒙特卡羅方法在模擬仿真時，收斂速度較慢，計算結(jié)果不穩(wěn)定，且受到采樣次數(shù)的限制；由于偽隨機數(shù)列的隨機性過強而均勻性不足，故使蒙特卡羅方法的應(yīng)用效果受到影響。因而有專家[10]提出用分布更加均勻地低偏差數(shù)列代替?zhèn)坞S機數(shù)列進行仿真模擬，稱為擬蒙特卡羅方法。目前還未見利用擬蒙特卡羅方法進行虛擬儀器測量不確定度的研究。本文采用擬蒙特卡羅方法對虛擬儀器不確定度評定，其核心是用低偏差的擬蒙特卡羅數(shù)列進行模擬仿真。

1偽隨機數(shù)列和擬隨機數(shù)列

1.1偽隨機數(shù)列

蒙特卡羅方法的基本思路是利用概率事件仿真技術(shù)產(chǎn)生服從特定概率分布的隨機數(shù)列來模擬虛擬儀器測量鏈中的隨機誤差源，從而獲取一系列測量結(jié)果，并用統(tǒng)計方法得出虛擬儀器的測量不確定度。仿真效果的關(guān)鍵在于隨機數(shù)列的質(zhì)量，最經(jīng)典的方法是用同余法產(chǎn)生均勻分布的隨機數(shù)，然后根據(jù)相應(yīng)的變換得到其他分布的隨機數(shù)列。生成的隨機數(shù)列要經(jīng)過隨機性檢驗，但該方法生成的隨機數(shù)列不是真正的隨機數(shù)列，稱其為偽隨機數(shù)列。偽隨機數(shù)列隨機性過強而均勻性不足，而虛擬儀器的不確定度源絕大部分屬于均勻分布，故仿真效果受到影響。

1.2 擬隨機數(shù)列

比偽隨機數(shù)列分布更為均勻的是低偏差數(shù)列，即擬隨機數(shù)列。擬蒙特卡羅方法正是采用擬隨機數(shù)列代替?zhèn)坞S機數(shù)列進行仿真模擬。擬隨機數(shù)列是基于確定性的點列而不是隨機點列，較常見的擬隨機數(shù)列有halton數(shù)列和sobol數(shù)列。分別用偽隨機數(shù)列和sobol數(shù)列產(chǎn)生1000個均勻分布的點列，如圖1和圖2所示，可見sobol數(shù)列比偽隨機數(shù)列更加均勻的分布于樣本空間。

此外，擬隨機數(shù)列對比偽隨機數(shù)列還有如下優(yōu)點：①偽隨機數(shù)列的收斂半徑為[ON-0.5]，而擬隨機數(shù)列[11]的收斂半徑為[OlogNdN]，可見收斂速度更快，仿真速度相應(yīng)提高。②蒙特卡羅模擬精度受偽隨機數(shù)列采樣次數(shù)的影響，如表1所示，對均勻分布和正態(tài)分布的偽隨機數(shù)采用不同的抽樣次數(shù)N，獲得其期望與方差，經(jīng)過10次這樣的統(tǒng)計，得到其期望與方差的變動量，可看出，抽樣次數(shù)越多，期望與方差的變動量越小。而擬蒙特卡羅方法采用的是擬隨機數(shù)列，擬隨機數(shù)由于采用確定性的數(shù)列，故在較小樣本下就能達到穩(wěn)定。

2 擬隨機數(shù)列生成與轉(zhuǎn)換

2.1 Sobol數(shù)列生成

擬隨機數(shù)列有多種，本文采用Sobol[12]數(shù)列進行模擬，Sobol數(shù)列是基于叫做“直接數(shù)”的數(shù)[di]而構(gòu)造的。設(shè)[qi]是小于[2i]的正奇數(shù)，則[di=qi/2i]，數(shù)[di]（同時[qi]）的生成借助于系數(shù)[ai]只為0或1的簡單多項式。多項式表示成：

[f（y）=yp+a1yp-1+…+ap-1y+ap] （1）

對于[i>p]，有遞歸公式：

[di=a1di-1⊕a2di-2⊕…⊕apdi-p⊕[di-p/2p]] （2）

這里[⊕]表示二進制按位異或。對于[qi]，對等的遞歸公式為：

[qi=2a1qi-1⊕22a2qi-2⊕…⊕2papqi-p⊕qi-p] （3）

則可以得到Sobol數(shù)列第[i]個數(shù)：

[si=e1d1⊕e2d2⊕e3d3⊕…] （4）

這里[e1]，[e2]，…是[i]的二進制表示形式。

2.2正態(tài)隨機數(shù)列生成

得到了[0，1]區(qū)間上均勻分布的擬隨機Sobol數(shù)列，就可以用現(xiàn)有的一些算法將均勻分布轉(zhuǎn)化為所需要的其它分布。轉(zhuǎn)化的主要方法就是做累計分布函數(shù)的逆變換。一般在虛擬儀器系統(tǒng)中用到最多的就是均勻分布和正態(tài)分布，這里介紹轉(zhuǎn)化為正態(tài)分布的方法。累計標準正態(tài)分布是標準正態(tài)分布密度函數(shù)的積分。累計標準正態(tài)分布函數(shù)如下：

[Yx=12π-∞xe-t22dt] （5）

累積標準正態(tài)分布的Y軸是服從單位均勻分布U[0，1]的。則從累積標準正態(tài)分布的Y軸，我們可以找到相應(yīng)的其X軸上的樣本值，即[xY]，這就是我們需要的仿真值。通常使用的正態(tài)分布的逆變換方法是Box-Muller算法、Moro算法等，詳見參考文獻[13]，這里不再贅述。

3 基于擬蒙特卡羅方法的虛擬儀器不確定度評定

3.1不確定度源

虛擬儀器一般包括傳感器、信號調(diào)理轉(zhuǎn)換、信號采集、計算機及虛擬儀器軟件，系統(tǒng)框圖如圖3所示。

從理論來說，虛擬儀器各個組成部分均會產(chǎn)生不確定度源，但計算機計算引入的浮點舍入誤差并不大，在絕大部分情況下浮點舍入誤差的影響遠遠低于虛擬儀器硬件產(chǎn)生的誤差，故這里主要考慮傳感器與A/D轉(zhuǎn)換引起的不確定度源，其主要不確定源信息分別如圖4和圖5所示。

3.2 評定流程

擬蒙特卡羅方法與蒙特卡羅方法的評定方法類似，評定步驟主要為：1.獲取所有不確定度源信息，確定其概率分布信息。2.確定樣本容量M，產(chǎn)生一系列跟不確定度源分布相同的擬隨機數(shù)。3.若已知誤差傳遞函數(shù)，則根據(jù)以下公式進行誤差合成

[Δy=?f?x1Δx1+?f?x2Δx2+…+?f?xnΔxn] （6）

其中，[?f/?xi]是誤差傳播系數(shù)，[Δxi]是各不確定度源誤差分量，這樣可以得到M次試驗的誤差，統(tǒng)計其標準差，就可以得到虛擬儀器系統(tǒng)的標準不確定度

若未知誤差傳遞函數(shù)，虛擬儀器系統(tǒng)不確定度源大多屬于此類，可通過建立誤差仿真模塊，結(jié)合GUM評定方法，進行多次測量，對測量結(jié)果統(tǒng)計分析，得到虛擬儀器的合成不確定度，其評定流程如圖6所示。

3.3評定實例

某數(shù)字稱重傳感系統(tǒng)由稱重傳感器YZC-B1（激勵電壓5V，輸出電壓0～10mV），A/D轉(zhuǎn)換器AD7190和計算機及軟件組成。傳感器與A/D轉(zhuǎn)換器的不確定度源分別如表2和表3所示。

1）基于GUM方法B類不確定度評定

根據(jù)GUM方法B類評定規(guī)定，如不確定度源屬于均勻分布，則其標準不確定度計算公式為：[uxi=a/3]；如不確定度源屬于三角分布，則其標準不確定度計算公式為：[uxi=a/6]；如不確定度源屬于正態(tài)分布，則其標準不確定度計算公式為：[uxi=a/3]。這里假設(shè)環(huán)境溫度為25℃，分別計算其標準不確定度，各結(jié)果見表2和表3。例如由非線性引起的不確定度為：

[u3=0.0016×25/3=0.0231FS%] （7）

由噪聲誤差引起的不確定為：

[u8=0.01LSB/3=3.3333μV] （8）

故由傳感器引起的不確定為：

[us=10mV100u21+…+u24=7.7530μV] （9）

由A/D轉(zhuǎn)換引起的不確定為：

[uad=u25+…+u210=12.0226μV] （10）

則合成標準不準定度為：

[u=u2s+u2ad=14.3057μV] （11）

根據(jù)被測重量估計值與電壓估計值的關(guān)系為：

[y=fx=kx] （12）

計算出[k]值為[0.001g/μV]，則被測量的B類合成不確定為0.0143[g]。

2）依據(jù)蒙特卡羅方法及擬蒙特卡羅方法評定

按圖4擬蒙特卡羅方法虛擬儀器不確定度評定流程，用偽隨機數(shù)列和擬隨機sobol數(shù)列分別對上述不確定度源概率分布進行抽樣，抽取服從相應(yīng)分布的隨機數(shù)，樣本容量M分別取10，50，100，500，1000和5000，根據(jù)

[u=fx1+…+x10] （13）

統(tǒng)計其標準差，得到蒙特卡羅方法[mc]與擬蒙特卡羅方法[qmc]合成不確定度的結(jié)果如表4所示。

由上表可見，根據(jù)采樣次數(shù)的不同，蒙特卡羅方法評定的結(jié)果的范圍在0.0126至0.0171之間，擬蒙特卡羅方法評定的結(jié)果在0.0144至0.0147之間，在樣本容量較小時，如[M=10]和[M=50]，采用蒙特卡羅方法仿真得到的不確定度波動較大，而擬蒙特卡羅方法在小樣本仿真時，也可達到很好的效果。統(tǒng)計同一不確定度源不同仿真次數(shù)的標準差，采用蒙特卡羅方法仿真的標準差為1.4353，而采用擬蒙特卡羅方法仿真的標準差為0.1355，可看出擬蒙特卡羅方法的穩(wěn)定性要遠高于蒙特卡羅方法。如圖7所示，采用GUM方法、蒙特卡羅方法和擬蒙特卡羅方法評定該系統(tǒng)的不確定度結(jié)果比較，可看出，當采樣次數(shù)M>1000時，三者不確定度評定的結(jié)果有較好的一致性。

4 結(jié)束語

本文旨在探索虛擬儀器不確定度評定的新方法，對比偽隨機數(shù)列，分析了擬隨機數(shù)列的均勻性、收斂性、穩(wěn)定性以及擬隨機數(shù)列的生成與轉(zhuǎn)換分布。數(shù)字稱重實驗的表明，擬蒙特卡羅方法與傳統(tǒng)不確定度方法評定的結(jié)果具有較好的一致性，而且可以有效克服GUM方法與蒙特卡羅方法在評定時的缺陷與不足。

參考文獻：

[1] 喬仁曉，孟曉風.虛擬儀器測量不確定度研究[J].電子測量與儀器學(xué)報，2007，21（6）：52-55.

[2]王中宇，張海濱.測量不確定度最大殘差系數(shù)的一種新算法[J].計量學(xué)報，2006，27（3）：201-205.

[3]詹惠琴，習友寶.基于數(shù)值算法的虛擬儀器不確定度評估研究[J].儀器儀表學(xué)報，2005，26（10）：

[4] Salvatore Nuccio，Ciro Spataro.Assessment of Virtual Instruments Measurement Uncertainty [J].Computer Standards&Interfaces，2001，23（2001）：39-46.

[5] Ghiani E，Locci N，Muscas C.Auto-Evaluation of the Uncertainty in Virtual Instruments[J].IEEE Transactions on Instrumentation and Measureme- nt，2004，53（3）：672-677.

[6]Ghiani E，Locci N，Muscas C.Evaluation of Uncertainty in Measurements Based on Digitized Data[J].Measurement，2002，32（4）：265-272.

[7]王中宇，葛樂矣.一種小樣本虛擬儀器測量不確定度評定新方法[J].計量學(xué)報，2008，29（4）：387-392.

[8]袁敏，劉桂雄.基于Monte Carlo法虛擬儀器測量不確定度評定[J].科學(xué)技術(shù)與工程，2009，9（23）：7122-7125.

[9]王偉，宋明順.蒙特卡羅方法在復(fù)雜模型測量不確定度評定中的應(yīng)用[J].儀器儀表學(xué)報，2008，29（7）：1446-1449.

[10]雷桂媛.關(guān)于蒙特卡羅及擬蒙特卡羅方法的若干研究[D].杭州：浙江大學(xué)，2003：2-14.

[11]盧秀玉.蒙特卡羅方法和擬蒙特卡羅方法解線性方程組[J].東華大學(xué)學(xué)報，2010，36（2）：224-228.

[12]楊旭娟.亞式期權(quán)定價的擬蒙特卡羅模擬[D].武漢：武漢理工大學(xué)，2007.

[13]黃美發(fā)，景暉.基于擬蒙特卡羅方法的測量不確定度評定[J].儀器儀表學(xué)報，2009，30（1）：120-125.