侯曉臣，孫偉，李建貴，李全勝

(1.新疆農(nóng)業(yè)大學計算機與信息工程學院,新疆 烏魯木齊 830052；2.新疆農(nóng)業(yè)大學林業(yè)研究所,新疆 烏魯木齊 830052)

焉耆盆地作為典型的干旱綠洲區(qū)，是巴音郭楞蒙古自治州乃至新疆重要的農(nóng)牧業(yè)基地之一.近年來，隨著綠洲城市的不斷發(fā)展，以及全球氣候變暖趨勢日益加劇[1-3]，焉耆盆地內(nèi)植被恢復、水資源供應及氣候變化之間的多重矛盾逐漸顯現(xiàn)，水碳資源的科學統(tǒng)籌與優(yōu)化管理是解決這一矛盾的有效途徑，其基本前提是對水碳耦合機制的深入研究與應用[4].隨著生態(tài)系統(tǒng)水碳關系研究的逐步深入，二者間的耦合關系已經(jīng)在不同的時空尺度上得到了證明[5-12].現(xiàn)階段，關于焉耆盆地水碳資源的研究，多集中于水文過程及植被分布變化[13-14]，有關水碳資源耦合的研究相對較少.

當前，渦度相關技術可提供大量的、長期的、連續(xù)的水碳通量相關數(shù)據(jù)[19-20]，水碳耦合機制及其定量化評估模型逐步成為研究熱點[5-16].自20世紀80年代開始，F(xiàn)arquhar等[17]提出的光合模型和Ball等[18]提出的光合與氣孔導度關系模型是水碳耦合模型研究的雛形.現(xiàn)有的水碳耦合模型可以劃分為兩種：(1)基于光合-氣孔-蒸騰機理構建的水碳耦合模型，主要有CEVSA模型[21]、BEPS模型[22]、IBIS模型[23]等；(2)集成現(xiàn)有水文模型和生態(tài)模型構建的集成模型[24]，主要有RHESSys模型[25]、DLEM模型[26]、WaSSI-C模型[27]等.其中，WaSSI-C模型(Water supply stress index_carbon model)融合了Hamon的潛在蒸散(PET)模型[28]、Mccabe的融雪模型[29]以及薩克拉門托模型[30]，已在我國和美國的部分地區(qū)取得了很好的應用效果[31-32].本研究將干旱綠洲區(qū)的焉耆盆地作為研究對象，針對研究區(qū)的自然地理和資源環(huán)境特征，在原WaSSI-C模型的基礎上增加冰川融化計算模塊，以提高WaSSI-C模型在研究區(qū)的適用性，進而模擬并分析研究區(qū)水碳資源特征，以期為焉耆盆地水碳資源綜合管理提供可靠的科學數(shù)據(jù)基礎和決策工具支撐.

1 材料與方法

1.1 研究區(qū)概況

焉耆盆地是新疆天山山脈中的一個斷陷盆地(E 85°50′～87°50′，N 41°40′～42°20′)，面積為7 231 km2；盆地地形西北高東南低，海拔1 047～1 409 m，盆地東南部的博斯騰湖為最低處；屬于暖溫帶大陸性干旱氣候，冬季嚴寒，夏季炎熱，氣溫年較差大，2000～2015年，年均溫為9.2 ℃，1月平均氣溫-9 ℃，7月平均氣溫24.4 ℃；整體降水量少，年降水量為50～80mm；光熱資源豐富，具有典型的干旱區(qū)綠洲氣候特征.盆地內(nèi)主要有4條河流流經(jīng)，即開羅河、黃水溝、清水河和孔雀河.在行政區(qū)劃上，焉耆盆地位于新疆巴音郭楞蒙古自治州內(nèi)，包含了焉耆、和靜、和碩和博湖縣的平原部分.

圖1 研究區(qū)示意圖Figure 1 Diagram of the study area

1.2 數(shù)據(jù)獲取與處理

