（1.河南农业大学理学院,河南 郑州 450002；2.河南农业大学信息与管理科学学院,河南 郑州 450046；3.河南农业大学农学院,河南 郑州 450046）

　　为探究太阳黑子活动与地区降雨量的关联性,采用连续小波变换方法,分析了1980年以来郑州地区(新郑站点)降雨量和太阳黑子数的数据变化,并对不同时段两者的相关性进行研究;进一步对两者进行交叉小波变换和小波相干谱分析,并根据两者关系通过太阳黑子数观测值对降雨量进行预测。结果表明:(1)不同时段降雨量和太阳黑子数的相关性存在正负差异的现象。降雨量的第一主周期尺度是21 a,在此主周期尺度下得到14 a主周期;太阳黑子数的第一主周期尺度是16 a,在此主周期尺度下得到11 a主周期(与经验值相符)。降雨量与太阳黑子数的主周期相差了3 a,因此导致两者的相关性在不同时段存在正负差异。(2)降雨量与太阳黑子数在1992—2008年的8~12 a时间尺度上关联性显著,且降雨量比太阳黑子数存在规律性的时间滞后,两者在2~4 a和7~10 a的时间尺度上关联性较好;降雨量比太阳黑子数分别在1991—2004年和2006—2013年这两个时间段上呈现规律性的滞后,两者在其他时间段的各时间尺度上关联性不明显。(3)根据延迟年数经验公式,由太阳黑子数观测值对降雨量进行预测,最近的降雨量峰年在2022年附近,与2021年郑州“7·20”特大暴雨的发生时间相近,而下一个降雨量极小年在2028年前后,且随着时间的推移极端旱涝情况可能会愈加严重。

　　河南省郑州市地处平原,农业作为其经济支柱,对降雨量的敏感度较高。2021年7月,河南省遭遇极端强降雨天气,称为“7·20”特大暴雨[1];7月17日8时至23日8时,郑州市累计降雨量400 mm以上面积达到5590 km2,600 mm以上面积达到2068 km2,降雨过程在郑州市气象观测记录历史上属于强度最强[2];7月20日郑州国家气象站出现624.1 mm最大日降雨量,与郑州平均年降雨量(640.8 mm)接近,是建站以来最大值(189.4 mm,1978年7月2日)的3.4倍,特别是7月20日16时至17时,郑州小时雨强达到201.9 mm,突破我国大陆气象观测记录历史极值(198.5 mm,1975年8月5日河南林庄)[3]。在地球视角下这次强降雨的成因,是在西太平洋副热带高压异常偏北、夏季风偏强等气候背景下,同期形成的2个台风汇聚输送的海上水汽与河南上空对流系统叠加,遇伏牛山、太行山地形抬升而形成的[4]。

　　然而,探究此事件是偶发性事件还是趋势性事件,并以此预测未来的极端强降雨情况,需要对降雨量进行时间序列分析[5]。从长期趋势探究降雨量的影响因素,太阳的影响作用不可忽略,太阳作为地球重要的能量来源,对大气温度、河流的径流量、地区的平均降雨量等都有影响[6]。太阳黑子数是反映太阳活动的重要数据指标,常作为评估太阳活动的代表被人们研究[7]。在太阳黑子活动周期与地球自然条件(气温、降雨量等)变化周期的相关性研究方面,赵新华等基于太阳黑子等数据运用小波变换和交叉相关分析等方法对太阳活动与地表温度的关系进行研究,发现太阳活动在百年的时间尺度上对地表温度的变化具有不可忽略的影响[8]。龚兰兰等运用一元线性回归分析、Mann-kendall非参数检验法以及Morlet小波变换及交叉小波变换等方法分析太阳黑子活动与安康市气温变化的关系,认为安康市近60年来平均气温变化周期与太阳活动相关[9]。在降雨量方面,郭天驰等运用互相关谱分析与小波变换法研究沈阳地区降雨量与太阳黑子活动,得出4个降水周期近似等于1个太阳黑子活动周期的认识[10]。王涛等基于小波变换方法分析了天山中西部降水与太阳黑子相对数变化序列,发现两者在准10 a尺度上关联性显著[11]。在数据处理方法上,采用连续小波变换(continuous wavelet transform, CWT)[12]与交叉小波变换(cross wavelet transform, XWT)[13]方法对数据进行处理,可以形成CWT谱与XWT谱,用于判断相应物理量扰动的高低能量区域的分布,虽然不能精确地分析到具体时间,但是该图谱可以图像化展示数据的规律性分布[14]。

