广义B样条

广义B样条（共7篇）

广义B样条 篇1

引言

基于广义线性模型的统计分析常常被用于生物、医学、金融工程等各个领域。该模型有两个基本的假定:一是在给定协变量X的条件下, 响应变量Y的条件分布属于某一指数分布族;二是通过某一变换, 该模型可以转化为线性回归模型。但是, 在实际问题分析中, 假定协变量X与响应变量Y满足某一线性结构的参数模型往往是不成立的。如下的广义变系数模型是经典广义线性模型的一个简单而有意义的推广。

其中β (u) = (β1 (u) , …, βP (u) ) T是未知的函数系数向量, g (·) 是一个已知的单调联系函数, X和U为协变量, Y为响应变量, ε为模型误差, 并且满足E (ε|X, U) =0。

注意到在模型 (1) 中, 回归系数可以随着协变量U的变化而变化, 因此该模型在表达协变量与响应变量的关系上, 具有更强的适应性和解释能力。另外, 显然当β (u) ≡β为某一未知参数向量时, 模型 (1) 就转化为经典的广义线性模型。因此, 模型 (1) 可以看成经典广义线性模型的一个推广形式。

目前, 模型在统计领域中的应用越来越受到人们的关注, 并且关于模型 (1) 的研究也已有大量的文献。比如Cai等研究了模型 (1) 的有效估计以及拟合优度检验问题。Lian在高维数据下, 考虑了模型 (1) 的变量选择问题。Kuruwita等则在联系函数未知的情况下, 来研究模型 (1) 的估计问题。在本文中, 我们主要考虑模型的B样条估计问题。结合基函数逼近技术以及拟似然方法, 对模型 (1) 提出了一个B样条估计方法.并结合R统计软件, 对所提出的估计方法给出了一个基于R软件的实现过程, 数据模拟表明所提出的B样条估计方法以及R软件的实现过程均是可行的。

与核估计以及局部多项式估计等局部非参数估计方法相比, B样条估计有其独特的优势。比如, B样条估计属于整体非参数估计方法, 在估计过程中可以充分利用样本的信息, 因而具有相对较高的估计精度。关于B样条估计的更多性质, 参见文献[6]。

一、方法论和主要结果

(一) B样条估计过程

设 (Xi, Ui, Yi) , i=1, …, n为来自模型 (1) 的一个容量为n的样本, 并且记B (u) = (B1 (u) , …, BL (u) ) T为阶数M的一组B样条基函数, 其中L=K+M, K为内部节点的个数, 那么βK (u) 可以近似表示成:

其中γk为基函数系数, 把 (2) 代入模型 (1) 可得:

其中进一步, 记μi=E (Yi|Wi) ≈g (WTiγ) , V (μi) =Var (Yi|Wi) ≈V (g (WTiγ) ) 。

那么关于γ的拟似然函数可以定义为:

其中定义为为最小化公式所得的解, 那么函数系数βK (u) 估计为

(二) R软件的实现过程

在统计分析领域, 常用的统计软件有SAS、SPSS、S-Plus以及R等。其中R软件是一套完整的数据处理、计算和制图软件系统。其功能包括:数据存储和处理系统、数组运算工具、完整连贯的统计分析工具、优秀的统计制图功能、简单而强大的编程语言等功能。在Tiobe公布的2011年11月编程语言排行榜上, R语言位列第二十七位, 市场占有率是0.5%, SAS是0.386%, 排在第三十五位。可见R软件是目前用户增长最快, 且已成为了统计学家和数据分析师所钟爱的统计软件之一。接下来, 我们利用R软件, 对上节所提出的估计方法给出一个具体的实现过程。

首先我们以3次B样条为例, 给出B样条基函数的程序代码, 其他阶数的B样条基函数程序可以类似地给出。程序中u为样条函数的自变量, m为样条函数的阶数, kk为内部节点个数, ul为自变量u的取值下限, uu为自变量u的取值上限。

接下来, 我们给出函数系数β (u) 估计的程序代码, 其中x和u为协变量, y为响应变量, m为B样条的阶数, kk为内部节点个数。在该程序中, 我们采用的是等间距节点。

二、模拟研究

接下来我们通过数字模拟实验来说明本文提出的估计方法以及算法的可行性。为实施模拟, 我们假定数据满足的广义变系数模型为变系数泊松回模型, 即响应变量Y服从泊松分布, 并且假定其条件均值具有如下结构:

其中β (u) =sin (2πu) 。另外, 在产生数据时, 假定协变量服从X～N (1, 1) 和U～U (0, 1) 。在估计过程中我们采用3次B样条基函数, 并且取等间距的内部节点, 内部节点个数取为K=2, 样本量取为n=200。该模拟过程的R程序代码如下:

运行上述程序代码, 则得出下页图1的模拟结果, 其中实线代表真实的函数曲线, 点虚线代表利用本文的方法给出的估计曲线。从图1可以看出, 估计曲线可以很好地拟合真实曲线。这也表明本文提出的估计方法以及给出的R程序代码是行之有效的。

参考文献

[1]Cai Z, Fan J Q, Li R.Efficient estimation and inferences for varying-coefficient models[J].Journal of the American Statistical Associ-ation, 2000, 95:888-902.

[2]Lian H.Variable selection for high-dimensional generalized varying-coefficient models[J].Statistica Sinica, 2012, 22:1563-1588.

[3]Kuruwita C N, Kulasekera K B, Gallagher C M.Generalized varying coefficient models with unknown link function[J].Biometrika, 2011, 98:701-710.

[4]Wu C O, Chiang C T.Kernel smoothing on varying coefficient models with longitudinal dependent variable[J].Statistica Sinica, 2000, 10:433-456.

[5]Fan J, Gijbels I.Local Polynomial Modeling and Its Applications[M].Chapman and Hall:London, 1996.

[6]Schumaker L L.Spline Functions[M].Wiley:New York, 1981.

B样条曲线的绘制算法及其实现 篇2

1、B样条曲线定义

B样条曲线由分段曲线组成, 具有Bezier曲线的优点, 以克服了其不具备局部性质的缺点。p次B样条曲线的数学方程[1]为:

其中Pi为控制顶点, 曲线有n-k+1段;Ni, p (u) 称为B样条基函数, 由节点矢量确定.B样条基函数Ni, p (u) 由德布尔-考克斯递推公式计算[1]:

2、计算曲线上点的算法

已知节点向量给定参数u, p次B样条曲线的基函数可由式 (2) 递推得到。其matlab算法如下:

算法1

其中i表示参数u所在的区间[Ui%, Ui+1], u为所求参数, p为B样条曲线的次数, U为节点向量, 返回值N表示在区间[Ui, Ui+1]上, 参数u对应的p+1个基函数值。由算法1容易求得参数u所对应的B样条曲线上的点。

3、参数选取

从算法1可知, 只要给出参数, 就可以计算得到B曲线上对应的点。通常的绘制方法是用在参数上均匀取点, 计算得到对应的曲线上的点, 然后连线得到直线段来近似代替原曲线。这种方法得到的直线与原曲线的误差较大。本文提出参数选取的方法, 就是使得直线段与原曲线的误差限制在一定的误差值tol内。主要思想是:用节点插入方法[1]把B样条曲线分解成几段Bezier曲线Qi (u) 。计算每段Qi (u) 到其两端点连线Q0Qp的最大距离是比较耗时的.由Bezier曲线的凸包性可知:

max (d (Qi, Q0Qp) , i=1, …, p-1) 表示顶点到边Q0Qp最大的距离, 是容易计算的。因此, 当max (d (Qi, Q0Qp) , i=1, …p-1) ≤tol时, 参数的分割终止。否则, 在该参数区间中点插入节点, 直至满足误差小于tol的要求时终止。算法过程如下:

算法2

1.利用节点插入方法的B样条曲线分解成几段Bézier曲线.

2.对每段Bézier曲线B (Q0, …, Qp)

用曲线B (Q0, …, Qp) 初始化栈;

while栈非空 do

曲线B出栈;

if B控制顶点误差<tol then保存参数u;

else在参数的中点插入节点, 分成两段Bézier曲线;

把两段Bézier曲线再存入栈;

End

4、数值例子

