基于C1自然鄰近插值的曲面擬合

聶志峰， 周慎杰， 王 凱， 孔勝利

（1. 山東大學機械工程學院，山東 濟南 250061；2. 山東科技大學機械電子工程學院，山東 青島 265510）

基于C1自然鄰近插值的曲面擬合

聶志峰1,2， 周慎杰1， 王 凱1， 孔勝利1

（1. 山東大學機械工程學院，山東 濟南 250061；2. 山東科技大學機械電子工程學院，山東 青島 265510）

以C0連續non-Sibsonian 插值作為三次單純形Bernstein-Bézier多項式的基坐標，構造C1連續自然鄰近插值函數。介紹了高階連續函數的構建原理和性質。將C1連續自然鄰近插值函數應用于曲面擬合場合，由于Voronoi圖能夠自動調整數據點分布不規則和密度不均勻在空間上的差異，即使對于散亂數據點，也能獲得較好的擬合結果。

計算機應用；曲面擬合；C1自然鄰近插值；Bernstein-Bézier多項式

對于規則網格數據點的插值問題，通常采用樣條曲線來構建插值函數。Lai和Wenston采用連續三次樣條函數計算了橢圓型偏微分方程問題[1]。對于不規則數據點的曲線曲面擬合，通常采用最小二乘法[2]，通過使誤差的平方和最小，得到一個線性方程組，求解線性方程組可以得到擬合曲線或者曲面。但是，如果離散數據量比較大、形狀復雜，需要分段擬合和平滑處理，實際操作中存在一定的困難。為了克服以上困難，文獻[3]利用插值與逼近相結合的曲面擬合思路，較好地解決原始數據點分布不均勻的曲面造型問題。移動最小二乘法對最小二乘法作了改進，能夠生成精度高﹑光滑性好的曲線曲面[4-5]。自然鄰近插值法是一種較新的插值理論，主要用于多變量不規則數據點的插值。由于Voronoi圖和對偶的 Delaunay三角化能夠自動調整數據點分布不規則和密度不均勻在空間上的差異，即使對于散亂數據點，也可以得到較好的曲面擬合結果。另外，自然鄰近插值法具有嚴格的插值特性，偏微分方程Garlerkin法可以直接施加本質邊界條件，這是移動最小二乘法所無法相比的。Tily和Brace以汽車發動機半規則實驗數據為研究對象，對C0連續Sibson插值函數和MATLAB interp2 插值函數進行了比較，發現當實驗數據標準化處理后，Sibson 插值函數具有更高的插值精度[6]。高階偏微分方程的數值計算要求插值函數至少 C1連續，對 C1連續插值函數的插值性能進行深入研究是有必要的。Farin將C0連續Sibson插值函數以自然鄰近坐標的形式代入三次單純形Bernstein-Bézier曲面多項式，構造出C1連續插值函數[7]。Sukumar和Moran提出了一種新的計算方法，將Bernstein基函數轉換為C1形函數，保證任意結點I的3個形函數分別與它的3個自由度、和直接相聯，從而適合于高階偏微分方程 Galerkin法[8]。與Sibson插值比較，non-Sibsonian插值具有計算簡單、計算效率高等優點[9]，所以本文將non-Sibsonian 插值函數作為三次單純形Bernstein-Bézier曲面多項式的基坐標，采用文獻[8]的計算方法，構造出新的連續自然鄰近插值函數，該函數具有二次完備性、對結點函數值和梯度值的插值特性以及在分析域的邊界退化為三次多項式等性質。為了考察新建連續插值函數的插值精度及其影響因素，作者將它應用到曲面擬合場合。

1 C1自然鄰近插值

1.1 non-Sibsonian自然鄰近插值

Voronoi 圖及其對偶的Delaunay三角化源自計算幾何，是由一組不規則的點定義的最基本的幾何結構。考慮平面上由m個離散結點組成的集合，結點I的Voronoi圖是指與結點I對應的區域，內任意點到結點I的距離都小于該點到其它結點J的距離，用數學公式表示為

當積分點x被引進集合N的一階Voronoi 圖時，根據空圓準則可以確定點x的自然鄰子。圖1是平面上7個結點的一階Voronoi 圖，點x共有4個自然鄰子，分別是結點1、2、3和4。

假設點x有n個自然鄰子，該點關于結點I的non-Sibsonian插值函數可以定義為

1.2 以non-Sibsonian插值作為自然鄰近坐標的C1自然鄰近插值函數

將 non-Sibsonian插值φ代入三次單純形Bernstein-Bézier曲面多項式，沿用文獻[8]的計算方法，構造出C1自然鄰近插值函數

對公式（3）作如下變換，得到

1.2.1 Bernstein-Bézier基函數

non-Sibsonian坐標φ下由積分點x的n個自然鄰子組成的n邊形的三次基函數為

對于三次Bernstein-Bézier曲面，基函數的多指標記法只會出現下面3種情況：和。表示數組中只有第α個元素為1，其余為0的多指標記法。

