基于擬蒙特卡洛方法的動態(tài)測量不確定度評定

漢澤西，邢靖虹(西安石油大學 電子工程學院，陜西 西安710065 )

0 引言

不確定度是對測量精度的定量表征[1]。當測量系統(tǒng)的特性隨時間發(fā)生變化時，測量系統(tǒng)不確定度將伴隨著工作時間的延續(xù)而發(fā)生變化與漂移。目前，國內(nèi)外對不確定度的表征與評定均依據(jù)1993年ISO公布的《測量不確定度指南》（簡稱GUM）。而我國也為了便于國際交流，于1999年發(fā)布了JJF059-1999《測量不確定度評定與表示》[2]。但二者均回避了關于動態(tài)不確定度的問題。隨著動態(tài)測量在整個測試領域里的比重日益增加，動態(tài)測量不確定度的評定與應用問題已成為現(xiàn)代不確定度理論的核心，關于動態(tài)不確定度的問題有待深入研究。

本文提出的動態(tài)測量不確定度評定方法，注重改進了隨機數(shù)的產(chǎn)生方法，在隨機過程理論以及應用蒙特卡洛方法評定測量不確定度的基礎之上，采用擬蒙特卡洛方法根據(jù)動態(tài)測量數(shù)據(jù)的特征，得到隨機變量的概率分布，產(chǎn)生空間分布更加均勻的擬隨機數(shù)，得到符合該隨機變量概率分布特性的隨機數(shù)值序列，作為動態(tài)測量系統(tǒng)的輸入進行仿真，并結合GUM闡述的不確定度評定方法，完成了動態(tài)測量的不確定度評定。該方法具備較高的計算精度，穩(wěn)定性和效率。

1 蒙特卡洛方法與擬蒙特卡洛方法

1.1 蒙特卡洛方法

蒙特卡洛（Monte Carlo，簡稱MC）方法，是一種統(tǒng)計試驗方法，它是采用基于“隨機數(shù)”的計算機隨機模擬方法。其基本思想是：當所求解問題是某種隨機事件出現(xiàn)的概率，或者是某個隨機變量的期望值時，通過某種“實驗”的方法，以這種事件出現(xiàn)的頻率估計這一隨機事件的概率，或者得到這個隨機變量的某些數(shù)字特征，并將其作為問題的解。但是，其生成算法具有確定性，這樣生成的序列在本質(zhì)上不是隨機的，通常稱之為偽隨機數(shù)。

偽隨機數(shù)的產(chǎn)生原理如下：首先基于一定方法，如線性同余法、乘同余法、逆同余法等，產(chǎn)生的服從均勻分布的隨機數(shù)，再根據(jù)相應的數(shù)學變換得到其他分布的隨機數(shù)。上述各方法均存在一定的不足，如高維不均勻性和長周期相關性現(xiàn)象，會導致仿真收斂速度慢及結果波動大等一系列問題。引起上述問題的一個主要原因是蒙特卡洛方法使用的偽隨機數(shù)隨機性過強而均勻性不足。利用蒙特卡洛方法解題主要步驟如圖1所示。

圖1 蒙特卡洛方法解題步驟

1.2 擬蒙特卡洛方法

傳統(tǒng)蒙特卡洛方法采樣往往會形成“空隙和簇”的現(xiàn)象，造成對采樣空間的搜索不充分。為了獲得分布更加均勻的采樣序列，人們提出了多種形式的確定性采樣方法，通常稱之為擬蒙特卡洛（Quasi-Monte-Carlo，簡稱QMC）方法。

擬蒙特卡洛方法是用精選的確定性的樣本點來代替蒙特卡洛采樣中的隨機性樣本點，其中確定性的樣本點是由比偽隨機數(shù)更加均勻地充滿采樣空間的低偏差序列通過某種變換而得的。擬蒙特卡洛方法計算的準確性及收斂速度主要取決于偏差。偏差用來度量點列在函數(shù)域上的均勻分布程度，點列分布越均勻，偏差就越小，收斂速度就加快，波動也就相應地減小，計算精度也隨之提高[3]。這里用Koksma-Hlawka不等式[4-5]給出擬蒙特卡洛方法的計算誤差：

式中：Q是空間 Id中包含的點集 x1, x2,… ,xN；對于任意屬于空間 Id的子集B， W*是Id上具有形式的所有子區(qū)間的集合；是d維勒貝格測度。

