基于分位数回归的天山1号冰川径流变化分析

孟 晴,范永辉

(天津师范大学 数学科学学院,天津 300387)

冰川径流是我国西北干旱地区重要的水资源供给[1],由于全球气候变暖,世界上的大多数冰川长期处于萎缩状态[2],而冰川的加速消融使高寒地区的水循环加剧,导致水资源匮乏.受人类活动的影响,冰川流域径流的变化可能会使得流域冰川灾害发生的频率增加.因此,通过对冰川区径流的长期变化进行分析,了解其变化趋势,有助于更好地对西北干旱区水资源进行合理的开发利用[3].

对于冰川径流的研究,大都是以模拟为主,学者们通过使用各种水文模型,将复杂的冰川水文过程进行简化来探索冰川径流的变化趋势及其影响因子[4-5].分位数回归模型可对不同分位点下的因变量随自变量的影响变化情况进行分析,其结果不易受到异常值的影响,且对于模型中的残差项并不要求其是正态的,与常规的线性回归模型相比,分位数回归具有更强的稳健性与适应性.该方法最早由Koenker和Bassett[6]提出,目前已被广泛的应用到经济、医学、教育、生物学等领域.近年来,学者们将这一方法应用推广到了水文气象领域,已经取得了一些成果.Mazvimavi[7]利用分位数回归方法分析了津巴布韦地区近100年的极端降雨变化,没有证据表明在1892～2000年期间干旱的频率和严重程度发生了变化;Wang等[8]使用二阶平稳性分析、变点分析和分位数回归法对美国东南部地区的降水情况进行了对比分析,其中分位数回归法的结果显示出降水在低(高)分位数有增加(减少)的趋势;王永珂等[9]采用分位数回归法研究了娘子关泉降水及径流变化,结果表明年降水序列在高分位点的下降趋势比低分位点更为显著,年径流量序列在不同分位水平下保持一致性,都随时间显著下降;李婧等[10]使用分位数回归模型和GAMLSS模型对渭河及珠江流域的洪水进行非一致性频率分析,表明分位数回归模型的拟合优度整体优于GAMLSS模型.

天山乌鲁木齐河源1号冰川(简称1号冰川)是世界冰川监测组织(WGMS)设立的全球重点观测的10条冰川之一,是中国西部干旱地区的代表性冰川[11].我国在1号冰川的径流量方面已经有了一定的研究成果,孙美平等[12]采用数理统计、Morlet小波分析和 Mann-Kendall 突变检验对1号冰川1959～2008年径流量序列进行分析,表明这50年来1号冰川融水径流显著增加,并于1993年发生突变;贾玉峰等[13]研究了1号冰川流域径流及其组分长期变化,发现1号冰川流域径流主要由冰川径流和非冰川区降水径流组成,分别占 70% 和 30%;何天豪等[14]将水稳定同位素信息与冰川流域水文模型FLEXG相耦合,建立了FLEXG-iso模型,并对1号冰川年径流进行模拟检验,结果发现模拟效果良好.为了分析1号冰川径流随时间的变化趋势,本文采用分位数回归方法,在选取的高、低分位数水平下对1989～2017年间1号冰川年径流量序列的变化特征进行研究,并进一步探讨其与气候变化之间的响应关系.

1 研究区概况与数据

1.1 研究区概况

天山乌鲁木齐河源1号冰川位于中国新疆维吾尔自治区天山乌鲁木齐河的源头(43°05′N,86°49′S)(图1),根据中国冰川编目资料,于1962年对1号冰川(编号5Y730C0029)进行了首次调查和测绘,其面积为1.95 km2,2017年退缩为1.54 km2[15].在全球变暖的情况下,冰川一直处于持续的消融状态,从而使得1号冰川在1993年分裂为东、西两个独立的分支.该冰川是典型的大陆性冰川,是中国观测时间最长,资料最为详尽的冰川,根据大西沟气象站数据可知,1959～2017年平均年降水量为468 mm,其中5～9月份降水占全年总降水量的90%[12],7～8月平均气温为5 ℃.

图1 天山乌鲁木齐河源1号冰川、水文站及气象站位置示意图Figure 1 Location diagram of glacier No. 1, hydrology station and weather station in the Urumqi River Source of Tianshan Mountains