　　本文在连续小波变换的基础上运用交叉小波变换和小波相干性来综合分析1980年以来郑州地区降雨量与太阳黑子数之间的关系,以期获得更加可靠的周期性规律与周期对应的时间分布情况。本文关于太阳黑子活动与郑州地区降雨量对应关系的研究有望对于预防和降低郑州地区因旱涝灾害所形成的巨大农业生产损失起到指导和帮助作用。

　　河南郑州地处华中地区,位于黄河下游,是中原核心地带,属于典型的北温带大陆性季风气候。郑州市全年平均气温15.6 ℃,8月份最热,月平均气温25.9 ℃,1月份最冷,月平均气温2.15 ℃。全年平均降雨量542.15 mm,夏天的降雨量约占全年的55%[15]。郑州市春季较短、沙尘较大且比较干燥,夏季降雨充沛,秋季日照时间较长,冬长寒冷少雨雪[16]。因为华中地区的年降雨量主要集中在每年的6至8月份[17],所以本文基于河南省郑州市新郑气象站站点数据,分析了1980—2021年6—8月份降雨量数据,构建了郑州地区42 a年平均降雨量的变化趋势。

　　一个太阳活动周期约为11 a,一个太阳活动周期开始和结束的标准是:当月平滑太阳黑子相对数达到极小值时,旧太阳周结束、新太阳周开始。太阳活动周期以1755—1766作为第一个周期开始计算,每个太阳周的长度不尽相同[18]。太阳黑子年平均数来源于比利时皇家天文台,网址为 。本文主要分析了1980—2021年期间的太阳黑子年平均数。

　　小波变换适用于解决地球物理时间序列的复杂周期问题。对郑州地区在大气层气流、经纬度、海拔高度、地形变化等影响下的多年降雨量数据,按照时间先后顺序排列成数列,小波变换可清晰地显示出潜藏在此时间序列中各种各样的变化周期,展现不同时间尺度下变化的趋势,因此可根据变换趋势对未来发展做预测。小波变换使用一簇小波函数系来趋近目标信号或函数,展现出正负交替的振荡波形式,其函数值在正负两部分的能量相等,因此,积分值等于零,用公式可以表示为[19]:

　　小波函数有很多种,在应用于时间序列分析时,需要得到平滑连续的小波振幅和位相信息。由于复值小波有虚部,可以较好地表达位相信息,通常选择复值小波进行分析。在复值小波中,Morlet小波在时频空间的局部性较好,因此常被用于地球物理时间序列分析。Morlet小波解析形式为[20]:

　　从Morlet小波变换系数的实部可以得到各种特征时间尺度在时间段中的分布和位相信息,而从其模中可以得到特征时间尺度信号的强弱信息。在应用于降雨量时间序列分析时,Morlet小波变换系数的实部可以表示不同周期尺度在考察时间段中的分布,以及降雨量峰谷值信息。离开时间尺度谈周期是没有意义的。Morlet小波变换系数的模可以表示不同时间尺度下变化周期性强弱在时间域中的分布。

　　小波方差跟随尺度因子a的变化而形成的图像称为“小波方差图”。小波方差在图像中能够显现出信号变化随尺度a的能量分布情况。所以,小波方差图能够用来测定信号中不同种时间尺度受到干扰的波动的相对强度,在时间序列中可以表示周期波动的强弱,分析其主要趋势下的周期,即为主周期。

　　以上连续小波变换可以分析单个时间序列的周期强度等信息,在连续小波变换的基础上,通过交叉小波变换可以分析两个时间序列在多个时间尺度方面的相互关系,揭示它们在时频空间中的共同功率和相对位相。两个时间序列x和y的交叉小波谱定义为[20]:

　　因小波在时间域未完全局部化,使得连续小波变换存在边缘伪影。为此,引入一个影响锥(Cone of Influence,COI),本文将影响锥设为由边缘不连续小波功率下降到边缘值的e2的区域。

　　交叉小波变换适用于突出两时间序列间耦合性比较显著的区域,对于两序列间交叉小波功率谱一致性不高的区域,可以利用小波相干谱来表征两序列间变化的一致程度.

　　小波位相角箭头方向可判定两时间序列在相应的时间、尺度下的位相关系。当两时间序列的交叉小波相关系数值较高时,可以通过小波位相角箭头方向与时间轴夹角估算后一序列相对前一序列的提前(滞后)程度。由此,本文通过交叉小波变换位相角来分析降雨量序列相对于太阳黑子数序列在某一时间段上的提前(滞后)程度。

　　对两个时间序列进行连续小波变换和交叉小波变换,可以揭示它们相同的高能区域和位相关系。而小波相干谱则用来度量时频空间中两个时间序列局部相关的密切程度,即使对于交叉小波功率谱中低能量值区,两者在小波相干谱中的相关性也可能显著。定义时间序列x和y的小波相关谱为[23]:

