CN103115601A - 轴类零件圆柱度的超差测定方法 - Google Patents

Abstract

一种轴类零件圆柱度的超差测定方法，该方法首先用三坐标测量机对轴的圆柱面测量四次并分别获取测点坐标，建立了最小区域圆柱度误差评定模型；用粒子群算法分别搜索四次测量轴与坐标平面xoy的交点坐标、轴线方向向量的实际值及最小区域圆柱度误差；构建参数矩阵并求其协方差阵，获取交点坐标及轴线方向向量的不确定度、相关不确定度并计算单个测点测量值的不确定度；执行自适应蒙特卡洛算法获得圆柱度误差不确定度值及其在设定的置信概率下的包含区间。本发明中蒙特卡洛算法次数不断增加，直至所需要的各种结果达到统计意义上的稳定，能够同时计算最小区域圆柱度误差、圆柱度误差不确定度及其包含区间，准确测定圆柱度超差的轴类零件。

Description

轴类零件圆柱度的超差测定方法

技术领域

本发明涉及一种圆柱度的超差测定方法，尤其涉及一种轴类零件圆柱度的超差测定方法，属于精密计量与计算机应用领域。 

背景技术

形状误差是评定机械零件的重要指标，其大小对产品质量及其使用寿命至关重要，在生产中要加以测量和控制。评定形状误差有多种方法，以圆柱度误差评定为例，评定方法有最小区域法、最小二乘法、最大内切法和最小外接法，但各种方法得出的结果都不相同，甚至差异很大，导致产品出现误收或误废，直接影响产品的质量和成本，因此国际标准ISO/1101和国家标准GB/T1958-2004都规定，形状误差值用包容实际被测要素且具有最小宽度E或最小直径φE的包容区域来表示（简称最小区域法），并以此为仲裁方法。以最小区域法评定形状误差，能够在不改变硬件设备的前提下，提高测量设备的检测精度。 

一个科学完整的测量结果，除了应给出被测量的最佳估计值外，还应同时给出测量结果的不确定度。不确定度必须正确评价，否则会导致工件的误收和误废。因此如何快速精确评定不确定度成为测定零件是否合格的关键。为此，国家计量技术规范JJF1059-1999《测量不确定度评定与表示》（简称GUM）规定了测量不确定度的评定与表示的通用规则，它对科学研究、工程技术、以及商贸中大量存在的测量结果的处理和表示，具有指导意义。由于基于GUM的测量不确定度评定方法容易受到直接测量量相关性的限制，特别是对于非线性模型，在使用GUM计算测量不确定度时计算过程中存在诸多近似，导致计算的精度降低；加之其严格的解析操作在实践中有时缺乏可操作性，给测量不确定度的评定带来不便。针对GUM中测量不确定度评定方法存在的不足，2008年国际标准化组织正式颁布了ISO/IEC导则98-32008（GUM）及其一系列补充标准，使不确定度的应用更加科学合理。 

轴类零件是机械产品的重要组成部分，其精度的高低对产品的质量及其使用寿命至关重要，衡量轴类零件形状误差大小的指标有轴线直线度、轴剖面素线直 线度、横截面圆度及轴的圆柱度误差，因为圆柱度误差能够同时反映轴的横截面圆度、轴剖面素线直线度和轴线直线度误差，因此被广泛应用于轴类零件形状误差评定。测量圆柱度误差设备有圆柱度仪、三坐标测量机（CMM）等，虽然圆柱度仪测量精度高，但因价格昂贵，对测量环境要求高而使其应用受到一定限制。目前在实验室、工厂常用CMM测量圆柱度误差，使用CMM测量圆柱度误差时得到的是一系列离散测点值，需要经过数据处理求解圆柱度误差，目前的三坐标测量机只是给出最小二乘拟合的圆柱度误差，既没有给出最小区域法的圆柱度误差（最小区域圆柱度误差），更没有对测量结果的不确定度进行评价。虽然近年来一直有学者致力于最小区域圆柱度误差求解研究，提出了多种计算方法，但是没有对最小区域圆柱度误差的测量不确定度进行评价；仅有极少数学者研究了采用GUM方法通过对圆柱度误差最小二乘评定模型求一阶导数近似计算圆柱度误差不确定度，由于最小二乘法本身提供的仅是形状误差的近似评价结果，并不保证解的最小区域性，按最小二乘法计算的结果比最小区域法求得结果大1.8%～30%，平均过估计为10%。因而不能准确测定轴类零件是否超差。 

综合上述分析，当前对相关领域研究工作存在的不足主要是：缺乏能够对轴类零件圆柱度是否超差进行准确测定的方法。 

发明内容

本发明提供一种能够提高误差识别精度的轴类零件圆柱度的超差测定方法。 

本发明的技术方案为：一种轴类零件圆柱度的超差测定方法，包括如下步骤： 

步骤1以测量平台上任意一点o为原点，做三条互相垂直的数轴x轴，y轴和z轴，建立测量空间直角坐标系oxyz，坐标平面xoy位于测量平台上，将被测轴置于测量空间直角坐标系oxyz中且被测轴的轴线L与oz轴平行，轴线L的理想方向向量为(p,q,1)，p,q和1分别为轴线L沿x，y和z方向的理想方向向量，轴线L与坐标平面xoy的理想交点为O′(a,b,0)，a,b和0分别为理想交点O′在测量空间直角坐标系oxyz下的坐标值，轴线L可以表示为， 

$\frac{x-a}{p}=\frac{y-b}{q}=\frac{z}{1},$

步骤2令j=1，j为测量序数 

步骤3使用三坐标测量机测得被测轴的圆柱面的测点P_ij(x_ij,y_ij,z_ij)，P_ij为 第j次测量的第i个测点，i为测点序号，i=1,2,…,n,n为测点数目且n为正整数，x_ij，y_ij和z_ij分别为测点P_ij在测量空间直角坐标系oxyz下的坐标值， 

