侯晓臣，孙伟，李建贵，李全胜

(1.新疆农业大学计算机与信息工程学院,新疆 乌鲁木齐 830052；2.新疆农业大学林业研究所,新疆 乌鲁木齐 830052)

焉耆盆地作为典型的干旱绿洲区，是巴音郭楞蒙古自治州乃至新疆重要的农牧业基地之一.近年来，随着绿洲城市的不断发展，以及全球气候变暖趋势日益加剧[1-3]，焉耆盆地内植被恢复、水资源供应及气候变化之间的多重矛盾逐渐显现，水碳资源的科学统筹与优化管理是解决这一矛盾的有效途径，其基本前提是对水碳耦合机制的深入研究与应用[4].随着生态系统水碳关系研究的逐步深入，二者间的耦合关系已经在不同的时空尺度上得到了证明[5-12].现阶段，关于焉耆盆地水碳资源的研究，多集中于水文过程及植被分布变化[13-14]，有关水碳资源耦合的研究相对较少.

当前，涡度相关技术可提供大量的、长期的、连续的水碳通量相关数据[19-20]，水碳耦合机制及其定量化评估模型逐步成为研究热点[5-16].自20世纪80年代开始，Farquhar等[17]提出的光合模型和Ball等[18]提出的光合与气孔导度关系模型是水碳耦合模型研究的雏形.现有的水碳耦合模型可以划分为两种：(1)基于光合-气孔-蒸腾机理构建的水碳耦合模型，主要有CEVSA模型[21]、BEPS模型[22]、IBIS模型[23]等；(2)集成现有水文模型和生态模型构建的集成模型[24]，主要有RHESSys模型[25]、DLEM模型[26]、WaSSI-C模型[27]等.其中，WaSSI-C模型(Water supply stress index_carbon model)融合了Hamon的潜在蒸散(PET)模型[28]、Mccabe的融雪模型[29]以及萨克拉门托模型[30]，已在我国和美国的部分地区取得了很好的应用效果[31-32].本研究将干旱绿洲区的焉耆盆地作为研究对象，针对研究区的自然地理和资源环境特征，在原WaSSI-C模型的基础上增加冰川融化计算模块，以提高WaSSI-C模型在研究区的适用性，进而模拟并分析研究区水碳资源特征，以期为焉耆盆地水碳资源综合管理提供可靠的科学数据基础和决策工具支撑.

1 材料与方法

1.1 研究区概况