WaSSI-C模型運行的輸入數(shù)據(jù)包含地形數(shù)據(jù)、氣象數(shù)據(jù)、葉面積指數(shù)(LAI)數(shù)據(jù)、土地覆被數(shù)據(jù)等，驗證數(shù)據(jù)包括徑流數(shù)據(jù)、蒸散(ET)數(shù)據(jù)和總生態(tài)系統(tǒng)生產(chǎn)力(GEP)數(shù)據(jù)等，具體情況如下.

(1)數(shù)字高程模型DEM為美國國家航空航天局(NASA)的ASTER-GDEM V2產(chǎn)品，利用ArcGIS軟件水文分析工具，經(jīng)填洼、流向流量分析、流域分水嶺提取等處理過程，將研究區(qū)劃分為42個水文響應單元.

(2)土地覆被數(shù)據(jù)來源于美國陸地衛(wèi)星5號(Landsat 5)，通過人工解譯，并利用ENVI軟件進行監(jiān)督分類，依據(jù)國土資源部組織修訂的國家標準《土地利用現(xiàn)狀分類》(GB/T 21010-2017)，將焉耆盆地劃分為建設用地、耕地、林地、草地、未利用地和水域6種土地類型，生成2005年焉耆盆地土地覆被數(shù)據(jù).

(3)氣象數(shù)據(jù)來源于國家氣象局，包括月尺度的降水和氣溫數(shù)據(jù)，利用ArcGIS軟件的地統(tǒng)計模塊，引入經(jīng)度、緯度和高程作為協(xié)變量，采用協(xié)同克里金插值工具得到1 km×1 km的降水和溫度柵格數(shù)據(jù)集.經(jīng)處理后可得2000～2015年不同月份的平均降水量和氣溫.

(4)水文數(shù)據(jù)為月徑流監(jiān)測數(shù)據(jù)，來源于焉耆盆地內(nèi)焉耆水文站(圖1).由于模型模擬結果為徑流深數(shù)據(jù)，因此利用ArcGIS軟件的曲面面積計算工具計算出流域曲面面積并將流量數(shù)據(jù)轉化為徑流深數(shù)據(jù)，從而用于模型的驗證.

此外，WaSSI-C模型用到的其他數(shù)據(jù)包括：MODIS葉面積指數(shù)數(shù)據(jù)(LAI)(MODIS15)、MODIS蒸散數(shù)據(jù)(ET)(MODIS16)和MODIS總初級生產(chǎn)力數(shù)據(jù)(GPP)(MODIS17)，來源于美國國家航空航天管理局網(wǎng)站，地址為http://modis.gsfc.nasa.gov.在生態(tài)系統(tǒng)尺度上，GPP與GEP近似相等，二者均表征生態(tài)系統(tǒng)的碳生產(chǎn)能力，因此本研究將MODIS的GPP數(shù)據(jù)作為總生態(tài)系統(tǒng)生產(chǎn)力(GEP)的驗證數(shù)據(jù).

1.3 模型評價方法

本研究采用的評價指標是決定系數(shù)(R2)和效率系數(shù)(NS)[31，33].決定系數(shù)指示了回歸直線對樣本數(shù)據(jù)的擬合程度，其值在0～1，決定系數(shù)越高說明模擬效果越好.效率系數(shù)常被用于評價水文模型模擬效果的優(yōu)劣，其值最大為1，最小為負無窮，且越接近于1，模型計算結果越好.

(1)

(2)

式(1)、(2)中，Qp為模擬值，Qm為觀測值，Qavg為觀測值的平均值，n為觀測數(shù)據(jù)個數(shù).

2 WaSSI-C模型與改進

2.1 WaSSI-C模型理論框架和運行機理

WaSSI-C模型是以月為時間尺度，以水文響應單元為空間尺度，旨在模擬水碳資源動態(tài)變化過程的生態(tài)系統(tǒng)水碳耦合模型，核心模塊是水量平衡模塊和生態(tài)系統(tǒng)生產(chǎn)力模塊[34]，其理論框架和邏輯結構見圖2.

2.2 水量平衡模塊

2.2.1 蒸散過程

