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

(1.中國鐵道科學研究院 研究生部，北京 100081； 2.中國鐵道科學研究院集團有限公司 鐵道科學技術研究發(fā)展中心，北京 100081)

1 引 言

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

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

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

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

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

φ(Zi,Zj,RZij) dZidZj

(1)

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

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

(2)

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

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

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

(3)

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

RZ=USVT

(4)

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

(5)

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

(1) 設k=1，按式(5)計算A的初使值A1，按式(3)計算目標函數(shù)初始值f(A1),設初始步長αk=0.5；

(2) 計算目標函數(shù)的梯度：

(6)

(3) 確定步長(k≥1)：

(7)

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

Ak+1=Ak-αkgk

(8)

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

4 輸入量相關的MCM實施步驟

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

4.1 MCM輸入

4.1.1 確定測量模型及輸入量概率分布

首先確定測量模型：

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

(9)

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

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

(1) 當?shù)玫絻蓚€輸入量的估計值Xi和Xj，使用了同一個測量標準、測量儀器、參考數(shù)據(jù)或采用了相同的具有相當大不確定度的測量方法時，則Xi和Xj明顯相關[19]。

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

(10)

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

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

4.1.3 MCM試驗次數(shù)

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

4.2 MCM傳播

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

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

(3) 對每個樣本向量(X1,r,X2,r, …,Xn,r)計算相應輸出量Yr的值

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

(11)

4.3 MCM輸出

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

5 輪軌縱向蠕滑率測量不確定度評定

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

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

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

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

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

5.1.1 輪軌縱向蠕滑率測量模型

(12)

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

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

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

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

(13)

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

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

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

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

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

迭代修正矩陣如下：

同時得到所需的線性變換矩陣為：

5.3 MCM傳播及結果報告

采用MCM對輪軌縱向蠕滑率不確定度進行評定，設定抽樣次數(shù)為106，評定結果如表2所示。

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

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

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

6 結 論

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

(2) 描述了考慮相關性的MCM評定測量不確定度的具體實施步驟。

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