平面度误差分析

平面度误差分析
摘要
在生产生活中，有时平面的平整度尤为重要，特别是在某些精密仪器的生产组装中。

我们用平面度误差来表示客观表面具有相对不平整高度对于理想平面的偏差，以Hmax-Hmin来表示，其中Hmax表示测量点至理想平面最大距离，Hmin 表示测量点至理想平面最小距离。

因此，对于测定平面是否满足使用要求，平面度的误差分析必不可少。

本文就平面度误差分析问题，通过建立数学模型，找出较为适宜精确的平面度误差计算方案。

对于问题一，由于平面度的分析是判定相对性问题，则需找出参考的理想平面，由此想到利用常见的最小二乘法可较为直接的利用数据点拟合出理想的平面，因此构建了模型一：利用最小二乘法拟合出理想平面后计算出各点至平面的最大最小距离，两者之和便是所求的平面度误差。

但在检验模型一过程中，我们发现利用最小二乘法求解理想平面时受最大最小点变化扰动较大，使所求矩阵很可能是病态矩阵，不能保证解在最小区域范围内，这往往不能求得最优解，从而使偏差大于实际情况，这限制了模型一的使用范围只适用于变化差异不大的数据点中。

加之题中所测数据多达2500个点，变化较多，计算量较大。

很可能模型一会造成较大误差。

因此，我们改进了利用最小二乘法拟合平面模型，将三维问题简化为二维问题，建立了模型二：将三维点投影到两坐标面上，在坐标面中进行直线拟合，间接求出理想平面，类似模型一，计算出平面度误差。

利用此模型较好的简化了整体计算量，避开了较大差异点，较好的减少了所计算误差，更加直观精确。

对于问题二，我们利用matlab均匀生成了2500个测试点，代入各模型中进行计算，由所求结果可以看出，模型二较于模型一适用范围更广，准确度较好，这较好的验证了结果的准确性，
一、问题重述
1.1问题背景
某工件的某部分是一个100cm*100cm的平面，制作完成后要检测此平面是否合格，也就是要判定这个平面是否足够“平”。

假设在这个平面上采集到2500个均匀分布的数据点（三维），不考虑测量误差，即假定数据都是足够精确的，产品合格的标准为平面度误差不超过ε=0.001mm。

1.请给出一个精确的、直观的、可行的检测方法并说明理由。

2.构造产品合格和产品不合格的数值例子，试验你提出的检测方法。

1.2有关信息
平面度误差是指被测实际表面相对其理想平面的变动量,理想平面的位
置应符合最小条件。

即变动量ε=Hmax-Hmin,Hmax指距离理想平面的最大值，Hmin指距离理想平面最小值。

二、问题分析
依据上述平面度误差定义，我们只需依据所测出的各个数据点，计算出所有点距理想平面的距离，找出最大值和最小值，那么平面度误差ε可由Hmax-Hmin 求得。

也就意味着我们只需根据测量点找出理想平面即可，因此，我们就将平面度误差测量问题转化成了求理想平面问题。

因此只需建立适当的模型，求解出相
可知Hmax-Hmin中结应的理想平面。

根据点到平面的距离公式H=
A2+B2+C2
果与D无关，因此求解可简化为求A,B,C三个未知数的解。

三、模型假设
1.在平面上所采集的点是均匀的，避免取点集中理想平面确定有较大偏差，导致平面度误差测定不准确。

2.工件材料均匀，不可变形，不随外界因素如光照，温度，湿度等变化而改变工件表面的平整度。

3.所测定的数据足够精确，不考虑测量误差。

四、符号说明
(x i,y i,z i)测量点的坐标；
ε平面度误差；
Q 测量点数据的误差平方和； n 测量点的个数；
A,B,C,D,a,b,c模型一中平面方程系数；
A1，B1，A2,B2模型二中两直线方程系数；
→m,→
n模型二中两直线的方向向量；
五、建模思路
要评定一个平面的平面度，必须要找到其参考的基准，即理想平面。