例采用本文所述算法编程, 对平面3次B样条曲线在不同误差下绘制。控制顶点为 (-3.5, 0) , (-2.0, -47) , (0, 10.7) , (2.0, -4.7) , (3.5, 0) ;节点向量为 (0, 0, 00, 0.5, 1, 1, 1, 1) .图1为tol=0.1时, 得到对应折线段数为16.图2为tol=0.01时, 得到对应折线段数为46.从图1和2看出, 在曲线较平的地方, 折线分割较少;在曲线弯曲的地方, 折线较密集。表1列出了误差tol取不同值下的折线段数, 折线段数随着tol取值的减少而增加。

5、结论

本文实现了一种绘制B样条曲线的方法, 提出在指定误差限内绘制平面或三维空间中的B样条曲线。数值例子显示本文方法优于等分参数区间方法, 程序是高效、实用的。按照文中提供的方法可推广到有理B样条曲线绘制, 由于篇幅有限, 本文不再赘述。

摘要：B样条曲线是广泛流行的形状数学描述方法。本文对B样条曲线原理进行研究, 并给出算法终止条件。最后用matlab语言设计和实现了B样条曲线.数值例子显示本文程序是稳定和实用的。

关键词：B样条曲线,算法,Matlab

参考文献

[1]Piegl, Tiller.The NURBS book[M].2nd ed..New York:SpringerVerlag, 1997.

[2]孙家广.计算机图形学 (第三版) [M].清华大学出版社, 2003.

[3]施法中.计算机辅助几何设计与非均匀有理B样条[M].高等教育出版社, 2001.

[4]张玉连.Bezier曲线的一种快速绘制算法及其实现[J].燕山大学学报, 2000;24 (3) .

代数曲线的有理二次B样条逼近 篇3

曲线的参数和代数表示在几何造型中各有优缺点,CAGD中同时采用曲线的两种表示法,两种形式间的转换成为必须解决的课题。参数曲线的代数化[1]总是可以实现的,但大部分代数曲线都不能精确参数化,近似参数化问题由此产生。

相关研究工作已取得一些有价值的成果。如线性逼近方法[2];多边形逼近技术[3],可变的B样条曲线逼近方案[4],特殊情形的有理逼近手段[5],分段有理逼近途径[6],区间三次Bézier曲线逼近算法[7],分段有理二次曲线逼近理论[8]。

经过对比分析不难发现,有些逼近算法的计算复杂度较高,有些计算必须构造数值方法才能解决,如“肩点”等问题的计算。这里提出的有理二次B样条逼近算法侧重于曲线段“中点”的信息,规避了复杂的“肩点”计算问题。同时,逼近曲线保持了原始曲线的一些重要几何性质,如凹凸性,单调性和G1连续性;递归调用逼近算法,可以将误差控制在指定的范围之内。

1 代数曲线分段

1.1 曲线的凹凸和单调区间分段

已知代数曲线

$\tilde{A}$:F(x,y)=0 x∈[a0,an] (1)

假定根据凹凸性将区间[a0,an]分割为n段

[a0,an]=[a0,a1]∪[a1,a2]∪…∪[an-1,an] (2)

这里,ai∈[a0,an] i=1,2,…,n-1表示曲线拐点的横坐标。

在凹凸区间

[ai,ai+1] i=0,1,…,n-1 (3)

上,第i段曲线表示为:

$\tilde{A}{}_i$:F(x,y)=0 x∈[ai,ai+1] (4)

设曲线段(4)在区间(3)上的极值点为ami,将区间(3)进一步分割(若已经单调,则不再分割)为如下形式:

[ai,ami]∪[ami,ai+1] i=0,1,…,n-1 (5)

则曲线段(4)表现为两段:

$\begin{array}{l}f(x,y)=0x\in [a{}_i,a{}_{mi}]\\ f(x,y)=0x\in [a{}_{mi},a{}_{i+1}]\end{array}(6)$

这样,所给代数曲线在每个分段区间上具有固定的凹凸性和单调性。

1.2 曲线段的三角形凸包假定分割后的某条代数曲线段为

A:F(x,y)=0 x∈[a,b] (7)

设它的两个端点为V0,V2,求出曲线段在两端点处切线,记为:

$\begin{array}{l}{\rm T}{}_0:a{}_{0}x+b{}_{0}y+c{}_0=0\\ {\rm T}{}_2:a{}_{2}x+b{}_{2}y+c{}_2=0\end{array}(8)$

容易证明,切线T0,T2必相交于一点。

记交点为V1。这样,由三点V0,V1,V2形成了曲线段A的三角形凸包(如图1)。

例如,给定一条平面代数曲线段(如图2(a)):

C:x6+y6+3x4y2+3x2y4-x4-y4+2x2y2=0 (9)

它的右四分之一(如图2(b))分割为四段(如图3(a)),各段曲线的三角形凸包如图3(b)所示,其中A,B,C,D四点坐标依次为$(0‚0)‚(\sqrt{30}/9‚\sqrt{6}/9)‚(1,0),(\sqrt{30}/9,-\sqrt{6}/9)$。

2 曲线段的有理二次B样条逼近

2.1 有理二次B样条曲线

设控制顶点为:

Qi(xi,yi) i=0,1,2

权因子为:

wi(>0) i=0,1,2

则有理二次B样条曲线[4,5,6,7]P(u)u∈[0,1]可以表示为:

$\{\begin{array}{l}x(u)=\frac{w{}_{0}x{}_0(1-u){}^2+w{}_{1}x{}_1(1+2u-2u{}^2)+w{}_{2}x{}_{2}u{}^2}{w{}_0(1-u){}^2+w{}_1(1+2u-2u{}^2)+w{}_{2}u{}^2}\\ y(u)=\frac{w{}_{0}y{}_0(1-u){}^2+w{}_{1}y{}_1(1+2u-2u{}^2)+w{}_{2}y{}_{2}u{}^2}{w{}_0(1-u){}^2+w{}_1(1+2u-2u{}^2)+w{}_{2}u{}^2}\end{array}(10)$

它的基本性质是:起始点P0=P(0)和终止点P1=P(1)分别位于Q0Q1,Q1Q2上,并和Q0Q1,Q1Q2相切。要确定与分段曲线A对应的B样条曲线段P(u),首先必须确定控制顶点和权因子。

根据Farin和Worsey提出的B样条曲线的参数化标准形式[8],当曲线首末两端的权因子为1时,曲线具有较好的参数化效果。这里取

w0=w2=1 (11)

2.2 控制多边形的确定

已知A的三角形凸包V0V1V2,设:

$V{}_i=(\tilde{x}{}_i,\tilde{y}{}_i)i=0,1,2$

取控制顶点 Q1=V1 (12)

则$x{}_1=\tilde{x}{}_1,y{}_1=\tilde{y}{}_1$,由V0=P0,V2=P1得:

$\frac{Q{}_0+w{}_{1}Q{}_1}{1+w{}_1}=V{}_0‚\frac{w{}_{1}Q{}_1+Q{}_2}{1+w{}_1}=V{}_2$

即有:

$x{}_0=(1+w{}_1)\tilde{x}{}_0-w{}_1\tilde{x}{}_1‚y{}_0=(1+w{}_1)\tilde{y}{}_0-w{}_1\tilde{y}{}_1(13)$

$x{}_2=(1+w{}_1)\tilde{x}{}_2-w{}_1\tilde{x}{}_1‚y{}_2=(1+w{}_1)\tilde{y}{}_2-w{}_1\tilde{y}{}_1(14)$

这样,控制顶点将借助于下文的权因子w1的确定而确定(如图4),曲线P(u)中的w1的自由变化对应代数曲线段A的逼近曲线族(如图5)。

2.3 “中点”到“中点”的最小距离逼近

这里,A上的“中点”是指SF(x(a+b)/2,y(a+b)/2);B样条曲线P(u)上的“中点”是指P(1/2),记作S1/2。“中点”到“中点”的最小距离逼近的条件:

$\begin{array}{l}\text{Μ}\text{i}\text{n}:f(w{}_1)=|S{}_{1/2}S{}_F|{}^2=\left(\frac{x{}_0+6x{}_{1}w{}_1+x{}_2}{2+6w{}_1}-x{}_{(a+b)/2}\right){}^2+\\ \left(\frac{y{}_0+6y{}_{1}w{}_1+y{}_2}{2+6w{}_1}-y{}_{(a+b)/2}\right){}^2(15)\end{array}$

即:

$\frac{\partial \text{f}}{\partial \text{w}{}_1}=0(16)$

将(13)、(14)代入(16)可以得到关于w1的方程,求出w1,并由(12)、(13) 、(14)确定控制顶点为:

Qi(xi,yi) i=0,1,2

最后由(10)确定P(u)。

2.4 逼近曲线的性质

1) 保持原有曲线的凹凸性不变;