　　其中，c1和c2是标准化常数,П是矩形函数,0.6是由经验确定的Morlet小波的尺度去相关长度因子[20]。

　　交叉小波变换和小波相干谱都可以得到两个时间序列之间位相信息,但交叉小波变换主要在时频域中展现两序列周期性强度一致的区域,能直观了解两序列间共有周期的强度,而小波相干谱侧重于展现两个序列周期性变化趋势的一致性,并不直接展现变化周期的强度关系。

　　1980—2021年郑州年平均降雨量与太阳黑子数统计结果如图1所示,在1990和2000年处,平均降雨量数据与太阳黑子数同时出现相对高值;在1986、1997、2008和2018年处,平均降雨量与太阳黑子数同时出现相对低值;在1991、2001和2014年处平均降雨量出现相对低值,而太阳黑子数则出现相对高值;在2009和2019年处,平均降雨量出现相对高值,而太阳黑子数则出现相对低值。以上初步对比表明,太阳黑子数与郑州降雨量的周期可能存在某种规律。

　　对郑州降雨量数据进行连续小波变换处理,由于郑州降雨量数据是有限时间序列,在时间序列的两端可能会产生“边界效应”。为消除或减小序列两端点附近的边界效应,需要先对两端数据进行对称延伸,在进行完小波变换后,再去掉两端延伸数据的小波变换系数,保留原数据序列时段内的小波系数。

　　利用MATLAB小波工具箱中的信号延伸(Signal Extension)功能,对郑州降雨量数据两端进行对称延伸,由于1980—2021年是42 a的数据,分别向前对称延伸21个单位,向后对称延伸21个单位;然后计算其复小波系数,删除延伸数据,计算Mrolet复小波系数的实部、模、模方和方差。下文将逐一进行绘图分析。

　　对复小波系数的实部数据,以年份为横坐标,时间尺度为纵坐标,绘制小波系数实部等值线所示。用颜色冷暖展现复小波系数的实部值,颜色较暖的区域代表降雨量较高的时期,颜色较冷的区域代表降雨量较低的时期。复小波系数实部等值线图展现了不同时间尺度下降雨量变化周期性强弱随时间分布。在图2中每个等值线圈都代表降雨量数据的一次振荡,颜色越深代表振荡越剧烈。两个颜色相反的等值线圈连在一起就是一个振荡周期(从右边区域能量密度条上也可看出数值从正值最大到负值最大)。等值线圈中心所在纵坐标表示了在此时间尺度下周期性较明显。从线圈中心横坐标能看出大致周期。由图2可以看出,郑州降雨量的变化存在5~8 a、16~22 a、29~35 a的周期变化特征。在29~35 a的时间尺度上存在2个完整的振荡周期;在16~22 a的时间尺度上存在3个完整的振荡周期;在5~8 a的时间尺度上振荡周期的次数较多且较复杂,1980—1996年期间的周期性不如1996—2020年期间的明显。

　　对前面所得郑州降雨量数据的复小波系数模值作图,如图3所示。同样用颜色冷暖表示复小波系数模值大小,其含义为不同时间尺度下变化周期性强弱在时间域中的分布,系数模值越大(对应颜色条上越暖),表明其所对应时段和尺度的周期性越强。从图3中可以看出,在1980—2005年期间,18~23 a的周期性比较明显;在2009—2021年期间,22~29 a的周期性最为显著;此外,在1997—2003年、1995—2005年期间也存在一些相对不明显的周期性。总体上在全时间段上都有一定的周期性,在2005年以后的19~35 a的时间尺度上的周期比较明显,此结果与小波系数实部等值线)相互呼应。

　　对前面所得郑州降雨量数据的复小波系数模的方差进行制图,得到小波方差曲线所示,横轴为周期,纵轴为方差。小波方差图可以用来确定时间序列变化过程中存在的主周期尺度,主周期尺度有多个,小波方差图中每个峰值都是一个主周期尺度。从图4上可知,郑州降雨量变化中有准8 a、21 a以及30 a的主周期尺度,分别为第三、第一和第二主周期尺度(方差值最高的极值点是第一主周期尺度,以此类推)。

　　从图4可以发现,第一主周期尺度准21 a的方差最大,代表此时间尺度上周期性最强,绘制降雨量21 a尺度主周期趋势图。取出前面所得郑州降雨量数据的复小波系数实部数据中纵坐标为21对应的那一行数据(图2)进行画图,得到21 a尺度的主周期趋势图,如图5所示,横坐标为时间,纵坐标为降雨量数据的复小波系数实部值。由图5可知,在21 a尺度下,降雨量的主周期为14 a。

