二次曲面的NURBS最优化表示方法研究

孔德明， 黄紫双， 杨 丹

(1.燕山大学电气工程学院，河北秦皇岛066004；2.根特大学 通信与信息处理系，比利时根特B-9000)

1 引 言

近年来，随着现代化工业的发展，逆向工程在机械零件加工和复杂形状的模型构建等相关领域中得到了广泛的应用。在许多机械零件的模型构建中二次曲面是最为常见、最为基础的几何构型，这些基础构型在制造过程上往往需要达到较高的加工精度。对于同时具有自由型曲面和二次曲面的零件，B样条方法[1]所包括的特例Bézier方法[2]不能精确表示除抛物面以外的二次曲面，只能给出近似的表示进而引入了误差问题。为了解决这个问题，具有多项式表达的有理B样条[3]被提出。其中非均匀有理B样条(non-uniform rational B-spline，NURBS)方法[4]利用非均匀的节点向量表达式构造有理B样条函数，为标准的解析结构和自由型曲面提供了统一的数学表示，适用于各种自由型曲面及组合式曲面模型的构建，并且具有更好的连续性和光顺性。在国际标准化组织(ISO)颁布的工业产品几何定义的STEP标准中，NURBS作为自由型曲线曲面的唯一表示方法在逆向工程中得到了广泛的应用。

目前曲面重构研究主要集中在自由曲面拟合，二次曲面的NURBS拟合方法还存在着拟合精度低、拟合过程复杂的问题。为利用NURBS方法得到高精度的二次曲面，Piegl L等[5，6]提出采用二阶周期NURBS曲线来精确表示球面或椭球面。该方法通过对周期曲线进行平移、旋转和缩放得到目标曲面，所需数据量少，拟合过程简单，拟合效果好。但该方法只能对周期曲线进行修改，仅适用于由一条拟合曲线通过旋转、缩放得到的旋转曲面和对称曲面，对于非旋转类的拼接曲面并不适用。马力全等[7，8]提出通过对数据点插值生成拟合曲面，并采用最小二乘极小化法对构造的曲面进行修正，该方法具有充分的灵活性，可以通过修改输入数据项来逼近曲面的形状，适用于任何二次曲面的拟合，但采用最小二乘逼近需要确定逼近曲面的合适次数才能满足精度要求，实现过程较复杂。

本文根据二次曲面的形状特征选取u、v参数化方向，改善了插值曲面的拟合效果；在避免多次进行控制点调整的前提下，求解得到NURBS拟合高质量二次曲面的最优控制点个数的合理选取范围，简化了拟合过程，为非对称的二次曲面拼接零件的构建提供了数据基础；并以实例曲面为对象，对比分析了本文方法和传统NURBS拟合方法[8,9]的拟合效果和拟合时间。相较于多项式拟合NURBS方法可以任意修改二次曲面拼接类零件的几何形状。

2 方 法

2.1 u、v参数化方向选取

NURBS曲面的参数化方向对应于曲面上参数增加的方向。曲面的参数化方向是不唯一的，选择合适的参数化方向能够改善曲面的拟合精度。曲面的u、v参数化方向选择步骤如下：

Step1：根据要求拟合的二次曲面的一般方程式中的系数求解曲面判别式。

Step2：由判别式结果和判别条件判定曲面的形状特征。

Step3：根据曲面的形状特征选择u、v参数化方向。

根据二次曲面的方程式推断曲面的形状特征，给出二次曲面的一般方程式式(1)，以及判别式(2)～式(4)。

F(x，y，z)=a11x2+a22y2+a33z2+a12xy+a23yz+

a31zx+a1x+a2y+a3z+a4=0

(1)

(2)

(3)

(4)

式中：(x，y，z)是二次曲面上的点坐标；Δ、δ是二次曲面的不变量矩阵，矩阵中的元素a是二次曲面方程系数；Δ0、δ0是位置判别式。

将二次曲面方程式中的各项系数分别代入判别式中进行求解，当判别式结果满足下面两个条件时，所需重构的二次曲面在xy面投影为四边形截面，其余条件下二次曲面在xy面投影则为圆形截面。

条件1：Δ>0，δ=0；

条件2：Δ=0，δ=0；Δ0≠0。

为了精确表示曲面在3个坐标轴平面内各点的曲率，获得较高的拟合精度。根据曲面的形状特征选定当二次曲面在xy面投影为圆形截面时，选择从z=zmax沿侧面曲线到z=0的方向为u方向，沿底面圆形的逆时针方向为v方向。当二次曲面在xy面投影为四边形截面时，选择u、v参数化方向为从原点沿底面四边形两个相邻边长到x=∞、y=∞的方向。

