胡沛然，陳少輝

(1.中國(guó)科學(xué)院地理科學(xué)與資源研究所 陸地水循環(huán)及地表過(guò)程重點(diǎn)實(shí)驗(yàn)室，北京 100101；2.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引言

人們?cè)谕ㄟ^(guò)各種手段獲取觀測(cè)數(shù)據(jù)時(shí)，觀測(cè)儀器記錄的要素一般為特定格式的離散資料。由于觀測(cè)儀器故障、惡劣的環(huán)境或其他不可避免的外部因素干擾，觀測(cè)資料會(huì)出現(xiàn)一定數(shù)量的缺測(cè)值或異常值，使得用戶無(wú)法正常使用觀測(cè)數(shù)據(jù)，這就需要用戶選取合適的數(shù)學(xué)插值方法來(lái)填補(bǔ)數(shù)據(jù)空缺并處理數(shù)據(jù)的異常值。另一方面，雖然現(xiàn)代的觀測(cè)技術(shù)和觀測(cè)手段得到了很大進(jìn)步，用戶可以獲取各種高質(zhì)量的觀測(cè)資料，但是我們的研究尺度往往會(huì)和觀測(cè)資料的尺度無(wú)法完全匹配，例如人們獲得的衛(wèi)星遙感圖像或氣象資料與實(shí)際研究區(qū)域的空間尺度相差很大，為了和實(shí)際的研究精度相匹配，就必須要對(duì)觀測(cè)數(shù)據(jù)進(jìn)行插值處理，才能滿足數(shù)據(jù)使用者的特定需求。

拉格朗日插值法最早于1779年由一位英國(guó)的數(shù)學(xué)家提出[1]，后來(lái)由于著名科學(xué)家拉格朗日在他的一本書(shū)中又提到了這個(gè)經(jīng)典的插值方法而得以命名[2]。拉格朗日插值法具有簡(jiǎn)潔明了的迭代過(guò)程和較高的計(jì)算精度，因而十分適合通過(guò)計(jì)算機(jī)編程來(lái)處理各種數(shù)據(jù)。近幾十年來(lái)，拉格朗日插值方法在全球定位系統(tǒng)(global positioning system，GPS)衛(wèi)星定位[3]、水文氣象模型[4]和計(jì)算機(jī)圖像視頻處理[5]等領(lǐng)域均有成熟的應(yīng)用。根據(jù)拉格朗日插值方法的理論可以得出，計(jì)算結(jié)果的精度會(huì)隨著多項(xiàng)式階數(shù)的增加而大幅提高，但階數(shù)的增加可能會(huì)造成龍格現(xiàn)象(runge phenomenon)。具體來(lái)說(shuō)，簡(jiǎn)單的一階拉格朗日插值不會(huì)帶來(lái)數(shù)值振蕩，然而計(jì)算結(jié)果會(huì)出現(xiàn)一定的數(shù)值誤差，無(wú)法達(dá)到用戶所需要的數(shù)值精度。二階甚至更高階的拉格朗日插值可以有效減小插值結(jié)果的誤差，同時(shí)也會(huì)導(dǎo)致插值區(qū)間上出現(xiàn)顯著的數(shù)值振蕩(即為龍格現(xiàn)象)。雖然拉格朗日插值具有諸多優(yōu)點(diǎn)，但高次插值所引發(fā)的龍格現(xiàn)象卻在很大程度上限制了該方法在實(shí)際中的廣泛應(yīng)用[6]。為了解決這個(gè)問(wèn)題，目前的主要方法有適當(dāng)增加計(jì)算區(qū)間邊緣插值點(diǎn)個(gè)數(shù)和分段使用低次插值多項(xiàng)式等。在近年的研究中，Zhang在海洋模型的鹽度計(jì)算中采用混合一階和二階的拉格朗日插值法[7]。Wu等提出了一種混合N階拉格朗日插值方法[8]，即首先使用N階拉格朗日插值法進(jìn)行插值處理，然后在結(jié)果誤差較大的地方使用低階拉格朗日插值法進(jìn)行修正，這種方法對(duì)提高插值精度和減小數(shù)值振蕩有顯著效果。熊歡歡等在GPS星歷插值計(jì)算中采用高階和低階組合的滑動(dòng)插值方法，用以提高插值區(qū)間邊緣的計(jì)算精度和減弱龍格現(xiàn)象[9]。壽媛等使用Multiquadric(MQ)擬插值法和分段線性插值避免高次插值時(shí)出現(xiàn)的龍格現(xiàn)象[10]。

