鄧興升 湯仲安

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

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

鄧興升1)湯仲安2)

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

式中

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

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

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

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

3 實(shí)例分析

3.1 插值精度比較

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

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

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

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

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

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

由圖3可知，Shepard反距離加權(quán)插值曲面與真實(shí)曲面相差較大，誤差最大值為3.183，最小值為-2.727。由圖4可知，本文方法插值曲面與真實(shí)曲面在區(qū)域中央誤差接近為0，沒有“上凸”或“下凹”現(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)，由不同方法分別計算100個測試點(diǎn)的插值精度。為了在同等條件下進(jìn)行比較，3種方法都采用了全部已知點(diǎn)進(jìn)行插值，表1為3種方法插值中誤差的比較結(jié)果。

圖2 二元三次多項式擬合誤差曲面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插值法和二元三次多項式十參數(shù)法精度沒有改善，甚至降低;而本文方法插值精度隨著已知點(diǎn)數(shù)量的增加不斷提高。

3.2 計算時間比較

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

表2 整體插值與移動插值的計算時間比較Tab.2 Comparison between the computation time as k= 3 564 and k=100

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

圖5 最鄰近100點(diǎn)移動插值與整體插值的格網(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á)閾值時即出現(xiàn)飽和現(xiàn)象，即插值精度不再提高，此時采用最鄰近k點(diǎn)插值與采用全部點(diǎn)插值的結(jié)果相同，因此沒有必要采用全部點(diǎn)插值，從而在保障精度的同時又節(jié)省計算時間。

4 結(jié)論

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

1 尤淑撐，嚴(yán)泰來.基于人工神經(jīng)網(wǎng)絡(luò)面插值的方法研究［J］.測繪學(xué)報，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 鄧興升，郭云開，花向紅.似大地水準(zhǔn)面格網(wǎng)雙二次多項式插值方法［J］.測繪學(xué)報，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é)報，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年生，高級工程師，博士，主要從事大地測量數(shù)據(jù)處理研究.E-mail：whudxs@163.com

P207

A