基于CASA模型和SEVI指數(shù)的福建省植被NPP遙感估算與分析

江 洪，虞嘉瑋，蔣世豪，黃貝瑩，李玉潔

（1.福州大學(xué)數(shù)字中國研究院，福建 福州 350108；2.福州大學(xué)空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實驗室，福建 福州 350108；3.福建省基礎(chǔ)地理信息中心，福建 福州 350003）

植被凈初級生產(chǎn)力（Net Primary Productivity，NPP）是指綠色植被在單位時間、單位面積上通過光合作用同化CO2并扣除自養(yǎng)呼吸消耗后的剩余有機(jī)物［1］，是評估陸地生態(tài)系統(tǒng)碳匯和調(diào)節(jié)過程的主要因子，在全球的碳平衡中發(fā)揮重要的作用［2］.目前，遙感反演植被NPP的模型主要有統(tǒng)計模型、過程模型和光能利用率模型.統(tǒng)計模型通常使用葉面積指數(shù)（Leaf Area Index，LAI）或者歸一化植被指數(shù)（Normalized Differential Vegetation Index，NDVI）進(jìn)行估算［3-4］.過程模型常見的有Biome-BGC［5］，BEPS［6］，TEM［7］，InTEC［8］等.光能利用率模型以太陽輻射為基礎(chǔ)，結(jié)合植被的吸收、轉(zhuǎn)化能力估算植被NPP，例如CASA［9］，SEBAL［10］，Nnpp［11］，C-Fix［12］，Glo-PEM［13］等模型.其中CASA模型綜合考慮了太陽輻射、水分、溫度、不同植被類型的影響.根據(jù)研究區(qū)的差異，通過參數(shù)的改進(jìn)，CASA模型可以較好地應(yīng)用于不同地區(qū).例如朱文泉［14］等利用誤差最小原則，根據(jù)中國的植被NPP實測數(shù)據(jù)模擬中國各種植被類型的最大光能利用率參數(shù)進(jìn)而估算中國的植被NPP.在山地丘陵等地形崎嶇的區(qū)域，遙感影像上光學(xué)輻射信號受地形影響會發(fā)生畸變，影響植被NPP估算精度.結(jié)合數(shù)字高程模型（Digital Elevation Model，DEM）數(shù)據(jù)對CASA模型中的參數(shù)進(jìn)行改進(jìn)可以達(dá)到一定的效果［15］，但是基于DEM的地形校正模型只對地形本影校正效果良好，對地形落影的校正效果欠佳［16］.近年提出的陰影消除植被指數(shù)（Shadow Eliminated Vegetation Index，SEVI）可以有效地消除地形本影和落影對山區(qū)植被信息的影響，并且無需DEM數(shù)據(jù)的支持［17］.因此，筆者將SEVI應(yīng)用到CASA模型，改進(jìn)其中的植被光合有效輻射吸收比率因子（Fraction of Photosynthetically Active Radiation，F(xiàn)PAR）計算模型，提高山區(qū)植被NPP遙感估算精度，并計算分析福建省2005年與2015年的植被NPP.

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況福建?。?3°33′～28°20′N、115°50′～120°40′E）地處我國東南沿海，下轄9個地級市和1個綜合實驗區(qū)（如圖1所示），陸地面積12.4×104km2，山地、丘陵占全省總面積的80%以上，地勢總體上西北高東南低.福建省是我國南方重點(diǎn)林區(qū)之一，氣候溫和，冬短夏長，降水充沛，自然條件優(yōu)越，植被種類多達(dá)5 000余種，森林覆蓋率達(dá)66.8%［18］，居全國第一.

圖1 研究區(qū)與驗證樣本示意圖

