基于校正的CASA模型NPP遙感估算及分析——以河西走廊為例

李傳華,曹紅娟,范也平,韓海燕,孫 皓,王玉濤

西北師范大學地理與環境科學學院, 蘭州 730070

植被凈初級生產力(Net Primary Productivity, NPP)是陸地生態系統碳循環的重要組成,也是揭示碳源/匯的關鍵環節,NPP估算不僅可以反映植被的物質生產能力,也是揭示大氣二氧化碳收支和氣候環境變化的重要指標[1-3]。

NPP估算模型較多,CASA模型其原型由Monteith提出,之后Potter和Field等[4]在Monteith基礎上改進,因結構簡單,參數可通過遙感獲得,廣泛應用于全球和區域尺度上的NPP估算[5-7]。CASA模型中需要輸入氣象數據,氣象數據一般來源于氣象站點,而氣象站點通常選擇在空曠平坦的地帶,因此,通過站點數據直接空間插值的氣象數據較難反應復雜地形地區的氣候狀況,在估算NPP時會產生較大的誤差。基于CASA模型估算NPP與氣候因子相關的參數有太陽輻射、溫度脅迫系數和水分脅迫系數。太陽輻射校正方法比較成熟[8]；溫度脅迫因子是通過氣溫得到的,大部分研究使用如Kriging、IDW和Spline等方法對氣溫進行空間插值,效果較差[9-10],有研究表明,在甘肅省海拔、經度和緯度對氣溫影響也很大[11]；水分脅迫系數在CASA模型中是實際蒸散量與潛在蒸散量的比值,這兩個變量是根據經驗公式計算,沒有考慮到地表植被差異,即相同的比值對不同植被類型的水分脅迫是一樣的,顯然,這種計算存在較大的不確定性。并且,原CASA模型中全球植被的最大光能利用率為0.389 gC/MJ[12],對于干旱區并不適用。目前,有很多學者都對CASA模型進行了校正[13-15],但對模型校正前后NPP估算的影響卻鮮有研究。

河西走廊位于中國西北地區,屬于干旱半干旱區,是國家絲綢之路經濟帶的重要通道,祁連山國家森林公園也位于此,區內景觀多樣,在自然環境和人類活動雙重作用下,該區生態系統非常脆弱,植被退化非常嚴重,荒漠化發展速度較快[16]。本文利用校正的CASA模型,以氣象數據、遙感數據、DEM數據以及土地利用數據等為基礎,針對CASA模型的以上缺點,利用海拔、經度和緯度對氣溫進行校正,并使用地表水分指數(Land Surface Water Index, LSWI)反演水分脅迫系數,最大光能利用率采用馮益明等[17]針對干旱區得到的不同植被類型上的相應值,估算了2015年河西走廊NPP。旨在提高基于CASA模型的NPP模擬精度,并著重分析校正前后的CASA模型對NPP估算的差異及原因,同時也為該區植物生產力評價,生態資產核算,以及生態補償機制建立提供科學參考。

1 研究區概況及數據

1.1 研究區概況

河西走廊地處甘肅省境內,位于黃河以西,南部被祁連山和阿爾金山、北部被合黎山、龍首山以及馬鬃山夾持,是古絲綢之路的交通要道,地理位置為37°17′—42°48′ N,92°12′—103°48′ E,行政區劃包括武威市、張掖市、金昌市、酒泉市和嘉峪關市(圖1)。河西走廊地勢南高北低,南部山地海拔大多在3000—3500 m以上；中部走廊平原海拔為1000—2600 m左右；北部山地斷續分布。河西走廊屬于溫帶大陸性氣候,光照強烈,年降水量大部分地區在35—200 mm,大部分地區的年蒸發量在1500 mm以上,年均溫5.8—9.3℃。該區域內植被類型主要有林地、耕地、草地和荒漠。河西走廊由三大流域組成,從東向西依次為石羊河流域、黑河流域和疏勒河流域,三大河流都發源于祁連山,區域內綠洲農業發達。

圖1 研究區位置及氣象站點空間分布圖Fig.1 Location of the study area and spatial distribution of meteorological stations

1.2 數據來源與預處理

1.2.1 遙感數據

