權重歸一化拉格朗日插值及其空間降尺度應用

胡沛然，陳少輝

(1.中國科學院地理科學與資源研究所 陸地水循環及地表過程重點實驗室，北京 100101；2.中國科學院大學，北京 100049)

0 引言

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

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

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

1 方法介紹

1.1 拉格朗日插值原理

拉格朗日插值法的基本原理是通過迭代計算求得一個特定的多項式，其在每個觀測點恰好取到給定的觀測值。插值方法的定義如下：

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

(1)

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

(2)

從上面的定義中可以看出，拉格朗日插值法原理簡單，僅涉及到了簡單的迭代計算，基函數的取值只與插值點坐標有關，非常適合理論分析和編程計算。然而在實際的應用過程中，插值多項式會隨著階數的增加而在插值區間邊緣產生強烈的數值不穩定的現象，也就是說插值多項式雖然在給定的觀測點取到了觀測值，但在插值區間邊緣附近卻會產生巨大的偏差，這被稱為龍格現象。如圖1所示，雖然多項式在插值區間的(-1，1)這部分精度很高，但是在區間的邊緣部分產生了大幅的數值振蕩。

圖1 插值中的龍格現象

1.2 拉格朗日插值法的一種改進方法

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

(3)

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

(4)

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

(5)

2 改進前后插值結果對比實驗

2.1 插值效果主觀對比

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

①考慮到計算精度和效率，采用五階拉格朗日插值法來進行滑動插值計算，即每次插值計算時，待插值點必須始終位于插值區間的中心。

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

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

④通過Fortran語言編寫程序來實現數據讀取和插值計算，并輸出圖像直觀對比插值方法改進前后的效果。

2.2 插值結果精度分析

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

圖2 樣本數據

圖3 重采樣數據

表1 實驗精度評價指標

2.3 樣本數據介紹

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

表2 數據集變量信息表

3 插值結果分析

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

圖4 2001年2月1日12時原始數據

圖5 2001年4月12日15時原始數據

圖6 2001年7月22日18時原始數據

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

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

圖7 2001年2月1日12時

圖8 2001年4月12日15時

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

圖9 2001年7月22日18時

圖10 2001年2月1日12時

圖11 2001年4月12日15時

圖12 2001年7月22日18時

圖13 傳統方法插值結果圖

圖14 改進方法插值結果

圖15 改進前后相對誤差頻率直方圖

表3 拉格朗日插值法精度評價表

圖16 改進方法的誤差空間分布

4 結束語

通過對上述兩個實驗的定性比較與定量分析，可以得到以下3點結論：

①原始數據在實際應用時為了和研究精度相匹配，必須綜合考慮計算精度和復雜度來選取合適的插值方法，在保證數值精度的前提下進行插值處理。

②拉格朗日插值方法的插值原理簡單，公式中僅包含了簡單的迭代運算，故十分適合編程來處理數據，當選取合適的階數進行插值時，插值的精度完全可以滿足用戶需求。

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