1.2 數(shù)據(jù)源及預(yù)處理遙感影像采用2005年和2015年的MODIS數(shù)據(jù)，包括空間分辨率250 m的MOD13Q1數(shù)據(jù)和空間分辨率500 m的MCD12Q1數(shù)據(jù)和MOD17A3H數(shù)據(jù)（來源于美國國家航天局，https：//ladsweb.modaps.esodis.nasa.gov/）.MOD13Q1數(shù)據(jù)時間分辨率為16 d，用于CASA模型中光合有效輻射（APAR）的計算.MCD12Q1數(shù)據(jù)時間分辨率為1 a；用于光能利用率ε計算中植被生態(tài)系統(tǒng)的劃分.MOD17A3H數(shù)據(jù)為NPP產(chǎn)品，時間分辨率為1 a，用于本文植被NPP估算結(jié)果的對比.氣象數(shù)據(jù)采用福建省所有氣象站點(diǎn)的月值觀測數(shù)據(jù)，包括降水、氣溫、濕度、日照時數(shù)等（來源于中國氣象數(shù)據(jù)網(wǎng)，http：//data.cma.cn/）.

MODIS數(shù)據(jù)的預(yù)處理包括數(shù)據(jù)格式轉(zhuǎn)換、投影轉(zhuǎn)換以及裁剪，并利用最大值合成法將年內(nèi)多景MOD13Q1數(shù)據(jù)合成12個月的月值數(shù)據(jù).DEM數(shù)據(jù)的處理包括格式轉(zhuǎn)換、投影轉(zhuǎn)換、鑲嵌以及裁剪.氣象數(shù)據(jù)的預(yù)處理包括：數(shù)據(jù)的清洗，用往年的平均值代替異常值或者缺失值；數(shù)據(jù)插值，采用克里金插值法將站點(diǎn)氣溫、降水、濕度等點(diǎn)數(shù)據(jù)轉(zhuǎn)換成柵格數(shù)據(jù).

2 研究方法

CASA模型主要通過光能利用率ε和植被吸收的光合有效輻射（Absorbed Photosynthetic Active Radiation，APAR）計算植被NPP，

其中，x表示像元，t表示時間；N PP(x，t)表示像元x在t月份的凈初值生產(chǎn)力，單位：g·m-2；ε(x，t)表示像元x在t月份的實際光能利用率，單位：g·MJ-1；APA R(x，t)表示像元x在t月份獲得的光合有效輻射，單位：MJ·m-2.

2.1 光能利用率的計算光能利用率指植被將實際吸收的太陽有效輻射轉(zhuǎn)換成自身有機(jī)物的效率，通過溫度脅迫因子T、水分脅迫因子W以及最大光能利用率εmax計算獲得，見式（2）.溫度脅迫因子的計算見式（3）和（4），水分脅迫因子的計算見式（5），最大光能利用率參數(shù)如表1所示［14］.

表1 最大光能利用率參數(shù)表 %

其中，ε(x，t)表示像元x在t月的光能利用率，Tε1(x，t)表示植被在低溫與高溫情況下受到的限制作用，由式（3）計算獲得；Tε2(x，t)表示植被所處的環(huán)境溫度，由公式（4）計算獲得；Wε(x，t)表示水分條件對光能利用率的影響，由公式（5）計算獲得；εmax表示在理想條件下植被轉(zhuǎn)換成A P AR的效率.

其中，Tε1(x，t)表示像元x在t月的第一溫度脅迫因子；T opt(x)表示像元x年內(nèi)生長的最適溫度，當(dāng)像元x在t月的溫度小于-10℃時，Tε1(x，t)=0［19］.

其中，Tε2(x，t)表示像元x在t月的第二溫度脅迫因子；T(x，t)表示像元x在t月的溫度，比最適溫度高10℃或者低13℃時，Tε1(x，t)等于最適溫度T o pt(x)時Tε2(x，t)的一半［19］.

其中，Wε(x，t)表示像元x在t月的水分脅迫因子；E ET(x，t)表示像元x在t月的實際蒸發(fā)量，采用區(qū)域?qū)嶋H蒸散模型求?。?0］；P ET(x，t)表示像元x在t月的潛在蒸散量，根據(jù)Thornthwaite植被-氣候關(guān)系模型求取［21］.

2.2 光合有效輻射的計算光合有效輻射由太陽輻射總量（SOL）以及植被光合有效輻射的吸收比率（FPAR）決定，

其中，SO L(x，t)表示像元x在t月的太陽輻射總量，單位MJ·m-2，由天文輻射總量和日照百分率線性擬合而成［22］，

其中，Q S為月太陽輻射總量；Q0為月天文輻射總量，起伏地形下的天文輻射總量通過DEM數(shù)據(jù)計算獲得，見式（8）和（9）［23］；s為日照百分率；a，b為經(jīng)驗系數(shù)，引用陳惠［24］等利用實測數(shù)據(jù)擬合得到的福建省的經(jīng)驗系數(shù)，如表2所示.

