鄧興升 湯仲安

(1)長(zhǎng)沙理工大學(xué)測(cè)繪工程系，長(zhǎng)沙 410004 2)湖南省測(cè)繪科技研究所，長(zhǎng)沙410004)

移動(dòng)格林基函數(shù)樣條二維插值算法研究*

鄧興升1)湯仲安2)

(1)長(zhǎng)沙理工大學(xué)測(cè)繪工程系，長(zhǎng)沙 410004 2)湖南省測(cè)繪科技研究所，長(zhǎng)沙410004)

針對(duì)用于插值的已知點(diǎn)較多時(shí)，插值計(jì)算需要解算大規(guī)模矩陣、計(jì)算耗時(shí)長(zhǎng)甚至無(wú)法解算的問(wèn)題，引入移動(dòng)曲面的思想，取插值點(diǎn)周邊最鄰近k個(gè)已知點(diǎn)進(jìn)行格林基函數(shù)二維樣條移動(dòng)插值，實(shí)例計(jì)算結(jié)果表示，該方法的插值精度高于Shepard插值法與多項(xiàng)式擬合法的精度。插值范圍大及測(cè)點(diǎn)數(shù)量眾多時(shí)，該方法仍可用，無(wú)需數(shù)據(jù)分區(qū)與光滑接邊，與整體插值相比可大大降低計(jì)算時(shí)間。

移動(dòng)格林基函數(shù);二維樣條;插值算法;整體插值;移動(dòng)插值

1 引言

在大地測(cè)量領(lǐng)域，需要由離散觀測(cè)的型值點(diǎn)來(lái)構(gòu)造自由曲面。由于區(qū)域廣、范圍大，與空間位置相關(guān)的屬性難以用一個(gè)數(shù)學(xué)函數(shù)來(lái)表達(dá)，因而需先建立規(guī)則的格網(wǎng)，然后在格網(wǎng)的基礎(chǔ)上進(jìn)行插值計(jì)算。目前國(guó)內(nèi)外常采用的插值方法有：多項(xiàng)式擬合、Shepard反距離加權(quán)、Hardy多面函數(shù)、Kriging、人工神經(jīng)網(wǎng)絡(luò)［1］、最小曲率樣條插值［2］等。

多項(xiàng)式擬合存在以下不足：1)復(fù)雜曲面難以由有限參數(shù)的多項(xiàng)式精確擬合;2)構(gòu)造的曲面容易產(chǎn)生多余擺動(dòng)且難以發(fā)現(xiàn);3)多項(xiàng)式插值存在“龍格現(xiàn)象”，曲面邊緣會(huì)產(chǎn)生畸變;4)大面積曲面必須分區(qū)擬合，接邊時(shí)容易產(chǎn)生錯(cuò)位;5)多項(xiàng)式次數(shù)及參數(shù)個(gè)數(shù)需要依賴經(jīng)驗(yàn)調(diào)整或多次試算;6)當(dāng)曲面復(fù)雜時(shí)，必須用高次多項(xiàng)式，而高次多項(xiàng)式振動(dòng)大、不穩(wěn)定，用它進(jìn)行預(yù)測(cè)效果很差。Shepard插值原理是反距離加權(quán)插值，即將插值函數(shù)定義為各型值點(diǎn)的加權(quán)平均值，權(quán)函數(shù)定義為與距離成反比或采用分段函數(shù)，反距離加權(quán)插值會(huì)引起曲面畸形及錯(cuò)開(kāi)的現(xiàn)象［3］。Hardy多面函數(shù)的光滑因子需人工主觀確定，且對(duì)擬合結(jié)果影響顯著，這種方法有嚴(yán)重缺陷甚至錯(cuò)誤，已被否定［4］。