　　与郑州降雨量数据处理方法一样,对太阳黑子数据进行连续小波变换处理,绘制如图6所示的太阳黑子数小波系数实部等值线,图中颜色较暖的区域代表太阳黑子活动较为剧烈的时期,颜色较冷的区域代表太阳黑子活动较为平缓的时期。由图6可知,太阳黑子数的小波系数实部图在拥有本文分析数据的时间段上存在14~18 a的周期变化特征,时间尺度表示变化周期范围在14~18 a(所有等值线的中心所在纵坐标)的时间尺度上存在约3.5个震荡周期(一个完整的等值线圈是半个振荡周期,相邻冷暖的两个等值线圈是一个振荡周期)。

　　同样,对太阳黑子数据的复小波系数的模值制图,如图7所示,图中用颜色冷暖展现复小波系数模值的大小,用于表征不同时间尺度变化周期性强弱在时间域中的分布,系数模值越大(对应颜色越暖),表明其所对应时段和尺度的周期性越强。由图7可知,在2002—2021年期间,18~24 a的周期比较明显;在1980—2002年期间,14~18 a的周期能量值较高,周期性最显著。总体来说,在分析数据的全时间段上都有一定的周期性,在2002年以后的14~18 a的时间尺度上的周期比较明显。以上结果也与小波系数实部等值线   太阳黑子数的小波系数模等值线

　　将所得的太阳黑子数据绘制如图8所示的小波方差图,小波方差图可以用来确定时间序列变化过程中存在的主周期尺度,即峰值所在横坐标。由图8可知,太阳黑子数据变化中仅有16 a的周期变化,是在1980—2021年时间段内唯一的主周期尺度,即第一主周期尺度。

　　绘制太阳黑子数在16 a尺度的主周期趋势图。取出前面所得太阳黑子数据的复小波系数实部数据中纵坐标为16对应的该行数据(图6)进行画图,如图9所示,可知在16 a尺度下,太阳黑子数的主周期为11 a。

　　由图5可知,降雨量以14 a为主周期进行变化。由图9可知,太阳黑子数以11 a为主周期进行变化。因图5和图9都计算了复小波系数实部,等价于进行了数据归一化处理,可将两图合并进行对比分析,如图10所示,横轴是年份,纵轴是两数据的复小波系数实部值。再将图10与图1进行对比,可以发现在去除不必要的扰动基础上,两者保持相同的变化趋势。

　　由于降雨量与太阳黑子数两个序列的主周期不同,相差3 a,所以从图10中可看到,降雨量与太阳黑子数的几个周期的波峰所对应时间也相应延后,1980—1981年两者的波峰几乎处于同一时期,下一个周期中,降雨量的波峰(1993年)相对太阳黑子数的波峰(1991年)延迟了2 a,第三个周期中,降雨量的波峰(2007年)相对太阳黑子数的波峰(2002年)延迟了5 a。由此可知,降雨量与太阳黑子数的波峰(或波谷)在延迟年数上随着时间的推移越来越大。对波峰延迟年数建立数列模型,波峰延迟年数列的通项经验公式为:

　　其中,n=1表示从1980年开始的第一个峰值周期,n=2表示1991—1993的第二个峰值周期,以此类推。

　　可以由太阳黑子数下一个峰值时间,通过波峰延迟年数模型推测出下一个降雨量峰值的时间段。然而,对比图1可知,精确到1 a的降雨量数据可能会由于复杂因素影响而有波动,但不影响对相邻区间范围内整体趋势的预测。例如,根据模型可知,在第4个峰值周期中,两者峰值延迟为9 a,由于太阳黑子数峰值在2013年,由此可计算此周期中降雨量峰值应在2022年附近,与2021年郑州“7·20”特大暴雨的发生时间接近。

　　类似的,分析降雨量与太阳黑子数的几个周期的波谷所对应时间,可以看出在1986和1987年两者的波谷相差1 a;在下一个周期中,降雨量的波谷(2000年)相对太阳黑子数的波峰(1996年)延迟了4 a;在第三个周期中,降雨量的波谷(2007年)相对太阳黑子数的波峰(2014年)延迟了7 a。由此可知,降雨量与太阳黑子数的波谷在延迟年数上随着时间的推移越来越大。对波谷延迟年数建立数列模型,波谷延迟年数列的通项公式为:

　　其中,n=1表示从1986年开始的第一个谷值周期,n=2表示1996—2000的第二个谷值周期,以此类推。

　　根据模型可知,在第4个谷值周期中,两者谷值延迟为10 a,由于太阳黑子数谷值在2018年,由此可计算下一个干旱期可能在2028年前后。由此,即可解释两个时间序列既有正相关又有负相关的现象。

　　另外,从图10中纵向对比两条曲线可以看出,太阳活动的强度随时间越来越弱,而降雨量变化幅度随时间越来越大。所以从整体趋势上来看,今后出现极端气候的情形可能愈加严重。