表2 各月經(jīng)驗系數(shù)a和b

其中，Q是地形起伏下日天文輻射量，單位：MJ·m-2；E sc是太陽常數(shù)；r0/r是當(dāng)天日地距離訂正系數(shù)；n為可照時角的離散數(shù)目，δ為太陽赤緯角；ωsi和ωri分別是可照時間段的起始和終止太陽時角（°）.

其中，u、v、w是與地形、地理有關(guān)的特征因子，φ、α、β分別為當(dāng)?shù)氐木暥?、坡度和坡?

2.3 FPAR改進(jìn)植被光合有效輻射吸收比率（FPAR）是CASA模型估算植被NPP過程中主要的驅(qū)動因子.目前基于植被指數(shù)反演FPAR的方法眾多［25］，常用的植被指數(shù)有NDVI和比值植被指數(shù)（RVI）等，由NDVI反演FPAR的計算公式

其中，F(xiàn)PA Rmax和F PA Rmin的取值分別為0.95和0.001，ND V I(i，max)和ND V I(i，min)的取值分別為第i種地類N DV I的95%及5%的百分位數(shù)［14］.

由于山區(qū)植被受地形陰影影響顯著，根據(jù)文獻(xiàn)［26］研究結(jié)果表明：未經(jīng)地形校正的NDVI反演得到的FPAR在本影和落影的數(shù)值遠(yuǎn)低于相鄰非陰影陽坡，相對陽坡的誤差大于70%；經(jīng)C模型或者FLAASH+C模型校正的NDVI反演得到的FPAR，本影的相對誤差小于10%，但是落影的相對誤差仍超過24%；通過SEVI反演的FPAR，本影、落影與陽坡的數(shù)值十分接近，相對誤差小于3%.因此采用能有效消除地形陰影影響的SEVI替代NDVI，改進(jìn)FPAR的反演模型

其中，SE VI(i，max)和SE V I(i，min)的取值分別為第i種地類SE VI的95%及5%的百分位數(shù).

SEVI計算公式［17］

其中，ρnir和ρred分別代表近紅外波段和紅光波段的反射率，f(Δ)為地形調(diào)節(jié)因子，通過相關(guān)系數(shù)法進(jìn)行優(yōu)化選取.以0.001為步長，令f(Δ)從0到1.000，計算f(Δ)對應(yīng)下的SEVI與RVI的相關(guān)系數(shù)r1，SEVI與SVI的相關(guān)系數(shù)r2，當(dāng)r1和r2接近或相等時，確定最優(yōu)的調(diào)節(jié)因子（f opt），計算公式

其中，x，y1，y2分別為遙感影像SEVI，RVI以及SVI的計算結(jié)果，n為SEVI，RVI以及SVI影像的像元數(shù).

2.4 FPAR結(jié)果驗證分別利用MODIS數(shù)據(jù)計算的NDVI和SEVI反演FPAR，選取相鄰陰陽坡的FPAR樣本，統(tǒng)計分析陰影區(qū)域和對應(yīng)陽坡區(qū)域的FPAR均值以及相對誤差，驗證不同植被指數(shù)計算的FPAR精度.相對誤差的計算公式

其中，F(xiàn)P A Rshadow為陰影處的植被光合有效輻射吸收比率，F(xiàn)PA Rsunny為相鄰非陰影陽坡的植被光合有效輻射吸收比率.

250 m空間分辨率的MODIS影像相較于文獻(xiàn)［26］中使用的30 m空間分辨率的Landsat8 OLI數(shù)據(jù)，難以區(qū)分本影和落影，但在地形復(fù)雜的山區(qū)仍有可以目視觀察的陰影.根據(jù)地理學(xué)第一定律，一般情況下，地物間的距離越接近，其相關(guān)性也越大；同時，也有相關(guān)研究表明地處亞熱帶的福建省山地的陰坡水熱條件充沛，陰坡與陽坡的天然林長勢狀態(tài)極為接近［27-28］.因此，在Google Earth影像參考下，通過目視解譯，手動選取了45組相鄰陰陽坡植被長勢接近的樣本，如圖1所示.