最小曲率樣條插值認(rèn)為插值點(diǎn)形成的曲面應(yīng)滿足整體曲率最小的原則，由于整體曲率最小在理論上很合理，生成的曲面具有光滑性，因而成為主流插值方法。然而該方法生成的曲面會(huì)在已知控制點(diǎn)之間產(chǎn)生多余的擺動(dòng)，Schweikert［5］提出張力樣條，以克服多余擺動(dòng)，Cline［6］實(shí)現(xiàn)了張力樣條。David［7］在最小曲率準(zhǔn)則的基礎(chǔ)上提出了格林基函數(shù)插值法，在構(gòu)造曲面方式上進(jìn)步較大，但當(dāng)型值點(diǎn)的分布極度不均衡時(shí)，插值結(jié)果欠穩(wěn)定，且同樣存在曲面多余擺動(dòng)的問(wèn)題［8，9］，文獻(xiàn)［8］采用張力樣條以克服該方法的多余擺動(dòng)，張力樣條得到廣泛認(rèn)同和應(yīng)用［9-11］。在已知點(diǎn)相對(duì)較密時(shí)，張力樣條會(huì)提高插值精度，但已知點(diǎn)分布不均勻及相對(duì)較稀時(shí)插值結(jié)果仍欠穩(wěn)定，張力樣條構(gòu)造曲線在某些情形下呈波浪形，產(chǎn)生抖動(dòng)現(xiàn)象，致使曲線失真［12］。

在大范圍插值問(wèn)題中，會(huì)碰到大量非規(guī)則分布的離散點(diǎn)，如果要解算大規(guī)模矩陣，計(jì)算量巨大，計(jì)算速度慢，隨著型值點(diǎn)的增加計(jì)算非常緩慢甚至無(wú)法進(jìn)行。因此本文引入移動(dòng)曲面的思想，以待插值點(diǎn)周邊最鄰近的若干點(diǎn)實(shí)施格林基函數(shù)二維樣條插值，從而減少計(jì)算時(shí)間;采用插值點(diǎn)周圍最鄰近的若干點(diǎn)進(jìn)行移動(dòng)插值，改變了型值點(diǎn)分布的均衡性，有利于提高插值結(jié)果的穩(wěn)定性。

2 移動(dòng)格林基函數(shù)樣條二維插值算法

曲面s(x)可表達(dá)為：

式中，x為已知點(diǎn)的位置參數(shù);g(x，xj)為格林函數(shù); xj為第個(gè)插值點(diǎn)的位置;wj是與xj相關(guān)聯(lián)的未知權(quán)值;T(x)為傾向參數(shù)，一般為常數(shù)，它代表不能由格林函數(shù)表達(dá)的那部分地區(qū)性傾向參數(shù)，在對(duì)曲面s(x)建模之前，首先要把T(x)從數(shù)據(jù)中移去，建模之后再恢復(fù)T(x)，從而不考慮對(duì)T(x)的建模。式(1)中權(quán)值wj是滿足

的解。其中格林函數(shù)g(xi，xj)表達(dá)式為

設(shè)某測(cè)區(qū)共有n個(gè)已知點(diǎn)，當(dāng)n=5 000時(shí)，在PC機(jī)(CPU：2G Hz，內(nèi)存1G)上進(jìn)行一次5 000階的矩陣求逆運(yùn)算需要4～5分鐘，如果格網(wǎng)達(dá)10 000個(gè)結(jié)點(diǎn)以上，則計(jì)算耗時(shí)漫長(zhǎng)。移動(dòng)插值認(rèn)為，插值只是局部行為，待插值點(diǎn)p0的值只與其周邊最鄰近的若干個(gè)點(diǎn)相關(guān)。本文以待插值點(diǎn)為圓心，變動(dòng)半徑使落入圓中的已知點(diǎn)數(shù)等于閾值k。移動(dòng)格林基函數(shù)樣條二維插值算法如下：

1)設(shè)測(cè)區(qū)共有n個(gè)已知點(diǎn)，插值點(diǎn)為 p0(x0，y0)，參與p0點(diǎn)插值計(jì)算的最鄰近點(diǎn)數(shù)為k，較小的k值(如k＜30)會(huì)影響計(jì)算精度;增加k值能提高插值精度，但達(dá)到某一閾值時(shí)(通常小于100)會(huì)出現(xiàn)飽和現(xiàn)象，即插值精度不再提高。根據(jù)p0(x0，y0)與已知點(diǎn)(xi，yi)計(jì)算其距離r0：