1.2.2 變換矩陣的構建

通常，可以把Bézier空間坐標分為頂點、切面和中心Bézier空間坐標三種。三種Bézier空間坐標與結點函數值和梯度值具有不同的構造關系：頂點Bézier空間坐標等于結點函數值，即

切面Bézier空間坐標為

中心Bézier空間坐標的最佳方案為[7]

式中

由公式（4）可以看出：變換矩陣的作用是把結點函數值或者梯度值映射到Bézier空間坐標上。利用公式（6）～（8），可以構建出變換矩陣。一旦變換矩陣構建完畢，Bernstein-Bézier基函數和變換矩陣相乘，可以得到形函數

圖3 結點網格及?形函數ψ3A-2(x)、ψ3A-1(x)

（1） 高階連續性

（2） 二次完備性

式中 a為切面Bézier空間坐標的質心，

c為頂點Bézier空間坐標的質心。

（3） 對結點函數值和結點梯度值的插值特性

（4） 分析域邊界積分點的插值函數是三次多項式

如果點x的兩個自然鄰子的 non-Sibsonian坐標為，把它們代入方程（12），可以得到

顯然，分析域邊界上的插值函數是三次多項式。

2 曲面擬合

3個曲面函數分別定義為

圖4 單位正方形離散方案

為了衡量曲面擬合質量，定義整體域?的相對誤差范Le

圖5 規則25結點的3個函數曲面與擬合曲面的對比（左為函數曲面，右為擬合曲面）

Le值越小，曲面擬合結果越好。3個函數在5種離散方案下的相對誤差范Le，見表1。

表1 曲面擬合相對誤差范Le

3 結 論

（2） 對于高次多項式函數（大于等于3次），曲面擬合質量下降了。離散結點數量對曲面擬合質量影響較大，加密結點可以提高連續插值函數的插值精度，得到質量較高的擬合曲面。由于Voronoi圖能夠自動調整數據點分布不規則和密度不均勻在空間上的差異，所以對于不規則和密度不均勻數據點的情況，也得到了較好的曲面擬合結果。

[1] Lai M J, Wenston P. Report on numerical experi-ments with bivariate C1cubic splines for numerical solution of partial differential equations [R]. Technical Report, Department of Mathematics, University of Georgia, 1998.

[2] 彭芳瑜, 周 濟, 周艷紅, 等. 基于最小二乘法的曲面生成算法研究[J]. 工程圖學學報, 1999, 20(3): 41-46.

[3] 彭芳瑜, 周云飛, 周 濟. 基于插值與逼近的復雜曲面擬合[J]. 工程圖學學報, 2002, 23(4): 87-96.

[4] 曾清紅, 盧德唐. 基于移動最小二乘法的曲線曲面擬合[J]. 工程圖學學報, 2004, 25(1): 84-89.

[5] Lancaster P, Salkauskas K K. Surfaces generated by moving least squares methods [J]. Mathematics of Computation, 1981, (37): 141-158.

[6] Tily R, Brace C J. A study of natural neighbour interpolation and its application to automotive engine test data [J]. Proceeding of I mech E, Part D: Journal of Automobile Engineering, 2006, (220): 1003-1017.

[7] Farin G. Surfaces over Dirichlet tessellations [J]. Computer Aided Geometric Design, 1990, (7): 281-292.

[8] Sukumar N, Moran B.natural neighbor interpo-lant for partial differential equations [J]. Numerical Methods for Partial Differential Equations, 1999, 15(4): 417-447.

[9] Sukumar N, Moran B, Semenov A Y, et al. Natural neighbor Galerkin methods [J]. International Journal for Numerical Methods in Engineering, 2001, (50): 1-27.

Surface Fitting Based on C1Natural Neighbor Interpolation

NIE Zhi-feng1,2, ZHOU Shen-jie1, WANG Kai1, KONG Sheng-li1
( 1. School of Mechanical Engineering, Shandong University, Jinan Shandong 250061, China; 2. College of Mechanical and Electronic Engineering, Shandong University of Science and Technology, Qingdao Shandong 265510, China)

C1natural neighbor interpolation can be realized when C0non-Sibsonian interpolation is introduced in the Bernstein-Bézier surface representation of a cubic simplex. The principle and properties of C1natural neighbor interpolation are described and the surface fitting is given. The Voronoi diagram can adjust the spatial discrepancy caused by irregular data and data of varying density, so even for the scattered data, C1natural neighbor interpolation can get accurate fitting results.

computer application; surface fitting; C1natural neighbor interpolation; Bernstein-Bézier polynomial

TP 391

A

：1003-0158(2010)01-0110-06

2008-08-01

國家自然科學基金資助項目（10572077）；山東省自然科學基金資助項目（Y2007F20）；山東科技大學“春蕾計劃”資助項目（2009AZZ021）

聶志峰（1970-），男，山東日照人，博士研究生，主要研究方向為數值計算方法。