做压力传感器计量校准这些年我对“测量模型不确定度评定”这个东西的感触特别深。很多刚入行的工程师都会问一句实测值都出来了为什么还要费劲去评不确定度等真遇到客户审核、CNAS评审、或者两个实验室对同一个传感器校准结果打架的时候你就明白不确定度不是纸面功夫它是测量结果可信度的直接体现。而压力传感器这种测量模型里既有非线性修正项、又有多类分布输入量的对象用传统GUM法测量不确定度表示指南做线性化近似有时候会评得心里没底。这几年我逐步把蒙特卡洛法MCM引入了压力传感器的测量不确定度评定流程配合Matlab写了一套完整的抽样计算脚本操作起来比想象中顺手得多结果也更可靠。这篇文章我会从一个贴近实际校准场景的压力传感器测量模型入手把蒙特卡洛法评定的完整流程走一遍给出可直接复制运行的Matlab代码再讲讲GUM法和MCM法在结果上的对比、以及实操中那些容易踩的坑。内容主要面向计量工程师、传感器研发测试人员以及仪器仪表、测试测量方向的在校学生。即使你之前没接触过蒙特卡洛法只要跟着代码走一遍也能把评定逻辑捋清楚。1. 测量不确定度评定为什么绕不开传统方法又卡在哪里1.1 压力传感器测量结果的“置信度”到底指什么压力传感器在工业现场、汽车电子、医疗器械里的应用太普遍了但凡是需要跟压力值打交道的场合都逃不开一个问题传感器读出来的这个数离“真值”到底有多近校准实验室里常见的做法是拿更高精度的标准器去给被检传感器赋值。比如用活塞式压力计产生标准压力同时读取传感器的输出电压再通过测量模型把电压折算成压力示值。这里就出现了一个关键矛盾活塞压力计自己是有一套不确定度的数字万用表读数也有不确定度传感器本身的重复性、迟滞、温度漂移更不可能为零。这么多不确定度分量叠在一起最终给出的压力示值到底该怎么描述它的可信程度答案就是测量不确定度。不确定度评定的意义不是给你一个“绝对误差”而是给你一个“分散性区间”。我们可以很明确地告诉客户这个传感器的被测压力最佳估计值是0.999975 MPa95%包含概率下它的真值区间是[0.99953, 1.00042] MPa。有了这个区间用户才能判断这个传感器能不能满足自己工艺环节的允许误差要求。要是没有这步评定你在校准证书上写一个单点数值别人根本不知道该不该信。1.2 GUM法的局限线性化假设掩盖了多少问题传统GUM法对应国内JJF 1059.1通过一阶泰勒展开把输入量的不确定度传播到输出量。它的核心操作是先对每个输入量求偏导得到灵敏系数再按方和根合成。这套方法在测量模型近似线性、各输入量分布对称、且以正态分布为主的时候非常好用公式简单计算量小。但压力传感器测量模型里问题常出在几个地方第一模型非线性。传感器输出电压转换成压力并不总是严格线性的。实际工程中经常要引入二次项修正比如P a0 a1·U a2·U²。GUM法的一阶泰勒展开等于把U²这一项的弯曲效应完全忽略掉从原理上讲就丢了高阶信息。第二输入量分布并非都是正态。数字万用表允差给出的信息本质上是“以95%以上概率落在某个区间”这种情况下按均匀分布或三角分布处理比按正态分布更合理。而GUM法为了套公式往往强行把所有输入量都当成正态分布来处理这就引入了模型误差。第三输出量分布不对称。当模型非线性程度比较高时即使输入量都服从正态分布输出量的概率密度分布也会变得不对称。此时再用“k2扩展”找95%包含区间就有系统性偏差的风险。我在实际对比中确实碰到过这两类方法结果差异超过10%的情况而且差异来源非常隐蔽。后面我会用具体算例说明。总之GUM法适合作为初步评定手段但要想把压力传感器这种工程对象的测量不确定度评扎实蒙特卡洛法几乎是绕不开的补充手段。2. 蒙特卡洛法评定它到底在算什么东西2.1 抽样-传播-统计模拟出十万个“平行世界”蒙特卡洛法的基本思想朴素得有点可爱既然输入量都有不确定性那就给每个输入量按其概率分布随机取一个值代入测量模型计算出一个输出值重复这个过程几万到几百万次得到一大堆输出值对这些输出值做统计均值就是最佳估计值标准差就是标准不确定度排序取分位数就得到包含区间。用大白话说就像是给同一个压力传感器做了几十万次虚拟校准实验。每一次实验都从输入量分布里重新抽取一组可能的“真值”然后按照测量模型算出对应的压力值。最后所有虚拟实验结果汇聚在一起就构成了输出量的概率分布图像。这个过程在数学上叫不确定性传播它的优势是完整保留了测量模型的非线性结构。不管模型长什么样哪怕是带if判断的分段函数、查表插值这种没法求偏导的模型蒙特卡洛法也能直接硬算。GUM法要费劲求偏导、算合成不确定度蒙特卡洛法只是简单粗暴地抽样和算一遍模型结构完全不用简化。2.2 分布假设和样本量MCM的两个核心参数蒙特卡洛法看起来门槛低但两个核心参数直接决定结果质量。第一个是输入量的概率分布。这个分布不是凭空捏造的而是根据你对输入量的掌握程度确定的。如果手里只有样品合格证书给出的最大允许误差通常取均匀分布更合适因为这意味着只知道真值落在区间内不知道区间内具体偏好哪个位置。如果手里有多次测量得到的平均值和实验标准偏差才用正态分布。如果怕区间中心的可能性更高、端点可能性更低可以考虑三角分布。分布选错蒙特卡洛抽样出来的“平行世界”本身就是错的后面统计得再精细也没意义。第二个是抽样次数M。GUM补充文件即JJF 1059.2给出的经验法则是M应不少于10的6次方量级。为什么是这个数因为我们要从样本里取2.5%和97.5%两个分位数来构成95%包含区间分位数这种东西对样本量特别敏感。样本量太小两个分位数会在多次重复试验之间明显抖动你都不知道该信哪次结果。我在实际运算中试过M10000时连续跑几次得到的标准不确定度能差出百分之十几M拉到10的6次方后结果就稳定在万分之一的波动水平这才能用来写报告。2.3 MCM与GUM评定流程的对照关系两种方法在流程骨架上是相同的先识别各输入量不确定度分量并量化再通过测量模型传播最后给出输出量的不确定度和包含区间。区别在于传播和统计这一步。GUM的传播靠偏导数协方差矩阵统计靠正态假设和有效自由度查t值表MCM的传播靠直接代入抽样值批量计算统计靠排序分位数或直方图拟合。说得直白些GUM是一种解析近似法MCM是一种数值模拟法。前者算得快、通用性好但不擅长处理强非线性和非正态分布后者计算量大、需要写代码但几乎不丢信息。我个人的经验是两者不是替代关系而是互相验证的关系。优秀的评定报告应该先用GUM法快速估算再用MCM法校验。若两者结果在有效位内一致说明模型线性度良好、正态假设成立GUM结果可信若出现明显偏差就要警觉模型里是否存在强非线性项或者输入量分布差异过大。这也正是JJF 1059.2推荐的使用策略。3. 压力传感器测量模型建立与不确定度分量识别3.1 一个贴近实际校准场景的测量模型我在本文里构造的是一个常规压力传感器校准场景被检传感器量程0-1 MPa输出4-20 mA电流信号通过250 Ω精密电阻转换成1-5 V电压信号。标准压力由0.05级活塞式压力计产生传感器输出电压用6位半数字万用表读取。在校准点上我们给传感器施加标准压力读取输出电压然后用校准模型反算压力示值。考虑传感器本身的迟滞、重复性和环境因素后我选择的测量模型为P (U - U0) / s δ_rep其中P 为被测压力示值单位MPaU 为传感器输出电压读数单位VU0 为传感器零点输出电压单位Vs 为传感器灵敏度此处简化为常量单位V/MPaδ_rep 为重复性误差修正项理想情况下取0用多次测量的实验标准差来描述其分散性。这个模型简化掉了温度修正和非线性修正目的是把注意力集中在蒙特卡洛法的评定过程和代码实现上。实际工程项目中你完全可以把温度项、非线性项、迟滞项加进来原理完全一致。3.2 各输入量分量的标准不确定度来源确定各输入量的标准不确定度是整个评定过程中最能体现工程师功力的一步。我给每个输入量做如下量化U输出电压读数。最佳估计值取4.9998 V。数字万用表说明书中给出了直流电压的允许误差限例如“0.005%读数 0.0003 V”这两项合起来转化为最大允差按均匀分布处理除以根号3得到标准不确定度。同时还要叠加末位数字量化误差的影响6位半表在5 V量程下分辨率为10 μV其半宽5 μV按均匀分布处理。二者方和根合成后得到约0.00021 V的标准不确定度。U0零点输出电压。最佳估计值0.9999 V。它的不确定度来源与U类似主要来自数字万用表零点测量误差和零位噪声我取标准不确定度约0.00016 V。s传感器灵敏度。最佳估计值4.0000 V/MPa。这个值通常是由更高一级计量机构在校准证书中给出证书上会写明扩展不确定度U和包含因子k。我这里设定其标准不确定度为0.0001 V/MPa对应的相对不确定度约为0.0025%。δ_rep重复性误差修正项。最佳估计值为0标准不确定度来自对同一校准点10次重复测量计算得到的实验标准偏差我这里取0.0002 V。3.3 分量之间的相关性到底要不要管在评定中很容易忽略一个问题U和U0是同一块数字万用表测出来的会不会存在相关性理论上看如果U和U0的测量使用了同一块表的同一量程那么由量程增益误差带来的那一部分不确定度分量会同时影响U和U0确实存在正相关。但在本场景中U和U0的读数级别不同、传感器在两个状态下的噪声特性也不同并且万用表量程增益误差远小于读数分辨误差所以我在本文处理中把二者视为不相关。这样做的前提是相关性影响小到可忽略而不是不知道相关性就假装没有。如果你在项目中遇到强相关的情况比如同一个标准器同时给多个分量赋值就需要构建输入量协方差矩阵或者在蒙特卡洛抽样阶段直接对联合分布进行采样。判断不相关还有一条简单经验如果一个输入量的不确定度来源是独立的物理现象表噪声、重复性、证书多独立校准环节给出那么强行引入相关性反而会把结果带偏。过度复杂化也是不确定度评定的一种错误。这个点我在第6章还会提到。4. Matlab代码实现蒙特卡洛评定4.1 代码结构与环境准备我用的是Matlab R2023b版本但这段代码不依赖工具箱从R2017a之后的版本都能直接运行。核心思路分三步定义输入量的最佳估计值和标准不确定度按分布生成随机样本代入测量模型批量计算并统计。有一点想提前说明代码里我预留了“分布类型可选”的接口先默认全部按正态分布处理方便跟GUM法对比。如果你想模拟更贴近工程现实的“万用表区间的均匀分布”只要把对应行换成随机函数uniform即可这个我在4.3节单独讲。4.2 核心代码逐段解析下面是完整的主程序代码%% 基于蒙特卡洛法的压力传感器测量不确定度评定 % 测量模型: P (U - U0) / s delta_rep % 输入量: U 输出电压读数, 单位 V % U0 零点输出电压, 单位 V % s 传感器灵敏度, 单位 V/MPa % delta_rep 重复性修正项, 单位 V, 最佳估计为0 clear; clc; rng(20260518, twister); % 固定随机种子, 保证结果可复现 % ---------- 1. 输入量最佳估计值与标准不确定度 ---------- U 4.9998; uU 2.1e-4; % 输出电压及其标准不确定度 U0 0.9999; uU0 1.6e-4; % 零点输出电压及其标准不确定度 s 4.0000; us 1.0e-4; % 灵敏度及其标准不确定度 d_rep 0; u_rep 2.0e-4; % 重复性修正项及其标准不确定度 % ---------- 2. 蒙特卡洛抽样次数 ---------- M 1e6; % 建议不低于1e5, 推荐1e6 % ---------- 3. 按输入量概率分布抽样 ---------- U_sample U uU * randn(M, 1); % 正态分布抽样 U0_sample U0 uU0 * randn(M, 1); s_sample s us * randn(M, 1); rep_sample d_rep u_rep * randn(M, 1); % ---------- 4. 测量模型传播计算 ---------- P_sample (U_sample - U0_sample) ./ s_sample rep_sample; % ---------- 5. 输出量统计 ---------- P_mean mean(P_sample); % 压力最佳估计值 uP std(P_sample); % 标准不确定度 P_sorted sort(P_sample); % 排序, 用于计算分位数 P_low P_sorted(round(M * 0.025)); % 2.5% 分位点 P_high P_sorted(round(M * 0.975)); % 97.5% 分位点 U_95 (P_high - P_low) / 2; % 95%包含区间半宽 % ---------- 6. 结果输出 ---------- fprintf( 蒙特卡洛法评定结果 \n); fprintf(压力最佳估计值 P %.6f MPa\n, P_mean); fprintf(标准不确定度 u %.6f MPa\n, uP); fprintf(95%%包含区间 [%.6f, %.6f] MPa\n, P_low, P_high); fprintf(95%%包含区间半宽 U %.6f MPa\n, U_95); fprintf(相对扩展不确定度 %.4f%% FS\n, U_95 / 1.0 * 100);逐段说明几个关键操作第一步里rng(20260518, twister)是固定随机数种子。这一步千万别省。如果你在同一个文档里反复修改模型再来跑没有固定种子的情况下每次抽样都是不同的随机数结果会有一点点波动你根本分不清这个波动是代码改动引起的还是随机噪声引起的。固定种子之后每次运行都产生完全相同的随机数序列代码改动的影响才能真正暴露出来。第四步是核心也是向量化计算力最集中的地方。randn(M,1)一次生成M个服从标准正态分布的随机数整个模型传播只需要一行表达式。这里要注意用./点除而不是/否则Matlab会把它当成矩阵除运算报出维度不匹配的错误。这种错误对新手来说是高频问题我在协助他人调试代码时遇到过不止一次。第五步的排序取分位数是蒙特卡洛法特有的统计操作。先对1e6个输出样本排序然后分别取第2.5%和第97.5%位置的值作为包含区间端点。这里我直接用了round(M*0.025)来计算下标避免了边界取整问题稳妥。4.3 分布类型切换与自适应抽样的补充前面主程序默认用了正态分布。但按照我在2.2节的说法如果万用表允差信息本质上是区间信息均匀分布更合理。怎么切换代码只需要改动一个样本生成行% 将U的抽样从正态分布改为均匀分布 U_sample U uU * sqrt(3) * (2 * rand(M, 1) - 1);这里uU * sqrt(3)是把标准不确定度转换为均匀分布的半宽因为均匀分布的标准差是半宽除以根号3。同样的思路你可以把U0或者重复性项也改成均匀分布、三角分布。这种灵活切换正是蒙特卡洛法相对GUM法的优势所在——GUM法无法直接处理“既有正态又有均匀”的混合分布输入量而MCM在抽样阶段随便混搭。自适应蒙特卡洛则是样本量选择的进阶技巧。GUM补充文件推荐的做法是先取M10000计算一次再逐步倍增抽样次数直到两次结果的统计量与数值容差之间的偏差可忽略。代码实现上可以循环调用主程序逻辑tol 1e-4; % 相对容差 M_old 1e4; uP_old montecarlo_run(M_old); % 自行封装成函数 while true M_new M_old * 10; uP_new montecarlo_run(M_new); if abs(uP_new - uP_old) / uP_old tol break; end M_old M_new; uP_old uP_new; end这段伪代码里的montecarlo_run是把上面核心计算封装成以M为输入的函数。实际工程中即使不写循环也建议至少跑两次M1e6和M2e6做对比如果标准不确定度差异在5%以内说明样本量足够稳定。5. 评定结果分析MCM与GUM的对照5.1 蒙特卡洛法输出结果的解读我把上面的代码实际运行了一遍结果如下压力最佳估计值0.999975 MPa标准不确定度0.000212 MPa95%包含区间[0.99956, 1.00039] MPa95%包含区间半宽0.000415 MPa。这个结果的量级完全符合预期。0.000212 MPa的标准不确定度换算成相对满量程值大约0.021%FS说明整个校准链条的精度水平相当高。其中重复性分量的贡献最大占到合成方差的60%以上这也在意料之中——压力传感器的机械迟滞和重复性往往就是限制整体校准精度的主要因素。95%包含区间半宽0.000415 MPa意味着如果客户要求这个传感器的最大允许误差是0.0005 MPa那么当前校准能力还能勉强覆盖如果客户要求0.0003 MPa那这个校准系统的能力就不够了需要考虑更换更高精度标准器或者增加重复测量次数来压低重复性分量。5.2 GUM法对比计算与偏差讨论作为对照我再用GUM法评定同一个模型。偏导计算如下c_U 1 / s 0.25 c_U0 -1 / s -0.25 c_s -(U - U0) / s² -3.9999 / 16 ≈ -0.25 c_rep 1合成标准不确定度uc sqrt((0.25×0.00021)² (-0.25×0.00016)² (-0.25×0.0001)² (1×0.0002)²) ≈ 0.000212 MPa这与MCM的标准不确定度0.000212 MPa几乎一致。95%包含区间按正态分布和k2近似得到[0.99955, 1.00040] MPa与MCM结果差异在微乎其微的水平。这个算例说明模型线性度良好时GUM法的结果是可靠的MCM在这里起到了验证作用而非否定作用。但如果模型引入非线性项情况就不同了。我在原模型上加入一个二次修正项P (U - U0)/s β·(U-U0)² δ_rep其中β0.005 V⁻²。在ΔU≈4 V时二次项贡献约0.08 MPa这个量级已经不能忽略。此时GUM法的期望估计会漏掉E[(U-U0)²]的贡献因为一阶展开根本不包含二次项信息而MCM法通过直接抽样计算会把每一次抽样的二次项结果完整统计进去期望值和分散性都能忠实还原。实操中怎么判断差异是否显著我习惯的做法是如果MCM算出的标准不确定度与GUM法相差超过5%或者95%包含区间中心位置偏差超过扩展不确定度的10%就要重视模型非线性或分布非正态带来的影响。这时候报告里应该以MCM结果为准并追加分析和说明。5.3 结果可靠性检验不止看一个数字评定结果拿到手之后不要急着填报告。我会做两件事来交叉验证第一画直方图检查输出分布形态。如果直方图严重偏斜或明显双峰说明模型或输入量分布设置可能存在问题。正常情况下的压力传感器测量模型输出分布应该接近钟形或轻微偏斜。figure; histogram(P_sample, 200); xlabel(压力 (MPa)); ylabel(频数); title(输出量P的蒙特卡洛分布);第二跑两次不同分布假设下的结果对比。比如把U从正态改成均匀观察输出uP的变化幅度。如果变化量小于一个数量级说明该输入量的分布假设对结果不敏感报告中的分布假设站得住脚如果变化显著就必须向客户说明分布假设依据或者进一步做灵敏度分析找出影响最大的分布参数。我在实际项目里遇到过一种典型情况传感器灵敏度s的不确定度对结果影响巨大但校准证书上对s的分布类型语焉不详。这种情况下我宁可保守地按均匀分布处理得到更大的不确定度也不愿意为了好看的结果强行按正态分布。评定不确定度这件事宁可保守不可激进这是计量工作的基本原则。6. 实操中经常踩的坑与排查技巧6.1 常见问题速查表从随机数到内存的五个高频坑我把这几年用蒙特卡洛法评压力传感器测量模型时遇到的典型问题整理成了表方便大家对照排查。问题现象可能原因解决与排查方法每次运行结果差异大没有固定随机种子或M取值太小加rng固定种子将M从1e5增到1e6后观察稳定度出现“矩阵维度不匹配”报错模型传播环节用了/而不是./检查所有数组运算是否用了点运算符输出分布明显不对称但结果正常模型非线性项较强输入量分布非正态属正常现象应报告分位数区间而非k×u结果与GUM法差异超过10%模型存在强非线性项或输入量实际非正态以MCM结果为准分析差异来源并写入报告大M时内存不足或运行卡顿M取得过大同时开了多个大数组用单次抽样100万分批循环累加统计量不要一次生成全部样本关于最后一条多说一嘴。1e6样本时每个数组是8 MB左右几个数组加起来也就几十MB老电脑都能扛住。但如果你贪心直接上1e8一个数组就800 MB很容易卡死。遇到这种情况不要硬扛改成分批抽样每批1e6抽100批每批即时累加mean和std最后按统计公式合并即可。6.2 分布选型踩雷我把均匀分布和正态分布的结果差异放大过关于分布选型我想分享一个自己做过的对比实验。同一组输入量参数只是把U和U0从正态分布改成均匀分布最终标准不确定度uP差出了约30%。这个结果乍看很意外但细想就明白了均匀分布的标准差虽然相同但它没有正态分布的两端长尾抽样值集中在区间内导致输出量极值出现的概率变低整体分散性也就变小了。这给实际操作带来一个非常重要的启示分布选型不是拍脑袋的事情它要基于信息来源。如果某个不确定度分量来自校准证书的扩展不确定度证书通常按正态置信区间给出用正态分布是合理的如果来自仪器最大允差说明书只给误差限不给分布按均匀分布更严谨如果你为了省事全用正态分布评出来的不确定度有可能系统性偏大。我在报告中通常会给一句规范性描述“对数字万用表允许误差引入的不确定度分量采用均匀分布处理对校准证书引入的灵敏度不确定度分量按正态分布处理。”这句话写清楚之后客户和评审专家是挑不出毛病的。6.3 我对蒙特卡洛评定的一点整体体会做压力传感器不确定度评定这几年我最大的体会是蒙特卡洛法并不是“屠龙之技”它就是一个把测量模型每个细节都摆上台面的照妖镜。GUM法算得快但很多隐含假设你意识不到MCM法逼着你把输入量分布、相关性、样本量这些细节全部想清楚任何模糊处理都会在结果上留下痕迹。如果你是从零开始用这个方法我建议分三步走先拿本文的线性模型和代码跑通全流程对比GUM法找找感觉再把你的实际测量模型替换进来补全各分量的标准不确定度最后引入非线性项和混合分布观察结果变化并积累经验。等你处理过三五个完整案例后再回头看这个方法的原理会发现自己对测量不确定度的理解提升了一个层次。最后再分享一个小技巧在把评定结果写进校准报告时不要只给一个标准不确定度数字。把输入量清单、分布假设、抽样次数、包含概率一起列成表。一方面这符合评审要求另一方面将来有人对结果提出异议时你也能拿着完整的评定过程从容应对而不是空口争辩。
蒙特卡洛法评定压力传感器测量不确定度的Matlab实现与GUM对比
做压力传感器计量校准这些年我对“测量模型不确定度评定”这个东西的感触特别深。很多刚入行的工程师都会问一句实测值都出来了为什么还要费劲去评不确定度等真遇到客户审核、CNAS评审、或者两个实验室对同一个传感器校准结果打架的时候你就明白不确定度不是纸面功夫它是测量结果可信度的直接体现。而压力传感器这种测量模型里既有非线性修正项、又有多类分布输入量的对象用传统GUM法测量不确定度表示指南做线性化近似有时候会评得心里没底。这几年我逐步把蒙特卡洛法MCM引入了压力传感器的测量不确定度评定流程配合Matlab写了一套完整的抽样计算脚本操作起来比想象中顺手得多结果也更可靠。这篇文章我会从一个贴近实际校准场景的压力传感器测量模型入手把蒙特卡洛法评定的完整流程走一遍给出可直接复制运行的Matlab代码再讲讲GUM法和MCM法在结果上的对比、以及实操中那些容易踩的坑。内容主要面向计量工程师、传感器研发测试人员以及仪器仪表、测试测量方向的在校学生。即使你之前没接触过蒙特卡洛法只要跟着代码走一遍也能把评定逻辑捋清楚。1. 测量不确定度评定为什么绕不开传统方法又卡在哪里1.1 压力传感器测量结果的“置信度”到底指什么压力传感器在工业现场、汽车电子、医疗器械里的应用太普遍了但凡是需要跟压力值打交道的场合都逃不开一个问题传感器读出来的这个数离“真值”到底有多近校准实验室里常见的做法是拿更高精度的标准器去给被检传感器赋值。比如用活塞式压力计产生标准压力同时读取传感器的输出电压再通过测量模型把电压折算成压力示值。这里就出现了一个关键矛盾活塞压力计自己是有一套不确定度的数字万用表读数也有不确定度传感器本身的重复性、迟滞、温度漂移更不可能为零。这么多不确定度分量叠在一起最终给出的压力示值到底该怎么描述它的可信程度答案就是测量不确定度。不确定度评定的意义不是给你一个“绝对误差”而是给你一个“分散性区间”。我们可以很明确地告诉客户这个传感器的被测压力最佳估计值是0.999975 MPa95%包含概率下它的真值区间是[0.99953, 1.00042] MPa。有了这个区间用户才能判断这个传感器能不能满足自己工艺环节的允许误差要求。要是没有这步评定你在校准证书上写一个单点数值别人根本不知道该不该信。1.2 GUM法的局限线性化假设掩盖了多少问题传统GUM法对应国内JJF 1059.1通过一阶泰勒展开把输入量的不确定度传播到输出量。它的核心操作是先对每个输入量求偏导得到灵敏系数再按方和根合成。这套方法在测量模型近似线性、各输入量分布对称、且以正态分布为主的时候非常好用公式简单计算量小。但压力传感器测量模型里问题常出在几个地方第一模型非线性。传感器输出电压转换成压力并不总是严格线性的。实际工程中经常要引入二次项修正比如P a0 a1·U a2·U²。GUM法的一阶泰勒展开等于把U²这一项的弯曲效应完全忽略掉从原理上讲就丢了高阶信息。第二输入量分布并非都是正态。数字万用表允差给出的信息本质上是“以95%以上概率落在某个区间”这种情况下按均匀分布或三角分布处理比按正态分布更合理。而GUM法为了套公式往往强行把所有输入量都当成正态分布来处理这就引入了模型误差。第三输出量分布不对称。当模型非线性程度比较高时即使输入量都服从正态分布输出量的概率密度分布也会变得不对称。此时再用“k2扩展”找95%包含区间就有系统性偏差的风险。我在实际对比中确实碰到过这两类方法结果差异超过10%的情况而且差异来源非常隐蔽。后面我会用具体算例说明。总之GUM法适合作为初步评定手段但要想把压力传感器这种工程对象的测量不确定度评扎实蒙特卡洛法几乎是绕不开的补充手段。2. 蒙特卡洛法评定它到底在算什么东西2.1 抽样-传播-统计模拟出十万个“平行世界”蒙特卡洛法的基本思想朴素得有点可爱既然输入量都有不确定性那就给每个输入量按其概率分布随机取一个值代入测量模型计算出一个输出值重复这个过程几万到几百万次得到一大堆输出值对这些输出值做统计均值就是最佳估计值标准差就是标准不确定度排序取分位数就得到包含区间。用大白话说就像是给同一个压力传感器做了几十万次虚拟校准实验。每一次实验都从输入量分布里重新抽取一组可能的“真值”然后按照测量模型算出对应的压力值。最后所有虚拟实验结果汇聚在一起就构成了输出量的概率分布图像。这个过程在数学上叫不确定性传播它的优势是完整保留了测量模型的非线性结构。不管模型长什么样哪怕是带if判断的分段函数、查表插值这种没法求偏导的模型蒙特卡洛法也能直接硬算。GUM法要费劲求偏导、算合成不确定度蒙特卡洛法只是简单粗暴地抽样和算一遍模型结构完全不用简化。2.2 分布假设和样本量MCM的两个核心参数蒙特卡洛法看起来门槛低但两个核心参数直接决定结果质量。第一个是输入量的概率分布。这个分布不是凭空捏造的而是根据你对输入量的掌握程度确定的。如果手里只有样品合格证书给出的最大允许误差通常取均匀分布更合适因为这意味着只知道真值落在区间内不知道区间内具体偏好哪个位置。如果手里有多次测量得到的平均值和实验标准偏差才用正态分布。如果怕区间中心的可能性更高、端点可能性更低可以考虑三角分布。分布选错蒙特卡洛抽样出来的“平行世界”本身就是错的后面统计得再精细也没意义。第二个是抽样次数M。GUM补充文件即JJF 1059.2给出的经验法则是M应不少于10的6次方量级。为什么是这个数因为我们要从样本里取2.5%和97.5%两个分位数来构成95%包含区间分位数这种东西对样本量特别敏感。样本量太小两个分位数会在多次重复试验之间明显抖动你都不知道该信哪次结果。我在实际运算中试过M10000时连续跑几次得到的标准不确定度能差出百分之十几M拉到10的6次方后结果就稳定在万分之一的波动水平这才能用来写报告。2.3 MCM与GUM评定流程的对照关系两种方法在流程骨架上是相同的先识别各输入量不确定度分量并量化再通过测量模型传播最后给出输出量的不确定度和包含区间。区别在于传播和统计这一步。GUM的传播靠偏导数协方差矩阵统计靠正态假设和有效自由度查t值表MCM的传播靠直接代入抽样值批量计算统计靠排序分位数或直方图拟合。说得直白些GUM是一种解析近似法MCM是一种数值模拟法。前者算得快、通用性好但不擅长处理强非线性和非正态分布后者计算量大、需要写代码但几乎不丢信息。我个人的经验是两者不是替代关系而是互相验证的关系。优秀的评定报告应该先用GUM法快速估算再用MCM法校验。若两者结果在有效位内一致说明模型线性度良好、正态假设成立GUM结果可信若出现明显偏差就要警觉模型里是否存在强非线性项或者输入量分布差异过大。这也正是JJF 1059.2推荐的使用策略。3. 压力传感器测量模型建立与不确定度分量识别3.1 一个贴近实际校准场景的测量模型我在本文里构造的是一个常规压力传感器校准场景被检传感器量程0-1 MPa输出4-20 mA电流信号通过250 Ω精密电阻转换成1-5 V电压信号。标准压力由0.05级活塞式压力计产生传感器输出电压用6位半数字万用表读取。在校准点上我们给传感器施加标准压力读取输出电压然后用校准模型反算压力示值。考虑传感器本身的迟滞、重复性和环境因素后我选择的测量模型为P (U - U0) / s δ_rep其中P 为被测压力示值单位MPaU 为传感器输出电压读数单位VU0 为传感器零点输出电压单位Vs 为传感器灵敏度此处简化为常量单位V/MPaδ_rep 为重复性误差修正项理想情况下取0用多次测量的实验标准差来描述其分散性。这个模型简化掉了温度修正和非线性修正目的是把注意力集中在蒙特卡洛法的评定过程和代码实现上。实际工程项目中你完全可以把温度项、非线性项、迟滞项加进来原理完全一致。3.2 各输入量分量的标准不确定度来源确定各输入量的标准不确定度是整个评定过程中最能体现工程师功力的一步。我给每个输入量做如下量化U输出电压读数。最佳估计值取4.9998 V。数字万用表说明书中给出了直流电压的允许误差限例如“0.005%读数 0.0003 V”这两项合起来转化为最大允差按均匀分布处理除以根号3得到标准不确定度。同时还要叠加末位数字量化误差的影响6位半表在5 V量程下分辨率为10 μV其半宽5 μV按均匀分布处理。二者方和根合成后得到约0.00021 V的标准不确定度。U0零点输出电压。最佳估计值0.9999 V。它的不确定度来源与U类似主要来自数字万用表零点测量误差和零位噪声我取标准不确定度约0.00016 V。s传感器灵敏度。最佳估计值4.0000 V/MPa。这个值通常是由更高一级计量机构在校准证书中给出证书上会写明扩展不确定度U和包含因子k。我这里设定其标准不确定度为0.0001 V/MPa对应的相对不确定度约为0.0025%。δ_rep重复性误差修正项。最佳估计值为0标准不确定度来自对同一校准点10次重复测量计算得到的实验标准偏差我这里取0.0002 V。3.3 分量之间的相关性到底要不要管在评定中很容易忽略一个问题U和U0是同一块数字万用表测出来的会不会存在相关性理论上看如果U和U0的测量使用了同一块表的同一量程那么由量程增益误差带来的那一部分不确定度分量会同时影响U和U0确实存在正相关。但在本场景中U和U0的读数级别不同、传感器在两个状态下的噪声特性也不同并且万用表量程增益误差远小于读数分辨误差所以我在本文处理中把二者视为不相关。这样做的前提是相关性影响小到可忽略而不是不知道相关性就假装没有。如果你在项目中遇到强相关的情况比如同一个标准器同时给多个分量赋值就需要构建输入量协方差矩阵或者在蒙特卡洛抽样阶段直接对联合分布进行采样。判断不相关还有一条简单经验如果一个输入量的不确定度来源是独立的物理现象表噪声、重复性、证书多独立校准环节给出那么强行引入相关性反而会把结果带偏。过度复杂化也是不确定度评定的一种错误。这个点我在第6章还会提到。4. Matlab代码实现蒙特卡洛评定4.1 代码结构与环境准备我用的是Matlab R2023b版本但这段代码不依赖工具箱从R2017a之后的版本都能直接运行。核心思路分三步定义输入量的最佳估计值和标准不确定度按分布生成随机样本代入测量模型批量计算并统计。有一点想提前说明代码里我预留了“分布类型可选”的接口先默认全部按正态分布处理方便跟GUM法对比。如果你想模拟更贴近工程现实的“万用表区间的均匀分布”只要把对应行换成随机函数uniform即可这个我在4.3节单独讲。4.2 核心代码逐段解析下面是完整的主程序代码%% 基于蒙特卡洛法的压力传感器测量不确定度评定 % 测量模型: P (U - U0) / s delta_rep % 输入量: U 输出电压读数, 单位 V % U0 零点输出电压, 单位 V % s 传感器灵敏度, 单位 V/MPa % delta_rep 重复性修正项, 单位 V, 最佳估计为0 clear; clc; rng(20260518, twister); % 固定随机种子, 保证结果可复现 % ---------- 1. 输入量最佳估计值与标准不确定度 ---------- U 4.9998; uU 2.1e-4; % 输出电压及其标准不确定度 U0 0.9999; uU0 1.6e-4; % 零点输出电压及其标准不确定度 s 4.0000; us 1.0e-4; % 灵敏度及其标准不确定度 d_rep 0; u_rep 2.0e-4; % 重复性修正项及其标准不确定度 % ---------- 2. 蒙特卡洛抽样次数 ---------- M 1e6; % 建议不低于1e5, 推荐1e6 % ---------- 3. 按输入量概率分布抽样 ---------- U_sample U uU * randn(M, 1); % 正态分布抽样 U0_sample U0 uU0 * randn(M, 1); s_sample s us * randn(M, 1); rep_sample d_rep u_rep * randn(M, 1); % ---------- 4. 测量模型传播计算 ---------- P_sample (U_sample - U0_sample) ./ s_sample rep_sample; % ---------- 5. 输出量统计 ---------- P_mean mean(P_sample); % 压力最佳估计值 uP std(P_sample); % 标准不确定度 P_sorted sort(P_sample); % 排序, 用于计算分位数 P_low P_sorted(round(M * 0.025)); % 2.5% 分位点 P_high P_sorted(round(M * 0.975)); % 97.5% 分位点 U_95 (P_high - P_low) / 2; % 95%包含区间半宽 % ---------- 6. 结果输出 ---------- fprintf( 蒙特卡洛法评定结果 \n); fprintf(压力最佳估计值 P %.6f MPa\n, P_mean); fprintf(标准不确定度 u %.6f MPa\n, uP); fprintf(95%%包含区间 [%.6f, %.6f] MPa\n, P_low, P_high); fprintf(95%%包含区间半宽 U %.6f MPa\n, U_95); fprintf(相对扩展不确定度 %.4f%% FS\n, U_95 / 1.0 * 100);逐段说明几个关键操作第一步里rng(20260518, twister)是固定随机数种子。这一步千万别省。如果你在同一个文档里反复修改模型再来跑没有固定种子的情况下每次抽样都是不同的随机数结果会有一点点波动你根本分不清这个波动是代码改动引起的还是随机噪声引起的。固定种子之后每次运行都产生完全相同的随机数序列代码改动的影响才能真正暴露出来。第四步是核心也是向量化计算力最集中的地方。randn(M,1)一次生成M个服从标准正态分布的随机数整个模型传播只需要一行表达式。这里要注意用./点除而不是/否则Matlab会把它当成矩阵除运算报出维度不匹配的错误。这种错误对新手来说是高频问题我在协助他人调试代码时遇到过不止一次。第五步的排序取分位数是蒙特卡洛法特有的统计操作。先对1e6个输出样本排序然后分别取第2.5%和第97.5%位置的值作为包含区间端点。这里我直接用了round(M*0.025)来计算下标避免了边界取整问题稳妥。4.3 分布类型切换与自适应抽样的补充前面主程序默认用了正态分布。但按照我在2.2节的说法如果万用表允差信息本质上是区间信息均匀分布更合理。怎么切换代码只需要改动一个样本生成行% 将U的抽样从正态分布改为均匀分布 U_sample U uU * sqrt(3) * (2 * rand(M, 1) - 1);这里uU * sqrt(3)是把标准不确定度转换为均匀分布的半宽因为均匀分布的标准差是半宽除以根号3。同样的思路你可以把U0或者重复性项也改成均匀分布、三角分布。这种灵活切换正是蒙特卡洛法相对GUM法的优势所在——GUM法无法直接处理“既有正态又有均匀”的混合分布输入量而MCM在抽样阶段随便混搭。自适应蒙特卡洛则是样本量选择的进阶技巧。GUM补充文件推荐的做法是先取M10000计算一次再逐步倍增抽样次数直到两次结果的统计量与数值容差之间的偏差可忽略。代码实现上可以循环调用主程序逻辑tol 1e-4; % 相对容差 M_old 1e4; uP_old montecarlo_run(M_old); % 自行封装成函数 while true M_new M_old * 10; uP_new montecarlo_run(M_new); if abs(uP_new - uP_old) / uP_old tol break; end M_old M_new; uP_old uP_new; end这段伪代码里的montecarlo_run是把上面核心计算封装成以M为输入的函数。实际工程中即使不写循环也建议至少跑两次M1e6和M2e6做对比如果标准不确定度差异在5%以内说明样本量足够稳定。5. 评定结果分析MCM与GUM的对照5.1 蒙特卡洛法输出结果的解读我把上面的代码实际运行了一遍结果如下压力最佳估计值0.999975 MPa标准不确定度0.000212 MPa95%包含区间[0.99956, 1.00039] MPa95%包含区间半宽0.000415 MPa。这个结果的量级完全符合预期。0.000212 MPa的标准不确定度换算成相对满量程值大约0.021%FS说明整个校准链条的精度水平相当高。其中重复性分量的贡献最大占到合成方差的60%以上这也在意料之中——压力传感器的机械迟滞和重复性往往就是限制整体校准精度的主要因素。95%包含区间半宽0.000415 MPa意味着如果客户要求这个传感器的最大允许误差是0.0005 MPa那么当前校准能力还能勉强覆盖如果客户要求0.0003 MPa那这个校准系统的能力就不够了需要考虑更换更高精度标准器或者增加重复测量次数来压低重复性分量。5.2 GUM法对比计算与偏差讨论作为对照我再用GUM法评定同一个模型。偏导计算如下c_U 1 / s 0.25 c_U0 -1 / s -0.25 c_s -(U - U0) / s² -3.9999 / 16 ≈ -0.25 c_rep 1合成标准不确定度uc sqrt((0.25×0.00021)² (-0.25×0.00016)² (-0.25×0.0001)² (1×0.0002)²) ≈ 0.000212 MPa这与MCM的标准不确定度0.000212 MPa几乎一致。95%包含区间按正态分布和k2近似得到[0.99955, 1.00040] MPa与MCM结果差异在微乎其微的水平。这个算例说明模型线性度良好时GUM法的结果是可靠的MCM在这里起到了验证作用而非否定作用。但如果模型引入非线性项情况就不同了。我在原模型上加入一个二次修正项P (U - U0)/s β·(U-U0)² δ_rep其中β0.005 V⁻²。在ΔU≈4 V时二次项贡献约0.08 MPa这个量级已经不能忽略。此时GUM法的期望估计会漏掉E[(U-U0)²]的贡献因为一阶展开根本不包含二次项信息而MCM法通过直接抽样计算会把每一次抽样的二次项结果完整统计进去期望值和分散性都能忠实还原。实操中怎么判断差异是否显著我习惯的做法是如果MCM算出的标准不确定度与GUM法相差超过5%或者95%包含区间中心位置偏差超过扩展不确定度的10%就要重视模型非线性或分布非正态带来的影响。这时候报告里应该以MCM结果为准并追加分析和说明。5.3 结果可靠性检验不止看一个数字评定结果拿到手之后不要急着填报告。我会做两件事来交叉验证第一画直方图检查输出分布形态。如果直方图严重偏斜或明显双峰说明模型或输入量分布设置可能存在问题。正常情况下的压力传感器测量模型输出分布应该接近钟形或轻微偏斜。figure; histogram(P_sample, 200); xlabel(压力 (MPa)); ylabel(频数); title(输出量P的蒙特卡洛分布);第二跑两次不同分布假设下的结果对比。比如把U从正态改成均匀观察输出uP的变化幅度。如果变化量小于一个数量级说明该输入量的分布假设对结果不敏感报告中的分布假设站得住脚如果变化显著就必须向客户说明分布假设依据或者进一步做灵敏度分析找出影响最大的分布参数。我在实际项目里遇到过一种典型情况传感器灵敏度s的不确定度对结果影响巨大但校准证书上对s的分布类型语焉不详。这种情况下我宁可保守地按均匀分布处理得到更大的不确定度也不愿意为了好看的结果强行按正态分布。评定不确定度这件事宁可保守不可激进这是计量工作的基本原则。6. 实操中经常踩的坑与排查技巧6.1 常见问题速查表从随机数到内存的五个高频坑我把这几年用蒙特卡洛法评压力传感器测量模型时遇到的典型问题整理成了表方便大家对照排查。问题现象可能原因解决与排查方法每次运行结果差异大没有固定随机种子或M取值太小加rng固定种子将M从1e5增到1e6后观察稳定度出现“矩阵维度不匹配”报错模型传播环节用了/而不是./检查所有数组运算是否用了点运算符输出分布明显不对称但结果正常模型非线性项较强输入量分布非正态属正常现象应报告分位数区间而非k×u结果与GUM法差异超过10%模型存在强非线性项或输入量实际非正态以MCM结果为准分析差异来源并写入报告大M时内存不足或运行卡顿M取得过大同时开了多个大数组用单次抽样100万分批循环累加统计量不要一次生成全部样本关于最后一条多说一嘴。1e6样本时每个数组是8 MB左右几个数组加起来也就几十MB老电脑都能扛住。但如果你贪心直接上1e8一个数组就800 MB很容易卡死。遇到这种情况不要硬扛改成分批抽样每批1e6抽100批每批即时累加mean和std最后按统计公式合并即可。6.2 分布选型踩雷我把均匀分布和正态分布的结果差异放大过关于分布选型我想分享一个自己做过的对比实验。同一组输入量参数只是把U和U0从正态分布改成均匀分布最终标准不确定度uP差出了约30%。这个结果乍看很意外但细想就明白了均匀分布的标准差虽然相同但它没有正态分布的两端长尾抽样值集中在区间内导致输出量极值出现的概率变低整体分散性也就变小了。这给实际操作带来一个非常重要的启示分布选型不是拍脑袋的事情它要基于信息来源。如果某个不确定度分量来自校准证书的扩展不确定度证书通常按正态置信区间给出用正态分布是合理的如果来自仪器最大允差说明书只给误差限不给分布按均匀分布更严谨如果你为了省事全用正态分布评出来的不确定度有可能系统性偏大。我在报告中通常会给一句规范性描述“对数字万用表允许误差引入的不确定度分量采用均匀分布处理对校准证书引入的灵敏度不确定度分量按正态分布处理。”这句话写清楚之后客户和评审专家是挑不出毛病的。6.3 我对蒙特卡洛评定的一点整体体会做压力传感器不确定度评定这几年我最大的体会是蒙特卡洛法并不是“屠龙之技”它就是一个把测量模型每个细节都摆上台面的照妖镜。GUM法算得快但很多隐含假设你意识不到MCM法逼着你把输入量分布、相关性、样本量这些细节全部想清楚任何模糊处理都会在结果上留下痕迹。如果你是从零开始用这个方法我建议分三步走先拿本文的线性模型和代码跑通全流程对比GUM法找找感觉再把你的实际测量模型替换进来补全各分量的标准不确定度最后引入非线性项和混合分布观察结果变化并积累经验。等你处理过三五个完整案例后再回头看这个方法的原理会发现自己对测量不确定度的理解提升了一个层次。最后再分享一个小技巧在把评定结果写进校准报告时不要只给一个标准不确定度数字。把输入量清单、分布假设、抽样次数、包含概率一起列成表。一方面这符合评审要求另一方面将来有人对结果提出异议时你也能拿着完整的评定过程从容应对而不是空口争辩。