將r0i按歸并排序算法從小到大排序，取得距離最小的前k(k≤n)個(gè)已知點(diǎn)的序號(hào)及坐標(biāo)(xj，yj，zj)，j=1，2，…，k。

2)設(shè)前k個(gè)已知點(diǎn)構(gòu)成的矩陣為X=［x1x2…xk］T，Y=［y1y2… yk］T，Z=［z1z2… zk］T。設(shè)k×k階方陣為

式中

3)計(jì)算插值點(diǎn)p0(x0，y0)的1×k階矩陣Gp：

式中d01，d02，…，d0k按式(6)計(jì)算，式(6)中的r0j即為p0至k個(gè)已知點(diǎn)的距離。采用第1)步計(jì)算的相應(yīng)結(jié)果，則所求p0點(diǎn)插值z(mì)p為：

4)重復(fù)1)～3)步，依次求出其他點(diǎn)pi(xi，yi)的插值z(mì)pi。

從以上步驟可知，移動(dòng)格林基函數(shù)樣條二維插值算法有如下特點(diǎn)：算法簡(jiǎn)潔;算法基于點(diǎn)與點(diǎn)之間的距離及其衍生函數(shù)，而與方向無(wú)關(guān);已知點(diǎn)數(shù)量任意多時(shí)，插值只在最鄰近的k個(gè)點(diǎn)中進(jìn)行，因而仍可行;無(wú)需對(duì)數(shù)據(jù)分區(qū)與光滑接邊。

3 實(shí)例分析

3.1 插值精度比較

以二元三次多項(xiàng)式十參數(shù)法、Shepard插值法作為對(duì)比方法。已知曲面為：

其中a1=0.4，b1=0.3，a2=-0.8，b2=0.2。實(shí)驗(yàn)中變更這些參數(shù)，結(jié)果類似。由式(10)產(chǎn)生若干已知點(diǎn)和測(cè)試點(diǎn)，均無(wú)誤差，因而插值結(jié)果只與不同方法自身有關(guān)。在(x，y)的定義域內(nèi)以2.9為間隔取64點(diǎn)作為已知點(diǎn)，以間隔2.13取100點(diǎn)作為測(cè)試點(diǎn)。已知點(diǎn)構(gòu)造的曲面如圖1所示。

圖1 試驗(yàn)的已知曲面Fig.1 Experimental surface(8×8 grid)

二元三次多項(xiàng)式十參數(shù)法模型如下：

其中a～j為模型參數(shù);x、y為點(diǎn)的坐標(biāo);z為與x、y對(duì)應(yīng)的屬性值。根據(jù)最小二乘法計(jì)算出模型參數(shù)，再計(jì)算測(cè)試點(diǎn)，并與真值比較，其差異如圖2所示。

由圖2可知，二元三次多項(xiàng)式擬合會(huì)出現(xiàn)規(guī)則的變形，曲面中間出現(xiàn)規(guī)則的“鞍狀上凸”，邊緣出現(xiàn)了“翹曲”，誤差最大值為 0.354，最小值為-0.427。Shepard插值采用了以距離倒數(shù)為權(quán)值的加權(quán)插值，測(cè)試點(diǎn)插值結(jié)果與真值比較，其差異如圖3所示。采用本文方法由64個(gè)已知點(diǎn)插值得到的曲面與真實(shí)曲面比較，其差異曲面如圖4所示。

由圖3可知，Shepard反距離加權(quán)插值曲面與真實(shí)曲面相差較大，誤差最大值為3.183，最小值為-2.727。由圖4可知，本文方法插值曲面與真實(shí)曲面在區(qū)域中央誤差接近為0，沒(méi)有“上凸”或“下凹”現(xiàn)象;但在曲面內(nèi)側(cè)邊緣有畸變，誤差最大值為0.326，最小值為-0.201。

