俱巖飛， 王俊彪， 常崇義， 趙澤平

(1.中國(guó)鐵道科學(xué)研究院 研究生部，北京 100081； 2.中國(guó)鐵道科學(xué)研究院集團(tuán)有限公司 鐵道科學(xué)技術(shù)研究發(fā)展中心，北京 100081)

1 引 言

基于不確定度傳播率的GUM(guide to the expression of uncertainty in measurement)不確定度評(píng)定方法存在一定的局限性，如輸入量需呈對(duì)稱分布,輸出量需呈近似正態(tài)分布或t分布,測(cè)量模型為線性模型或可用線性模型近似的模型等[1,2]。而GUM的補(bǔ)充文件1[3]提出的基于分布傳播的蒙特卡洛方法(Monte Carlo method，MCM)打破了GUM法不確定度評(píng)定的局限，可以有效解決非線性測(cè)量模型、分布不對(duì)稱等情況下的不確定度評(píng)定問題。文獻(xiàn)[4]～[6]對(duì)MCM評(píng)定測(cè)量不確定度進(jìn)行了深入研究。文獻(xiàn)[7]～[11]研究表明：當(dāng)輸入量的概率分布為對(duì)稱分布、輸出量的概率分布近似為正態(tài)分布或t分布，測(cè)量模型為線性模型或可用線性模型近似表示的情況下，GUM法與MCM法得到的不確定結(jié)果基本一致；但當(dāng)輸入模型非線性或者輸入量為均勻分布、反正弦分布等，兩種方法法得到的輸出量分布形式有很大差異，且GUM法的評(píng)定結(jié)果較為保守；當(dāng)輸入量的分布區(qū)間不對(duì)稱時(shí)，會(huì)導(dǎo)致GUM法的評(píng)定結(jié)果的區(qū)間發(fā)生偏離?？偟膩碚f，MCM法比GUM法適用范圍更廣、可信度更高。

在實(shí)際的測(cè)量不確定度評(píng)定過程中，往往會(huì)遇到輸入量之間相關(guān)性不可忽略的情況。文獻(xiàn)[12]基于Cholesky分解產(chǎn)生多維正態(tài)分布抽樣序列；文獻(xiàn)[13]以最小二分法為基礎(chǔ)，給出了具有相關(guān)性的2個(gè)隨機(jī)變量的模擬方法；文獻(xiàn)[14]、[15]基于Cholesky因子線性-非線性變換產(chǎn)生了具有任意邊緣分布和相關(guān)系數(shù)的多維隨機(jī)抽樣序列，然而僅限于各輸入量之間相關(guān)系數(shù)矩陣正定即可以進(jìn)行Cholesky分解的情況。目前未見對(duì)于輸入量相關(guān)系數(shù)矩陣非正定的MCM評(píng)定測(cè)量不確定度的研究。

本文將提出一種基于Barzilai-Borwein梯度法[16]的迭代修正算法，對(duì)非正定的相關(guān)系數(shù)矩陣進(jìn)行修正，同時(shí)得到線性變換矩陣。結(jié)合Nataf變換[17]梳理出考慮相關(guān)性的MCM法評(píng)定測(cè)量不確定度的實(shí)施步驟，并對(duì)高速輪軌關(guān)系試驗(yàn)臺(tái)輪軌縱向蠕滑率測(cè)量不確定度進(jìn)行評(píng)定。

2 基于Nataf變換產(chǎn)生n維相關(guān)非正態(tài)隨機(jī)變量

對(duì)于相關(guān)系數(shù)矩陣為RX的n維非正態(tài)相關(guān)隨機(jī)變量X=[X1,X2,…,Xn]T，可由相關(guān)系數(shù)矩陣為RZ的標(biāo)準(zhǔn)正態(tài)向量Z=[Z1,Z2,…,Zn]T通過Nataf逆變換產(chǎn)生，其中RXij與RZij滿足如下函數(shù)關(guān)系：

φ(Zi,Zj,RZij) dZidZj

(1)