步骤4建立最小区域圆柱度误差的目标函数，所述目标函数为： 

E_j=f(a_j,b_j,p_j,q_j)=min(max(R_ij)-min(R_ij))=min(R_maxj-R_minj) 

其中， $R_{\mathrm{ij}}=\sqrt{{[x_{\mathrm{ij}}-(p_j{z}_{\mathrm{ij}}+a_j)]}^2+{[y_{\mathrm{ij}}-(q_j{z}_{\mathrm{ij}}+b_j)]}^2}$

式中a_j,b_j分别表示第j次测量的轴线L与坐标平面xoy的交点的x坐标值、y坐标值，p_j,q_j分别表示第j次测量的轴线L沿x和y方向的方向向量，R_ij为第j次测量的测点P_ij到轴线L的距离，E_j为第j次测量得到被测轴的最小区域圆柱度误差，R_maxj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最大值，所述最大值对应的测点坐标为（x_maxj，y_maxj，z_maxj），R_minj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最小值，所述最小值对应的测点坐标为（x_minj，y_minj，z_minj）， 

步骤5使用粒子群算法求解步骤4所述的目标函数，获得第j次测量的轴线L的a_j,b_j,p_j,q_j值及最小区域圆柱度误差E_j，如果j>4，则转入步骤7，否则，进入步骤6， 

步骤6令j=j+1，返回步骤3， 

步骤7计算E_j的平均值获得被测轴的最小区域圆柱度误差

分别计算a_j,b_j,p_j及q_j的平均值

及

由步骤5得到的四组轴线L的a_j,b_j,p_j,q_j值构建参数矩阵Ir， 

$\mathrm{Ir}=\left(\begin{array}{cccc}{a}_1& b_1& p_1& q_1\\ a_2& b_2& p_2& q_2\\ a_3& b_3& p_3& q_3\\ a_4& b_4& p_4& q_4\end{array}\right),$

步骤8求参数矩阵Ir的协方差阵V，获取轴线L的a,b,p,q估计值的不确定度u_a,u_b,u_p,u_q及相关不确定度u_ab,u_ap,u_aq,u_ba,u_bp,u_bq,u_pa,u_pb,u_pq,u_qa,u_qb,u_qp, 

$V=\mathrm{cov}\left(\mathrm{Ir}\right)=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right]$

其中，cov(Ir)表示对参数矩阵Ir求协方差阵，u_a,u_b,u_p,u_q分别为轴线L的a,b,p,q估计值的不确定度，u_ab,u_ap,u_aq分别为a估计值与b估计值,p估计值,q估计值之间的相关不确定度，u_ba,u_bp,u_bq分别为b估计值与a估计值,p估计值,q估计值之间的相关不确定度，u_pa,u_pb,u_pq分别为p估计值与a估计值,b估计值,q估计值之间的相关不确定度，u_qa,u_qb,u_qp分别为q估计值与a估计值,b估计值,p估计值之间的相关不确定度，其中， 

u_ab=u_ba,u_ap=u_pa,u_aq=u_qa,u_bp=u_pb,u_pq=u_qp， 

步骤9根据三坐标测量机精度、测量条件、环境因素计算单个测点测量值的不确定度u₀， 

步骤10执行自适应蒙特卡洛算法，求解圆柱度误差不确定度及其设定置信概率下的包含区间， 

步骤10.1令指针h=1， 

步骤10.2执行蒙特卡洛算法，得到圆柱度误差

r为蒙特卡洛算法所产生模拟量组数的序号，r=1,2,...,M,M为所产生模拟量的组数，M=10⁴，并将 

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.2.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方

为方差，产生服从正态分布的随机数，

表示均值为*、方差为

的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布  $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量 

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为μ、协方差阵为V的多变量正态分布，其中， 

$\mathrm{\μ}=\left(\begin{array}{c}\stackrel{\‾}{a}\\ \stackrel{\‾}{b}\\ \stackrel{\‾}{p}\\ \stackrel{\‾}{q}\end{array}\right),$ $V=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right],$

再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.2.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差并将得到的10⁴个圆柱度误差

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3令h=h+1，执行蒙特卡洛算法，得到圆柱度误差 

并将

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方

为方差，产生服从正态分布

的随机数，

表示均值为*、方差为

的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布  $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量 

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为 _μ、协方差阵为V的多变量正态分布，再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.3.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差

并将得到的10⁴个圆柱度误差

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.4计算排序后的圆柱度误差的平均值

标准不确定度 

100d%包含区间的左端点

右端点

d为置信概率，d=0.95， 

${\stackrel{\‾}{e}}^{\left(h\right)}=\frac{1}{{10}^4}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{e}_r^{\left(h\right)}$

$u^2\left({\stackrel{\‾}{e}}^{\left(h\right)}\right)=\frac{1}{{10}^4-1}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{(e_r^{\left(h\right)}-{\stackrel{\‾}{e}}^{\left(h\right)})}^2$

${\mathrm{\β}}_{\mathrm{left}}^{\left(h\right)}=e_{250}^{\left(h\right)},$ ${\mathrm{\β}}_{\mathrm{right}}^{\left(h\right)}=e_{9750}^{\left(h\right)},$

步骤10.5计算

平均值的标准偏差计算 

平均值

的标准偏差

$\stackrel{\‾}{e}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\stackrel{\‾}{e}}^{\left(g\right)},g=\mathrm{1,2},...,h,$ g为圆柱度误差平均值 

的序号， 

$s_{\stackrel{\‾}{e}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\stackrel{\‾}{e}}^{\left(g\right)}-\stackrel{\‾}{e})}^2}$

$\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}u\left({\stackrel{\‾}{e}}^{\left(g\right)}\right)$

$s_{\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{(u\left({\stackrel{\‾}{u}}^{\left(g\right)}\right)-\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right))}^2}$

步骤10.6计算

