侯曉臣,孫 偉1,,李建貴,李全勝
(1.中國農業科學院農業信息研究所,北京 100081; 2.新疆農業大學計算機與信息工程學院,新疆 烏魯木齊 830052;3.新疆農業大學林業研究所,新疆 烏魯木齊 830052)
塔里木河位于西北干旱區,是我國最大的內陸河,被譽為“生命河”和“母親河”[1],在新疆發展戰略中占據著重要地位。塔里木河干流自身不產流,徑流補給全部來自源流,為純耗散型河流,上游區的水碳資源狀況對于中下游地區乃至整個干流區生態系統的可持續發展有重要影響。近幾年隨著經濟的快速發展,塔里木河干流上游區的資源開發強度不斷加大,致使塔里木河干流流域的生態環境正逐步惡化,表現為干流徑流量逐年遞減,越往下游遞減趨勢越明顯[2];植被衰退,土地荒漠化趨勢加劇,下游綠色走廊瀕臨消失[3]。水資源短缺和土地荒漠化成為阻礙區域發展的兩大主要因素,如何有效協調二者的可持續性發展已成為研究區當前亟需解決的問題之一。水碳資源的科學統籌與優化管理是解決這一問題的有效途徑,其基本前提是對水碳耦合機制的深入研究與應用[4]。張永強和朱治林等[5-6]發現水分利用效率在一天內呈現由高到低的趨勢,且最高值并不是出現在正午時分。郭維華和蘇培璽等[7-8]發現植物的水碳通量變化具有較高的一致性。冠層尺度上,生態系統的水碳通量與太陽輻射具有較強的一致性[8]。Baldocchi等[9]認為陸地生態系統的年總生態系統生產力和蒸散發之間保持著很好的一致性。Beer等[10]發現,世界范圍內40%以上林地的總初級生產力與降水量密切相關。Zhao等[11]發現,陸地凈初級生產力的變化趨勢與區域干旱程度的變化特征具有很強的相關性。目前,在塔里木河流域干流區關于水、碳資源的動態變化研究很多[12-13],關于水碳耦合關系的研究仍然相對較少,沒有形成比較科學成熟的體系。
現階段,渦度相關技術可提供大量的、長期的、連續的水碳通量相關數據,水碳耦合模型已逐步成為研究熱點。現有的水碳耦合模型可以劃分為2種:(1)基于光合-氣孔-蒸騰機理構建的水碳耦合模型,主要有CEVSA模型[14]、BEPS模型[15]、IBIS模型[16]等;(2)集成現有水文模型和生態模型構建的集成模型,主要有RHESSys模型[17]、DLEM模型[18]、WaSSI-C模型[19]等。其中,WaSSI-C模型(Water Supply Stress Index_Carbon model)是集成了水分供需計算模塊(Water Supply Stress Index,WASSI)和水碳循環模塊的月尺度生態系統水碳耦合集成模型,主要融合了Hamon的潛在蒸散(PET)模型[20]、Mccabe G J等[21]的融雪模型以及薩克拉門托模型,且如地形起伏、土壤濕度、植被狀況和水分劃分等關鍵水碳要素都得到了充分的考慮,已被應用于我國和美國的部分地區,并且取得了很好的應用效果[22-23]。
本文的研究對象是典型干旱區塔里木河干流上游區,由于該區域降水稀少,水資源主要來源于高山冰川消融。因此,冰川消融模擬對于WaSSI-C模型在研究區的應用具有關鍵性的作用,而原WaSSI-C模型的融雪模塊忽視了高山冰川消融對區域的水分補充作用,嚴重影響了研究區的模擬精度,因此在原模型融雪模塊的基礎上補充冰川消融計算過程,提高了模型在研究區的適用性。本研究探討WaSSI-C模型在塔里木河干流上游區的適用性,以期為研究區的水碳資源統籌管理提供有效的分析工具。
塔里木河干流流經塔克拉瑪干沙漠北部邊緣,總長1 321 km,分上、中、下游3個區段,其中肖夾克至英巴扎為上游段,河長495 km;英巴扎至恰拉為中游段,河長398 km;恰拉至臺特瑪湖為下游段,河長428 km,其水量補充主要來源于阿克蘇河、和田河、葉爾羌河和開都-孔雀河,形成了塔里木河“四源一干”的基本格局。基于已獲取的DEM數據,利用ArcGIS的流域分析工具提取出塔里木河干流上游流域范圍(80°45′~85°10′E,40°25′~41°30′N),并將其作為本研究的研究區域。該區總面積約為1.5萬km2,地勢平坦,海拔908~1 045 m,屬于典型的暖溫帶內陸干旱荒漠性氣候,干旱少雨,晝夜溫差大,年降水量在20~50 mm之間,年均潛在蒸散量在1 800~2 900 mm之間。近幾年,隨著人類活動強度的不斷加大,大片的林、灌、草地被開發為農田,耕地面積迅速增長,引發了農業用水與生態用水的矛盾[24],進而使得水碳資源矛盾日益突出。因此,探討WaSSI-C模型在該地區的適用性,對掌握研究區的水、碳資源動態變化及內在關系,進而有效平衡二者的協調發展具有重要的現實意義。
WaSSI-C模型運行需要基礎數據包括輸入數據和驗證數據。輸入數據包括數字高程模型(DEM)、氣象數據、葉面積指數(LAI)數據和土地利用類型數據等;驗證數據包含徑流數據、蒸散(ET)數據和總初級生產力(GEP)數據等。數據獲取和處理過程如下:
(1)DEM數據來源于美國國家航空航天局(NASA)的ASTER-GDEM V2產品(http://gdem.ersdac.jspacesystems.or.jp/),空間分辨率為30 m,應用ArcGIS對原始數據進行拼接、投影變換、提取等處理后,可得模型所需的DEM數據。
(2)基于研究區的Landsat 5遙感影像,利用ENVI軟件進行監督分類和人工目視解譯等預處理,得到2005年塔里木河干流上游區的土地利用類型數據。
(3)流域氣象數據包含降水和氣溫數據,基于國家氣象局2000-2015年的月降水和氣溫站點數據,使用ArcGIS的協同克里金插值,以經緯度和高程作為協變量,可得連續分布的30 m×30 m的降水和氣溫數據。
(4)基于英巴扎水文站2000-2015年的月徑流量數據,利用ArcGIS計算研究區面積,可將徑流量數據轉化為徑流深數據。
(5)葉面積指數(LAI)數據、蒸散(ET)數據和總初級生產力(GEP)數據使用MODIS數據,經拼接、投影變換等預處理操作后可得研究區1 km×1 km的相應數據。
本研究采用相關系數(R2)和效率系數(NS)對模型在塔里木河干流上游區的適用性進行分析和評價。其中,相關系數(R2)的計算公式為:
(1)
式中,O為觀測值,S為模擬值,n為觀測數據的個數,R2可在EXCEL中利用線性回歸法自動求得,0≤R2≤1,其值越小表示數據吻合程度越差,當R2=1時證明模擬效果達到最優,數據完全吻合。效率系數(NS)的計算公式如下:
(2)

WaSSI-C模型是以月為時間尺度,以水文響應單元為空間尺度的,旨在模擬水碳資源動態變化過程的生態系統水碳耦合模型,其核心模塊是水量平衡模塊和生態系統生產力模塊[19]。水量平衡模塊是生態系統生產力模塊的計算基礎,主要用來模擬流域生態系統的降水、蒸發、滲透、積雪和徑流過程,包含了蒸散量計算過程、冰雪累積與融化過程、土壤水分及滲透計算過程、徑流計算過程等;生態系統生產力模塊又稱為碳循環模塊,它的計算反映了水碳的耦合過程,該模塊利用經驗回歸模型進行關系構建,旨在模擬流域內碳的收益和損失過程,見圖1。
2.2.1 蒸散過程
PEThamon=0.1651×n×k×ρw
(3)
式中,PEThamon為Hamon潛在蒸散(mm);n為月時長(d);k為日晝長(12 h);ρw為月均飽和蒸汽密度(g·m-3)。
PAET=a×PEThamon+b×LAI+c×P×PEThamon
(4)
式中,PAET為不考慮土壤水分狀況條件下的植被蒸散潛力(mm);LAI為月均葉面積指數;P為降水量(mm);a、b、c為經驗參數。
2.2.2 融雪過程
Pin=Prain+SnowM
(5)
式中,Pin為總來水量(mm);Prain為降雨量(mm);SnowM為地表積雪融化量(mm)。
2.2.3 水循環過程

圖1 WaSSI-C模型的理論框架和邏輯結構Fig.1 The theoretical framework and logical structure of the WaSSI-C model
薩克拉門托土壤濕度計算模型和McCabe的融雪模型是WaSSI-C模型水循環計算過程的思想來源,基于土壤水分含量在垂直方向上的差異性,將其劃分為上下兩層,進而模擬土壤水分的循環過程。水循環過程旨在模擬研究區產流過程和實際蒸散過程,在其計算過程中充分考慮了自然地形和植被狀況對土壤水分循環的影響,提高了模型模擬結果的準確性。地表徑流產生后直接匯入河網;經線型水庫的調節后,地下徑流和壤中流依次匯入河網;地表地下產流和壤中流相加并扣除一定時段內的蒸散發后,得到總入流量,經河網調節后,得到最終徑流量。實際蒸散水分主要來源于土壤上層張力水和自由水及下層張力水,具體計算過程如下:
當PAET (6) 當UZTWC ET=PAET (7) 當PAET>UZTWC+UZFWC時, ET=UZTWC+UZFWC+ (8) 式中,ET為實際蒸散(mm);PAET為不考慮土壤水分狀況條件下的植被蒸散潛力(mm);UZTWC為土壤上層張力水含量(mm);UZTWM為土壤上層張力水最大容量(mm);UZFWC為土壤上層自由水含量(mm);LZTWC為土壤下層張力水含量(mm);LZTWM為土壤下層張力水最大容量(mm)。 生態系統生產力模型旨在模擬生態系統的碳循環動態過程,其構建思路主要是基于全球通量網絡中的GEP和ET數據,利用線性回歸模型建立回歸方程。因此回歸方程斜率代表了基于GEP的水分利用效率,基于此可評價和分析研究區的水碳特征。 GEP=k×ET (9) 式中,GEP為總生態系統生產力(g·m-2);ET為實際蒸散(mm);k為經驗參數。 研究區地處中國西部內陸干旱區,遠離海洋,降水稀少,其水資源主要來源于高山冰雪融化。然而在WaSSI-C模型的融雪過程中,區域總來水量主要由降雨量和地表積雪融化量組成,這顯然與研究區實際情況不符。因此,在融雪過程中加入冰川融化量計算模塊[25-26],以提高模型在研究區的適用性。 SnowI=f×(T+9)t (10) Pin=Prain+SnowM+SnowI (11) 式中,SnowI為冰川融化量(mm);T為氣溫(℃);f,t為經驗參數;Pin為總來水量(mm);Prain為降雨量(mm);SnowM為地表積雪融化量(mm)。 基于已獲取的實測徑流數據和MODIS的蒸散(ET)數據和總生態系統生產力(GEP)數據,將2000-2009年作為WaSSI-C模型的率定期,2010-2015年作為驗證期,利用決定系數(R2)和效率系數(NS)分別對2個時期的總徑流模擬效果、蒸散模擬效果和總生態系統生產力模擬效果進行定量評價,進而分析WaSSI-C模型在研究區的適用性。 模型在塔里木河干流上游區的總徑流模擬結果在2個時期均呈現出良好的模擬效果。在率定期,總徑流模擬值與實測值的R2為0.72,NS為0.71,總體上吻合效果較為理想,僅在徑流峰值處存在差異略大的情況(圖2)。在驗證期,總徑流模擬值與實測值的R2為0.68,NS為0.67,與率定期相似;在徑流量較小時,模擬徑流與實測徑流基本完全吻合,當徑流量較大時,模擬效果較差(圖3)。綜合模型在率定期和驗證期的徑流模擬評價結果可以發現,模型總體上可以很好地反映研究區的徑流動態變化,但除個別年份(2009年)外,不同時期的總徑流模型均存在著徑流峰值低估的問題。造成這樣問題的原因可能是:(1)研究區產流主要以超滲產流為主,而模型的水循環模擬中兼有超滲產流和蓄滿產流,這在一定程度上使得模擬結果在徑流峰值處被低估[27];(2)補充的冰川融化量計算模塊受溫度限制過大,忽略了其他因素對冰川融化的影響。另外,趙銳鋒等[2]研究了1957-2005年塔里木河干流年徑流量的變化趨勢,計算出英巴扎水文站的年均徑流量為28.47×108m3,將其轉化為徑流深為185 mm,并指出塔里木河干流徑流量表現出遞減的變化趨勢。楊鵬等[28]指出塔里木河徑流量在2006年后持續減少,2009年是塔里木河在過去50 a里水量最少的一年,而2010年為豐水年。本研究模擬所得2000-2015年的年均徑流深為150 mm,與趙銳峰等的研究略有差異,這可能是由于研究時間尺度不同所造成的,在年際變化上與楊鵬等的研究具有一致性,進而從側面驗證了徑流模擬結果的可靠性。 模型在率定期(2000-2009年)和驗證期(2010-2015年)的蒸散模擬在不同時期都取得了良好的模擬效果。在率定期,模型模擬蒸散值與MODIS蒸散值對比的R2和NS均為0.60,不同月份平均蒸散模擬值與MODIS值整體保持較好的一致性。但也存在差異較大的情況,在5月和11月模擬值略高于MODIS值,而在7月模擬值要明顯低于MODIS值(圖4)。在驗證期,模型模擬蒸散值與MODIS蒸散值對比的R2為0.64,NS為0.56,與率定期相似,不同月份模擬蒸散與MODIS蒸散整體上具有一定的一致性,但在部分月份差異略大,即4月、5月和11月模擬值高于MODIS值,7月模擬值要明顯低于MODIS值(圖5)。綜合模型在2個時期的蒸散模擬評價結果可以發現,模型總體上可以很好地反映研究區的蒸散變化特征,但也存在部分地區和月份模擬蒸散與MODIS蒸散差異略大的問題,其原因可能是:(1)MODIS蒸散在無植被區域為無效值,本研究在處理時將其設置為0,進而導致其計算蒸散有所下降;(2)MODIS的蒸散是基于飽和水汽壓差計算的,忽略了土壤水分對實際蒸散的影響,可能會造成對實際蒸散值的高估。目前,有關塔里木河干流上游區的實際蒸散估算研究多為塔里木河干流區乃至新疆范圍內的實際蒸散研究。因此,利用前人研究成果對塔里木河干流上游區實際蒸散難以進行準確驗證,只能基于區域自然地理特征對其進行比較驗證。代超[29]基于區域水量平衡計算出塔里木河干流區年平均陸地蒸發深為162.3 mm;阿迪來·烏甫等[30]在研究新疆地表蒸散時空分布與變化趨勢時指出塔里木盆地邊緣地區年實際蒸散在48~248 mm之間。本研究模擬所得塔里木河干流上游區年均實際蒸散為217 mm;各水文響應單元的年均實際蒸散值在110~360 mm之間;考慮到塔里木河干流上中下游的差異性及研究區的地理位置,可以認為模型模擬所得實際蒸散值是較為可靠的。 圖2 率定期月總徑流的WaSSI-C模擬值與觀測值的對比Fig.2 Comparison between the simulated values by WaSSI-C and the observed values of themonthly total runoff during calibration period 圖3 驗證期月總徑流的WaSSI-C模擬值與觀測值的對比Fig.3 Comparison between the simulated values by WaSSI-C and the observed values of the monthly total runoff during validation period 模型的GEP模擬效果能夠較好地再現研究區的GEP動態變化過程。在率定期,模型模擬GEP值與MODIS的GEP值對比的R2為0.68,NS為0.66,不同月份GEP模擬值與MODIS值整體保持較好的一致性,僅在個別月份(4-7月)有較大的差異(圖6)。在驗證期,模型模擬蒸散值與MODIS蒸散值對比的R2為0.61,NS為0.69,不同月份模擬GEP與MODIS的GEP計算結果整體比較一致,僅在部分月份(4-5月、7月)差異略大(圖7)。綜合模型在2個時期的GEP模擬評價結果可以發現,在率定期和驗證期的兩個評價指標均大于0.6,處于可接受范圍,但在部分月份GEP模擬值與MODIS的GEP值存在較大差異,這可能是由于:(1)MODIS的GEP計算時利用了輻射數,而本文使用了溫度, 二者間盡管有較強的相關性, 但并不完全等價[31];(2)MODIS基于全球氣象數據來計算GEP,而本文利用基于觀測點的插值氣象數據,二者間的差異增加了模型的不確定性。陳福軍等[32]利用陸地生態系統碳循環模型( CASA 模型) 對中國陸地生態系統近30 a的NPP時空變化特征進行了研究,指出中國西部荒漠生態系統NPP 年均值在100 g·m-2以下。盧玲等[33]指出新疆南部年均NPP累積量在0~100 g·m-2之間。本研究在2000-2015年的年均GEP模擬值為88 g·m-2,與陳福軍和盧玲等[32-33]的研究基本一致。 本研究在充分研究WaSSI-C模型理論框架和運行機理的基礎上,針對塔里木河干流上游區自然地理和資源環境特征,通過增加冰川融化計算過程對模型進行了適用性改進,進而探討了修正后的WaSSI-C模型在塔里木河干流上游區的適用性,并得到了如下結論: 1)修正后的WaSSI-C模型能夠很好地應用于研究區。統計評價指標R2和NS均顯示出WaSSI-C模型對研究區的水、碳模擬效果良好。雖然模擬結果和驗證數據間有著一定的不一致性,但這可能是由于輸入數據和驗證數據本身可能具有的誤差以及模型的不確定性對模擬效果評價的影響導致的,而且,模型計算結果在徑流、蒸散和生產力方面的數據范圍與相關研究大體一致[2,28-30,32-33]。因此可以認為WaSSI-C模型經適用性改進后能夠很好地應用于研究區。 圖4 率定期不同月份平均蒸散的WaSSI-C模擬值與MODIS值的對比Fig.4 Comparison of WaSSI-C and MODIS values of average monthly evapotranspiration during calibration period 圖6 率定期不同月份平均總生態系統生產力(GEP)的WaSSI-C模擬值與MODIS值的對比Fig.6 Comparison of WaSSI-C and MODIS values of average monthly gross ecosystem productivity during calibration period 圖5 驗證期不同月份平均蒸散的WaSSI-C模擬值與MODIS值的對比Fig.5 Comparison of WaSSI-C and MODIS values of average monthly evapotranspiration during validation period 圖7 驗證期不同月份平均總生態系統生產力(GEP)的WaSSI-C模擬值與MODIS值的對比Fig.7 Comparison of WaSSI-C and MODIS values of average monthly gross ecosystem productivity during validation period 2)根據WaSSI-C模型在研究區的3個模擬結果變量驗證可以發現:模型對于徑流的模擬效果最優,對于ET和GEP 的模擬效果略差。具體表現為率定期和驗證期的ET和GEP對比驗證的R2和NS均低于徑流對比驗證的R2和NS值;在模擬結果與驗證數據的對比曲線上,模擬徑流與實測徑流呈現出更好的一致性。這可能是由于缺乏實測的ET和GEP數據用于模型的率定和驗證,進而增加了模型在ET和GEP模擬中的不確定性。因此,在未來除繼續對模型進行適用性改進外,尋求更高質量的ET和GEP數據將成為提升WaSSI-C模型在研究區適用性的關鍵因素之一。2.3 生態系統生產力模塊
2.4 WaSSI-C模型的適應性改進
3 結果與分析
3.1 總徑流模擬結果評價與分析
3.2 蒸散模擬結果評價與分析


3.3 碳生產力模擬結果評價與分析
4 結論與討論