式中：μi、μj、σi、σj為Zi和Zj的期望和標(biāo)準(zhǔn)差；φ(Zi,Zj,RZij)為相關(guān)系數(shù)為RZij的二維標(biāo)準(zhǔn)正態(tài)分布聯(lián)合概率密度函數(shù)。

式(1)所示積分方程很難直接求解。文獻(xiàn)[17]通過大量的數(shù)值實(shí)驗(yàn)，得出了各種分布類型對(duì)應(yīng)的RXij和RZij之間的函數(shù)關(guān)系；文獻(xiàn)[14]采用了插值擬合方法進(jìn)行求解RZij；文獻(xiàn)[18]根據(jù)Weierstrass逼近定理，將RXij表示為RZij的多項(xiàng)式函數(shù)，求解得出RXij和RZij之間的關(guān)系，其中，當(dāng)Xi和Xj都服從均勻分布時(shí)，RXij和RZij滿足以下函數(shù)關(guān)系：

(2)

相關(guān)系數(shù)矩陣為RZ的標(biāo)準(zhǔn)正態(tài)向量Z=[Z1,Z2,…,Zn]T需要n維標(biāo)準(zhǔn)正態(tài)向量經(jīng)過線性變換產(chǎn)生。若RZ為正定矩陣，則所需的線性變換矩陣可通過對(duì)RZ進(jìn)行Cholesky分解產(chǎn)生；若RZ為非正定矩陣，則不能進(jìn)行Cholesky分解。

3 非正定相關(guān)系數(shù)矩陣迭代修正

由于計(jì)算誤差等原因，RZ可能出現(xiàn)接近“0”的負(fù)特征值。部分文獻(xiàn)提到用奇異值分解來處理非正定矩陣，事實(shí)上，對(duì)于含有負(fù)特征值的對(duì)稱矩陣，奇異值分解得到的左奇異矩陣和右奇異矩陣并不相等，即非正定矩陣不能寫成一個(gè)矩陣與其轉(zhuǎn)置乘積的形式。因此，非正定矩陣無法通過奇異值分解產(chǎn)生線性變換矩陣。本文通過迭代修正得到與原非正定矩陣盡可能接近的正定矩陣，進(jìn)而得出線性變換矩陣，并使得產(chǎn)生的相關(guān)多維隨機(jī)變量盡可能符合預(yù)期的相關(guān)性要求。

(3)

為了減少迭代步數(shù)，A的初始值應(yīng)該使AAT盡可能地接近RZ。對(duì)RZ進(jìn)行奇異值分解：

RZ=USVT

(4)

由于非正定矩陣奇異值分解產(chǎn)生的左右奇異矩陣U和V非常接近，因此可以設(shè)定A的初始值為:

(5)

以下為基于Barzilai-Borwein梯度法的矩陣迭代修正算法的具體實(shí)施步驟：

(1) 設(shè)k=1，按式(5)計(jì)算A的初使值A(chǔ)1，按式(3)計(jì)算目標(biāo)函數(shù)初始值f(A1),設(shè)初始步長(zhǎng)αk=0.5；

(2) 計(jì)算目標(biāo)函數(shù)的梯度：

(6)

(3) 確定步長(zhǎng)(k≥1)：

(7)

(4) 計(jì)算下一次的迭代值：

Ak+1=Ak-αkgk

(8)

(5) 按式(3)計(jì)算此步迭代的目標(biāo)函數(shù)值f(Ak+1)；

4 輸入量相關(guān)的MCM實(shí)施步驟

既有的MCM只針對(duì)輸入量相互獨(dú)立的情況，本節(jié)將考慮輸入量的相關(guān)性，簡(jiǎn)要梳理出含相關(guān)性的MCM實(shí)施步驟。首先通過對(duì)標(biāo)準(zhǔn)正態(tài)隨機(jī)變量的概率密度函數(shù)(PDF)離散抽樣，經(jīng)過Nataf逆變換產(chǎn)生符合輸入量概率分布及相關(guān)性的多維隨機(jī)變量；然后由測(cè)量模型傳播輸入量的概率分布，計(jì)算輸出量的PDF的離散抽樣值；最后由輸出量的離散分布數(shù)值獲取輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度和包含區(qū)間。

