四川省不同SPEI指数计算方法适用性评价

康银红 王嘉驰 宋 鑫,2 武剑飞 贺 帅

(1.四川农业大学水利水电学院, 雅安 625014; 2.四川省绵阳市涪城区农业农村局, 绵阳 621000;3.中国农业科学院农田灌溉研究所, 新乡 453002)

0 引言

干旱是一个地区由于长期降水与蒸散收支不平衡导致水分异常短缺而产生的自然灾害,具有发生频率高、影响范围广、持续时间长等特点,不仅会给当地农业等部门带来直接经济损失,也会对水资源、土地资源及社会资源造成间接影响[1]。因此,加强干旱研究有助于减轻干旱造成的危害。目前,干旱研究大多基于干旱指数进行[2-4]。干旱指数能够有效反映地表的干湿情况,是监测、评估干旱的重要参数,也是旱灾研究的基础。因此,开展干旱指数适用性评估,对干旱预警及资源保护具有重大意义[5]。

目前干旱评价指标主要有帕默尔干旱指数(PDSI)[6]、标准化降水指数(SPI)[7]、标准化降水蒸散指数(SPEI)[8]等。其中PDSI指数虽然考虑了温度对干旱的影响,适用全球变暖背景下的干旱研究,但其时间尺度固定,难以表征短期干旱[9]。相较于PDSI指数,SPI指数能够在不同空间和时间尺度上评估干旱,但其仅考虑了降水作用,并未考虑蒸散的影响。SPEI指数综合了前两种指数的优点,不仅可以反映多时间尺度的干湿情况,还通过引入参考作物腾发量(ET0)评估了其他气象因素对干旱的影响,能更加客观地描述地表干湿变化,适用于气候变暖背景下干旱特征的分析,在国内外应用广泛[10-12]。

由于ET0的计算精度会影响SPEI指数的干旱监测能力,因此在利用SPEI监测干旱前,需要考虑ET0的计算精度。Penman-Monteith(PM)法[13]是联合国粮农组织(FAO)推荐使用的ET0计算方法,具有较高的权威性和准确性。由于PM法计算ET0时,需要大量气象参数,对站点的要求较高,所以在某些资料匮乏的地区难以得到有效应用[14-16]。为此,学者们提出了许多ET0简化算法,但这些方法都有特定的适用条件,计算ET0时,需要根据不同地区实际情况选用[17]。目前,针对不同ET0计算方法适用性的研究较多,然而考虑不同ET0算法对SPEI影响的研究还相对较少。因此,有必要厘清不同ET0计算方法的适用性及其对SPEI的影响。

四川省水资源分布差异显著,区域降水不均、地势复杂,容易发生区域性、季节性干旱。近年来,由于受到全球气候持续变暖的影响,四川省干旱强度及频率有增长趋势[18]。特别是2022年7月28日—8月22日,持续高温事件综合强度为1961年有完整气象观测记录以来最强,部分地区旱情严重,农作物减产,群众饮水困难。因此需要选择适宜的SPEI指数对区域降水不均、地势复杂、容易发生干旱的四川省进行准确的干旱监测与评估,并较为准确地研究当地的干旱时空分布情况。本文以PM法计算的ET0为标准,通过误差分析、空间插值等方法比较7种不同方法的ET0计算结果,并选取其中表现较好的3种方法计算得到其对应的SPEI指数,通过理论验证、泰勒图等方法找出适合四川省不同时间尺度和不同区域的SPEI指数,为全球气候变暖背景下四川省的干旱监测和防控提供理论依据。

1 研究区概况

研究区为四川省(26°03′～34°19′N,97°21′～108°33′E),地形地貌复杂多样,东部为盆地,西部为青藏高原,地势西高东低。根据四川省地形特征将研究区分为川西高原、川西南山地和川中盆地(图1)。该区地处亚热带,受地理位置及地形条件影响,区域气候差异显著。川中盆地属于亚热带季风气候,气温高、湿度大;川西南山地为亚热带半湿润气候区,气温高、雨量少;川西高原为典型的高原气候,气候立体变化明显。省内水资源丰富,年平均降水量在1 000 mm以上,但区域水资源分布差异显著,降水多集中在夏季,冬春季雨量稀少。

图1 四川省地形图Fig.1 Topographic map of Sichuan Province

2 数据与研究方法

2.1 数据来源

