白 杰，胡紅波

(中國計(jì)量科學(xué)研究院，北京 100029)

1 引 言

自GUM實(shí)施以來，國內(nèi)外學(xué)者做了很多研究與討論[1～4]，包括從統(tǒng)計(jì)學(xué)不同角度(經(jīng)典統(tǒng)計(jì)與貝葉斯統(tǒng)計(jì))分析其特點(diǎn)與不足之處，這對后期GUM的修訂以及相關(guān)補(bǔ)充標(biāo)準(zhǔn)文件的發(fā)布實(shí)施起到了很大的作用[5～8]。

在實(shí)際的計(jì)量校準(zhǔn)工作中，除了直接對被測量進(jìn)行分析以外，回歸分析也是一種常用的數(shù)據(jù)分析方法。比如在振動加速度計(jì)量校準(zhǔn)中，要通過回歸分析從測量得到的數(shù)據(jù)中計(jì)算出表征運(yùn)動數(shù)學(xué)模型的參數(shù)，即幅度和相位等信息。有些儀器，其輸入輸出關(guān)系本身就是一個線性回歸模型，比如風(fēng)速傳感器輸出電壓與風(fēng)速的關(guān)系。對于簡單的線性回歸模型Y=θ0+θ1X，其參數(shù)至少包含回歸直線的截距θ0與斜率θ1參數(shù)，因其無法轉(zhuǎn)換一個確定的單一被測量模型[9]，無法直接依據(jù)GUM系列文件的規(guī)定進(jìn)行處理。本文依據(jù)GUM不確定度評估的準(zhǔn)則，針對計(jì)量中回歸模型參數(shù)計(jì)算的問題，在測量數(shù)據(jù)的基礎(chǔ)上采用最小二乘與貝葉斯推斷參數(shù)計(jì)算的方法計(jì)算了參數(shù)值以及相應(yīng)的不確定度。對于多參數(shù)計(jì)算的問題，通常需采用馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo，MCMC)法進(jìn)行數(shù)值計(jì)算[10]。本文也采用了MCMC方法，并且與最小二乘方法進(jìn)行了比較，說明2種方法的相同與不同之處。

2 測量模型與統(tǒng)計(jì)模型

2.1 GUM與測量模型

依據(jù)GUM系列文件進(jìn)行不確定度評估，首先需建立測量模型。測量模型表示被測量θ與各輸入量(Y,X)之間的函數(shù)關(guān)系，為了便于說明，假定測量模型如式(1)所示線性關(guān)系，實(shí)際測量模型可依據(jù)測量過程建立。

θ=φ(Y,X)=αY+βX

(1)

式中：θ為被測量；Y為通過測量得到的數(shù)據(jù)；X為通過其他方式得到的被測輸入量；α、β為固定常系數(shù)。通過測量得到一組數(shù)據(jù)Y后，依據(jù)GUM即可確定被測量θ的最佳估計(jì)值與不確定度。另外，若能依據(jù)測量過程以及其他信息確定測量方程(1)中所有輸入量的聯(lián)合分布fxy(x,y)，且能對其進(jìn)行抽樣，則能夠采用GUM S1中的方法確定被測量的近似概率密度函數(shù)，最終得到被測量的最佳估計(jì)值以及任意置信概率下區(qū)間。

2.2 貝葉斯推斷與統(tǒng)計(jì)模型

貝葉斯推斷主要是指利用貝葉斯定理計(jì)算參數(shù)的方法與過程，其出發(fā)點(diǎn)是數(shù)據(jù)，目的是為了得到參數(shù)的后驗(yàn)分布。貝葉斯定理形式如式(2)所示。

π(θ|y)∝L(y|θ)π(θ)

(2)

式中：π(θ|y)為參數(shù)的后驗(yàn)分布；L(y|θ)為似然函數(shù)，即測量數(shù)據(jù)的分布；π(θ)為參數(shù)的先驗(yàn)分布，即在測量數(shù)據(jù)得到之前對參數(shù)的認(rèn)識。