3 結(jié)果與分析

3.1 FPAR估算結(jié)果估算結(jié)果顯示采用NDVI反演的FPAR會造成陰影區(qū)域的數(shù)值偏低（如圖2a所示），而采用SEVI反演的FPAR在陰陽坡的數(shù)值接近（如圖2b所示）.樣本FPAR統(tǒng)計結(jié)果顯示采用NDVI反演的FPAR在陰影區(qū)域的相對誤差為7.85%，而采用SEVI反演的FPAR在陰影區(qū)域的相對誤差降至0.53%（如表3所示）.這說明采用MOD13Q1數(shù)據(jù)反演復(fù)雜地形山區(qū)地表參數(shù)時，有必要進(jìn)行地形校正，這與文獻(xiàn)［26］中基于Landsat8 OLI數(shù)據(jù)研究得到的結(jié)論一致.

圖2 樣本F P A R折線圖

表3 FPAR在陰影、陽坡樣本處的均值M以及相對誤差E

3.2 植被NPP估算結(jié)果2005年福建省植被NPP估算結(jié)果平均值為861.9 g·m-2·a-1，不同植被類型NPP均值大小順序依次為常綠闊葉林、針闊混交林、常綠針葉林、灌叢、草原和農(nóng)用地，均值分別為1 082 g·m-2·a-1、908.2 g·m-2·a-1、863.7 g·m-2·a-1、852.7 g·m-2·a-1和711.2 g·m-2·a-1.2015年福建省植被NPP均值為855.7 g·m-2·a-1，不同植被類型年均NPP均值依次為1 107.9 g·m-2·a-1、875.5 g·m-2·a-1、831.7 g·m-2·a-1、796.7 g·m-2·a-1、764.8 g·m-2·a-1和758.0 g·m-2·a-1，大小順序與2005年保持一致.

由于植被NPP實測結(jié)果較難獲取，采用交叉驗證法分析本文估算結(jié)果的合理性.查閱文獻(xiàn)表明已有多位學(xué)者對2005年福建省植被NPP進(jìn)行了估算（如表4所示），均值范圍在788.3～811.5 g·m-2·a-1，常綠闊葉林、針闊混交林、常綠針葉林、灌叢、草原和農(nóng)用地的NPP范圍分別為710.7～1 600 g·m-2·a-1、639.0～1 100 g·m-2·a-1、1 022.2～1 300 g·m-2·a-1、800～900 g·m-2·a-1、700～800 g·m-2·a-1、600～700 g·m-2·a-1.基 于MOD17A3H數(shù)據(jù)估算的2001～2010中國東南部植被NPP均值為511.4 g·m-2·a-1，其中福建省的年均植被NPP范圍在500～1 000 g·m-2·a-1，不同植被類型NPP大小順序為常綠闊葉林、針闊混交林、常綠針葉林、草地、農(nóng)田植被［29］.福建、廣東、廣西等長江以南各省區(qū)，常綠闊葉林實測NPP范圍在910～1 340 g·m-2·a-1［30］.利用CASA模型估算的武夷山區(qū)一帶亞熱帶森林年均NPP超過500 g·m-2·a-1［31］.BEPS模型估算的2004年福建省森林生態(tài)系統(tǒng)的NPP全年均值為578.97 g·m-2·a-1，闊葉林為780.0 g·m-2·a-1［32］.本文估算結(jié)果高于2004年BEPS模型估算的結(jié)果，主要是由于模型不同、年份不同，以及本文研究涵蓋的植被類型更多，因此，結(jié)果更高.對比其他學(xué)者2005年的研究結(jié)果，除針闊混交林和灌木外，本文估算的各植被類型NPP結(jié)果與已有研究結(jié)果吻合.本文估算的福建省植被NPP結(jié)果高于Biome-BGC模型的估算結(jié)果［33］，也略高于其他學(xué)者的CASA模型估算結(jié)果［33-35］，可能的原因是SEVI相對NDVI提高了FPAR的估算值，進(jìn)而提高了植被NPP的估算結(jié)果.

表4 2005年福建省植被NPP估算結(jié)果比較 g·m-2·a-1

