雅鲁藏布江流域极端降水模拟及预估

2021-05-24高佳佳

冰川冻土 2021年2期

高佳佳，杜军

（1.中国气象局成都高原气象研究所拉萨分部，西藏拉萨 850000； 2.西藏高原大气环境科学研究所，西藏拉萨 850000；3.西藏高原大气环境研究重点实验室，西藏拉萨 850000）

0 引言

极端降水是全球最受关注、影响最大的自然灾害之一，是短期气候预测研究的重点［1］。IPCC［2］曾指出，随着全球气候变暖，内陆地区的极端降水事件频率呈现出增加趋势。近年来，区域性洪水、干旱、高温、雨雪冰冻等极端事件频发，尤其是20世纪80 年代以后，频繁的极端事件给生态环境、经济发展和人民生活造成了严重影响［3-4］。因此，研究区域性极端降水事件对科学认识气候变化背景下水循环的时空演变，把握气候异常对极端降水的影响规律，评估水资源管理及区域水资源安全具有十分重要的理论和实践意义。

极端气候事件归根到底是气候极值问题，气候极值是极端事件产生的必要条件，极端事件发生发展的预测首先要考虑极值的分布规律。国内外部分学者研究表明，中国西部的极端降水天数呈增加趋势，尤其是西北地区［5-8］。Fischer 等［9］分析了珠江的极端降水分布特征，并估算了极端降水指数。Hong 等［10］认为海河流域的极端降水主要发生在38°N，大部分站点的降水呈现出减少趋势。程炳岩等［11］、江志红等［12］研究认为广义帕累托分布函数（generalized Pareto distribution，GPD）在重庆、中国东部的日降水模拟中具有更高的拟合度；李占玲等［13］基于GPD 函数分析了黑河流域的极端降水频率特征，得出该流域20 世纪60 年代发生的极端降水次数最多，90 年代以后次数较少。Eylon 等［14］运用极值理论分析了巴拿马运河的极端降水分布特征，并估算了极端降水的重现期及相应的置信区间。刘彩红等［15］运用CMIP5 模式指出，青海高原的降水有极端化的趋势，极端降水频次增加，强度增大。韩国军［16］、游庆龙等［17］运用统计方法计算出青藏高原极端降水大部分呈增加趋势，且逐年平均降水强度和逐年连续降水天数均有所增加，90 年代以来增加明显。

雅鲁藏布江发源并流经西藏高原，地理位置特殊，是世界上海拔最高的大河之一，平均海拔4 000 m 以上，是全球气候主要变化区与敏感区。流域沿岸为西藏主要农、牧业生产区，其洪涝和干旱的频繁发生导致了水资源分布不均，进而影响了流域的用水矛盾和生态环境的恶化，而极端事件的发生是对区域气候、环境变化的重要响应。目前对气候极值进行定量评估的方法以气候动力模式为主，从概率论角度对极端气候事件及可预测性研究并不多，尤其是预测方法。IPCC报告中特别强调统计方法对极端气候事件的重要性［2］，因此，本文基于广义帕累托方法（GPD）的分布参数模型，针对超出阈值的数据作为样本数据来建模，从气候极值的分布规律出发，揭示极端降水的发生发展规律，探索极端降水的可预测性，从而更好地预估极端事件，为提高防灾减灾能力提供科学依据。

1 研究区与数据

1.1 研究区域

雅鲁藏布江（简称雅江）全部在中国境内，横贯西藏高原南部，干流全长约2.1×103km，流域面积2.4×105km2，雅鲁藏布江干流河谷沿东西向的断裂带发育，流域呈东西向的狭长带，支流多而短小，较大支流有拉萨河、雅鲁藏布江帕隆藏布、易贡藏布、拉喀藏布、尼泽曲、年楚河等。干流在仲巴县里孜以上为上游，是高寒河谷地带。里孜到米林县派区为中游，支流众多，流量增大，河谷展宽，气候温和，水利条件较好，是西藏农业最发达的地区。下游位于林芝一带。截止国境线，年径流总量为1.1×1010m3，洪水由强降水形成，持续时间较长。流域东部地区降水量充足，年平均降水量超过900 mm，达到半湿润地区水平，西部地区降水量少，年均降水量不足100 mm，为干旱地区，整个流域的降水量从东至西呈现递减的趋势。