评定的
方法不同，其求出来的平面度误差也不相同。

所以解决此题的关键是用合理的方法确定理想平面，使得所测得的平面误差度最小，利用各个点到理想平面的距离平方和的最小值来确定理想平面，即拟合出唯一的最小二乘平面。

用在平面两侧离平面最远的两点沿平面法向量方向的距离作为平面度误差。

最小二乘法常被用于尺寸与公差的评价。

该问题用此法时，所求的结果唯一且利用数学软件较易实现，便于计算机处理数据。

但是用最小二乘法时，由于本题处理的数据较多，所拟合的平面的形位和其计算的误差和理论的理想平面会有一定差距，不能得出更好的平面解。

在对问题的分析中我们得出，当理想平面的法向量确定后，平面度误差就可以用所有测量点在沿法线方向上的相距最大距离来表示。

这样可以避开对平面位置的确定的影响，只需要找出平面的法向量，其实误差度只与平面方程的两个参量有关，简化了模型，降低了时间复杂度，提高运行效率。

沿着这个思路，我们建立了简单且直观的双直线拟合模型，通过将所有测量点往xOz和yOz平面投影，分别拟合出最小二乘直线，将其叉乘，求出理想平面的法向量，再求出平面度误差。

六、模型的建立与求解
6.1模型一的建立
一般平面方程为Ax+By+Cz+D=0，两边同除D移项得ax+by+c=z,a,b,c为变量,对于点到平面距离H=
a2+b2+12
,(xi,yi,zi)为第i点坐标值，利用最小二乘法,要求出理想平面，则每一点到理想平面坐标的值平方和最小，Q=(zi−z)2,有Q=(axi+byi+c−zi)2，当有理想平面时，Q有最小，根据极小值存在的条
件有ðQ
ða =0,ðQ
ðb
=0,ðQ
ðc
=0,代入Q即有方程组：
ðQ
ða
=2xi(axi+byi+c−zi)=0
ðQ
ðb
=2yi(axi+byi+c−zi)=0
ðQ
ðc
=2(axi+byi+c−zi)=0
展开化简有
a xi 2
+b xiyi +c xi = xizi a xiyi +b yi 2+c yi = yizi a xi +b yi +nc = zi
将此式转化为矩阵形似有a b c *xi 2xiyi xi xiyi yi 2yi xi yi 1=xizi
yizi zi
；
令xi 2xiyi xi xiyi yi 2yi xi yi 1
=A ,xizi
yizi zi =B,则矩阵化为A*abc=B,abc 为a,b,c 构成的列向
量，abc=B/A,利用矩阵的算法可以很容易将系数a,b,c 算出，将算法编码利用matlab 导入点的数据，可以轻松的计算出系数，进而得出平面方程而平面误差度ε=Hmax-Hmin,由于H=zi-z,得出ε=max
a 2+
b 2+12
-min
a 2+
b 2+12。

6.1.1模型一的求解与结果：
将数据点导入matlab 中，编写相应算法（见附件1），在matlab 中实现相应数据理想平面的拟合图像如下：
再利用ε=Hmax-Hmin，计算出每个点的Hi,比较得出Hmax和Hmin,进而算出该实际平面的平面度误差ε计算结果如下所示：
ε=0.0080
理想平面表达式
z= p(1)*x+p(2)*y+p(3) （p=[4.20E-061.25E-060.00642]）
6.2模型二的建立 6.2.1双直线拟合模型：
上述模型一的主要思想是将所测的2500个点进行多元线性拟合，从而得到一个理想平面。

但此模型有如下三点缺陷：
（1）本题中测量平面度误差需要考虑2500个测量点，用二元最小二乘法拟合计算量大且耗时多。

（2）用二元最小二乘法可能得到不可逆的矩阵，从而对结果产生较大的影响。

（3）个别偏离较大的点可能会使最小二乘法产生较大的偏差，使得我们无法得到更好的解。