2.2 最优控制点个数选取范围分析

理论上曲面拟合的控制点越多且分布均匀，拟合的精度越高。但控制点的个数过多又会延长计算时间，考虑到实验时间与设计过程将控制点选取分为两个方面进行分析：一是选择最优的控制点分布，在最优的分布条件下确定利用NURBS拟合二次曲面的最优控制点个数的最小取值；二是限定程序运行时间，求解较短运行时间内利用NURBS拟合二次曲面的最优控制点个数的最大取值。最优控制点个数的选取范围求解步骤如下：

Step1：根据3个不重合控制点决定一条拟合曲线的特性确定二次曲面所需的基本控制点个数。

Step2：选取最优位置分布。采用极坐标法将传统直角坐标系的参数值用极坐标参数表示，并代入二次曲面一般方程式进行联合解算得到二次曲面的极坐标方程式。代入均分角度值求得坐标数据在曲面上均匀选取数据点。

Step3：改变极坐标方程式中角度参数θ和φ的大小对数据点个数进行增减。

Step4：利用积累弦长法对数据点进行参数化处理，对得到的节点参数值u′k和v′l采用取平均值法选定合适的节点矢量，将曲面反算问题分解为曲线反算问题求解线性方程组得到控制点，对得到的控制点进行NURBS曲面重构。

Step5：求解最优控制点个数的合理选取范围。利用差值绝对值和均方根误差(RMSE)对不同控制顶点数下生成的各个重构曲面进行定量分析，根据定量分析的结果曲线求解NURBS拟合二次曲面的最优控制点个数的最小取值。结合程序运行时间关系曲线求解NURBS拟合二次曲面的最优控制点个数最大取值。

2.2.1 基本控制点个数确定

NURBS曲面是由多条NURBS曲线在u、v方向上多次构建而成的，是由(m+1)×(n+1)个控制点构成的控制网格。NURBS曲面的方程式定义[4]为

(5)

式中：{Pi，j}是曲面拟合的控制网格；{ωi，j}为权因子均取为1；{Ni，p(u)}和{Nj，q(v)}分别表示节点矢量U和V上的p次和q次非有理B样条基函数；u，v是节点值。

(6)

式中：r=p+m+1；s=q+n+1；U和V中包含的节点总个数分别为(r+1)和(s+1)。

利用NURBS进行曲面拟合通常先获得其低阶曲面，再通过升阶、修改控制点的方式逐步改善曲面的光顺性和拟合精度[9]。因此本文选取p=q=2，以u向为例，结合式(6)可知节点矢量U能够允许的最少节点个数r+1=6，即u方向上拟合3阶NURBS曲线最少需要满足3个控制点。而3个不重合的控制点任意拟合的二次NURBS曲线如图1所示，C(u)是一条具有凸包性的无规则的弧线。同理对于在u方向、v方向上的NURBS曲线结构较为复杂的二次曲面，可以采用外接多边形法将其分别看作g段和h段弧线相连接而成，则曲线所需的基本控制网格顶点数为

(7)

式中：λ为闭合曲线的个数。

图1 3个不重合控制点拟合曲线Fig.1 Fitting curve of three non-coincidence control points

2.2.2 数据点个数范围的选取

相较于直接设计控制点的分布，先确定在曲面上最优分布的数据点再反算求解控制点的方法更加直观简便。最优分布即在曲面上均匀选取数据点，由于在直角坐标系下的均分坐标轴参数法不能实现在弧线上的均匀取点，因此本文采用极坐标取点法,其示意图如图2所示。

图2 极坐标法取点示意图Fig.2 Diagram of polar coordinate method

图2中，点P是选取的曲面上的点，点P′是点P在xy面上的投影。l是点P距离原点的距离，θ是点P在xy平面上的投影点P′距离x轴的角度，φ是点P距离z轴的角度。

将二次曲面一般方程式的直角坐标轴参数xp、yp、zp用极坐标方程表示：

(8)

将式(8)中的坐标轴参数替换方程代入二次曲面的一般方程式，得到式(9)二次曲面极坐标方程：

(a11cosθcos2θsin2φ+a22sin2θsin2φ+a33cos2φ+
a12sinθcosθsin2φ+a23sinθsinφcosφ+
a31cosθsinφcosφ)l2+(a1cosθsinφ+
a2sinθsinφ+a3cosφ+a4)l=0