PEThamon=0.1651×n×k×ρw

(3)

PAET=a×PEThamon+b×LAI+c×P×PEThamon

(4)

式(3)中，PEThamon為Hamon潛在蒸散；n為月時長(d)；k為日晝長(12 h)；ρw為月均飽和蒸汽密度(g/m3)；式(4)中，PAET為不考慮土壤水分狀況條件下的植被蒸散潛力；LAI為月均葉面積指數(shù)；P為降水量(mm)；a、b、c為經(jīng)驗參數(shù).

2.2.2 融雪過程

Pin=Prain+SnowM

(5)

式(5)中，Pin為總來水量(mm)；Prain為降雨量(mm)；SnowM為地表積雪融化量(mm).

2.2.3 水循環(huán)過程 薩克拉門托土壤濕度計算模型和Mccabe的融雪模型是WaSSI-C模型水循環(huán)計算過程的思想來源，基于土壤水分含量在垂直方向上的差異性，將其劃分為上下兩層，進而模擬土壤水分的循環(huán)過程.水循環(huán)過程主要包括蒸散發(fā)過程、地表地下產(chǎn)流過程和壤中流過程.其中，土壤上下層的張力水是蒸散發(fā)過程的主要水分來源，地表徑流產(chǎn)生后直接匯入河網(wǎng)，經(jīng)線性水庫的調(diào)節(jié)后，地下徑流和壤中流依次匯入河網(wǎng)，地表地下產(chǎn)流和壤中流相加并扣除一定時段內(nèi)的蒸散發(fā)后，得到總入流量，經(jīng)河網(wǎng)調(diào)節(jié)后，得到最終徑流量.在該過程中，自然地形和植被狀況對土壤水分循環(huán)的影響得到了充分的考慮，提高了模型模擬結果的準確性.

圖2 WaSSI-C模型的理論框架和邏輯結構Figure 2 The theoretical framework and logical structure of the WaSSI-C model

2.3 生態(tài)系統(tǒng)生產(chǎn)力模塊

生態(tài)系統(tǒng)生產(chǎn)力模型旨在模擬生態(tài)系統(tǒng)的碳循環(huán)動態(tài)過程，其構建思路主要是基于全球通量網(wǎng)絡中的GEP和ET數(shù)據(jù)，利用線性回歸模型建立回歸方程.因此回歸方程斜率代表了基于GEP的水分利用效率，基于此可評價和分析研究區(qū)的碳循環(huán)特征.

GEP=k×ET

(6)

式(6)中，GEP為總生態(tài)系統(tǒng)生產(chǎn)力；ET為實際蒸散；k為經(jīng)驗參數(shù).

2.4 WaSSI-C模型的適應性改進

焉耆盆地處中國西部內(nèi)陸干旱區(qū)，遠離海洋，降水稀少，其水資源主要來源于高山冰雪融化.然而在WaSSI-C模型的融雪過程中，區(qū)域總來水量主要由降水量和地表積雪融化量組成，這顯然與研究區(qū)實際情況不符.因此，在融雪過程中加入冰川融化量計算模塊[35-36]，以提高模型在焉耆盆地的適用性.

SmowI=f×(T×9)t

(7)

Pin=Prain+SmowM+SnowI

(8)

式(7)中，SnowI為冰川融化量；f，t為經(jīng)驗參數(shù).

3 結果與分析

3.1 WaSSI-C模型模擬結果評價與適用性分析

本研究將2000～2009年作為WaSSI-C模型的率定期，2010～2015年作為驗證期，利用決定系數(shù)(R2)和效率系數(shù)(NS)，結合水文站的實測徑流數(shù)據(jù)和MODIS的ET和GEP產(chǎn)品數(shù)據(jù)[37-38]，分別對兩個時期的總徑流模擬效果、蒸散模擬效果和總生態(tài)系統(tǒng)生產(chǎn)力模擬效果進行定量評價.

