滤棒检测中圆度误差评定方法的选择

'滤棒检测中圆度误差评定方法的选择邓春宁（福建龙岩卷烟厂龙岩364000）摘要本文对滤棒圆度误差评定进行了研究，重点研究了圆度误差评定的最小二乘方法，指出只有满足“小误差假设”和“小偏差假设”，最小二乘通用算法的评定结果才是严格意义上最小二乘解。针对这一结论，提出了一种由最小二乘迭代算法得到的精确的最小二乘解。关键词滤棒检测圆度误差误差评定方法0引言形状误差测量就是将被测形状与理想形状进行比较，从而确定并用数值描述实际形状与理想形状的差异。每一个形状误差参数的确定过程包括测量和评定两个阶段。测量阶段是根据被测形状误差项的定义，获得统一坐标系下采样面、采样线和采样点上的轮廓坐标值；评定阶段则是根据定义对统一的坐标值进行处理，求得具体的形状误差值。圆度误差作为一种形状误差是卷烟滤棒生产过程中一个重要的检测指标，其误差评定可以采用多种方法，实质是一个非线性最优化问题。1圆度误差的几何特性根据国际标准ISO1101和国标GB1850的定义，圆度公差带是指在半径差为公差t的两同心圆之间的区域。因此，圆度误差的量值大小反映在圆周的半径方向上，即误差具有径向性。圆度误差的第二个几何特性是圆度误差变化具有“周期性”。因为圆形零件横截面的实际轮廓形状是一个复杂的封闭曲线轮廓，轮廓上各点径向误差的大小不同，而且在圆周的一圈上以2π为周期连续变化。圆度误差的周期性可以用富氏级数来表示，在极坐标系中：∞∞∞ρ(θ)=r+acosiθ+bsiniθ=r+csin(iθ+a)（1）0∑1∑i0∑iii=1i=0i=122式中，ρ(θ)为极径，r为富氏级数常量，a,b为富氏级数的系数，c=a+b，011iiia=arctgab。iii式（1）的实际意义可以看成由一个平均半径为r的圆和若干个周期不同的形状圆误差0波形迭加而成。i=1，级数展开中的acosθ+bsinθ项决定了平均半径圆，偏心量是1122e=c=a+b;i=2，级数展开式中的acos2θ+bsin2θ项反映在极坐标系中是一11122 个椭圆；i=3，级数展开式中acos3θ+bsin3θ项反映在极坐标系中是一个三棱圆。依33次类推，n次谐波，反映在极坐标系中是一个n边棱圆。根据国标GB1850定义，偏心影响、表面粗糙度和表面波纹度（高次谐波分量）的影响均应从ρ(θ)中剔除，所以圆度误差函数nΔr(θ)=∑(acosiθ+bsiniθ)（2）i1i=2通常n取值为7~9。由圆度误差的几何特性可知，圆度误差的评定关键在于如何根据最小条件或最小条件的近似来确定理想圆的圆心位置。2圆度误差的评定方法2.1圆度误差评定方法种类圆度误差评定主要有4种方法：①最小包容区域法（亦称最小半径法）以包含实际轮廓、且半径差为最小的两同心圆的圆心为理想圆心，这样确定的理想圆心满足最小条件。最小条件的判定法则是：两个同心包容圆应与实际圆的轮廓分别至少有两个切点，4个切点相间分布。②最小外接圆法最小外接圆法用实际轮廓的最小外接圆圆心来近似理想圆圆心。最小外接圆的判定法则为：在外接圆上至少有3点与外接圆相切，且组成锐角三角形，或者只存在两点与外接圆相切，且两点的连线过圆心。这两个法则称为锐角三角形法则和对径法则。③最大内接圆法最大内接圆法用实际轮廓的最大内接圆圆心来近似理想圆圆心，最大内接圆的判定法则与最小外接圆相似。④最小二乘法最小二乘法用实际轮廓的最小二乘圆圆心近似理想圆圆心。上述4种方法只有最小包容区域法严格符合最小条件，形状误差评定结果精度最高。但是最小包容区域中心的求解过程实质是一个复杂的、带约束条件的非线性最优化问题，这类问题的求解通常采用两种方法：一是在一定的假设前提下，利用线性评定模型并结合单纯形法、置换法求解；二是线性逼近、迭代求解。前者存在模型上的误差，评定精度受到限制；后者求解过程复杂，计算量大，而且还存在解的收敛性和稳定性问题，应用受到一定的限制。最小二乘法虽然是一种近似方法，但其评定结果优于最小外接圆和最大内接圆法，非常接近于最小包容区域法（评定误差通常小于5%），广泛应用于工程实践。最小外接圆法和最大内接圆法相当于与轴、孔实际装配时的情况，具有实际工程意义。通常最小外接圆法用于评定轴类零件的圆度误差，最大内接圆法用于评定孔类零件的圆度误差。2.1最小二乘法给定任意的圆度测量数据及测量点的直角坐标(xy)(i=1,2,?N)，可以得到如下误差ii方程22(x−u)+(y−u)=R+e(i=1,2,??N)（3）i1i2i 运用最小二乘原理，令准则函数N22Ja=∑(x1−u1)+(y1−u2)−R（4）i=1使J→min的(,)uu及R为最小二乘圆的圆心及半径。由于e为(,)uu及R的非线性函a12i12数，不易直接求解，故将准则函数J改为如下准则函数：aN2222Jb=∑{}(x1−u1)+(y1−u2)−R（5）i=1222令u=R−u−u，则式(5)可写为如下形式：312N222Jb=∑[]2xiu1+2yiu2+u3−(xi1+yi)（6）i=1这样参数（u,u）和R的估计问题可由线性最小二乘法直接求解。这种算法称为圆度误差12评定的最小二乘通用算法，被人们普遍采用。由于J有别于J，它所得到的理想圆不是严ba格的最小二乘圆。由式（3）、（5）可得N222222Jb=∑[](xi−u1)+(yi−u2)+R.(xi−u1)+(yi−u2)−Ri=1N2222=∑(2R+ei)(xi−u1)+(yi−u2)−Ri=12Ne22122=4R∑1+.(xi−u1)+(yi−u2)−R（7）i−12r2ei⎛ei⎞考虑到ei与R相比是微量，1+≈⎜1+⎟，Jb等价于如下准则函数：2R⎝2R⎠Ne2i22Jc=∑1+.(xi−u1)+(yi−u2)−R（8）i=12R比较式（4）和（8）可知，通用算法可以看成是对J进行加权平均处理后的结果，相应的权a 值为1+eR(i=1,2,?N)。eR越小，基于J和J的计算结果越接近。因此，当e与R11abi相比是微量时，即满足“小误差假设”，通用算法是适用的。圆度误差测量数据通常采用极坐标表示，此时误差方程可表示为22R=ucosa+usina+R−(usina−ucosa)+Σi1i2i1i2ii(i=1,2,??N)（9）式中，Σ(i=1,2,??N)为残余误差。同样，为了求得圆心坐标u,u及半径R，可以极i12N2小化如下准则函数：Jb=∑εii=1N222=∑[R1−u1cosai−u2sinai−R−(u1sinai−u2cosai)]（10）i=1显然，ε是u、u及R的非线性函数，不易直接求解。为此，定义如下准则函数：i12NN22222Jpdi==−−∑∑εi2(Ruisincaiua2osi).+εεii（11）式中，ii==11222p=2R−(usina−ucosa)+ε（12）111211由式（9）可得N222Je=∑[]R1−u1cosai−u2sinai+R−(u1sinai−u2cosai)i=1222.R−ucosa−usina−R−(usina−ucosa)i1i2i1i2iN22222=∑(Ri+u1+u2−R−2Ricosaiu1−2Risinaiu2)（13）i=1令222K=R−u−u（14）12则有N22Jb=∑(Ri−K−2Ricosaiu1−2Risinaiu2)（15）i=1 由式（15）可以进行线性最小乘求解，确定理想圆参数。采用准则函数J的最小二乘e算法称为极坐标下的圆度误差评定的最小二乘通用算法。由式（11）可以看出，J对J进行了加权处理（其权值为p），因此，通过J而求edie得的结果与通过J求得的结果两者之间存在一定的偏差，只有当uu、、ε较小时，即满d12i足“小偏差假设”和“小误差假设”，权值p才基本相等，两种算法圆度误差评定结果基本i一致。由以上分析可知，只有满足“小误差假设”和“小偏差假设”时，选用准则函数J和Jbe的通用算法所得到的圆度误差评定结果才与严格意义上最小二乘结果基本一致。当以上两个假设不成立时，通用算法结果可能会出现较大偏差。由于准则函数J（J）是非线性的，无法直接求解，本文提出一种迭代方法来逼近最ba小二乘解。将极坐标下的准则函数J改写为以下形式：bn12JN=∑(Ri−Rˆ)（16）2Ni=1式中，22Rˆ=ucosa+usina+R−(usina−ucosa)（17）i1i2i1i2iT定义θ=[]u,uR为最小二乘圆的圆心坐标和半径，则由式（16）可得到递推公式1212(k+1)J(θ)=kJ(θ)+[R−Rˆ(θ)]（18）k+1kk+1k+12不失一般性，令θˆ表示极小化J而得到的圆心和半径的估计值，在θˆ处作一阶泰勒展开来kkk逼近Rˆ(θ)，得k+1TRˆ(θ)≈Rˆ(θ)+Φ(θ−θˆ)（19）k+1k+1k+1kdRˆk+1式中，Φ=，k+1dθθ=θˆk将式（19）代入式（18），并对θ求导可得 T(k+1)J(θ)≈kJ(θ)−Φ{R−Rˆ(θ)−Φ(θ−θˆ)}（20）k+1kk+1k+1k+1k+1k又将J(θ)在θˆ处作一价泰勒展开，得kkj(θ)=j(θˆ)+j(θˆ)(θ−θˆ)（21）kkkkk将式（21）代入式(20)，并注意θˆ使J(θ)达到极小值，即J(θˆ)=0，所以kkkT(k+1)J(θ)≈kJ(θˆ)(θ−θˆ)−Φ{R−Rˆ(θˆ)−Φ(θ−θˆ)}（22）k+1kkkk+1k+1k+1kk+1k在式（22）中令θ=θˆ，则J(θ)=0，整理上式可得k+1k+1−1Tθˆ≈θˆ+[]kj(θˆ)+ΦΦΦ[R−Rˆ(θˆ)]（23）k+1kkkk+1k+1k+1K+1K+1K上式即为参数估计的递推形式，为了便于计算，令−1Tp=[]kj(θˆ)+ΦΦ（24）k+1kkk+1k+1对式（20）求导可得T(k+1)j=kj(θˆ)+ΦΦ（25）k+1kkk+1k+1由式（24）、（25）用矩阵求逆，可得TpΦpkk+1kp=p−（26）k+1kT1+ΦpΦk+1kk+1由此，可归纳出最小二乘圆参数估计的递推计算公式：⎧Σ=R−Rˆ(θˆ)kkkk−1⎪T⎪pk−1ΦkΦkpk−1⎨pk=pk−1−T（27）1+ΦpΦ⎪kk−1k⎪θˆ=θˆ+pΦΣ⎩kk−1kkk由以上递推公式，一次递推过程往往不能得到最小二乘圆的精确参数，需要经过多次迭代才能得到符合要求的参数。由以上分析，可以得到最小二乘圆求解递推算法步骤如下：①取θˆ为合适的初值，p=102~104I(I为3阶单位矩阵)；00②k=k+1；③由式（27）递推计算ε、p、θˆ；kkk④转步骤②，直到计算θˆ，重复步骤②~④，直到θˆ收敛。0N 2.3最小外接圆法和最大内接圆法最小外接圆法和最大内接圆法的算法过程和判别法则基本相同，因此，下面只给出最小外接圆法的算法实现。为提高算法效率，对圆度误差测量数据作以下预处理：首先采用最小二乘通用算法确定理想圆的粗略圆心0(x，y)和半径R；然后以圆心0(x，y)为极点，重新计算各采样点00000的极径R和极角a，并将极径偏差R−R存入数组r(n)，对应的极角数组a(n)也作相应iii0交换。最小外接圆法的基本算法过程：首先按锐角三角形法则搜索3点组，要求3点组所确定的外接圆半径最小，且包容圆的实际轮廓，如果搜索结束仍然不能找到满足条件的3点组，则按对径法则搜索两点组，要求以两点连线为直径的外接圆半径最小，且包容圆的实际轮廓。为了缩短平均搜索时间，循环选点应从极径偏差数组r(n)的大端开始（最大内接圆法则从小端开始选点），并且设置常数C，循环选点只选择极径偏差大于r(c)的采样点（极径偏差小于r(c)的点不可能在最小外接圆上，因此不必选择）。3实验结果为了检验圆度误差评定算法性能，设计了如下仿真试验。仿真数据用极坐标表示，由以下公式产生：9r(θ)=r0+c1sin(θ+a1)+∑(aicosiθ+bisiniθ)（28）i=2式中，r为理想圆的半径，c,a为圆心偏移量和偏移方向，理论圆度误差为Δr(θ)的极大011值减去极小值。Δr(θ)为9Δr(θ)=∑(aicosiθ+bisiniθ)（29）i=2式（29）中c、b(i=2,?,9)确定了圆度误差的形式和大小，可以随机产生，也可以根据具11体应用任务设计。如香烟滤棒圆截面变形主要是椭圆变形，设计仿真数据时，应满足a>>a、b>>b(i≠2)。2121共进行了两组仿真试验，实验数据随机产生，r=100，a为[0,π]上均匀分布的随机数，01第一组仿真试验c、a、b(i=2,?,9)为[−1,1]上均匀分布的随机数，第二组仿真试验11i c、a、b(i=2,?,9)为[]−5,5上均匀分布的随机数。两组实验各生成1000组数据，每组数11i由圆周上的180个采样点组成（均匀分布）。对每组数据分别用最小二乘算法、最小二乘迭代算法、最小外接圆算法、最大内接圆算法进行圆度误差评定。统计每种算法的平均相对评定误差δ%和相对评定误差大于8%的次数n。实验结果如表1所示。m表1圆度误差评定仿真试验结果最小二乘迭最小二乘能最小外接圆最大内接圆评定算法代算法用算法算法算法δ0.911.455.125.05第一组n0119192mδ0.974.015.325.35第二组n0289597m表2香烟滤棒圆度误差评定结果（单位：pixels）最小二乘迭最小二乘通最小外接圆最大内接圆评定算法代算法用算法算法算法评定结果10.0912.1112.6712.71备注最小二乘迭代算法得到的圆半径为152.4pixesl分析实验结果，可以得出以下结论：①最小二乘迭代算法、最小二乘通用算法、最小外接圆算法、最大内接圆算法等4种圆度误差评定算法中，最小二乘迭代算法精度最高，最小外接圆算法和最大内接圆算法精度最低；②当满足“小误差假设”和“小偏差假设”时，最小二乘通用算法能达到较高的评定精度，当小误差假设”和“小偏差假设”不满足时，最小二乘通用算法会出现较大的评定误差；③当圆度变形为椭圆或三棱圆时，采用最小二乘迭代算法、最小二乘通用算法、最小外接圆算法、最大内接圆算法等4种圆度误差评定算法，结果非常接近，最大相对评定误差均不超过2%。圆度误差评定结果如表2所示。第一作者邓春宁，男，1974年4月生，工程师，龙岩卷烟厂计量科。'