平均值

的标准偏差

以及计算 

平均值的标准偏差

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}})}^2}$

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}})}^2}$

步骤10.7根据单个测点测量值的不确定度u₀计算数值容差δ， 

将u₀表示为C×10^l的形式，C为两位有效十进制整数，l是负整数，l=-6,-5,-4或-3，数值容差δ取为， 

$\mathrm{\δ}=\frac{1}{2}{10}^l$

步骤10.8如果

及

中的任何一个值都不大于δ/5，则进入步骤11，否则，返回步骤10.3， 

步骤11以标准不确定度的平均值

为圆柱度误差不确定度u(e)，以100d%包含区间的左端点平均值右端点平均值

构造包含区间 

当被测轴最小区域圆柱度误差Em大于100d%包含区间右端点

则被测轴为圆柱度超差的轴类零件。 

本发明的有益效果在于： 

建立了圆柱度误差最小区域评定模型，采用粒子群算法求解最小区域圆柱度误差，极大地提高了圆柱度误差计算精度和效率；采用自适应蒙特卡洛算法计算轴类零件圆柱度误差不确定度及设定的置信概率下的包含区间，其中蒙特卡洛算法次数不断增加，直至所需要的各种结果达到统计意义上的稳定，克服了对于非线性模型在使用GUM计算测量不确定度时计算过程中存在诸多近似，导致计算的 精度降低和蒙特卡洛方法需要事先设定蒙特卡洛模拟实验次数，不能保证结果是否稳定的不足。该方法能够同时计算圆柱度误差，圆柱度误差不确定度和设定的置信概率下的包含区间，不仅算法可靠，优化效率高，而且能够准确测定被测轴的圆柱度是否超差。 

附图说明

图1轴类零件圆柱面测量图。 

图2最小区域圆柱度误差示意图。 

图3本发明的流程图。 

图4用粒子群算法搜索最小区域圆柱度误差进化过程图。 

图5自适应蒙特卡洛算法流程图。 

图6基于自适应蒙特卡洛算法的圆柱度误差频率分布图。 

具体实施方式

一种轴类零件圆柱度的超差测定方法，其步骤如下： 

步骤1以测量平台上任意一点o为原点，做三条互相垂直的数轴x轴，y轴和z轴，建立测量空间直角坐标系oxyz，坐标平面xoy位于测量平台上，将被测轴置于测量空间直角坐标系oxyz中且被测轴的轴线L与oz轴平行，轴线L的理想方向向量为(p,q,1)，p,q和1分别为轴线L沿x，y和z方向的理想方向向量，轴线L与坐标平面xoy的理想交点为O′(a,b,0)，a,b和0分别为理想交点O′在测量空间直角坐标系oxyz下的坐标值，轴线L可以表示为， 

$\frac{x-a}{p}=\frac{y-b}{q}=\frac{z}{1}$

步骤2令j=1，j为测量序数 

步骤3使用三坐标测量机测得被测轴的圆柱面的测点P_ij(x_ij,y_ij,z_ij)，P_ij为第j次测量的第i个测点，i为测点序号，i=1,2,…,n,n为测点数目且n为正整数，n取等于或大于60的整数，x_ij，y_ij和z_ij分别为测点P_ij在测量空间直角坐标系oxyz下的坐标值， 

步骤4建立最小区域圆柱度误差的目标函数，所述目标函数为： 

E_j=f(a_j,b_j,p_j,q_j)=min(max(R_ij)-min(R_ij))=min(R_maxj-R_minj) 

其中， 

$R_{\mathrm{ij}}=\sqrt{{[x_{\mathrm{ij}}-(p_j{z}_{\mathrm{ij}}+a_j)]}^2+{[y_{\mathrm{ij}}-(q_j{z}_{\mathrm{ij}}+b_j)]}^2}$

式中a_j,b_j分别表示第j次测量的轴线L与坐标平面xoy的交点的x坐标值、y坐标值，p_j,q_j分别表示第j次测量的轴线L沿x和y方向的方向向量，R_ij为第j次测量的测点P_ij到轴线L的距离，E_j为第j次测量得到被测轴的最小区域圆柱度误差，R_maxj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最大值，所述最大值对应的测点坐标为（x_maxj，y_maxj，z_maxj），R_minj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最小值，所述最小值对应的测点坐标为（x_minj，y_minj，z_minj）， 

步骤5使用粒子群算法求解步骤4所述的目标函数，获得第j次测量的轴线L的a_j,b_j,p_j,q_j值及最小区域圆柱度误差E_j，如果j>4，则转入步骤7，否则，进入步骤6， 

步骤5.1随机产生粒子的初始位置和初始速度 

选择粒子大小popsize为20的种群，以1×4维的实数向量为种群中的第k个粒子的位置pos_k，k=1,2,…,20，第k个粒子的位置表示为pos_k=(A_1k,A_2k,A_3k,A_4k)，其中A_1k,A_2k,A_3k,A_4k分别为对应a_j,b_j,p_j,q_j的可能取值，以另一1×4维的实数向量为种群中的第k个粒子的速度，表示为v_k=(B_1k,B_2k,B_3k,B_4k)，其中B_1k,B_2k,B_3k,B_4k分别为对应粒子在a_j,b_j,p_j,q_j上的飞行速度， 

在

[-0.05，0.05]和[-0.05，0.05]数值区域内随机产生20个粒子的a_j,b_j,p_j,q_j，

为第j次测量所有测点x_ij坐标的平均值，

为第j次测量所有测点y_ij坐标的平均值，以产生的A_1k,A_2k,A_3k,A_4k作为第k个粒子初始位置

为第k个粒子在第t代的位置，令t=1，第k个粒子初始位置

进入粒子迭代，并根据随机产生的粒子初始位置

