許金鑫, 由 強(qiáng)

(1. 中國計量科學(xué)研究院, 北京 100029; 2. 清華大學(xué) 電機(jī)工程與應(yīng)用電子技術(shù)系, 北京 100084)

1 引 言

多項式最小二乘擬合是數(shù)據(jù)處理中常用的曲線擬合方法。測量不確定度評定的指導(dǎo)性文件(guide to the expression of uncertainty in measurement, GUM)給出了一階線性最小二乘擬合的不確定度評估方法[1],然而,對于更高階的多項式最小二乘擬合,卻很難從GUM給出的方法中推出其不確定度的理論計算方法。一些其他的評估直線或者平面擬合不確定度的方法,并不適用于最普遍的情況[2～5]。Hibbert在2006年提出了一種相對完備的方法,給出了線性、拋物線和加權(quán)線性擬合的不確定表達(dá)式[6],然而,三階以上的多項式擬合仍然很難從這種方法中得到不確定度的表達(dá)式。蒙特卡羅法可以用來評估任意階次的最小二乘擬合的不確定度[7～10]。但這種方法只適用于真實(shí)的測量函數(shù)已知的情況,如算法仿真測試[8,9]和測量超精密自由曲面[7]等。因此,采用一種通用的方法來給出任意階次多項式擬合不確定度的表示式是很有必要的。

本文提出了一種采用矩陣方程表示的任意階次多項式最小二乘擬合的不確定度計算方法。根據(jù)這個方法,可以進(jìn)一步得到測量數(shù)據(jù)的最佳擬合階次以及擬合后積分計算的不確定度。同時,本文通過仿真進(jìn)一步研究了該方法的特性,驗(yàn)證了其查找最佳擬合階次的有效性。

2 原理分析

2.1 擬合系數(shù)的方差

m組測量數(shù)據(jù)的n階多項式最小二乘擬合可以看成是求解n+1個參數(shù)的線性方程組,方程組個數(shù)為m個,一般情況下測量數(shù)據(jù)個數(shù)m大于n+1,所以,這個方程組為超定方程組,可以寫為：

(1)

(2)

這個超定方程組可以由廣義逆矩陣求解,寫成：

(3)

式中Pinv(A)或者K都為矩陣A的廣義逆矩陣。由于測量存在誤差,Ys并不是對應(yīng)測量點(diǎn)的真實(shí)結(jié)果,假設(shè)對應(yīng)測量點(diǎn)的真實(shí)結(jié)果為Y,可以認(rèn)為測量結(jié)果滿足的真實(shí)的多項式關(guān)系參數(shù)為α0,…,αn。與式(3)類似,可以得到真實(shí)參數(shù)表達(dá)式：

α=Pinv(A)Y=KY

(4)

測量結(jié)果與真實(shí)結(jié)果之間的誤差可以表示為：

εi=ysi-yi

(5)

式中ysi和yi分別為Ys和Y中的第i個元素。將式(3)展開并代入式(5),消掉Ys中的元素,得到：

(6)

(7)

對式(7)兩邊同時計算方差,得到：

(8)

式(8)中,由于αi為真實(shí)的多項式關(guān)系參數(shù),為固定值,其方差為0。同時,假設(shè)每次測量的誤差是獨(dú)立的且同分布的,式(8)可以化簡為：

(9)

式中σ為測量值和真實(shí)值之間的標(biāo)準(zhǔn)偏差。然而,真實(shí)值在測量過程中是無法得到的,所以,這里采用擬合殘差的標(biāo)準(zhǔn)偏差s來代替σ。擬合殘差為測量結(jié)果與擬合后得到結(jié)果的差,對于m個測量數(shù)據(jù)點(diǎn),采用n階擬合時,得到的m個擬合殘差數(shù)據(jù)的自由度為m-n-1,所以可以得到s的表達(dá)式為：

(10)

式中Δyi為擬合殘差。用s來代替σ后,式(9)進(jìn)一步化簡為：

(11)

2.2 擬合系數(shù)的協(xié)方差矩陣

(12)

由于每次測量的誤差是獨(dú)立的且同分布的,可以得到誤差之間的協(xié)方差為：

(13)

將式(13)代入式(12)中并采用擬合殘差的標(biāo)準(zhǔn)偏差s來代替σ,式(12)可以簡化為：

(14)

結(jié)合式(11)和式(14),可以得到擬合參數(shù)的協(xié)方差矩陣C可以表示為：

C=s2KTK

(15)

這里注意在計算式(15)時,需要對測量點(diǎn)數(shù)據(jù)進(jìn)行歸一化：

(16)