3.1.1 總徑流模擬結果評價與適用性分析 率定期總徑流模擬值與實測值的R2和NS分別為0.80和0.77，模擬值與觀測值基本吻合，但也應注意到在低徑流期和高徑流期模擬值與觀測值仍有一定差距(圖3)；驗證期總徑流模擬值與實測值的R2和NS分別為0.77和0.69，擬合線斜率為0.58，對比曲線的吻合程度相對率定期略低，但基本吻合，與率定期相似，當徑流量較小或較大時，模擬值與觀測值存在較大差異(圖4).造成模型在低徑流期和高徑流期模擬效果有所下降的原因可能是模型融雪過程中冰川融化計算模塊受溫度限制過大，而忽略了其他水文要素對徑流的影響.

圖3 率定期月總徑流(mm·month-1)的WaSSI-C模擬值與觀測值的對比曲線圖Figure 3 Graph of comparison between the simulated values by WaSSI-C and the observed values of the monthly total runoff (mm·month-1)during calibration period

圖4 驗證期月總徑流(mm·month-1)的WaSSI-C模擬值與觀測值的對比曲線圖Figure 4 Graph of comparison between the simulated values by WaSSI-C and the observed values of the monthly total runoff (mm·month-1) during validation period

3.1.2 蒸散(ET)模擬結果評價與適用性分析 由于用于模型蒸散對比驗證的MODIS16數(shù)據(jù)缺乏對水域下墊面蒸散值的計算結果值，因此在本研究中僅對除博斯騰湖區(qū)域(水文響應單元41)外其他區(qū)域的蒸散模擬結果值進行對比驗證.率定期模型的模擬蒸散值與MODIS的蒸散值的R2和NS分別為0.82和0.79，在4月、7月和8月，MODIS蒸散值略高于模擬蒸散值，但差異值并不大，而在其他月份MODIS和模型計算的蒸散值非常接近(圖5)；驗證期模型的模擬蒸散值與MODIS的蒸散值的R2和NS分別為0.80和0.76，在6～8月，MODIS蒸散值略高于模型得到的蒸散值，但差異值較小，而在其他月份MODIS和模型計算的蒸散值基本完全吻合(圖6).造成模型在蒸散量較高時模擬值低于MODIS蒸散值的原因可能是估算方法的差異，MODIS基于飽和水汽壓差估算蒸散，忽略了土壤水分對實際蒸散的限制，而本研究考慮了土壤水分的限制.

圖5 率定期焉耆盆地不同月份平均蒸散的WaSSI-C模擬值與MODIS值的對比Figure 5 Comparison of WaSSI-C and MODIS values of average evapotranspiration in different months of Yanqi Basin during calibration period

圖6 率定期焉耆盆地不同月份平均蒸散的WaSSI-C模擬值與MODIS值的對比Figure 6 Comparison of WaSSI-C and MODIS values of average evapotranspiration in different months of Yanqi Basin during validation period

3.1.3 總生態(tài)系統(tǒng)生產(chǎn)力(GEP)模擬結果驗證 率定期GEP的模擬值與MODIS值的R2和NS分別為0.88和0.84，GEP模擬值與MODIS值總體上呈現(xiàn)出較好的一致性，僅在個別月份(2月、11月和6～8月)有較小的差異(圖7)；驗證期GEP的模擬值與MODIS值的R2和NS分別為0.87和0.82，僅在6～8月和10月有較小差異，在其他月份MODIS和模型計算的GEP值吻合程度很高(圖8).造成模型在個別月份的GEP模擬值與MODIS值之間存在較大差異的原因可能是：(1)MODIS的GEP計算時利用了輻射數(shù),而本研究使用了溫度，二者間盡管有較強的相關性，但并不完全等價[42].(2)MODIS基于全球氣象數(shù)據(jù)來計算GEP，而本研究利用基于觀測點的插值氣象數(shù)據(jù)，二者間的差異增加了模型的不確定性.

圖7 率定期焉耆盆地不同月份平均總生態(tài)系統(tǒng)生產(chǎn)力(GEP)的WaSSI-C模擬值與MODIS值的對比Figure 7 Comparison of WaSSI-C and MODIS values of average total ecosystem productivity in different months of Yanqi Basin during calibration period

