含相关性的测量不确定度拟蒙特卡罗评定方法-论文.pdf

'第35卷第6期仪器仪表学报Vo1．35No．62014年6月ChineseJournalofScientificInstrumentJun．2014含相关性的测量不确定度拟蒙特卡罗评定方法凌明祥，李会敏，黎启胜，李明海(中国工程物理研究院总体工程研究所绵阳621900)摘要：针对测量不确定度评定中存在相关性分量以及蒙特卡罗法收敛速度和精度较低的问题，提出能处理不确定度源相关且兼顾均匀性和随机性的复杂系统测量不确定度拟蒙特卡罗评定方法。该方法通过产生均匀分布的独立随机化Sobol序列，并基于Cholesky因子线性一非线性变换产生具有任意边缘分布和相关系数的多维随机抽样序列，建立测量模型进行概率分布传播得到被测埴的估计值和包含区间。算例和实际测量表明，该方法具有比蒙特卡罗法更高的精度和效率，而且除保持一般拟蒙特卡罗法的优良低偏差性外更能真实反映测量不确定度传播的随机过程，尤其适用于测量不确定度源服从非正态分布且相关的情况。对采用拟蒙特卡罗法评定测量不确定度的一些问题进行了讨论和总结。关键词：测量不确定度；蒙特卡罗法；相关性；随机变量抽样；Sobol序列中图分类号：TB9文献标识码：A国家标准学科分类代码：460．40QuasiMonteCarlomethodforthemeasurementuncertaintyevaluationconsideringcorrelationLingMingxiang，LiHuimin，LiQisheng，LiMinghai(InstituteofSystemEngineering，ChinaAcademyofEngineeringPhysics，Mianyang621900，China)Abstract：Aimingattheproblemsthatthereexistcorrelationcomponentsandtheconvergencespeed&precisionofMonteCarlomethod(MCM)arelowinmeasurementuncertaintyevaluation，aquasiMonteCarlomethod(QMCM)isproposedforthemeas—urementuncertaintyevaluationofcomplicatedsystem，whichcandealwithcorrelateduncertaintysourcesandconsiderbothuni-formityandrandomness．ThroughgeneratingindependentrandomizedSobolsequencewithuniformdistribution，thismethodgener—．atesmulti—dimensionalrandomvariablesamplingsequencewitharbitrarymarginalprobabilitydistributionandtargetcorrelationCO-efficientbasedonimprovedlinear—nonlineartransformingtechniquewithCholeskyfactor．Themeasurementmodelisbuilt，theprobabilitydistributionofthemeasurandispropagatedandcalculated，andtheestimationandcoverageintervalofthemeasurandareobtainedfromitssamplingsequence．Theexampleandactualmeasurementillustratethattheproposedmethodhashigherpreci·sionandefficiencycomparedwithMCM；besideskeepingtheexcellentlowdeviationperformanceofQMCM，thismethodcanreallyreflecttherandomprocessofmeasurementuncertaintypropagation．Theproposedmethodisespeciallysuitableforthecaseswheretheuncertaintysourcesarecorrelatedandhavenon—normaldistribution．SomeissuesaboutusingQMCMformeasurementuncer—taintyevaluationarealsodiscussedandsummarized．Keywords：measurementuncertainty；MonteCarlomethod(MCM)；correlation；randomvariablesampling；Sobolsequence负参数，对计量、检测以及比对实验具有重要意义。JJF(测量不确定度评定与表示》⋯以及《Guidetothe1引言expressionofuncertaintyinmeasurement)(GUM)提供了一套基于不确定度传播律的测量不确定度评定测量不确定度是表征赋予被测量值分散性的非方法。然而对于测量模型明显呈非线性、不易求灵收稿日期：2014-01ReceivedDate：2014-01基金项目：国家重大科学仪器设备开发专项(2011YQ13O047)资助项目 1386仪器仪表学报第35卷敏度系数、不确定度分量相差较大以及被测量(输关系确定时，输出与输入的估计值可表示为：出量)概率密度分布较大程度偏离正态分布或t分Y=／1，2，⋯，)(1)布的情况，GUM方法的精度和可操作性差。对被测量l，的系统误差进行修正以及对随机误差补针对复杂测量，基于蒙特卡罗法的测量不确定偿，得到l，的测量模型为：度评定具有独特优势，适用范围更广，该方面的研究Y=f(x。，⋯，Xn)+∑A+∑(2)已有较多文献介绍，而且较为充分深入。然而即：这些研究均没有考虑不确定度分量之间的相关性，Y=X)，其中：=(一，，A一，△，占一，不能适用于不可忽略相关性的场合。针对以上不6)足，国内外人员进行了部分研究，但都是基于式中：△和6分别是被测量l，的系统误差修正值和随机Cholesky分解变换得到相关的多维正态分布抽样序误差补偿值。列后进行被测量的概率分布传播，目前未见能处理测量所得的值只是对被测量的估计，测量过程中不确定度分量为非正态分布且相关的方法报道，而的系统效应和随机效应均会导致测量不确定度，式具有非正态分布且相关的不确定度源在实际测量中(2)的测量模型考虑了被测量l，的由系统效应和随机更加常见，例如由同一仪器的示值误差或相同基准效应引入的不确定度源。其中，△和6的值可以而引入的不确定度分量一般视为相关的均匀分布¨。且经常为零，但必须考虑系统误差修正和随机效应引蒙特卡罗计算精度方面，文献[16]采用均匀性更好入的不确定度分量。需注意的是测量模型中各分量的Hahon序列代替伪随机数进行被测量的概率分布的单位应统一，不确定度来源尽量考虑全面且不传播，进行不确定度评定，弥补了蒙特卡罗法使用伪重叠。随机数均匀性不足而导致计算结果不够稳定以及收敛较慢的不足。然而，6维以上Hahon序列的高度2．2概率分布传播相关性会对多不确定度源的被测量概率分布传播造若确定了测量模型中各输入量(包括各修正和成较大误差。采用Hahon、Faure、Sobol序列等随机补偿量)、△和6的概率分布及相关性，即可确定数的拟蒙特卡罗法虽然在多重积分计算等领域获得各输入量的随机离散抽样序列，进而南测量模型式巨大成功，但若直接将其用于测量不确定度评定，不(2)传播输入量的分布，计算输出量y的概率密度能真实反映被测量的概率分布传播的随机性，尤其函数离散抽样值，得到】，的最佳估计值、标准不确定是采用多维拟随机序列来传播具有多个不确定度源度以及包含区间。整个计算流程如图1所示。的被测量概率分布时容易失误。本文将提供一套考虑不确定度源具有任意边缘概率分布和相关系数的测量不确定度拟蒙特卡罗评定方法。该方法的核心在于拟随机数的产生和随机化以及如何生成具有任意边缘概率分布和相关系数的多维不确定度源随机分量抽样序列。2评定方法框架及流程蒙特卡罗是一种数值计算方法，其解决问题的思路是先建立所求解的概率模型，使所求问题的解正好是所建概率模型的数学期望或其他特征量，然后多次模拟，统计出某事件发生的概率，求出要估计的参数。拟蒙特卡罗法与蒙特卡罗法的区别在于使用的随机数不同，前者采用低偏差的拟随机数，均匀性更好却牺牲了随机性。图1(拟)蒙特卡罗测量不确定度评定流程基于拟蒙特卡罗法的测量不确定度评定包括测量模型建Fig．1Uncertaintyevaluationprocedurewith(quasi)MCM立和被测量概率分布传播2个阶段。2．1测量模型建立1)拟蒙特卡罗输入被测量y由l"1个输入量。、，⋯，通过函数确定测量模型输入量个数Ⅳ及各输入分量服从 第6期凌明祥等：含相关性的测量不确定度拟蒙特卡罗评定方法1387的概率分布，并明确其中的相关性。首先产生Ⅳ维的性质：即任意长的子序列都能均匀地填充整个函相互独立的[0，1]区间均匀分布拟随机序列，每一数空间。鉴于Sobol序列高维(260维以上)才出现维序列的长度为拟蒙特卡罗试验次数。并对该Ⅳ退化的优势，本文选择Sobol序列作为所有测量不确维序列进行随机化使其成为均匀(低偏差)性好的定度源概率分布抽样序列的种子序列，该序列为独立随机抽样序列。不能直接将拟随机序列用于被(0，1)区间的均匀分布。产生Sobol序列的过程测量的概率分布计算，原因是拟随机序列自身是相如下：关的确定性点列，不能反映测量不确定度传播的随取Ⅳ个系数为0或1的本原多项式，各多项式的阶机性，造成输出量的估计不符合实际情况。数不同。本原多项式为：以上述方法得到的Ⅳ维独立拟随机序列作为种z)=+c1+⋯+Cp-1+c。(5)子，通过一定变换可得到服从特定概率分布的测量模式中：P为阶次；c为系数，取0或1。型输入量的离散抽样值，该离散抽样序列是Ⅳ维矢量对应每个多项式，取P个正奇数作为初值m一，X，i=1，⋯，Ⅳ，每一维由M个离散值组成且分别服从m，m<2‘，共Ⅳ组。按式(6)得到前P个vj值：测量模型中各输入量的概率分布类型。如果某些测=／2；=1，⋯，P(6)量不确定度源分量相关，则相应的离散抽样序列之间对于i>P，采用如下递归公式产生m和：应满足一定的相关系数。具体抽样方法可归纳为fm=2clmo2c2m①⋯①2cPm0m表1。L=c1①C2①⋯①cP①／2表1由独立拟随机序列生成各种分布抽样序列的方法Table1Themethodsofgeneratingprobabilitydistribufion(7)samplingsequencesfromindependentquasi式中：符号。表示2进制按位异或。randomsequences则Sobol序列的第n个数。为：0=bl1ob2v2ob3v30⋯(8)式中：⋯b3bb是i的2进制表示形式。置n=n+1，直至得到长度为M的拟随机数。Sobol序列只有在非常高维(260维以上)时相邻维之间才存在高度相关性，而测量不确定度源最多达数十个，因此不用考虑多维Sobol序列之间的Cholesky分解相关正态分布抽样序列相关性影响。但每一维Sobol序列仍然为低偏差的独立正态分布确定性点列，不能直接用于模拟多不确定度源的随线性一非线性变换具喜和机传递过程，需对Sobol序列进行随机化。本文采用随机抽取交换法对Sobol序列进行随机置乱，基本思想如下：2)拟蒙特卡罗传递和输H{设生成的Sobol序列为，⋯，，在区间[1，M]将得到的服从特定概率分布和相关系数的Ⅳ个离散中产生一个随机整数，记为r，表示抽取到，将与抽样矢量代入测量模型计算，输出个被测量l，的离散最后一个元素交换。为了避免元素被重复抽取，值Y，，r=l，⋯，，将Y，绘成直方图得到被测量y的近似第二次抽取时，从。，⋯，。中抽取一个，将，与概率密度函数并分别用平均值和标准偏差作为，，的最佳估计值和标准不确定度：一交换。依次类推，当第i次抽取时，从。，⋯，iM⋯中抽取，即在[1，M—i+1]中产生一个随机整Y==1∑Y(3)数r，将与⋯交换。抽取次即可生成均匀性厂一1————一好的随机抽样离散值。u()，)=“()√(y一)(4)图2显示了5000个服从U(0，1)均匀分布的二被测量包含区间的确定按JJF1059．2—2012《用维伪随机数和随机化Sobol序列及其直方图，可明显蒙特卡罗法评定测量不确定度》提供的方法进行。看出随机化Sobol序列更均匀地分布在(0，1)区间，用于测量不确定度传播时计算精度和收敛速度更好。服从具体分布的独立随机抽样序列可基于随机化3拟随机数产生及独立随机序列抽样方法Sobol序列根据分布函数逆变换方法和Moro方法获得。拟蒙特卡罗法中拟随机数的选取应该具有这样 1390仪器仪表学报第35卷表2采用不同方法的测量不确定度评定结果(表中M为实验次数)Table2Measurementuncertaintyevaluationresultswithdiferentmethods(MdenotesMonteCarlotrialnumber)由表2可知，由于非线性影响，GUM法得不到被测量的无偏估计值，本文方法的被测量估计值明显优于蒙特卡罗法，且随着实验次数增加，估计值接近理论真萋值。GUM2考虑了不确定度传播的高阶非线性项，标准不确定度结果与本文方法、理论值较为接近，而相同实验霉次数下本文方法的标准不确定度计算结果明显优于MCM，说明了本文方法采用低偏差序列作为所有抽样序列种子的优越性，相同计算精度下相比蒙特卡罗法本文方法只需更少的实验次数。此外，GUM方法得到的被测被测量／10。量95％包含区间与MCM、本文方法均存在较大差异，原图4基于QMCM的被测量概率密度直方图及因是GUM是假定被测量服从正态分布或t分布来计算基于3种方法的被测量包含区间包含区间，由图4可知，被测量概率分布明显呈非对称分Fig．4ProbabilitydensityhistogramwithQMCMandthe布，说明了GUM法的局限性。通过该例可得出以下4点coverageintervalswiththreekindsofmethods结论：1)测量模型明显呈非线性时，采用二阶或高阶Tay—lor展开的GUM法可以得到正确的标准不确定度；6实际应用2)输出量偏离正态分布或t分布时，虽然GUM法能得出正确的标准不确定度，但扩展不确定度仍然会出现精密离心机用于提供一个高精度的加速度基准作为错误，应根据不确定度分量的概率分布理论分析或通过惯性加速度计的输入，从而实现对惯性加速度计的标校。(拟)蒙特卡罗方法来计算或验证包含区间；精密离心机旋转时，受动不平衡和主轴回转误差等因素3)蒙特卡罗法与拟蒙特卡罗法的计算精度均与影响，转盘会偏离水平位置在垂直面上下摆动，产生动态实验次数有关，可采用自适应蒙特卡罗来确定合适的俯仰失准角，该角度将使待标校的加速度计输入重力分实验次数。然而，拟蒙特卡罗方法采用低偏差序列来量，影响加速度计的标校精度。因此必须对动态俯仰失代替伪随机数，计算精度和稳定性更高，但是需要将准角进行实时测量并对精密离心机输出的加速度值进行拟随机数随机化并且要避免高维拟随机序列的相关补偿。性带来的干扰；当精密离心机转速稳定后，动态俯仰失准角也趋于4)采用GUM法的不确定度传播律可以很方便地处稳定，可近似看作静态测量，其测量原理如图5所示。在理测量模型输入量之间的相关性。然而，对于含相关性离心机定位平台正下方安装电容测微仪，转盘下表面的输入量的测量问题，当测量模型明显呈非线性、灵敏度系局部区域机械加工精度达到txm量级。缓慢旋转离心数不易求解或者输出量的概率分布偏离正态分布或t分机，电容测微仪在局部区域采集30个点求平均值后作为布时，本文方法是一种有效的尝试。传感器测头与转盘下表面之间的间距，J．。当精密离心机 第6期凌明祥等：含相关性的测量不确定度拟蒙特卡罗评定方法1391处于工作状态时，采用定点测量方法，使电容测微仪在局表3测量结果及分布类型(仅供本文方法的应用例)部区域采集30个点求平均值后作为离心机工作状态时Table3Measurementresultsanddistributiontypes对应的传感器测头与转盘下表面之间的间距。则动态(Onlyillustratefortheproposedmethod)俯仰失准角可由式(21)计算得到：项目L1L2R^r，，估计值98．17ml16．23m0．99541921m入=arctan()=arctaI(!)(21)标准不确定度0．53Dxm0．67umO．1Ixm一1mm』』分布类型均匀分布均匀分布正态分布式中：A是转盘动态俯仰失准角，是精密离心机工包含因子√1作半径，是精密离心机以非常低的转速旋转时转盘下相关性l和，J2相关：p一0．59，R独立表面与电容测微仪测头之间的距离，：是精密离心机工作状态时转盘下表面与测微仪测头之间的距离。取蒙特卡罗试验次数M=10。，按照本文所提方法进行动态俯仰失准角的测量不确定度评定，结果如图6所示。当工作半径的测量不确定度在0．1m～1mm变化时，动态俯仰失准角的测量不确定度评定结果不变，原因是测量模型式(21)中工作半径的灵敏度系数极小，即使其测量不确定度较大，对动态俯仰失准角的测量不确定度贡献也很小可以忽略不计。说明了所提测量方法的优越性。此外，所提测量方法的输人量之间的相关性有利于降低传递到被测量的不确定度。由于测量模型的非线性4，如3果3梧采2用4{G2将UM1婪法霉1进行O测量不确定度评定，虽然相关性处理很简单，但灵敏度系数求解复杂，而且输出量并不服从正态分布或t分布(见图6)，计算动态俯仰失准角的扩展不确定度时存在一定误差。图5动态俯仰失准角测量原理Fig．5Measurementprincipleofdynamicpitchingmisalignmentangle电容测微仪采用米铱cs系列，量程500Ixrn，测量信号采用低噪声屏蔽电缆传输到PXI机箱接口，PXI板卡内部时钟驱动24位A／D采集卡进行信号的变速率采集。机室温度控制在20±0．5℃，接地电阻不大于0．1n。分别在精密离心机缓慢旋转和工作状态(旋转角速度设置为10rad／s)下重复采样20次，取平均值分别作为传感器测头与转盘下表面之间的间距。和。由于图6动态俯仰失准角测量结果概率密度直方图和：来自同一个传感器和数采系统，传感器和数采系Fig．6Probabilitydensityhistogramofthemeasurement统的误差将同时影响，J，和L：。在评定动态俯仰失准角resultofdynamicpitchingmisalignmentangle的测量不确定度时必须考虑输入量之间的相关性。为此，采用式(22)计算￡和之间的相关系数。测量结果如表3所示。由于电容测微仪轴心很难与加速度计质7结论心对齐，因此工作半径引入的测量不确定度可能很大，根本文探讨了测量不确定度源存在相关性时的不确定据实际经验一般不超过1mm。度拟蒙特卡罗评定方法。包括评定流程和步骤、拟随机∑(～五。)(。一己)数产生和随机化、具有特定概率分布的相关随机变量抽p旦—一一(22)样序列产生方法。本文方法尤其适用于不确定度源为非式中：n为重复性测量次数，、分别为。、的第i次正态分布且相关的测量过程。算例表明本文方法比一般重复测量结果，、：分别为、的n次重复测量结果蒙特卡罗方法的收敛速度和精度更优。平均值。采用本文方法对精密离心机的动态俯仰失准角的测 1392仪器仪表学报第35卷量不确定度进行了评定。相对于GUM方法，本文方法能[S]．2013．处理输入量相关、测量模型复杂、输出量不服从正态分布JJF1059．2—2012．EvaluationandMonteCarlomethod的测量问题。forevaluationofmeasurementuncertainty[S]．2013．修正了产生非正态分布且相关的随机抽样序列的非[8]GUERRASIOR．Measurementuncertaintyofisotopo—loguefractionsinfluxomicsdeterminedviamassspee—线性变换公式，构造目标函数采用搜索方法对相关系数矩阵进行求解。本文方法适用于随机变量概率分布函数trometry[J]．AnalBioanalChem，2013，405：5133—5146．能求逆的情况，否则应考虑采用其他方法来产生含相关[9]GERDW．MIVHAELK．Evaluationofmeasurementun—性的随机变量抽样序列。certaintyanditsnumericalcalculationbyaMonteCarlo参考文献method[J]．MeasurementScienceandTechnology，2008，19：084009—1-4．JJF1059．1—2012．测量不确定度评定与表示[10]HALLBD．Evaluatingmethodsofcalculatingmeasure—[s]．2013．mentuncertainty[J]．Metrologia，2008，24：L5一L8．JJF1059．1—2012．Evaluationandexpressionofuncer—WENXL，ZHAOYB．AdaptiveMonteCarloandGUMtaintyinmeasurement[S]．2013．methodsfortheevaluationofmeasurementuncertaintyof[2]ISO／IEC98—3，Uncertaintyofmeasurement一胁3：thecylindricityerror[J]．PrecisionEngineering，2013，37：Guidetotheexpressionofuncertaintyinmeasurement856—864．[S]．1999．[12]王辉，沈钺．相关性条件下合成不确定度的MC方[3]张亦弛，齐晓辉，林茂六．基于协方差矩阵的VNA／法[J]．电测与仪表，2007，44(10)：37—39．NVNA／LSNA测量不确定度传播规律与分析方法研究WANGH，SHENY．Combineduncertaintybythetheory[J]．仪器仪表学报，2012，33(11)：2495—2502．0fMCintheconditionoftherelativityjJ1．ElectricalZHANGYCH，QIxH，LINML．Covariance—matrixMeasurement＆Instrumentation，2007，44(10)：37—39．basedmeasurementuncertaintypropagationlawandanal—ysismethodforVNA’NVNAandLANA『J]．Chinese[13]WENLLW，ZHAOYB．AdaptiveMonteCarloandGUMmethodsfortheevaluationofmeasurementuncer—JournalofScientificInstrument，2012，33(11)：taintyofcylindricityeror[J]．PrecisionEngineering，372495—2502．(2013)：856—864．[4]陈怀艳，曹芸，韩洁．基于蒙特卡罗法的测量不确定[14]COXMG，SIEBERTBR．TheuseofaMonteCarlo度评定[J]．电子测量与仪器学报，2011，25(4)：methodforevaluatinguncertaintyandexpandeduncer—301—308．CHENHY，CAOY，HANJ．Evaluationofuncertaintytainty．Metrologia，2006(43)：S178—188．inmeasurementbasedonaMonteCarlometll0dlJI．[15]李慎安．JJF1059—1999(~量不确定度评定与表示》JournalofElectronicMeasurementandInstrument，201l，讨论之廿四：相关系数与灵敏系数问的联系[J]．T25(4)：301—308．业计量，2009，19(5)：35—36．[5]赵玉华，袁锋，丁振良，等．基于合作目标的姿态测量LiSHAN．JJF1059—1999，‘EvaluationandExpressionofUncertaintyinMeasurement’discussing：Therelation—系统建模及精度的蒙特卡罗估计[J]．仪器仪表学shipbetweencorrelationcoeficientandsensitivityeoeffi—报，2010，31(8)：1873—1877．ZHAOYH，YUANF，DINGZHL，eta1．Modelingofcient[J]．IndustrialMeasurement，2009，19(5)：theattitudemeasurementsystembasedoncooperationtar—35—36．getanditsaccuracyestimationwithMonte-Carlosimula-[16]黄美发，景晖，匡兵．基于拟蒙特卡罗方法的测量不tion[J]．ChineseJournalofScientificInstrument，2010，确定度评定[J]．仪器仪表学报，2009，30(1)：31(8)：1873—1877．120．125．[6]王伟，宋明顺，陈意华，等．蒙特卡罗方法在复杂模型HUANGMF，JINGH。KUANGB，eta1．Measure—测量不确定度评定中的应用[J]．仪器仪表学报，mentuncertaintyevaluationbasedonquasiMonte—Carlo2008，29(7)：1446—1449．methodJ]．ChineseJournalofScientificInstrument，WANGW，SONGMSH，CHENYH，eta1．Application2009，30(1)：120—125．0fMonte—Carlomethodintuncertaintyevalu一[17]DIMOVI．GEORGIEVR．OSTROPMSKYT．Ad—ationofcomplicatedmodel[J]．ChineseJournalofScien—vancedalgorithmsformultidimensionalsensitivitystud—tificInstrument，2008，29(7)：146—149．iesoflarge—scaleairpollutionmodelsbasedonSobol[7]JJF1059．2—2012．用蒙特卡罗法评定测量不确定度sequences[J]．ComputersandMathematicswithAp— '