為了利用貝葉斯定理，首先需建立數(shù)據(jù)的統(tǒng)計(jì)模型，也稱為數(shù)據(jù)觀測方程[7]。對于式(1)所示的測量方程，相應(yīng)的觀測方程如式(3)所示。

(3)

式中ε為測量過程中的噪聲。若相互獨(dú)立的測量過程得到的數(shù)據(jù)y=(y1,y2,…,yN)，且噪聲ε為均值為0、方差為σ2的正態(tài)分布，即ε～N(0,σ2)。則式(2)所示的似然函數(shù)為如式(4)所示。

(4)

參數(shù)先驗(yàn)分布的確定也是利用貝葉斯定理的一個關(guān)鍵條件。在計(jì)量校準(zhǔn)的實(shí)際工作中，對某一個確定校準(zhǔn)對象，往往都會有較為明確的經(jīng)驗(yàn)信息，比如前面做過同樣型號的校準(zhǔn)，或者往年對同一對象校準(zhǔn)的數(shù)據(jù)等。如何由先驗(yàn)信息轉(zhuǎn)化為先驗(yàn)分布以及先驗(yàn)分布的各種取法可參考文獻(xiàn)[11]。對式(3)所示的測量方程而言，輸入量x與被測量θ可根據(jù)依據(jù)GUM S1推薦的確定輸入量分布的方法根據(jù)先驗(yàn)信息確定其合適的分布分別為π(x)與π(θ)。對被測量若無先驗(yàn)信息，則可取π(θ)∝C，C為常數(shù)。對測量過程中的噪聲方差，通常取其無信息先驗(yàn)形式，即π(σ2)∝(σ2)-1。由此可得被測量θ的后驗(yàn)分布如式(5)所示。

π(θ|y)∝?L(y|θ,σ2,x)π(σ2)π(x)π(θ) dxdσ2

(5)

由此也可以看出，應(yīng)用貝葉斯推斷進(jìn)行參數(shù)不確定度分析，不需要再進(jìn)行A類與B類方法的區(qū)分。

另外，若式(3)所示的觀測方程α、β未知，同樣可應(yīng)用貝葉斯推斷對其進(jìn)行估計(jì)。假定參數(shù)α、β的先驗(yàn)分布π(α,β)，則所有參數(shù)的聯(lián)合后驗(yàn)分布如式(6)所示。

π(α,β)π(x)π(θ)π(σ2) dx

(6)

對式(6)所示的聯(lián)合分布進(jìn)行積分，則可得到關(guān)于任何參數(shù)的后驗(yàn)分布。比如參數(shù)α、β的聯(lián)合后驗(yàn)分布為式(7)所示。

π(α,β|y)=?π(θ,α,β,σ2|y) dθdσ2

(7)

3 回歸模型參數(shù)的確定

對于式(1)或式(3)所示的模型，均可以通過變換將其等效為統(tǒng)一的關(guān)于參數(shù)的線性測量模型。假定某一計(jì)量校準(zhǔn)過程，在給定的一組設(shè)定試驗(yàn)條件X=(X1,X2,…,Xp)下，被校對象的輸出Y如式(8)所示。

Y=β0+β1X1+…+βpXp+ε

(8)

為簡化表示，將式(8)寫成矩陣形式，即Y=Xβ+ε。Y=(y1,y2,…,yN)T為N組不同設(shè)定試驗(yàn)條件下被校儀器的輸出；X=(1,xi1,xi2,…,xip)(i=1,2,…,N)為一個N×(p+1)維的矩陣；β=(β0,β1,…,βp)T為待確定的參數(shù)；ε為N維噪聲向量。

3.1 最小二乘算法

S(β)=argmin(Y-Xβ)T(Y-Xβ)

(9)

=σ2(XTX)-1

(10)

(11)

在得到參數(shù)估計(jì)值協(xié)方差矩陣后，其對角線上元素即為估計(jì)參數(shù)的方差，即

(12)

3.2 貝葉斯推斷參數(shù)估計(jì)

對于式(8)所示的回歸模型，為了便于分析參數(shù)計(jì)算過程，通常假定測量噪聲ε～N(0,σ2I)。在得到N組不同設(shè)定試驗(yàn)條件下被校儀器的輸出Y后，可得似然函數(shù)為：