计算粒子初始位置的目标函数值选取初始位置目标函数值最小的粒子的 位置作为第一代全局最佳粒子位置gbest^t，t=1；第k个粒子初始位置作为第k个粒子第一代的局部最佳粒子位置进入粒子迭代，t=1，k=1,2,…,20，在[-0.05，0.05]数值区域内随机产生20个粒子的B_1k,B_2k,B_3k,B_4k作为初始速度第k个粒子初始速度

进入粒子迭代， 

第k个粒子至第t代以前搜索到的最优位置称为粒子k第t代的局部最佳粒子位置

整个粒子群至第t代以前搜索到的最优位置称为第t代的全局最佳粒子位置gbest^t， 

步骤5.2采用浓缩因子法修改粒子速度 

第k个粒子在迭代的第t代采用如下浓缩因子法修改速度： 

$v_k^{t+1}=K(v_k^t+c_1{\mathrm{rand}}_1^t({\mathrm{pbest}}_k^t-{\mathrm{pos}}_k^t)+c_2{\mathrm{rand}}_2^t({\mathrm{gbest}}^t-{\mathrm{pos}}_k^t))$

式中

分别为第k个粒子在第t代的速度和位置,

和分别为在第t代随机产生的1×4维向量，向量中的每一元素在[0,1]区间随机产生，c₁，c₂为加速因子，分别决定第k个粒子向局部最佳粒子

和全局最佳粒子gbest^t方向飞行的相对拉力，K为浓缩因子，c₁，c₂满足

为加速因子的和，加速因子c₁，c₂和浓缩因子K取值分别为2.05，2.05和0.73， 

步骤5.3用步骤5.2得到修改后的速度

改变粒子位置 

在迭代的第t代，将第k个粒子位置修改为： 

${\mathrm{pos}}_k^{t+1}={\mathrm{pos}}_k^t+v_k^{t+1}\mathrm{\Δt}$

△t是时间步长，设置为1， 

步骤5.4计算粒子位置改变后的所有粒子目标函数值

计算第k个粒子位置改变为

后的粒子目标函数值

步骤5.5更新局部最佳粒子位置

如果位置改变后第k个粒子的目标^函数值

小于未改变前第k个粒子局部最佳位置的目标函数值

则用

更新第k个粒子的第t代的局部最佳粒子位置

作为第k个粒子的第t+1代的局部最佳粒子

位置，否则，第k个粒子的局部最佳粒子位置

作为第t+1代的局部最佳粒子位置 

步骤5.6更新全局最佳粒子位置 

找出位置改变后所有粒子

目标函数值

最小的粒子mpos，如果粒子mpos的目标函数值E_j(mpos)小于未改变前全局最佳粒子位置的目标函数值E_j(gbest^t)，则用mpos更新全局最佳粒子位置gbest^t，作为第t+1代的全局最佳粒子位置gbest^t+1，否则，第t代的全局最佳粒子位置gbest^t作为第t+1代的全局最佳粒子位置gbest^t+1， 

步骤5.7令t=t+1，如果t=201，则进入步骤6，否则，重复步骤5.3～5.6， 

当算法达到终止条件时，全局最佳粒子位置gbest²⁰⁰对应参数a_j,b_j,p_j,q_j的最优值，全局最佳粒子位置gbest²⁰⁰的目标函数值E_j(gbest²⁰⁰)即为搜索到的圆柱度误差最小区域解E_j， 

步骤6令j=j+1，返回步骤3， 

步骤7计算E_j的平均值获得被测轴的最小区域圆柱度误差

分别计算a_j,b_j,p_j及q_j的平均值

及

由步骤5得到的四组轴线L的a_j,b_j,p_j,q_j值构建参数矩阵Ir， 

$\mathrm{Ir}=\left(\begin{array}{cccc}{a}_1& b_1& p_1& q_1\\ a_2& b_2& p_2& q_2\\ a_3& b_3& p_3& q_3\\ a_4& b_4& p_4& q_4\end{array}\right),$

步骤8求参数矩阵Ir的协方差阵V，获取轴线L的a,b,p,q估计值的不确定度u_a,u_b,u_p,u_q及相关不确定度u_ab,u_ap,u_aq,u_ba,u_bp,u_bq,u_pa,u_pb,u_pq,u_qa,u_qb,u_qp,其中u_a,u_b,u_p,u_q分别为轴线L的a,b,p,q估计值的不确定度，u_ab,u_ap,u_aq分别为a估计值与b估计值，p估计值,q估计值之间的相关不确定度，u_ba,u_bp,u_bq分别为b估计值与a估计值,p估计值,q估计值之间的相关不确定度，u_pa,u_pb,u_pq分别为p估计值与a估计值,b估计值,q估计值之间的相关不确定度，u_qa,u_qb,u_qp分别为q估计值与a估计值,b估计值,p估计值之间的相关不确定度，利用Matlab语言提供的函数cov()对参数矩阵Ir求协方差阵V， 

$V=\mathrm{cov}\left(\mathrm{Ir}\right)=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right]$

其中，u_ab=u_ba,u_ap=u_pa,u_aq=u_qa,u_bp=u_pb,u_pq=u_qp， 

步骤9根据三坐标测量机精度、测量条件、环境因素计算单个测点测量值的不确定度u₀， 

采用三坐标测量机对轴的圆柱度误差测量时，影响单个测点测量值的不确定度的主要因素有： 

重复性引起的不确定度u_re，温度引起的不确定度u_T，CMM的漂移和迟滞引起的不确定度u_pc， 

上述不确定度因素彼此独立，单个测点测量值的不确定度为， 

$u_0=\sqrt{u_{\mathrm{re}}^2+u_T^2+u_{\mathrm{pc}}^2}$

步骤10执行自适应蒙特卡洛算法，求解圆柱度误差不确定度及设定置信概率下的包含区间， 

步骤10.1令指针h=1， 

步骤10.2执行蒙特卡洛算法，得到圆柱度误差

