不同氣象插值方法精度評估及對草地NPP估算的影響

朱玉果, 杜靈通, 謝應忠,3, 劉 可, 宮 菲, 丹 楊, 王 樂

(1.寧夏大學 西北土地退化與生態系統恢復省部共建國家重點實驗室培育基地, 銀川 750021;2.寧夏大學 西北退化生態系統恢復與重建教育部重點實驗室, 銀川 750021; 3.寧夏大學 農學院, 銀川 750021)

植被凈初級生產力(Net Primary Productivity,NPP)是單位面積、單位時間內綠色植物通過光合作用所獲取的有機物質總量中扣除自養呼吸后的干物質總量[1]。其反映了在自然環境條件下植被對CO2的固定能力，不僅是植被活動的重要表征，而且是調節生態過程的重要依據，在生態系統的質量狀況、生產能力評估等方面發揮著重要作用[2]，自19世紀植被凈初級生產力被提出之后便受到了全球學者的廣泛關注[3-5]。草地生態系統是全球分布最廣的生態系統，準確估算其NPP對分析草地植被在全球氣候變化中的生態價值及全球碳循環收支平衡具有重要意義[6]。遙感技術是從空間上估算NPP的有效手段[7]，其中CASA模型應用最為廣泛，它是一種充分考慮環境條件和植被本身特征的光能利用率模型，可從不同空間尺度估算NPP[8-10]。CASA模型在估算NPP時需要輸入遙感、氣象、輻射和土地利用類型等數據，而氣象數據則需經過空間插值才可與遙感空間數據匹配，但不同插值方法獲取的氣象要素空間數據精度不一，這對CASA模型的運算產生一定的影響，進而影響到NPP的估算精度。

氣象要素的空間插值以盡可能建立接近實際空間分布特征的數學模擬方程為核心，根據區域內可得的氣象要素樣本值以一定的插值函數模型進行模擬，獲取氣象要素空間柵格數據[11]。目前應用較多的插值方法有反距離權重(IDW)、克里金、最近鄰、樣條函數和Anusplin等插值法，現有研究表明，不同的插值方法存在著各自的優劣[12-16]。也有學者對不同氣象插值方法在草地NPP估算中的可靠性進行了評價，但不同類型草地NPP的估算精度可能對不同氣象插值方法的敏感度存在差異，特別是在草地類型多樣化的農牧交錯帶，因此需進一步研究氣象數據插值方法對不同類型草地NPP估算的影響。

位于農牧交錯帶的寧夏有天然草地3.01萬km2，依靠草地發展的草畜產業是寧夏今后打造國家級生態草牧業試驗示范區的基礎，而準確估算草地NPP及草地載畜能力對制定生態草牧業發展戰略具有重要的指導意義。為此，本研究根據寧夏境內及周邊14個氣象站點2000—2015年的氣象數據，對利用反距離權重(IDW)、樣條函數及Anusplin這3種主流插值法獲取的寧夏氣象要素空間插值數據進行對比，利用交叉檢驗法對插值精度進行驗證，并探討氣象要素不同插值方法對不同類型草地NPP估算精度的影響，采用實測數據及MOD17A3 NPP數據進行對比驗證，以期優選出最適宜于寧夏地區NPP估算的氣象要素空間插值方法。

1 數據與方法

1.1 研究區概況與數據

寧夏位于104°17′—107°39′E，35°14′—39°23′N，全區屬典型溫帶大陸性氣候，年降水量183～677 mm，主要降水集中于7月、8月，全區總面積5.18萬km2，南北狹長，地勢南高北低。寧夏草原面積廣闊，以干草原、荒漠草原為主，主要分布在中部風沙干旱區的鹽池、靈武、同心等縣，另外，還有賀蘭山、南華山、西華山、月亮山、六盤山和云霧山等地的山地草原，天然草原是寧夏生態系統的重要組成部分，對構建黃河中游上段的生態保護屏障意義重大。