(9)

通过代入极坐标参数和二次曲面一般方程式中的各项系数进行联合解算，求得点P距离原点的距离l值代入式(8)，进而得到在曲面上均匀选取数据点的坐标数据。

在逆向工程中，常通过进行多次控制点调整以用于在给定允差内对曲线曲面形状进行表示[10]。本文为避免控制点调整给简单二次曲面拟合带来耗时长、计算复杂等问题，从而对二次曲面拟合的最优控制点个数的选取范围进行求解。在均匀选取数据点的基础上，逐渐增加控制点个数生成拟合曲面进行定量分析，判定在某一控制点个数的基础上，增加控制点的个数对误差变化的影响，若误差变化的大小限定在一个很小的范围区间内，且拟合精度满足给定阈值条件，此控制点个数即为NURBS拟合二次曲面的最优控制点个数的最小取值。在最小取值的基础上选取控制点个数进行拟合，选取个数越多拟合曲面的精度越高，但程序运行时间越长。为了节省实验时间，需要结合程序运行时间与控制点个数的关系曲线，给定时间的限定值，求解最优控制点个数的最大取值。

2.2.3 NURBS反算的一般过程

NURBS曲面反算就是指构造一张p×q次NURBS曲面，使其插值于给定的呈拓扑矩形阵列的数据点集{Qk，l}(k=0，1，…，m；l=0，1，…，n)，表示成求解未知控制点Pi，j(i=0，1，…，m+p-1；j=0，1，…，n+q-1)的一个线性方程组。但曲面的线性方程组过于庞大，给求解及计算机上实现带来困难，一般将曲面反算问题分解为分别沿u向和v向的两次曲线反算问题。曲面反算时一般使两个参数化方向上的曲线的首末端点分别与首末数据点Qk，l一致，且数据点将分别依次与二次曲面的节点一一对应[11]。为确定与数据点相对应的节点值，需要先对数据点进行参数化处理。本文选择积累弦长参数化法，这种参数化法可以反映数据点按弦长的分布情况，适用于构造任意阶次的非均匀有理B样条节点矢量，且使所得插值曲线具有较好的光顺性[12]。给出弦长参数化法公式(仅以u向为例[13]):

(10)

(11)

式中：d为总弦长;u′k为Qk决定的节点参数值。

接着采用式(12)取平均值的方法选定合适的节点矢量U，这样就可以建立一个系数矩阵为(m+1)×(m+1)的线性方程组。

(12)

待求的非均匀有理B样条二次插值曲面的线性方程组可写为，

(13)

将式(13)改写为类似于非均匀有理B样条曲线方程的表达式：

(14)

式中u方向上的控制顶点被另一方向上的控制曲线替代：

(15)

首先将式(11)和式(12)计算得到的节点参数值u′k和节点ui依次代入式(16)计算出节点矢量U上的B样条基函数，然后同已知的数据点集{Qk，l}一起代入式(14)就得到了这些控制曲线上m+1个点Ri(v′l)，i=0，1，…，m。这些点又作为数据顶点代入式(15)，运用与节点矢量U上B样条基函数相同的运算公式求得节点矢量V上的B样条基函数，就反算出了曲面上数据点对应的控制点。

(16)

式中

3 具体实现

3.1 3种曲面的u 、v参数化方向选取

根据二次曲面的方程式确定u、v参数化方向。首先给出球面的解析方程式如下：

x2+y2+z2-a2=0

(17)

式中：(x，y，z)是曲面上的点；a是球面半径。

将方程式(17)中各项系数的值代入判别式，求解得到Δ=-a2<0，δ=1>0，Δ0=-3a2<0，δ0=3>0。由判别式结果确定球面的xy面投影为圆形截面，则选择如图3所示的u、v参数化方向。u向为由球面顶点沿侧面弧线指向底面的方向，其侧面圆弧线具有相同的曲率；v向为沿半球底面的逆时针方向。

图3 球面u 、v参数化方向示意图Fig.3 Diagram of u ，v parameter directions of sphere

将圆锥面的解析方程式变换成如下形式：

c2(x2+y2)-a2z2=0

(18)

式中：a是圆锥面的高；c是底面半径。

将方程式(18)中各项系数代入判别式，更易求得Δ=0，δ=-a2c2<0，Δ0=0，δ0=c2-2a2c2不满足所给出的两个判别条件，因此选择如图4所示的u、v参数化方向，u向为由圆锥面顶点沿侧面母线指向底面z=0的方向，v向为沿圆锥面在xy面投影圆形截面的逆时针方向。