3.3 植被NPP空間分布特征福建省年均植被NPP的空間分布總體呈現(xiàn)內(nèi)陸向沿海遞減的趨勢，如圖3所示.按照行政區(qū)劃統(tǒng)計各地級市年均植被NPP，結(jié)果與圖3相似（如表5所示）：2005年植被NPP最大的為龍巖市，其次為三明市、南平市、寧德市、漳州市、福州市、莆田市、泉州市和廈門市；2015年全省植被NPP第一、二位和第八、九位保持不變，第三到七位分別為漳州市、南平市、寧德市、莆田市和福州市.

圖3 福建省2005年和2015年年均植被NPP空間分布

表5 不同行政區(qū)的年均植被NPP及變化

以100 m為步長，對福建省0～1 800 m海拔劃分為18個等級，統(tǒng)計不同海拔等級植被NPP.結(jié)果顯示，植被NPP隨海拔變化的趨勢在2005年和2015年基本一致.低海拔區(qū)域（500 m以下）植被NPP數(shù)值較低，隨海拔升高NPP逐步增加；中海拔區(qū)域（500～1200 m）植被NPP達(dá)到峰值并穩(wěn)定，在700 g·m-2·a-1附近小幅度波動；高海拔區(qū)域（1 200 m以上）植被NPP逐步下降，如圖4所示.這一結(jié)果表明，福建省植被NPP的空間分布與海拔高度相關(guān)性強(qiáng)，沿海城市的平均海拔處于第一級，人類活動密集，森林覆蓋率低，因此平均植被NPP較低；而內(nèi)陸城市的平均海拔高度處于第二級到第四級，森林覆蓋率高，人類對森林的干擾較少，因此平均植被NPP較高.

圖4 不同海拔等級對應(yīng)的植被NPP

3.4 植被NPP月度分布特征福建省月均植被NPP與溫度、降水、日照百分率3個水熱因子統(tǒng)計結(jié)果表明福建省植被NPP月平均值在時間序列上的分布呈現(xiàn)倒“V”型特征，如圖5所示.2005年1～2月，全省植被NPP均值在20 g·m-2以下；從3～7月，植被NPP均值總體上升，7月達(dá)到最大值148.95 g·m-2；8～12月，全省植被NPP均值穩(wěn)步下降，12月再次低至20 g·m-2·以下.2015年福建省植被NPP變化趨勢與2005年一致，1～6月份的植被NPP均值高于2005年同月的結(jié)果，8月份全省植被NPP達(dá)到最大值123.84 g·m-2·，9～12月植被NPP均值逐漸減少.

圖5 2005年和2015年月均植被NPP、溫度、降水和日照百分率

2005年月平均溫度1月最低，為8.01℃，7月出現(xiàn)最大值27.95℃，總體呈現(xiàn)先逐步上升，后緩慢下降的趨勢.月均降水在6月出現(xiàn)最大值401.41 mm，12月出現(xiàn)最小值23.81 mm，上半年波動明顯，下半年逐步下降.月均日照百分率也出現(xiàn)明顯的波動特征，7月出現(xiàn)最大值58.65%，2月出現(xiàn)最小值11.82%.與2005年相比，2015年月平均溫度走勢保持不變，月均降水和日照百分率的走勢區(qū)別較大，月均降水在8月出現(xiàn)最大值350.98 mm，1月出現(xiàn)最小值47.20 mm.日照百分率在10月出現(xiàn)最大值41%，12月出現(xiàn)最小值17.29%.綜合月均植被NPP和水熱因子的走勢，可以看出植被NPP月度分布特征與溫度因子的走勢最為相似，相關(guān)性分析結(jié)果顯示，二者的相關(guān)系數(shù)在2005年和2015年分別達(dá)到0.96和0.95.

3.5 2005～2015年植被NPP變化2005～2015年福建省植被NPP變化結(jié)果顯示，福建省植被NPP總體呈現(xiàn)略有減少的趨勢.結(jié)合圖6和表5可以看到，植被NPP減少的城市主要為南平市、寧德市和三明市，其中寧德市減少27.45 g·m-2·a-1，減幅達(dá)3.54%.植被NPP增加的為福州市、莆田市、泉州市、廈門市、漳州市和龍巖市，其中漳州市增加117.95 g·m-2·a-1，增幅達(dá)15.78%.