為了考察已知點(diǎn)數(shù)量增多的情況下各種方法的插值精度，將取點(diǎn)間隔適當(dāng)縮小，在(x，y)的定義域內(nèi)以間隔1.8取得144點(diǎn);以間隔1.5取得225點(diǎn);以間隔1.3取得289點(diǎn)，由不同方法分別計(jì)算100個(gè)測(cè)試點(diǎn)的插值精度。為了在同等條件下進(jìn)行比較，3種方法都采用了全部已知點(diǎn)進(jìn)行插值，表1為3種方法插值中誤差的比較結(jié)果。

圖2 二元三次多項(xiàng)式擬合誤差曲面Fig.2 Error surface of binary cubic polynomial fitting

圖3 Shepard插值誤差曲面Fig.3 Error surface of Shepard’s interpolation

圖4 格林基函數(shù)樣條插值誤差曲面Fig.4 Error surface of spline interpolation based on Green’s function

表1 不同已知點(diǎn)數(shù)量的情況下3種方法的中誤差比較Tab.1 Comparison among the mean square errors with three methods with different number of data points

從表1可知，Shepard方法在本例中的插值中誤差較大;格林基函數(shù)樣條二維插值法在本例中的插值精度最高;隨著已知點(diǎn)數(shù)量的增加，Shepard插值法和二元三次多項(xiàng)式十參數(shù)法精度沒(méi)有改善，甚至降低;而本文方法插值精度隨著已知點(diǎn)數(shù)量的增加不斷提高。

3.2 計(jì)算時(shí)間比較

已對(duì)某曲面均勻觀測(cè)了3 564個(gè)測(cè)點(diǎn)，欲根據(jù)觀測(cè)數(shù)據(jù)生成60×60的正方形格網(wǎng)，插值方法采用格林基函數(shù)二維樣條插值，本例對(duì)采用全部3 564點(diǎn)整體插值與采用最鄰近100點(diǎn)移動(dòng)插值的計(jì)算時(shí)間進(jìn)行比較，結(jié)果見(jiàn)表2。算法采用C++語(yǔ)言編程實(shí)現(xiàn)，測(cè)試計(jì)算機(jī)的配置為CPU：2G Hz，內(nèi)存1G。

表2 整體插值與移動(dòng)插值的計(jì)算時(shí)間比較Tab.2 Comparison between the computation time as k= 3 564 and k=100

由表2可知，采用最鄰近100點(diǎn)移動(dòng)插值可以大大縮短計(jì)算時(shí)間。為了考查采用最鄰近100點(diǎn)插值時(shí)的精度有沒(méi)有受到影響，因該曲面模型未知，插值結(jié)果無(wú)法與真值比較，我們將其與整體插值的格網(wǎng)進(jìn)行比較，兩格網(wǎng)之間的差異如圖5所示，差異最大值為2.7 mm，最小值為-3.3 mm，中誤差僅為0.4 mm，表明插值結(jié)果非常接近。

圖5 最鄰近100點(diǎn)移動(dòng)插值與整體插值的格網(wǎng)差異Fig.5 Interpolated Grids’differences between the results from moving interpolation and global interpolation

試驗(yàn)進(jìn)一步增加k值，當(dāng)k到達(dá)閾值時(shí)即出現(xiàn)飽和現(xiàn)象，即插值精度不再提高，此時(shí)采用最鄰近k點(diǎn)插值與采用全部點(diǎn)插值的結(jié)果相同，因此沒(méi)有必要采用全部點(diǎn)插值，從而在保障精度的同時(shí)又節(jié)省計(jì)算時(shí)間。

4 結(jié)論

針對(duì)插值過(guò)程中求解大規(guī)模矩陣?yán)щy的問(wèn)題，引入移動(dòng)曲面的思想，以待插值點(diǎn)周圍最鄰近的若干點(diǎn)進(jìn)行移動(dòng)格林基函數(shù)二維樣條插值，算法簡(jiǎn)潔，易于編程使用。實(shí)例中該方法插值精度高于二元三次多項(xiàng)式擬合和Shepard反距離加權(quán)插值的精度，且隨著已知點(diǎn)數(shù)量的增加能不斷地提高插值精度。已知點(diǎn)數(shù)量任意多時(shí)，整體插值難以進(jìn)行，移動(dòng)插值只在最鄰近的k點(diǎn)中進(jìn)行因而仍可實(shí)施，插值精度與整體插值相同，并且無(wú)需對(duì)數(shù)據(jù)分區(qū)與光滑接邊，減小過(guò)程復(fù)雜性，極大地縮短計(jì)算時(shí)間。