式中：CB是B的特征函數(shù)，F(xiàn)(B;Q)是計算滿足xi∈B(1≤i≤N)的那些點的個數(shù)的函數(shù)。

2 基于擬蒙特卡洛方法的動態(tài)測量不確定度評定

2.1 動態(tài)測量不確定度原理

為了解決傳統(tǒng)誤差評定方法的不足，人們引入不確定度的概念。結合1993年《測量不確定度指南》（GUM），將動態(tài)測量不確定度定義為：表征在相同的測量條件下（相對零時刻）合理地賦予被測量值的分散性，與測量結果相聯(lián)系的參數(shù)。動態(tài)測量不確定度是對測量結果的質(zhì)量的動態(tài)評價，具有時變性、隨機性、關聯(lián)性、動態(tài)性等特點。當實際測量系統(tǒng)的特性隨時間發(fā)生變化時，不但會使測量結果的估計值與被測量真值的差異逐漸增大，而且一系列隨機因素和未知因素對系統(tǒng)測量不確定性的影響增強，使測量標準差隨時間發(fā)生變化，即s=s(t)，則由測量標準差所決定的測量不確定度的量值隨時間延續(xù)不斷增大。從而，總不確定度也是一個隨時間變化的量，即測量不確定度 U= U(t)。圖2所示為被測量的估計值漂移及測量不確定度隨時間的變化趨勢[7]。

圖2 測量結果漂移和測量標準差變化圖

動態(tài)測量不確定度的評定通過不斷引入新數(shù)據(jù)，除去舊數(shù)據(jù)，跟蹤測量不確定度的變化規(guī)律，從而實現(xiàn)動態(tài)測量不確定度分析。

2.2 低偏差點集

低偏差點集的生成致力于點集的均勻性，偽隨機數(shù)序列側重于點集的隨機性。目前統(tǒng)計學家已經(jīng)提出了各種不同的低偏差點集，雖然構造點集的方法有差異，但這些點集序列都比偽隨機數(shù)序列在單位超立方體中分布更均勻。本文采用基于Halton序列的擬隨機數(shù)產(chǎn)生算法[8-9]：

若q為一質(zhì)數(shù)，則任意自然數(shù)k均有唯一的q進制表示：

對于 有 ai∈ {0 ,1,…, q ?1}，并且qr≤≤k≤≤qr+1。對任何可用下式表示的整數(shù)，且k≥1，令：

則q( k ) ∈ [ 0 , 1 ] 稱為k關于基q的根式逆運算。

令pi(1≤i≤m)為m個互不相同的質(zhì)數(shù)，則點集：

稱為Halton序列。

圖3和圖4分別表示偽隨機數(shù)序列與Halton序列抽樣點同為200時在單位面積上的分布。

圖3 單位面積上的偽隨機數(shù)序列

圖4 單位面積上的Halton序列

由圖3和圖4可以看出，Halton序列抽樣點的分布更加均勻，偽隨機數(shù)序列分布明顯有密集區(qū)和大的空白區(qū)，因此要填滿整個空間就需要抽取大量的樣本點，降低了計算效率。而低偏差序列在概率空間上分布較為規(guī)整，這使得擬蒙特卡洛方法比蒙特卡洛方法收斂速度更快且穩(wěn)定性更高。

2.3 評定步驟

本文提出的基于擬蒙特卡洛方法的動態(tài)測量不確定度評定步驟如下：

（1）建立動態(tài)測量過程的模型方程 ；

（2）確定m時刻動態(tài)測量序列不確定度分量的概率密度函數(shù) ；

（3）對各分量概率分布進行抽樣，確定抽樣樣本容量M。產(chǎn)生滿足指定概率密度函數(shù) 的擬隨機數(shù)；

（5）求得m時刻標準不確定度

重復執(zhí)行（2）～（5）步，得到系統(tǒng)動態(tài)測量不確定度。

由上述步驟可見，動態(tài)不確定度評定的QMC方法和MC法的根本區(qū)別就在于第（3）步中以各種不同低偏差點集取代偽隨機數(shù)序列，這樣不但繼承了MC方法程序容易實現(xiàn)的優(yōu)點，而且以其高精度和確定性的特點改進了傳統(tǒng)方法的缺陷，節(jié)省了大量計算成本。