图4 圆锥面u 、v参数化方向示意图Fig.4 Diagram of u ，v parameter directions of cone

同样给出椭球面的变换方程式如下：

b2c2x2+a2c2y2+a2b2z2-a2b2c2=0

(19)

式中：a、b、c分别是椭球面的长半轴、中半轴和短半轴。由判别式矩阵求解得到Δ=-a6b6c6<0，δ=a4b4c4>0，Δ0=-a4b4c4(a2+b3+c2)<0，δ0=a2b2c2(a2+b3+c2)>0 不满足条件。则椭球面的u向为由椭球面顶点沿侧面弧线指向底面的方向，v向为沿椭球底面的逆时针方向，见图5。

图5 椭球面u 、v参数化方向示意图Fig.5 Diagram of u，v parameter direction of ellipsoid

3.2 最优控制点个数范围求解

3.2.1 球面选取范围最小值求解

首先定义拟合曲面是一个由(m+1)×(n+1)个控制网格顶点构成的单位球面，设定阶数为3。根据球面标准方程和式(1)可知a11=a22=a33=1，a4=-1，代入式(9)中得到球面上任意一点到原点的距离值l=1。根据设定的控制顶点个数，求解在球面上均匀选取数据点的极坐标参数φi和θj：

(20)

将l值和极坐标参数值分别代入式(8)中，得到均匀选取数据点的极坐标方程式(21)。对式(21)中得到的数据点坐标进行反算得到最优的控制点分布。

(21)

最优控制点个数的求解是通过选取不同控制顶点数生成各个重构曲面，在最优的控制点分布条件下，计算4×104个重构球面坐标数据与标准球面坐标数据之间的差值绝对值进行定量分析。表1给出了拟合球面与标准球面之间的差值平均值和均方根误差，其数值大小随着控制顶点数的增加而减小，根据精度要求设定二者的阈值分别为10-3、0.05。

为了求解拟合球面的最优控制点个数，本文将表1中拟合球面与标准球面之间差值绝对值的平均值及变化率与控制顶点数作为状态变量进行多项式拟合[14]，得到两者之间的关系曲线如图6(a)所示，拟合关系式见式(22)。

表1 球面不同控制顶点数拟合结果分析Tab.1 Analysis of fitting results of different control vertices on sphere dm

图6 球面差值平均值及其变化率与控制点个数的关系曲线Fig.6 Curves of relationship between average of differences，rate of change and number of control points of sphere

当控制顶点数增加时，差值绝对值的平均值随之减小，拟合球面与标准球面之间的误差减小，拟合精度提高。

fs(Ks)=1.79e-0.432 6Ks+0.020 78e-0.037 82Ks

(22)

式中：fs(Ks)是差值绝对值的平均值；Ks为球面的控制点个数。

求解差值平均值曲线的切矢量即对式(22)进行求导，得到差值平均值的变化率和控制顶点数的关系曲线如图6(b)所示，关系式见式(23)。图6(b)中变化率随着选取控制顶点数的不断增加逐渐趋近于0，设定变化率阈值为10-5，当变化率的值小于阈值时，就认为变化率大小几乎不再发生变化；当差值平均值与变化率都小于阈值，这时的控制顶点数就是选取的最小控制顶点数。

f′s(Ks)=-0.774 3e-0.432 6Ks-7.859×10-4e-0.037 82Ks

(23)

将阈值分别代入式(22)、(23)联合解算得到球面控制点个数的选取条件：

Ks>115.395 1

(24)

当选取的控制点个数大于等于116时，利用NURBS方法可以得到满足高精度要求的拟合球面且误差变化率在一个极小的范围内波动。此时拟合球面的最优控制点个数的最小取值对应的控制顶点数是7×20。

不同控制顶点数下的拟合曲面如图7所示，曲面上浅色的点为均匀选取的数据点。其中，当选取的控制顶点数为3×4少于基本控制顶点数时(见图7(a))，v向上的控制顶点不足以拟合一个整圆，拟合球面无法形成一个类似球面。图7(b)是在基本控制定点数3×8下重构得到的类似球面，其与标准球面之间的差值平均值超过所设阈值，不满足精度要求。图7(c)是最优控制顶点最小取值7×20下得到的拟合球面，其与标准二次曲面之间的差值平均值和均方根误差均满足阈值要求，具有较高的拟合精度。图7(d)是控制顶点数为51×200的拟合曲面，其表面选取的大量数据点构成了一个光滑的标准球面。