因此，我们可以分别进行两次一元线性拟合来简化模型一，从而得到模型二——双直线拟合模型。

模型二的基本思想是将所测的所有的点分别投影于xOz 平面和yOz 平面，然后分别对这两个面上的点用一元最小二乘法拟合出两条直线。

假设这两条直线分别为z =A 1x +B 1与z =A 2y +B 2,从而分别得到两条直线的
方向向量→m =(1,0,A 1),→
n
=(0,1,A 2)。

由立体几何公理可知，空间中相交的两条
直线构成一个平面；又由叉乘知识可知该平面的法向量→r =→m ×→
n
=(-A1,-A2,1)。

设该平面的方程为z =Ax +By +E ,又因平面法向量可知，所以只有C 未知。

然而我们由平面度的定义可知常数E 对平面度ε并无影响，故无需算出E 。

6.2.2一元线性回归的推导：
对于投影在xOz 平面上的点M ix (x i ,0,z i )来说，要求出它们所拟合出的理想直线z =A 1x +B 1，只需每一点到理想直线坐标的值平方和最小。

设Q =
(z i −z )2=211()i i A x B z +-∑，当有理想直线时，Q 有最小，根据极小值存在的条件有
11
0,0Q Q
A B ∂∂==∂∂,代入Q 即有方程组：
∂Q
∂A 1
=2x i (A 1x i +B 1−z i )=0∂Q
∂B 1=2 (A 1x i +B 1−z i )=0 展开化简可为：
211i i i i A x B x x z +=∑∑∑11i i
A x n
B z +=∑∑
将此式转化为矩阵
1
A B 1
2i
x ∑i
x
∑i
x ∑n
=
i i
x z ∑i
z
∑
设A=
1
A B 1
,B=
2i
x
∑i
x
∑i
x
∑n
,C =
i i
x z ∑i
z
∑ ;
所以我们有A B=C .又由C/B 可知A ，从而得知A 1与B 1的值，从而得知该直线的方程。

同理，对于投影在yOz 平面上的点M ix (x i ,0,z i )，亦可求得A 2与B 2的值。

6.2.3平面方程的推导：
对于xOz 平面上的直线z =A 1x +B 1,可知其法向量→
r 1
=(A1,0,-1)。

设该直线
的方向向量→m =(s,0,t),又由→r 1.→m
=0可求得→
m =(1,0,A 1)。

同理可知，yOz 平面上
的直线z =A 2y +B 2的方向向量→
n
=(0,1,A 2)。

由叉乘知识可得上述两条直线所构成的平面的法向量
→r =→m ×→n = →e 1→e 210→
e 3A 1
01A 2 = −A 1 ,−A 2,1 ,所以可知该平面方程为z =A 1x +A 2y +E 。