圖8 率定期焉耆盆地不同月份平均總生態(tài)系統(tǒng)生產(chǎn)力(GEP)的WaSSI-C模擬值與MODIS值的對比Figure 8 Comparison of WaSSI-C and MODIS values of average total ecosystem productivity in different months of Yanqi Basin during validation period

3.2 焉耆盆地水碳資源空間分布特征

3.2.1 徑流空間分布特征 焉耆盆地各水文響應單元年均徑流深在36～220 mm，在空間分布上，總體呈現(xiàn)出“西高東低，北高南低”的特點(圖9).研究區(qū)水分補給主要來源于高山冰雪融水和降水.在盆地西部和北部(水文響應單元1～7、9、11～13、16～18、22、24、26～31、34～35、37)由于海拔較高，受冰雪融水補給最多，因此具有較高的產(chǎn)流量(97～220 mm/y)；而在盆地東部和南部(水文響應單元8、10、14～15、19、20～23、32～33、36、38～42)，其海拔較低，受冰雪融水補給較少，且由于上游用水截留的影響，導致在這一地區(qū)年均產(chǎn)流量較低(36～128 mm).另外，研究區(qū)水汽主要來自于西北方向[45]，造成盆地西北部的降水量較東南部更高，客觀上加劇形成了這一特征.

圖9 焉耆盆地年均徑流空間分布Figure 9 Spatial distribution of annual average evapotranspiration in Yanqi Basin

3.2.2 蒸散空間分布特征 焉耆盆地陸地表面各水文響應單元年均蒸散在15～416 mm，博斯騰湖區(qū)域年均蒸散在416～776 mm，在空間分布上呈現(xiàn)出“中間高，四周低”的特點(圖10).這可能是受到海拔和土地覆被類型的影響.在盆地西部和北部(水文響應單元1、2、4、5、7、9、11、12、15～18、20～22、24～28、31、34、35)由于海拔較高，導致其氣溫略低，且土地覆被類型多為未利用地，可用于蒸散的土壤水分較少，進而使得其年均蒸散量較低(15～282 mm)；在盆地南部和東南部(水文響應單元37～39、42)，主要是由于其產(chǎn)流量較低且未利用地比重較高，導致其年均蒸散量較低(15～282 mm)；而在盆地中部(水文響應單元3、6、8、10、13、14、19、23、30、32、33、36、40、41)，由于海拔較低，氣溫較高，且土地覆被類型以耕地和水域為主，涵養(yǎng)水源能力略差，進而導致其年均蒸散量較高(281～776 mm).

圖10 焉耆盆地年均蒸散空間分布Figure 10 Spatial distribution of annual average runoff in Yanqi Basin

3.2.3 總生態(tài)系統(tǒng)生產(chǎn)力空間分布特征 研究區(qū)各水文響應單元年均總生態(tài)系統(tǒng)生產(chǎn)力在20～161 g C/m2之間，空間分布上呈現(xiàn)出與蒸散相似的特征，即中間高，四周低(圖11).在盆地四周(水文響應單元1，2、4、5、7、9、11～13、15～18、20～22、24～31、33～39、42)，土地覆被類型多為未利用地，植被稀疏，使得其總生態(tài)系統(tǒng)生產(chǎn)力相對較低(20～118 g C/m2)；而在盆地中部(水文響應單元3、6、8、10、14、19、23、32、40)，土地覆被類型多為林地、耕地和草地，植被覆蓋度高，因此總生態(tài)系統(tǒng)生產(chǎn)力相對較高(117～161 g C/m2).由于模型在水域的總生態(tài)系統(tǒng)生產(chǎn)力值設定為0，因此本研究暫不研究博斯騰湖區(qū)域(水文響應單元41)的總生態(tài)系統(tǒng)生產(chǎn)力.

圖11 焉耆盆地年均總生態(tài)系統(tǒng)生產(chǎn)力空間分布Figure 11 The spatial distribution of the annual total ecosystem productivity in the Yanqi Basin