焉耆盆地是新疆天山山脉中的一个断陷盆地(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；光热资源丰富，具有典型的干旱区绿洲气候特征.盆地内主要有4条河流流经，即开罗河、黄水沟、清水河和孔雀河.在行政区划上，焉耆盆地位于新疆巴音郭楞蒙古自治州内，包含了焉耆、和静、和硕和博湖县的平原部分.

图1 研究区示意图Figure 1 Diagram of the study area

1.2 数据获取与处理

WaSSI-C模型运行的输入数据包含地形数据、气象数据、叶面积指数(LAI)数据、土地覆被数据等，验证数据包括径流数据、蒸散(ET)数据和总生态系统生产力(GEP)数据等，具体情况如下.

(1)数字高程模型DEM为美国国家航空航天局(NASA)的ASTER-GDEM V2产品，利用ArcGIS软件水文分析工具，经填洼、流向流量分析、流域分水岭提取等处理过程，将研究区划分为42个水文响应单元.

(2)土地覆被数据来源于美国陆地卫星5号(Landsat 5)，通过人工解译，并利用ENVI软件进行监督分类，依据国土资源部组织修订的国家标准《土地利用现状分类》(GB/T 21010-2017)，将焉耆盆地划分为建设用地、耕地、林地、草地、未利用地和水域6种土地类型，生成2005年焉耆盆地土地覆被数据.

(3)气象数据来源于国家气象局，包括月尺度的降水和气温数据，利用ArcGIS软件的地统计模块，引入经度、纬度和高程作为协变量，采用协同克里金插值工具得到1 km×1 km的降水和温度栅格数据集.经处理后可得2000～2015年不同月份的平均降水量和气温.

(4)水文数据为月径流监测数据，来源于焉耆盆地内焉耆水文站(图1).由于模型模拟结果为径流深数据，因此利用ArcGIS软件的曲面面积计算工具计算出流域曲面面积并将流量数据转化为径流深数据，从而用于模型的验证.

此外，WaSSI-C模型用到的其他数据包括：MODIS叶面积指数数据(LAI)(MODIS15)、MODIS蒸散数据(ET)(MODIS16)和MODIS总初级生产力数据(GPP)(MODIS17)，来源于美国国家航空航天管理局网站，地址为http://modis.gsfc.nasa.gov.在生态系统尺度上，GPP与GEP近似相等，二者均表征生态系统的碳生产能力，因此本研究将MODIS的GPP数据作为总生态系统生产力(GEP)的验证数据.

1.3 模型评价方法

本研究采用的评价指标是决定系数(R2)和效率系数(NS)[31，33].决定系数指示了回归直线对样本数据的拟合程度，其值在0～1，决定系数越高说明模拟效果越好.效率系数常被用于评价水文模型模拟效果的优劣，其值最大为1，最小为负无穷，且越接近于1，模型计算结果越好.

(1)

(2)

式(1)、(2)中，Qp为模拟值，Qm为观测值，Qavg为观测值的平均值，n为观测数据个数.

2 WaSSI-C模型与改进

2.1 WaSSI-C模型理论框架和运行机理

WaSSI-C模型是以月为时间尺度，以水文响应单元为空间尺度，旨在模拟水碳资源动态变化过程的生态系统水碳耦合模型，核心模块是水量平衡模块和生态系统生产力模块[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为月均叶面积指数；P为降水量(mm)；a、b、c为经验参数.

2.2.2 融雪过程

Pin=Prain+SnowM

(5)

式(5)中，Pin为总来水量(mm)；Prain为降雨量(mm)；SnowM为地表积雪融化量(mm).

2.2.3 水循环过程 萨克拉门托土壤湿度计算模型和Mccabe的融雪模型是WaSSI-C模型水循环计算过程的思想来源，基于土壤水分含量在垂直方向上的差异性，将其划分为上下两层，进而模拟土壤水分的循环过程.水循环过程主要包括蒸散发过程、地表地下产流过程和壤中流过程.其中，土壤上下层的张力水是蒸散发过程的主要水分来源，地表径流产生后直接汇入河网，经线性水库的调节后，地下径流和壤中流依次汇入河网，地表地下产流和壤中流相加并扣除一定时段内的蒸散发后，得到总入流量，经河网调节后，得到最终径流量.在该过程中，自然地形和植被状况对土壤水分循环的影响得到了充分的考虑，提高了模型模拟结果的准确性.

图2 WaSSI-C模型的理论框架和逻辑结构Figure 2 The theoretical framework and logical structure of the WaSSI-C model

2.3 生态系统生产力模块

生态系统生产力模型旨在模拟生态系统的碳循环动态过程，其构建思路主要是基于全球通量网络中的GEP和ET数据，利用线性回归模型建立回归方程.因此回归方程斜率代表了基于GEP的水分利用效率，基于此可评价和分析研究区的碳循环特征.

GEP=k×ET

(6)

式(6)中，GEP为总生态系统生产力；ET为实际蒸散；k为经验参数.

2.4 WaSSI-C模型的适应性改进

焉耆盆地处中国西部内陆干旱区，远离海洋，降水稀少，其水资源主要来源于高山冰雪融化.然而在WaSSI-C模型的融雪过程中，区域总来水量主要由降水量和地表积雪融化量组成，这显然与研究区实际情况不符.因此，在融雪过程中加入冰川融化量计算模块[35-36]，以提高模型在焉耆盆地的适用性.

SmowI=f×(T×9)t

(7)

Pin=Prain+SmowM+SnowI

(8)

式(7)中，SnowI为冰川融化量；f，t为经验参数.

3 结果与分析

3.1 WaSSI-C模型模拟结果评价与适用性分析

本研究将2000～2009年作为WaSSI-C模型的率定期，2010～2015年作为验证期，利用决定系数(R2)和效率系数(NS)，结合水文站的实测径流数据和MODIS的ET和GEP产品数据[37-38]，分别对两个时期的总径流模拟效果、蒸散模拟效果和总生态系统生产力模拟效果进行定量评价.

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数据缺乏对水域下垫面蒸散值的计算结果值，因此在本研究中仅对除博斯腾湖区域(水文响应单元41)外其他区域的蒸散模拟结果值进行对比验证.率定期模型的模拟蒸散值与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 总生态系统生产力(GEP)模拟结果验证 率定期GEP的模拟值与MODIS值的R2和NS分别为0.88和0.84，GEP模拟值与MODIS值总体上呈现出较好的一致性，仅在个别月份(2月、11月和6～8月)有较小的差异(图7)；验证期GEP的模拟值与MODIS值的R2和NS分别为0.87和0.82，仅在6～8月和10月有较小差异，在其他月份MODIS和模型计算的GEP值吻合程度很高(图8).造成模型在个别月份的GEP模拟值与MODIS值之间存在较大差异的原因可能是：(1)MODIS的GEP计算时利用了辐射数,而本研究使用了温度，二者间尽管有较强的相关性，但并不完全等价[42].(2)MODIS基于全球气象数据来计算GEP，而本研究利用基于观测点的插值气象数据，二者间的差异增加了模型的不确定性.

图7 率定期焉耆盆地不同月份平均总生态系统生产力(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 率定期焉耆盆地不同月份平均总生态系统生产力(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，在空间分布上，总体呈现出“西高东低，北高南低”的特点(图9).研究区水分补给主要来源于高山冰雪融水和降水.在盆地西部和北部(水文响应单元1～7、9、11～13、16～18、22、24、26～31、34～35、37)由于海拔较高，受冰雪融水补给最多，因此具有较高的产流量(97～220 mm/y)；而在盆地东部和南部(水文响应单元8、10、14～15、19、20～23、32～33、36、38～42)，其海拔较低，受冰雪融水补给较少，且由于上游用水截留的影响，导致在这一地区年均产流量较低(36～128 mm).另外，研究区水汽主要来自于西北方向[45]，造成盆地西北部的降水量较东南部更高，客观上加剧形成了这一特征.

图9 焉耆盆地年均径流空间分布Figure 9 Spatial distribution of annual average evapotranspiration in Yanqi Basin

3.2.2 蒸散空间分布特征 焉耆盆地陆地表面各水文响应单元年均蒸散在15～416 mm，博斯腾湖区域年均蒸散在416～776 mm，在空间分布上呈现出“中间高，四周低”的特点(图10).这可能是受到海拔和土地覆被类型的影响.在盆地西部和北部(水文响应单元1、2、4、5、7、9、11、12、15～18、20～22、24～28、31、34、35)由于海拔较高，导致其气温略低，且土地覆被类型多为未利用地，可用于蒸散的土壤水分较少，进而使得其年均蒸散量较低(15～282 mm)；在盆地南部和东南部(水文响应单元37～39、42)，主要是由于其产流量较低且未利用地比重较高，导致其年均蒸散量较低(15～282 mm)；而在盆地中部(水文响应单元3、6、8、10、13、14、19、23、30、32、33、36、40、41)，由于海拔较低，气温较高，且土地覆被类型以耕地和水域为主，涵养水源能力略差，进而导致其年均蒸散量较高(281～776 mm).

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

3.2.3 总生态系统生产力空间分布特征 研究区各水文响应单元年均总生态系统生产力在20～161 g C/m2之间，空间分布上呈现出与蒸散相似的特征，即中间高，四周低(图11).在盆地四周(水文响应单元1，2、4、5、7、9、11～13、15～18、20～22、24～31、33～39、42)，土地覆被类型多为未利用地，植被稀疏，使得其总生态系统生产力相对较低(20～118 g C/m2)；而在盆地中部(水文响应单元3、6、8、10、14、19、23、32、40)，土地覆被类型多为林地、耕地和草地，植被覆盖度高，因此总生态系统生产力相对较高(117～161 g C/m2).由于模型在水域的总生态系统生产力值设定为0，因此本研究暂不研究博斯腾湖区域(水文响应单元41)的总生态系统生产力.

图11 焉耆盆地年均总生态系统生产力空间分布Figure 11 The spatial distribution of the annual total ecosystem productivity in the Yanqi Basin

4 讨论

本研究的径流模拟结果年际变化在率定期呈现出先增后减的变化趋势，且在2002年达到最高，这与麦麦提吐尔逊·艾则孜等[14]关于焉耆盆地径流变化特征的研究结果保持一致，近16 a平均年径流量132 mm，与毕雪丽[39]的研究结果略有差异，造成这一问题的原因可能是：(1)本研究计算了2000～2015年的平均径流量，毕雪丽计算的是1961～2008年的平均径流量，时间跨度的不同可能会引起计算结果的差异.(2)补充的冰川融化计算模块受温度限制过大，高径流期模拟值与观测值的差异大于低径流期，进而使得模拟径流值总体上有一定高估.WaSSI-C模型在各水文响应单元年均蒸散模拟值为15～416 mm，焉耆盆地年均蒸散模拟值为220 mm左右(博斯腾湖区域除外)，与徐永明等[40]、郭玉川等[41]在焉耆盆地的蒸散量估算结果差异较小.WaSSI-C模型在焉耆盆地的年均GEP模拟值在80 g C/m2左右，与陈福军等[43]利用陆地生态系统碳循环模型( CASA模型)的估算结果基本一致，年内GEP模拟值在0.11～20 g C/m2，略低于王辉等[44]的研究结果，这可能是由于：(1)模型基于水文响应单元计算GEP均值，进而降低了数据间的差异性；(2)模型的GEP计算依赖于实际蒸散，而实际蒸散的计算受气温影响较大，进而导致在气温较低和较高时期，GEP模拟值不够准确；(3)王辉等的研究对象为整个开孔河流域，本研究仅针对焉耆盆地，研究范围的差异可能导致研究结果的差异.尽管模拟结果和验证数据间有着一定的不一致性，但这可能是由于输入数据和验证数据本身可能具有的误差以及模型的不确定性对模拟效果评价的影响导致的，而且，模型计算结果在径流、蒸散和生产力方面的数据范围与相关研究大体一致.因此可以认为WaSSI-C模型经适用性改进后能够很好地应用于焉耆盆地.

另外，焉耆盆地水碳资源的空间分布的模拟结果与徐永明、郭玉川等[33-34]关于焉耆盆地的蒸散空间分布呈现出中部高于四周的特征极为相似，与王辉等[37]关于开孔河流域NPP高值在其研究区内呈带状分布，并向四周扩散衰减的研究结果保持一致.由此可见，改进后WaSSI-C模型对于焉耆盆地的水碳资源空间分布的模拟结果具有较高的可信度.

5 结论

本研究基于WaSSI-C模型理论框架和运行机理的深入探讨，针对焉耆盆地自然地理和资源环境特征，通过增加冰川融化计算过程对模型进行了适用性改进，并分析了焉耆盆地2000～2015年的水碳资源空间分布特征，得到了如下结论.

1) 改进后的WaSSI-C模型能够很好地应用于焉耆盆地.统计评价指标R2和NS均显示出WaSSI-C模型对研究区的水、碳模拟效果良好.

2) 焉耆盆地的水碳资源在空间分布上具有明显的不均衡性.总体上，径流空间分布呈现“西高东低、北高南低”的特点；蒸散和总生态系统生产力均呈现“中间高，四周低”的分布特点.具体表现为：盆地西部和北部具有高产流量、低蒸散量和低碳生产力；中部具有低径流量、高蒸散量和高碳生产力；南部和东南部具有低径流量、低蒸散量和低碳生产力.