2) 保持原有曲线的单调性不变;

3) 在曲线的分段点处保持G1连续性[6];

4) 逼近曲线段的“中点”到原始曲线段的“中点”距离最小。

3 分段逼近算法

3.1 逼近误差

在几何学中,曲线段C和它的近似表示Ca之间的误差经常使用Hausdorff距离来描述

$e(C,C{}_a)=dis(C,C{}_a)=\underset{{\rm P}\in C}{{\displaystyle \max }}\underset{Q\in C{}_a}{{\displaystyle \min }}d({\rm P},Q)$

这种距离便于理论分析但不便于计算。

设A为代数曲线段,P(u)为逼近曲线段,下面给出一种易于计算的误差概念。

定义1 在参数域[0,1]上任意取一点

ui=i/n 0≤i≤n n∈N

记P(u)上的对应点为(x(ui),y(ui)),A上的对应点为$(x(u{}_i),\widehat{y}(u{}_i))$,定义逼近误差为

$e(A,{\rm P}(u))=\underset{i}{{\displaystyle \max }}$$\widehat{y}(u{}_i)-y(u{}_i)|(17)$

3.2 算法描述及其说明

输入代数曲线$\tilde{A}$以及误差限δ,输出误差小于δ的分段逼近的有理二次B样条曲线。

1) 由(1)至(6)式对原始代数曲线进行分段(涉及拐点和极值点的计算);

2) 计算曲线段A两端切线T0,T2及其交点V1,确定三角形凸包V0V1V2;

3) 根据(13)、(14)、(16)式确定权因子w1;

4) 由(12)、(13)、(14)式确定Q0,Q1,Q2;

5) 依据(10)式计算逼近曲线P(u);

6) 若逼近误差e(A,P(u))<δ,则计算过程中止;否则,采用逐步二分分段区间的方法,递归调用上述分段逼近算法,直至满足给定的误差要求。

3.3 算法的收敛性

定理1 在上述算法的实施过程中,A和P(u)之间的Hausdorff距离收敛于0。即如果设二分k次后的Hausdorff距离为Ek(A,P(u)),则

$\underset{k\to +\infty }{{\displaystyle \lim }}E{}_k(A,{\rm P}(u))=0(18)$

证明:如图6,由代数曲线的分段方法及其凸包形成过程可知,△V0V1V2是确定的。设:

∠V1V0V2=α,∠V1V2V0=β 0<α,β≤π/2

记边V0V2的高为h,则:

h=V0V1sinα h=V2V1sinβ

记“中点”SF到边V0V2的距离为d。设经过SF的切线分别和V0V1、V2V1相交于P1,P2,自P1,P2分别向V0SF、V2SF引垂线,记P1H1=d11,P2H2=d12。

又设∠P1V0SF=α1,∠P2V2SF=β1 0<α1,β1<π/2。

显然E0(F,P(u))<d<h,经过在“中点”处一次分割后,

d11=P1V0sinα1=a1V1V0sinα1 0<a1<1

d12=P2V2sinβ1=b1V1V2sinβ1 0<b1<1

经过k次分割后,

dk1=a1a2…akV1V0sinαk 0<ai<1,i=1,2,…,k

dk2=b1b2…bkV1V2sinβk 0<bj<1,j=1,2,…,k

设$a=\underset{i}{{\displaystyle \max }}\{a{}_i\}‚b=\underset{i}{{\displaystyle \max }}\{b{}_i\}$,显然0<a,b<1,于是

$\begin{array}{l}\underset{k\to +\infty }{{\displaystyle \lim }}d{}_{k1}<\underset{k\to +\infty }{{\displaystyle \lim }}a{}^k|V{}_{1}V{}_0|\sin \alpha {}_k=0\\ \underset{k\to +\infty }{{\displaystyle \lim }}d{}_{k2}<\underset{k\to +\infty }{{\displaystyle \lim }}b{}^k|V{}_{1}V{}_2|\sin \beta {}_k=0\end{array}$

从而,$\underset{k\to +\infty }{{\displaystyle \lim }}E{}_k(F,{\rm P}(u))=0$。

4 实验结果

逼近算法的求解(参见(15)或(16)式)不存在实质性的困难,算法容易实施。

图3中曲线C的右四分之一分段一步有理二次B样条曲线的逼近结果参见图7,图8为误差函数图形,一步逼近误差

$e{}_1(A,{\rm P}(u))=\underset{i}{{\displaystyle \max }}$$\widehat{y}(u{}_i)-y(u{}_i)$<0.015

5 结 论

代数曲线的有理二次B样条逼近算法计算简单,保持了原始曲线的一些重要几何性质,通过算法的递归调用,可以将逼近误差控制在指定的范围之内。

参考文献

[1] Sederberg T,Zheng J.Algebraic methods for computer aided geometric design.In:Farin G,Hoschek J,Kim M s.Handbook of Computer Aided Geometric Design.North-Holland,Amsterdam,2002.

[2] Ihm I,Naylor B.Piecewise linear approximations of digitized space curves with applications.In:Patrikalakis,N.Scientific Visualization of Physical Phenomena.Springer-Verlag,Berlin,1991:545-569.

[3]Farouki R.Hierarchical segmentations of algebraic curves and some ap-plications.In:Lyche T,Sshumaker LL.Math.Methods in Computer Ai-ded Geometric Design.Academic Press,Boston,MA,1989:239-248.

[4] Pottmann H,Leopoldseder S,Hofer M.Approximation with active B-spline curves and surfaces.In:Coquillart S,Shum H,Hu S M.Proc.of Pacific Graphics 2002.IEEE press,Los Alamitos,CA,2002:8-25.

[5] Sederberg T,Zhao J,Zundel A.Rational approximation of algebraic curves.In:Strasser W,Seidel H.Theory and Practice of Geometric Modeling.Springer-Verlag,Berlin,1989:33-54.

[6] Bajaj C,Xu G L.Piecewise rational approximations of real algebraic curves.J Computational,1997,15(1):55-71.

[7]Chen F,L.Interval parametrization of planar algebraic curves,In:Li ZM,Sit W.Proceeding of the 6th Asian Symposium on Computer Mathe-matics.World Scientific,Singapore,2003:64-76.

[8] Xiao Shan Gao,Ming Li.Rational quadratic approximation to real algebraic curves.Computer Aided Geometric Design,2004,21:805-828.

[9] Les Piegl ,Wayne Tiller,The NURBS Book [M].New York:Springer,1997.

函数的分段有理二次B样条插值 篇4

科学和工程计算中,函数的近似表示一直是一个重要课题。近似方法一般可归结为插值、逼近和拟合三种基本类型,经历长期发展,函数逼近方法[1,2,3]十分丰富。诸如Taylor公式、Lagrange-Newton插值、多项式样条方法、最佳平方逼近、最小二乘拟合、连分式插值、有理样条表示等等,它们在函数近似表示方面均发挥着不同程度的作用。

有理B样条曲线[4,5,6,7]是计算机辅助几何设计中常用的一元函数的参数形式,它具有良好的几何性质以及比较成熟的算法。文献[8]在讨论代数曲线的分段逼近时,采用“单权因子”的标准形式的Bézier,提出了代数曲线的“肩点”到有理Bézier曲线的“肩点”的最小距离逼近算法,以解决代数曲线的参数化问题。考虑到传统的插值方法较少考虑函数的几何性质,同时也很少见到以参数方程给出的插值函数,这里将函数的有理B样条形式引入函数近似表示领域,给出一种基于有理二次B样条曲线的插值算法。