r为蒙特卡洛算法所产生模拟量组数的序号，r=1,2,...,10⁴，并将

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差  $e_r^{\left(1\right)},r=1,2,...,{10}^4,$

步骤10.2.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方

为方差，产生服从正态分布的随机数，

表示均值为*、方差为的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布  $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量 

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为μ、协方差阵为V的多变量正态分布，其中， 

$\mathrm{\μ}=\left(\begin{array}{c}\stackrel{\‾}{a}\\ \stackrel{\‾}{b}\\ \stackrel{\‾}{p}\\ \stackrel{\‾}{q}\end{array}\right),$ $V=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right],$

再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.2.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差

并将得到的

采用Matlab语言提供的函数sort()实现从小到大排序，且相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3令h=h+1，执行蒙特卡洛算法，得到圆柱度误差 

并将

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方

为方差，产生 服从正态分布

的随机数，

表示均值为*、方差为

的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布  $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量 

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为μ、协方差阵为V的多变量正态分布，再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.3.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差

并将得到的

采用Matlab语言提供的函数sort()排序实现从小到大排序，且相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.4计算排序后的圆柱度误差

的平均值标准不确定度 

100d%包含区间的左端点

右端点

${\stackrel{\‾}{e}}^{\left(h\right)}=\frac{1}{{10}^4}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{e}_r^{\left(h\right)}$

$u^2\left({\stackrel{\‾}{e}}^{\left(h\right)}\right)=\frac{1}{{10}^4-1}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{(e_r^{\left(h\right)}-{\stackrel{\‾}{e}}^{\left(h\right)})}^2$

${\mathrm{\β}}_{\mathrm{left}}^{\left(h\right)}=e_{250}^{\left(h\right)},$ ${\mathrm{\β}}_{\mathrm{right}}^{\left(h\right)}=e_{9750}^{\left(h\right)},$

步骤10.5计算

平均值的标准偏差

计算 

平均值

的标准偏差

$\stackrel{\‾}{e}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\stackrel{\‾}{e}}^{\left(g\right)},g=\mathrm{1,2},...,h,$ g为圆柱度误差平均值 

的序号， 

$s_{\stackrel{\‾}{e}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\stackrel{\‾}{e}}^{\left(g\right)}-\stackrel{\‾}{e})}^2}$

$\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}u\left({\stackrel{\‾}{e}}^{\left(g\right)}\right)$

$s_{\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{(u\left({\stackrel{\‾}{u}}^{\left(g\right)}\right)-\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right))}^2}$

步骤10.6计算

平均值

的标准偏差

以及计算 

平均值

的标准偏差

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}})}^2}$

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}})}^2}$

步骤10.7根据单个测点测量值的不确定度u₀计算数值容差δ， 

将u₀表示为C×10^l的形式，C为两位有效十进制整数，l是负整数，l=-6,-5,-4或-3，数值容差δ取为， 

$\mathrm{\δ}=\frac{1}{2}{10}^l$

步骤10.8如果

及

中的任何一个值都不大于δ/5，则进入步骤11，否则，返回步骤10.3， 

步骤11以标准不确定度的平均值

为圆柱度误差不确定度u(e)，以100d%包含区间的左端点平均值

右端点平均值

构造包含区间 

当被测轴最小区域圆柱度误差Em大于100d%包含区间右端点

则被测轴为圆柱度超差的轴类零件。 

以下结合附图对本发明做进一步的说明： 

1、在圆柱面上获取测点P_ij(x_ij,y_ij,z_ij)(i=1,2,…,n,n为测点数目，n=60，测量序数j=1,2,3,4)，见附图1。 

为了证实本方法的正确性，采用Miracle NC454三坐标测量机对一直径36.2mm，长90mm，公差为0.012mm轴的圆柱面实测四次，获得四组测量数据。 

2、初始化算法参数 

加速因子c₁，c₂和浓缩因子K取值分别为2.05，2.05和0.73。 

3、随机产生粒子的初始位置和初始速度 

在用粒子群算法求解圆柱度误差时，粒子的种群大小popsize取值为20，目标函数值E_j的大小由a_j,b_j,p_j,q_j四个参数决定，粒子k的初始位置pos_k=(A_1k,A_2k,A_3k,A_4k)的四个分量A_1k,A_2k,A_3k,A_4k分别在 

[-0.05，0.05]和[-0.05，0.05]区间上随机产生；粒子k的初始速度v_k=(B_1k,B_2k,B_3k,B_4k)的四个分量B_1k,B_2k,B_3k,B_4k均在[-0.05,0.05]上随机产生。 

4、根据已建立最小区域圆柱度误差评定模型及测量数据，计算粒子初始位置的目标函数值，选取粒子初始位置目标函数值最小的粒子作为第一代全局最佳粒子位置gbest^t，t=1；第k（k=1,2,…,popsize）个粒子初始位置作为粒子k第一代局部最佳粒子位置

进入粒子迭代，t=1。 

5、采用浓缩因子法修改粒子速度、改变粒子的位置 

粒子k在迭代的第t代采用如下浓缩因子法修改速度： 

$v_k^{t+1}=K(v_k^t+c_1{\mathrm{rand}}_1^t({\mathrm{pbest}}_k^t-{\mathrm{pos}}_k^t)+c_2{\mathrm{rand}}_2^t({\mathrm{gbest}}^t-{\mathrm{pos}}_k^t))$

粒子k的位置根据来改变，△t取值为1。 

6、计算改变位置后所有粒子的目标函数值 

由步骤4所述公式计算位置改变后所有粒子的目标函数值

7、更新局部最佳粒子位置 

如果位置改变后粒子k的目标函数值

小于未改变前该粒子局部最佳位置的目标函数值则用

更新局部最佳粒子位置

8、更新全局最佳粒子位置 

如果位置改变后粒子k的目标函数值

小于未修改前全局最佳粒子位置的目标函数值E_j(gbest^t)，则用

更新全局最佳粒子位置gbest^t。 