图7 球面不同控制顶点的拟合结果Fig.7 Fitting results of different control vertices on sphere

3.2.2 其他曲面选取范围最小值求解

定义一个底面半径为1，高为1.5的圆锥面，以及一个长半轴长为1，中半轴长为0.8和短半轴长为0.6的椭球面。同理球面在最优控制点分布的情况下，选取不同的控制顶点数对圆锥面和椭球面进行重构和定量分析。

对比表1～表3中的数据发现，在相同控制顶点数下，形状特征较简单的曲面得到的NURBS拟合曲面精度更高。

通过对表2、表3中控制顶点数与差值平均值进行多项式拟合，得到NURBS拟合圆锥面和椭球面控制点个数的选取条件：

(25)

式中：Kc和Ke分别为圆锥面和椭球面的控制点个数。可知，圆锥面满足高精度拟合要求的最少控制点个数为43，控制网格顶点数为5×12。椭球面满足高精度拟合要求的最小控制点个数为167，控制网格顶点数为8×24。

表2 圆锥面不同控制顶点数拟合结果分析Tab.2 Analysis of fitting results of different control vertices on cone dm

表3 椭球面不同控制顶点数拟合结果分析Tab.3 Analysis of fitting results of different control vertices on ellipsoid dm

如图8所示，利用NURBS方法得到控制顶点数为5×12的拟合圆锥面和控制顶点数为8×24的拟合椭球面，其两者与标准曲面之间的差值平均值和均方根误差均低于所设阈值，具有较高拟合精度。

图8 其他曲面最少控制顶点下的拟合结果Fig.8 Fitting results of the least control vertices on other surfaces

3.2.3 3种曲面选取范围的最大值求解

考虑到程序运行总时间最优，对3种曲面不同控制顶点数下进行拟合的程序运行时间进行实验分析。为保证误差不影响实验结果，求取100次各个重构曲面的程序运行时间的平均值，获得3个二次曲面在不同控制顶点下进行NURBS拟合的程序运行时间与控制点个数的关系曲线如图9所示。

图9 运行时间及与控制点个数关系曲线Fig.9 Curves of relationship between runtime and number of control points

关系多项式为

(26)

由式(26)求解可知，采用本文方法在3种曲面各自的最优控制顶点数下对曲面进行拟合时所需的程序运行总时间均在0.8 s左右。由图9可知:3种曲面的关系曲线在1 s内走势趋于稳定，斜率较小；1 s以后的时间与控制点个数之间的增长关系呈类似指数增长的趋势，斜率较大。为节省实验时间并求得合适的控制点个数选取范围的上限，本文将运行时间限定在1 s以内，求解式(26)得到3个二次曲面允许的最大控制点个数如下：

(27)

综合实例结果可知，根据时间和精度的要求，在均匀选取数据点的前提下，利用NURBS拟合3种二次曲面的最优控制点个数的选取范围是167～677。考虑其他二次曲面的结构复杂程度的不同，在已求得的选取范围基础上允许20%的波动差异，则利用NURBS拟合一般二次曲面的最优控制点个数的合理选取范围是201～541。

3.3 对比分析

利用MATLAB软件采用传统NURBS方法[8，9]对3种二次曲面进行拟合，各曲面选取的控制点个数为上述求得的最小控制顶点数。从差值平均值和均方根误差两个方面对比传统方法与本文方法的拟合效果，对比参数见表4。然后，保证曲面拟合精度一定，比较两种方法所需的程序执行总时间，对比结果如表5所示。从表5中数据可知，与传统NURBS方法相比，本文方法在拟合效果上均有不同程度的提高，运行时间短。

表4 两种方法拟合效果的对比结果统计Tab.4 Statistics of comparative results of fitting effect of two methods dm

表5 两种方法拟合时间的对比结果统计Tab.5 Statistics of comparative results of fitting time of two methods

4 结 论

本文从参数化方向和控制点个数选取范围出发，提出了一种高效拟合二次曲面的方法。实例证明：与传统方法相比，本文方法能够有效改善整体曲面的拟合效果；求解得到的最优选取范围为201～541，选取该范围内的控制点个数对二次曲面进行拟合，无需进行多次控制点调整即可获得高精度的拟合曲面，保证了测量精度，节省了时间；为利用NURBS对二次曲面与自由曲面的拼接曲面进行加工和设计提供了二次曲面拟合的数据基础和理论支持，为后续的曲面升阶、控制点修改奠定基础，对低阶标准解析曲面构建提供了技术借鉴。