基于齿面误差模型的成形磨齿机关键几何误差筛选方法
技术领域
[0001] 本
发明属于数控机床误差分析与
精度控制技术领域，具体的为一种基于齿面误差模型的成形磨齿机关键几何误差筛选方法。
背景技术
[0002] 大规格精密
齿轮广泛用于
风电、舰船、航空航天和重型机械等高端装备的动
力传递系统中，直接影响着装备的服役性能和可靠性。随着对装备的服役寿命、安全性和可靠性的要求日益增高，必须对制齿装备的精度性能进行控制和增强。五轴数控成形磨齿机是进
行大规格齿轮精加工的典型装备，由于生产效率高、易于齿面修形等优点得以广泛使用。但是该机床加工精度受多源误差协同影响，包括机床几何误差、热误差、力致
变形误差、伺服控制误差等。其中，机床几何误差是最重要的误差源，由机床各运动部件的制造及安装误差引起，对加工精度存在时不变、可重复的系统性影响。
[0003] 为高效率地提高成形磨齿机的加工精度，建立几何误差与加工误差间的定量映射关系模型至关重要。目前，国内外专家学者大多将刀具与
工件间的相对
位姿误差视为加工
误差，并致力于研究几何误差对该位姿误差的影响。然后考虑机床几何误差的随机特征和
耦合效应，利用敏感性分析方法，筛选出影响该加工误差的关键误差项。然而，由于实际加工中存在运动干涉等问题，对于复杂的齿轮
螺旋齿面加工而言，刀具与工件间的相对位姿
误差并不足以准确地描述和替代加工误差。显然，基于该位姿误差的关键几何误差筛选结
果也并不完全可靠。另外，即使误差模型已经成功建立，但如何准确有效地筛选出关键误差也是一大难题。
发明内容
[0004] 有鉴于此，本发明的目的在于提供一种基于齿面误差模型的成形磨齿机关键几何误差筛选方法，可以同时分析误差项个体效应和与其它误差项间的耦合效应对机床加工误
差的影响，从而识别出影响齿面精度的关键误差，为后续加工误差精确补偿提供可靠的理
论依据，以高效低耗地提高磨齿精度。
[0005] 为达到上述目的，本发明提供如下技术方案：
[0006] 一种基于齿面误差模型的成形磨齿机关键几何误差筛选方法，包括如下步骤：
[0007] 步骤一：基于机床
运动链和成形共轭磨削理论建立齿面误差模型
[0008] 步骤11：根据五轴数控成形磨齿机结构定义几何误差
[0009] 对五轴数控成形磨齿机进行结构分析和运动分析，得到影响磨齿精度的
位置无关几何误差和位置相关几何误差；
[0010] 步骤12：
砂轮坐标系和齿轮坐标系之间的位姿变换
[0011] 基于机床运动链，考虑几何误差，利用矩阵连乘法得到齿轮坐标系与砂轮坐标系间的位姿变换关系；
[0012] 步骤13：根据成形磨齿原理构建齿面误差模型
[0013] 根据实际成形共轭磨削理论，建立磨削
接触条件，利用二分法求解实际成形磨削接触线，并通过多条接触线拟合构建出齿面误差模型；
[0014] 步骤二：基于改进Sobol法的全局敏感性分析方法，筛选影响齿轮磨削精度的关键几何误差
[0015] 步骤21：根据齿面误差模型，采用一阶敏感指数表示单参数对模型输出方差的贡献率，采用高阶敏感指数表示多个参数耦合效应对模型输出方差的贡献率；同时考虑单参
数的个体效应和与其它参数的耦合效应，采用全局敏感指数表示该参数对模型输出方差的
综合贡献率；
[0016] 步骤22：采用改进Sobol法，根据几何误差测量数据确定取值区间并设计Sobol采样序列，基于Monte-Carlo估计值计算一阶敏感指数和全局敏感指数。
[0017] 进一步，所述步骤11中，五轴数控成形磨齿机包括X轴、Y轴、Z轴三根直线运动轴以及A轴、C轴两根旋转运动轴；
[0018] 由安装偏差引起的位置无关几何误差有11项，包括X轴、Y轴和Z轴之间的3项直线度误差和8项安装误差，分别为：
[0019] Z轴与X轴间的直线度误差SZX；
[0020] Y轴与X轴间的直线度误差SYX；
[0021] Y轴与Z轴间的直线度误差SYZ；
[0022] A轴在y方向上的安装位置误差δAy；
[0023] A轴在z方向上的安装位置误差δAz；
[0025] A轴绕z方向的安装姿态误差γAY；
[0026] C轴在x方向上的安装位置误差δCx；
[0027] C轴在y方向上的安装位置误差δCy；
[0028] C轴绕x方向的安装姿态误差αCY；
[0029] C轴绕y方向的安装姿态误差βCX；
[0030] 由制造
缺陷及运动磨损引起的位置相关几何误差有30项，每个运动轴包含6项误差，其中3项为位移误差，3项为
角度误差，分别为：
[0031] X轴在xyz方向上的位移误差分别为δx(X),δy(X)和δz(X)；
[0032] X轴在xyz方向上的角度误差分别为εx(X),εy(X)和εz(X)；
[0033] Y轴在xyz方向上的位移误差分别为δx(Y),δy(Y)和δz(Y)；
[0034] Y轴在xyz方向上的角度误差分别为εx(Y),εy(Y)和εz(Y)；
[0035] Z轴在xyz方向上的位移误差分别为δx(Z),δy(Z)和δz(Z)；
[0036] Z轴在xyz方向上的角度误差分别为εx(Z),εy(Z)和εz(Z)；
[0037] A轴在xyz方向上的位移误差分别为δx(A),δy(A)和δz(A)；
[0038] A轴在xyz方向上的角度误差分别为εx(A),εy(A)和εz(A)；
[0039] C轴在xyz方向上的位移误差分别为δx(C),δy(C)和δz(C)；
[0040] C轴在xyz方向上的角度误差分别为εx(C),εy(C)和εz(C)。
[0041] 进一步，所述步骤12中，五轴数控成形磨齿机的运动链为：齿轮G-C轴-底座R-X轴-Z轴-A轴-Y轴-砂轮W；以X,Y,Z,A,C分别表示X轴、Y轴、Z轴、A轴和C轴的数控指令，则各轴的运动变换矩阵分别为：
[0042]
[0043]
[0044] MCR＝MWY＝I
[0045] 理想情况下，利用矩阵连乘法得到齿轮坐标系与砂轮坐标系间的位姿变换关系为：
[0046]
[0047] 同理，几何误差引起的误差运动也可采用齐次变换矩阵表示，则由位置相关几何误差引起的误差矩阵可表示为：
[0048]
[0049] 其中，N表示N轴的数控指令； 表示N轴的位置相关几何误差引起的误差矩阵；此外，由位置无关几何误差引起的各轴误差矩阵为：
[0050]
[0051]
[0052] 由于存在几何误差，砂轮坐标系与齿轮坐标系之间的实际位姿变换关系为：
[0053]
[0054] 进一步，所述步骤13中，在成形磨齿加工中，忽略砂轮修整误差，则磨削砂轮轴向廓形与齿轮端面廓形完全相同；若以η表示砂轮轴向廓形参数，则砂轮轴向廓形的坐标矢量可表示为：
[0055] rq(η)＝[xq(η),0,yq(η),1]T
[0056] 在砂轮坐标系中，若将其绕砂轮轴旋转，形成的轨迹曲面就是砂轮回转曲面；以φ表示砂轮回转参数，则砂轮的坐标矢量为：
[0057] rw(η,φ)＝Mwq(φ)rq(η)，其中，
[0058] 在砂轮坐标系中，砂轮的单位法矢为：
[0059]
[0060] 根据步骤12得到的齿轮坐标系与砂轮坐标系间的理想位姿变换矩阵，可求得齿轮坐标系中理想砂轮为：
[0061]
[0062]
[0063] 考虑几何误差对砂轮的位姿影响，并以向量G表示所有几何误差时，可求得实际砂轮为：
[0064]
[0065]
[0066] 由于砂轮参数已知，齿轮与砂轮的共轭接触条件可表述为：从齿轮坐标系原点向砂轮回转面上的一点作径矢rg，如果这一点绕齿轮轴线kg作螺旋运动时的线速度矢量与其
在回转面上的法线ng垂直，则该点就是磨削接触点：
[0067]
[0068] 因此，理想接触条件和实际接触条件为：
[0069]
[0070]
[0071] 其中，γ为砂轮安装角；a为砂轮与齿轮中心距； 为
工作台旋转带动齿轮绕自身轴线旋转角度，相应的，自旋转的砂轮沿齿轮轴向向上移动 形成砂轮与齿轮间的相对
螺旋运动，其中p为螺旋参数；
[0072] 砂轮安装角γ、砂轮与齿轮中心距a均为常数，螺旋加工参数 仅在几何误差存在下才影响接触线的形状；同时机床几何误差向量G仅与机床运动轴位置有关，因此当 为常
数 时，接触条件f＝0只与砂轮轴向廓形参数η和砂轮回转廓形参数φ相关；
[0073] 由于η的范围已知，因此可将其等间距分成n个离散数值，令η＝ηj(j＝1,2,3,...,n)，然后根据接触条件利用二分法求得相应的φj，最后将(ηj,φj)代入rg，则可求得接触点的坐标向量 和单位法矢 由于φ是由η与接触条件联合求得，故φ可表示为η的函数，则由n个接触点拟合而成的第k条接触线可表示为：
[0074]
[0075]
[0076] 由于加工齿面可视作是由λ条接触线构成，则理想齿面和实际齿面可分别表示为：
[0077]
[0078]
[0079]
[0080]
[0081] 因此，齿面误差模型为：
[0082]
[0083]
[0084] 进一步，所述步骤21中，将齿面误差模型简化为一个多输入多输出的参数系统：
[0085] PE＝f(G)
[0086] m项几何误差作为输入，6项齿面位姿误差作为输出，其中，G＝[x1,x2,x3,…,xm]T表示几何误差向量，PE＝[δx,δy,δz,εx,εy,εz]T表示齿面位姿误差向量；
[0087] 由于几何误差输入是相互独立的且均匀分布于m维单位超空间，因此，以δx＝f(G)为例，f(G)可以分解为以下形式：
[0088]
[0089] 其中，f0是常数，fi(xi)表示参数xi的函数，fij(xi,xj)表示xi和xj的联合函数；
[0090] 由于各输入参数是相互独立的，各子函数也是相互
正交的，因此每个子函数的积分应等于0，也即：
[0091]
[0092] 基于条件期望，各子函数可定义为：
[0093] f0＝E(δx)
[0094] fi(xi)＝E(δx|xi)-f0
[0095] fij(xi,xj)＝E(δx|xi,xj)-f0-fi-fj
[0096] 其中，f0表示所有输入参数的总效应，fi(xi)表示xi的个体效应，fij(xi,xj)表示xi和xj的综合效应；同理，高阶函数项也可类似定义；由于输入参数相互独立，δx＝f(G)是平方可积的，则：
[0097]
[0098] 其中，Ωm表示m维单位参数空间；由于上式左边表示f(G)的方差，右边的各项表示各子函数的方差，因此方差分解表达式为：
[0099]
[0100] 其中，
[0101]
[0102]
[0103] G～i表示除xi外的所有输入参数，G～ij表示除xi和xj外的所有输入参数；
[0104] 为计算单参数对模型输出方差的贡献率，定义一阶敏感指数的计算表达式为：
[0105]
[0106] 同时，为衡量参数xi1,xi2,…,xip的耦合效应对模型输出方差的贡献率，定义高阶敏感指数的计算表达式为：
[0107]
[0108] 所有敏感指数必须满足以下隐含条件：
[0109]
[0110] 此外，同时考虑参数xi的个体效应和与其它参数的耦合效应，该参数对模型输出方差的综合贡献率可以用全局敏感指数表示：
[0111]
[0112] 其中，由于参数xi和xj的耦合效应在STi和STj中被重复计算，因此所有参数的全局敏感指数之和会大于1。
[0113] 进一步，所述步骤22中，基于改进Sobol法的敏感指数计算流程如下：
[0114] (1)根据输入参数的概率分布和Sobol序列，生成一个N×2m维的随机
采样矩阵HN×2m；
[0115] (2)根据HN×2m，构建输入矩阵AN×m,BN×m和ABiN×m；以H的前m列构造矩阵A，后m列构造矩阵B，同时生成衍生矩阵ABi，ABi中除第i列被矩阵B的第i列替换外，与矩阵A完全相等；
[0116] (3)将输入矩阵的任意一行作为模型的一组输入参数，则共有N(m+2)组输入数据；将每组输入数据分别代入模型δx＝f(G)中，运行模型后，可得相应的输出值f(A)，f(B)和f(ABi)；最后，利用Monte-Carlo估计值来计算一阶敏感指数和全局敏感指数：
[0117]
[0118] 本发明的有益效果在于：
[0119] 本发明的基于齿面误差模型的成形磨齿机关键几何误差筛选方法，首先，基于机床运动链和齐次变换矩阵，建立机床各部件坐标系间的位姿变换关系，得到砂轮与齿轮坐
标系间的带误差的位姿变换；然后，基于共轭磨削理论，考虑螺旋齿面成形磨削过程，建立实际齿面误差模型，以确切反映几何误差对齿面误差的映射关系；然后，根据几何误差测量数据确定取值区间，设计Sobol采样序列，通过齿面误差模型总方差的分解和Monte-Carlo
积分法计算每个几何误差项引起的模型方差及其占比，定量获得每项几何误差的各阶敏感
指数和全局敏感指数，该全局敏感性分析方法可以同时分析误差项个体效应和与其它误差
项间的耦合效应对机床加工误差的影响，从而识别出影响齿面精度的关键误差，为后续加
工误差精确补偿提供可靠的理论依据，以高效低耗地提高磨齿精度，具有采样样本数低、计算高效等优点。
附图说明
[0120] 为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：
[0121] 图1为五轴数控成形磨齿机的结构示意图；
[0122] 图2为利用五轴数控成形磨齿机成形磨齿加工的原理示意图；
[0123] 图3为δx误差方向上的敏感指数；
[0124] 图4为δy误差方向上的敏感指数；
[0125] 图5为δz误差方向上的敏感指数；
[0126] 图6为εx误差方向上的敏感指数；
[0127] 图7为εy误差方向上的敏感指数；
[0128] 图8为εz误差方向上的敏感指数。
具体实施方式
[0129] 下面结合附图和具体
实施例对本发明作进一步说明，以使本领域的技术人员可以更好的理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
[0130] 本实施例基于齿面误差模型的成形磨齿机关键几何误差筛选方法，包括如下步骤：
[0131] 步骤一：基于机床运动链和成形共轭磨削理论建立齿面误差模型
[0132] 根据齐次变换矩阵，建立机床相邻部件间的位姿变换关系；基于机床运动链，利用矩阵连乘法得到齿轮坐标系与砂轮坐标系间的位姿变换关系；再考虑实际成形共轭磨削理论，建立磨削接触条件，利用二分法求解实际成形磨削接触线，并通过多条接触线拟合构建出齿面误差模型。
[0133] 步骤11：根据五轴数控成形磨齿机结构定义几何误差
[0134] 对五轴数控成形磨齿机进行结构分析和运动分析，得到影响磨齿精度的位置无关几何误差和位置相关几何误差。
[0135] 具体的，如图1所示，五轴数控成形磨齿机包括X轴、Y轴、Z轴三根直线运动轴以及A轴、C轴两根旋转运动轴；成形磨齿机由于运动轴存在安装偏差和制造缺陷，多轴联动制齿时难免产生加工误差。一般而言，由安装偏差引起的位置无关几何误差(PIGE)有11项，包括三直
线轴间的3项直线度误差和两旋
转轴的8项安装误差，其中每个
旋转轴有2项安装位置误差和2项安装姿态误差，分别为：
[0136] Z轴与X轴间的直线度误差SZX；
[0137] Y轴与X轴间的直线度误差SYX；
[0138] Y轴与Z轴间的直线度误差SYZ；
[0139] A轴在y方向上的安装位置误差δAy；
[0140] A轴在z方向上的安装位置误差δAz；
[0141] A轴绕y方向的安装姿态误差βAZ；
[0142] A轴绕z方向的安装姿态误差γAY；
[0143] C轴在x方向上的安装位置误差δCx；
[0144] C轴在y方向上的安装位置误差δCy；
[0145] C轴绕x方向的安装姿态误差αCY；
[0146] C轴绕y方向的安装姿态误差βCX；
[0147] 由制造缺陷及运动磨损引起的位置相关几何误差(PDGE)有30项，每个运动轴包含6项误差，其中3项为位移误差，3项为角度误差，分别为：
[0148] X轴在xyz方向上的位移误差分别为δx(X),δy(X)和δz(X)；
[0149] X轴在xyz方向上的角度误差分别为εx(X),εy(X)和εz(X)；
[0150] Y轴在xyz方向上的位移误差分别为δx(Y),δy(Y)和δz(Y)；
[0151] Y轴在xyz方向上的角度误差分别为εx(Y),εy(Y)和εz(Y)；
[0152] Z轴在xyz方向上的位移误差分别为δx(Z),δy(Z)和δz(Z)；
[0153] Z轴在xyz方向上的角度误差分别为εx(Z),εy(Z)和εz(Z)；
[0154] A轴在xyz方向上的位移误差分别为δx(A),δy(A)和δz(A)；
[0155] A轴在xyz方向上的角度误差分别为εx(A),εy(A)和εz(A)；
[0156] C轴在xyz方向上的位移误差分别为δx(C),δy(C)和δz(C)；
[0157] C轴在xyz方向上的角度误差分别为εx(C),εy(C)和εz(C)。
[0158] 所有几何误差的定义和编号见表1。其中，Y轴主要用于刀具对中和砂轮修整，与成形磨齿过程相对独立，因此暂时忽略其几何误差对磨齿精度的影响，不进行编号。
[0159] 表1五轴数控成形磨齿机几何误差定义
[0160]
[0161] 步骤12：砂轮坐标系和齿轮坐标系之间的位姿变换
[0162] 基于机床运动链，利用矩阵连乘法得到齿轮坐标系与砂轮坐标系间的位姿变换关系。
[0163] 五轴数控成形磨齿机的运动链为：齿轮G-C轴-底座R-X轴-Z轴-A轴-Y轴-砂轮W；以X,Y,Z,A,C分别表示X轴、Y轴、Z轴、A轴和C轴的数控指令，则各轴的运动变换矩阵分别为：
[0164]
[0165]
[0166] MCR＝MWY＝I
[0167] 理想情况下，利用矩阵连乘法得到齿轮坐标系与砂轮坐标系间的位姿变换关系为：
[0168]
[0169] 同理，几何误差引起的误差运动也可采用齐次变换矩阵表示，则由位置相关几何误差引起的误差矩阵可表示为：
[0170]
[0171] 其中，N表示N轴的数控指令； 表示N轴的位置相关几何误差引起的误差矩阵；此外，由位置无关几何误差引起的各轴误差矩阵为：
[0172]
[0173]
[0174] 由于存在几何误差，砂轮坐标系与齿轮坐标系之间的实际位姿变换关系为：
[0175]
[0176] 步骤13：根据成形磨齿原理构建齿面误差模型
[0177] 根据实际成形共轭磨削理论，建立磨削接触条件，利用二分法求解实际成形磨削接触线，并通过多条接触线拟合构建出齿面误差模型。
[0178] 如图2，成形磨齿的加工原理为：首先，旋转A轴并移动X轴，使得砂轮安装角为γ，砂轮与齿轮中心距为a；然后，工作台旋转带动齿轮绕自身轴线旋转 同时自旋转的砂轮沿齿轮轴向向上移动 形成砂轮与齿轮间的相对螺旋运动，完成单个螺旋齿面的加工，其
中p为螺旋参数。成形磨齿是典型的线接触加工，接触线可视作加工齿面和磨削砂轮的基本组成单元。
[0179] 在成形磨齿加工中，忽略砂轮修整误差，则磨削砂轮轴向廓形与齿轮端面廓形完全相同；若以η表示砂轮轴向廓形参数，则砂轮轴向廓形的坐标矢量可表示为：
[0180] rq(η)＝[xq(η),0,yq(η),1]T
[0181] 在砂轮坐标系中，若将其绕砂轮轴旋转，形成的轨迹曲面就是砂轮回转曲面；以φ表示砂轮回转参数，则砂轮的坐标矢量为：
[0182] rw(η,φ)＝Mwq(φ)rq(η)，其中，
[0183] 在砂轮坐标系中，砂轮的单位法矢为：
[0184]
[0185] 根据步骤12得到的齿轮坐标系与砂轮坐标系间的理想位姿变换矩阵，可求得齿轮坐标系中理想砂轮为：
[0186]
[0187]
[0188] 考虑几何误差对砂轮的位姿影响，并以向量G表示所有几何误差时，可求得实际砂轮为：
[0189]
[0190]
[0191] 由于砂轮参数已知，齿轮与砂轮的共轭接触条件可表述为：从齿轮坐标系原点向砂轮回转面上的一点作径矢rg，如果这一点绕齿轮轴线kg作螺旋运动时的线速度矢量与其
g
在回转面上的法线n垂直，则该点就是磨削接触点：
[0192]
[0193] 因此，理想接触条件和实际接触条件为：
[0194]
[0195]
[0196] 其中，γ为砂轮安装角；a为砂轮与齿轮中心距； 为工作台旋转带动齿轮绕自身轴线旋转角度，相应的，自旋转的砂轮沿齿轮轴向向上移动 形成砂轮与齿轮间的相对
螺旋运动，其中p为螺旋参数；
[0197] 砂轮安装角γ、砂轮与齿轮中心距a均为常数，螺旋加工参数 仅在几何误差存在下才影响接触线的形状；同时机床几何误差向量G仅与机床运动轴位置有关，因此当 为常
数 时，接触条件f＝0只与砂轮轴向廓形参数η和砂轮回转廓形参数φ相关；
[0198] 由于η的范围已知，因此可将其等间距分成n个离散数值，令η＝ηj(j＝1,2,3,...,n)，然后根据接触条件利用二分法求得相应的φj，最后将(ηj,φj)代入rg，则可求得接触点的坐标向量 和单位法矢 由于φ是由η与接触条件联合求得，故φ可表示为η的函数，则由n个接触点拟合而成的第k条接触线可表示为：
[0199]
[0200]
[0201] 由于加工齿面可视作是由λ条接触线构成，则理想齿面和实际齿面可分别表示为：
[0202]
[0203]
[0204]
[0205]
[0206] 因此，齿面误差模型为：
[0207]
[0208]
[0209] 步骤二：基于改进Sobol法的全局敏感性分析方法，筛选影响齿轮磨削精度的关键几何误差
[0210] 改进Sobol法是一种基于方差分解的全局敏感性分析方法，适用于非线性高阶系统，具有采样要求低、计算效率高等优点。它通过计算参数及参数间耦合效应对模型输出方差的贡献率来定量评价参数敏感性。
[0211] 步骤21：根据齿面误差模型，采用一阶敏感指数表示单参数对模型输出方差的贡献率，采用高阶敏感指数表示多个参数耦合效应对模型输出方差的贡献率；同时考虑单参
数的个体效应和与其它参数的耦合效应，采用全局敏感指数表示该参数对模型输出方差的
综合贡献率。
[0212] 将齿面误差模型简化为一个多输入多输出的参数系统：
[0213] PE＝f(G)
[0214] m项几何误差作为输入，6项齿面位姿误差作为输出，其中，G＝[x1,x2,x3,…,xm]T表示几何误差向量，PE＝[δx,δy,δz,εx,εy,εz]T表示齿面位姿误差向量；
[0215] 由于模型输出是标量的，因此多输出的分析可以分解为多个独立的单输出分析过程；由于几何误差输入是相互独立的且均匀分布于m维单位超空间，因此，以δx＝f(G)为例，f(G)可以分解为以下形式：
[0216]
[0217] 其中，f0是常数，fi(xi)表示参数xi的函数，fij(xi,xj)表示xi和xj的联合函数；
[0218] 由于各输入参数是相互独立的，各子函数也是相互正交的，因此每个子函数的积分应等于0，也即：
[0219]
[0220] 基于条件期望，各子函数可定义为：
[0221] f0＝E(δx)
[0222] fi(xi)＝E(δx|xi)-f0
[0223] fij(xi,xj)＝E(δx|xi,xj)-f0-fi-fj
[0224] 其中，f0表示所有输入参数的总效应，fi(xi)表示xi的个体效应，fij(xi,xj)表示xi和xj的综合效应；同理，高阶函数项也可类似定义；由于输入参数相互独立，δx＝f(G)是平方可积的，则：
[0225]
[0226] 其中，Ωm表示m维单位参数空间；由于上式左边表示f(G)的方差，右边的各项表示各子函数的方差，因此方差分解表达式为：
[0227]
[0228] 其中，
[0229]
[0230]
[0231] G～i表示除xi外的所有输入参数，G～ij表示除xi和xj外的所有输入参数，其他高阶子函数的方差也可类似定义；
[0232] 为计算单参数对模型输出方差的贡献率，定义一阶敏感指数的计算表达式为：
[0233]
[0234] 同时，为衡量参数xi1,xi2,…,xip的耦合效应对模型输出方差的贡献率，定义高阶敏感指数的计算表达式为：
[0235]
[0236] 所有敏感指数必须满足以下隐含条件：
[0237]
[0238] 此外，同时考虑参数xi的个体效应和与其它参数的耦合效应，该参数对模型输出方差的综合贡献率可以用全局敏感指数表示：
[0239]
[0240] 其中，由于参数xi和xj的耦合效应在STi和STj中被重复计算，因此所有参数的全局敏感指数之和会大于1。
[0241] 步骤22：采用改进Sobol法，利用Monte-Carlo估计值计算一阶敏感指数和全局敏感指数。由于上述等式中含有大量多维积分，采用Monte-Carlo积分法进行计算。传统的
Sobol法直接采用计算机内部随机函数生成采样序列，会产生周期性重复。为克服这个问
题，改进的Sobol采用Sobol采样序列计算多维积分。该序列致力于在参数概率空间以巧妙
的随机方法生成均匀分布的采样值。
[0242] 基于改进Sobol法的敏感指数计算流程如下：
[0243] (1)根据输入参数的概率分布和Sobol序列，生成一个N×2m维的随机采样矩阵HN×2m。
[0244] (2)根据HN×2m，构建输入矩阵AN×m,BN×m和ABiN×m。以H的前m列构造矩阵A，后m列构造矩阵B，同时生成衍生矩阵ABi，ABi中除第i列被矩阵B的第i列替换外，与矩阵A完全相等；
[0245] (3)将输入矩阵的任意一行作为模型的一组输入参数，则共有N(m+2)组输入数据；将每组输入数据分别代入模型δx＝f(G)中，运行模型后，可得相应的输出值f(A)，f(B)和f(ABi)；最后，利用Monte-Carlo估计值来计算一阶敏感指数和全局敏感指数：
[0246]
[0247] 下面结合成形磨齿机的关键几何误差筛选仿真算例对本发明做进一步说明。
[0248] 步骤一：基于机床运动链和成形共轭磨削理论建立齿面误差模型
[0249] 五轴数控成形磨齿机的运动链为：G-C-R-X-Z-A-Y-W。齿轮坐标系与砂轮坐标系间的位姿变换矩阵可以计算为：
[0250]
[0251] 其中， 其余的齐次变换矩阵为：
[0252]
[0253]
[0254] 而由于几何误差的存在，实际的位姿变换矩阵为：
[0255]
[0256] 其中，
[0257]
[0258]
[0259]
[0260]
[0261]
[0262] 同理，根据共轭磨削理论：
[0263] f＝(kg×rg+pkg)·ng＝0
[0264] 利用二分法，可以求得由n个接触点拟合而成的第k条接触线可表示为：
[0265]
[0266]
[0267] 由于加工齿面可视作是由λ条接触线构成，则理想齿面和实际齿面可分别表示为：
[0268]
[0269]
[0270]
[0271]
[0272] 因此，齿面误差模型为：
[0273]
[0274]
[0275] 步骤二：基于改进Sobol法的全局敏感性分析方法
[0276] 改进Sobol法是一种基于方差分解的全局敏感性分析方法，适用于非线性高阶系统，具有采样要求低、计算效率高等优点。它通过计算参数及参数间耦合效应对模型输出方差的贡献率来定量评价参数敏感性。
[0277] 现将齿面误差模型简化为一个多输入多输出的参数系统，m项几何误差作为输入，6项齿面位姿误差作为输出。由于输入几何误差的取值范围显著影响敏感性分析结果，因
此，首先应确定各项几何误差的概率分布范围。对于该五轴数控成形磨齿机而言，测定结果表明，位移误差分布于[0,20]μm，角度误差分布
[0278] 于[0,29]mdeg。以砂轮中心对应的某刀位点而言，以x＝709.7784mm，y＝-17.7358mm，z＝50mm为例，几何误差的Sobol随机采样序列的统计值如表2。
[0279] 表2输入几何误差的Sobol采样序列统计值
[0280]
[0281] 考虑计算耗时和收敛速率，将砂轮中心的磨削轨迹离散化为11个刀位点，也即11条磨削接触线，构建最终的加工齿面。因此，共计11×50×(33+2)次Monte-Carlo仿真将被
执行，以计算对于全齿面的误差敏感指数。
[0282] 因此，基于齿面误差模型，利用改进Sobol法进行全局敏感性分析，在6个误差方向上的各项几何误差的一阶敏感指数和全局敏感指数计算结果如图3-图8所示。
[0283] 根据图3-图8，可以看出虽然Si和STi两种敏感指数的分布具有一定的相似性，但由于全局敏感指数还考虑了几何误差之间的耦合效应，它反映几何误差对齿面误差的影响更为全面有效。因此，将STi作为关键误差的评价指标。表3列出了在每个误差方向上筛选出的关键几何误差。
[0284] 表3各误差方向上的关键几何误差
[0285]
[0286] 根据图3-8所示，在33项几何误差中，共筛选出13项关键误差，包括δx(X)，εy(X)，εz(X)，δx(Z)，εz(Z)，δx(A)，εy(A)，βAZ，δx(C)，εy(C)，εz(C)，δCx和βCX。因此，在后续的关键误差精确补偿中，这13项误差应该重点关注进行补偿消减。若需要着重提高某误差方向上的加工精度，则相应的多项关键几何误差应当进行补偿消减。但同时必须注意，单方向的关键误差补偿，可能导致其它方向上的加工精度降低。因此，必须对各方向上的精度进行必要的平衡，以防止总体加工精度下降。
[0287] 以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。
本技术领域的技术人员在本发明
基础上所作的等同替代或变换，均在本发明
的保护范围之内。本发明的保护范围以
权利要求书为准。