氣象數據包括全區及周邊省區在內14個氣象站點(圖1)2000—2015年的月平均氣溫、月總降水量以及太陽總輻射數據，這些數據均來自中國氣象數據網(http:∥data.cma.cn/)。遙感影像數據包括MOD13A3 NDVI月值數據及MOD17A3 NPP年值數據，其空間分辨率均為1 km，均來自于NASA數據分發網站(https:∥ladsweb.modaps.eosdis.nasa.gov/)，采用MRT軟件進行格式轉換及投影坐標系的轉換。數字高程模型(DEM)為地理空間數據云(http:∥www.gscloud.cn/)發布的90 m分辨率的SRTM數據。草地類型圖由1∶1 200 000的紙質草地分類圖矢量化所得，寧夏草地共分為10種類型。以上所有數據均轉換成橫軸莫卡托投影，空間分辨率均重采樣為1 km。草地初級生產力實測數據為寧夏15個市縣1981—2010年監測到的所屬市縣各類型草地凈初級生產力的加權平均值，數據來自文獻[17]。

MOD17A3 NPP數據是基于BIOME-BGC模型的覆蓋全球的NPP數據。該模型屬于生態過程模型，通過模擬生態系統內的光合、呼吸等生理活動及植物組織的營養物質的傳遞與循環等生理生態過程估算植被總初級生產力。在BIOME-BGC模型中植被NPP是植被總初級生產力與植被呼吸消耗之差。模型輸入參數分為初始化文件、氣象數據及生態生理指標，與CASA模型相比，其中加入了土壤有效深度、土壤顆粒組成、大氣CO2濃度年際變化及冠層消光系數、葉氮在羧化酶中的百分比含量等眾多植被生理指標，但由于部分模型參數獲取困難，因此在中小尺度區域的NPP估算中實現較困難[18-19]。而該模型在國內外認可度高，MOD17A3 NPP數據在全球范圍內廣泛應用于不同生態系統的碳循環研究，因此本文以全球MOD17A3 NPP數據為基準，對比不同插值算法下的CASA模型的估算精度[20-24]。

圖1 寧夏氣象站點及草地類型分布

1.2 插值方法

1.2.1 反距離權重法(IDW) 反距離權重法(Inverse Distance Weighted, IDW)認為插值點的權重與樣點的距離成反比，即以插值站點和樣本點之間的距離作為權重進一步加權平均，樣點與插值點的距離越近，其實測值對插值點的影響越大。反距離權重是對距離進行加權平均，因此樣點越密集，其模擬效果越好，插值效果最佳，算法如下[25]：

(1)

式中：Z為待模擬的插值點的柵格值；Z(xi)為第i(i=1,2,3,…,n)個氣象站的實測值；n為樣點數；di為插值點到第i個站點的距離。

1.2.2 樣條函數法(Spline) 樣條函數法采用最小化表面總曲率的數學函數來模擬未知點，即通過實測樣點生成恰好經過輸入點的平滑表面，通常有規則樣條函數和張力樣條函數兩種方法。本研究采用規則樣條函數方法(Regularized)進行氣象數據插值，采用此方法權重越大，擬合表面越光滑，算法如下[26]：

(2)

(3)

T(x,y)=a1+a2x+a3y

(4)

式中：Z為待模擬的插值點的柵格值；n為樣點數；λi為線性方程確定的系數；γi為模擬插值點到第i點的距離；τ2為權重系數；k0為修正貝塞爾函數;c=0.577215；a為線性方程系數。以上兩種插值方法均在ArcMap軟件中進行模擬。

1.2.3 Anusplin插值法 Anusplin插值法基于普通薄盤及局部薄盤樣條函數插值理論，利用集成Anusplin軟件包進行空間插值，該插值法可以將海拔等協變量引入模型以提高插值精度，算法如下[27]：

Zi=f(xi)+bTyiρ+ei(i=1,…,N)

(5)

式中：Zi表示空間i點的因變量；xi為d維樣條獨立變量;f(xi)是模擬xi的未知光滑函數；yi為ρ維獨立協變量；b是yi的ρ維系數；ei是自變量隨機誤差，其期望值為0。式中的函數f(xi)及b系數通過最小二乘法估計得出：

(6)

式中：Jm(f)是測度函數，用來監測函數f(xi)的粗糙度，其定義為函數f(xi)的m階偏導，在Anusplin插值法中的應用為樣條次數，在本研究中經過試驗選擇2次樣條進行插值；ρ為正的光滑參數，作為數據保真與曲面粗糙度之間的平衡[28]。

1.3 誤差檢驗