所需的气象数据主要包括四川省34个气象站点1967—2016年的平均相对湿度、日照时数、日最高气温、日最低气温、日平均气温、日降水量、高度10 m处风速等。数据来自中国气象数据网中国地面气候资料日值数据集(V3.0)(http:∥data.cma.cn),选取的气象站点资料质量控制良好,对资料中部分气象站点的缺失数据通过同时期的临近站点插补处理。四川省历年来的实际旱情数据和相关描述来源于中国气象局国家气候中心(https:∥cmdp.ncc-cma.net/)、《中国气象灾害大典》(四川卷)[19]、《四川水旱灾害》[20]等相关统计资料。

2.2 研究方法

2.2.1ET0计算方法

本文采用8种常用的ET0计算方法,分别为标准综合算法PM法[13];辐射法中的FAO-24Radiation(FAO-Ra)法[21]、Priestley-Taylor(PT)法[22]以及Makkink(MK)法[23];温度法中的Hargreaves-Samani(HS)法[24]和Blaney-Criddle(BC)法[25];物质传输法中的WMO(World Meteorological Organization)法[26]和Rohwer(Ro)法[26]。各计算方法见表1。通过各方法计算得到1967—2016年的日ET0,然后将1—12月的日ET0求和作为当年ET0,各区域的ET0由相应区域内各站点的平均值来表示。

表1 ET0计算方法Tab.1 Calculation methods of reference crop evapotranspiration

2.2.2SPEI计算方法

利用四川省34个气象站点的逐月降水量和ET0计算1、3、12个月时间尺度的SPEI(SPEI1、SPEI3、SPEI12),12个月时间尺度的SPEI用来进行年际分析,其中第12个月的SPEI被用来监测当年的干旱状况。以PM公式为例,所对应的SPEI表示为SPEI_PM。

SPEI计算过程如下[8,27]:

首先计算逐月降水量与参考作物腾发量的差值Di,计算式为

Di=Pi-ET0

(1)

式中i——月数Pi——逐月降水量,mm

ET0——月参考作物腾发量,mm

再采用三参数的log-logistic概率密度函数对Di数据序列进行拟合,计算式为

(2)

式中α——尺度参数

β——形状参数

γ——位置参数

3个参数可以通过线性矩方法拟合获得。通过概率密度函数得到Di数据序列的累计概率函数,计算式为

(3)

对Di数据序列的累计概率函数进行正态标准化处理,并求得SPEI指数。计算式为

P=1-F(x)

(4)

(5)

(6)

式中系数c0、c1、c2、d1、d2、d3分别为2.515 517、0.802 853、0.010 328、1.432 788、0.189 269、0.001 308。

2.2.3精度评价指标

为了探究ET0方法及对应的SPEI在各区域的差异,以PM/SPEI_PM为标准,采用均方根误差(RMSE)、相对误差(RE)评价其他方法的精度,两个指标的值越小,方法的精度越高。为了验证水分亏缺量是否满足log-logistic概率分布,选用Kolmogorov-Smirno(K-S)检验比较经验分布与假设分布的拟合优度,选用纳什效率系数(Nash-Sutcliffe efficiency coefficient,NSE)比较经验分布与假设理论分布的接近程度。NSE越接近于1,表示二者越接近。

3 结果与分析

3.1 ET0计算方法年值

四川省各区域ET0计算结果见图2。由图2可知,在川西高原,PM法多年平均ET0为931.89 mm,最大值为977.82 mm,最小值为882.21 mm。FAO-Ra法、PT法、BC法、HS法均高估了ET0,其中FAO-Ra法偏差较大,多年ET0变化范围为1 067.21～1 195.46 mm,相比PM法偏高166.33～234.92 mm,PT法与PM法表现出较好的吻合度,多年均值偏差为72.92 mm,偏差较小;HS法、BC法均值则分别偏大104.26、71.30 mm;MK法、WMO法均低估了ET0,其中WMO法偏低较多,偏差范围为418.57～483.04 mm,MK法则偏低88.78～151.09 mm;Ro法年际变化不稳定,多年均值比PM法偏低103.66 mm。

图2 四川省各区域ET0计算结果Fig.2 Calculation results of ET0 in each region of Sichuan Province