根据函数曲线的凹凸性和单调性,对函数进行分段,定义分段函数的“凸包”,并由“凸包”来确定控制多边形;重点介绍基于“肩点”等“关键点”信息的“双权因子”插值算法;相对于“单权因子”逼近算法[8],误差有明显改善;最后以数值实验说明了算法的有效性。

1 函数的分段

1.1 函数的几何特征分段

已知函数:

$\tilde{F}$:y=f(x) x∈[a0,an] (1)

设ai∈[a0,an],i=1,2,…n-1为函数拐点的横坐标,将定义域[a0,an]分割为:

[a0,an]=[a0,a1]∪[a1,a2]∪…

∪[an-2,an-1]∪[an-1,an] (2)

在区间:

[ai,ai+1] i=0,1,…,n-1 (3)

上,第i段函数表示为:

$\tilde{F}{}_i$:y=f(x) x∈[ai,ai+1] (4)

设分段函数式(4)在区间式(3)上的极值点为ami,将区间式(3)进一步分割(若已经单调,则不再分割)为如下形式:

[ai,ami]∪[ami,ai+1],i=0,1,…n-1 (5)

则分段函数式(4)进一步分割为如下两段:

y=f(x) x∈[ai,ami] x∈[ami,ai+1] (6)

这样,所给函数在每个分段区间上具有固定的凹凸性和单调性。

1.2 分段函数的三角形凸包

假定分割后的一段曲线为:

F:y=f(x) x∈[a,b] (7)

设它的两个端点为V0(xa,ya)、V2(xb,yb),求出两端点处切线,记为:

T0:a0x+b0y+c0=0 T2:a2x+b2y+c2=0 (8)

容易证明,切线T0、T2一定相交于一点,记交点为V1(x1,y1)。这样,由三点V0、V1、V2形成了分段函数曲线F的三角形凸包(如图1所示)。

例如,给定函数(如图2所示)。

C:y=x3+x2ex-sin(2πx)+1 x∈[-0.5,0.5] (9)

它有一个拐点(-0.0077,1.0484)和两个极值点(-0.2538,2.0333)、(0.2296,0.0866),将函数分割为四段,各段及其三角形凸包如图3所示。

2 函数分段二次插值与逼近

2.1 “单权因子”标准有理二次Bézier曲线逼近

设控制顶点为Qi(xi,yi),i=0,1,2,权因子为wi(>0),i=0,1,2,则有理二次Bézier曲线P(u),u∈[0,1]可以表示为:

$\{\begin{array}{l}x(u)=\frac{w{}_{0}x{}_0(1-u){}^2+2w{}_{1}x{}_{1}u(1-u)+w{}_{2}x{}_{2}u{}^2}{w{}_0(1-u){}^2+2w{}_{1}u(1-u)+w{}_{2}u{}^2}\\ y(u)=\frac{w{}_{0}y{}_0(1-u){}^2+2w{}_{1}y{}_{1}u(1-u)+w{}_{2}y{}_{2}u{}^2}{w{}_0(1-u){}^2+2w{}_{1}u(1-u)+w{}_{2}u{}^2}\end{array}$

按照Farin和Worsey提出的Bézier曲线的参数化标准形式[9],取w0=w2=1,为确定式(2)中余下的权因子w1,文献[8]中采用“肩点”到“肩点”距离最小逼近方法。具体方案是:已知凸包$V{}_i(\tilde{x}{}_i,\tilde{y}{}_i),i=0,1,2$;设F的“肩点”为SF(xF,yF),P(u)的“肩点”为:

$\begin{array}{l}S{}_{{\rm P}}={\rm P}(1/2)=S(w{}_1)=(x{}_{1/2},y{}_{1/2})\\ =(\frac{x{}_0+2w{}_{1}x{}_1+x{}_2}{2(1+w{}_1)},\frac{y{}_0+2w{}_{1}y{}_1+y{}_2}{2(1+w{}_1)})\end{array}$

记d2(SF,SP)=(xF-x1/2)2+(yF-y1/2)2,由$\frac{\partial d{}^2}{\partial w{}_1}=0$来确定权因子w1,结果为:

$w{}_1=\frac{1}{2}\cdot \frac{(x{}_0+x{}_2-2x{}_F)+\alpha (y{}_0+y{}_2-2y{}_F)}{(x{}_F-x{}_1)+\alpha (y{}_F-y{}_1)}$

其中,$x{}_i=\tilde{x}{}_i,y{}_i=\tilde{y}{}_i,i=0,1,2$;$\alpha =\frac{y{}_0+y{}_2-2y{}_1}{x{}_0+x{}_2-2x{}_1}$。

2.2 “双权因子”有理二次均匀B样条插值

2.2.1 有理二次均匀B样条曲线

设控制顶点为Qi(xi,yi),权因子为wi(>0),i=0,1,2;则有理二次均匀B样条曲线[4,5,6,7]P(u),u∈[0,1]可以表示为:

考虑到计算复杂性等因素,这里选取:

w0=1 (11)

并称由此对应的B样条曲线为“准标准形式”。

2.2.2 控制顶点的确定

已知F的三角形凸包V0V1V2,设$V{}_i=(\tilde{x}{}_i,\tilde{y}{}_i),i=0,1,2$;取控制顶点Q1=V1,则:

$x{}_1=\tilde{x}{}_{1}y{}_1=\tilde{y}{}_1$ (12)

由P(0)=V0,P(1)=V2得$\frac{Q{}_0+w{}_{1}Q{}_1}{1+w{}_1}=V{}_0,\frac{w{}_{1}Q{}_1+w{}_{2}Q{}_2}{w{}_1+w{}_2}=V{}_2$,即有:

x0=(1+w1)x～0-w1x～1

y0=(1+w1)y～0-w1y～1 (13)

x2=(w1+w2)x～2-w1x～1w2

y2=(w1+w2)y～2-w1y～1w2 (14)

这样,控制顶点将借助于下文的权因子w1、w2的确定而确定(如图4所示),曲线P(u)中的w1、w2的自由变化对应曲线段A的插值曲线族(如图5所示)。

2.3 有理二次B样条插值

选取F上的“关键点”,如“肩点”[9]SF(xF,yF)(如图6所示),它指的是F上到线段V0V1具有最大距离的点,“中点”${\rm M}(\frac{(a+b)}{2},f(\frac{(a+b)}{2}))$,“n等分点”等,同时选取B样条曲线P(u)上的“关键点”,如“肩点”SP(xP,yP)、参数域中点P(1/2)(记为S1/2)、参数域的“n等分点”等。这时,插值条件可以在“关键点”间的对等关系中选取一种,如SP=SF;S1/2=M等等,这里我们选取。

S1/2=SF (15)

考虑到式(10)、(11)两式有:

$\frac{x{}_0+6x{}_{1}w{}_1+x{}_{2}w{}_2}{1+6w{}_1+w{}_2}=x{}_F$;$\frac{y{}_0+6y{}_{1}w{}_1+y{}_{2}w{}_2}{1+6w{}_1+w{}_2}=y{}_F$ (16)

将式(13)、式(14)代入式(16)可以得到关于w1、w2的方程组,求出w1、w2,并由式(12)、式(13) 、式(14)确定控制顶点为Qi(xi,yi),i=0,1,2,最后由式(10)确定P(u)。

2.4 插值函数的性质

1) 保持原有函数的凹凸性不变;

2) 保持原有函数的单调性不变;

3) 在函数的分段点处保持G1连续性[6];

4) 插值曲线通过原始曲线的某个“关键点”。

3 分段插值算法

3.1 插值误差

在几何学中,曲线段C和它的近似表示Ca之间的误差经常使用Hausdorff距离[9]来描述:

$e(C,C{}_a)=dis(C,C{}_a)=\underset{{\rm P}\in C}{{\displaystyle \max }}\underset{Q\in C{}_a}{{\displaystyle \min }}d({\rm P},Q)$

这种距离便于理论分析,但不便于计算。

设F为分段待插函数,P(u)为分段插值函数,下面给出一种易于计算的误差概念。

定义1 在参数域[0,1]上的任意取一点$u{}_i=\frac{i}{n},0\le i\le n,n\in {\rm N}$,记P(u)上的对应点为(x(ui),y(ui)),F上的对应点为(x(ui),yi),定义插值误差为:

$e(F,{\rm P}(u))=\underset{i}{{\displaystyle \max }}|y{}_i-y(u{}_i)|$ (17)

3.2 算法描述

输入原始函数$\tilde{F}$以及误差限δ,输出误差小于δ的分段有理二次B样条曲线P(u)。

1) 由式(1)至式(6)对原始函数进行分段;

2) 计算分段函数曲线F两端切线T0,T2及其交点V1,确定三角形凸包V0V1V2;

3) 计算F的“肩点”SF以及P(u)上的点S1/2;

4) 根据式(12)、式(13)、式(14)、式(16)计算权因子w1,w2和控制顶点为Qi(xi,yi),i=0,1,2;

5) 依据式(10)计算插值曲线P(u);

6) 若插值误差e(F,P(u))<δ,则计算过程中止;否则,以“肩点” SF为界,将F一分为二,递归调用上述插值算法,直至满足给定的误差要求。

3.3 算法的收敛性

定理1 在算法的实施过程中,F和P(u)之间的Hausdorff距离收敛于0。即如果经过k次分割后的Hausdorff距离为Ek(F,P(u)),则:

$\underset{k\to +\infty }{{\displaystyle \lim }}E{}_k(F,{\rm P}(u))=0$ (18)

证明 如图7所示,由函数的分段方法及其凸包形成过程可知,ΔV0V1V2是确定的。设∠V1V0V2=α,∠V1V2V0=β,0<α,β≤π/2,记边V0V2的高为h,则h=|V0V1|sinα,h=|V2V1|sinβ,记“肩点”SF到边V0V2的距离为d。设经过SF的切线分别和V0V1、V2V1相交于P1、P2,自P1、P2分别向V0SF、V2SF引垂线,记P1H1=d11,P2H2=d12。

又设∠P1V0SF=α1,∠P2V2SF=β1,0<α1,β1<π/2。显然E0(F,P(u))<d<h,经过在“肩点”处一次分割后:

d11=|P1V0|sinα1=a1|V1V0|sinα1 0<a1<1

d12=|P2V2|sinβ1=b1|V1V2|sinβ1 0<b1<1

经过k次分割后:

$d{}_{k1}=a{}_{1}a{}_2\cdots a{}_k\left|V{}_{1}V{}_0\right|\sin \alpha {}_{k}0<a{}_i<1i=1,2‚\cdots ‚k$

$d{}_{k2}=b{}_{1}b{}_2\cdots b{}_k\left|V{}_{1}V{}_2\right|\sin \beta {}_{k}0<b{}_j<1j=1,2‚\cdots ‚k$

设$a=\underset{i}{{\displaystyle \max }}\{a{}_i\}$,$b=\underset{j}{{\displaystyle \max }}\{b{}_j\}$,显然0<a,b<1,于是:

$\underset{k\to +\infty }{{\displaystyle \lim }}d{}_{k1}<\underset{k\to +\infty }{{\displaystyle \lim }}a{}^k\left|V{}_{1}V{}_0\right|\sin \alpha {}_k=0$

$\underset{k\to +\infty }{{\displaystyle \lim }}d{}_{k2}<\underset{k\to +\infty }{{\displaystyle \lim }}b{}^k\left|V{}_{1}V{}_2\right|\sin \beta {}_k=0$

从而,$\underset{k\to +\infty }{{\displaystyle \lim }}E{}_k(F,{\rm P}(u))=0$。

4 数值实验

图8为函数C分四段(从左至右)“双权因子”的一步插值效果,图9为其误差函数。一步插值误差$e{}_1(F,{\rm P}(u))=\underset{i}{{\displaystyle \max }}\left|y{}_i-y(u{}_i)\right|<6\times 10{}^{-3}$。

图10为函数C分四段(从左至右)“单权因子”的一步逼近效果,图11为其误差函数。一步插值误差$e{}_1(F,{\rm P}(u))=\underset{i}{{\displaystyle \max }}\left|y{}_i-y(u{}_i)\right|<15\times 10{}^{-3}$。

显然,插值算法的误差不到逼近算法的一半,同时误差分布比较均匀,近似效果明显优于逼近算法。

5 结 语

基于函数的单调性和凹凸性分段的函数“双权因子”有理二次B样条插值算法,充实了插值方法的理论内涵,具有一定的实用价值。算法可以以任意精度近似表示已知函数。误差估计方法易于理解,计算简单。“双权因子”插值误差远小于“单权因子”逼近误差,而且误差分布比较均匀。但误差缩小的理论原因有待进一步探讨。

参考文献

[1] Wang Nengchao.Concise tutorial of computational approach[M].High Educational Press,2004.

[2] Wang Renhong.Numerical approximation[M].High Educational Press,1999.

[3] Yang Fengxiang,Zhai Ruicai, Sun Jing. Numerical Analysis [M]. Tinajin University Press, 1996.

[4] Zhong Xinxiong. Modeling Technology of Free Formed Curves and Surfaces[M]. Beijing, Science Press,2000.

[5] Shi Fazhong.Computer Aided Geometric Design an Non Uniform Rational Basic Spline[M]. Beijing,Higher Educational Press,2001.

[6] Les Piegl,Wayne Tiller. The NURBS Book [M]. New York: Springer,1997.

[7] Gerald Farin. Curves and Surfaces for Computer Aided Geometric Design [M].New York:ACADEMIC PRESS, 2000.

[8] Gao Xiao-Shan, Li Ming. Rational quadratic approximation to real algebraic curves[J].Computer Aided Geometric Design,2004,21:805-828.

广义B样条 篇5

关键词：逆向工程,轮廓约束点,B样条,曲面拟合

0引言

由截面测量数据进行曲面重构是科学计算可视化(visualization in scientific computing,VISC) 的重要内容之一,在医学可视化、无损探伤和逆向工程等领域应用广泛[1]。B样条具有凸包性、连续性与光滑性等诸多优点,故在曲面重构中被广泛采用[2,3]。对于数目庞大的测量数据点,若采用直接插值的方法获取B样条曲面,将导致节点和控制顶点数量过大,给后续的曲面光顺与修改、数据交换等造成负担[4,5]。因此,在满足数据轮廓精度要求的前提下,如何使曲面的数学模型包含尽量小的数据量,一直是该领域研究热点问题之一。

文献[6-11]对B样条曲面拟合中的相关技术,包括散乱点数据参数化和数据压缩等进行了研究。早期有代表性的插值方法[6]是先插值各截面数据点列构造B样条曲线,再通过节点插入获得可相容的曲线簇,最后对曲线簇的控制顶点放样得到插值曲面。这种方法较直观,但由于引入了节点插入操作,势必会引起控制顶点数量的膨胀[2]。为此,Piegl等[7]构造了具有非退化特性的参数区间,以增加节点选择的灵活性,使所得的控制顶点数目明显减少;在此基础上,Park[8]采用能量极小化方法进一步压缩了曲面模型的数据量。 王文珂等[10]借鉴传统蒙皮算法的思想,提出了一种由无序B样条曲线来拟合B样条曲面的算法, 可得到满足误差要求的光滑拟合曲面,但当曲线数目较多时该算法会遇到数值稳定性问题。综合上述情况,对B样条曲面拟合相关技术展开进一步研究仍然非常必要。

本文提出一种面向截面测量数 据的拟合 算法。首先采用弦高法剔除原始数据点中的坏点, 以得到较光滑的线形数据点列;然后求出离散曲率值,将超过平均曲率值的局部曲率极大值点定义为曲率优势点,并作为曲线轮廓的初始约束点, 得到初始的插值B样条曲线;再计算出偏差分布情况,并将超过允许误差的局部偏差极大值点增补为轮廓约束点,直到满足拟合精度要求;然后以约束点数最多的曲线为准,将其余曲线的约束点补齐以得到拓扑矩形的数据点阵;最后计算出公用的节点矢 量,获得最终 的B样条单向 蒙皮曲面。对一组截面测量数据的曲面拟合实验验证了算法的可行性。

1截面B样条曲线的获取

对于各截面测量数据,先用弦高法对截面数据进行消噪处理,然后根据点列的曲率分布情况, 确定曲率优势点并将优势点作为型值点,计算出各截面的初始B样条曲线。