采用交叉檢驗進行插值精度誤差檢驗，即選取總樣本中30%的站點為檢驗樣點，剩余的作為訓練樣點進行空間插值，本研究選取寧夏境內的陶樂、中寧、鹽池及固原4個站點作為檢驗樣點(圖1)。假設檢驗站點的值未知，先通過插值算法獲取其模擬值，再通過計算檢驗站點實測值與模擬值之間的平均誤差(MAE)、相對平均誤差(MRE)和均方根誤差(RMSE)等指標來評估插值效果，具體公式如下：

(7)

(8)

(9)

式中：n為檢驗樣本站點的個數；Zai為第i個站點的實測值；Zλi為估測值。

1.4 NPP估算模型

CASA(Carnegie-Ames-Stanford Approach)模型建立在植物光合作用過程和光能利用率的概念上，是一個充分考慮環境條件和植被特征的NPP估算模型，計算公式如下[5]：

NPP=APAR×ε

(10)

式中：APAR為光合有效輻射，ε為光能轉化率。植被吸收的光合有效輻射APAR由太陽總輻射(SOL)及光合有效輻射吸收比率(FPAR)估算所得；FPAR由以歸一化植被指數(NDVI)為基礎的FPARNDVI和以比值植被指數(SR)為基礎的FPARSR加權平均獲得；光能轉化率ε由溫度脅迫因子、水分脅迫因子及最大光能利用率估算獲得，本研究參考朱文泉等[8]的研究結果，即草地的光能轉化率ε取值為0.542。

2 結果與分析

2.1 氣象要素空間插值結果與檢驗

2.1.1 氣象要素空間插值結果 采用3種不同的插值方法空間化的寧夏年均溫為-3～13℃，氣溫北高南低，高值均出現在寧夏中北部地區，而西南部為低值區(圖2—3)。3種插值方法都能模擬出寧夏氣溫的基本空間分布特征，均表現出由北向南降水量遞增的趨勢。由于Anusplin插值法將高程作為協變量，根據氣溫直減率和2次樣條函數進行氣溫空間插值，故其插值結果在局部特征上比樣條函數和IDW插值法的結果更為細膩，即地溫隨海拔起伏而變化明顯，特別是南部六盤山、南華山、中部羅山以及北部賀蘭山的高海拔區，其年均溫較低；而在海拔較低的北部引黃灌區、清水河河谷等地，年均氣溫形成高值區，這與寧夏的實際情況較為相符。而樣條函數及IDW插值法插值結果只反映出了寧夏年均氣溫自南至北的梯度變化特征，將賀蘭山山區插值成與引黃灌區相近的溫度特征，這與實際情況不符。此外，IDW插值法獲取的氣溫空間分布還在銀川、同心、海原等氣象站點附近出現了“牛眼效應”，即以氣象站為中心出現區域高值或低值中心。

3種不同的插值方法獲取的年總降水量空間變化特征基本相近，均表現出由北向南降水量遞增的趨勢，Anusplin插值法及樣條函數插值法獲取的降水量空間上由北向南遞增梯度基本一致，形成了有規律的降水遞增梯度線，而IDW插值法仍然出現了較為明顯的“牛眼效應”。

圖2年均溫在不同插值法下的空間分布

2.1.2 插值精度檢驗 為評估3種插值方法的精度，研究選擇10個氣象站數據進行插值，選擇4個氣象站的數據進行驗證。通過不同插值算法空間氣象要素估算，在ArcGIS中提取4個驗證站點的插值結果，通過平均誤差(MAE)、相對平均誤差(MRE)和均方根誤差(RMSE)等誤差分析指標來對比其插值精度(表1)。從擬合度來看，Anusplin和樣條函數的插值結果明顯優于反距離權重插值結果，而Anusplin的插值精度又略高于樣條函數插值結果。3種插值法的空間插值效果及插值精度交叉檢驗結果顯示，Anusplin插值法在氣溫與降水的空間模擬值的相對平均誤差(MRE)、平均誤差(MAE)及均方根誤差(RMSE)均小于其他兩種插值方法，其模擬精度較高。

2.2 CASA模型估算結果及驗證