在川西南山地,PM法多年平均ET0为1 036.46 mm,最大值为1 132.88 mm,最小值为963.02 mm。BC法在川西南山地计算值最高,多年ET0变化范围为1 245.12～1 322.58 mm,相比PM法偏高158.46～322.44 mm,偏差较大;FAO-Ra法、HS法、PT法与PM法年际变化趋势一致,但FAO-Ra法出现了较大的偏差,多年均值比PM法偏高155.01 mm,HS法、PT法多年均值则分别偏高121.71、86.44 mm;MK法、WMO法同样低估了ET0,其中WMO法偏低最多,为392.82～546.74 mm,MK法则偏低 98.25～198.18 mm;Ro法年际变化不稳定,多年均值比PM法偏低75.23 mm。

在川中盆地,各方法计算得到的ET0差异显著。PM法计算得到的多年平均ET0为866.38 mm,最大值为953.31 mm,最小值为795.47 mm。PT法、HS法、BC法均比PM法计算值高,尤其BC法多年ET0变化范围在1 369.22～1 451.76 mm之间,相比PM法偏高469.30～594.48 mm,精度最差。PT法、HS法与PM法年际变化基本一致,多年平均ET0较PM法分别偏高123.73、203.68 mm;辐射法中的FAO-Ra法与PM法计算值最为接近,变化趋势一致,其多年平均ET0为876.75 mm,偏差范围在-19.86～60.15 mm之间,计算结果最优;MK法、WMO法的ET0均比PM法偏低,其中WMO法低估程度明显较大,偏差范围为363.35～504.46 mm,在该地区适用性最差,MK法则偏低155.10 mm;Ro法年际变化不稳定,多年均值比PM法偏低93.39 mm。

3.2 ET0计算方法误差

利用RMSE来评价四川省各简便方法年ET0差异。各方法在不同区域误差变化见图3。从图3来看,在川西高原,PT法和BC法表现出较好的计算精度,RMSE均在76.93 mm以下;HS法、Ro法、MK法的RMSE在102.09～116.62 mm之间,FAO-Ra法的RMSE高达200.59 mm;WMO法计算精度最差,RMSE为448.25 mm。川西南山地PT法和Ro法的RMSE低于99.11 mm,计算结果较好,其次为HS法、FAO-Ra法、MK法,RMSE变化范围在125.13～157.09 mm之间,BC法、WMO法均表现出较大的误差,其中WMO法的RMSE为479.90 mm,精度最差。在川中盆地,FAO-Ra法RMSE最小,仅为18.66 mm,其次为Ro法、PT法、MK法、HS法,RMSE变化范围为123.62～204.74 mm,WMO法和BC法误差较大,RMSE均在454.08 mm以上。

图3 四川省各区域ET0计算方法的RMSEFig.3 RMSE of ET0 calculation methods in each region of Sichuan Province

利用ArcGIS中的克里金插值法对四川省34个气象站点年ET0的RE进行空间插值(图4),分析不同方法在四川省的空间适用性。

图4 四川省各区域ET0的RE空间分布Fig.4 RE spatial distributions of ET0 in each region of Sichuan Province

由图4可以看出,在川西高原,PT法的适用性最好,大部分气象站点的RE介于-3.8%～14.2%之间,只有高原北部的色达、若尔盖及红原一带的RE相对较大,在16.8%～19.3%间变化。而WMO法的RE最差,在整个川西高原地区都低估较多,各站点低估区间为-65.1%～-25.7%,其中北部的若尔盖、红原、色达及石渠地区更是低估-58.4%以上,适用性最差。在川西南山地,PT法的适用性仍然最好,除了昭觉和雷波的RE在15.2%左右外,其余站点RE仅在-1.3%～13.3%变化;而WMO法在整个川西南地区都低估较多,各站点相对低估-64.5%～-28.7%,计算精度最差。在川中盆地,FAO-Ra法的适用性最佳,除了巴中的RE高于6.4%外,其余各站点的RE都相对较小,介于-2.2%～3.3%之间,而WMO法在整个川中盆地都低估较多,除了万源和雅安外,其余各站点均低估-62.6%～-48.1%,评价效果最差。