1.1数据预处理

扫描数据点中一般会包含失真点,它们不能反映实际轮廓且对曲线光顺性存在较大影响,称作噪声点或坏点。在截面扫描数据点中,如果某个数据点偏离其相邻的两个点,且偏离程度超过预设值,则可视为坏点。关于剔除坏点的研究报道较多,针对截面测量数据的线形分布特点,本文采用文献[12]中的弦高差方法来去除坏点,对于考察点Pi,连接两个相邻点得到弦线Pi-1Pi+1,如果点Pi到该弦的距离大于允许误差,则将该点剔除。图1所示是某截面数据预处理效果。

下面以图1中的数据点为例,阐述获取截面B样条曲线的具体步骤。

1.2轮廓约束点的确定

1.2.1离散曲率计算

本文采用文献[13]给出的近似法,将数据点Pi处的曲率Ki近似为3个相邻点Pi-1、Pi和Pi+1所在圆弧的曲率。图1中数据点的曲率分布情况如图2所示。

1.2.2轮廓约束点的遴选

经预处理之后的密集数据点,其曲率分布能大体反映出轮廓走势。图1中的点列共包含202个数据点,为了压缩拟合曲线的数据量,现通过曲率值来确定初始轮廓约束点,以下简称约束点。

首先考虑曲率极大值法。即数据点Pi处的曲率值Ki,需同时满 足条件Ki> Ki-1和Ki>Ki+1。找出曲率极大值所对应的数据点,并将其作为约束点。则图1中的数据点包含67个约束点,其拟合情形如图3所示。由此得到数据点与拟合曲线间的最大偏差值为0.144mm,平均偏差值为0.051mm,均方差为0.028mm,偏差分布情况见图4。

从图3所示的约束点分布情况来看,该方法较有效地压缩了数据量,但即使在相对平直的局部区域,也分布了比较密集的约束点,显然没有必要。

再来考察曲率临界值法。事先设定曲率临界值Kcrt,将曲率值满足Ki> Kcrt的数据点作为约束点。不失一般性,此处将离散曲率的平均值Kavg设为临界值,即Kcrt= Kavg。对于图1中的数据点,用此方法确定的约束点为64个,所得拟合曲线如图5所示。最大拟合偏差值为0.969mm,平均偏差值为0.100mm,均方差为0.229mm,偏差分布情况见图6。

由图5可见,约束点可提取轮廓弯曲的局部区域,但由于该区域的离散曲率值一般都大于平均值,故约束点非常密集,而相对平直的区域约束点过少,这种分布显然也是不可取的。

上述两种确定约束点的方法,各有优劣。为利用二者的优点,本文选取超过平均曲率值的曲率极大值点作为约束点。即同时满足以下两个条件:1Ki> Ki-1且Ki> Ki+1;2Ki> Kavg。此处将符合上述条件的数据点称为曲率优势点,并将其作为轮廓约束点,对于开曲线,则需将首尾两个端点直接当作约束点。

采用曲率优势点法,则图1中的数据点只包含31个约束点,其拟合情形如图7所示。由此得到最大拟合偏差值为0.918mm,平均偏差值为0.164mm,均方差为0.261mm,偏差分布情况见图8。

与图3和图5相比,图7中的约束点数目得到了进一步的缩减,同时能在一定程度上保证数据点的轮廓精度。下面将插值于曲率优势点的曲线作为初始曲线。

诚然,由曲率优势点所确定的曲线,尽管大大压缩了数据量,但并不能保证拟合精度。若初始曲线存在拟合偏差超过预设值的局部区域,则需进行局部优化。

1.3拟合精度的优化

由图8可以看出,只有少量数据点处的拟合偏差值为零,因此偏差曲线一般会出现多个峰值, 记允许误差为[δ],若初始曲线的某些峰值大于 [δ],则直接将偏差峰值对应的数据点增设为约束点。也就是说,新的约束点将由两部分构成,即曲率优势点和少数偏差峰值数据点。

以图7中的初始 拟合曲线 为例 , 设 [ δ ] =0.5mm,对照图8可知,需要增加2个约束点,新曲线的拟合效果如图9所示。相应的拟合偏差最大值为0.235mm,平均偏差值为0.060mm,均方差为0.057mm,分布情况见图10。

可见,增加约束点之后,曲线局部区域的拟合精度得到了明显改善,而此时的轮廓约束点数为33个。该方法与曲率极值法和曲率临界值法相比,显然更加合理。当然,若偏差峰值仍大于[δ], 则需采用相同方法增加约束点,直到获得满足要求的B样条曲线。

2B样条曲面的构造

利用上述方法,可以获得满足精度要求的各截面B样条曲线,但是这些曲线是孤立的,缺乏相容性。各截面数据轮廓曲率变化各异,因而各节点矢量也不尽相同。由于B样条曲面具有张量积特性,需要解决曲线之间的相容性问题。为此,下面设法获取各截面曲线公用的节点矢量。

2.1公用节点矢量的获取

2.1.1约束点的参数化

为计算公用节点矢量,首先必须保证各曲线具有相同数量的约束点。前面得到的孤立B样条曲线簇,尽管已满足拟合精度要求,其拟合偏差曲线仍会存在多个峰值。现将各峰值点对应的数据点,按其偏差值由大到小的顺序存入各自数组, 称作备用约束点数组。

记第i条曲线的现有约束点数目为Ci,Ci的最大值记为Cmax。为确保各曲线具有相同的节点数目,现从各备用数组中依次取出必要数目的约束点,对于第i条曲线,需取出Cmax-Ci个数据点作为约束点。至此,所有曲线的约束点数目都为Cmax。

接下来实现各约束点的参数化。为尽量考虑型值点的分布对曲线形状的影响,此处采用平均弦长参数化方法。设共有k条截面曲线,对于每一条曲线,依次连接各约束点得到Cmax-1段弦线, 记第i条曲线的第j个弦长值为Li,j,则曲线簇第j段的弦长平均值Lj为

因此,用来计算节点矢量的总弦长d为

约束点的参数可以表示为

而其余参数

2.1.2节点矢量的计算

记曲线次数为p,m =Cmax+p。下面计算曲线簇公用的节点矢量U。其中,u0=u1= … =up=0,um-p=um-p+1= … =um=1,中间的节点采用平均值法,用下式求出:

这种方法得到的控制顶点与约束点的数目相等,避免了非奇异矩阵的产生,便于后续控制顶点的计算。

2.2控制网络的反算

下面根据文献 [ 14 ] 中给出的方法反求曲线的控制顶点P 。 通过各型值点参数和节点矢量U求出基函数)。 记n=Cmax-1 , 第k条截面曲线的第s个型值点为Qk , s, 则对应的控制顶点Pk , j可通过下式求解 :

图11所示的是B样条曲线簇的矩形拓扑控制网络示意图。

至此,曲面的节点矢量和控制顶点均已求出, 即B样条曲面已被定义。

3仿真实验与分析

为验证本文算法的可行性与有效性,下面以一组截面点云数据为例,用一张双三次B样条曲面对其进行拟合。图12所示的是经预处理之后的截面数据点,共21条线性点列,包含2557个数据点,分布在40mm×50mm×12mm的空间里, 平均每条点列包含约120个数据点。

按照上述步骤,通过遴选曲率优势点获得了初始拟合曲线,并以 [δ]=0.5mm为逼近偏差约束条件,对局部区域进行了精度优化,最终得到了V向B样条曲线簇(即各曲线的参数u为常数), 图13所示的是其中的一条曲线。单条曲线包含的控制顶点数目最多为22个,最少为13个。

为满足曲线簇的相容性要求,将所有曲线的轮廓约束点都增加至22个。然后运用平均弦长参数法,得到公用节点矢量为

所得的控制网络为21行×22列,即共有462个控制顶点,如图14所示。最终重构的B样条蒙皮曲面见图15。与原数据点数目相比,减少了2095个,即减少81.97%。

4结论

(1)通过曲率优势点法所确定的轮廓优势点, 可在提取截面数据基本形状的前提下,大幅缩减曲面模型的数据量,为后续的数据运算奠定基础。

(2)将初始拟合的偏差峰值对应点增设为B样条曲线型 值点,可快速有 效地提高 局部逼近 精度。