L(Y|β,σ2)=(2πσ2)-N/2·

(13)

對于參數(shù)β，假定其先驗(yàn)分布為π(β)且互不相關(guān)，即π(β)=π(β0)π(β1)…π(βp)，噪聲方差的先驗(yàn)分布為π(σ2)，可得參數(shù)的后驗(yàn)分布為：

π(β,σ2|Y)∝L(Y|β,σ2)π(β)π(σ2)

(14)

有了所有參數(shù)的聯(lián)合后驗(yàn)分布，則可以通過積分得到感興趣參數(shù)的后驗(yàn)分布，從而計(jì)算參數(shù)的最佳估計(jì)值以及相應(yīng)置信概率下的區(qū)間。從式(14)可以看出，當(dāng)校準(zhǔn)對象的回歸模型參數(shù)較多(大于3個)時(shí)，通常難以通過積分得到單個參數(shù)概率分布的解析解，需采用MCMC算法來對后驗(yàn)分布進(jìn)行數(shù)值計(jì)算。

4 實(shí) 例

4.1 加速度計(jì)的振動校準(zhǔn)

加速度計(jì)，特別是壓電加速度計(jì)，在科學(xué)研究與工程實(shí)踐中廣泛應(yīng)用于機(jī)械結(jié)構(gòu)振動沖擊測量以及狀態(tài)監(jiān)測。頻率響應(yīng)和幅值線性度是加速度計(jì)2個關(guān)鍵的指標(biāo)。對于加速度計(jì)頻率響應(yīng)的校準(zhǔn)，一般采用振動校準(zhǔn)法。圖1是在設(shè)定頻率100 Hz，且振動加速度幅度大約50 m/s2下時(shí)，激光干涉儀測量得到的振動臺運(yùn)動的加速度波形數(shù)據(jù)。

圖1 振動臺運(yùn)動的加速度波形Fig.1 Acceleration waveform of vibrator

為了得到圖1所示序列的幅度和相位等信息，我們采用回歸的方法。由于振動臺按正弦規(guī)律運(yùn)動，故可設(shè)定其加速度波形為式(15)所示，即：

a(t)=A0cos(ω0t+θ0)+C0

(15)

式中：A0與θ0為關(guān)注且需確定的參量；ω0為設(shè)定的振動臺工作頻率；C0為未知的直流分量。對其按式(8)形式進(jìn)行展開并對參數(shù)做適當(dāng)變換可得[13]：

a(t)=A0cos(ω0t)cos(θ0)-A0sin(ω0t)sin(θ0)+C0

=A1cos(ω0t)-A2sin(ω0t)+C0

(16)

(17)

為了利用貝葉斯定理對參數(shù)進(jìn)行推斷，需設(shè)定參數(shù)的先驗(yàn)分布π(β,σ2)，假定各參數(shù)之間是相互獨(dú)立的，即π(β,σ2)=π(β)π(σ2)=π(A1)π(A2)π(C0)π(σ2)。直觀觀測圖1所示加速度波形圖再結(jié)合參數(shù)A1與A2的形式，設(shè)定π(A1)π(A2)～Unif[-60,60]，即為從-60到60之間的均勻分布，同時(shí)π(C0)～Unif[-10,10]，考慮到方差為尺度參數(shù)且必須大于0，選擇π(σ2)～χ2(3)，該范圍包括了實(shí)際振動加速度校準(zhǔn)過程中所有方差可能取到的值。利用MCMC數(shù)值計(jì)算的方法，對式(16)模型參數(shù)進(jìn)行貝葉斯推斷，得到主要參數(shù)的后驗(yàn)分布樣本、概率分布圖以及95%的區(qū)間(陰影部分),如圖2所示。

圖2 參數(shù)A1與A2的后驗(yàn)分布抽樣值與概率分布Fig.2 Samples and posterior distribution of parameter A1 and A2