9、进化代数t大于200时，算法终止 

用上述算法搜索该实例四组测量数据最小区域圆柱度误差的优化过程见附图4，由图可见，在约140代时就已搜索到最小区域圆柱度误差，a_j,b_j,p_j,q_j的优化结果见表1，由表1求得被测轴的最小区域圆柱度误差

表1参数及圆柱度误差的优化结果 

由表1可见对用粒子群算法计算得到的最小区域圆柱度误差明显小于最小二乘法得到的圆柱度误差，提高了圆柱度误差评定精度；另外对同一圆柱面进行四 次测量求出轴线的a_j,b_j,p_j,q_j值及最小区域圆柱度误差明显不同，说明圆柱度误差测量过程中确实存在不确定度，对圆柱度误差不确定度测定是非常必要的。 

10、根据三坐标测量机精度、测量条件、环境因素计算单个测点测量值的不确定度值u₀

采用Miracle NC454三坐标测量机对轴实测时，影响单个测点测量值的不确定度的因素有： 

(1)重复性引起的不确定度u_re：在Miracle NC454三坐标测量机上对某一点重复测量30次，求30次测量值的平均值及标准差，得到重复性引起的不确定度u_re为1.48μm。 

(2)温度引起的不确定度u_T： 

测量过程中经观察温度偏离标准温度(20℃)的最大值为1℃，工件的温度膨胀系数为1.1μm/(100mm×℃)，由温度引起的不确定度极限值α_T为： 

该不确定服从均匀分布，由温度引起的不确定度u_T: 

$u_T=0.40/\sqrt{3}=0.23\mathrm{\μm}$

(3)CMM的漂移和迟滞引起的不确定度u_pc

由CMM的漂移和迟滞引起的不确定度为0.10μm。 

上述不确定度因素彼此独立，单点测量的不确定度为 

$u_0=\sqrt{u_{\mathrm{re}}^2+u_T^2+u_{\mathrm{pc}}^2}=\sqrt{{1.48}^2+{0.23}^2+{0.10}^2}=1.50\mathrm{\μm}=0.0015\mathrm{mm}$

11、执行自适应蒙特卡洛算法，求解圆柱度误差不确定度值及包含区间11.1指针h=1，首次执行蒙特卡洛算法，得到圆柱度误差E_r ⁽¹⁾,r=1,2,...,10⁴,用三坐标测量机测量圆柱面求解圆柱度误差E_j时，其目标函数值由a_j,b_j,p_j,q_j四个输入量通过步骤4所述函数关系f来确定，以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj，最小值对应的测点坐标为x_minj，y_minj，z_minj，采用Matlab语言提供的函数randn(ξ,η)产生ξ×η维均值为0，方差为1的标准正态分布阵列，用Matlab语言编写X1=(x_maxj,y_maxj,z_maxj,x_minj,y_minj,z_minj)′+u₀*randn(6,10⁴) 产生服从正态分布

的6×10⁴维模拟阵列X1，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，阵列X1的第一行到第六行分别表示服从正态分布  $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 的随机数，阵列X1的每一列表示一组蒙特卡洛算法的模拟量

由表1求得参数矩阵： 

$\mathrm{Ir}=\left(\begin{array}{cccc}-0.0333& 0.0126& 0.0048& -0.0020\\ -0.0345& 0.0147& 0.0047& -0.0020\\ -0.0346& 0.0155& 0.0048& -0.0020\\ -0.0327& 0.0131& 0.0047& -0.0020\end{array}\right)$

求参数矩阵Ir的协方差阵V： 

$V=\mathrm{cov}\left(\mathrm{Ir}\right)=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right]={10}^{-5}*\left[\begin{array}{cccc}0.0862& -0.1126& -0.0012& 0\\ -0.1126& 0.1836& 0.0005& 0\\ -0.0012& 0.0005& 0.0003& 0\\ 0& 0& 0& 0\end{array}\right]$

由协方差阵可见a,b,p,q之间相关，需要采用多变量正态分布N(μ,V)产生a,b,p,q的模拟量

其中， 

$\mathrm{\μ}=\left(\begin{array}{c}\stackrel{\‾}{a}\\ \stackrel{\‾}{b}\\ \stackrel{\‾}{p}\\ \stackrel{\‾}{q}\end{array}\right)=\left(\begin{array}{c}-0.0338\\ 0.0140\\ 0.0047\\ -0.0020\end{array}\right),$ ${V=10}^{-5}*\left[\begin{array}{cccc}0.0862& -0.1126& -0.0012& 0\\ -0.1126& 0.1836& 0.0005& 0\\ -0.0012& 0.0005& 0.0003& 0\\ 0& 0& 0& 0\end{array}\right]$

多变量正态分布随机数由下述方法产生： 

(1)对矩阵V进行乔里斯基（Cholesky）因子分解V=RR^T，得到下三角阵R 

$R=\left[\begin{array}{cccc}0.0009& 0& 0& 0\\ -0.0012& 0.0006& 0& 0\\ -0.0000& -0.0000& 0.0001& 0\\ 0& 0& 0& 0\end{array}\right],$

(2)由randn(4,10⁴)函数产生4×10⁴维标准正态阵列Z， 

(3)由X2=μ1^T+R^TZ生成4×10⁴维正态阵列X2，其中上角标T表示对矩阵求转置，1表示1×10⁴维的全1矩阵，X2中每一列表示一组蒙特卡洛算法的模拟量

$X=\left[\begin{array}{c}{X}_1\\ X_2\end{array}\right]$ 组成10×10⁴维阵列，每一列表示一组蒙特卡洛算法的模拟量 

由所产生的10⁴组模拟量代入步骤4所述目标函数计算公式得到10⁴个圆柱度误差

将得到

采用Matlab语言提供的函数sort()排序，得到排序后圆柱度误差

11.2令h=h+1，执行蒙特卡洛算法，得到圆柱度误差