广义B样条 篇6

近年来,无人驾驶技术备受关注,各大研究机构和企业争相推出各自的无人驾驶平台。无人车作为未来智能交通的主要主体也逐渐融入到我们的日常生活中,比如自主巡航[1]和自动泊车等等。然而,为了使其更好地服务于我们,需要进一步提高其智能化水平,而路径规划作为连接环境感知和运动规划的桥梁,是无人车智能化水平的重要体现[2]。

由于受到自身动力学和运动学模型的约束,车辆的路径规划问题除过要严格满足端点状态约束之外,还要求其中间状态满足运动系统的微分约束。由于实现简单,并且高阶多项式曲线能够很好地满足运动系统的微分约束,生成高阶平滑的路径,所以很多路径规划系统选择使用基于多项式曲线的方法生成路径。B样条曲线是一种典型的多项式曲线,且因为其所有的中间状态均是由控制点加权生成,所以其能够完全满足端点状态约束。综合考虑无人车路径规划的要求和实现复杂度,在仅已知初始位姿和目标位姿的情况下,本文选择B样条曲线生成路径,重点讲述分步规划模型,即路径簇生成、最大曲率约束、碰撞检测以及最优评价四个过程,并通过Matlab仿真对本文方法进行了验证。

2 问题描述

本节分别描述了无人车路径规划问题和B样条曲线。

2.1 路径规划问题描述

路径规划得到的是一条从初始位置到目标位置的路径,即二维平面内一条从初始位置点到目标位置点的曲线,曲线上的每一个点表示车在行驶过程中的一个状态。考虑到实现方便,本文将路径描述成离散点序列{Sstart,S1,∙∙∙,Sn,Sgoal},如图1所示,序列中每一个点Si(xi,yi,θi)表示车的一个状态,其中(xi,yi)表示此时刻车辆的位置,θi表示车辆的航向,Sstart和Sgoal分别表示车辆的初始状态和目标状态。图1中的圆(xobs1,yobs1,robs1)表示环境中的障碍物,(xobs1,yobs1)表示障碍物的位置信息,robs1表示障碍物的半径。

2.2 B样条曲线

如果给定m+n+1个控制点Pi(i=0,1,∙∙∙,m+n),就可以构造m+1段n次B样条曲线,其可以表示为公式1:

其中,Pi,n(t)表示第i个n次B样条曲线片段,n表示B样条曲线的次数,t为控制参数,其取值范围为[0,1],Pi+k为控制点,Fk,n(t)为B样条基。依次连接全部n阶B样条曲线段就组成了整条B样条曲线。

3 本文算法

本节重点讲述基于B样条曲线的路径规划方法和基于该方法生成路径的过程。

3.1 基于B样条曲线的路径规划方法

选择三阶B样条曲线生成几何路径,即每四个控制点生成一个路径片段,然后通过片段的拼接就可以实现从初始状态到目标状态的路径规划,下面着重讲述基于六控制点的三阶B样条曲线生成满足车辆端点位姿约束路径的方法,如图2所示。

依据初始状态选择控制点P0,P1,P2。当P0,P1,P2三个控制点共线,并且P1为线段的中点时,生成的B样条曲线与线段相切于点P1。所以选择无人车的初始位置为控制点P1,将控制点P0和P2选在初始航向角θstart所在直线上,并关于控制点P1对称。如是,即能满足车辆的初始位姿约束;

依据目标状态选择控制点P3,P4,P5。当P3,P4,P5三个控制点共线,并且P4为线段的中点时,生成的B样条曲线与线段相切与点P4。所以选择无人车的目标位置为控制点P3,将控制点P3和P5选在目标航向角θgoal所在的直线上,并关于控制点P3对称。如是,即能满足车辆的目标位姿约束。

3.2 分步路径规划

本小节将以图3所给定的场景为例,讲述最优路径生成的过程。

3.2.1 路径簇生成

在选定控制点P1和P4之后,通过选择不同的控制点P2和P3,从而得到多组控制点,进而得到多条路径。将控制点选择的极限定为线段相等,但是不能出现交叉。将这个范围等间隔量化,各取四个点,共组成16组控制点,得到16条路径,如图3(a)中的蓝色曲线所示。

3.2.3 最大曲率约束

理论上,车辆的最小转弯半径是其本身属性[3],只取决于车辆的轴距L和最大前轮转角θmax。那么,车辆可行驶路径的最大曲率也是固定的,假设无人车可行驶路径的最大曲率以此为约束条件,在路径簇中选择满足最大曲率约束的路径,如图3(b)所示,图中绿色曲线表示不满足最大曲率约束的路径。

3.2.4 碰撞检测

碰撞检测的目的是保证车辆在规划路径上行驶而不与障碍物发生碰撞。采取的碰撞检测的方法很简单,因为环境中的障碍物采用圆来描述,所以只要判断路径上每一点到圆心的距离与障碍物半径的关系,就能确定其是否发生碰撞。由两点间距离公式:

如果d大于障碍物半径,则不发生碰撞;如果d小于障碍物半径,则发生碰撞。图3(c)中的蓝色曲线表示既满足最大曲率约束,又不与障碍物碰撞的路径。

3.2.2最优路径

路径要求的侧重点不同,优化的目标函数也可以有多种选择,常用的目标函数有最短和最平滑等。其中,路径最短可以抽象成优化问题:

路径最平滑可以抽象成优化问题:

式中,tr aoptimal为最优路径,为第i条路径的长度,为第i条路径上所有点处的曲率平方之和。图3(d)中的红色曲线即为得到的最短可行驶路径。

如是,就能得到满足车辆运动学约束,并且无碰撞的最优路径。

4 结论

本文选择使用B样条曲线解决无人车路径规划问题,并建立了基于B样条曲线的分步规划模型。仿真结果表明,使用基于B样条曲线的路径规划方法能够很好地解决简单障碍物场景中无人车的路径规划问题,并且因为路径生成过程简单,所以该方法常常表现得十分高效,能够完全满足无人车路径规划系统对算法实时性的要求。

参考文献

[1]Vahidi A,Eskandarian A.Research advances in intelligentcollision avoidance and adaptive cruise control[J].IEEETransactions on Intelligent Transportation Systems,2003,4(3):143-153.

[2]Siegwart R,Nourbakhsh I R,Scaramuzza D.Introduction toautonomous mobile robots[M].US:MIT press,2011.

一种自适应B样条插补算法研究 篇7

在复杂形状零件的高速加工中,采用B样条曲线插补指令插补的数控文件比通常采用的微小直线段插补的数控文件小10~100倍,而且避免了以直代曲,因而提高了工件加工精度,改善了表面质量。而且,B样条曲线插补刀轨在允许的加工方向变化范围内,无需频繁加减速,提高了加工速度。故研究B样条曲线插补器具有重要意义,是数控技术发展的需要。

研究人员在样条曲线插补方面已经做了一些工作,如Bedi[1]将参数增量设为一个恒定的常量,由此得出了一种均匀参数插补算法。然而采用这种算法时,在每一个相同的周期内走过的路径段长度大小不等,因而沿曲线的进给速度也不相等,存在着较大波动。为了对插补过程中沿曲线的进给速度进行控制,降低进给速度的波动,Huang和Yang[2]给出了基于欧拉和一阶泰勒展开式的近似插补算法,显示了比恒定参数插补优良的特性,但他们的设计没有加入误差控制机制。Yeh等[3]设计了一种自适应的参数曲线插补算法,保证进给速度在大部分加工过程中稳定,当加工误差超过了允许误差的范围时,进给速度作自适应调整,这种方法能提高加工速度,缺点是忽略了机床加减速能力对进给速度的影响。

本文在以上研究的基础上,设计了一种自适应B样条曲线插补器,能参照插补精度和机床加减速能力要求,对进给速度进行调整。同时算法的改进使插补计算速度更快,插补精度也得到提高。

2 插补点参数计算

B样条曲线的方程可以表示如下:

曲线插补的过程就是用时间序列{t1,t2,…,ti,…,tn}对参数变量进行分割以形成参数序列{u1,u2,…,ui,…,un},从而得到插补点序列{p"(u1),p"(u2),…,p"(ui),…,p"(un)}的过程[4]。从计算的角度看,参数曲线插补的关键就是如何根据插补点p"(ui)处进给速度vi计算出下一个插补点的参数ui+1,从而得到下一个插补点p"(ui+1)。