在陸面模型進(jìn)行徑流模擬時(shí)，研究者發(fā)現(xiàn)在對(duì)驅(qū)動(dòng)模型的氣象強(qiáng)迫數(shù)據(jù)進(jìn)行空間降尺度的過(guò)程中，原始的拉格朗日插值方法雖然易于編程且計(jì)算速度快，但插值效果卻呈現(xiàn)出較大波動(dòng)。增加區(qū)間邊緣插值點(diǎn)個(gè)數(shù)和使用混合高低階插值方法在一定程度上可以有效減弱高階拉格朗日插值的龍格現(xiàn)象，但這些方法會(huì)涉及復(fù)雜的判斷和較高的計(jì)算量，并不適合編程處理大量數(shù)據(jù)。針對(duì)這個(gè)問(wèn)題，本文提出了一種相對(duì)簡(jiǎn)便的改進(jìn)方法，在滿足一定精度的情況下，快速有效地減少了高階拉格朗日插值方法所帶來(lái)的數(shù)值振蕩。

1 方法介紹

1.1 拉格朗日插值原理

拉格朗日插值法的基本原理是通過(guò)迭代計(jì)算求得一個(gè)特定的多項(xiàng)式，其在每個(gè)觀測(cè)點(diǎn)恰好取到給定的觀測(cè)值。插值方法的定義如下：

對(duì)某個(gè)多項(xiàng)式，已知有k+1個(gè)給定點(diǎn)：(x0，y0)，…，(xk，yk)，假設(shè)任意兩個(gè)xj都互不相同，那么拉格朗日插值多項(xiàng)式為：

(1)

把式中每個(gè)lj(x)叫做基函數(shù)，lj(x)具體的表達(dá)式為[11]：

(2)

從上面的定義中可以看出，拉格朗日插值法原理簡(jiǎn)單，僅涉及到了簡(jiǎn)單的迭代計(jì)算，基函數(shù)的取值只與插值點(diǎn)坐標(biāo)有關(guān)，非常適合理論分析和編程計(jì)算。然而在實(shí)際的應(yīng)用過(guò)程中，插值多項(xiàng)式會(huì)隨著階數(shù)的增加而在插值區(qū)間邊緣產(chǎn)生強(qiáng)烈的數(shù)值不穩(wěn)定的現(xiàn)象，也就是說(shuō)插值多項(xiàng)式雖然在給定的觀測(cè)點(diǎn)取到了觀測(cè)值，但在插值區(qū)間邊緣附近卻會(huì)產(chǎn)生巨大的偏差，這被稱為龍格現(xiàn)象。如圖1所示，雖然多項(xiàng)式在插值區(qū)間的(-1，1)這部分精度很高，但是在區(qū)間的邊緣部分產(chǎn)生了大幅的數(shù)值振蕩。

圖1 插值中的龍格現(xiàn)象

1.2 拉格朗日插值法的一種改進(jìn)方法

拉格朗日插值方法在高階插值時(shí)會(huì)出現(xiàn)顯著的龍格現(xiàn)象，在模擬一個(gè)較為平穩(wěn)的函數(shù)時(shí)，插值多項(xiàng)式的數(shù)值不穩(wěn)定性更加強(qiáng)烈，究其原因是由基函數(shù)計(jì)算的插值權(quán)重所致，為此提出了一種權(quán)重歸一化方法。數(shù)據(jù)歸一化的主要思想是通過(guò)一系列數(shù)值計(jì)算，把原始數(shù)據(jù)映射到特定的范圍內(nèi)。數(shù)據(jù)的歸一化計(jì)算可以將數(shù)據(jù)按比例縮放，使其落入一個(gè)規(guī)定的區(qū)間內(nèi)(例如0～1或-1～1之間)。歸一化后既可以保留原始數(shù)據(jù)的某些特性，又可以消除數(shù)據(jù)的異常波動(dòng)，在某種程度上可以提高插值精度。歸一化有多種方法，其中一種常用的簡(jiǎn)單歸一化方法定義為：對(duì)于正項(xiàng)序列x1，x2，…，xn，進(jìn)行如下變換：

(3)