從2種方法計(jì)算的結(jié)果可以看出，在參數(shù)先驗(yàn)分布取弱信息先驗(yàn)分布且數(shù)據(jù)量較多的情況下，2種方法計(jì)算的結(jié)果基本一致；同時(shí)從圖3也可以看出，2個參數(shù)基本相互獨(dú)立，這也與前面計(jì)算的結(jié)果一致。但當(dāng)數(shù)據(jù)量較少時(shí)，先驗(yàn)分布的選取對參數(shù)計(jì)算的結(jié)果有較大的影響[14]。

圖3 參數(shù)聯(lián)合分布Fig.3 Joint probability distribution for parameters

4.2 加速度計(jì)的沖擊校準(zhǔn)

加速度計(jì)的幅值線性度指標(biāo)，通常采用沖擊的方式校準(zhǔn)。將被校加速度計(jì)安裝在標(biāo)準(zhǔn)沖擊校準(zhǔn)臺上，測量輸入加速度計(jì)的沖擊加速度量值a(單位m/s2)，同時(shí)記錄加速度計(jì)的輸出電壓量y(單位mV)。設(shè)定加速度計(jì)輸入輸出是一個線性關(guān)系，即y=b+ka。由1組測量的數(shù)據(jù)即可確定直線截距b與系數(shù)k的值，再依據(jù)相關(guān)規(guī)程即可計(jì)算出被校加速度計(jì)在試驗(yàn)范圍內(nèi)的線性度指標(biāo)。表1為1典型的加速度計(jì)在沖擊加速度范圍1.0×102～5.0×104m/s2的試驗(yàn)結(jié)果。

表1 加速度計(jì)沖擊校準(zhǔn)結(jié)果Tab.1 Calibration results of acceleromter by shock

表1中的脈沖持續(xù)時(shí)間為沖擊校準(zhǔn)的一個試驗(yàn)條件，此處不是我們關(guān)注的對象。與振動校準(zhǔn)不同，沖擊校準(zhǔn)在不同的范圍有不同的測量不確定度，也就是測量噪聲模型為ε～N(0,σ2Ω)，Ω為N×N非奇異正定且對稱矩陣。故應(yīng)采用加權(quán)最小二乘算法進(jìn)行參數(shù)的計(jì)算。此時(shí)式(9)變?yōu)槿缡?18)所示[7,15]。

S(β)=argmin(Y-Xβ)TW(Y-Xβ)

(18)

(19)

利用貝葉斯推斷首先確定參數(shù)的先驗(yàn)分布。對于斜率，由同類加速度計(jì)校準(zhǔn)經(jīng)驗(yàn)數(shù)據(jù)可知該加速度計(jì)的靈敏度范圍在0.5到1.5之間，可設(shè)定π(k)～N(1,0.252)；對于截距，通常沒有1個明確的范圍，故可設(shè)定為無信息先驗(yàn)分布，即π(b)∝C，C為常數(shù)。利用MCMC算法，計(jì)算得到2個參數(shù)的后驗(yàn)分布以及后驗(yàn)聯(lián)合分布分別如圖4所示。

圖4 斜率與截距的后驗(yàn)分布Fig.4 Posterior distribution for slope and intercept

5 小 結(jié)

通過測量數(shù)據(jù)計(jì)算參數(shù)的值，即將參數(shù)表示為數(shù)據(jù)的函數(shù)，本文所述的最小二乘法與貝葉斯推斷法(包括參數(shù)的極大似然估計(jì)等)都是應(yīng)用較為廣泛的方法。其中貝葉斯推斷的最大的優(yōu)點(diǎn)在于可以利用計(jì)量校準(zhǔn)中的經(jīng)驗(yàn)信息或者是以前校準(zhǔn)的數(shù)據(jù)，并且以1種符合邏輯的方式并入測量數(shù)據(jù)的分析。本文對上述2種方法通過實(shí)際校準(zhǔn)中的例子進(jìn)行了比較，最小二乘法計(jì)算非常簡潔方便，可以迅速的得到參數(shù)的值；而貝葉斯推斷則計(jì)算較為復(fù)雜，通常需要采用數(shù)值計(jì)算的方式，但其能充分利用計(jì)量校準(zhǔn)中的經(jīng)驗(yàn)和歷史數(shù)據(jù)等信息。實(shí)際工作中可根據(jù)情況選擇合適的方法。