3 應用實例分析

用相同測量條件下，30次慣性陀螺的1s采樣原始輸出脈沖數(shù)測量值為例，計算測量數(shù)據(jù)的動態(tài)不確定度，將蒙特卡洛方法與擬蒙特卡洛方法進行了分析對比。本次實驗暫不考慮環(huán)境、噪聲及系統(tǒng)內(nèi)部因素等對測量結果產(chǎn)生的影響，旨在突出隨機數(shù)產(chǎn)生方法改變后，擬蒙特卡洛方法與蒙特卡洛方法用于動態(tài)測量不確定度評定的不同效果。原始測量數(shù)據(jù)見表1。

表1 慣性陀螺的1s采樣原始脈沖測量數(shù)據(jù)（單位：脈沖/秒）

現(xiàn)取動態(tài)測量序列長度為20，即第一次計算取表1中序號為1~20的數(shù)據(jù)，第二次取表1中序號為2~21的數(shù)據(jù)，以此類推。表2為任意選取一組動態(tài)測量序列，進行不同樣本容量下的仿真計算結果和相應不確定度。表3給出傳統(tǒng)蒙特卡洛方法和本文方法計算得到的測量結果及相應不確定度。

表2 不同樣本容量的仿真結果

表3 不同方法得到的測量結果及標準不確定度

由表2中的數(shù)據(jù)可分析得出，當在不同樣本容量情況下進行仿真時，蒙特卡洛方法收斂速度比擬蒙特卡洛方法慢，且計算結果也不如擬蒙特卡洛方法穩(wěn)定。從表3中數(shù)據(jù)可以看出，本文方法評定的動態(tài)測量不確定度的精度優(yōu)于蒙特卡洛方法。

4 結束語

將數(shù)論網(wǎng)格中的低偏差點集引入到動態(tài)不確定度分析領域，結合測量不確定度評定指南《GUM》形成了動態(tài)測量不確定度評定的擬蒙特卡洛方法，并以慣性陀螺原始脈沖數(shù)為例，進行了數(shù)值算例分析。實驗結果表明，該方法的計算精度和效率優(yōu)于傳統(tǒng)的基于蒙特卡洛方法的不確定度評定方法，而且能夠提供確定的誤差估計，從而有效改善蒙特卡洛方法的缺陷。對于復雜動態(tài)系統(tǒng)的動態(tài)不確定度評定，擬蒙特卡洛方法能夠減少大量的有限元分析次數(shù)，從而大幅度提高計算效率，可得到優(yōu)于偽隨機數(shù)更好的結果，收斂速度快，結果穩(wěn)定，適用性強，為動態(tài)系統(tǒng)的不確定度評定提供了一條新思路。

[1] 高玉英,陳曉懷.應用蒙特卡羅方法計算動態(tài)測量的不確定度[J].黑龍江科技學院學報,2006,16(6).

[2] 漢澤西,肖志紅,董浩.現(xiàn)代測試技術[M].北京:機械工業(yè)出版社,2006:169-171.

[3] 戴鴻哲,王偉.結構可靠性分析的擬蒙特卡羅方法[J].航空學報,2009,30(4):2-5.

[4] 張艷林,朱麗莎,張義民,等.機械強度可靠性靈敏度分析的擬蒙特卡羅法[J].東北大學學報,2010,31(11):2-4.

[5] Ronald Kleiss,Achilleas Lazopoulos.Error in Monte Carlo,quasi-error in Quasi-Monte Carlo[J].Computer Physics Communications,2006,175:93-115.

[6] 孟軻.擬蒙特卡羅方法在Variance-Gamma模型里的應用[D].北京:清華大學,2009:11-12.

[7] 謝少鋒,陳曉懷,張勇斌.測量系統(tǒng)不確定度分析及其動態(tài)性研究[J].計量學報,2002,23(3):1-2.

[8] B IPM/JCGM.Evaluation of measurement data-supp lement 1 to the‘guide to the expression of uncertainty in measurement’- propagation of distributions using a monte carlo method[J].D raft Supp lement 1,2006:4-12.

[9] Xavier Legrand.A Quasi-Monte Carlo Integration Method Applied to the Computation of the Pollaczek Integral[J].IEEE Transactions On Power Delivery,2008,23(3):3-8.