式中：X為由測量點(diǎn)x1，…，xm構(gòu)成的矢量；mean(X),max(X)和min(X)分別為x1，…，xm的平均值、最大值和最小值。歸一化后,可以使得矩陣A中的元素相差不大,從而使得在計算廣義逆矩陣K時避免截斷誤差。

2.3 擬合結(jié)果的不確定度

在擬合區(qū)間內(nèi)的任意點(diǎn)x0所對應(yīng)的擬合結(jié)果y0可表示為：

(17)

根據(jù)不確定度的合成法則,y0的合成不確定度可以表示為：

u2(n,y0)=(xn0)2Var(αn)+…+(x0)2Var(α1)+

xn0xn-10Cov(αn,αn-1)+…+

x20x0Cov(α2,α1)

(18)

將式(18)改為為矩陣表達(dá)：

u2(n,y0)=X0CXT0

(19)

式中矩陣C為擬合參數(shù)的協(xié)方差矩陣,可由式(15)求得。矩陣X0為：

X0=[xn0xn-10…x01]

(20)

將式(15)代入式(19)中并開根號,就可以得到擬合結(jié)果y0的不確定度的表達(dá)式：

(21)

由式(21)可知,y0的不確定度u(y0)由3個變量決定：一是擬合殘差的標(biāo)準(zhǔn)偏差s,擬合殘差越小,擬合結(jié)果的不確定度越小；二是矩陣K,矩陣K由擬合階次和測量點(diǎn)決定,與測量結(jié)果無關(guān),不同的測量點(diǎn)和擬合階次決定了不同的擬合不確定度；三是矩陣X0,為擬合點(diǎn)的位置,也就是說不同的擬合點(diǎn)對應(yīng)的擬合不確定度不同。

2.4 最佳擬合階次

由第2.3節(jié)可知,同一條擬合曲線上,每一個擬合點(diǎn)所對應(yīng)的擬合不確定度都不相同。為了得到最佳的擬合階次,就需要定義一個擬合曲線的平均擬合不確定度。在擬合區(qū)間內(nèi),等間隔的取一定數(shù)量的采樣點(diǎn),可以得到每個采樣點(diǎn)所對應(yīng)的擬合不確定度u(y0i)。將每個采樣點(diǎn)看成是獨(dú)立的一次測量,則可以得到采樣點(diǎn)平均的擬合不確定度為：

(22)

式中l(wèi)為采樣點(diǎn)的個數(shù)。當(dāng)采樣間隔足夠小時,式(22)則可以寫為積分的表達(dá)式：

(23)

式中[a,b]是擬合的區(qū)間。對于不同擬合階次所得到的擬合曲線,都可以得到一個平均的擬合不確定度。而平均擬合不確定度最小的那個曲線對應(yīng)的擬合階次,認(rèn)為是最佳的擬合階次。

2.5 擬合曲線積分結(jié)果的不確定度

在一些場合中,需要對擬合后的曲線做積分,而積分結(jié)果的不確定度同樣與擬合有關(guān)。n階多項式擬合后的積分結(jié)果可以表示為：

(24)

與從式(17)到式(21)的推導(dǎo)相類似,根據(jù)不確定度的合成法則,式(23)積分結(jié)果的不確定度可以表示為：

u2int(n,x1,x2)=X12NCNTXT12

(25)

同樣,代入將式(15)代入式(25)中并開根號可以得到：

(26)

其中矩陣N和X12分別為：

(27)

(28)

3 實(shí)例仿真分析

3.1 仿真模型

假設(shè)一個實(shí)際的物理系統(tǒng)所滿足的真實(shí)的多項式關(guān)系如下：

(29)

式(29)描述的原函數(shù)在實(shí)際測量過程中是無法得知的,這里采用這樣的假設(shè)模型來對本文中描述的方法進(jìn)行仿真分析。在式(29)定義的區(qū)間等間隔取樣x1,…,xm各點(diǎn),由式(29)可直接得到取樣點(diǎn)所對應(yīng)的真實(shí)值y1,…,ym。

在真實(shí)值上疊加一些隨機(jī)的誤差,模擬實(shí)際測量時的誤差。隨機(jī)的誤差滿足矩形分布,其幅值相對于測量值約為10-2量級,表示為：

(30)

式(29)中e也可以看成是測量值誤差的最大值。增加誤差后,測量得到的值可以表示為：

ysi=yi+rand(2e)-e

(31)

其中rand(2e)表示從0到2e之間產(chǎn)生矩形分布的隨機(jī)數(shù)。所以,由式(31)產(chǎn)生的測量值ysi是從yi-e到y(tǒng)si+e之間分布的,其標(biāo)準(zhǔn)測量不確定度為：

(32)

3.2 擬合不確定度計算

由式(21)可得到每個點(diǎn)的擬合不確定度。將式(21)看成是2部分的乘積,一部分是s,可由式(10)計算；另一部分寫成t,表示為：

