《1、 引言》
1、 引言
气候条件和土地利用状况的变化已被证实对特定流域内水资源的时空分布有显著影响[1‒3]。这些变化通过改变降水、蒸散（ET）和相应的水流状态来影响相关的水文循环[4]，这可能会导致可利用水资源的减少进而带来水资源冲突的增加[5‒6]。如前所述，为理解分析气候变化和人为干扰活动下的水文循环，已经开展了许多研究[1‒7]。然而，由于水循环相关的复杂性，特别是对于水循环的核心过程即径流，研究面临许多挑战。一般来说，径流的识别和预测被认为是水文循环分析的必要前提。变化条件下的径流预测被认为是有助于水设施管理、水库调节和水资源规划的一个重要过程。因此，识别气候和土地利用影响下水文过程的动态特征变化对于水资源可用性评估和适应策略制定至关重要[8]。
此前，学者采用耦合了气候和土地利用变化预测的水文模型来评估多种情景下气候变化对水文过程的潜在影响[9‒11]。Wilby和Harris [12]利用两组参数、4组全球气候模式（GCM）输出以及两种水文模型评估了两种排放情景下气候因子对泰晤士河流域旱季现象的影响。Ahiablame等[13]设计了几个土地利用情景和气候模式，以分析可能的径流响应情况。Morán-Tejeda等[14]将一个区域气候模型和多个假设的土地利用情景与山区流域的水文模型相结合。一些研究还表明，气候变化是导致流量变化的主要因素之一[15]。水文模型通常用历史数据进行率定，然后通过未来气候和土地利用变化的预测作为驱动进行研究。随着环境变化的显著性增强，未来环境因子也将导致水文模型参数的动态变化。与气候变化对土壤特性的间接影响相比，土地利用变化的影响要显著得多，因为它直接而快速地改变了土壤的孔隙度和连通性[16]。对于土地利用变化，Merz等[17]指出可以认为模型参数是5年内不变的，而Luo等[18]则指出水文模拟应该考虑不同季节模型参数的不同，然后再考虑土地利用变化通过改变土壤特性对水文过程的间接影响。
通常，GCMs被广泛用作未来气候变化预测的主要工具。然而，由于气候系统的复杂性，GCMs在气候预测方面存在着多种不确定性。GCMs通过降尺度到合适尺度可以表示空间异质性[19]，通常采用动态降尺度和统计降尺度两种常用方法对GCMs进行预估[20]。统计降尺度旨在数据充足的情况下，建立历史GCMs与观测值之间的统计关系，在水文相关领域研究应用广泛[21‒23]。由于多种降尺度方法和GCMs选择的累积影响，水文预报的不确定性增加[23‒24]。多个GCM集成已被证明优于单个GCM，通常采用长短期记忆（long short-term memory, LSTM）和多模型集成均值（multi-model ensemble mean, MMEM）两种方法对GCMs进行集成[25‒28]。然而在区分“好”和“坏”的度量标准上几乎没有一致意见，但多个气候模式的集合通常优于任何单个模式，因为这可以提高未来气候预测的可靠性[29‒31]。使用机器学习模型或多个GCMs的平均值进行集成可以减少未来水文响应的预测偏差[32]。
此外，由于土地空间特征和利用类型的多样性，预测土地利用变化的模型较少。目前常用的方法有元胞自动机（CA）模型、马尔可夫（Markov）模型和土地利用转换及其效应（CLUE-S）模型。CLUE-S模型在小规模尺度应用中得到了广泛应用，但模型的类型和网格数量有限[33]。马尔可夫模型主要用于时间尺度的土地利用预测，但不能反映土地利用变化的空间分布，该模型通常与其他模型一起使用，以发挥其优势[34]。CA模型能够识别复杂系统的空间动态演化过程，但难以有效反映社会经济因素等宏观影响。因此，结合CA模型强大的空间计算能力和Markov模型的长期预测能力的CA-Markov模型具有很强的优势，且被广泛应用于土地利用预测[35]。此外，CA-Markov模型还可以考虑自然、社会、经济、行政中心、交通和气候等条件对土地利用变化的影响。然而，CA-Markov模型的研究大多集中在大城市，流域尺度尤其是山区的土地利用预测研究较少。
安宁河流域是长江上游的一条支流，是中国农牧业的交错带。由于其地形独特，地形起伏，海拔高，坡度陡，山谷生态系统极其脆弱。降水季节和区域分布不均匀，雨季暴雨强度大，旱季降水少，蒸散发强烈[36]。已有大量研究利用多种方法和技术对长江流域的气候变化进行了分析[37]，然而针对长江脆弱的干旱或半干旱地区评估气候和土地利用共同影响下的水文过程更加重要。
随着变化环境复杂性和不确定性的增加，安宁河流域水文过程研究面临着以下挑战：①地形复杂且脆弱地区进行水文模拟的稀缺性；②众多气候变化情景中寻找适用的气候预测方法；③未来土地利用的时空变化特征；④简化的水文模型和非线性水文过程之间的差异[30]。综上，时空变化尺度下长江干暖河谷地区气候和土地利用影响的研究是长期水文预测必须解决的重要科学问题。
因此，本研究的目的是：①评估MIKE SHE模型在山区的适用性；②通过LSTM和MMEM方法，运用31个GCMs对安宁河流域的降水和温度进行集成预报，以提高预报精度；③预测自然、社会和经济条件下的土地利用变化；④研究气候和土地利用变化对径流量和实际蒸散发的潜在影响。开发的模型结合了基于31个GCMs集成的气候预测模型和基于CA-Markov的土地利用预测模型，不仅可以解析气候和土地利用预测的不确定性，还可以研究水文过程的动态变化和范围。这对于处理随机不确定性和复杂变化条件下的水资源优化管理具有重要意义，为长江上游和下游地区的水资源和环境可持续性适应措施提供了基准信息。
《2、 研究区概况》
2、 研究区概况
横断山区干暖河谷具有独特的气候、地理和水文特征，基于全球气候变化和当地的物理条件，它具有典型的非地带性气候。Zhang等[38]发现，干暖河谷的气候变化与全球气候变化存在不一致现象，20世纪50年代以来气温下降，而降水增加。作为中国典型的生态脆弱区和重要的农业生产基地，变化条件下的水文分析有助于作物系统、生物资源和生态恢复的决策制定。
安宁河流域位于中国西南干热河谷，是长江上游的二级支流（图1）。其经度范围为102°06′E~102°10′E，纬度范围为26°38′N~29°02′N，海拔为900~4750 m。流域全长337 km，流经冕宁、西昌、德昌和米易，面积约11 150 km2。多年平均降水量为1240 mm，平均气温为17~19 °C，有明显的干湿季节，雨季从6月到10月，雨季降水量占全年降水量的90%以上，强降雨加上山坡陡峭，通常会导致山洪暴发，通常干季降水量少，蒸发量大。无论是雨季还是干季，降水都表现出明显的空间差异性，即山上降水多，山谷降水少。
《图1》
图1 安宁河流域图。DEM：数字高程模型；金沙江：长江上游。
作为凉山州的母亲河，安宁河流域是14个少数民族的聚集地。它是四川省第二大河谷平原，是长江上游生态脆弱的地区，也是中国西部经济发展的重要地区。多年来，由于全球气候和土地利用变化的影响，出现了许多生态环境问题，如水资源短缺、土地荒漠化、水质恶化等，已成为西南地区生态经济发展的严重障碍。社会经济活动相对集中，水利水电工程规划建设较多，灌区规模和城市规模较大，为流域系统管理带来的不确定性和难度加大。因此，分析气候和土地利用变化的水文响应对长江流域水资源合理利用、构建生态屏障和可持续发展具有重要意义。
《3、 方法》
3、 方法
《3.1 模型集成框架》
3.1 模型集成框架
研究主要目标是利用水文模型预测气候和土地利用变化影响下的水文响应。图2展示了水文响应预测的流程图。首先，根据气象水文资料的可获得性和连续性，利用1977—1983年数据对MIKE SHE进行率定，利用1984—1986年数据进行验证，并输入1980年的土地利用数据和湾滩站的日尺度流量数据进行率定和验证；随后，进一步利用2005年的土地利用数据和2006—2007年的气象水文数据对MIKE SHE进行验证，保持敏感性参数不变，模型结果表现良好。然后，利用LSTM和MMEM方法对31个GCMs集成4组气象数据，对未来气候变化进行了预测，其中利用1961—2000年的历史气象数据对LSTM方法进行训练，并利用2001—2010年的数据进行了测试。随后，利用CA-Markov模型结合多准则评价（MCE）和空间多元回归（LogisticReg）模块预测了未来3组土地利用数据。最后，我们利用4组气象预测数据和3组土地利用预测数据来驱动MIKE SHE模型以评估安宁河流域未来时期的水文响应。
《图2》
图2 气候和土地利用变化影响下的水文响应预测流程图。NSE：纳什效率系数；R：相关性系数。
《3.2 水文模型构建》
3.2 水文模型构建
选用基于物理机制的分布式水文模型的典型代表MIKE SHE模型，探讨气候和土地利用变化对水文过程的影响[39]。该模型主要考虑了植物截流、蒸散发、坡面流、非饱和区和饱和区、地表水和地下水交换及融雪径流等水文过程，由于这些参数具有物理意义，因此可以从测量数据中估计得出[40]。MIKE SHE模型主要在以下4个方面应用较广：①探索模型的适用性及模型构建过程[41]；②参数敏感性分析[42]、校准与验证[43]及模型不确定性分析[44]；③在数据不足地区建立分布式水文模型的可行性分析[45]；④全球变暖背景下水文对气候和土地利用变化的响应[46]。MIKE SHE水文模型在多领域、多尺度上应用广泛，尤其在湿润和干旱地区具有很大的应用潜力[47]。此外，不同于SWAT、VIC和TOPMOEL等分布式水文模型，研究人员可以根据研究区不同的气候条件、地质特征和可获得数据选择MIKE SHE中适合的模块进行建模[48]。
《3.2.1. 水文模型的输入数据》
3.2.1. 水文模型的输入数据
水文模拟和分析中的空间异质性主要通过数字高程模型（DEM）、气象数据、土地利用数据、土壤数据和水文数据等进行表征（表1）。气象数据选取出水口上游冕宁、西昌和米易三个气象站，并从国家气象中心获取了1961—2010年的降水、高温和低温站点数据。本研究中使用了联合国粮食及农业组织（FAO）的Penma-Monteith方法来估算潜在蒸散发[49]。使用《中华人民共和国水文年鉴》中湾滩水文站的日尺度流量观测值来率定和验证MIKE SHE。DEM和土地利用数据是从中国科学院的数据平台获取的，而土壤数据是从世界土壤数据库获取的。ArcGIS工具用于处理和获取研究区的空间数据信息，如流向、汇流、河网、土壤和土地利用数据等。
《表1》
表1 MIKE SHE模型输入数据
| Data type | Name | Data Source |
|---|---|---|
| DEM | ASTER GDEM （30 m × 30 m） | GS Cloud |
| Climate data | Precipitation | National Meteorological Center |
| Temperature | National Meteorological Center | |
| Reference evapotranspiration | Penman´s equation | |
| Vegetation | Land use map （30 m × 30 m） | Resource & Environment Cloud |
| Leaf area index | FAO and field collection | |
| Root depth | FAO and field collection | |
| Soil | Surface and sectional type （1 km × 1 km） | Harmonized World Soil Database |
| Hydrology | Streamflow | The People´s Republic of China Hydrologic Yearbook |
此外，水文模型的静态参数很难反映未来变化环境下的水文特征。植被覆盖是流域最关键的环境因素之一，它直接影响到降雨在蒸发和径流之间的分布。本研究采用植被叶面积指数（LAI）和根系深度（RD）作为模型参数动态变化的表征因子，并对植被生长期和非生长期进行划分。随后，基于降水和温度等气候变化条件，建立并分析LAI与RD之间的函数关系，估算未来每10年LAI与RD的动态参数。水文模型的动态参数能够建立模型参数与植被变化的响应关系。
《3.2.2. 率定和验证》
3.2.2. 率定和验证
参数校准是调整模型将模拟结果与观测数据进行比较的过程，如果所有参数都进行调整则会发生过参数化现象。在校准过程中。主要敏感性参数根据文献综述获得[50‒51]。MIKE SHE模型中的自动率定工具被用来进行参数校准（表2）。
《表2》
表2 MIKE SHE敏感性参数
| Parameter | Initial value | Value range | Final value |
|---|---|---|---|
| Manning coefficient | 20 | 5‒50 | 25 |
| Detention storage | 2 | 0‒20 | 4 |
| Horizontal hydraulic conductivity | 1.00 × 10-4 | 1.00 × 10-5‒2.00 × 10-3 | 1.25 × 10-3 |
| Vertical hydraulic conductivity | 1.00 × 10-5 | 1.00 × 10-6‒2.00 × 10-4 | 1.25 × 10-4 |
| Specific yield | 0.20 | 0.10‒0.50 | 0.11 |
| Storage coefficient | 1.00 × 10-4 | 1.00 × 10-7‒1.00 × 10-3 | 2.22 × 10-5 |
| Drainage level | -0.5 | -2.0 ‒ -0.1 | -1.2 |
| Time constant | 1.00 × 10-6 | 1.00 × 10-7‒1.00 × 10-5 | 3.45 × 10-7 |
纳什效率（NSE）系数和相关性系数（R）用来评估MIKE SHE模型的表现。纳什效率（NSE）系数对模拟值和观测值的均值和方差都很敏感。NSE和R系数的表达式如下：
（1）
（2）
式中，
NSE值范围为0~1，越接近1表示模拟效果越好，模拟精度越高。一般而言，当模拟值和测量值的相关性系数R大于0.6时，表明该模型效果较好，当NSE值在0.50~0.95之间时，表明模拟效果更好。
《3.3 多气候模式的气候预测》
3.3 多气候模式的气候预测
在本研究中，从CMIP5网站下载了31个1961—2050年月尺度的大气环流模型（GCMs）（表3）。运用Liu和Zuo [52]开发的NWAI-WG统计降尺度模型获取每个气象站的日尺度数据。数据集包括降水、最高温度和最低温度三个变量。政府间气候变化专门委员会报告（IPCC5）指出，温室气体排放是未来气温和降水不确定性的主要来源，因此代表中等和较高排放情景的RCP4.5和RCP8.5被用于气候变化影响研究[53]。
《表3》
表3 CMIP5中31个GCMs的基本信息
| Model | Region | Resolution | Model | Region | Resolution |
|---|---|---|---|---|---|
| ACCESS1.0 | Australia | 1.87° × 1.25° | GISS-E2-R | the United States | 2.50° × 2.00° |
| ACCESS1.3 | Australia | 1.87° × 1.25° | GFDL-CM3 | the United States | 2.50° × 2.00° |
| BCC-CSM1.1 | China | 2.80° × 2.80° | GFDL-ESM2G | the United States | 2.50° × 2.00° |
| BCC-CSM1.1（m） | China | 1.10° × 1.10° | GFDL-ESM2M | the United States | 2.50° × 2.00° |
| BNU-ESM | China | 2.80° × 2.80° | HadGEM2-AO | Republic of Korea | 1.87° × 1.25° |
| CanESM2 | Canada | 2.80° × 2.8° | INM-CM4 | Russia | 2.00° × 1.50° |
| CCSM4 | the United States | 1.25° × 0.94° | IPSL-CM5A-MR | France | 1.27° × 2.50° |
| CESM1（BGC） | the United States | 1.25° × 0.94° | IPSL-CM5B-LR | France | 1.89° × 3.75° |
| CESM1（CAM5） | the United States | 1.25° × 0.94° | MIROC5 | Japan | 1.40° × 1.40° |
| CESM1（WACCM） | the United States | 2.50° × 1.90° | MIROC-ESM | Japan | 2.80° × 2.80° |
| CMCC-CM | Europe | 0.75° × 0.75° | MPI-ESM-LR | Germany | 1.87° × 1.86° |
| CMCC-CMS | Europe | 1.86° × 1.87° | MRI-CGCM3 | Japan | 1.10° × 1.10° |
| EC-EARTH | Europe | 1.10° × 1.10° | NorESM1-M | Norway | 2.50° × 1.90° |
| FIO-ESM | China | 2.80° × 2.80° | NorESM1-ME | Norway | 2.50° × 1.90° |
| GISS-E2-H | The United States | 2.50° × 2.00° | MIROC-ESM-CHEM | Japan | 2.80° × 2.80° |
| GISS-E2-H-CC | The United States | 2.50° × 2.00° |
用大气环流模式进行大空间尺度上的气候预测存在显著的不确定性和模拟偏差，多个GCMs的集成通常被认为优于单一GCM [54]。采用LSTM和MMEM方法进行气候变化预测，可以提高多个GCMs的预测精度。MMEM是一种传统的利用多模式集合均值进行评价的方法，LSTM是一种时间循环神经网络，适用于分析和预估时间序列中间隔和延迟非常长的重要事件，已被证明比传统RNN更好地捕捉时间序列数据之间的关系[26]。通过引入LSTM方法检测时间信号，建立31个GCMs与各气象站观测数据之间的非线性相关关系，然后将长期记忆信息传播到未来时期。LSTM的主要目的是通过GCMs集成提高降水和温度的精度和质量。
LSTM的特殊结构是存储单元，用于存储信息和传播不同时间步长的输出。单元状态和输出由称为门的结构管理，包括遗忘门、输入门和输出门。门由带乘法运算的sigmoid函数（神经网络中的激活函数）组成。sigmoid函数从0到1不等，表示可以传递或丢弃多少信息。具体地说，首先，一个称为遗忘层的sigmoid层决定要从先前的单元状态丢弃多少信息，例如，0表示将丢弃所有信息；其次，一个称为输入门的sigmoid层结合双曲正切函数（即神经网络中的一个值为-1~1的激活函数），创建新的单元状态；最后，输出门根据单元状态输出决定输出的部分。
本研究利用31个GCMs分别建立了温度和降水预估的LSTM模式。首先，温度LSTM模型的输入包含31个平均温度，分别由上述31个GCM最高和最低日温度计算得出，建立1961—2010年日平均温度计算值与观测值的配对数据集。然后利用1961—2010年31个GCMs的月尺度降水和观测降水建立降水预报配对数据集。模型的训练时间为1961—2000年的40年，测试时间确定为2001—2010年的10年。将2021—2050年31个GCMs的温度或降水输入已开发的LSTM模型，分别预估未来时段内的温度或降水。评价标准NSE和R被用来表示预估的准确性。
如果能探测到1961—2010年31个GCMs的降尺度数据与观测数据之间的信号，就可以很好地实现2021—2050年的气候预测。利用LSTM进行深度学习和记忆，提高了预测精度。利用多个大气环流模式集成，以1961—2010年作为基准期，利用LSTM的记忆单元预测2021—2050年的降水和温度数据。相对而言，MMEM是直接对31个GCMs进行平均获得2021—2050年的降水和气温数据。
《3.4 CA-Markov模型》
3.4 CA-Markov模型
研究采用元胞自动机（CA）与马尔科夫（Markov）相结合的CA-Markov模型反映土地利用格局的时空动态变化。CA-Markov模型综合了多标准评价和多目标决策支持系统的优点，定义了土地利用类型间的转换规则。
马尔可夫模型是一种预测土地利用类型和随机转换过程的最优控制理论方法[35]。基于贝叶斯的土地利用预测如下：
（3）
式中，C(t)和C(t+1)分别是t和t+1时刻的状态；
（4）
CA模型是一种能够将时间、空间和状态进行离散的网格动力学模型，具有综合考虑自然、社会和经济条件等历史趋势模拟土地利用时空演变的能力。CA模型的计算方法如下[54]:
（5）
式中，t和t+1代表时间；C是细胞状态的集合；D代表某个邻域中所有细胞；f是细胞状态转换规则。对于本研究，D包含影响因素（即海拔、坡度、降水量、温度、河流、道路、人口密度等），f表示邻域内每个单元的组合规则，C(t+1)表示后一时刻土地利用状况受前一时刻C(t)的影响。
具体地，基于IDRISI软件中的CA-Markov模型，以1995年、2005年和2015年的土地利用地图为基础，预测了2025年、2035年和2045年的土地利用变化。具体过程如下：①计算1995—2005年和2005—2015年安宁河流域土地利用类型的转移面积矩阵和转移概率矩阵；②选择海拔、坡度、降水量、温度、河流等自然因素和道路、人口密度等人为因素作为土地利用预估的驱动因素，模块中的MCE和LogisticReg用于为每种土地利用类型创建过渡适宜性图集；③基于1995年和2005年的空间分布和驱动因子，模拟了2015年土地利用的空间变化，并利用2015年实际解译数据对模拟结果进行了验证，确保了模型的科学性和可靠性；④基于土地利用变化与各驱动因素之间的联系，分别采用模块中的三种方法对2025年、2035年和2045年的土地利用进行预测。
《4、 结果与讨论》
4、 结果与讨论
《4.1 MIKE SHE模型表现》
4.1 MIKE SHE模型表现
利用DEM，土地利用图，土壤数据，降雨，风速，日照时数，日最高、最低气温等数据对MIKE SHE模型进行校正。将流域出水口的日尺度流量与实测流量进行比较，以评估模型在变化环境下的表现。首先，输入1980年土地利用图，利用1977—1983年的湾滩水文站日尺度流量对MIKE SHE进行率定，再利用1984—1986年的逐日流量进行验证。然后，利用2005年土地利用图对2006—2007年的日尺度流量进一步验证（图3）。对于日尺度模拟，率定期的R和NSE值分别为0.85和0.72，第一个验证期的R和NSE值分别为0.87和0.73，第二个验证期的R和NSE值分别为0.84和0.70。各时段模拟值与实测值变化趋势一致，表明湾滩水文站径流模拟表现较好。对于月尺度模拟，率定期内R值和NSE值分别为0.88和0.79，两个验证期的相关系数R分别达到0.91和0.85，NSE效率系数分别达到0.86和0.74。结果表明，MIKE SHE能有效模拟安宁河流域的日尺度流量，月尺度的模拟效果优于日尺度，整体来说在山区具有较好的适用性。
《图3》
图3 MIKE SHE水文模型的日尺度模拟值与观测值的比较。（a）率定期（1977—1983年）；（b）验证期（1984—1986年和2006—2007年）。
从图3中可以看出，春冬季节流量相对较小，夏秋季节流量相对较大。从峰值模拟来看，某些时刻模拟值小于观测值。由于输入土地利用数据的变化，水文模型在第二个验证期表现不佳，说明土地利用变化也会影响水文过程。另外，由于研究区内气象站点稀少，降水等气象资料不能完全代表每个计算单元，这将对模型的模拟精度产生一定影响。
为了验证MIKE SHE模拟其他水文要素的适用性，对1977—1986年的实际蒸散发进行了模拟，并将其转化为研究区平均值。计算出研究区年平均实际蒸散发量约为580~720 mm·a-1，与国家青藏高原数据中心最新研究的范围一致，表明西南流域（含长江）实际蒸散发量应显著高于400 mm·a-1 [55]。2021—2050年的实际蒸散发将在接下来进一步分析和讨论。
《4.2 用于气候预估的多气候模式集成的表现》
4.2 用于气候预估的多气候模式集成的表现
研究运用LSTM和MMEM方法对RCP4.5和RCP8.5情景下的31个GCMs分析2021—2050年的气候变化情况，利用1961—2000年的历史数据和实测数据进行训练，用2001—2010年的数据检验模型表现。表4为1961—2010年两种方法下日尺度气温和月尺度降水量的表现。R和NSE系数表明，LSTM在多个GCMs集成方面表现较优，能够较好地反映基准期内安宁河流域的日温度，然而，仍然存在某些极端值无法捕捉到。由于90%的日尺度降水量为零或极低，LSTM不能很好地预测绝大部分的极值，也有研究表明，日尺度降水量变化不规则且波动较大，神经网络模型通常难以准确预测极高或极低的数据[56]。因此，将波动较大的日尺度降水转化为月尺度降水以消除零值，LSTM对月尺度降水则有较强的预测能力。此外，R和NSE的系数表明，月尺度降水预测比日尺度气温预测表现更好。
《表4》
表4 1961—2010年三个气象站日尺度气温和月尺度降水预测的准确性检验
| Station | Period | Temperature | Precipitation | ||||||
|---|---|---|---|---|---|---|---|---|---|
| LSTM | MMEM | LSTM | MMEM | ||||||
| NSE | NSE | NSE | NSE | ||||||
| Mianning | 1961‒2000 | 0.89 | 0.79 | 0.84 | 0.69 | 0.92 | 0.83 | 0.76 | 0.76 |
| 2000‒2010 | 0.85 | 0.72 | 0.83 | 0.68 | 0.90 | 0.79 | 0.82 | 0.80 | |
| Xichang | 1961‒2000 | 0.88 | 0.76 | 0.81 | 0.65 | 0.92 | 0.84 | 0.74 | 0.74 |
| 2000‒2010 | 0.81 | 0.64 | 0.79 | 0.62 | 0.89 | 0.76 | 0.79 | 0.75 | |
| Miyi | 1961‒2000 | 0.86 | 0.73 | 0.82 | 0.66 | 0.92 | 0.84 | 0.71 | 0.71 |
| 2000‒2010 | 0.84 | 0.71 | 0.81 | 0.65 | 0.88 | 0.76 | 0.76 | 0.74 | |
通过将降水量预估结果与实测值进行比较，发现LSTM方法对31个GCMs的预估数据总体上比MMEM方法表现好（图4）。LSTM预估得到的月平均降水和日平均气温与观测数据相关性整体较强，能够较好地预测基准期的气温和降水。对于西昌多年月平均降水和温度来说，MMEM预估的降水和温度曲线更接近观测值，特别是5~9月，然而对于1961—2010年的月尺度降水和月尺度气温来说，R和NSE表明LSTM预估的气候数据在三个站点的表现优于MMEM。结果表明，先进的深度学习方法LSTM在处理时间信号方面比传统方法MMEM更有效。
《图4》
图4 分别由MMEM和LSTM方法预测的基准期月平均降水和温度（1961—2010）。Pre_O：降雨观测值；Pre_M：MMEM预估的降雨；Pre_L：LSTM预估的降雨；Tem_O：温度观测值；Tem_M：MMEM预估的温度；Tem_L：LSTM 预估的温度。
《4.3 未来气候变化预估》
4.3 未来气候变化预估
以历史时期（1961—2010年）为基准期，采用LSTM和MMEM方法分别分析了2021—2050 年气温变化趋势（图5），并生成了4组气候变化情景（RCP4.5_L、RCP8.5_L、RCP4.5_M和RCP8.5_M）。从月平均气温的角度看，不同的集成方法可以观测到不同的变化趋势。不同气候情景下LSTM在10~4月的预测值普遍高于历史值，而在5~9月则低于历史值。冕宁、西昌气象站MMEM值与历史气温相比偏高，而米易站MMEM值与历史气温相比偏低。在RCP4.5和RCP8.5两种温室气体排放情景下均表现出不同程度的气温上升趋势，其中RCP8.5情景比RCP4.5情景下的温度略高。冕宁、西昌在RCP4.5和RCP8.5情景下的年平均气温均较基准期升高，后者较前者升高约1.4 ℃，这与IPCC5报道的结果一致。流域温度的升高必然导致安宁河流域源头冰雪的快速融化，进而影响产汇流机制发生变化，这将给当地水资源开发管理带来新的机遇和挑战。
《图5》
图5 预估期（2021—2050年）的月平均温度。Hist：历史时期1961—2010年；RCP4.5_L：由LSTM预估的RCP 4.5情景下的温度；RCP8.5_L：由LSTM预估的RCP 8.5情景下的温度；RCP4.5_M：由MMEM预估的RCP 4.5情景下的温度；RCP8.5_M：由MMEM预估的RCP 8.5情景下的温度。
图6展示了两种方法（MMEM和LSTM）及两种情景（RCP4.5和RCP8.5）下2021—2050年预估的月平均降水量的变化趋势。在冕宁站和西昌站，LSTM预估的降水量高于MMEM。未来一段时间内，两个情景下的月降水量与5~10月的历史降水量相比有较大波动，其余月份变化不大。结果还表明，RCP4.5和RCP8.5情景下的年平均降雨较历史时期均有所增加，且RCP4.5比RCP8.5情景增加约1.9 mm，尤其是在5~10月更为显著。
《图6》
图6 预估期（2021—2050年）的月平均降雨量。
《4.4 土地利用变化预估》
4.4 土地利用变化预估
1980年、1995年、2005年和2015年土地利用主要分为林地、耕地、草地、水域、建筑用地五大类型（表5）。根据遥感影像资料的可获得性，运用研究区1980年的土地利用图进行率定（1977—1983年）和验证（1984—1986年）。本研究主要利用1995年、2005年和2015年的土地利用数据生成转换矩阵预估每10年的土地利用数据，前15年时间间隔（即1980—1985年）对未来土地利用预测影响不大。从表 5可以看出，安宁河流域主要土地利用类型是森林，其次是农田和草地。结果表明，近几十年来耕地面积持续减少，这可能是由于90年代以来国家实施“退耕还林”等自然保护工程造成的，水域变化的原因可能是由于90年代以来流域上游大桥水库的建设和运行造成面积的增加引起的。
《表5》
表5 流域不同土地利用类型及1980—2015年变化率
| Type | Land-use area (km2) | Change rate （%） | |||||
|---|---|---|---|---|---|---|---|
| 1980 | 1995 | 2005 | 2015 | 1980‒1995 | 1995‒2005 | 2005‒2015 | |
| Cropland | 2673 | 2661 | 2636 | 2613 | -0.45 | -0.94 | -0.87 |
| Forest | 5990 | 5978 | 5978 | 6003 | -0.20 | 0 | 0.42 |
| Grassland | 2233 | 2248 | 2254 | 2225 | 0.67 | 0.27 | -1.29 |
| Water | 110 | 112 | 128 | 127 | 1.82 | 14.29 | -0.78 |
| Built-up | 57 | 64 | 67 | 95 | 12.28 | 4.69 | 41.79 |
在综合考虑土地利用变化驱动因素的基础上，选取海拔、坡度、降水量、气温、河流、道路、人口密度等因子，对不同土地利用类型的转化进行约束和限制，生成适宜性图像。采用MCE和LogisticReg代表性模块，结合CA-Markov模型，对3组土地利用预估驱动因子与土地利用的关系进行了评价。为了检验CA-Markov模型预估各土地类型数量和空间分布的准确性，根据1995年和2005年的土地利用数据预估了2015年的土地利用。然后利用IDRISI软件中Crosstab模块对1995年、2005年和2015年3组土地利用预估图和实际土地利用图进行比较，得到Kappa系数分别为0.96、0.89和0.91。结果表明，CA-Markov模型总体上是一个可靠的模拟模型，对未来土地利用预估具有较好的适用性和较高的精度。
基于2015年预估图与观测图校准的满意度，运用马尔可夫模型中考虑多个驱动因素的3种算法（即情景1：MCE；情景2：LogisticReg；情景3：Markov），预估2025年、2035年和2045年土地利用变化，每幅土地利用影像代表10年的土地利用状况（即2025年代表2021—2030年；2035年代表2031—2040年；2045年代表2041—2050年）。从图7可以看出，建设用地总体上有3种扩展情况，主要集中在住宅区和工业区以及交通用地等，尤其是中下游区域。主要趋势是由自然用地类型向人工用地类型转变，受区域经济发展、人口增长等因素影响，这种转变是合理可行的，并且所有情况都可能在未来发生。伴随着土地利用多样化和空间复杂化发展，必须加强对农田、森林和水体的保护，实现土地合理利用。
《图7》
图7 未来三个时期的土地利用预估情景。（a）MCE方法；（b）LogisticReg方法；（c）Markov方法。
《4.5 水文响应评价》
4.5 水文响应评价
《4.5.1. 径流对未来气候变化和土地利用的响应》
4.5.1. 径流对未来气候变化和土地利用的响应
在上述4组气候预测数据和3组土地利用预测数据的基础上，利用MIKE SHE模型对未来时期的水文过程进行了评价。以12种情景的交叉组合作为水文模型的输入，分别预估了气候和土地利用变化下2021—2050年的径流变化趋势以及蒸散发时空变异性。图8给出了不同情景下模拟的年平均流量，可以看出LSTM预估的径流比MMEM预估的高。我们可以看到，在3种土地利用预测下，径流没有明显的变化，但对气候变量的变化更加敏感。气候变化是历史时期径流变化的主要原因，这与前人的研究一致。从多年平均排放量来看，RCP4.5情景下的模拟值高于RCP8.5情景下的模拟值，特别是2031—2040年。尽管RCP4.5和RCP8.5温度没有显著差异，但由于降水变化，RCP4.5和RCP8.5的未来流量较高。由此也说明降水量对水文变化影响较大，MMEM降水预测数据较低，因此径流预测值远低于LSTM。结果表明，年径流量随降水量的增加而显著增加，随气温的升高而减少。
《图8》
图8 未来不同情景下的径流变化。（a）LSTM预估的气候变化；（b）MMEM预估的气候变化；LU1、LU2和LU3：3种土地利用情景。
为了进一步分析未来流量变化的多年月平均分布，将预估期（2021—2050年）RCP4.5和RCP8.5情景下的多年月均径流量与基准期（1977—1986年）进行了比较（图9）。我们可以清楚看到，月径流量从1~4月明显增加，而从8~11月则急剧下降。显著增加的原因可能是由于温度的轻微升高，这将不可避免地导致安宁河流域源头冰川和积雪的融化。如图9所示，5~9月的月平均流量比其他月份的变化大，这表明在此期间流量预测的不确定性较高。可以从两个方面来解释9~10月不同排放情景下流量略有下降趋势的原因：一是由于特殊地理条件下坡陡谷深，校验的水文模型在洪水期峰值表现不佳；二是虽然降水量在这一时期略有增加，但实际蒸散发也会受到其他因素的影响，如气温升高、地表类型、云量、日照时数、风速等。径流可能在雨季尤其是特别是5~8月增加，应防范全球环境变化尤其是气温变化情况下带来的冰川和积雪融化，导致洪涝灾害等事件的发生。与基准期相比，旱季的径流量特别是2~4月可能增加，然而11月份径流量的减少表明干旱风险将会增加。
《图9》
图9 预估期（2021—2050年）相对于基准期（1977—1986年）的月平均径流量变化。（a）RCP4.5；（b）RCP8.5。
《4.5.2. 实际蒸散发的时空尺度变化》
4.5.2. 实际蒸散发的时空尺度变化
为了更好地评估变化环境下的水文过程，对实际蒸散发进行了分析。本研究基于Yan和Smith公式算法，通过MIKE SHE的双层UZ/ET模型输出实际蒸散发结果[57]。图10显示了不同气候和土地利用情景下实际年蒸散发的空间变异性。通过MIKE SHE的模拟结果，对比研究了2021—2050年两种气候情景和3种土地利用情景下实际蒸散发的动态变化。实际蒸散发随温度升高而增大，由于RCP8.5情景下的温度较高，实际蒸散发的增幅[图10（c）和（d）]高于RCP4.5情景下的增幅[图10（a）和（b）]。研究表明，实际蒸散发对气候变化比对土地利用变化更为敏感。
《图10》
图10 RCP4.5和RCP8.5情景下预估期（2021—2030年、2031—2040年和2041—2050年）相对于基准期（1977‒1986年）实际蒸散量的变化。（a）、（b）分别代表RCP4.5情景下由LSTM和MMEM预估的气候情景；（c）、（d）分别代表RCP8.5情景下由LSTM和MMEM预估的气候情景。
图11展现了每种情景下每10年的平均径流量和平均实际蒸散量的比较。RCP4.5情景下的年平均径流量高于RCP8.5情景，除RCP8.5情景下的第一种土地利用情景外，其他情景下的实际蒸散发均表现出一致的变化趋势。水分平衡组成成分包括降水、林冠截留、地表径流、非饱和和饱和流等。显然，降水变化与径流和实际蒸散发的变化密切相关。并且在一定温度范围内，气温升高会影响蒸散发，降低林冠阻力，其中，在干旱和半干旱地区，蒸散量受到土壤水分正向影响[11]。由于未来气候和土地预估的不确定性及地形条件的复杂性，实际蒸散量与径流量之间没有绝对一致或相反的关系。
《图11》
图11 预估期每10年平均径流量和实际蒸散量的变化比较。（a）RCP4.5；（b）RCP8.5；C1、C2：两个气候情景；L1、L2、L3：三个土地利用情景；P1、P2、P3：分别代表2021—2030年、2031—2040年、2041—2050年。
《5、 结论》
5、 结论
本研究揭示了安宁河干暖河谷流域气候和土地利用变化对水文过程的影响。为解决单一气候环流模式产生的不确定性问题，本研究采用了先进的LSTM方法和传统的MMEM方法对多个GCM进行集成，获得了4组气候预测数据。综合考虑了海拔、坡度、降水、温度、河流、道路和人口密度等驱动因素，生成了3组土地利用预测数据。在数据可获得性和地形异质性的约束条件下，研究了复杂变化条件下安宁河流域水文过程的动态变化和范围。
结果表明，MIKE SHE水文模型可有效应用于青藏高原干暖河谷流域进行水文模拟。相比传统的MMEM方法，LSTM方法在日尺度温度和月尺度降水量预测方面表现更好，然而当降水为零时，LSTM对日尺度降水的预测难度较大。RCP8.5情景下未来气温的预测略高于RCP4.5情景，然而6~10月未来降水的预测略低于RCP4.5。年径流量随降水量的增加而显著增加，随气温的升高而减少。径流和实际蒸散发对气候变化比土地利用变化更敏感。径流量可能在雨季的5~8月份增加，而在旱季的11月份减少。实际蒸散发随温度升高而增大，且在RCP8.5情景下增幅较大。随着未来的环境变化，应采取有效措施应对5~8月的洪涝灾害以及旱季特别是11月的干旱灾害。作为长江上游生态环境脆弱区的重要支流，研究复杂地形区水文的时空变化和随机水文过程具有十分重要的意义。本研究可为应对环境变化下的洪涝干旱风险提供有效的指导，有助于长江上游地区和下游地区的可持续发展。
不断变化的环境对生态脆弱地区和长江下游地区的水资源和环境可持续发展提出了新的挑战，有望开发一个具有更多动态参数的水文模型以探索未来复杂和不确定的变化环境。
京公网安备 11010502051620号