4.1 MCM輸入

4.1.1 確定測(cè)量模型及輸入量概率分布

首先確定測(cè)量模型：

Y=f(X1，X2，…，Xn)

(9)

然后根據(jù)有關(guān)專家的研究成果或經(jīng)驗(yàn)、貝葉斯定理、最大熵原理等[2～5]來確定各個(gè)輸入量的PDF。

4.1.2 確定輸入量之間的相關(guān)系數(shù)矩陣及線性變換矩陣

(1) 當(dāng)?shù)玫絻蓚€(gè)輸入量的估計(jì)值Xi和Xj，使用了同一個(gè)測(cè)量標(biāo)準(zhǔn)、測(cè)量?jī)x器、參考數(shù)據(jù)或采用了相同的具有相當(dāng)大不確定度的測(cè)量方法時(shí)，則Xi和Xj明顯相關(guān)[19]。

對(duì)Xi，Xj兩個(gè)輸入量同時(shí)觀測(cè)的n組測(cè)量數(shù)據(jù)，相關(guān)系數(shù)的估計(jì)值按下式計(jì)算：

(10)

(2) 當(dāng)某個(gè)輸入量由其它輸入量通過數(shù)學(xué)公式得出時(shí)，則這個(gè)輸入量和其它輸入量相關(guān)，可根據(jù)輸入量之間的數(shù)學(xué)關(guān)系推導(dǎo)出輸入量之間的相關(guān)系數(shù)矩陣。本文中的黏著輪半徑根據(jù)軌道輪半徑、軌道輪轉(zhuǎn)速和黏著輪轉(zhuǎn)速求出。

若相關(guān)系數(shù)矩陣RZ為正定矩陣，則對(duì)RZ進(jìn)行Cholesky分解產(chǎn)生線性變換矩陣；否則，根據(jù)本文第3節(jié)迭代修正算法，確定線性變換矩陣。

4.1.3 MCM試驗(yàn)次數(shù)

試驗(yàn)次數(shù)M至少應(yīng)大于1/(1-p)的104倍，通常取 105～106倍，也可采用自適應(yīng)蒙特卡洛方法[3]選擇M。

4.2 MCM傳播

(1) 從標(biāo)準(zhǔn)正態(tài)分布的PDF中抽取M個(gè)樣本值wi,r(i=1，2，…,n;r=1，2，…，M)。

(2) 基于Nataf逆變換將wi變換成相關(guān)系數(shù)矩陣為RX，邊緣分布為輸入量概率分布的隨機(jī)變量Xi，從而得到M個(gè)樣本值Xi,r(i=1，2，…，n；r=1，2，…，M)。

(3) 對(duì)每個(gè)樣本向量(X1,r,X2,r, …,Xn,r)計(jì)算相應(yīng)輸出量Yr的值

Yr=f(X1,r,X2,r,…,Xn,r) ,
r=1,2,…,M

(11)

4.3 MCM輸出

將M個(gè)輸出值嚴(yán)格按遞增次序排序，得到輸出量的PDF的離散表示，進(jìn)而獲取輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度和包含區(qū)間。

5 輪軌縱向蠕滑率測(cè)量不確定度評(píng)定

輪軌滾動(dòng)接觸蠕滑理論是研究輪軌接觸幾何關(guān)系的主要理論，而輪軌蠕滑率作為輪軌之間相對(duì)滑動(dòng)的一種度量參數(shù)，是輪軌接觸蠕滑理論的重要研究?jī)?nèi)容[20]。

高速輪軌關(guān)系試驗(yàn)臺(tái)作為試驗(yàn)速度能夠達(dá)到500 km/h的輪軌試驗(yàn)臺(tái)，可以模擬干燥、潮濕、撒砂等條件下的輪軌界面環(huán)境，可進(jìn)行輪軌黏著、蠕滑、脫軌、磨耗、疲勞、制動(dòng)、噪聲等試驗(yàn)。本節(jié)將基于上文所述方法對(duì)高速輪軌關(guān)系試驗(yàn)臺(tái)縱向蠕滑率試驗(yàn)結(jié)果的不確定度進(jìn)行評(píng)定。