本文針對(duì)拉格朗日插值方法中的龍格現(xiàn)象，發(fā)展了一種基于歸一化思想的插值改進(jìn)方法，通過(guò)對(duì)插值基函數(shù)進(jìn)行歸一化計(jì)算將其轉(zhuǎn)換在(0，1)的范圍內(nèi)，減少了多項(xiàng)式在區(qū)間內(nèi)的數(shù)值振蕩。相比于目前主流的解決思路，新的改進(jìn)方法在插值計(jì)算時(shí)簡(jiǎn)便易行，精度在可接受范圍內(nèi)。具體步驟如下：

(4)

②將新的基函數(shù)ωj(x)帶入多項(xiàng)式(1)中，得到改進(jìn)后的插值多項(xiàng)式如下：

(5)

2 改進(jìn)前后插值結(jié)果對(duì)比實(shí)驗(yàn)

2.1 插值效果主觀對(duì)比

為了檢驗(yàn)插值方法改進(jìn)的效果，本文設(shè)計(jì)了一組主觀對(duì)比實(shí)驗(yàn)。我們從中國(guó)區(qū)域地面氣象要素?cái)?shù)據(jù)集中選取一片區(qū)域(經(jīng)緯度范圍為114.85°E～119.05°E，35.65°N～38.25°N)，并以2001年2月1日12時(shí)、2001年4月12日15時(shí)、2001年7月22日18時(shí)3個(gè)時(shí)刻的氣溫(2 m)、氣壓、向下短波輻射和向下長(zhǎng)波輻射數(shù)據(jù)作為此次實(shí)驗(yàn)的樣本數(shù)據(jù)。具體的插值對(duì)比方法為：

①考慮到計(jì)算精度和效率，采用五階拉格朗日插值法來(lái)進(jìn)行滑動(dòng)插值計(jì)算，即每次插值計(jì)算時(shí)，待插值點(diǎn)必須始終位于插值區(qū)間的中心。

②使用傳統(tǒng)拉格朗日插值方法將樣本數(shù)據(jù)(空間分辨率為0.1°)內(nèi)插成0.02°的空間分辨率，插值前后的空間范圍不變。

③使用本文改進(jìn)的拉格朗日插值法將樣本數(shù)據(jù)內(nèi)插成0.02°的空間分辨率，插值前后的空間范圍不變。

④通過(guò)Fortran語(yǔ)言編寫(xiě)程序來(lái)實(shí)現(xiàn)數(shù)據(jù)讀取和插值計(jì)算，并輸出圖像直觀對(duì)比插值方法改進(jìn)前后的效果。

2.2 插值結(jié)果精度分析

上文2.1中的對(duì)比實(shí)驗(yàn)只是通過(guò)插值結(jié)果與原始數(shù)據(jù)的定性圖形對(duì)比來(lái)主觀判斷插值效果。為了準(zhǔn)確地定量比較改進(jìn)前后插值方法精度的變化，我們?cè)谥袊?guó)區(qū)域地面氣象要素?cái)?shù)據(jù)集中取出某一時(shí)刻的地面向下長(zhǎng)波輻射數(shù)據(jù)作為定量對(duì)比實(shí)驗(yàn)的樣本數(shù)據(jù)，樣本數(shù)據(jù)尺寸為172列乘126列(經(jīng)緯度范圍為104.25°E～121.45°E，23.45°N～36.05°N)。將樣本數(shù)據(jù)(圖2)重采樣至分辨率更低的水平(圖3)，隨后再將其分別使用原始和改進(jìn)后的拉格朗日插值法，按照2.1中的五階滑動(dòng)插值重新內(nèi)插為0.1°的空間分辨率。將插值后得到的2組數(shù)據(jù)分別與樣本數(shù)據(jù)做誤差分析，畫(huà)出相對(duì)誤差頻率直方圖進(jìn)行比較，并選取3個(gè)常用的誤差指數(shù)作為本實(shí)驗(yàn)插值精度的評(píng)價(jià)指標(biāo)，具體見(jiàn)表1，表中的Xi與Zi分別為樣本數(shù)據(jù)和插值計(jì)算后的數(shù)據(jù)。最后對(duì)比原始樣本數(shù)據(jù)和使用改進(jìn)方法插值得到的結(jié)果，進(jìn)而分析改進(jìn)方法的誤差空間分布特征。

圖2 樣本數(shù)據(jù)

圖3 重采樣數(shù)據(jù)

表1 實(shí)驗(yàn)精度評(píng)價(jià)指標(biāo)