4 討論

本研究的徑流模擬結果年際變化在率定期呈現(xiàn)出先增后減的變化趨勢，且在2002年達到最高，這與麥麥提吐爾遜·艾則孜等[14]關于焉耆盆地徑流變化特征的研究結果保持一致，近16 a平均年徑流量132 mm，與畢雪麗[39]的研究結果略有差異，造成這一問題的原因可能是：(1)本研究計算了2000～2015年的平均徑流量，畢雪麗計算的是1961～2008年的平均徑流量，時間跨度的不同可能會引起計算結果的差異.(2)補充的冰川融化計算模塊受溫度限制過大，高徑流期模擬值與觀測值的差異大于低徑流期，進而使得模擬徑流值總體上有一定高估.WaSSI-C模型在各水文響應單元年均蒸散模擬值為15～416 mm，焉耆盆地年均蒸散模擬值為220 mm左右(博斯騰湖區(qū)域除外)，與徐永明等[40]、郭玉川等[41]在焉耆盆地的蒸散量估算結果差異較小.WaSSI-C模型在焉耆盆地的年均GEP模擬值在80 g C/m2左右，與陳福軍等[43]利用陸地生態(tài)系統(tǒng)碳循環(huán)模型( CASA模型)的估算結果基本一致，年內(nèi)GEP模擬值在0.11～20 g C/m2，略低于王輝等[44]的研究結果，這可能是由于：(1)模型基于水文響應單元計算GEP均值，進而降低了數(shù)據(jù)間的差異性；(2)模型的GEP計算依賴于實際蒸散，而實際蒸散的計算受氣溫影響較大，進而導致在氣溫較低和較高時期，GEP模擬值不夠準確；(3)王輝等的研究對象為整個開孔河流域，本研究僅針對焉耆盆地，研究范圍的差異可能導致研究結果的差異.盡管模擬結果和驗證數(shù)據(jù)間有著一定的不一致性，但這可能是由于輸入數(shù)據(jù)和驗證數(shù)據(jù)本身可能具有的誤差以及模型的不確定性對模擬效果評價的影響導致的，而且，模型計算結果在徑流、蒸散和生產(chǎn)力方面的數(shù)據(jù)范圍與相關研究大體一致.因此可以認為WaSSI-C模型經(jīng)適用性改進后能夠很好地應用于焉耆盆地.

另外，焉耆盆地水碳資源的空間分布的模擬結果與徐永明、郭玉川等[33-34]關于焉耆盆地的蒸散空間分布呈現(xiàn)出中部高于四周的特征極為相似，與王輝等[37]關于開孔河流域NPP高值在其研究區(qū)內(nèi)呈帶狀分布，并向四周擴散衰減的研究結果保持一致.由此可見，改進后WaSSI-C模型對于焉耆盆地的水碳資源空間分布的模擬結果具有較高的可信度.

5 結論

本研究基于WaSSI-C模型理論框架和運行機理的深入探討，針對焉耆盆地自然地理和資源環(huán)境特征，通過增加冰川融化計算過程對模型進行了適用性改進，并分析了焉耆盆地2000～2015年的水碳資源空間分布特征，得到了如下結論.

1) 改進后的WaSSI-C模型能夠很好地應用于焉耆盆地.統(tǒng)計評價指標R2和NS均顯示出WaSSI-C模型對研究區(qū)的水、碳模擬效果良好.

2) 焉耆盆地的水碳資源在空間分布上具有明顯的不均衡性.總體上，徑流空間分布呈現(xiàn)“西高東低、北高南低”的特點；蒸散和總生態(tài)系統(tǒng)生產(chǎn)力均呈現(xiàn)“中間高，四周低”的分布特點.具體表現(xiàn)為：盆地西部和北部具有高產(chǎn)流量、低蒸散量和低碳生產(chǎn)力；中部具有低徑流量、高蒸散量和高碳生產(chǎn)力；南部和東南部具有低徑流量、低蒸散量和低碳生產(chǎn)力.