考虑非正定相关性的洛评定方法

俱岩飞， 王俊彪， 常崇义， 赵泽平

(1.中国铁道科学研究院 研究生部，北京 100081； 2.中国铁道科学研究院集团有限公司 铁道科学技术研究发展中心，北京 100081)

1 引 言

基于不确定度传播率的GUM(guide to the expression of uncertainty in measurement)不确定度评定方法存在一定的局限性，如输入量需呈对称分布,输出量需呈近似正态分布或t分布,测量模型为线性模型或可用线性模型近似的模型等[1,2]。而GUM的补充文件1[3]提出的基于分布传播的蒙特卡洛方法(Monte Carlo method，MCM)打破了GUM法不确定度评定的局限，可以有效解决非线性测量模型、分布不对称等情况下的不确定度评定问题。文献[4]～[6]对MCM评定测量不确定度进行了深入研究。文献[7]～[11]研究表明：当输入量的概率分布为对称分布、输出量的概率分布近似为正态分布或t分布，测量模型为线性模型或可用线性模型近似表示的情况下，GUM法与MCM法得到的不确定结果基本一致；但当输入模型非线性或者输入量为均匀分布、反正弦分布等，两种方法法得到的输出量分布形式有很大差异，且GUM法的评定结果较为保守；当输入量的分布区间不对称时，会导致GUM法的评定结果的区间发生偏离。总的来说，MCM法比GUM法适用范围更广、可信度更高。

在实际的测量不确定度评定过程中，往往会遇到输入量之间相关性不可忽略的情况。文献[12]基于Cholesky分解产生多维正态分布抽样序列；文献[13]以最小二分法为基础，给出了具有相关性的2个随机变量的模拟方法；文献[14]、[15]基于Cholesky因子线性-非线性变换产生了具有任意边缘分布和相关系数的多维随机抽样序列，然而仅限于各输入量之间相关系数矩阵正定即可以进行Cholesky分解的情况。目前未见对于输入量相关系数矩阵非正定的MCM评定测量不确定度的研究。

本文将提出一种基于Barzilai-Borwein梯度法[16]的迭代修正算法，对非正定的相关系数矩阵进行修正，同时得到线性变换矩阵。结合Nataf变换[17]梳理出考虑相关性的MCM法评定测量不确定度的实施步骤，并对高速轮轨关系试验台轮轨纵向蠕滑率测量不确定度进行评定。

2 基于Nataf变换产生n维相关非正态随机变量

对于相关系数矩阵为RX的n维非正态相关随机变量X=[X1,X2,…,Xn]T，可由相关系数矩阵为RZ的标准正态向量Z=[Z1,Z2,…,Zn]T通过Nataf逆变换产生，其中RXij与RZij满足如下函数关系：

φ(Zi,Zj,RZij) dZidZj

(1)

式中：μi、μj、σi、σj为Zi和Zj的期望和标准差；φ(Zi,Zj,RZij)为相关系数为RZij的二维标准正态分布联合概率密度函数。

式(1)所示积分方程很难直接求解。文献[17]通过大量的数值实验，得出了各种分布类型对应的RXij和RZij之间的函数关系；文献[14]采用了插值拟合方法进行求解RZij；文献[18]根据Weierstrass逼近定理，将RXij表示为RZij的多项式函数，求解得出RXij和RZij之间的关系，其中，当Xi和Xj都服从均匀分布时，RXij和RZij满足以下函数关系：

(2)

相关系数矩阵为RZ的标准正态向量Z=[Z1,Z2,…,Zn]T需要n维标准正态向量经过线性变换产生。若RZ为正定矩阵，则所需的线性变换矩阵可通过对RZ进行Cholesky分解产生；若RZ为非正定矩阵，则不能进行Cholesky分解。

3 非正定相关系数矩阵迭代修正

由于计算误差等原因，RZ可能出现接近“0”的负特征值。部分文献提到用奇异值分解来处理非正定矩阵，事实上，对于含有负特征值的对称矩阵，奇异值分解得到的左奇异矩阵和右奇异矩阵并不相等，即非正定矩阵不能写成一个矩阵与其转置乘积的形式。因此，非正定矩阵无法通过奇异值分解产生线性变换矩阵。本文通过迭代修正得到与原非正定矩阵尽可能接近的正定矩阵，进而得出线性变换矩阵，并使得产生的相关多维随机变量尽可能符合预期的相关性要求。

(3)

为了减少迭代步数，A的初始值应该使AAT尽可能地接近RZ。对RZ进行奇异值分解：

RZ=USVT

(4)

由于非正定矩阵奇异值分解产生的左右奇异矩阵U和V非常接近，因此可以设定A的初始值为:

(5)

以下为基于Barzilai-Borwein梯度法的矩阵迭代修正算法的具体实施步骤：

(1) 设k=1，按式(5)计算A的初使值A1，按式(3)计算目标函数初始值f(A1),设初始步长αk=0.5；

(2) 计算目标函数的梯度：

(6)

(3) 确定步长(k≥1)：

(7)

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

Ak+1=Ak-αkgk

(8)

(5) 按式(3)计算此步迭代的目标函数值f(Ak+1)；

4 输入量相关的MCM实施步骤