由于四川省不同地区存在气候差异,导致各方法在不同地区应用时差异较大。HS法在纬度较低、平均温度高、辐射强的川西南山地ET0最大,这与HS法主要考虑平均温度和昼夜温差影响,同时利用太阳天顶辐射进行计算有关。BC法随着区域海拔降低,其数值逐渐偏大,这与区域海拔降低,气温逐渐升高有关。PT法忽略了空气动力学项,仅考虑了辐射项,在四川省3个区域均高估了ET0,这与符娜等[28]对云南地区进行研究得出PT法高估ET0的结论一致。由于忽略了土壤热量,MK法在四川省3个区域均低估了ET0。RAHIMIKHOOB等[29]对伊朗北部8个气象站点研究也发现,MK法在所有站点均低估了ET0,左德鹏等[30]、李志[31]也得出类似结论。在湿润的川中盆地,FAO-Ra法的ET0计算值最接近PM法,这是因为FAO-Ra法在湿润区的计算效果较优。NANDAGIRI等[32]评价了印度4个不同气候区的气象站点逐日ET0,发现FAO-Ra法在湿润区的计算效果优于干旱区。樊军等[33]的研究也表明FAO-Ra法在半干旱地区估算ET0较差。Ro法与WMO法同属于物质传输法,其计算精度对使用的参数变化非常敏感[26]。尽管Ro、WMO法在3个区域计算结果一致偏小,但Ro法偏差较小,WMO明显偏差较大,这说明Ro法的初始参数在四川省具有较好的适用性,而WMO法的参数需要根据天气情况进行调整。

从各项误差指标来看,MK、Ro法在四川省3个区域的表现较为稳定,PT法在川西南山地和川西高原均表现出较低的误差,WMO法在各区域适用性最差,综合考虑,选取PT、MK、Ro 3种简便方法计算的ET0进行SPEI分析。

3.3 基于不同ET0计算方法的多尺度SPEI适应性

3.3.1SPEI指数理论验证

计算SPEI指数的前提是降水量和蒸散量之差(Di)服从三参数对数逻辑分布。因此,在计算指数之前,本文通过K-S检验和NSE检验探究四川省各气象站点的Di是否服从以上概率分布,结果如表2所示。

表2 Di分布检验结果Tab.2 Distribution test results of Di

由表2可知,采用三参数的对数逻辑分布函数对Di进行拟合,通过K-S检验和NSE评估拟合效果,结果发现34个站点的Di均通过K-S检验,且NSE基本都在0.95以上,表明三参数对数逻辑分布函数在四川省各月份不同时间尺度下的Di拟合效果都较好,可以用于SPEI指数的多时间尺度计算。

3.3.2时间变化

对四川省不同区域1967—2016年SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro的年值进行比较(图5),其中红色图例表示实际旱情。

图5 四川省1967—2016年12个月时间尺度SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro年际变化曲线Fig.5 Temporal variations of SPEI_PM, SPEI_PT, SPEI_MK and SPEI_Ro at 12-month timescale in Sichuan Province during 1967—2016

由图5可知,在四川省的3个区域中,SPEI_PT、SPEI_MK、SPEI_Ro具有与SPEI_PM相似的变化情况,其变化值均处于-2～2之间。从四川省各区域的SPEI计算结果来看,基于不同ET0方法计算的年SPEI差异较小,变化趋势基本一致,4种方法均显示1967—2016年,川西高原、川西南山地的SPEI有增加的趋势,川中盆地则存在SPEI减小的趋势。虽然SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro最小值不同,但时间分布一致。在有实际旱情的年份,最小值均低于0,因此4种SPEI指数都能够识别历史干旱事件。总的来说,在同一区域,不同方法之间会有差异,但SPEI几乎都维持在同一水平,干湿情况交替出现。

3.3.3误差分析

图6为四川省1967—2016年间SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro的泰勒图,其中点A为观测点,代表SPEI_PM。由图6可知,在四川省各区域,12个月时间尺度下的SPEI_PT、SPEI_MK、SPEI_Ro均与SPEI_PM有较好的相关系数,这是因为SPEI是通过历史同期数据的异常来监测干旱,主要与ET0变化情况有关,而受ET0本身的影响较小。从标准偏差来看,SPEI_PT、SPEI_MK、SPEI_Ro在川西高原、川西南山地差别不大,都较接近1,而在川中盆地,SPEI_Ro的标准偏差明显大于SPEI_PT、SPEI_MK,达到1.14。从RMSE来看,SPEI_MK与SPEI_PM的RMSE最小,其次依次为SPEI_PT、SPEI_Ro。综合来看,在四川省全域,SPEI_MK与SPEI_PM有更好的相关性。

图6 四川省12个月时间尺度SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro的泰勒图Fig.6 Taylor diagrams of SPEI_PM, SPEI_PT, SPEI_MK and SPEI_Ro at 12-month timescale in Sichuan Province