(33)

s只與測量的結(jié)果和擬合的結(jié)果有關(guān),所以這里將其稱為y軸相關(guān)量；t只與采樣的點(diǎn)和計算點(diǎn)位置有關(guān),這里將其稱為x軸相關(guān)量。y軸相關(guān)量反應(yīng)了測量結(jié)果和擬合結(jié)果之間的偏差對不確定度的影響。x軸相關(guān)量反應(yīng)了不同擬合階次和擬合點(diǎn)位置對不確定度的影響。下面分別對這2部分進(jìn)行計算和分析。

在仿真計算中,一共模擬產(chǎn)生了21個測量數(shù)據(jù),分別進(jìn)行了從1階到18階擬合。其中,y軸相關(guān)量的計算結(jié)果如圖1所示,總體上,隨擬合階次升高,y軸相關(guān)量減小。圖1中的放大示圖可以看出,從10階到12階擬合,以及從13階到16階擬合出現(xiàn)的小幅度的上升,由式(10)可知,這是由于自由度的減小造成的。當(dāng)采用13階擬合時,y軸相關(guān)量的值最小,與實(shí)際的原函數(shù)為5階相差較大。

圖1 仿真得到的y軸相關(guān)量與擬合階次的關(guān)系Fig.1 The relationship between y-related part of fitting uncertainty and fitting order in this simulation

x軸相關(guān)量t的計算結(jié)果如圖2所示,對圖2中得到的不同擬合階次下的t求平均值就可以得到如圖3所示的曲線。

圖2 不同擬合階次下的擬合曲線Fig.2 The fitted value with respect to different fitting order

圖3 x軸相關(guān)量tave隨擬合階次的變化Fig.3 The average of x-related part tave with respect to the fitting order

(34)

式(34)實(shí)際上就是式(22)除以s得到的結(jié)果。由圖2和圖3可以得到如下3個結(jié)論：1) 在擬合曲線的中間部分,x軸相關(guān)量較小,而在擬合曲線的兩頭,x軸相關(guān)量變大；2) 隨著擬合階次的增加,tave單調(diào)增加。在11階以下時,tave基本為線性增加,當(dāng)超過12階時,tave急速增加,這也是為什么一般不選取高階多項式擬合的原因。3) 隨著擬合階次的增加,x軸相關(guān)量在擬合區(qū)間兩邊增加的速度遠(yuǎn)大于中間區(qū)域。從圖2中可以看出,5階到15階的擬合曲線都與實(shí)際測量點(diǎn)符合的很好,但是到了18階的時候,在最右端出現(xiàn)了比較大的波動,這就是所謂的龍格現(xiàn)象[11],在該位置,x軸相關(guān)量也相應(yīng)的出現(xiàn)了極大值。

s與tave的乘積就可以得到由式(23)所定義的曲線的平均擬合不確定度,再除以測量數(shù)據(jù)的平均值,就可以得到相對的擬合不確定度。圖4顯示了擬合曲線的相對平均不確定度隨擬合階次的變化。在低階的擬合階段,相對平均不確定隨擬合階次的增加快速減小,主要是因?yàn)閟隨擬合階次1階到5階快速從5.1×10-1快速減小到了2.3×10-3,而tave再線性地緩慢增加。擬合階次從5階到15階變化時,s減小的速度變緩,tave增加的速度超過了s減小的速度,使得整體的曲線隨擬合階次緩慢增加,從而形成了一個極小值點(diǎn)。當(dāng)擬合階次進(jìn)一步增加時,tave增加的速度大大提升,從而使得整體的擬合不確定度快速增加。所以,在這個仿真模型中,通過本文的計算方法,最終得到的最佳擬合階次為5階,與假設(shè)的原函數(shù)階次相同,擬合的平均不確定度為2.3×10-3。同時也可以看出,擬合的平均不確定度是小于測量的不確定度的,所以說通過最小二乘擬合,是可以一定程度上減小測量引入的不確定度的。

圖4 擬合相對不確定度隨擬合階次的變化Fig.4 Relative uncertainty with respect to the fitting order

4 結(jié) 論

在多項式最小二乘擬合中,每一點(diǎn)的擬合不確定度都是不同的。中間點(diǎn)的擬合不確定度要遠(yuǎn)小于兩邊的擬合不確定度。同時,測量的不確定和測量點(diǎn)會影響到曲線的擬合不確定度。曲線的擬合不確定度可以看成是x軸相關(guān)量t和y軸相關(guān)量s的乘積。當(dāng)擬合階次增加時,由于擬合殘差降低,s減小,然而t會增加,所以存在不確定度最小的擬合階次,為最佳擬合階次。