图1 研究区域位置示意Fig.1 Location of study area

1.2 数据来源

选取流域内9 个气象站近50 年（1959—2017年）5—9 月逐日降水资料作为研究对象。部分站点开始于1978年和1979年。

2 研究方法

2.1 GPD分布函数和密度函数

GPD 可以直接利用历年的原始数据，人为设置阈值，在设置好阈值后，以此为标准来抽取每一年超过此阈值的极大或极小值，即“超门限峰值POT”（peaks over threshold），可以提高估算精度［18］。具体为：

GPD的分布函数为：

其相应的密度函数为：

式中：ξ为门限值；a为尺度参数；k为形状参数，自变量x的取值，取决于k的值，当k< 0时，ξ 0 时，GPD 为指数分布；当k=1 时，GPD 在区间[ξ，ξ+a]上均匀分布。

表1 雅鲁藏布江流域内各站点资料长度Table 1 Basic information of each stations in Yarlung Zangbo River basin

2.1.1 阈值的选取

GPD 模型的核心在于确定阈值，它是正确估计参数的前提。如果阈值选取的过高，会使得超额数据量太少，导致估计出来的参数方差很大；如果阈值选取的过低，则不能保证函数的收敛性，所估计参数有较大偏差。本文主要使用Hill 图估计、百分位法和年交叉率法来确定各站点日降水量阈值。

Hill 图法是基于Hill 估计量的一种阈值图形法［19］，由点（k，1/H（k，n））构成的曲线，通过观察图中尾部指数稳定的区域来选择阈值。其定义为：

百分位法是将该站点的日降水量从小到大排序，并计算相应的累计百分位，某一百分位所对应数据的值就称为这一百分位的百分位数，文中分别计算了90、95、97、98、99百分位。

年交叉率法是假定超阈值降水极值出现次数服从泊松分布，以一年为时段所量度，极值超过阈值的次数，λ=n/T，其中，λ为年交叉率，T为资料总年数。

2.1.2 参数估计

对GPD 分布进行参数估计的方法有很多，本文主要通过极大似然估计法对参数进行估计。该方法具有很强的灵活性，可以适应不同模型的需求，且统计特性良好，能够综合各种有关信息到统计推断中去。模型估计时，样本变异可能会导致模型参数的不确定性，由于极大似然方法具有渐进正态性，容易给出估计值及其标准误差（标准误差是参数不确定性或变异性的度量之一）。该方法唯一的缺点是计算时迭代繁琐。具体方法见文献［20］。

2.1.3 重现期计算

某一指定重现期T时间的降水量分位计算公式：

2.2 独立性和稳定性检验

GPD 拟合需要超阈值数据序列满足平稳性的条件，因此，拟合之前需对超阈值序列进行稳定性检验。文章中使用Mann-Kendall（M-K 检验）对序列的变化趋势和突变点进行检验。M-K 检验［21］是气象学、气候学中经常用来进行突变检验的一种非参数检验方法，它不要求样本符合一定的分布。即给定显著性水平α=0.05，则统计量的临界值为±1.96。统计量大于0，表示序列呈上升趋势；反之，表明呈下降趋势，大于或小于±1.96，表示上升或下降趋势明显。该方法能有效区分某一自然过程是处于自然波动还是存在确定的变化趋势，常用于气候变化影响下的降水、干旱频次趋势检测。

2.3 拟合度检验

