朱躲萍,葉 輝,王軍邦,趙烜嵐,左 嬋,蘆光新,張法偉,李英年
1 九江學院旅游與地理學院,九江 332005 2 中國科學院地理科學與資源研究所 生態系統網絡觀測與模擬重點實驗室 生態系統大數據與模擬中心,北京 100101 3 長江大學園藝園林學院,荊州 434000 4 青海大學農牧學院,西寧 810016 5 中國科學院西北高原生物研究所,西寧 810008
土地利用變化和溫室氣體排放被認為是全球氣候變暖的主要因素,并分別從其生物地球物理過程及生物地球化學過程開展研究[1—3]。表征下墊面反射短波輻射通量與入射短波輻射通量比值的地表反照率,是影響地球能量平衡的關鍵變量[4—5],是生物地球物理過程中植被變化對氣候的反饋機制之一[1, 6—8]。
在植被變化對氣候的反饋機制研究中,對地表反照率的影響存在不同的解釋[3]。模型模擬研究表明,若全球尺度毀林,其凈氣候效應是通過改變地表反照率與蒸散而抑制由于毀林產生的溫室氣體排放導致的氣候變暖[9],這也得到了觀測數據的支持[10]。過去30年全球植被綠度增加減緩了全球陸地氣候變暖增溫的12%,其中地表反照率降低貢獻了6%[2]。而在青藏高原的研究發現,該地區植被綠度增加,盡管地表反照率降低使地表吸收的輻射增多,但同時也使地表蒸散增強,最終導致了減緩氣候變暖的降溫效應[11]。然而草地退化與恢復作為土地利用變化的一種漸變方式,所導致的地表反照率的變化及其影響機制,特別是其輻射效應,仍然缺乏研究[2, 12—15]。
作為長江源、黃河源和瀾滄江源發源地的青海三江源地區,其分布最廣的草地生態系統不僅維持著多種生態功能,也為當地牧民畜牧業發揮重要的生活和生產功能;其敏感而又脆弱的生態系統,是我國重要的生態安全屏障之一[16],在全球氣候變化研究領域中占據著重要地位。本文基于2001—2018年中分辨率成像光譜儀(MODIS)地表反照率等數據[17],研究三江源區土地覆蓋與利用變化對地表反照率的時空變化的影響,及其地表輻射溫度效應,以深入理解區域土地利用與覆被變化的生物地球物理過程與機制,為適應全球氣候變化、生態系統管理及科學決策提供支撐。
三江源地處青藏高原腹地(圖1),地理范圍約為31°39′—36°12′N,89°45′—102°23′E,海拔為3335—6564 m,總面積為36.3萬km2[16],是長江、黃河和瀾滄江的發源地,素有“中華水塔”之美稱。三江源具有冷熱兩季交替、干濕兩季分明、年溫差小、日溫差大、日照時間長、輻射強烈、無明顯四季區分的氣候特征,為典型的高原大陸性氣候[18]。

圖1 三江源區地理位置及主要土地覆蓋類型
根據2018年30 m空間分辨率土地利用數據一級分類(圖1)[19],三江源的植被以草地為主,占全區面積的71.2%,其中高、中、低覆蓋草地分別占全區的5.7%、24.9%、40.6%;未利用地占16.7%,森林僅占全區面積的4.3%(表1)。長江源園區、黃河源園區以及瀾滄江源園區的草地面積占比分別為85.3%、89.4%和50.8%,未利用地面積分別占40.4%、3.4%和3.6%,森林為0.01%、0.30%和0.90%。其中,長江源園區、黃河源園區以及瀾滄江源園區草地高、中、低覆蓋草地面積占比分別為1.1%、42.6%、41.6%,0.3%、23.6%、65.5%和2.5%、12.7%、35.6%。

