基于移动最小二乘法的曲线曲面拟合_曾清红
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
3.1 曲线拟合 使 用 移 动 最 小 二 乘 (MLS) 法 对 离 散 数 据 点 x=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]， y=[0, 4, 5, 14, 15, 14.5, 14, 12, 10, 5, 4]进行曲线拟 合，拟合过程采用线性基函数 p(x)=[1, x]T，图 3 给 出了拟合结果。图 3(a)~(c)三组图分别给出了权函 数为全局定常、紧支定常和紧支光滑时，MLS 拟 合所得到的结果曲线。 （1）图 3（a）使用的权函数在整个区域内都 是常数， 这时移动最小二乘法相当于传统的最小二 乘法，其行为相当于是一个线性拟合。
收稿日期：2002－10－21 基金项目：国家自然科学基金资助项目（10102020） ；国家 973 资助项目（G1999032805） 作者简介：曾清红（1978-） ，男，湖南长沙人，硕士生，主要研究领域为大规模科学计算及科学计算可视化。
第 １ 期 曾清红等：基于移动最小二乘法的曲线曲面拟合 ・８５・
其中
（4） （5）
示，图 2 给出了三次样条函数曲线
2 3 2 − 4s + 4s 3 2 4 3 4 w( s ) = − 4s + 4s − s 3 3 0
A( x ) = ∑ w( x − x I ) p( x I ) p T ( x I )
I =1
n
（6） （7） （8）
0.4
拟合曲线 拟合 曲 线
10 8 6 4 2 0
0.3
0.2
0.1
0.0 0.0 0.2 0.4 0.6 0.8 1.0
0.0
0.2
0.4
0.6
0.8
1.0
x
x
来自百度文库
（c）紧支光滑权函数 图3 权函数与拟合曲线
3.2 曲面拟合 假设区域Ω :{－3.0≤x≤3.0, －3.0≤y≤3.0}上的 函数 f ( x, y ) 定义为
B( x ) = [ w( x − x1 ) p ( x1 ), w( x − x2 ) p( x2 ), Λ , w( x − xn ) p( xn )] y = [ y1 , y 2 , Λ , y n ]
T
1 (s ≤ ) 2 1 ( < s ≤ 1) （12） 2 ( s > 1)
影响区域应该包含足够多的节点，以使得式（6） 中 A(x)可逆。如果式（1）使用的是线性基函数， 则曲线拟合的影响区域内应至少包含不重叠的两 个节点， 曲面拟合的影响区域内应至少包含不在同 一直线上的 3 个节点。
第 １ 期 曾清红等：基于移动最小二乘法的曲线曲面拟合 ・８７・
2.0
16 14
1.5
12
拟合曲线 拟合 曲 线
0.0 0.2 0.4 0.6 0.8 1.0
权函数 权函数
10 8 6 4 2
1.0
0.5
0.0 0 0.0 0.2 0.4 0.6 0.8 1.0
x 的 函 数 。
待 求 系 数 ， 它 是 坐 标
p( x ) = [ p1 ( x ), p2 ( x ),Λ , pm ( x )]T 称为基函数， 它
是一个 k 阶完备的多项式，m 是基函数的项数。例 如，对于二维问题， 线性基
需要注意的是即使基函数 p(x)为多项式，式 （ 9 ）中的 f (x) 也不再是多项式。如果基函数
x
smax
权函数
X
图1
影响区域
图2
三次样条权函数
2 移动最小二乘拟合流程
使用移动最小二乘法进行曲线曲面拟合的基 本思想是先将拟合区域网格化，然后利用公式（9）
表1
求出网格点上节点值， 最后连接网格节点形成拟合 曲线(曲面)。表 1 给出了使用移动最小二乘法进行 曲线曲面拟合的程序设计流程。
基于移动最小二乘法的曲线曲面拟合程序设计流程
x − xI
2
的增加单调递减。权函数还应具有一定
的光滑性，因为拟合函数会继承权函数的连续性： 如果权函数 w(x－xI)是 C1 阶连续的，则拟合函数也 是 C 阶连续的。常用的权函数是样条函数，记 s = x－xI，s =
1
s smax
则三次样条权函数如式（12）所
∂J = A( x )α ( x ) − B( x ) y = 0 ∂α α ( x ) = A −1 ( x ) B ( x ) y
２００４ 年 工 程 图 学 学 报 ２００４ 第 １ 期 JOURNAL OF ENGINEERING GRAPHICS No.1
基于移动最小二乘法的曲线曲面拟合
p ∈ C r ， 权 函 数 w∈Cs ， f ∈ C min( r , s ) 。
1.2 权函数
则拟合函数
p( x ) = [1, x, y ]T , m = 3 （2a） 2 2 T 二次基 p( x ) = [1, x, y , x , xy , y ] , m = 6
（2b） 范 式 ( 向 量
考 虑 下 面 的 加 权 离 散 L2
x = [ x1 , x 2 Λ , xn ] 的 L2 范式 x 2 = ∑ i =1 xi J = ∑ w( x − x I )[ f ( x ) − y I ]2 =
I =1 n
)
∑ w( x − x I )[ p T ( x I )α ( x ) − y I ]2
f ( x ) = ∑ α i ( x ) pi ( x ) = p T ( x )α ( x )
i =1
m
（1）
Ö IShepard ( x ) =
∑
w( x − x I ) w( x = x J ) J =1
n
（11）
式中
α ( x ) = [ a 1 ( x ), a 2 ( x ), Λ , a m ( x )] T 为
其中
Ö k ( x ) 称为形函数， k 表示基函数的阶数 Ö k ( x) = [ Φ1k ,Φ2k , Λ ,Φnk ] = p T ( x ) A−1 ( x ) B( x )
（10）
1 移动最小二乘法
1.1 拟合函数的建立 在拟合区域的一个局部子域上， 拟合函数 f (x) 表示为
如果 k = 0 ，则基函数 p ( x ) = {1} ，这时的形函数 为 Shepard 函数
曾清红， 卢德唐
（中国科学技术大学力学和机械工程系，合肥 230027）
摘 要：建立了一种基于移动最小二乘(Moving Least-Squares， MLS)法的曲线曲 面拟合方法。这种方法对传统的最小二乘(LS)法的作了比较大的改进，使生成的曲线曲面具 有精度高、光滑性好等许多优点。详细介绍了移动最小二乘法的原理、应用和特点，并且给 出了使用移动最小二乘法进行曲线曲面拟合的程序设计流程。 最后给出了曲线拟合和空间散 乱数据曲面拟合算例，将拟合结果与最小二乘拟合结果作了比较，分析了 MLS 拟合曲线曲 面的光滑性和拟合质量，表明了该方法的优越性和有效性。 关 键 词：计算机应用；移动最小二乘法；曲线曲面拟合；权函数 中图分类号：TP 391 文献标识码：A 文 章 编 号：1003-0158(2004)01-0084-06
（2）图 3（b）中权函数在一个小的子域内等 于常数，在其它地方等于零。虽然这个权函数也具 有紧支性，但是不光滑，它导致的拟合曲线就是将 所有的数据点连线，即相当于分段线性插值。 （3）图 3（c）中权函数是三次样条函数，不 但具有紧支特性，而且光滑。即使使用的是线性基 函数，但由于继承了权函数的连续性，拟合曲线仍 然具有很好的光滑性。 图 3 说明了 MLS 中权函数的重要性。在最后 一组图中，仅使用线性基函数，就获得很好的拟合 结果， 这是传统最小二乘法和其它拟合方法无法做 到的。
（1）将拟合区域网格化； （2）对每个网格点 x 进行循环： (a) 确定网格点 x 的影响区域的大小； (b) 确定包含在 x 的影响区域内的节点； (c) 计算形函数 Ö k ( x ) ； (d) 计算网格点 x 处的节点值。 （3）结束网格点循环； （4）连接网格点形成拟合曲线(曲面)。
3 算 例
f ( x , y ) = 2(1 − x ) 2 exp[ − x 2 − ( y + 1) 2 ] − x 10 ( x 3 − y 5 ) exp[ − x 2 − y 2 ] − （13） 5 1 exp[ − ( x + 1) 2 − y 2 ] 3
Belytschko 将其应用于无网格方法中[4]。移动最小 二乘法与传统的最小二乘法相比， 有两个比较大的 改进： （1）拟合函数的建立不同。这种方法建立拟 合函数不是采用传统的多项式或其它函数， 而是由 一个系数向量 a(x)和基函数 p(x)构成，这里 a(x)不 是常数，而是坐标 x 的函数。 （2）引入紧支（Compact Support）概念，认 为点 x 处的值 y 只受 x 附近子域内节点影响， 这个 子域称作点 x 的影响区域，影响区域外的节点对 x 的取值没有影响。在影响区域上定义一个权函数 w(x)，如果权函数在整个区域取为常数，就得到传 统的最小二乘法。 这些改进能够带来许多优点， 减缓或解决传统 曲线曲面拟合过程中存在的困难。 可以取不同阶的 基函数以获得不同的精度， 取不同的权函数以改变 拟合曲线(曲面)的光滑度，这是其它拟合方法无法 做到的。
将式(5)代入式(1)， 就可以得到 MLS 拟合函数
f ( x ) = ∑ΦIk ( x ) y I = Ö k ( x ) y
I =1
n
（9）
・８６・ 工 程 图 学 学 报 ２００４ 年
x
x
（a）全局定常权函数
1.0
16 14
0.8
12
拟合曲线 拟合 曲 线
0.0 0.2 0.4 0.6 0.8 1.0
权函数 权函数
10 8 6 4
0.6
0.4
0.2 2 0.0 0 0.0 0.2 0.4 0.6 0.8 1.0
x
x
（b）紧支定常权函数
0.7
16 14 12
0.6
0.5
权函数 权函数
n 1/ 2 2
权函数在移动最小二乘法中起着非常重要的 作用。移动最小二乘法中的权函数 w(x－xI)应该具 有紧支性， 也就是权函数在 x 的一个子域内不等于 零，在这个子域之外全为零，这个子域称为权函数 的支持域(即 x 的影响区域)。一般选择圆形作为权 函数的支持域(见图 1)，其半径记为 smax。 由于权 函数的紧支性， 只有这些包含在影响区域内的数据 点对点 x 的取值有影响。 权函数 w(x－xI)应该是非负的，并且随着
I =1
n
（3）
n 是影响区域内节点的数目， f ( x ) 是拟合 函数， yI 是 x = x I 处的节点值， yI = y ( x I ) ； w( x − xI ) 是 节 点 x I 的 权 函 数 。 为 确 定 系 数 α ( x ) ，式（3）应该取极小值。式（3）对 α 求导
式中 得
曲线曲面拟合是一种古老而常用的技术， 在工 程、实验、统计和计算机图形等方面有着广泛的应 用。通常，经过测量或者采集可以得到一组离散数 据点(xi, yi), i=1, 2, …, n， 这里 xi 为坐标值。 由于这 些数据点并非完全精确， 而且函数 y= f (x)的表达式 预先无法知道，需要在给定的函数类 f 上根据这些 离散数据作出逼近曲线(曲面)。因为离散数据有误 差，并不要求逼近曲线(曲面)经过数据点，而只是 要求逼近曲线(曲面) f (x)的误差的某个指标达到最 小。传统的曲线(曲面)拟合方法一般使用最小二乘 法[1], [2]，通过使误差的平方和最小，得到一个线性 方程组，求解线性方程组就可以得到拟合曲线(曲 面)。如果离散数据量比较大、形状复杂，还需要 进行分段(分块)拟合和平滑化，这在实际中往往带 来一定的困难。 笔者使用移动最小二乘 (MLS) 法建立一种新 的曲线(曲面)拟合方法，这种方法能够克服以上困 难 ， 并 且 还 具 有 许 多 其 它 优 点 。 Lancaster 和 Salkauskas 最先在曲面生成中使用了 MLS[3]，后来
收稿日期：2002－10－21 基金项目：国家自然科学基金资助项目（10102020） ；国家 973 资助项目（G1999032805） 作者简介：曾清红（1978-） ，男，湖南长沙人，硕士生，主要研究领域为大规模科学计算及科学计算可视化。
第 １ 期 曾清红等：基于移动最小二乘法的曲线曲面拟合 ・８５・
其中
（4） （5）
示，图 2 给出了三次样条函数曲线
2 3 2 − 4s + 4s 3 2 4 3 4 w( s ) = − 4s + 4s − s 3 3 0
A( x ) = ∑ w( x − x I ) p( x I ) p T ( x I )
I =1
n
（6） （7） （8）
0.4
拟合曲线 拟合 曲 线
10 8 6 4 2 0
0.3
0.2
0.1
0.0 0.0 0.2 0.4 0.6 0.8 1.0
0.0
0.2
0.4
0.6
0.8
1.0
x
x
来自百度文库
（c）紧支光滑权函数 图3 权函数与拟合曲线
3.2 曲面拟合 假设区域Ω :{－3.0≤x≤3.0, －3.0≤y≤3.0}上的 函数 f ( x, y ) 定义为
B( x ) = [ w( x − x1 ) p ( x1 ), w( x − x2 ) p( x2 ), Λ , w( x − xn ) p( xn )] y = [ y1 , y 2 , Λ , y n ]
T
1 (s ≤ ) 2 1 ( < s ≤ 1) （12） 2 ( s > 1)
影响区域应该包含足够多的节点，以使得式（6） 中 A(x)可逆。如果式（1）使用的是线性基函数， 则曲线拟合的影响区域内应至少包含不重叠的两 个节点， 曲面拟合的影响区域内应至少包含不在同 一直线上的 3 个节点。
第 １ 期 曾清红等：基于移动最小二乘法的曲线曲面拟合 ・８７・
2.0
16 14
1.5
12
拟合曲线 拟合 曲 线
0.0 0.2 0.4 0.6 0.8 1.0
权函数 权函数
10 8 6 4 2
1.0
0.5
0.0 0 0.0 0.2 0.4 0.6 0.8 1.0
x 的 函 数 。
待 求 系 数 ， 它 是 坐 标
p( x ) = [ p1 ( x ), p2 ( x ),Λ , pm ( x )]T 称为基函数， 它
是一个 k 阶完备的多项式，m 是基函数的项数。例 如，对于二维问题， 线性基
需要注意的是即使基函数 p(x)为多项式，式 （ 9 ）中的 f (x) 也不再是多项式。如果基函数
x
smax
权函数
X
图1
影响区域
图2
三次样条权函数
2 移动最小二乘拟合流程
使用移动最小二乘法进行曲线曲面拟合的基 本思想是先将拟合区域网格化，然后利用公式（9）
表1
求出网格点上节点值， 最后连接网格节点形成拟合 曲线(曲面)。表 1 给出了使用移动最小二乘法进行 曲线曲面拟合的程序设计流程。
基于移动最小二乘法的曲线曲面拟合程序设计流程
x − xI
2
的增加单调递减。权函数还应具有一定
的光滑性，因为拟合函数会继承权函数的连续性： 如果权函数 w(x－xI)是 C1 阶连续的，则拟合函数也 是 C 阶连续的。常用的权函数是样条函数，记 s = x－xI，s =
1
s smax
则三次样条权函数如式（12）所
∂J = A( x )α ( x ) − B( x ) y = 0 ∂α α ( x ) = A −1 ( x ) B ( x ) y
２００４ 年 工 程 图 学 学 报 ２００４ 第 １ 期 JOURNAL OF ENGINEERING GRAPHICS No.1
基于移动最小二乘法的曲线曲面拟合
p ∈ C r ， 权 函 数 w∈Cs ， f ∈ C min( r , s ) 。
1.2 权函数
则拟合函数
p( x ) = [1, x, y ]T , m = 3 （2a） 2 2 T 二次基 p( x ) = [1, x, y , x , xy , y ] , m = 6
（2b） 范 式 ( 向 量
考 虑 下 面 的 加 权 离 散 L2
x = [ x1 , x 2 Λ , xn ] 的 L2 范式 x 2 = ∑ i =1 xi J = ∑ w( x − x I )[ f ( x ) − y I ]2 =
I =1 n
)
∑ w( x − x I )[ p T ( x I )α ( x ) − y I ]2
f ( x ) = ∑ α i ( x ) pi ( x ) = p T ( x )α ( x )
i =1
m
（1）
Ö IShepard ( x ) =
∑
w( x − x I ) w( x = x J ) J =1
n
（11）
式中
α ( x ) = [ a 1 ( x ), a 2 ( x ), Λ , a m ( x )] T 为
其中
Ö k ( x ) 称为形函数， k 表示基函数的阶数 Ö k ( x) = [ Φ1k ,Φ2k , Λ ,Φnk ] = p T ( x ) A−1 ( x ) B( x )
（10）
1 移动最小二乘法
1.1 拟合函数的建立 在拟合区域的一个局部子域上， 拟合函数 f (x) 表示为
如果 k = 0 ，则基函数 p ( x ) = {1} ，这时的形函数 为 Shepard 函数
曾清红， 卢德唐
（中国科学技术大学力学和机械工程系，合肥 230027）
摘 要：建立了一种基于移动最小二乘(Moving Least-Squares， MLS)法的曲线曲 面拟合方法。这种方法对传统的最小二乘(LS)法的作了比较大的改进，使生成的曲线曲面具 有精度高、光滑性好等许多优点。详细介绍了移动最小二乘法的原理、应用和特点，并且给 出了使用移动最小二乘法进行曲线曲面拟合的程序设计流程。 最后给出了曲线拟合和空间散 乱数据曲面拟合算例，将拟合结果与最小二乘拟合结果作了比较，分析了 MLS 拟合曲线曲 面的光滑性和拟合质量，表明了该方法的优越性和有效性。 关 键 词：计算机应用；移动最小二乘法；曲线曲面拟合；权函数 中图分类号：TP 391 文献标识码：A 文 章 编 号：1003-0158(2004)01-0084-06
（2）图 3（b）中权函数在一个小的子域内等 于常数，在其它地方等于零。虽然这个权函数也具 有紧支性，但是不光滑，它导致的拟合曲线就是将 所有的数据点连线，即相当于分段线性插值。 （3）图 3（c）中权函数是三次样条函数，不 但具有紧支特性，而且光滑。即使使用的是线性基 函数，但由于继承了权函数的连续性，拟合曲线仍 然具有很好的光滑性。 图 3 说明了 MLS 中权函数的重要性。在最后 一组图中，仅使用线性基函数，就获得很好的拟合 结果， 这是传统最小二乘法和其它拟合方法无法做 到的。
（1）将拟合区域网格化； （2）对每个网格点 x 进行循环： (a) 确定网格点 x 的影响区域的大小； (b) 确定包含在 x 的影响区域内的节点； (c) 计算形函数 Ö k ( x ) ； (d) 计算网格点 x 处的节点值。 （3）结束网格点循环； （4）连接网格点形成拟合曲线(曲面)。
3 算 例
f ( x , y ) = 2(1 − x ) 2 exp[ − x 2 − ( y + 1) 2 ] − x 10 ( x 3 − y 5 ) exp[ − x 2 − y 2 ] − （13） 5 1 exp[ − ( x + 1) 2 − y 2 ] 3
Belytschko 将其应用于无网格方法中[4]。移动最小 二乘法与传统的最小二乘法相比， 有两个比较大的 改进： （1）拟合函数的建立不同。这种方法建立拟 合函数不是采用传统的多项式或其它函数， 而是由 一个系数向量 a(x)和基函数 p(x)构成，这里 a(x)不 是常数，而是坐标 x 的函数。 （2）引入紧支（Compact Support）概念，认 为点 x 处的值 y 只受 x 附近子域内节点影响， 这个 子域称作点 x 的影响区域，影响区域外的节点对 x 的取值没有影响。在影响区域上定义一个权函数 w(x)，如果权函数在整个区域取为常数，就得到传 统的最小二乘法。 这些改进能够带来许多优点， 减缓或解决传统 曲线曲面拟合过程中存在的困难。 可以取不同阶的 基函数以获得不同的精度， 取不同的权函数以改变 拟合曲线(曲面)的光滑度，这是其它拟合方法无法 做到的。
将式(5)代入式(1)， 就可以得到 MLS 拟合函数
f ( x ) = ∑ΦIk ( x ) y I = Ö k ( x ) y
I =1
n
（9）
・８６・ 工 程 图 学 学 报 ２００４ 年
x
x
（a）全局定常权函数
1.0
16 14
0.8
12
拟合曲线 拟合 曲 线
0.0 0.2 0.4 0.6 0.8 1.0
权函数 权函数
10 8 6 4
0.6
0.4
0.2 2 0.0 0 0.0 0.2 0.4 0.6 0.8 1.0
x
x
（b）紧支定常权函数
0.7
16 14 12
0.6
0.5
权函数 权函数
n 1/ 2 2
权函数在移动最小二乘法中起着非常重要的 作用。移动最小二乘法中的权函数 w(x－xI)应该具 有紧支性， 也就是权函数在 x 的一个子域内不等于 零，在这个子域之外全为零，这个子域称为权函数 的支持域(即 x 的影响区域)。一般选择圆形作为权 函数的支持域(见图 1)，其半径记为 smax。 由于权 函数的紧支性， 只有这些包含在影响区域内的数据 点对点 x 的取值有影响。 权函数 w(x－xI)应该是非负的，并且随着
I =1
n
（3）
n 是影响区域内节点的数目， f ( x ) 是拟合 函数， yI 是 x = x I 处的节点值， yI = y ( x I ) ； w( x − xI ) 是 节 点 x I 的 权 函 数 。 为 确 定 系 数 α ( x ) ，式（3）应该取极小值。式（3）对 α 求导
式中 得
曲线曲面拟合是一种古老而常用的技术， 在工 程、实验、统计和计算机图形等方面有着广泛的应 用。通常，经过测量或者采集可以得到一组离散数 据点(xi, yi), i=1, 2, …, n， 这里 xi 为坐标值。 由于这 些数据点并非完全精确， 而且函数 y= f (x)的表达式 预先无法知道，需要在给定的函数类 f 上根据这些 离散数据作出逼近曲线(曲面)。因为离散数据有误 差，并不要求逼近曲线(曲面)经过数据点，而只是 要求逼近曲线(曲面) f (x)的误差的某个指标达到最 小。传统的曲线(曲面)拟合方法一般使用最小二乘 法[1], [2]，通过使误差的平方和最小，得到一个线性 方程组，求解线性方程组就可以得到拟合曲线(曲 面)。如果离散数据量比较大、形状复杂，还需要 进行分段(分块)拟合和平滑化，这在实际中往往带 来一定的困难。 笔者使用移动最小二乘 (MLS) 法建立一种新 的曲线(曲面)拟合方法，这种方法能够克服以上困 难 ， 并 且 还 具 有 许 多 其 它 优 点 。 Lancaster 和 Salkauskas 最先在曲面生成中使用了 MLS[3]，后来