通过对超阈值序列进行平稳性检验后，即可对GPD 进行拟合，拟合结果需通过Kolmogorov-Simimov（K-S 检验）。一般在K-S 检验中，先计算需要做比较的两组观察数据的累积分布函数，然后求这两个累积分布函数的差的绝对值中的最大值D。最后通过查表以确定D值是否落在所要求对应的置信区间内。若D值落在了对应的置信区间内，说明被检测的数据满足要求。

3 结果与分析

3.1 阈值的选取

通过Hill 图尾部特征的稳定性来选取阈值（图2），由于Hill图具有较强的主观性，从而会导致选取的阈值不同。以中游地区的拉孜、拉萨和墨竹工卡站为例，拉孜站点的Hill 图尾部指数趋于稳定大致位于65～70 附近，对应的降水量为7.9～15.4 mm；拉萨和墨竹工卡站的尾部指数均位于160左右时趋于稳定，对应的降水量分别为8.7～18.5 mm 和9.9～20.3 mm。同样的，下游地区的林芝站在80～120位置时趋于稳定，对应的降水量为14.4～17.6 mm。为更好的确定阈值，在此基础上，我们结合百分位法和年交叉率法进行阈值确定（表2）。

图2 各站点Hill图指数分布Fig.2 The distribution of index of Hill plot for each stations

表2 各站点的百分位阈值（mm）和年交叉率Table 2 Precipitation threshold selection in GPD and the average annual occurrence number for each stations

与Hill 图相比，百分位法可以更精确的确定阈值。以拉孜站为例，Hill 图显示的阈值为7.9～15.4 mm，而对应的是该站点的95 百分位；同样拉萨站、墨竹工卡站的Hill 图确定的阈值均位于93 百分位和95 百分位。总体而言，Hill 图确定的阈值要小于百分位确定的阈值。根据前人研究结果得出，年交叉率为1～2 时对应的阈值可作为GPD 拟合分析时的参考阈值。结合表1 得出，当各站点日降水量达到99 百分位时，年交叉率均稳定在1.5 附近，因此我们确定99 百分位时的阈值为最佳阈值。

3.2 超阈值序列的平稳性检验

根据阈值的平稳性和稳定性要求，我们用M-K法在显著性水平α=0.05 条件下检测超阈值降水序列的变化趋势和突变。结果显示，下游地区的林芝和米林站的统计量略大于1.96，分别为2.3 和2.8，其余大部分站点均位于临界区域内，通过了显著性水平检验。通过时间序列曲线（UFk曲线）可以看出（图3），米林、墨竹工卡、南木林站呈显著增长趋势，说明三个站点的超阈值序列的日降水量呈逐渐增加趋势；其余站点呈下降趋势（此部分只给出拉孜、墨竹工卡、日喀则、米林站的趋势变化图，其余图表省略）。各站点超阈值序列没有明显的突变性。

图3 流域内四个站点的时间序列统计量变化图Fig.3 Time series statistics of four sites in the basin

3.3 模型拟合效果检验

通过对流域内各站点进行GPD 拟合（图4），并使用K-S检验看其是否符合已知理论分布函数。大部分站点的统计值小于0.01显著性水平，接受原假设。说明雅江流域各站点之间虽存在差异，但由GPD 拟合曲线可知，理论频数和实测频数基本相符。另外，虽然西藏地区大部分台站的观测记录起始年代不一致，如拉萨站、拉孜站，资料开始时间分别为1955 年、1977 年，但从图4 中可看出，数据满足方程需求，拟合结果表明资料长度并不影响降水极值的统计推断，且资料年限越长拟合结果越好。

图4 各站点累积频率和实际频率分布曲线对比Fig.4 The distribution of cumulative frequency and empirical frequency over the observation stations

3.4 雅江流域极端降水的分布特征

从各站点的极端降水阈值分布来看，林芝地区的阈值最大，阈值最小的站点为拉孜站。这与我们的观测事实一致，林芝地区日降水量大，连续降雨日数长，得到的阈值就大。