表1 2018年三江源以及三個園區土地利用與土地覆被類型及面積/(×104 km2)
2.1.1地表反照率
采用2001—2018年每16天0.05°(約5 km)空間分辨率MODIS的地表反照率產品(MCD43C3),是基于分別搭載在Terra和Aqua兩個衛星上的傳感器數據合成的數據產品[20]。該產品包括太陽時中午的黑空短波反照率、白空反照率、可見光以及近紅外波段反照率,精度評估基于3個Ameriflux通量塔及53個FLUXNET塔點觀測的反照率數據,均方根誤差小于0.013[21]。本研究選擇了黑空短波反照率數據進行了分析研究,數據下載自美國地質調查局官網。
2.1.2歸一化植被指數
歸一化植被指數(Normalized Difference Vegetation Index,NDVI),是利用2001—2018年250 m空間分辨率的MODIS波段反射率數據產品(MOD09Q1)計算而得。在計算得到NDVI后,通過Timesat軟件進行了SG濾波[22—23],以消除云等對NDVI數據的影響。
2.1.3氣象數據
以中國氣象科學數據共享網(http://date.cma.cn/)提供的研究區及周邊氣象臺站觀測的降水量、最高和最低氣溫日值數據,采用氣象數據空間插值軟件ANUSPLIN,經空間插值得到三江源8天時間分辨率250 m空間分辨率的柵格數據,用于分析。ANUSPLIN軟件是基于普通薄膜和局部薄膜樣條函數多變量數據進行內插的工具,能同時進行多表面空間插值[24],插值時以90 m空間分辨率SRTM數字高程模型(DEM),經空間重采樣為250 m的DEM數據作為協變量。
2.1.4土地利用和覆蓋(LUCC)數據
本研究所用2018年三江源地區LUCC數據,來源于中國科學院資源環境科學數據中心(http://www.resdc.cn)。該數據空間分辨率為30 m,是根據土地資源的利用屬性和自然屬性,以及地貌位置等信息通過人工目視解譯的方式進行分類;其二級類型包括水田、旱田、有林地、疏林地、灌叢、高覆蓋草地、中覆蓋草地、低覆蓋草地、水域、建設用地、濕地、荒漠及未利用地等[19]。
2.1.5地表反照率地面觀測數據
以通量觀測塔測定的上行和下行短波輻射,計算得到的反照率作為地表反照率地面觀測數據,對地表反照率遙感數據進行了精度評價。通量塔測定數據來自中國通量觀測網絡(ChinaFLUX)海北高寒草地長期定位觀測站(37.6°N,101.3°E)。通量塔所提供的數據為2015—2017年每30分鐘的數值,在精度評價時,以每天11:00—15:00時范圍內數值計算均值,作為每日地表反照率,然后計算月平均值,與對應時段MCD43C3月平均數據進行對比分析。
2.2.1地表反照率數據質量控制
對于16天和0.05°時空分辨率MCD43C3地表反照率數據產品,經投影轉換、范圍裁剪、數據空間重采樣以及缺值插補,處理為8天250 m時空分辨率數據,以與植被指數及氣象數據時空匹配。針對本研究所用的黑空短波輻射反照率數據,根據其質量標識數據(QA)進行了如下質量控制:選擇QA小于等于3的值作為有效值,對QA大于等于4的像元值,作為缺值進行插補。當QA等于4,地表反照率值小于三倍有效值的標準差時,保留原值,否則以QA小于3的像元有效值的平均值進行插補[25]。
2.2.2太陽短波輻射
到達地面的太陽總輻射的計算,一般需要考慮地形、云量、大氣透明度、太陽光線入射角以及太陽高度角等因素[26],這里根據文獻采用下式計算:
(1)
式中,Ra為地面晴空太陽總輻射(MJ m-2d-1),N為可照時數,n為實照時數,n/N為日照百分率,as和bs的值由地區確定[27]。
2.2.3歸一化差植被指數(NDVI)
歸一化差植被指數(NDVI),MODIS波段反射率數據計算得到[25],計算公式如下:
(2)
式中,NR為近紅外波段,RR為可見光波段,即近紅外波段和可見光波段之差除以近紅外與可見光波段和。
2.2.4隨機森林回歸算法
隨機森林算法最早是在2001年由Breiman提出,主要用于解決分類和回歸的問題,是通過集成學習思想將多棵決策樹集成的一種機器學習算法[28]。隨機森林算法能提高預測能力,具有抗擬合能力,模型的泛化能力強,參數優化簡單[29—30]。本研究基于地表反照率逐像元的年際變化影響因子分析,選擇NDVI、最高氣溫、最低氣溫以及降水4個變量,構建了地表反照率年際變化的隨機森林回歸模型,對變量重要性進行評估,來判斷像元尺度對地表反照率年際變化影響因子。本文通過Python中RandomForestRegressor庫實現了隨機森林回歸模型,變量重要性X按式(3)計算:
X=∑(E1-E2)/Ntree
(3)
式中,E1為每一棵樹對應的袋外數據誤差,E2為對袋外數據的變量加入噪聲干擾,重新排序和再次計算的袋外數據誤差,Ntree隨機森林測試樣本。
2.2.5一元線性回歸分析
本研究采用最小二乘法擬合2001—2018年地表反照率的變化趨勢,并進行像元尺度和分土地覆蓋類型的年際變化分析,計算給出統計顯著性水平P值。
2.2.6地表反照率變化的輻射溫度效應計算
在陸面能量收支對氣候變化的貢獻計算中,不考慮各種反饋過程,僅考慮地表反照率對輻射平衡的影響,輻射平衡狀態時[31],計算公式如式(4):
(4)
式中,RS是太陽短波輻射(MJ m-2d-1),a是地表反照率,ε是比輻射率(灰體系數,取1.0),δ是斯蒂芬-玻耳茲曼常數(為5.67×10-8W m-2K-4)。根據式(4),草地覆蓋度變化引起的輻射溫度效應,采用下式計算[27]:
(5)
式中是草地覆蓋度變化引起的輻射溫度效應,Q為研究時段內平均年太陽短波輻射,a1、a2分別為前后兩個研究時段不同覆蓋度草地的平均地表反照率。本文中分別以2001—2004年期間和2005—2018年期間為草地覆蓋度變化前后時段,分別組合計算高、中、低覆蓋度草地變化的輻射溫度效應。
2.2.7植被覆蓋度
根據李苗苗等像元二分模型計算植被覆蓋度[32]:
(6)
式中,NDVIS生長季平均歸一化差植被指數NDVIveg和NDVIsoil分別為植被完全覆蓋時和無植被覆蓋時NDVIS值,以土地利用單元內頻率為95%和5%的NDVIS值分別確定。
如圖2所示,衛星遙感的地表反照率數據產品MCD43C3,在以海北站通量塔測定的下行和上行短波輻射計算的地表反射率進行精度評價時,二者之間存在顯著的季節變化一致性,復相關系數R2為0.79,均方根誤差為0.022,顯著性水平P小于0.001,表明該數據產品在三江源區域具有較高的數據精度和可靠性。

圖2 基于2015—2017年海北通量塔觀測數據對MCD43C3產品的精度評價
如圖3所示,2001—2018年三江源區平均反照率空間分布具有較大的空間變化,總體上呈自東南向西北逐漸增大的變化格局。區域平均值為0.163±0.027,主要集中在0.12—0.18之間,占總面積的76.4%;部分地表反照率大于0.18的地區,占總面積的21%,主要分布在三江源東北部和中西部;最高值達到0.56,分布于研究區西部地區戈壁、沙地以及裸土地等未利用地;研究區東南部地表反照率在0.12以下。長江源、黃河源和瀾滄江源三個園區的地表反照率差異較大,以長江源園區最高(0.177±0.036),瀾滄江源園區(0.156±0.002)略高于黃河源園區(0.153±0.037)。
在圖3中,三江源區各土地覆蓋類型的地表反照率,未利用地的最大,為(0.176±0.006),林地的最小(0.136±0.02);不同覆蓋草地間,低覆蓋草地的地表反照率高于高覆蓋草地,具體地高、中、低覆蓋草地地表反照率分別為(0.154±0.003)、(0.155±0.003)、(0.163±0.004)。值得注意的是長江源園區的未利用地的地表反照率,為(0.182±0.008),相比黃河源園區(0.162±0.006)和瀾滄江源園區(0.156±0.004)更高,因為未利用地主要集中在長江源園區,而黃河源園區和瀾滄江源園區主要分布著低覆蓋草地(表1),在瀾滄江源園區的森林(0.9%)的地表反照率較長江源園區(0.01%)和黃河源園區(0.3%)更低,可能是因為林地的冠層會影響地表反照率,冠層可以吸收可見光波段,其面積的增加會使其對反照率的吸收增加,導致地表反照率下降。

圖3 2001—2018三江源生長季年均地表反照率空間分布和三江源及三個園區各LUCC類型生長季平均地表反照率柱狀圖
如圖4所示,2001—2018年三江源區生長季地表反照率變化整體呈每年不顯著下降趨勢。全區平均下降速率為(0.152±0.763)%/10a(P=0.47),趨勢速率的空間變化范圍為-0.010—0.007之間。呈上升和下降趨勢的地區分別占全區總面積的41.42%和58.58%,其中,顯著下降區域占總面積的8.4%,主要分布在三江源區的東北部。呈顯著增加的面積僅占全區1.9%,零散分布在三江源東部。長江源園區、黃河源園區以及瀾滄江源園區平均地表反照率的年際變化為統計不顯著趨勢下降(P>0.05),每10年分別下降(0.078±0.900)%、(0.215±0.740)%、(0.152±0.450)%的速率下降。
如表2中,三江源不同土地覆蓋與利用類型生長季地表反照率均以不同趨勢下降。其中高、中、低覆蓋草地分別以-0.0003%/10a、-0.0013%、0.0001%/10a速率下降,而林地和未利用地的下降速率分別為每10年0.0003%和0.0014%,但下降趨勢均統計不顯著(P>0.05)。

表2 三江源植被類型年際變化趨勢、相關系數及其顯著性
應用隨機森林以2001—2018年三江源區域氣候因子(最高氣溫、最低氣溫、降水)和植被因子(NDVI)為自變量,以地表反照率為應變量,自變量的權重如圖5所示,NDVI、降水、最高氣溫、最低氣溫為主導因子的區域,分別占全區面積的40.7%、16.4%、29.0%和13.8%,表明植被是影響地表反照率年際變化的主要因素,且各自變量的權重值分別為(0.31±0.149)、(0.22±0.106)、(0.21±0.098)、(0.26±0.131),如圖6所示,表明植被因子對地表反照率的年際變化具有最高的貢獻。NDVI的權重值,在三個園區間相比時,在黃河源園區貢獻最高(0.319±0.148)、其次是瀾滄江源園區(0.316±0.151),在長江源園區貢獻最低(0.283±0.131)。

圖5 2001—2018年三江源生長季地表反照率年際變化的貢獻因子權重:NDVI、降水、最低氣溫和最高氣溫

圖6 2001—2018年三江源區及長江源園區、黃河源園區、瀾滄江源園區(NDVI、降水、最低氣溫和最高氣溫)各因子對地表反照率年際變化的貢獻頻數分布
以2001—2004年期間為基準,采用2001—2018年生長季平均太陽短波輻射,以消除太陽輻射變化的影響,定量分析2005—2018年期間生長季地表反照率變化的輻射溫度效應,結果如圖7,期間植被覆蓋度的變化全區平均值為(0.01±0.08),即增加了1%[33],增加區域占全區面積的57.1%;地表反照率變化整體上呈輻射增溫效應,增溫(0.11±0.42)℃,輻射溫度變化值主要集中在-0.3—0.6℃之間,占全區面積的84.8%;空間上表現為中西部以及東北部為輻射增溫效應,東南和西北部為輻射降溫效應,分別占全區面積54.86%和45.14%。
空間疊加分析表明,植被覆蓋度增加區域既存在輻射增溫效應,也存在輻射降溫效應,分別占全區32.6%和24.5%。輻射增溫效應主要分布在東北部和黃河源園區,植被覆蓋度增加,地表反照率降低,地表吸收太陽凈輻射增多,進而對氣候產生了變暖的正反饋[34—35];植被覆蓋度增加伴隨著的輻射降溫效應主要分布在東南部,主要土地覆被類型為森林,說明這些地方森林綠度增加導致了地表反照率上升,從而導致凈輻射吸收減少表現出輻射降溫效應,已有研究也發現森林變綠會導致地表反照率增加[25]。
三江源42.9%的區域植被覆蓋度下降,也是存在輻射增溫和輻射降溫效應,分別占全區面積的21.7%和21.2%,主要分布于三江源中西部和西部地區。植被覆蓋度下降導致的地表反照率下降表現出的輻射增溫效應,可能是生長季土壤冰層融化、土壤水分增多,而土壤水分增多會引起地表反照率下降[36—37],從而凈輻射吸收增多表現為輻射增溫效應,西北部植被覆蓋度下降導致地表反照率上升,而地面吸收凈輻射減少表現為輻射降溫效應。在黃河源園區植被覆蓋度的增加導致了輻射增溫效應,值為(0.24±0.37)℃,而長江源園區和瀾滄江源園區植被變化引起的輻射溫度效應不顯著。

圖7 2005—2018年相比于2001—2004年基準期三江源生長季植被覆蓋度變化及地表反照率變化的輻射溫度效應
三江源生長季2015—2018年與2001—2004年高覆蓋草地漸變到低覆蓋草地NDVI值之差以及低覆蓋草地過渡到高覆蓋草地NDVI值之差與相應時段地表反照率變化之差引起的輻射溫度效應關系均呈正相關,如圖8所示,上升趨勢分別為0.060(R=0.34,P<0.001)和0.058(R=0.30,P<0.001),其中低覆蓋草地漸變到高覆蓋草地引起的草地覆蓋度增加(草地覆蓋度變化大于0)導致的輻射增溫效應(大于0℃),占選取隨機點總數的57.3%,高覆蓋草地過渡到低覆蓋草地引起的草地覆蓋度降低(草地覆蓋度變化小于0)導致輻射降溫效應(小于0℃),占選取隨機點總數的56.3%,從而得出生長季草地覆蓋度增加主要導致的地表反照率減少的輻射溫增效應,而植被覆蓋度變低主要導致的地表反照率增加的輻射降溫效應,這與地表反照率上升和下降引起的地面吸收的凈輻射下降和上升導致的輻射溫度效應相一致。

圖8 2005—2018年與2001—2004年三江源生長季高低覆蓋草地對應的NDVI值之差與輻射溫度效應散點回歸圖
土地覆被變化導致地表輻射特性參數改變進而影響地表輻射收支,其中作為關鍵物理參數的地表反照率的較小改變,將影響地氣系統的能量收支,進而影響到區域氣候[38]。三江源區域地表反照率在空間上呈現從東南到西北上升的分布特征,空間上差異性明顯,這與趙之重等人得出的三江源地表反照率空間分布規律具有較好的一致性[39];其中,地表反照率較低的東南部,分布著森林[40],而森林較其他覆蓋類型具有較低的地表反照率[41—44]。三江源區生長季地表反照率變化整體呈現不顯著下降趨勢,且相對較為穩定,這與2001—2010年期間整體不顯著增加趨勢的研究結論基本一致[39];對整個青藏高原的研究表明,約占77.5%的區域地表反照率變化趨勢統計顯著,僅有20.1%(2.4%)的區域呈顯著降低(增加)趨勢[39]。三江源區域地表反照率年際變化的主導因子是植被因子,其中NDVI對地表反照率權重值最大占全區比接近全區面積的一半,且NDVI對地表反照率影響最大的區域與地表反照率高于0.15的區域(東北部和中西部)相重合,這些區域主要是沙地等低植被覆蓋度的未利用地,即與三江源區域地表反照率與NDVI負相關的結論一致[39],而降水、最高氣溫和最低氣溫對地表反照率的影響區域分布相對較零散。
除植被和氣候因子外,短波反照率的變化還受云覆蓋和積雪覆蓋及其變化的影響[44]。本研究對反照率數據質量(QA)做了控制處理,把填充值在25%以上的(QA=4、5)像元用臨近像元數值平均值代替,進而達到去云處理作用[25];對于地表積雪對地表反照率的影響[45],本研究選擇了6—8月作為生長季,以減弱積雪對地表反照率的影響。已往研究指出三江源部分區域存在多年凍土[46],而對此還需要進一步開展觀測和機理模擬研究和分析。
本研究進一步分析了三江源區域植被覆蓋度和地表反照率變化的地表輻射溫度效應,表明植被恢復可能會導致輻射增溫效應,從而對區域氣候變暖產生正反饋。但值得注意的是,本研究是基于衛星遙感觀測數據的統計分析,在未考慮短波輻射本身年際變化前提下,僅僅計算了由地表反照率變化導致的輻射溫度效應,所得到的結論值得參考和進一步思考。當然,今后研究需進一步考慮由地表凈輻射增加后,地表輻射在潛熱和顯熱間的分配,這就需要基于陸表過程模式的模擬和多源數據整合的方法從機制過程方面開展研究,進而更準確地確定其輻射溫度效應。
三江源生長季地表反照率在空間上呈現從東南到西北上升的分布特征,空間上差異性明顯,全區平均地表反照率為0.163±0.027,其中76.4%區域的值集中在0.12—0.18。2001—2018年生長季三江源地表反照率年際變化趨勢范圍為-0.010—0.007,顯著下降區域占總面積的8.4%,顯著上升區域占總面積1.9%;總體上三江源地表反照率以每10年(0.152±0.763)%速率不顯著下降,而植被是這種變化的主導因素。三江源生長季植被覆蓋度增加整體上導致輻射增溫效應,升溫(0.11±0.42)℃。由植被覆蓋度的增大導致的輻射增溫占總面積的32.6%,導致的輻射降溫占總面積的24.5%,同時覆蓋度下降導致的輻射增溫效應占總面積的21.7%,導致的輻射降溫效應占總面積的21.2%。且草地覆蓋度增加,引起地表反照率下降,導致輻射增溫效應;而草地覆蓋度降低,引起地表反照率上升,導致輻射降溫效應。今后需基于陸表過程模式模擬和多源數據整合方法從機制和過程方面開展研究,并考慮積雪覆蓋、永久凍土等對地表輻射收支的影響;本研究基于衛星遙感數據的統計分析,有助于對草地退化與恢復過程中的生物地球物理效應的理解和進一步思考。