計量中回歸模型參數(shù)值及其不確定度評估

白 杰，胡紅波

(中國計量科學研究院，北京 100029)

1 引 言

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

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

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

2.1 GUM與測量模型

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

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

(1)

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

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

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

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

(2)

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

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

(3)

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

(4)

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

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

(5)

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

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

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

(6)

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

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

(7)

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

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

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

(8)

為簡化表示，將式(8)寫成矩陣形式，即Y=Xβ+ε。Y=(y1,y2,…,yN)T為N組不同設(shè)定試驗條件下被校儀器的輸出；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ù)估計值協(xié)方差矩陣后，其對角線上元素即為估計參數(shù)的方差，即

(12)

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

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

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

(13)

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

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

(14)

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

4 實 例

4.1 加速度計的振動校準

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

圖1 振動臺運動的加速度波形Fig.1 Acceleration waveform of vibrator

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

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

(15)

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

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

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

(16)

(17)

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

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

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

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

4.2 加速度計的沖擊校準

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

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

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

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

(18)

(19)

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

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

5 小 結(jié)

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