相较于川西高原和川西南山地,川中盆地各方法计算的SPEI更具敏感性。对川中盆地不同时间尺度下SPEI_PT、SPEI_MK、SPEI_Ro的误差进行分析,如图7所示。由图7可知,1、3、12个月时间尺度下,SPEI_PT、SPEI_MK与SPEI_PM都具有较好的相关性,而SPEI_Ro与SPEI_PM的相关性最差。从RMSE来看,1、3、12个月时间尺度SPEI_PT、SPEI_MK与SPEI_PM的RMSE大都集中在0.05～0.20区间内,而12个月时间尺度SPEI_Ro与SPEI_PM的RMSE达到0.44。综合来看,大部分情况下,SPEI1相关性比SPEI3、SPEI12好。1个月时间尺度下,SPEI_MK与SPEI_PM有更好的相关性,相关系数达到0.99,RMSE仅为0.15,SPEI_Ro与SPEI_PM相关性最差。

图7 川中盆地1、3、12个月时间尺度SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro的泰勒图Fig.7 Taylor diagrams of SPEI_PM, SPEI_PT, SPEI_MK and SPEI_Ro at 1-, 3- and 12-month timescales in central Sichuan basin

3.4 周期分析

由于川中盆地各方法计算的SPEI更具敏感性,因此对川中盆地不同时间尺度下SPEI_PT、SPEI_MK、SPEI_Ro的小波实部及方差进行分析,其对比结果如图8、9所示。图8中等值线的正负代表偏湿润和偏干旱时期,图9中小波方差反映SPEI时间序列的能量波动随时间尺度的变化特征,峰值对应的时间尺度为SPEI变化的主要周期。

图8 川中盆地1、3、12个月时间尺度SPEI小波分析Fig.8 Wavelet analysis of SPEI at 1-, 3- and 12-month timescales in central Sichuan basin

图9 川中盆地1、3、12个月时间尺度SPEI小波方差Fig.9 Wavelet variances of SPEI at 1-, 3- and 12-month timescales in central Sichuan basin

由图8可以看出,在时间尺度10～15 a上,各方法的SPEI1周期振荡变化明显,贯穿整个研究时段,经历了5次“偏干旱-偏湿润”的交替变化,2016年等值线小于0且向负值中心变化,说明未来一段时间将处于偏干旱时期。在时间尺度5～10 a、15～20 a上,各方法的SPEI3周期振荡变化明显。结合图9发现,SPEI3的小波方差峰值出现在9 a的3.5和 19 a的5.5,分别对应第2主周期、第1主周期。在时间尺度10～15 a上,各方法的SPEI12周期振荡变化明显,经历了5次“偏干旱-偏湿润”的交替变化。在1个月时间尺度上,SPEI_PT表现出与SPEI_PM较好的一致性,两者第1主周期都位于13 a,小波方差差距最小,SPEI_Ro则表现出与SPEI_PM较为不同的周期性,其第1主周期位于11 a,与SPEI_PM差距最大。在3、12个月时间尺度上,各方法的SPEI都具有相似的第1主周期。总体而言,SPEI_PT和SPEI_PM具有最相似的周期振荡变化,SPEI_Ro和SPEI_PM的周期差距最大。

4 结论

(1)PT法等7种ET0算法的空间变异性十分明显。在川西高原和川西南山地,PT法与PM法有较好的相关性,适用性较强;在川中盆地,辐射法中的FAO-Ra法适用性最好。

(2)结合历史干旱事件分析,SPEI_PM、SPEI_PT、SPEI_MK、SPEI_Ro 4种指数都能够识别出四川省大多数干旱事件,具有较好的干旱监测能力,其中SPEI_MK监测的干湿情况最接近SPEI_PM。

(3)大多数情况下,SPEI1相关性比SPEI3、SPEI12好。1个月时间尺度下,SPEI_MK与SPEI_PM有更好的相关性,而SPEI_Ro与SPEI_PM相关性最差。

(4)SPEI_PT和SPEI_PM具有最相似的周期振荡变化,SPEI_Ro和SPEI_PM的周期差距最大。

(5)四川省SPEI_MK监测的干湿情况最接近SPEI_PM,数据缺失的条件下,SPEI_MK可以作为SPEI_PM的替代指数,这为四川省的干旱监测和防控提供理论依据,能够减少和预防干旱带来的灾害。由于本文仅研究了ET0计算方法及相应SPEI指数在四川省的适用性,下一步可深入探究不同SPEI指数监测的气象干旱是否会在干旱传播过程中对水文干旱、农业干旱造成影响。