基于移动最小二乘法的曲线曲面拟合论文阅读笔记
论文地址:http://www.cnki.com.cn/Article/CJFDTotal-GCTX200401016.htm
一、Problem Statement
传统的曲线(曲面)拟合方法一般使用最小二乘法, 通过使误差的平方和最小, 得到一个线性方程组,求解线性方程组就可以得到拟合曲线(曲面)。如果离散数据量比较大、形状复杂,还需要进行分段(分块)拟合和平滑化,这在实际中往往带来一定的困难。
二、Direction
建立了了一种基于移动最小二乘(Moving Least-Squares, MLS)法的曲线曲面拟合方法。
- 拟合函数的建立不同。这种方法建立拟合函数不是采用传统的多项式或其他函数,而是由一个系数向量 α ( x ) \alpha(x) α(x)和基函数 p ( x ) p(x) p(x)构成,这里 α ( x ) \alpha(x) α(x)不是常数,而是坐标 x x x的函数。
- 引入紧支(Compact Support)概念,认为点 x x x处的值 y y y只受 x x x附近子域内节点影响,这个子域称作点 x x x的影响区域。影响区域上定义一个权函数 w ( x ) w(x) w(x),如果权函数在这个区域内取为常数,就得到传统的最小二乘法。
三、Method
使用移动最小二乘法进行曲线曲面拟合的基本思想是先将拟合区域网格化,然后求出网格点上节点值,最后连接网格节点形成拟合曲线(曲面)。
1. 拟合函数的建立
在拟合区域的一个局部子域上,拟合函数 f ( x ) f(x) f(x)表示为:
f ( x ) = ∑ i = 1 m α i ( x ) p i ( x ) = p T ( x ) α ( x ) \begin{equation}f(x) = \sum_{i=1}^m \alpha_i(x)p_i(x)= p^T(x)\alpha(x) \end{equation} f(x)=i=1∑mαi(x)pi(x)=pT(x)α(x)
式中, α ( x ) = [ α 1 ( x ) , α 2 ( x ) , Λ , a m ( x ) ] T \alpha(x)=[\alpha_1(x), \alpha_2(x), \Lambda, a_m(x)]^T α(x)=[α1(x),α2(x),Λ,am(x)]T为待求系数,它是坐标 x x x的函数。 p ( x ) = [ p 1 ( x ) , p 2 ( x ) , Λ , p m ( x ) ] T p(x)=[p_1(x), p_2(x), \Lambda, p_m(x)]^T p(x)=[p1(x),p2(x),Λ,pm(x)]T称为基函数,它是一个 k k k阶完备的多项式, m m m是基函数的项数。例如:
空间维度 | 一阶形式 | 二阶形式 |
---|---|---|
1维: x = ( x ) x=(x) x=(x) | p 1 = 1 , p 2 = x p_1 = 1, p_2 = x p1=1,p2=x | p 1 = 1 , p 2 = x , p 3 = x 2 p_1 = 1, p_2 = x, p_3=x^2 p1=1,p2=x,p3=x2 |
2维: x = ( x , y ) x=(x,y) x=(x,y) | p 1 = 1 , p 2 = x , p 3 = y p_1 = 1, p_2 = x, p_3=y p1=1,p2=x,p3=y | p 1 = 1 , p 2 = x , p 3 = y , p 4 = x y , p 5 = x 2 , p 6 = y 2 p_1 = 1, p_2 = x, p_3=y, \\p_4=xy, p_5=x^2, p_6=y^2 p1=1,p2=x,p3=y,p4=xy,p5=x2,p6=y2 |
3维: x = ( x , y , z ) x=(x,y,z) x=(x,y,z) | p 1 = 1 , p 2 = x , p 3 = y , p 4 = z p_1 = 1, p_2 = x, p_3=y, p_4=z p1=1,p2=x,p3=y,p4=z | p 1 = 1 , p 2 = x , p 3 = y , p 4 = z , p 5 = x y , p 6 = x z , p 7 = y z p 8 = x 2 , p 9 = y 2 , p 10 = z 2 p_1 = 1, p_2 = x, p_3=y, p_4=z, \\p_5=xy, p_6=xz,p_7=yz\\p_8=x^2, p_9=y^2, p_{10}=z^2 p1=1,p2=x,p3=y,p4=z,p5=xy,p6=xz,p7=yzp8=x2,p9=y2,p10=z2 |
不同阶的基函数以获得不同的精度,取不同的权函数以改变拟合曲线(曲面)的光滑度。
MLS中,需要通过调整系数 α ( x ) \alpha(x) α(x),使得节点附近采样点取值与拟合函数在采样点取值之间的差的加权平方和最小。由此建立最优化模型,对向量 α ( x ) = [ α 1 ( x ) , α 2 ( x ) , Λ , a m ( x ) ] T \alpha(x)=[\alpha_1(x), \alpha_2(x), \Lambda, a_m(x)]^T α(x)=[α1(x),α2(x),Λ,am(x)]T计算加权L2范数 ∣ ∣ x ∣ ∣ 2 = ( ∑ i = 1 n x i 2 ) 1 2 ||x||_2 = (\sum_{i=1}^nx_i^2)^{\frac{1}{2}} ∣∣x∣∣2=(∑i=1nxi2)21, 如下所示:
J = ∑ I = 1 n w ( x − x I ) [ f ( x ) − y I ] 2 = ∑ I = 1 n w ( x − x I ) [ p T ( x ) α ( x ) − y I ] 2 \begin{equation}J = \sum_{I=1}^nw(x-x_I)[f(x)-y_I]^2=\sum_{I=1}^nw(x-x_I)[p^T(x)\alpha(x) -y_I]^2\end{equation} J=I=1∑nw(x−xI)[f(x)−yI]2=I=1∑nw(x−xI)[pT(x)α(x)−yI]2
其中, n n n是影响区域内节点的数目,如下图所示。
f ( x ) f(x) f(x)是拟合函数, y I y_I yI是 x = x I x=x_I x=xI处的节点值, y I = y ( x I ) y_I=y(x_I) yI=y(xI)。 w ( x − x I ) w(x-x_I) w(x−xI)是节点 x I x_I xI的权函数(后续讲解)。
为了确定系数 α ( x ) \alpha(x) α(x),式子 ( 2 ) (2) (2)应该取最小值,因此对 α ( x ) \alpha(x) α(x)求导可得,
∂ J ∂ α = A ( x ) α ( x ) − B ( x ) y = 0 \begin{equation} \frac{\partial J}{\partial \alpha} = A(x)\alpha(x)-B(x)y = 0 \end{equation} ∂α∂J=A(x)α(x)−B(x)y=0
α ( x ) = A − 1 ( x ) B ( x ) y \begin{equation} \alpha(x) = A^{-1}(x) B(x) y \end{equation} α(x)=A−1(x)B(x)y
其中
A
(
x
)
=
∑
i
=
1
n
w
(
x
−
x
I
)
p
(
x
I
)
p
T
(
x
I
)
A(x) = \sum_{i=1}^n w(x-x_I)p(x_I)p^T(x_I)
A(x)=i=1∑nw(x−xI)p(xI)pT(xI)
B ( x ) = [ w ( x − x I ) p ( x I ) , w ( x − x 2 ) p ( x 2 ) , Λ , w ( x − x n ) p ( x n ) ] B(x)=[w(x-x_I)p(x_I),w(x-x_2)p(x_2),\Lambda,w(x-x_n)p(x_n)] B(x)=[w(x−xI)p(xI),w(x−x2)p(x2),Λ,w(x−xn)p(xn)]
y T = [ y 1 , y 2 , Λ , y n ] y^T=[y_1, y_2, \Lambda, y_n] yT=[y1,y2,Λ,yn]
最后,将式子 ( 4 ) (4) (4)带入式子 ( 1 ) (1) (1)中,可以得到MSL拟合函数:
f
(
x
)
=
∑
i
=
1
m
α
i
(
x
)
p
i
(
x
)
=
p
T
(
x
)
α
(
x
)
=
∑
i
=
1
m
p
T
(
x
)
A
−
1
(
x
)
B
(
x
)
y
=
∑
I
=
1
n
ϕ
I
k
(
x
)
y
I
=
O
¨
k
(
x
)
y
f(x) = \sum_{i=1}^m \alpha_i(x)p_i(x)= p^T(x)\alpha(x) \\= \sum_{i=1}^m p^T(x) A^{-1}(x) B(x) y \\=\sum_{I=1}^n \phi_I^k(x)y_I=\ddot{O}^k(x)y
f(x)=i=1∑mαi(x)pi(x)=pT(x)α(x)=i=1∑mpT(x)A−1(x)B(x)y=I=1∑nϕIk(x)yI=O¨k(x)y
其中
=
O
¨
k
(
x
)
=\ddot{O}^k(x)
=O¨k(x)称为形函数,
k
k
k表示基函数的阶数。
O ¨ k ( x ) = [ ϕ 1 k , ϕ 2 k , Λ , ϕ n k ] = p T ( x ) A − 1 ( x ) B ( x ) \ddot{O}^k(x) = [\phi_1^k, \phi_2^k, \Lambda, \phi_n^k] = p^T(x) A^{-1}(x) B(x) O¨k(x)=[ϕ1k,ϕ2k,Λ,ϕnk]=pT(x)A−1(x)B(x)
需要注意的是,即使基函数 p ( x ) p(x) p(x)为多项式,MLS拟合函数的 f ( x ) f(x) f(x)也不再是多项式。
2. 权函数
上面公式 ( 2 ) (2) (2)中,有一个权函数 w ( x − x I ) w(x-x_I) w(x−xI)。MLS认为点 x x x处的值 y y y只受 x x x附近子域内节点影响,这个子域称作点 x x x的影响区域。 影响区域上定义一个权函数 w ( x ) w(x) w(x),如果权函数在这个区域内取为常数,就得到传统的最小二乘法。权函数在 x x x的一个子域内不等于零,在这个子域之外全为零,这个子域称为权函数的支持域(即 x x x的影响区域)。一般选择圆形作为权函数的支持域,其半径记为 S m a x S_{max} Smax,如上图1所示。由于权函数的紧支性,只有这些包含在影响区域内的数据点对点 x x x 的取值有影响。
- 权函数 w ( x − x I ) w(x-x_I) w(x−xI)应该是非负的,并且随着 ∣ ∣ x − x I ∣ ∣ 2 ||x-x_I||_2 ∣∣x−xI∣∣2的增加单调递减。
- 权函数具有一定的光滑性,因为拟合函数会继承权函数的连续性。
- 如果权函数 w ( x − x I ) w(x-x_I) w(x−xI)是 C 1 C^1 C1阶连续的,则拟合函数也是 C 1 C^1 C1阶连续的。
- 影响区域应该包含足够多的节点,以使得 A ( x ) A(x) A(x)可逆。
四、Performance
不同权函数的表现:
● 图 3( a)使用的权函数在整个区域内都是常数,这时移动最小二乘法相当于传统的最小二乘法,其行为相当于是一个线性拟合。
● 图 3( b)中权函数在一个小的子域内等于常数, 在其它地方等于零。 虽然这个权函数也具有紧支性, 但是不光滑, 它导致的拟合曲线就是将所有的数据点连线,即相当于分段线性插值。
● 图 3( c)中权函数是三次样条函数,不但具有紧支特性, 而且光滑。 即使使用的是线性基函数, 但由于继承了权函数的连续性, 拟合曲线仍然具有很好的光滑性。
图 3 说明了 MLS 中权函数的重要性。在最后一组图中, 仅使用线性基函数, 就获得很好的拟合结果,这是传统最小二乘法和其它拟合方法无法做到的。
曲面拟合情况:
五、Conclusion
移动最小二乘法使生成的曲线曲面具有精度高、 光滑性好等许多优点。
- 不需要事先确定拟合函数的类型;
- 由于紧支概念的引入,不需要进行分段(分块)拟合;
- 不需要求解线性方程组,从而避免了求解时方程组的系数矩阵可能是病态的情况;
- 精度比较高,能够捕捉到数据的剧烈变化;
- 只要选择合适的基函数和权函数,可以得到足够光滑的拟合曲线(曲面)。
移动最小二乘法进行曲线曲面拟合一个比较大的缺点就是计算过程比其它拟合方法要稍微复杂, 计算量也相对要大一些。