1 引言
中国地处东亚季风区， 受地理位置、 地形地貌及气候特征等因素的影响， 气象灾害种类多， 气象灾害频发， 其中干旱与洪涝是我国最主要的气象灾害（马柱国和任小波， 2007； 黄荣辉等， 2012； He et al， 2016a； 庞轶舒等， 2021）， 并且干旱是我国影响范围最广、 造成经济损失最严重的自然灾害。统计资料表明（张庆云等， 2008； 王劲松等， 2012）， 干旱造成农作物受灾面积占农作物受灾总面积的55%， 洪涝占27%， 每年我国三分之二的领土及超过一半的人口都受到各种干旱、 洪涝事件的影响。受全球变暖的影响， 我国洪涝与干旱呈现频发多发之势， 并且干湿转换（旱涝转换/临界转换）越来越频繁（王邵武和叶瑾琳， 1995； 秦大河和Stocker， 2014）。干湿转换事件包括干转湿以及湿转干两种， 是一种典型的降水异常事件。相比于某次干旱或者洪涝事件， 干湿转换事件往往能够导致巨大的灾害（Wu et al， 2006； 王革丽等， 2007； Li et al， 2015）。我国的长江流域、 淮河流域及西南地区干湿转换灾害发生较多， 例如2011年长江中下游地区的旱涝急转给人们的生命财产带来了重大的损失（封国林等， 2012； 沈柏竹等， 2012； 谢五三等， 2021）。
干湿转换往往带来区域性洪涝、 干旱或者水环境恶化等问题（王映思等， 2021）。若在干湿转换前能够进行有效的预警， 政府和相关部分便有足够的时间进行管理规划， 从而尽可能地减小灾害受损情况。因此， 研究我国干湿转换的早期预测预警方法非常重要。
然而干湿转换涉及复杂的非线性过程， 现有的气候模式对可能发生的干湿转换的预测能力非常有限（卞建春和杨培才， 2003； 李建平和丁瑞强， 2009； 万仕全等， 2014）。准确预测干湿转换所引发的区域性旱涝事件是气象部门面临的长期挑战， 因此， 基于现有预测理论和技术水平， 想要准确预测干湿转换仍是一项异常艰巨的任务。已有研究表明， 许多复杂非线性系统具有临界阈值， 即在系统行为或结构发生急剧变化的地方（又可称之为临界点）， 系统将在此迅速地从一种状态转换到另一种状态（Scheffer et al， 2009； Fisher and Scheffer， 2011）。这种临界转换现象广泛存在于自然界和人类社会的各个领域（Venegas et al， 2005； Lenton et al， 2012； Gopalakrishnan et al， 2016； Spielmann et al， 2016）， 并且这种转换往往都是灾难性的， 如气候系统的突变可能造成文明的衰落或朝代的兴亡（Alley et al， 2003； He et al， 2016b）。因此， 亟需开展非线性系统的临界转换的机理及早期预警信号的研究。
幸运的是， 21世纪初， 科学家们发现临界慢化理论在研究复杂动力系统的早期预警信号方面展示了重要潜力（Carpenter and Brook， 2006； Guttal and Jayaprakash， 2008； Scheffer et al， 2009）。临界慢化理论最早被应用于统计物理学， 它是由一种动力系统相态转换形成的分散涨落现象， 这种分散涨落往往表现为扰动的恢复速率变慢、 要素变化幅度增大、 恢复到旧相位的能力变小、 涨落的持续时间拉长等现象（于渌和郝柏林， 1984）。Scheffer et al（2009）通过气候系统、 生态系统、 医学等多个真实系统的研究指出系统向临界点趋近时， 会出现系统扰动后的恢复速率变慢、 方差增大、 自相关系数增大等临界慢化现象， 这为复杂动力系统突变的早期预警信号研究指明了新方向。晏锐等（2011）基于临界慢化理论对水氡浓度进行研究， 结果表明在2008年汶川地震发生前夕出现了方差和自相关系数增大的早期预警信号。近年来， 吴浩等（2012， 2013， 2014）、 Wu et al （2015）、 Tong et al（2014）以及颜鹏程等（2012）将临界慢化理论用于分析我国各区域气温、 太平洋年代际振荡指数以及阿留申低压强度等气候资料， 研究气候突变的早期预警信号。鉴于干湿转换涉及复杂的非线性过程， 可以将这类状态转换看作一种系统的突变， 前人基于临界慢化原理的突变早期预警信号的研究能够为我国干湿转换的早期预警信号的研究提供有效的途径和方法。
鉴于此， 本文拟基于临界慢化理论对我国不同区域干湿转换的早期预警信号进行研究。首先利用滑动t检验（MTT）突变检测方法对我国不同区域干湿转换进行检测， 进而通过表征现临界慢化现象的方差和自相关系数的分析， 探究干湿转换的早期预警信号， 从而为我国旱涝事件的预测提供重要参考， 为防灾减灾工作提供科技支撑。
2 资料来源与方法介绍
2.1 资料来源
所用资料为中国气象信息中心提供的我国753个台站1961 -2020年逐日降水和蒸发资料。经过质量控制和剔除存在缺测5年及以上的站点， 剩余589个站点（图1 ）， 根据研究的需要， 利用589个站的降水和蒸发资料计算了各站K干旱指数， 并采用低通滤波方法提取了K干旱指数的年代际趋势信息。
K干旱指数是降水的相对变率与蒸发的相对变率之比（王劲松等， 2007）， 从该指数的定义上可以看出它不仅考虑了降水在干旱形成中的重要作用， 也考虑到了温度、 风速和太阳辐射等对干旱的影响， 有明确的物理意义。从监测干旱的程度、 范围和持续时间的角度， 来对比分析K干旱指数与其他常用干旱指数的区域适应性， 结果表明， K干旱指数在我国干旱监测应用中有良好的效果， 可以很好地表征干湿变化（王劲松等， 2013； 吴哲红等， 2013； 王素萍等， 2015， 2020）。K干旱指数也是干旱灾害风险评估中表征致灾因子的一个良好指标（郭晓梅， 2018）。为细化干旱监测的时间尺度， Wang et al（2018）又进一步建立了多时间尺度的K干旱指数， 使得K干旱指数在干旱监测和干湿变化的表征方面有了更好的效果（Uttaruk and Laosuwan， 2019）。因此， 选取K干旱指数作为表征干湿变化的指标。
2.2 方法介绍
2.2.1 方差和自相关系数
方差是描述样本中数据相对于均值 偏离程度的特征量， 自相关系数是描述同一变量不同时刻相互之间相关性的特征量（魏凤英， 1999）。本文通过滑动计算序列的方差和自相关系数， 从而寻找K干旱指数序列状态转换的早期预警信号。具体做法概述如下（图2 ）：
图2 滑动计算方差、 自相关系数示意图（a）滑动计算方差， 其中L1， L2， L3， …Ln…代表长度相同的各个窗口（ML）， S1， S2， S3， …Sn…代表对应窗口长度数据的均方差， L为序列总长度， MT为滑动步长； （b） 滑动计算自相关系数， L1（L12）， L2（L22）， L3（L32）， …Ln（Ln2）…代表长度相同的各个窗口（ML）， α1代表L1和L12的自相关系数， α2代表L2和L22的自相关系数， …αn代表Ln和Ln2的自相关系数， LT表示滞后时间， L、 MT意义同（a） Fig.2 We calculated the variance and autocorrelation coefficient by sliding the window. （a） We calculated variance by sliding the window， L1， L2， L3， …， Ln … denote windows of the same length， S1， S2， S3， …Sn … represent the variances of the corresponding windows， L is the total length of the sequence， and MT is the sliding step； （b） we calculated the autocorrelation coefficient by sliding the window， L1（L12）， L2（L22）， L3（L32）， …， Ln（Ln2） … represent windows of the same length， a1 denotes the autocorrelation coefficients of L1 and the L12， a2 denotes the autocorrelation coefficient of the L2， while L22， … the an refer to the autocorrelation coefficients of the Ln and Ln2， LT represents the lag time， and L and MT have the same meaning as those in Fig.2（a） |
滑动计算方差： 以K干旱指数序列为例， K干旱指数序列总长度为L； 从K干旱指数序列中选取子序列， 称之为滑动窗口， 以ML表示； 然后滑动截取相同长度的K干旱指数子序列， 滑动长度称为滑动步长， 以MT表示； 选定ML和MT后， 通过滑动得到多条长度相同的子序列L1， L2， L3， ···Ln···， 计算L1， L2， L3， ···Ln···序列的均方差， 得到一条方差序列S1， S2， S3， ···Sn···。本文在滑动计算方差时， 将ML选为120个月（10年）， MT选为3个月。
滑动计算自相关系数： 同样， 以K干旱指数序列为例， 将子序列L1， L2， L3， ···Ln···往后滑动一个LT（滞后时间）长度得到新的子序列L12， L22， L32， ···Ln2···， 分别计算L1与L12、 L2与L22、 L3与L32、 ···Ln···与···Ln2···的相关系数， 得到一条相关系数序列α1， α2， α3， ···αn···， 因为相关系数计算的子序列都来自K干旱指数总序列， 因此称之为自相关系数。本文在滑动计算自相关系数时， 将ML选为120个月（10年）， MT选为3个月， LT选为1个月。
2.2.2 临界慢化现象与自相关增大、 方差增大的关系
当系统趋近于临界状态、 系统控制参数趋近于临界阈值时， 临界慢化现象表现为序列的自相关和方差系数增大（Carpenter and Brook， 2006； Scheffer et al， 2009）。详细描述可见参考文献（吴浩， 2014）。本文通过计算表征临界慢化特征的方差和自相关系数来寻找复杂动力系统趋于临界点（干湿转换）的早期预警信号。
2.2.3 滑动 t 检验法
滑动t检验法（MTT）是通过考察两组样本平均值的差异是否显著来检验突变的（魏凤英， 1999）。需要注意的是MTT方法子序列时段的选择带有人为性， 本文利用该方法并且多次变动子序列长度来检测不同干旱分区代表站的K干旱指数序列的干湿转换时间。
3 计算结果与讨论
3.1 不同干旱分区的确定
研究表明， 区域干旱化是全球变化的主要特征， 是全球变化区域研究的重要内容。干旱往往表现出明显的区域性特征， 而REOF方法的主要优点在于能更好地反映要素场变率的局地特征， 因此刘晓云等（2012）基于REOF方法选取通过99%显著性检验的值作为分区临界值， 从而保证了各分区之间相互重叠或遗漏的部分尽可能地少， 进而将我国划分为16个干旱异常区域。使用该方法能很好地反映出干旱的区域性特征。本文考虑到其中的12区、 14区和15区为常年干旱区（中国气象局， 2007）， 因此剔除了这3个区域， 将文献中的13区修改为12区、 16区修改为13区， 从而得到了本文的13个干旱异常区（图3 ， 本文仅开展对中国大陆地区的干旱分析， 未对海南省、 台湾省及诸多岛屿进行分析）。为探究13个干旱异常区干湿转换特征， 本文从各个分区中挑选出一个代表站（选取荷载大于0.6， 且1961 -2020年没有数据缺测的站点， 图1 中红色圆点表示）， 通过代表站的K干旱指数研究干旱异常区的干湿转换特征。
图3 我国13个干旱分区图（改绘自刘晓云等， 2012）海南省和台湾省未分析 Fig.3 13 Drought zoning of China （Moditied from Liu et al， 2012） No regional analysis was conducted for Hainan and Taiwan Provinces |
3.2 我国不同干旱分区干湿转换早期预警信号研究
3.2.1 我国不同干旱分区干湿转换检测
图4 近60年陕西榆林站的K干旱指数变化特征及其干湿转换检测图（c）和（d）中虚线表示显著性水平0.05的临界值 Fig.4 Variation characteristics of K drought index and the detection of dry-wet transition at Yulin station in Shaanxi Province in recent 60 years. In Fig.4 （c） and （d）， the dotted lines indicate the critical values at significance level of 0.05 |
表1 不同干旱分区代表站干湿转换时间Table 1 Dry-wet transition time of representative stations in different drought zones |
| 分区 | 站名 | 站号 | 经度/（°E） | 纬度/（°N） | 省份 | 突变年份 |
|---|---|---|---|---|---|---|
| 1 | 榆林 | 53646 | 109.70 | 38.23 | 陕西 | 1985年 |
| 2 | 天门 | 57483 | 113.17 | 30.67 | 湖北 | 2004年 |
| 3 | 赣州 | 57993 | 114.95 | 25.85 | 江西 | 2004年 |
| 4 | 白城 | 50936 | 122.83 | 45.63 | 吉林 | 1995年 |
| 5 | 佛坪 | 57134 | 107.98 | 33.53 | 陕西 | 1981年 |
| 6 | 兴仁 | 57902 | 105.18 | 25.43 | 贵州 | 2002年 |
| 7 | 福海 | 51068 | 87.47 | 47.12 | 新疆 | 1985年 |
| 8 | 甘孜 | 56146 | 100.00 | 31.62 | 四川 | 1989年 |
| 9 | 桦甸 | 54273 | 126.75 | 42.98 | 吉林 | 1990年 |
| 10 | 兖州 | 54916 | 116.85 | 35.57 | 山东 | 1989年 |
| 11 | 佳木斯 | 50873 | 130.28 | 46.82 | 黑龙江 | 2010年 |
| 12 | 吕泗 | 58265 | 121.60 | 32.07 | 江苏 | 2001年 |
| 13 | 都兰 | 52836 | 96.42 | 36.43 | 青海 | 1995年 |
3.2.2 我国不同干旱分区干湿转换早期预警信号-方差信号
根据检测到各区域干湿转换的时间， 本节将利用表征临界慢化现象的方差系数研究干湿转换的早期预警信号。
3.2.3 我国不同干旱分区干湿转换早期预警信号-自相关系数信号
临界慢化理论表明， 系统趋近临界点时会发生自相关系数增大的现象。本节将利用表征临界慢化现象的自相关系数研究干湿转换的早期预警信号。
综上所述， 通过计算序列方差和自相关系数指示的K干旱指数干湿转换早期预警信号的时间基本一致， 进一步验证了基于临界慢化理论研究不同干旱分区干湿转换早期预警信号的可行性。
4 结论
针对干湿转换的复杂非线性特征， 结合非线性科学的相关理论和方法， 本文基于K干旱指数对我国13个干旱区域近60年的干湿转换及早期预警信号进行了研究， 论证了表征临界慢化现象的方差和自相关系数作为突变早期预警信号的可靠性。主要结论如下：
（1） 基于我国589个站的降水及蒸发资料， 通过REOF方法可以将我国划分为13个干旱分区。各干旱分区方差贡献率较大， 满足REOF分区要求。
（2） 通过研究我国13个干旱分区代表站的K干旱指数序列， 检测到各个区域都存在较明显的干湿转换事件， 文章对各个分区特征最明显的干湿转换事件进行研究。
（3） 我国各区域代表站K干旱指数序列干湿转换前10年内都出现了不同程度的方差以及自相关系数增大的临界慢化现象， 表明基于临界慢化理论研究旱涝状态转换早期预警信号的可行性。我国长江中下游、 华南、 新疆北部地区干湿转换早期预警信号时间提前较少（1~3年）， 而我国东北、 华北、 华东及青藏高原等地早期预警信号提前时间较多（7年左右）。
本文的研究有利于深入认识我国干湿转换事件和寻找其干湿转换的早期预警信号， 该方法的广泛应用可以为我国的旱涝预测及防灾减灾提供科技支撑。需要指出的是， 尽管基于K干旱指数的研究表明临界慢化现象导致的动力学的方差增大、 自相关系数增大可能是干湿转换发生前的一个早期预警信号， 但对干湿转换发生前后的环流特征、 外强迫变化特征、 临界慢化现象与状态转换的强度/突变幅度之间的关系等问题仍然需要进一步研究。