既有的MCM只针对输入量相互独立的情况，本节将考虑输入量的相关性，简要梳理出含相关性的MCM实施步骤。首先通过对标准正态随机变量的概率密度函数(PDF)离散抽样，经过Nataf逆变换产生符合输入量概率分布及相关性的多维随机变量；然后由测量模型传播输入量的概率分布，计算输出量的PDF的离散抽样值；最后由输出量的离散分布数值获取输出量的最佳估计值、标准不确定度和包含区间。

4.1 MCM输入

4.1.1 确定测量模型及输入量概率分布

首先确定测量模型：

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

(9)

然后根据有关专家的研究成果或经验、贝叶斯定理、最大熵原理等[2～5]来确定各个输入量的PDF。

4.1.2 确定输入量之间的相关系数矩阵及线性变换矩阵

(1) 当得到两个输入量的估计值Xi和Xj，使用了同一个测量标准、测量仪器、参考数据或采用了相同的具有相当大不确定度的测量方法时，则Xi和Xj明显相关[19]。

对Xi，Xj两个输入量同时观测的n组测量数据，相关系数的估计值按下式计算：

(10)

(2) 当某个输入量由其它输入量通过数学公式得出时，则这个输入量和其它输入量相关，可根据输入量之间的数学关系推导出输入量之间的相关系数矩阵。本文中的黏着轮半径根据轨道轮半径、轨道轮转速和黏着轮转速求出。

若相关系数矩阵RZ为正定矩阵，则对RZ进行Cholesky分解产生线性变换矩阵；否则，根据本文第3节迭代修正算法，确定线性变换矩阵。

4.1.3 MCM试验次数

试验次数M至少应大于1/(1-p)的104倍，通常取 105～106倍，也可采用自适应蒙特卡洛方法[3]选择M。

4.2 MCM传播

(1) 从标准正态分布的PDF中抽取M个样本值wi,r(i=1，2，…,n;r=1，2，…，M)。

(2) 基于Nataf逆变换将wi变换成相关系数矩阵为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的离散表示，进而获取输出量的最佳估计值、标准不确定度和包含区间。

5 轮轨纵向蠕滑率测量不确定度评定

轮轨滚动接触蠕滑理论是研究轮轨接触几何关系的主要理论，而轮轨蠕滑率作为轮轨之间相对滑动的一种度量参数，是轮轨接触蠕滑理论的重要研究内容[20]。

高速轮轨关系试验台作为试验速度能够达到500 km/h的轮轨试验台，可以模拟干燥、潮湿、撒砂等条件下的轮轨界面环境，可进行轮轨黏着、蠕滑、脱轨、磨耗、疲劳、制动、噪声等试验。本节将基于上文所述方法对高速轮轨关系试验台纵向蠕滑率试验结果的不确定度进行评定。

为了较为全面地评定轮轨蠕滑率试验过程中的不确定度，综合评价轮轨蠕滑率的试验结果，均匀选取了蠕滑试验中蠕滑率从0增加至20%再减小为0的整个试验过程中9组轨道轮和黏着轮转速试验值，见表1。

表1 轨道轮和黏着轮转速试验值Tab.1 Speed test value of track wheel and adhesive wheel r/min

5.1 轮轨纵向蠕滑率测量模型及各输入量的概率分布

5.1.1 轮轨纵向蠕滑率测量模型

(12)

式中：nr和nw分别为轨道轮和黏着轮的转速；Rr和Rw分别为轨道轮和黏着轮接触点处的半径。

5.1.2 设定输入量的概率密度函数

(3) 轨道轮半径的测量值为1 489.80 mm，测量精度为±0.01 mm，则将Rr设定为均匀分布U(1 489.79， 1 489.81)；

(4) 黏着轮半径根据试验过程中纵向蠕滑率为0时，轮轨接触点处轨道轮和黏着轮线速度相等来求出，即：

(13)

当轨道轮和黏着轮速度为300 km/h时，设定纵向蠕滑率ζx=0，利用蒙特卡洛模拟，计算得出Rw服从参数为：下限a=439.65 mm，众数c=439.81 mm，上限b=439.97 mm的三角形分布。

5.2 相关系数矩阵处理

根据式(13)，通过编程求得Rw、Rr、nw、nr之间的相关系数矩阵为：

由于RZ非正定，采用本文第3节中所述迭代修正算法对其进行迭代修正，目标函数值随迭代次数增加的变化曲线如图1所示。

图1 目标函数与迭代次数的关系曲线图Fig.1 The relationship between the objective function f(A) and the number of iterations

迭代修正矩阵如下：

同时得到所需的线性变换矩阵为：

5.3 MCM传播及结果报告

采用MCM对轮轨纵向蠕滑率不确定度进行评定，设定抽样次数为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梯度法，用于修正非正定相关系数矩阵的迭代修正算法。解决了基于Nataf逆变换产生服从任意边缘概率分布的相关多维随机变量过程中，相关系数矩阵非正定时，无法产生线性变换矩阵的问题。

(2) 描述了考虑相关性的MCM评定测量不确定度的具体实施步骤。

(3) 采用本文所述迭代修正算法并基于Nataf逆变换的MCM对高速轮轨试验台轮轨纵向蠕滑率不确定度进行了评定，结果表明了本文所述迭代修正算法有效且具有很快的收敛速度。