2.2.1 草地NPP估算結果 在利用CASA模型估算區域NPP時，需要輸入空間氣象要素作為模型驅動變量，然而前文已述及不同插值方法獲取的氣象要素空間插值結果精度不同，每種插值方法對NPP的估算結果有何影響，需要進一步深入研究。為此，本節選擇以插值誤差最小的Anusplin插值法獲取的空間氣象要素和以IDW插值法獲取的空間氣象要素分別來驅動CASA模型，并進行NPP估算，并將二者估算的NPP結果與美國航空航天局(NASA)發布的全球MOD17A3 NPP數據進行對比分析，評估不同氣象插值方法對草地NPP估算的影響。寧夏草地NPP估算結果如圖4所示，兩種不同氣象插值數據驅動下的NPP估算值與MOD17A3 NPP數據在空間趨勢上基本一致，均表現出南部山區NPP較高，北部NPP較低，而中部干旱帶的NPP最低，這與寧夏草地的分布格局相符合。寧夏南部的六盤山、南華山地區主要以山地草原為主，黃土丘陵區主要以典型草原為主，加之南部山區降水較北部豐沛，故草地NPP為全區最高，而中部干旱帶則以干草原、荒漠草原為主，其NPP也自然最低。基于Anusplin插值法的CASA模型與基于IDW插值法的CASA模型草地NPP估算值在量級上較為接近，草地均值分別為149.42 gC/(m2·a)與150.45 gC/(m2·a)，而MOD17A3 NPP數據的草地NPP均值為147.65 gC/(m2·a)，略低于CASA模型估算值。在草地NPP值域范圍上，基于Anusplin插值氣象要素估算的NPP值為45.06～807.83 gC/(m2·a)，而基于IDW插值氣象要素估算的NPP的值介于47.16～733.63 gC/(m2·a)，二者的值域范圍均大于MOD17A3 NPP數據的值域范圍47.28～586.66 gC/(m2·a)，特別是NPP像元最大值相差較大。

圖3 年降水量在不同插值法下的空間分布表1 氣溫降水量插值誤差對比

為對比不同氣象插值方法在不同類型草地NPP估算中的應用效果，本研究分別計算不同插值方法驅動下的寧夏全區各草地類型的平均NPP(圖5)。從3種模型估算的NPP平均值來看，兩種基于CASA模型估算的NPP結果相近，與MOD17A3 NPP數據存在較大差異。在灌叢草原、低濕地草甸和山地草甸中，MOD17A3 NPP數據明顯比CASA模型估算的NPP偏低；在荒漠草原、荒漠化類草原、干荒漠類、草甸草原、灌叢草原和沼澤類草地中，MOD17A3 NPP數據明顯比CASA模型估算的NPP偏高；而3種模型對干草原的NPP估算中非常接近。從NPP值來看，荒漠草原類的草地NPP最低，而山地草甸類的草地NPP最高。

2.2.2 與全區實測數據對比 基于兩種不同氣象插值方法和CASA模型估算了近16 a寧夏草地年平均NPP，并用寧夏境內15個縣市多年NPP實測均值進行對比。草地實測值為地上草地產草量干重，而CASA模型的估算值NPP是地上和地下生物量的總和，因此在驗證中參照前人文獻中有關干物質產量到NPP的轉換關系及不同草地類型地下與地上生物量的比例系數的研究[29-30]，將CASA模型估算的草地地上NPP值求出，并將其與草地實測地上生物量進行對比(圖6)。基于不同氣象插值方法驅動CASA模型估算的NPP值均與實測NPP存在良好的線性關系，且其相關系數R2均在0.82以上，說明CASA模型在寧夏草地NPP估算中效果較為理想，實用性較強。基于Anusplin插值法估算的草地NPP值與實測值相關性達0.86，優于IDW插值法估算的草地NPP值，說明Anusplin氣象插值法更適宜于驅動CASA模型。此外，對比兩種插值方法在NPP估算中的誤差發現，基于Anusplin插值法估算的草地NPP其相對平均誤差(MRE)、均方根誤差(RMSE)分別為0.23，2.35，小于基于反距離權重插值法估算的草地NPP的誤差(MRE=0.23，RMSE=2.37)，這說明基于Anusplin插值獲取的氣象要素驅動的CASA模型所估算NPP的精度略高，更接近實測數據，能夠反映寧夏實際的草地NPP分布狀況。

圖4寧夏草地NPP空間分布

圖5不同模型模擬的寧夏各草地類型平均NPP對比