圖6 2005～2015年福建省植被NPP變化結(jié)果

4 討論

采用SEVI指數(shù)改進(jìn)CASA模型中FPAR參數(shù)的反演，降低了山區(qū)陰影和陽坡的FPAR相對誤差.利用改進(jìn)的CASA模型估算福建省2005年和2015年的植被NPP，估算結(jié)果略高于其他學(xué)者基于NDVI反演的結(jié)果，可能是由于SEVI降低了崎嶇山區(qū)地形效應(yīng)的影響從而提高了FPAR的估算結(jié)果，進(jìn)而提高植被NPP的反演結(jié)果.

福建省植被NPP的空間分布特征在2005年和2015年保持了相同的趨勢，龍巖市、南平市和三明市均高于800 g·m-2·a-1，而沿海地區(qū)的廈門市、泉州市和福州市均低于750 g·m-2·a-1，說明植被NPP的空間分布特征較為穩(wěn)定.這種空間分布特征可能與森林覆蓋率等因素有關(guān)，龍巖市、南平市和三明市森林覆蓋率較高，同時山區(qū)人類活動較少，植被生長狀態(tài)好.廈門市、泉州市和福州市地處福建省沿海，城鎮(zhèn)化水平高，人口密度大，植被覆蓋率較低.不同海拔對應(yīng)的植被NPP差異明顯，海拔較低時，人類活動明顯，對植被的生長影響較大，而在中海拔區(qū)域，較少的人類活動以及合適的地理環(huán)境均有利于植被的生長，植被NPP達(dá)到最大值，而過高的海拔不適宜植被的生長.因此，植被NPP在海拔1 200 m以上出現(xiàn)明顯的下降趨勢.

福建省植被NPP在2005年和2015年表現(xiàn)出相同的月序特征，夏季（6～8月）出現(xiàn)最大值，冬季（12～2月）出現(xiàn)最小值，主要是由于夏季平均氣溫高、日照充足、降水少，適合植被的生長，而冬季的水熱條件較差，植被生長緩慢，導(dǎo)致植被NPP數(shù)值較小.

2005～2015年福建省植被NPP總體呈現(xiàn)略有減少的趨勢.根據(jù)中科院地理科學(xué)與資源研究所公布的中國土地利用現(xiàn)狀遙感監(jiān)測數(shù)據(jù)，統(tǒng)計了福建省2005年到2015年土地利用類型的變化，發(fā)現(xiàn)有林地、灌木林、疏林地的占比均下降，而城鄉(xiāng)、工礦、居民用地占比均上升.由于林地的植被NPP高于建設(shè)用地的植被NPP，因此契合福建省植被NPP略有減少的趨勢.此外，統(tǒng)計空間分辨率500 m的MODIS NPP產(chǎn)品（MOD17A3H）也表明2005年到2015年福建省行政區(qū)域內(nèi)的植被NPP均值由925.91 g·m-2·a-1略降為922.59 g·m-2·a-1.還有其他相關(guān)研究也得出福建省近年來年均植被NPP呈現(xiàn)下降趨勢的結(jié)論［29，36-37］.

5 小結(jié)

采用SEVI指數(shù)改進(jìn)CASA模型中FPAR參數(shù)的計算，可以提高CASA模型對山區(qū)植被NPP的估算精度.福建省2005年和2015的植被NPP反演結(jié)果表明：2005～2015年福建省植被NPP出現(xiàn)小幅度的減少，但是各植被類型NPP以及時空分布特點(diǎn)保持不變，不同植被類型NPP均值排序依次為常綠闊葉林＞針闊混交林＞常綠針葉林＞灌叢＞草地＞農(nóng)用地.

福建省植被NPP空間分布上呈現(xiàn)內(nèi)陸高沿海低的特征，并且隨海拔升高呈現(xiàn)先增加后減少的特點(diǎn)，在500～1 200 m的中海拔區(qū)域達(dá)到最大值.植被NPP月度分布特征為7月和8月達(dá)到最大值，而12月和1月出現(xiàn)最小值，與月平均溫度的走勢相關(guān)性強(qiáng)，2005年和2015年月均植被NPP和月均溫度的相關(guān)系數(shù)分別為0.96和0.95.