并将得到的

采用Matlab语言提供的函数sort()排序，得到排序后圆柱度误差 

11.3计算排序后的圆柱度误差

的平均值

标准不确定度

95%包含区间的左端点右端点

11.4计算

平均值

的标准偏差

计算平均值 的标准偏差

计算

平均值

的标准偏差

以及计算 

平均值

的标准偏差

11.5将单个测点测量值的不确定度u₀表示为u₀=15×10^-4，则数值容差δ为： 

$\mathrm{\δ}=\frac{1}{2}{10}^{-4}$

即δ等于0.00005， 

11.6当h=156时，

及

中的任何一个值都不大于δ/5=0.00001，则得到圆柱度误差的不确定度为0.0045mm及95%置信概率的包含区间[0.0096,0.0185]，由156×10⁴组模拟量获得的圆柱度误差频率分布见图6。 

12.因被测轴的圆柱度公差为0.012mm，而求得的最小区域圆柱度误差为0.0146mm，误差大于公差，如果不考虑测量不确定度及置信概率下的包含区间，则认为该轴超差，判断该轴为不合格零件。事实上测量过程中存在不确定度，圆柱度误差0.0146落在95%置信概率包含区间[0.0096,0.0185]中，不能判定该零件不合格，需由供需双方协商决定该零件是否合格，只有当被测轴最小区域圆柱度 误差大于包含区间右端点0.0185mm则认为该轴圆柱度超差，判断该轴为不合格零件。由此可见能否准确测定圆柱度是否超差将直接关系到产品的质量和性能。 

Claims (1)

1.一种轴类零件圆柱度的超差测定方法，其特征在于，具体步骤如下：

步骤1以测量平台上任意一点o为原点，做三条互相垂直的数轴x轴，y轴和z轴，建立测量空间直角坐标系oxyz，坐标平面xoy位于测量平台上，将被测轴置于测量空间直角坐标系oxyz中且被测轴的轴线L与oz轴平行，轴线L的理想方向向量为(p,q,1)，p,q和1分别为轴线L沿x，y和z方向的理想方向向量，轴线L与坐标平面xoy的理想交点为O′(a,b,0)，a,b和0分别为理想交点O′在测量空间直角坐标系oxyz下的坐标值，轴线L可以表示为，

$\frac{x-a}{p}=\frac{y-b}{q}=\frac{z}{1},$

步骤2令j=1，j为测量序数

步骤3使用三坐标测量机测得被测轴的圆柱面的测点P_ij(x_ij,y_ij,z_ij)，P_ij为第j次测量的第i个测点，i为测点序号，i=1,2,…,n,n为测点数目且n为正整数，x_ij，y_ij和z_ij分别为测点P_ij在测量空间直角坐标系oxyz下的坐标值，

步骤4建立最小区域圆柱度误差的目标函数，所述目标函数为：

E_j=f(a_j,b_j,p_j,q_j)=min(max(R_ij)-min(R_ij))=min(R_maxj-R_minj)

其中， $R_{\mathrm{ij}}=\sqrt{{[x_{\mathrm{ij}}-(p_j{z}_{\mathrm{ij}}+a_j)]}^2+{[y_{\mathrm{ij}}-(q_j{z}_{\mathrm{ij}}+b_j)]}^2}$

式中a_j,b_j分别表示第j次测量的轴线L与坐标平面xoy的交点

的x坐标值、y坐标值，p_j,q_j分别表示第j次测量的轴线L沿x和y方向的方向向量，R_ij为第j次测量的测点P_ij到轴线L的距离，E_j为第j次测量得到被测轴的最小区域圆柱度误差，R_maxj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最大值，所述最大值对应的测点坐标为（x_maxj，y_maxj，z_maxj），R_minj为第j次测量的n个测点P_ij到轴线L的距离R_ij中的最小值，所述最小值对应的测点坐标为（x_minj，y_minj，z_minj），

步骤5使用粒子群算法求解步骤4所述的目标函数，获得第j次测量的轴线L的a_j,b_j,p_j,q_j值及最小区域圆柱度误差E_j，如果j>4，则转入步骤7，否则，进入步骤6，

步骤6令j=j+1，返回步骤3，

步骤7计算E_j的平均值获得被测轴的最小区域圆柱度误差

分别计算a_j,b_j,p_j及q_j的平均值及

由步骤5得到的四组轴线L的a_j,b_j,p_j,q_j值构建参数矩阵Ir，

$\mathrm{Ir}=\left(\begin{array}{cccc}{a}_1& b_1& p_1& q_1\\ a_2& b_2& p_2& q_2\\ a_3& b_3& p_3& q_3\\ a_4& b_4& p_4& q_4\end{array}\right),$

步骤8求参数矩阵Ir的协方差阵V，获取轴线L的a,b,p,q估计值的不确定度u_a,u_b,u_p,u_q及相关不确定度u_ab,u_ap,u_aq,u_ba,u_bp,u_bq,u_pa,u_pb,u_pq,u_qa,u_qb,u_qp,

$V=\mathrm{cov}\left(\mathrm{Ir}\right)=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right]$

其中，cov(Ir)表示对参数矩阵Ir求协方差阵，u_a,u_b,u_p,u_q分别为轴线L的a,b,p,q估计值的不确定度，u_ab,u_ap,u_aq分别为a估计值与b估计值,p估计值,q估计值之间的相关不确定度，u_ba,u_bp,u_bq分别为b估计值与a估计值,p估计值,q估计值之间的相关不确定度，u_pa,u_pb,u_pq分别为p估计值与a估计值,b估计值,q估计值之间的相关不确定度，u_qa,u_qb,u_qp分别为q估计值与a估计值,b估计值,p估计值之间的相关不确定度，其中，

u_ab=u_ba,u_ap=u_pa,u_aq=u_qa,u_bp=u_pb,u_pq=u_qp，