由插补过程可知曲线参数u为时间t的函数,即u=u(t),当t=ti+1时,由Taylor公式可得:

其中T=ti+1-ti,即插补周期,HOT为高次项,由于插补周期一般都很短(一般1~4ms),高次项几乎等于0。

由式(3)便可以由当前插补点的参数ui,速度v(ui)得到下一个插补点ui+1。

3 可控弦高误差的B样条曲线插补

利用一阶Taylor近似插补算法可以保证进给速度恒定,但是一直存在弦高误差,而当前曲线曲率半径过小且进给速度过大时,弦高误差有可能过大。为了使弦高误差限制在允许误差范围内,进给速度应该在插补过程中根据曲率半径的变化进行自动调整,这就是控制弦高误差的自动调节进给速度的插补[4]思想。通常用圆弧近似的方法确定弦高误差和进给速度及曲率半径之间的关系。

3.1 弦高误差计算

在插补过程中,由于每两个相邻的插补点之间的距离很小,所以可以将B样条曲线中两个插补点之间那段微曲线弧近似地看作一段圆弧。如图1所示,Li为两个插补点p(ui)和p(ui+1)之间的距离,!′为插补过程中的弦高误差,Ri为插补点p(ui)处的曲率半径,则弦高误差!′,可以用!代替,由图1可得:

由式(4)可知,在插补周期T一定的时候,vi越大,弦高误差越大。

3.2 B样条曲线曲率半径计算

设B样条曲线表示为p(u),则该曲线上参数u=ui点的曲率半径为:

根据德布尔算法,k次B样条曲线上一点处的r阶导矢(rp"()u)可按如下递推公式计算[5]:

其中

利用德布尔算法的递推式(6),即可以求出3次B样条曲线上数据点的一阶、二阶导矢,再通过式(5)可计算出该点的曲率半径Ri。

3.3 自适应可控弦高误差插补

为了保证参数曲线曲面的加工精度和质量,在切削加工过程中要求进给速度的变化比较平缓,同时保证加工弦高误差在允许范围内。由式(4)可知,在曲线插补时,插补精度分别与插补周期、曲率半径和进给速度有关。在给定插补周期的条件下,要得到尽可能高的插补精度,则要求进给速度尽可能的小。曲率半径越小的区域,要求进给速度也越小,这就是在加工到曲率半径小的区域必须减速的原因,否则就会超过加工要求的误差。

自适应可控弦高误差插补的目的就是插补过程中的进给速度能够根据当前曲线的曲率进行调整,从而保证零件的加工精度。

由式(4)可得:

所以当零件的加工精度ε一定的时候,则第i个插补点处的进给速度vi,必须满足下式

在调整插补进给速度的过程中,还必须要考虑机床的机械特性,以及机床的实际承受能力,比如加工过程中最大法向加速度,机床允许的最大进给速度等。

下面讨论加工过程中的进给速度和法向加速度之间的关系。设机床最大允许法向加速度为aNmax,则第i个插补点进给速度vi与法向加速度ai N间的关系由圆周运动可知

设机床最大允许进给速度为vmax,给合式(8)和式(10)可得:

本文采用T形加减速控制,在曲线开始插补时,设计曲线的起始速度为v0,然后以加速度a来加速,则第i个插补点处的插补进给速度vi=v0+i·a,当时,以后进给的插补点采用为插补进给速度。

3.4 前瞻加减速处理

可控弦高误差插补方法在理论上保证了曲线的轮廓精度,但是在计算过程中没有考虑机床允许的最大切向加速度aTmax的影响,在实际加工中并不能保证加工精度。

由图2可以得:

由式(7)可知,由于A、C点曲率半径较大,则速度较大;而B点曲率半径小,则速度较小。如果B样条曲线曲率变化过大,则很容易造成aB<-aTmax,aA>aTmax,这就很容易造成机床振动,影响加工精度。要避免这种情况,可以设计一个实时前瞻处理模块,主要有两个功能:检测速度变化敏感区、进行加减速处理。设机床的最大加减速能力为aTmax,如果机床遭遇小曲率半径的区域需要加减速处理,存在两种情况:一种是减速处理,一种是加速处理[4]。

在减速的情况下,出现减速度超过机床的加减速能力,则由当前插补点i依次向前一个插补点i-1进行速度调节,直到减速度不超过机床的最大减速能力为止。

考虑极限情况,也就是进给速度由机床的最大允许速度vmax减小到0的情况。此时需要的最小减速处理距离Smin=v2max/2aTmax,的情况。因此要保证插补过程中在曲率较大处插补点位置上的进给速度能及时的降下来,必须预先分析n个插补周期,而这n个插补周期的插补长度之和要不小于Smin,只有这样才能保证以后插补进给速度的调整不影响当前插补点处的进给速度,并且插补进给速度在适当的时候能及时的降下来。

在加速阶段,如果有(vi+1-vi)/T>aTmax,则调整速度为

在处理的过程中,需要构建一个已处理的插补点序列结构数组P[i]{u,v,s,Pos},其中:u和v分别为当前点参数和速度;s为当前插补点到下一个插补点的距离,P[i].s=v×T;Pos为当前点坐标。存入新插补点过程中还需要用变量Length记录数组P中第一个插补点与最后一个插补点的插补长度总和。

4 插补仿真实验

4.1 B样条曲线实例

采用三次B样条曲线对本文设计的B样条插补器进行插补仿真计算,使用Matlab2006开发环境进行编程。曲线特征和插补参数为:

控制点{(0,15),(3.2,4.6),(10.4,-4.8),(19.6,17.8),(31.1,13.1),(41.4,1.1),(50,40.9),(61.6,-8),(69.4,27.2),(76.1,19),(80,8)};节点矢量{0,0,0,0,0.13,0.28,0.34,0.41,0.57,0.74,0.87,1,1,1,1};机床最大允许速度为100mm/s,最大允许切向和法向加速度为1250mm/s2,插补周期为4ms,最大允许弦高误差为0.001mm。根据前文所介绍算法生成该曲线插补轮廓如图3所示。

4.2 仿真结果

插补进给速度比较如图4所示,其中虚线是自适应可控弦高误差插补进给速度曲线,实线是经加减速前瞻调整后的进给速度曲线,由图4可看出调整后的速度曲线不会大于弦高误差插补计算出的速度,且速度变化均匀,因此能保证良好的加工精度。由图5可见,加减速前瞻处理以前,在小曲率半径区域,最大的加减速度可以达到2800mm/s2,加入前瞻处理功能,超过规定的加速度(1250mm/s2)的区域,进给加速度调整,加速度处理的区域变宽,加速度的轮廓曲线得到改善。

5 结语

本文成功地设计了一个实时的自适应B样条插补算法,通过一个自适应速度调整和实时前瞻的加减速处理,保证了一个比较稳定的进给速度,同时保证了加工的精度误差,可以满足机床的加减速反应能力。本插补算法将应用在笔者实验室开发的基于DSP的运动控制器中。

摘要：针对复杂轮廓零件数控加工的实际要求,设计了一种在可控弦高误差插补算法基础上加入前瞻功能的自适应非均匀B样条插补算法。该方法能够根据曲线的形状,自适应地调整进给速度,保持速度平稳,并通过前瞻模块在速度变化敏感区域对加减速进行处理,满足机床加减速要求。

关键词：三次B样条曲线,插补器

参考文献

[1]S.Bedi,I.Ali,N.Quan.Advanced techniques for CNC machines[J].Trans.Of the ASME J.of Engineering for Industry,1993(115):329-336.

[2]J.T.Huang,D.C.H.Yang.A generalized interpolator for command generation of parametric curves in computer controlled machines[J].ASME Japan/USA Symposium on Flexible Automation,1992(1):393-399.

[3]YEH S S,HSU P L.Adaptive-feedrate interpolation for parametric curves with a confined chord error[J].Computer-Aided Design,2002,34(3):229-237.

[4]康书杰.NURBS曲线实时插补技术研究[D].南京:南京航空航天大学,2007.