尺度参数主要是描述极值分布的变率，尺度参数越大，极值波动范围越大，表明打破极端降水的记录值也越大。整体而言，雅江流域的尺度参数由下游向上游是逐渐减小的，平均值为5.95。由表3可知，下游地区的尺度参数最大，约为7.00，表明这一区域的极端降水变化幅度很大。从气候背景来看，该地区位于高原季风区，受印缅槽和西风带影响，季节性降水较大［22-23］，5—9 月的降水总量可达600 mm，是西藏地区夏季降水量最大的区域，因此可能出现的破极端降水记录值要高于其他地区。尺度参数最小的区域位于流域中上游地区，受地理位置和大气环流影响，雨期短，降水量少，且连续降水日数也少，降水极值的范围比较小，区域打破极端降水的记录值要比下游地区低。这与前人研究结果一致［24］，我国的干旱地区大部分位于非季风区，降水极值范围较小，破纪录的可能性较季风区小。

表3 流域内各站点的GPD模型参数估计及检验Table 3 Estimation and validation of parameters in GPD model

形状参数作为模型的第二个重要参数，不同的形状有不同的尾部分布特征，它表示该区域极端降水的破纪录率。由表3 看出，形状参数正值区主要位于拉孜地区，这些地区发生破纪录降水事件的可能性比其他地区大。正是由于该区域降水日数少，所以一旦有降水过程，就可能会打破降水极值。而下游地区的形状参数为负值，说明这些区域的降水发生破纪录的概率偏小。因为夏季，降雨越频繁的区域，极值变率大，较均值离散程度大，则破纪录的概率较小。有降水的地区，不是有较大的形状参数就是有较大的尺度参数，不可能两个参数都大。形状参数较大的地区，称为“形状参数主导区”，该区域多持续降水，极端气候事件较少；尺度参数较大的区域称为“尺度参数主导区”，该区域降水较少，且多变，极端气候事件较多［23］。

3.5 雅江流域极值重现期预估

极值模型建立最重要的目的之一就是预测极端事件的重现期或重现水平。如表4 所示，根据重现期公式可以得出，从5 年一遇和10 年一遇的极端降水值来看，雅江流域除拉孜站外，其他地区降雨极值均超过30 mm，日喀则地区的降水极值达50 mm，其中拉萨、泽当、墨竹工卡和南木林站5 年一遇和10年一遇的极端降水量分别在40 mm 左右。相关研究指出［26-27］，当温度上升为2 ℃时，青藏高原的强降水距平百分率平均增多44.5%～59.5%，大值区出现在山南附近，这与我们的研究结论相符合。各站点在15年一遇的极端降水值之后，极值水平的增长变得非常缓慢，其中林芝地区的增长最缓慢，以0.32 mm·a-1的速率增长；日喀则地区的降雨极值增长率最快，约0.72 mm·a-1。西藏地区极端降水频率一般为每年4.3 次，强度在20 mm·d-1以上，林芝地区为极端降雨的高值区，且沿雅江一线极端降水的频次呈增加趋势［28-29］。

表4 雅江流域各站点日降水量极大值重现水平（单位：mm）Table 4 The maximum daily precipitation of flood season reappeared in Yarlung Zangbo River（unit：mm）

为检验各站点不同重现期水平的合理性，将其带入各站点逐日降水序列中进行验证。以拉萨和拉孜站为例，拉萨站5 年重现期水平值为39 mm，在1967—2017 年所有逐日降水中，共有10 次大于39 mm的降水过程，平均每5年一次。拉孜站5年重现期水平值为35.1 mm，在1977—2017 年期间，共有9 次大于35.1 mm 的降水过程，平均每4.4 年一次。由GPD 拟合计算出的极端降水重现期水平基本符合实际。值得注意的是这里的“重现期”并不意味着经过T年之后一定会出现的“周期”，它是概率意义上的“统计周期”。