最后只需将所测点M i (x i ,y i ,z i )代入可求得理论值z ii ，平面度误差ε=max z i −z ii −min ⁡( z i −z ii 。

6.2.4模型二的求解与结果：
将测量数据写入表格中，取其中20个点作为测试模型数据点，如下：
x-z 与拟合曲线图像：
拟合直线结果：
z=p(1)*y+p(2) p=[ 4.2e-006 0.00672]
y-z与拟合曲线图像
拟合直线结果：
z=q(1)*y+q(2) q=[1.25e-006 0.00705] 求得理想平面方程：
z=p(1)*x+q(1)*y+p(2)
如图：
求得平面度误差为：
ε=0.00713
七、模型的检验
7.1问题二求解
对于问题二，我们利用matlab生成均匀的的2500个点，使每个点坐标的x 值均匀在[0,100]中，y坐标均匀取在[0,100]中，z分别随机取在[0,0.08]和[0,0.012]范围内,这样就构成了两组2500个点的均匀数据点（见附件1），一组处于误差范围内，另一组处于误差范围之外，为了检验所建立的模型，将数据分别导入不同模型中可分别计算出不同模型的算出的平面度误差，拟合的理想平面系数及计算结果如下表所示：
7.1.1 z分别随机取在[0,0.012]范围内
2500点的拟合平面
7.1.1 z分别随机取在[0,0.008]范围内
散点图
拟合平面图
由上表比较不同模型所算出的平面度误差大致一致，说明模型建立大致无差
错，这与我们预期是相符合的，但是，由于数据点多达2500点，因此在程序运
行过程中明显能够感觉到改进之后的双直线拟合法比单纯的利用最小二乘法进
行拟合快的多，明显提高了效率，
在模型运用中，我们显然也应该知道模型所适用的范围，
由模型一的性质可
知，模型一不适用于具有差异较大点的判断，在模型检验中因此我们特意在第一
组数据基础上添加了几个差异相当大的点，以相同的方法，利用两种模型计算结
果如下表所示：
上表对比表可知，在改变仅仅几个数据点，平面大体不变的情况下，以模型一计算得出的结果和原结果差异甚大，以模型二计算结果改变较小，而以实际如在一平面边角处磕碰等可造成点的数据差异较大，但这对平面的平面度误差影响较小，因此在一定情况下，模型二比模型一更加合理。

八、模型的评价与推广
8.1 模型的评价
由表中的数据可知，结果基本符合事实，因此所建立的模型基本上是可行有效的。

我们所建立的模型采用较为一般简单的方法，切实可行。

由上文可知模型一适用范围较小，但在一般情况下，在严格平面制造工序中，由于工艺精度要求，初检时，出现明显凸起或凹下（也就是较大差异）的平面大部分都会被直接废弃，因此在测量检验时，大部分平面已是较为平整，所以在平面测量点中出现较大差异点的情况较小，那么也就是说模型一仍是可适用的。

我们以最常见的最小二乘法为基础建立了较为简单且可行的模型，有效的解决了平面度误差计算问题，只需将所测量的数据导入模型中，利用计算机可十分简便的给出平面度误差结果，且结果能与预期较好的吻合。

虽然最小二乘法在平面度误差评定上具有一定的优越性，但在某些精密性仪器平面中，精度要求相当高，若此时只用最小二乘法评定，则明显是达不到要求的，因此找出一种更加精密的模型是有意义的，但由于时间和水平有限，我们在尝试过程中始终不能建立一个较为完备的模型，但在查阅资料的过程中，我们知道利用粒子群算法可以较好的优化方案，由于编程水平不够，我们不能很好的编写出相应的代码实现算法，但这确实是改进方向，也可以从遗传算法入手，同样可以很好的提高结果的精确度。

上述所提供的方法都是建立在以数学算法找出理想平面的优化，同样我们可以考虑在几何上更加直观的解决问题，譬如以一束光入射，经过平面反射，由于平面的不平整度，必然会发生较于理想平面的偏转，测量出偏转角度，以对应关系则可得出相应点的高度差，进而得出平面度误差，当然这只是理论思路，可行性还有待后续检验实现。

8.2 模型的推广
由于模型的建立是基于最常见简单的最小二乘法，因此，模型的推广也相对简单，我们的模型是不光对于平面而言，对于比如圆，球体，圆柱等的曲面同样适用，对应于平面度，相应为圆度，球度等，下面针对圆度，球度利用最小二乘法模型进行推广，其他类似可知。

8.2.1圆度模型建立
最小二乘法评定圆度误差方法是按测得的轮廓坐标值，出该轮廓的最小二乘方圆，并以该圆为理想圆（评定基准圆）来评定圆度误差的一种方法，如图
最小二乘方圆也称平均圆。

当被测轮廓上各点至某一圆的距离平方和为最小时，该圆即为最小二乘方圆。

若轮廓上某点至该圆圆心有最大距离Rmax ，另一点有最小距离Rmin ，则圆度误差值ε=Rmax-Rmin, 被测轮廓用最小二乘方圆评定其圆度误差时，主要是确定最小二乘方圆的圆心和半径。

在图中，被测轮廓上各采样点为（xi,yi ）,极坐标值为（ri,θi ）,θi =2πn ×i i =1,2,3….n
设最小二乘方圆的圆心为O ′,坐标为(a,b),半径为R,坐标原点O 到O ′距离为c,c 2=a 2+b 2。

利用最小二乘法可以相似于式1推出
a =1n xi
b =1n yi R =1n
ri 利用极坐标转化得
a =2n ricosθi a =2n risinθi R =1n
ri
以ri 表示点到回转中心的距离，那么最小偏差ei =ri −R −acos θi −bsinθi , 点到最小二乘圆圆心距离为Ri = (a −xi )2−(b −yi )2
Ri 的最大差或ei 中的最大值，即为按最小二乘方圆评定的圆度误差值。

按最小二乘方圆判定时，都用计算法求得圆度误差值。

8.2.2球度模型建立
最小二乘法评定球度误差方法是求出该轮廓的最小二乘法球，并以该球为理想球（评定基准球）来评定球度误差。

当被测轮廓上各点至某一球的距离平方和为最小时，该圆即为最小二乘方球。

若轮廓上某点至该球球心有最大距离Rmax ，另一点有最小距离Rmin ，
则球度
误差值ε为：ε＝Rmax-Rmin
坐标原点到最小二乘法拟合出的理想球面的球心距离为d,理想球体球心为（a,b,c）那么d2=a2+b2+c2。

同理利用最小二乘法推出
a=1
n
xi
b=1
n
yi
c=1
n zi
, R=1
n
ri
转化为极坐标表示有
a=2
n
risinφicosθi
b=2
n
risinφisinθi
c=2
n ricosφi
, R=1
n
ri
由采样点至最小二乘法球球心的距离为Ri=(a−xi)2+b−yi2+c−zi2, 误差值即为Ri中最大差值，表示为球度。

附录
在建立数学模型求解问题时，使用了windows 7，处理器 Intel core i5环境下的matlab 7.0软件的命令窗口进行操作。

附件一：程序代码
1.用最小二乘法求理想平面计算结果并绘制代码
p=[x y ones(length(x),1)]\z %求系数
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y )),'v4'); %插值绘图
figure, surf(X,Y,Z);
hold on;
m=0:10:300;
n=0:12:480;
[mm,nn]=meshgrid(m,n);
o=p(1)*mm+p(2)*nn+p(3);
surf(mm,nn,o);
2.求平面误差度代码
A=[x',y',ones(size(x'))];
B=z';
abc=A\B;
a=abc(1);
b=abc(2);
c=abc(3);
%至此,方程可写为a*x+b*y+c-z=0
F=(a*x'+b*y'+c-z')/sqrt(a^2+b^2+1^2);
s=max(F)-min(F)
3.双直线拟合代码
m=polyfit(x,z,1)
n=polyfit(y,z,1)
p=cross(m,n)
m=0:2:100;
n=0:2:100;
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y )),'v4'); %插值绘图
figure, surf(X,Y,Z);
hold on;
[mm,nn]=meshgrid(m,n);
o=p(1)*mm+p(2)*nn+p(3);
surf(mm,nn,o);
F=( p(1)*mm+p(2)*nn+p(3)-z')/sqrt(p(1)^2+p(2)^2+1^2);
s=max(F)-min(F)
4.2500点拟合代码和图像绘制代码
figure
[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y )),'v4'); %插值
pcolor(X,Y,Z);
shading interp %伪彩色图
figure,contourf(X,Y,Z) %等高线图
figure,surf(X,Y,Z);%三维曲面
scatter3(x,y,z) %散点图
r=[x y ones(length(x),1)]\z %求系数
F=( r(1)*mm+r(2)*nn+r(3)-z')/sqrt(r(1)^2+r(2)^2+1^2);
s=max(F)-min(F)
hold on;
[mm,nn]=meshgrid(m,n);
o=p(1)*mm+p(2)*nn+p(3);
surf(mm,nn,o);。