為了較為全面地評(píng)定輪軌蠕滑率試驗(yàn)過程中的不確定度，綜合評(píng)價(jià)輪軌蠕滑率的試驗(yàn)結(jié)果，均勻選取了蠕滑試驗(yàn)中蠕滑率從0增加至20%再減小為0的整個(gè)試驗(yàn)過程中9組軌道輪和黏著輪轉(zhuǎn)速試驗(yàn)值，見表1。

表1 軌道輪和黏著輪轉(zhuǎn)速試驗(yàn)值Tab.1 Speed test value of track wheel and adhesive wheel r/min

5.1 輪軌縱向蠕滑率測(cè)量模型及各輸入量的概率分布

5.1.1 輪軌縱向蠕滑率測(cè)量模型

(12)

式中：nr和nw分別為軌道輪和黏著輪的轉(zhuǎn)速；Rr和Rw分別為軌道輪和黏著輪接觸點(diǎn)處的半徑。

5.1.2 設(shè)定輸入量的概率密度函數(shù)

(3) 軌道輪半徑的測(cè)量值為1 489.80 mm，測(cè)量精度為±0.01 mm，則將Rr設(shè)定為均勻分布U(1 489.79， 1 489.81)；

(4) 黏著輪半徑根據(jù)試驗(yàn)過程中縱向蠕滑率為0時(shí)，輪軌接觸點(diǎn)處軌道輪和黏著輪線速度相等來求出，即：

(13)

當(dāng)軌道輪和黏著輪速度為300 km/h時(shí)，設(shè)定縱向蠕滑率ζx=0，利用蒙特卡洛模擬，計(jì)算得出Rw服從參數(shù)為：下限a=439.65 mm，眾數(shù)c=439.81 mm，上限b=439.97 mm的三角形分布。

5.2 相關(guān)系數(shù)矩陣處理

根據(jù)式(13)，通過編程求得Rw、Rr、nw、nr之間的相關(guān)系數(shù)矩陣為：

由于RZ非正定，采用本文第3節(jié)中所述迭代修正算法對(duì)其進(jìn)行迭代修正，目標(biāo)函數(shù)值隨迭代次數(shù)增加的變化曲線如圖1所示。

圖1 目標(biāo)函數(shù)與迭代次數(shù)的關(guān)系曲線圖Fig.1 The relationship between the objective function f(A) and the number of iterations

迭代修正矩陣如下：

同時(shí)得到所需的線性變換矩陣為：

5.3 MCM傳播及結(jié)果報(bào)告

采用MCM對(duì)輪軌縱向蠕滑率不確定度進(jìn)行評(píng)定，設(shè)定抽樣次數(shù)為106，評(píng)定結(jié)果如表2所示。

表2 輪軌縱向蠕滑率不確定度評(píng)定結(jié)果Tab.2 Uncertainty evaluation results of wheel/rail longitudinal creep rate (%)

其中，縱向蠕滑率估計(jì)值為-0.029%的蒙特卡洛模擬結(jié)果如圖2所示。

圖2 蠕滑率估計(jì)值為-0.029%的MCM模擬結(jié)果概率密度直方圖Fig.2 Probability density histogram of MCM simulation results with an estimated creep rate of -0.029%

6 結(jié) 論

(1) 本文提出了一種基于Barzilai-Borwein梯度法，用于修正非正定相關(guān)系數(shù)矩陣的迭代修正算法。解決了基于Nataf逆變換產(chǎn)生服從任意邊緣概率分布的相關(guān)多維隨機(jī)變量過程中，相關(guān)系數(shù)矩陣非正定時(shí)，無法產(chǎn)生線性變換矩陣的問題。

(2) 描述了考慮相關(guān)性的MCM評(píng)定測(cè)量不確定度的具體實(shí)施步驟。

(3) 采用本文所述迭代修正算法并基于Nataf逆變換的MCM對(duì)高速輪軌試驗(yàn)臺(tái)輪軌縱向蠕滑率不確定度進(jìn)行了評(píng)定，結(jié)果表明了本文所述迭代修正算法有效且具有很快的收斂速度。