步骤9根据三坐标测量机精度、测量条件、环境因素计算单个测点测量值的不确定度u₀，

步骤10执行自适应蒙特卡洛算法，求解圆柱度误差不确定度及其设定置信概率下的包含区间，

步骤10.1令指针h=1，

步骤10.2执行蒙特卡洛算法，得到圆柱度误差

r为蒙特卡洛算法所产生模拟量组数的序号，r=1,2,...,M,M为所产生模拟量的组数，M=10⁴，并将按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.2.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方为方差，产生服从正态分布

的随机数，

表示均值为*、方差为

的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布 $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为μ、协方差阵为V的多变量正态分布，其中，

$\mathrm{\μ}=\left(\begin{array}{c}\stackrel{\‾}{a}\\ \stackrel{\‾}{b}\\ \stackrel{\‾}{p}\\ \stackrel{\‾}{q}\end{array}\right),$ $V=\left[\begin{array}{cccc}{u}_a^2& u_{\mathrm{ab}}& u_{\mathrm{ap}}& u_{\mathrm{aq}}\\ u_{\mathrm{ba}}& u_b^2& u_{\mathrm{bp}}& u_{\mathrm{bq}}\\ u_{\mathrm{pa}}& u_{\mathrm{pb}}& u_p^2& u_{\mathrm{pq}}\\ u_{\mathrm{qa}}& u_{\mathrm{qb}}& u_{\mathrm{qp}}& u_q^2\end{array}\right],$

再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.2.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差

并将得到的10⁴个圆柱度误差

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3令h=h+1，执行蒙特卡洛算法，得到圆柱度误差

并将

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.3.1选取任意一次测量获得的测点坐标，分别以步骤4得到的最大值对应的测点坐标x_maxj，y_maxj，z_maxj、步骤4得到的最小值对应的测点坐标x_minj，y_minj，z_minj为均值，以步骤9得到的单个测点测量值的不确定度u₀的平方为方差，产生服从正态分布

的随机数，

表示均值为*、方差为

的正态分布，*分别为x_maxj，y_maxj，z_maxj，x_minj，y_minj和z_minj，并从服从正态分布 $N(x_{\max j},u_0^2),N(y_{\max j},u_0^2),N(z_{\max j},u_0^2),N(x_{\min j},u_0^2),N(y_{\min j},u_0^2),N(z_{\min j},u_0^2)$ 中产生6×10⁴维模拟阵列X1，每一列表示一组蒙特卡洛算法的模拟量

产生服从多变量正态分布N(μ,V)的随机向量，N(μ,V)表示期望值向量为μ、协方差阵为V的多变量正态分布，再从服从多变量正态分布N(μ,V)中产生4×10⁴维模拟阵列X2，每一列表示一组蒙特卡洛算法的模拟量

步骤10.3.2由所产生的10⁴组模拟量

由步骤4所述的目标函数计算得到10⁴个圆柱度误差

并将得到的10⁴个圆柱度误差

按照从小到大顺序排序，如果几个圆柱度误差相等，则相等的圆柱度误差任意排序，得到排序后的圆柱度误差

步骤10.4计算排序后的圆柱度误差

的平均值

标准不确定度

100d%包含区间的左端点

右端点

d为置信概率，d=0.95，

${\stackrel{\‾}{e}}^{\left(h\right)}=\frac{1}{{10}^4}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{e}_r^{\left(h\right)}$

$u^2\left({\stackrel{\‾}{e}}^{\left(h\right)}\right)=\frac{1}{{10}^4-1}\underset{r=1}{\overset{{10}^4}{\mathrm{\Σ}}}{(e_r^{\left(h\right)}-{\stackrel{\‾}{e}}^{\left(h\right)})}^2$

${\mathrm{\β}}_{\mathrm{left}}^{\left(h\right)}=e_{250}^{\left(h\right)},$ ${\mathrm{\β}}_{\mathrm{right}}^{\left(h\right)}=e_{9750}^{\left(h\right)},$

步骤10.5计算

平均值

的标准偏差

计算

平均值

的标准偏差

$\stackrel{\‾}{e}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\stackrel{\‾}{e}}^{\left(g\right)},g=\mathrm{1,2},...,h,$ g为圆柱度误差平均值的序号，

$s_{\stackrel{\‾}{e}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\stackrel{\‾}{e}}^{\left(g\right)}-\stackrel{\‾}{e})}^2}$

$\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}u\left({\stackrel{\‾}{e}}^{\left(g\right)}\right)$

$s_{\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right)}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{(u\left({\stackrel{\‾}{u}}^{\left(g\right)}\right)-\stackrel{\‾}{u}\left(\stackrel{\‾}{e}\right))}^2}$

步骤10.6计算

平均值

的标准偏差

以及计算

平均值

的标准偏差

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{left}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{left}})}^2}$

${\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}=\frac{1}{h}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}$

$s_{{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}}}=\sqrt{\frac{1}{h(h-1)}\underset{g=1}{\overset{h}{\mathrm{\Σ}}}{({\mathrm{\β}}_{\mathrm{right}}^{\left(g\right)}-{\stackrel{\‾}{\mathrm{\β}}}_{\mathrm{right}})}^2}$

步骤10.7根据单个测点测量值的不确定度u₀计算数值容差δ，

将u₀表示为C×10^l的形式，C为两位有效十进制整数，l是负整数，l=-6,-5,-4或-3，数值容差δ取为，

$\mathrm{\δ}=\frac{1}{2}{10}^l$

步骤10.8如果

及中的任何一个值都不大于δ/5，则进入步骤11，否则，返回步骤10.3，

步骤11以标准不确定度的平均值

为圆柱度误差不确定度u(e)，以100d%包含区间的左端点平均值右端点平均值

构造包含区间

当被测轴最小区域圆柱度误差Em大于100d%包含区间右端点则被测轴为圆柱度超差的轴类零件。

Images