对于1号冰川的观测始于1959年,此前已有研究表明,1号冰川水文点冰川径流量约占总径流量70%[16],故本文使用1号冰川水文点的径流数据来进行相关分析.1号冰川水文点位于1号冰川的下游约300 m处,海拔3 659 m,其控制面积为3.34 km2,2017年冰川覆盖率为46.7%[15].随着全球气温升高,1号冰川的融水径流量呈不断增加趋势,尽管融水径流量的增加会给西北干旱地区提供丰富的水资源,但同时也存在着弊端.若径流量持续增加且没有减少的趋势,那么很有可能会导致一系列自然灾害的发生,因此精确研究冰川径流的变化有助于更好地预防自然灾害.

1.2 数据来源

径流数据采用的是天山乌鲁木齐河源1号冰川水文点的实测数据,来源于国家冰川冻土沙漠科学数据中心(http://www.crensed.ac.cn/portal/).由于夏季气温升高,融水径流会显著增加,故而选取大西沟气象站强消融期(7～8月份)气温和年降水量作为气象数据,时间范围为1989～2017年.该气象站海拔为3 539 m,距1号冰川的水平距离约2 km[17],且与河源区其他几个气象站的观测资料有着较好的一致性,可以很好地反映出河源区冰川区的气候变化情况[12].

2 研究方法

2.1 分位数回归方法

分位数回归方法不仅可以对被解释变量的条件期望进行分析,还可以分析解释变量与被解释变量分布的中位数、分位数之间的关系.分位数回归依据因变量的条件分位数对自变量进行回归,可得到不同分位点下的回归模型,其中,中位数回归是分位数回归的特殊情况,用对称权重解决残差最小化问题,而其他条件分位数回归则需要用非对称权重解决残差最小化[18].

假设随机变量Y的分布函数为F(y)=P(Y≤y),Y的τ分位数定义为满足F(y)≥τ的最小y值,即

Q(τ)=inf{y∶F(y)}≥τ,0<τ<1.

(1)

对于Y的一组随机样本y1,y2,…,yn,样本均值回归是使残差平方和最小,即

(2)

样本分位数回归是使加权残差绝对值之和最小,即

(3)

上式可等价表示为损失函数的形式,即

(4)

其中:ρτ(u)为损失函数,其定义为ρτ(u)=(τ-I(u<0))u,I(.)为示性函数.

2.2 分位数回归模型的参数估计及检验

对于一般的线性回归模型Y=Xβ+ε(ε为随机误差项),在给定X=x的条件下,Y的条件均值函数E(Y|X=x)=x′β的参数估计值为

(5)

同样地,一般的分位数回归模型为

Yi=X′iβ(t)+εi,i=1,2,…,n.

(6)

其条件分位数函数为Q(τ|Xi=xi)=x′iβ(τ),可通过求解

(7)

分位数回归模型的参数估计大致分为直接优化法与参数化法两大类,直接优化法又包括单纯形算法、内点算法及平滑算法.三种方法各有优缺点,可根据不同的情形选择合适的算法.单纯形算法相较于其他两种算法来说有更强的稳定性,故而在数据不多的情况下,使用单纯形算法来进行参数估计效果会更好一些.

为探究所建立的回归模型是否真正刻画了因变量与不同的自变量之间的依赖关系,需对回归系数进行假设检验,通常是使用t统计量来判断回归系数的显著性.

3 1号冰川径流量的回归分析

3.1 冰川径流量变化趋势分析

为了解不同的分位数下1号冰川径流量随时间的变化趋势,建立如下的分位数回归方程

Y=X′β(τ)+e

(8)

其中:Y为1号冰川径流量,X表示年份,β(τ)为对变量X进行参数估计的第τ分位数的系数,e为随机扰动项.

使用R软件对上述回归方程进行分位数回归,分位点选取0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,显著性水平a=0.05.对变量的回归参数进行估计和检验,具体结果如表1所示.

表1 模型(8)分位数回归参数估计及检验结果Table 1 Quantile regression parameter estimation and test results of model (8)

从表1可以看到,年径流量序列在各个分位点下的回归系数均为正值,表明1号冰川径流量在逐年增加.在不同的分位数水平下,年径流量序列的显著性有所差别,年径流量序列在中低分位点下有显著的回归趋势,但在高分位点下没有呈现出显著的回归趋势.

3.2 气温和年降水量对冰川径流的显著性分析

为探究冰川径流与气象因子(气温、降水量)之间的响应关系,本文利用1989～2017年间1号冰川水文点的径流量Y以及大西沟气象站的7～8月均温X1与年降水量X2,建立1号冰川径流量的分位数回归方程,即

Yt=Xt1β1(t)+Xt2β2(t)+et

(9)

其中:β1(τ),β2(τ)分别代表对气温和年降水量这两个变量进行参数估计的第τ分位数的系数,et为随机扰动项.

同样的选取0.1～0.9这9个分位点,使用统计软件计算出参数的回归系数以及检验参数的显著性.

从表2可看出,变量气温X1在低分位点的回归系数值是低于高分位点的,说明在气候条件相同的情况下,在低分位点下的径流要比在高分位点下的径流上升坡度小.观察变量X1各个分位点下的P值可知,大部分分位点下的P值均是小于0.05的,表明气温对冰川径流的影响比较显著,只有在分位点为0.2,0.3,0.4时是不显著性的.观察变量X2各个分位点下的P值可知,降水量对于冰川径流的影响并不显著,只有在低分位点0.1时呈现出显著的回归趋势,这是由于降水发生时时常伴随着降温,且冰川区的降水主要以固态为主[19],故而降水量对于冰川径流的影响较小.降水量的回归系数在各个分位点下没有明显的规律性,这可能是由于降水对于冰川径流的影响较为复杂.由此可以得出,气温是影响冰川径流变化的主要气象因子,降水对径流的影响不是十分明显.

表2 模型(9)的分位数回归参数估计及检验结果Table 2 Quantile regression parameter estimation and test results of model (9)

3.3 冰川径流量拟合结果

为进一步分析各个变量对径流量的影响,对径流量随时间、气温及年降水量的分位数回归进行了曲线拟合,从0～1间每间隔0.05选取一个值作为分位点τ,即为0.05,0.1,0.15,…,0.9,0.95,其中τ=0.5,时所作的回归直线是中位数回归.此外,为了便于比较,还使用普通最小二乘法拟合了一条线性回归直线,结果如图2～4所示.

图2 1989～2017年径流量随时间的变化趋势Figure 2 Variation trend of annual runoff over time during 1989～2017

从图2可以看出,年径流量序列在各个分位点下的回归直线都呈现出上升趋势,表明随着时间的迁移,径流量也在不断增加.图中的中位数回归直线与线性回归直线相对处于平行状态,在中位点0.5附近的几个分位点的回归直线比较集中,相反高分位点和低分位点的回归直线较为分散,可以看到回归直线的左右侧间隙相差不大,说明数据点总体上是均匀的.此外,分位点为0.9时的回归直线比较平缓,说明在此分位点下,冰川径流随时间的变化并不明显.

观察年径流量随温度变化的回归线(图3)可知,中位数回归直线波动较大,说明在分位点取0.5时,7～8月均温对冰川径流的影响较大,随着气温的升高,径流量显著增加.其他分位点处的回归直线与线性回归直线处于一种相对平衡的状态,径流量随着气温的升高而有所增加.

图3 年径流量随温度变化回归线Figure 3 Regression line of annual runoff variation with temperature

观察年径流量随年降水量变化的回归线(图4)可知,大多数分位点下的回归直线呈上升趋势,分位点较高时的回归直线较为平缓,径流量随年降水量变化并不明显,在分位点是0.95时,径流量呈现出明显的下降趋势,分位点是0.7时,回归直线也呈下降趋势,说明在这两个分位点下,随着年降水量的逐步增加,径流量会呈现出减少趋势.

图4 年径流量随年降水量变化回归线Figure 4 Regression line of annual runoff with annual precipitation

4 结 论

本文使用分位数回归法,分析了天山乌鲁木齐河源1号冰川1989～2017年间径流趋势变化以及其对气候变化的响应关系,并构建了径流与气象因子(气温、降水)之间的分位数回归模型,得出如下结论:

1) 从模型的参数估计及检验可知,在不同分位数下,冰川年径流量呈逐年增加趋势,冰川径流与夏季气温呈明显的正相关关系,与降水的相关关系不稳定.气温和降水对冰川径流的显著性并不一致,在绝大多数分位点下,气温与冰川径流有明显的响应关系,相反,降水对冰川径流的影响并无明显的显著关系,只有在低分位点(0.1)的情形下,气温和降水同时对径流量有显著关系.

2) 1号冰川年径流量序列在较高分位点处的回归直线呈微弱收敛趋势,在中分位点处的回归直线呈明显的收敛趋势,而在较低分位点处的回归直线呈发散趋势.从整体来看,年径流量序列在各个分位点下的回归直线都呈现出上升趋势,表明径流量在逐年增加.

3) 在各个分位点下,径流量随温度的增加而有所增加,尤其在分位点是0.5时,随着气温升高,年径流量序列呈显著增加趋势.年径流量序列在大多数分位点下随降水量的增加有增加的趋势,但在高分位点处反而出现了略微的降低趋势,这可能与冰川区的降水主要以固态形式有关.