2.2.3 不同草地類型NPP估算結果對比 以全球MOD17A3 NPP數據為基準，對比分析CASA模型的估算誤差，對CASA模型在寧夏草地NPP估算中的可靠性做進一步的研究。通過提取不同NPP圖像中的對應像元值，逐像元制作散點圖(圖7)，對比發現基于Anusplin和IDW插值法獲取的氣象要素空間數據驅動CASA模型而估算的草地NPP值均高于MOD17A3 NPP數據值，線性擬合斜率分別為1.34，1.30，與MOD17A3 NPP數據的擬合優度(R2)分別為0.83，0.82，即基于兩種氣象插值方法驅動的CASA模型在寧夏草地NPP估算中的結果相近。二者的擬合優度(R2)僅僅相差不到0.01，盡管優勢微弱，但結果仍然顯示基于Anusplin插值法獲取的氣象要素空間數據驅動CASA模型的效果要好于基于IDW插值的氣象要素空間數據。由此可以看出，在寧夏草地NPP估算中，選擇精度較高的Anusplin氣象要素插值法，能夠提高草地NPP的估算精度。

基于以上研究，選擇估算精度整體較高的基于Anusplin插值的CASA模型與基于BIOME-BGC模型在不同草地類型的NPP估算中進行深入對比研究。研究分析了2000—2014年基于CASA模型的NPP估算結果與MOD17A3 NPP數據間的誤差(表2)。總體來看，10類草地的總體相關性分析結果顯示，基于Anusplin插值氣象要素估算的NPP與MOD17A3 NPP數據的相關系數為0.88(p<0.01，n=150)。從不同草地類型上來看，在干草原、灌叢草原、干荒漠類草原及荒漠草原的估算中CASA模型與MOD17A3產品值相關系數最高，在沼澤類草原估算中效果欠佳，相關系數未通過顯著性檢驗，其他類型草地相關性整體在0.7以上，均通過了p<0.01的顯著性檢驗。從估算誤差分析來看，CASA模型在干草原的估算誤差最小，荒漠草原估算誤差次之，這兩類草原總面積接近寧夏草原面積的80%，代表了寧夏大部分的草地生產力；在沼澤類和山地草甸的估算誤差最大，這兩類草地在寧夏面積較小，在區域尺度的模型估算中產生較大誤差屬正常情況。

圖6寧夏草地多年年均NPP模擬值與實測值相關性

圖7不同插值算法NPP模擬逐像元值散點圖

表2 不同類型草地CASA模型估算NPP與MOD17A3 NPP數據的誤差對比

注：**為p<0.01，*為p<0.05，n=15。

3 結 論

(1) 在寧夏這樣的氣象站點稀疏地區，除反距離權重插值法出現了“牛眼效應”之外，其他插值法均能較好模擬寧夏自南向北的空間分布特征；交叉檢驗顯示，Anusplin氣象要素插值相比于傳統插值算法明顯提高了氣象數據插值精度，其中氣溫插值表現最為明顯。

(2) 采用基于不同氣象要素插值方法的CASA模型對寧夏草地NPP進行估算，其估算結果均能反映寧夏草地NPP的空間分布格局；在與多年草地NPP實測值及MOD17A3 NPP數據草地估算值的對比驗證中發現，在全區草地NPP估算中，基于Anusplin插值的CASA模型的模擬值與實測值相關性較高，其模擬值可靠，具有一定的科學性。

(3) 對模型與MOD17A3 NPP數據對比研究顯示，插值精度較高的Anusplin插值法驅動CASA模型估算的NPP與MOD17A3 NPP數據的相關性高，因此，提高氣象要素的插值精度在一定程度上能提高CASA模型NPP的估算精度。

(4) 以MOD17A3 NPP數據為真值檢驗CASA模型在寧夏不同草地NPP估算的估算精度分析表明，CASA模型在干草原、灌叢草原、干荒漠類草原及荒漠草原的估算中精度較高，而在沼澤類和山地草甸的估算效果欠佳，其估算精度有待提高。

在CASA模型估算精度的對比中，由于實測數據為各個市區的多年草地NPP均值的限制，無法從草地類型上進一步對比驗證，此為本研究的一大遺憾。但為補充驗證數據源在研究中引入了MOD17A3 NPP數據進行對比驗證，前提假設MOD17A3 NPP數據不存在誤差，這是一個有限條件假設。因此在模型驗證上需注重實測數據的收集，以在后期研究中豐富模型的可行性分析。此外，寧夏草地的時空格局變遷及人類活動、自然因素與NPP波動之間的內在聯系需進一步做重點探討。