2.3 樣本數(shù)據(jù)介紹

中國(guó)區(qū)域地面氣象要素?cái)?shù)據(jù)集(China meteorological forcing dataset)，是一套以現(xiàn)有的多種再分析資料融合氣象觀測(cè)數(shù)據(jù)制作而成的氣象要素?cái)?shù)據(jù)集。此數(shù)據(jù)集由中國(guó)科學(xué)院青藏高原所研制(寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心http：//westdc.westgis.ac.cn/)，適合驅(qū)動(dòng)陸面過(guò)程模型進(jìn)行模擬計(jì)算。數(shù)據(jù)集在空間上覆蓋全國(guó)，時(shí)間范圍上覆蓋了1979—2010年，時(shí)空分辨率為3 h和0.1°，數(shù)據(jù)集包含了表2中的7個(gè)變量[12]。

表2 數(shù)據(jù)集變量信息表

3 插值結(jié)果分析

圖4、圖5、圖6分別為實(shí)驗(yàn)區(qū)域3個(gè)時(shí)刻的近地面氣溫、氣壓、向下短波輻射、向下長(zhǎng)波輻射樣本數(shù)據(jù)。樣本數(shù)據(jù)的空間分辨率為0.1°×0.1°，顏色越深表示像元的數(shù)值越低。從圖中可以清晰看出，原始數(shù)據(jù)集的空間分辨率較低，4個(gè)變量在空間上的分布極不連續(xù)，呈現(xiàn)出明顯的斑塊狀分布。在使用這種低分辨率氣象數(shù)據(jù)驅(qū)動(dòng)陸面模型進(jìn)行高分辨率模擬前，必須采用一定的數(shù)學(xué)方法對(duì)原始數(shù)據(jù)進(jìn)行降尺度插值計(jì)算。

圖4 2001年2月1日12時(shí)原始數(shù)據(jù)

圖5 2001年4月12日15時(shí)原始數(shù)據(jù)

圖6 2001年7月22日18時(shí)原始數(shù)據(jù)

圖7、圖8和圖9是使用傳統(tǒng)拉格朗日插值法將3個(gè)時(shí)刻的數(shù)據(jù)內(nèi)插成空間分辨率為0.02°后輸出得到的圖像。從圖中可以看出插值結(jié)果很不理想，在數(shù)據(jù)的邊緣處效果最差。其余部分雖然變化趨勢(shì)大致同樣本數(shù)據(jù)相同，但仍可以看出有明顯的數(shù)值波動(dòng)。出現(xiàn)這種現(xiàn)象是由于拉格朗日插值法在高階插值時(shí)會(huì)在插值區(qū)間邊緣產(chǎn)生龍格現(xiàn)象，尤其是在模擬類(lèi)似于氣象強(qiáng)迫數(shù)據(jù)這種總體平穩(wěn)的數(shù)據(jù)時(shí)龍格現(xiàn)象更加強(qiáng)烈。另一方面，在數(shù)據(jù)起始和結(jié)束的邊緣區(qū)域，數(shù)據(jù)無(wú)法滿足滑動(dòng)插值時(shí)待插值點(diǎn)始終位于插值區(qū)間中心的要求，故數(shù)值誤差在邊緣處最大。

圖10、圖11和圖12是使用本文改進(jìn)后的拉格朗日插值法將樣本數(shù)據(jù)內(nèi)插成空間分辨率為0.02°輸出得到的圖像。從輸出結(jié)果可以清晰看出，改進(jìn)后的插值結(jié)果不僅在變化趨勢(shì)上完全與樣本數(shù)據(jù)吻合，且基本可以有效避免原始插值方法中的缺陷，消除了明顯的數(shù)值誤差。

圖7 2001年2月1日12時(shí)

圖8 2001年4月12日15時(shí)

圖13和圖14為定量對(duì)比實(shí)驗(yàn)2.2輸出的傳統(tǒng)方法和改進(jìn)方法插值結(jié)果。從插值效果上看，原始插值方法產(chǎn)生了大量的數(shù)值波動(dòng)，未能準(zhǔn)確模擬出原始數(shù)據(jù)。而使用改進(jìn)方法插值得出的結(jié)果無(wú)論在整體趨勢(shì)上還是數(shù)值精度上，都可以較為真實(shí)地反映原始數(shù)據(jù)的空間分布特征。圖15為傳統(tǒng)方法和改進(jìn)方法的誤差頻率分布直方圖。從直方圖可以看出，傳統(tǒng)方法的相對(duì)誤差的頻率極值出現(xiàn)在0.1附近，相對(duì)誤差的范圍在-0.2到0.5之間，其概率分布無(wú)明顯規(guī)律。而方法改進(jìn)后的相對(duì)誤差大致呈現(xiàn)出正態(tài)分布，頻率極值出現(xiàn)在0附近，相對(duì)誤差的范圍分布在±0.2以內(nèi)，即改進(jìn)后的誤差范圍相比于傳統(tǒng)方法縮小了大約一半。從表3中的精度指標(biāo)數(shù)值可以看出，使用改進(jìn)方法得出的插值結(jié)果平均絕對(duì)誤差(MAE)從傳統(tǒng)方法的19.771下降為5.610，平均相對(duì)誤差(MRE)從8.339下降到2.284，而均方根誤差(RMSE)降幅最大，由34.926下降為8.584。通過(guò)定量的精度評(píng)估可以得出，使用改進(jìn)的拉格朗日插值方法后，數(shù)據(jù)的插值精度得到了大幅提高。圖16為使用改進(jìn)方法插值后數(shù)據(jù)與原始數(shù)據(jù)的相對(duì)誤差分布圖。通過(guò)對(duì)比分析，可以看出改進(jìn)后的拉格朗日插值方法在原始數(shù)據(jù)數(shù)值平穩(wěn)的地方插值精度較高，相對(duì)誤差大致控制在10%以內(nèi)，能夠較為準(zhǔn)確地反映出原始數(shù)據(jù)的空間分布特征。但是在數(shù)據(jù)邊緣部分和數(shù)據(jù)變化較大的地方插值精度較低，其相對(duì)誤差最高達(dá)到了30%，故改進(jìn)后的拉格朗日插值法不適合用于波動(dòng)劇烈數(shù)據(jù)的插值計(jì)算。

圖9 2001年7月22日18時(shí)

圖10 2001年2月1日12時(shí)

圖11 2001年4月12日15時(shí)

圖12 2001年7月22日18時(shí)

圖13 傳統(tǒng)方法插值結(jié)果圖

圖14 改進(jìn)方法插值結(jié)果

圖15 改進(jìn)前后相對(duì)誤差頻率直方圖

表3 拉格朗日插值法精度評(píng)價(jià)表

圖16 改進(jìn)方法的誤差空間分布

4 結(jié)束語(yǔ)

通過(guò)對(duì)上述兩個(gè)實(shí)驗(yàn)的定性比較與定量分析，可以得到以下3點(diǎn)結(jié)論：

①原始數(shù)據(jù)在實(shí)際應(yīng)用時(shí)為了和研究精度相匹配，必須綜合考慮計(jì)算精度和復(fù)雜度來(lái)選取合適的插值方法，在保證數(shù)值精度的前提下進(jìn)行插值處理。

②拉格朗日插值方法的插值原理簡(jiǎn)單，公式中僅包含了簡(jiǎn)單的迭代運(yùn)算，故十分適合編程來(lái)處理數(shù)據(jù)，當(dāng)選取合適的階數(shù)進(jìn)行插值時(shí)，插值的精度完全可以滿足用戶需求。

③拉格朗日插值方法在高階計(jì)算時(shí)會(huì)發(fā)生明顯的龍格現(xiàn)象，尤其是在插值類(lèi)似于氣象數(shù)據(jù)這種相對(duì)平穩(wěn)而又存在局部小幅波動(dòng)的數(shù)據(jù)時(shí)尤為顯著。如果采用滑動(dòng)拉格朗日插值法時(shí)，數(shù)據(jù)邊緣也會(huì)產(chǎn)生極大的誤差。本文提出的基于歸一化思想的改進(jìn)方法，通過(guò)重新計(jì)算拉格朗日插值基函數(shù)后，因龍格現(xiàn)象產(chǎn)生的誤差被有效減弱。改進(jìn)的拉格朗日插值方法在處理數(shù)值上相對(duì)平穩(wěn)的數(shù)據(jù)時(shí)效果比較理想，但不適用于處理數(shù)值波動(dòng)劇烈的數(shù)據(jù)。因此在陸面模型氣象驅(qū)動(dòng)數(shù)據(jù)空間降尺度時(shí)，可選用本文提出的朗格拉日插值改進(jìn)方法。