1 尤淑撐，嚴(yán)泰來(lái).基于人工神經(jīng)網(wǎng)絡(luò)面插值的方法研究［J］.測(cè)繪學(xué)報(bào)，2000，29(1)：30-34.(You Shucheng and Yan Tailai.A study on artificial neural network based on surface interpolation［J］.Acta Geodaetica et Cartograpphica Sinica，2000，29(1)：30-34)

2 Briggs I C.Machine contouring using minimum curvature［J］.Geophysics，1974，39(1)：39-48.

3 鄧興升，郭云開(kāi)，花向紅.似大地水準(zhǔn)面格網(wǎng)雙二次多項(xiàng)式插值方法［J］.測(cè)繪學(xué)報(bào)，2009，38(1)：35-40.(Deng Xingsheng，Guo Yunkai and Hua Xianghong.Quasigeoid grid network biquadratic interpolation method［J］.Acta Geodaetica et Cartograpphica Sinica，2009，38(1)：35-40)

4 Jekeli C.Hardy’s multiquadric-biharmonic method for gravity field predictions［J］.Computers＆Mathematical Applications，1994，28(7)：43-46.

5 Schweikert D G.An interpolating curve using a spline intension［J］.Jour Math Physics，1966，45(2)：312-317.

6 Cline A.Scalar and planar-value curve fitting using splines under tension［J］.Comm ACM.，1974，(17)：218-223.

7 Sandwell D T.Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data［J］.Geophys Res Lett.，1987，14(2)：139-142.

8 Wessel P and Bercovici D.Interpolation with splines in Tension：A Green’s function approach［J］.Mathematical Geology，1998，30(1)：77-93.

9 Smith W H F and Wessel P.Gridding with continuous curvature splines in tension［J］.Geophysics，1990，55(3)：293 -305.

10 Mitásová H and Hofierka J.Interpolation by regularized spline with tension：II.Application to terrain modeling and surface geometry analysis［J］.Mathematical Geology，1993，25(6)：657-669.

11 Mitásová H and Mitás L.Interpolation by regularized spline with tension：I.Theory and implementation［J］.Mathematical Geology，1993，25(6)：641-655.

12 鄧曙光，李婉.曲線光滑的張力樣條插值法VC實(shí)現(xiàn)［J］.工程地球物理學(xué)報(bào)，2005，2(5)：387-390.(Deng Shuguang and Li Wan.Realization of smoothing curve with tension spline interpolation under visual C++［J］.Chinese Journal of Engineering Geophysics，2005，2(5)：387 -390)

STUDY ON TWO-DIMENSIONAL SPLINE INTERPOLATION BASED ON MOVING GREEN FUNCTION

Deng Xingsheng1)and Tang Zhongan2)

(1)Department of Surveying Engineering，Changsha University of Science＆Technology，Changsha 410004 2)Hunan Research Institute of Surveying and Mapping，Changsha 410004)

When the data coverage is dense，some algorithms need to solve large size matrix，thus the computation time is proportional approximately to the cube of the number of data constraints，it makes the process very slow.Focusing on this problem，the moving curvature is introduced in interpolation.Only the nearest data points are chosen for interpolating by two-dimensional spline based on the moving Green’s function.The examples show that the interpolation accuracy of the proposed method is higher than that of two other methods.No matter how many data points there are，this method can be implemented fast.It is not necessary to split the data into subsets which can be modeled individually，or to blend the subsets together into a final model.Comparing with the global solution，this algorithm can greatly reduce the computation time.

moving Green’s function;two-dimensional spline;interpolation algorithm;global interpolation;moving interpolation

1671-5942(2011)06-0069-04

2011-05-06

湖南省自然科學(xué)基金(10JJ3090)

鄧興升，男，1971年生，高級(jí)工程師，博士，主要從事大地測(cè)量數(shù)據(jù)處理研究.E-mail：whudxs@163.com

P207

A