本文2015年的遙感數據來自美國國家航空航天局(https://search.earthdata.nasa.gov/search),產品是MOD13Q1(從中選取歸一化植被指數——NDVI(Normalized Difference Vegetation Index),為16日合成數據,空間分辨率250 m)、MOD09A1(8天合成的地表反射率數據,選取近紅外波段(841—875 nm)和短波紅外波段(1628—1652 nm),空間分辨率為500 m)和MOD17A3(選取NPP,分辨率為1 km,時間分辨率為1年)。利用MRT批處理工具對遙感數據進行了拼接、格式轉換和重投影,并在ArcGIS中完成掩膜,使用最大合成法(Maximum Value Composite, MVC)將MOD13Q1和MOD09A1產品分別合成月數據。

1.2.2 氣象數據

2015年的氣象數據來源于中國氣象局氣象數據中心(http://data.cma.cn/site/index.html),剔除數據缺失和有異常值的站點,最終選取33個氣象站點的月數據集,包括：日照百分率、月平均氣溫和月降水量數據。

1.2.3 土地利用數據

100 m分辨率的土地利用數據來源于中國科學院資源環境科學數據中心(http://www.resdc.cn)。本文參照馮益明等[17]針對甘肅省提出的土地利用/覆蓋類型的劃分,將研究區的土地利用/覆蓋類型劃分為9種類型：耕地、林地、草地、建筑用地、水域、沼澤地、裸沙地、戈壁和其他用地,并將其分辨率重采樣為250 m。

1.2.4 其他數據

河西走廊行政區劃圖來源于國家基礎地理信息中心(http://ngcc.sbsm.gov.cn)。DEM數據獲取自地理空間數據云(http://www.gscloud.cn),在DEM基礎上提取各柵格單元的海拔、坡度和坡向信息。

2 研究方法

CASA模型估算植被NPP的基本思想是利用植被獲取太陽輻射,加上植被自身利用的情況,從而估算出植被凈生長狀況[18]。模型中所估算的NPP可以由植被吸收的光合有效輻射(APAR)和實際光能利用率(ε)兩個因子來表示[1,4],公式如下：

NPP(x,t)=APAR(x,t)×ε(x,t)

(1)

式中,x代表單個像元,t表示月份,APAR(x,t)則表示像元x在t月吸收的光合有效輻射(gC/m2),ε(x,t)表示單個像元x在t月的實際光能利用率(gC/MJ)。

2.1 APAR(x,t)的估算

APAR(x,t)是通過植被所能吸收的太陽有效輻射與植被對入射光合有效輻射的吸收比例來確定,公式如下：

APAR(x,t)=SOL(x,t)×FPAR(x,t)×0.5

(2)

式中,SOL(x,t)表示t月在單個像元x處的太陽總輻射(gC/m2),FPAR(x,t)表示植被層對入射光合有效輻射的吸收比例,常數0.5表示植被所能利用的太陽有效輻射(波長為0.4—0.7 μm)占太陽總輻射的比例。

由于研究區輻射站點較少,直接利用站點太陽輻射數據插值誤差較大,故本文SOL(x,t)的估算采用和清華等[8]將天文輻射作為起始值計算西部地區太陽總輻射的公式,

Q=QA(a+bs)

(3)

式中,Q為太陽總輻射；a、b的取值分別為0.185、0.595；s表示日照百分率；QA表示天文輻射,鑒于篇幅有限,這里不再贅述,詳細計算請參考文獻[8]。

植被對太陽有效輻射的吸收比例取決于植被類型和植被覆蓋狀況[19]。研究證明,由遙感數據得到的歸一化植被指數(NDVI)能很好地反映植被覆蓋狀況[4]。本文利用NASA-MOD15算法中設計的NDVI-FPAR查找表(Mynenietal.,1999)計算每月截取的光合有效輻射比例(FPAR),選取的是MODIS-NDVI數據,其數據質量和分辨率都高于原模型中運用的NOAA/AVHRR NDVI數據[20-21],計算式如下：

(4)

2.2 光能利用率的估算

光能利用率是在一定時期單位面積上生產的干物質中所包含的化學潛能與同一時間投射到該面積上的光合有效輻射能之比。環境因子如氣溫、土壤水分狀況以及大氣水汽壓差等會通過影響植被的光合能力而調節植被的NPP,計算公式如下：

ε(x,t)=Tε1(x,t)×Tε2(x,t)×Wε(x,t)×εmax

(5)

式中,Tε1(x,t)和Tε2(x,t)表示低溫和高溫對光能利用率的脅迫作用；Wε(x,t)為水分脅迫系數,反映水分條件對光能利用率的影響；εmax是理想條件下的最大光能利用率(gC/MJ),取值因植被類型而不同,取值參考文獻[17]。

2.2.1 溫度脅迫系數估算

Tε1(x,t)的估算：其反映在低溫和高溫時植被內在的生化作用對光合的限制而降低凈初級生產力。

Tε1(x,t)=0.8+0.02×Topt(x)-0.0005×[Topt(x)]2

(6)

式中,Topt(x)表示某一區域一年內NDVI值達到最高時的當月平均氣溫(℃),即植被生長的最適溫度。

Tε2(x,t)的估算：表示環境溫度從最適溫度Topt(x)向高溫或低溫變化時植被光能利用率逐漸變小的趨勢,這是因為低溫和高溫時高的呼吸消耗必將會降低光能利用率,生長在偏離最適溫度的條件下,其光能利用率也一定會降低,計算公式如下：

Tε2(x,t)=1.184/{1+exp[0.2×(Topt(x)-10-T(x,t))]}×1/{1+exp[0.3×(-Topt(x)-10+T(x,t))]}

(7)

式中,T(x,t)表示月均溫,Topt(x)同上。當月均溫T(x,t)比最適溫度Topt(x)高10℃或低13℃時,該月的Tε2(x,t)值等于月平均溫度T(x,t)為最適溫度Topt(x)時Tε2(x,t)值的一半。

郭婧等提出AMMRR插值方法(“多元回歸分析+殘差法”)能充分體現復雜多變的地形特點[22]。因此,本文通過經度,緯度和海拔3個變量,利用AMMRR方法對氣溫進行插值,該插值方法不僅考慮了海拔對氣象要素的影響,而且還利用站點數據與插值數據之間的差值進行了修正,提高了插值精度[23-24]。

2.2.2水分脅迫系數Wε(x,t)的計算

水分脅迫系數Wε(x,t)反映了植被所能利用的有效水分條件對光能利用率的影響,隨著環境中有效水分的增加,Wε(x,t)逐漸增大,CASA模型中它的取值范圍為0.5(在極端干旱條件下)到1(非常濕潤條件下)。

由于Wε(x,t)與植被水分含量直接相關,地表水分指數(LSWI),是描述植被葉片水分含量的一項指標[25-27]。短波紅外(SWIR)對植被水分含量敏感,近紅外和短波紅外波段已經被用來獲取對水分敏感的植被指數(LSWI),計算公式如下：

(8)

式中,LSWI計算公式見(9),LSWImax表示單個像元生長期內LSWI的最大值,使用MVC方法估算。

(9)

這里,ρnir和ρswir分別表示MOD09A1數據的近紅外和短波紅外波段。

Xiao等[25]在VPM模型中首次運用該方法,目前已有許多學者將LSWI引入CASA模型來計算Wε(x,t)[28-29]。值的注意的是,上述公式計算出來的Wε(x,t)的取值范圍為0(非常濕潤)到1(極端干旱),與CASA模型中Wε(x,t)的取值正好相反,為了適應模型算法,Bao等[15]提出了修改算法,并證明該方法適用于干旱半干旱區。因此,本文采用地表反射率產品MOD09A1計算Wε(x,t),這不僅增強了數據的準確性,也提高了數據的分辨率,同時遙感數據也包含了地形信息[30],公式如下：

Wε(x,t)= (1-(1+LSWI)/(1+LSWImax)) +0.5

(10)

3 結果與討論

3.1 精度評價

對NPP的驗證可采用實測值驗證法和相對驗證法兩種方法,考慮到研究區范圍較大且植被類型復雜,NPP實測值短時間內較難獲取,故本文采用與NPP產品比較和與前人的研究成果比較的方法。

MOD17A3產品已經被驗證適用于區域和全球尺度的NPP研究[31-32],也有很多學者將產品與估算值進行比較來判斷估算值的可靠性[33-36]。由于2015年的MOD17A3產品未更新,只能使用2014年產品,本文采用同樣的方法估算出2014年NPP,兩者隨機抽樣比較,結果見圖2。可以看出,校正后的NPP與MOD17A3產品的相關性較強,R2為0.803(P<0.01),明顯高于校正前的R2(0.681)。

圖2 校正前后NPP與MOD17A3產品的比較Fig.2 Comparison between NPP of correction before and after and NPP products of MOD17A3

本文也與前人的相關研究進行了對比,河西走廊NPP的估算鮮有研究,因而采用干旱區其他區域的NPP均值進行比較。高志強等[37]利用CASA模型估算的2000年西北地區年均NPP為102.52 gC m-2a-1;裘駿一[38]利用CASA模型估算中衛市沙坡頭自然保護區植被的年均NPP為60—135 gC m-2a-1；焦偉等[39]利用CASA模型估算了西北干旱區2000—2014年的植被NPP,年均值為191.63 gC m-2a-1；本文基于校正的河西走廊年均NPP為151.51 gC m-2a-1,與前人的研究結果相近。各植被類型的比較見表1,分析發現,草地、耕地和荒漠的模擬結果與前人研究的結果相對較好。林地、水域和沼澤與朱文泉和焦偉的結果相差較大。

造成差異的原因較多：一是數據源的差異；二是研究區空間尺度和空間分辨率的差異,可能存在混合像元；三是研究時段不同；四還存在算法的差異等。總的來說,校正后的NPP與MOD17A3產品相關性更高,與前人的計算結果更加接近,由此判斷,校正后的NPP估算結果更加準確。

3.2 河西走廊NPP的時空分布特征

2015年河西走廊NPP空間分布見圖3,NPP自東南向西北遞減,介于0—1699.31 gC m-2a-1之間,年均值為151.51 gC m-2a-1,總量為34.29 TgC/a,高值區主要分布在祁連山和綠洲,低值區分布在北部荒漠,主要植被類型NPP見表1。

表1 河西走廊年均NPP與文獻值的比較/(gC m-2 a-1)

圖3 2015年NPP空間分布圖/(gC m-2 a-1)Fig.3 Spatial distribution of NPP in 2015

從流域來看,石羊河流域NPP總量為13.30 TgC/a,均值為307.19 gC m-2a-1,黑河流域總量為13.53 TgC/a、均值為243.59 gC m-2a-1,疏勒河流域總量為7.45 TgC/a,均值為58.44 gC m-2a-1。石羊河流域中上游以耕地和草地為主,NPP值高,下游荒漠和戈壁面積大,NPP值低。黑河流域上游祁連山的東南部及南部NPP較高,植被類型主要為森林和草地,中游綠洲、荒漠和戈壁相間分布。疏勒河流域上游祁連山區的植被類型主要為高山草甸,NPP較高,下游荒漠、戈壁廣布,NPP較低,年累積NPP高值區主要集中在綠洲,范圍較小。

3.3 校正結果的比較

3.3.1 校正前后平均氣溫比較

圖4 2015年7月平均氣溫校正對比Fig.4 Comparison of mean temperature correction in July, 2015

3.3.2 校正前后水分脅迫系數比較

圖5是2015年1月和7月河西走廊校正前后水分脅迫系數空間分布。河西走廊除祁連山區外屬于干旱區,植被基本都處于水分脅迫狀態。校正前1月和7月的水分脅迫系數在0.50—0.58之間,差異很小,全區域均處于高度水分脅迫狀態,與祁連山區植物生長的水分脅迫程度并不相符,王海軍等[44]的研究結果表明祁連山區的降水具有明顯的季節性,冬季降水均在13 mm以下,而在夏季降水量最高可達247 mm,因而在夏季祁連山區的水分脅迫程度明顯降低。在干旱半干旱地區植被生產力的限制因子隨海拔會發生轉變,低海拔地區植被主要受到低降水導致的干旱脅迫,同時較高的溫度導致蒸散增加,進一步加劇了植被的水分脅迫[45]。校正后1月和7月的水分脅迫系數分布在0.5—1之間,更能反映河西走廊制約植被生長的水分脅迫情況,區別出地形變化對水分脅迫的影響差異,同時還能較好的識別綠洲等人類活動區。1月份雖然降水較少,但植被停止生長并不需要太多水分,7月份降水較多,但植被在生長季需水也最多,因此,水分脅迫在生長季(7月)比1月份更高(除祁連山外)是合理的。值的一提的是,水分脅迫系數校正前后的計算原理差異很大,前者是實際蒸散量與潛在蒸散量的比值,只與降水、太陽凈輻射和溫度有關,與植被無關,后者是根據地表反射率反演水分指數而得,與植被有關,即同樣的水分狀況,針對不同的植被和不同的生長階段是不一樣的,因此,校正更能反映植被的水分脅迫程度,這種情況在中下游綠洲區表現尤為明顯。位于圖西北部有部分高值,這是由于MOD09A1影像中的云量所致。

圖5 2015年1月和7月水分脅迫系數校正對比Fig.5 Comparison of moisture stress coefficients in January and July, 2015

3.3.3 校正前后NPP比較

圖6 2015年NPP校正對比Fig.6 NPP contrast used correction in 2015

2015年河西走廊校正前后的NPP空間分布基本一致,也存在較為明顯的差異,見圖6。可以看出,地形平坦的下游荒漠區吻合較好,地形起伏較大的南部高山區,綠洲區存在明顯差異,校正前前者低估,后兩者高估。整體來看,差值介于-526.891—296.5 gC m-2a-1之間,在山區和綠洲相差較大,平原地區差異較小。

校正對南部高山區的溫度脅迫和水分脅迫均有影響,由圖4、5可以看出,校正前氣溫高估和水分低估,氣溫高估導致蒸發量增加,增強該地區的水分脅迫,這個效應與校正前水分低估疊加,加劇水分脅迫導致該區域NPP低估。綠洲區地勢平坦,校正對氣溫的影響較小,基于地表反射率的水分脅迫校正能準確反應植被的需水狀態,主要植被類型是耕地,在生長季水分需求大,校正后水分脅迫程度更高,因此校正前NPP存在高估。此現象可進一步拓展到人類活動區域,即改變了原始植被類型的區域,校正對水分脅迫系數的影響非常大,進而對NPP估算產生影響。除綠洲區外的其他荒漠區地形起伏較小,植被分布稀少,吻合較好。

Sun等[46]利用地形校正的BEPS(Boreal Ecosystem Productivity Simulator)模型模擬了武陵山區NPP,得出校正后NPP低于校正前的地區主要分布在山區,與本研究的高山地區校正前NPP低估相反,其原因是武夷山區植被的主要脅迫因子是溫度,校正后溫度降低使得溫度脅迫加強,校正后NPP減少。平地區校正前后NPP偏差較小,在10 gC m-2a-1以內,高山地區NPP偏差可能會達到500 gC m-2a-1,與本文的結論一致。Chen等[47]對秦嶺西南部山區NPP研究發現地形因子對NPP的影響很大,尤其是在山區,也與本文的結果一致。

3.4 討論

3.4.1 校正對月NPP的影響

當前，中央電視臺大型綜藝節目《朗讀者》的熱播，對中小學階段的學生產生了重要的影響。很多央視的“名嘴”加入朗讀行列，對小學語文課本的朗讀，可以說是給當前的小學生很好的示范和帶頭作用。在小學語文教學中開展朗讀訓練，朗讀的過程同時也是閱讀的過程，久而久之有助于促進學生語文思維能力的提升，讓學生獲得思想情感的熏陶，讓學生樂在其中，享受其中。

圖7 校正對月NPP的影響Fig.7 Impact of correction on monthly NPP

從2015年河西走廊校正前后各月NPP變化曲線(圖7)可以看出,各月NPP差異并不相同,校正后的NPP除7月份外均高于校正前,5月份NPP校正的影響相差最大(低估14.03 gC/m2)。春季(3、4和5月)是植被開始生長期,對水熱的變化最敏感,地形對NPP的影響較大,校正前NPP低估了21.27 gC/m2。夏季(6、7和8月)生物量大量累積,也是水熱需求量最大的時期,校正前低估了6.75 gC/m2。秋季(9、10和11月)氣溫下降,植被基本停止生長,此時植被對水熱的需要是維持性的,校正對NPP影響較小；冬季(12、1和2月)為非生長季,校正對NPP影響可忽略不計。

為了展示校正對月NPP的影響細節(圖8),本文在肅南裕固族自治縣西部截取高低估和吻合區域相間分布的一部分地區(圖1),海拔介于1737—5332 m,地形起伏較大,降水為140—216 mm。校正在非生長季的1、2、3、11和12月對NPP的影響較小,大多處于-10—10 gC/m2,而在生長季的4、5、6、7、8、9和10月,校正對NPP的影響較大,介于-290—161 gC/m2之間。校正前6、7、8月份出現較多高估區,該區域海拔較高,氣溫較低降水比較充沛,主要脅迫因子為溫度,校正后溫度降低,因此校正前是高估的。

圖8 校正前后月NPP變化Fig.8 Changes of monthly NPP correction before and after

3.4.2 校正對地形因子與NPP關系的影響

地形能夠導致局部的水熱再分配,尤其是干旱地區,水熱組合差異影響植物生長甚至造成植被種類不同[48]。本文由DEM提取出海拔、坡度、坡向信息,利用等差分級法(海拔級差為30 m[49]、坡度級差為3°[50]、坡向分為9級[51]),校正對NPP與地形關系的影響見圖9。

圖9 校正對地形因子與NPP關系的影響Fig.9 Effect of correction on the relationship between terrain factors and NPP

校正對NPP在不同海拔范圍內的影響各異(圖9)。在2200 m以下荒漠所占面積大,地勢平坦,校正對NPP的影響較小,與上文平原區變化不大結論相符。2200—3500 m范圍內,部分綠洲區處于該高程帶,校正前NPP高于校正后,其原因是基于LSWI的校正更能反映其水分脅迫狀況。3500—5000 m海拔較高,氣溫下降,植被蒸騰作用減弱,校正后水分更能滿足植被生長需求,校正前低估了NPP。海拔5000 m以上的面積非常小,而且都是冰川雪地,NPP基本為0,校正的影響可忽略不計。

校正對坡度與NPP關系的影響也較大,主要體現在土壤含水量上,坡度與土壤含水量呈負相關,即隨坡度的增大而降低[52]。小于15°屬于緩坡,坡度對土壤水分持有量的影響不顯著[53],所以校正對該區間影響不大(圖9)；大于15°區間,校正后NPP減少,其原因是在重力作用下土壤含水量流失的坡度效應,并且,坡度越大兩者的差異越大,因此,隨坡度增加校正對NPP的影響也越大。

坡向導致地表接受的太陽輻射不同,使得不同坡向水熱有顯著差異[54],從而影響植物生長。各個坡向的NPP變化見圖9,東坡NPP最高,南坡NPP最小。校正后北坡,東北坡、西北坡和西坡高于校正前,平地、東坡、東南坡、南坡、西南坡低于校正前。北坡差異最大,校正前低估了7.12 gC m-2a-1,其次是東南坡,校正前高估了3.38 gC m-2a-1,其他坡向的差值都低于2 gC m-2a-1。坡向校正前不能夠識別水熱在坡向上的變化,校正后陽坡太陽輻射強烈氣溫高蒸發強烈,水分減少,導致NPP減少,陰坡氣溫低,蒸發較弱,NPP增加,與Sun等[46]在武陵山區的研究結果一致。

4 結論

(1)溫度脅迫系數基于地理因子回歸校正,水分脅迫系數基于遙感數據校正能有效提高CASA模型的估算精度。河西走廊NPP在原CASA模型中存在高估,校正后NPP總量為34.29 TgC/a,原CASA模型高估了0.23 TgC/a。

(2)校正對高海拔和地形起伏較大區域,以及人類活動地區的NPP估算影響較大。在河西走廊,前者主要是溫度脅迫的影響造成,后者主要受水分脅迫的影響。

(3)綠洲區校正結果表明,在人類活動區域,基于遙感數據的水分脅迫計算比原模型中基于蒸散量的計算更加可靠,后者存在高估。

(4)校正對月份,海拔、坡度以及坡向與NPP關系的影響也較大。校正對生長季的影響大于非生長季；海拔上,2200—3500 m范圍內,校正前高估NPP,在3500—5000 m,校正前低估了NPP；坡度上,校正前高估了NPP,坡度越大兩者的差異越大；坡向上,校正前高估了陽坡NPP,低估了陰坡NPP。

本研究的不足之處,一是地表反射率產品存在有云現象,導致部分區域估算不準確;二是太陽輻射是通過克呂格插值得到,精度有限;三是缺乏NPP實地采樣驗證,今后都有待深入研究。