1 引言
在全球变暖背景下， 天气气候灾害的突发性、 频发性与持续性呈加剧趋势， 其成因亦变得日趋复杂（丁一汇， 2013）。青藏高原是中国主要的牧业区， 地广人稀， 生态环境十分脆弱， 雪灾是高原最主要、 影响最广、 破坏力最大的气象灾害， 长期以来严重制约着当地牧业生产的发展和民族地区经济的振兴（高懋芳等， 2011）。雪灾发生时， 由于降雪量大， 地面积雪深， 加之气温低， 雪面冻结， 往往造成大批牲畜死亡， 使牧区财产遭受重大损失， 并对其生命安全构成威胁（梁潇云等， 2002）， 如1985年10月7-20日青南牧区25万平方千米的地区发生雪灾， 唐古拉一带积雪厚度达50~100 cm， 气温骤降至-42~-24 ℃， 雪灾造成193万（只）牲畜死亡， 使3000多户牧民绝畜， 经济损失1.2亿元（温克刚等， 2007）； 又如1997年9月至1998年4月西藏6地市普降暴雪， 使全区40个县、 457个乡（镇）、 3104个村， 904779人受灾， 受灾牲畜占全区牲畜存栏总数近一半， 其中那曲地区因灾致2.6万人返贫， 因病及交通事故死亡63人， 2万多人患雪盲， 直接经济损失达6亿元（温克刚等， 2008）。因此， 加强高原雪灾机理研究， 对藏区的生产生活和生态安全屏障的建设有着重要的意义。
青藏高原积雪异常对周边区域及东亚天气气候的影响机理一直是科学家关注的焦点， 众多学者也做了大量研究（Yanai et al， 1992； 范广洲等， 1997； 郑益群等， 2000； Qian et al， 2003； 于乐江等， 2008； Yu et al， 2016）， 但对于青藏高原积雪本身异常变化原因的研究相对较少。韦志刚（2001）研究表明， 高原冬春总体积雪的异常主要是由11月至翌年3月的积雪异常决定的， 其中高原积雪负异常至少持续4个月， 正异常一般为2~4个月。积雪异常包含了许多因子的作用信息， 马丽娟（2008）指出青藏高原积雪厚度与北大西洋涛动（NAO）有显著正相关， NAO 偏强时高原积雪普遍偏多， 影响的显著区域在高原中部， 由西南向东北延伸呈带状分布。梁潇云等（2002）认为， 来自孟加拉湾地区的水汽为高原东部牧区强降雪天气过程提供条件， 主要降雪区与水汽通量辐合区吻合。很多海域海温与青藏高原积雪异常亦有着密切的关系。Shaman et al（2005）通过大气驻波遥相关机制提出ENSO事件在对流层激发出的准静止大气罗斯贝波， 引起冬季高原上的风暴活动增多、 对流活动增强、 降雪增加。强El Niño年， 冬春季节南支西风扰动的活跃， 将引起高原以及华南的多雨雪； 张涛涛（2013）分析显示大雪灾和特大雪灾多发生在厄尔尼诺年份和正常年份， 拉尼娜年份相对较少。徐小玉等（2016）认为赤道中东太平洋的厄尔尼诺事件是通过沿200 hPa副热带西风急流传播的一列准静止罗斯贝波对高原地区对流活动及雪深产生影响。Yuan et al（2012）也研究发现印度洋偶极子IOD（Indian Ocean Dipole）与秋季青藏高原积雪呈显著正相关。刘华强等（2006）通过RegCM模式并引入次网格尺度地形重力波拖曳模拟得出， 在排除ENSO的影响后， 印度洋偶极子IOD依然与秋季（初冬）青藏高原积雪相关显著。可以看出， 海温也是积雪异常的主要因子之一， 积雪通过“大气桥”受到海洋的强迫和调控作用。然而以往研究均局限于单一海区的影响， 并没有把多个海区放到同一平面上， 分离它们各自的作用。因此， 有必要在青藏高原敏感且强烈气候变化的背景下， 进一步对影响高原雪灾变化的海温因子进行细致地讨论。
本研究基于青藏高原气象台站逐日观测积雪资料及Hadley中心海温月平均资料， 利用广义平衡反馈分析与主成分分析（GEFA-EOF）相结合的最优反馈模分析方法， 探讨青藏高原雪灾频数对关键区海温异常的响应特征， 综合分析不同海域海温异常对雪灾的独立贡献， 以期为高原雪灾预测提供一定理论基础。
2 资料来源和方法介绍
2.1 资料来源
由于高原地区尤其是西藏地区的气象台站建站较晚， 考虑到资料的稳定性， 本文选取了青海、 西藏1978 -2014年冬半年（当年10月至次年3月）72个气象台站逐日积雪观测资料， 海温数据采用了分辨率为 1°×1°的Hadley 中心海温月平均数据， 大气环流资料为500 hPa高度场、 600 hPa湿度场和风场NCEP/NCAR 月平均再分析数据。
根据青海省气象灾害地方标准（DB63/T372-2011）的权威定义和西藏气候中心对雪灾的研究成果（黄晓清等， 2013）， 当积雪深度D≥2 cm且持续日数T>10天， 或积雪深度D≥5 cm且持续日数T≥5 天即达到雪灾标准， 该标准符合当地实际并长期在气候业务中应用。本文基于此标准进行雪灾统计， 目的是采用时间尺度来探究单站雪灾过程的次数， 不考虑强度分级。当单站积雪状态达到以上两个标准的任何一个， 记为一次雪灾， 即研究中基于单站空间尺度和日时间尺度来计算雪灾发生的频数。对于跨月的雪灾过程， 规定以某站某月内有一次雪灾过程的起始日为准， 雪灾频数计入本月。
2.2 广义平衡反馈分析方法（GEFA）
GEFA方法是近年来在研究海气耦合作用的过程中发展起来的， 用于将海洋对大气的反馈作用与大气对海洋的强迫作用分离， 是观测资料分析和动力模式分析联结在一起的有效统计工具（Wen et al， 2010）。该方法通过计算气候时间尺度上区域大气对下垫面异常的响应系数， 来反映该下垫面异常对大气的独自贡献。具体计算原理为：
设某地在气候时间尺度上大气异常变化为xt， 则：
根据GEFA原理， x（t）可用两部分线性表示， 一部分来自下垫面海温异常Yt的反馈作用； 另一部分是大气内变化 Nt， 即：
其中：
式中： bij为区域i大气对第j个下垫面异常的响应系数， 反映了该下垫面异常对大气的独自贡献。考虑到气候时间尺度上大气内变化可以看作白噪音， 后期大气无法影响前期下垫面异常， 即：
由方程（3）可得反馈系数阵：
式中： Cxy（τ）是 x（t）与 y（t-τ）的协方差阵； Cyy（τ）为y（t）与y（t-τ）的自协方差阵； τ为大气对下垫面的响应时间， 文中取τ=1（江志红等， 2013）。
B的估计精度、 稳定性与样本长度密切相关， 在给定样本情况下， GEFA会产生样本误差， 魏红成等（2013）研究表明， 当样本数大于130 时， 不同海区的响应系数将趋于稳定， 本文选用37年冬半年共计222个样本进行统计分析， 保证了反馈系数的估计精度和稳定性。同时研究中采用EOF方法建立主要海洋模态空间场， 可有效减少强迫场噪音， 并将运算结果进行蒙特卡洛检验（施能等， 2004）。
3 结果与分析
3.1 雪灾频数时空特征
统计1978 -2014年冬半年青藏高原雪灾频数， 由图1（a）可知， 青藏高原年平均雪灾频数为50.2次， 总体有减少的趋势， 气候倾向率为3.4次·（10a）-1， 其中20世纪80年代末期到90年代中期为雪灾的高发期， 1998年后雪灾显著减少， 1998 -2014年平均雪灾频数较1978 -1997年偏少22.4次。段安民等（2016）研究认为， 青藏高原1998年开始趋暖， 并有加速增暖的趋势， 而气温是积雪变化的主导因子（姜琪等， 2020）， 温度与雪灾相关性比降雪更高（黄晓清等， 2018）， 近40年冬春季气温升高导致了积雪深度、 积雪日数的减少（马荣， 2018）， 对高原雪灾趋势变化起着决定性作用。从冬半年各月来看， 雪灾集中发生在1月和2月， 平均发生次数分别达到28.8次和29.1次， 其次12月份， 达到20.3次， 而10月发生次数最少， 平均仅为3.5次（图略）。空间上雪灾频数有较大的区域性差异［图1（b）］， 总体表现出自北向南递增的分布形式， 高值区主要在喜马拉雅山脉北坡及嘉黎地区， 各站雪灾累计次数均在80次以上， 其中帕里发生雪灾次数最多， 达到105次； 其次是错那（100次）。位于青南牧区的清水河， 由于海拔高， 温度低， 雪灾累计发生次数达到50次， 是青海出现雪灾最多的地区， 而少值区主要集中在柴达木盆地及青海东部农业区， 雪灾发生频数小于10次。
3.2 海温异常对雪灾变率强迫作用的诊断
选取关键海区是进行GEFA 诊断的基础。将原始数据减去三次多项式拟合趋势线去除了雪灾频数与全球海温的变化趋势， 得到1978 -2014年冬半年青藏高原雪灾频数年际变化时间序列与全球同期海温距平场（SSTA）的相关分布（图 2）。从图2中可以看到， 热带太平洋180°E以东为显著正相关、 120°E以西的热带印度洋和南印度洋为显著正负相关分布、 南太平洋大部分地区以负相关为主， 其中西部负相关显著， 东部局部正相关明显。可见， 以上4个海区的海温异常都与高原雪灾频数有一定的关系， 但由于各海区温度间可能并非独立， 相关分析无法分离不同海区对雪灾异常的独立贡献， 因此引入GEFA方法分离不同海区海温异常对高原雪灾频数的影响。
图2 1978 -2014年青藏高原冬半年雪灾频数年际变化时间序列与全球 SSTA场的相关分布阴影为通过显著水平为95%的显著性检验 Fig.2 Correlation distribution between the interannual variation time series of snow disaster frequency on the Qinghai-Xizang Plateau and the global SSTA field in winter half year from 1978 to 2014.The shadow indicates the area that has passed the significance test at 95% level |
基于冬半年雪灾频数与全球 SSTA场的相关分布， 类似文献（Wen et al， 2010； 江志红等， 2013）的处理首先是把热带和中纬度大洋划分为4个非重叠的子海区， 分别为： 热带太平洋（TP， 20°S -20°N， 120°E -80°W）， 热带印度洋（TI， 20°S -20°N， 35°E -120°E）， 南太平洋（SP， 60°S -20°S， 120°E -80°W）和南印度洋（SI， 60°S -20°S， 35°E -120°E）。鉴于海温异常的一些主要物理模态可以通过 EOF 模来定义， 对这4个海区冬半年各月海温异常（SSTA）进行 EOF 分解， 以得到各海区的海温模态为强迫场。由表1可知， 各海区前三模态的累计方差贡献率为60.7%~89.3%， 可以体现各海域SSTA的主要特征。其中， 热带太平洋第一模态（TP1）主要表现为传统的El Niño海温异常空间特征， 第二空间模态（TP2）的海温异常分布呈“－＋”经向分布， 为热带太平洋东西向偶极型， 第三模态（TP3）中赤道中太平洋为负值， 其他区域为正值， 呈中太平洋La Nina 模态； 热带印度洋第一模态（TI1）表现为热带印度洋海域均为正值， 即海区一致模， 第二空间模态（ TI2）呈热带印度洋偶极子分布， 第三模态（TI3）海温呈东西经向反位相分布； 南印度洋第一空间模态（SI1）基本上为横跨副热带海区“－＋”温度异常分布， 第二模态SI2呈“－＋－”空间型， 第三模态（SI3）为印度洋中部正位相； 南太平洋第一空间模态（SP1）呈“＋－”纬向分布， 为南太平洋南北向偶极型， 第二模态（SP2）呈“－＋－”的三极子分布， 第三模态（SP3）为东西反位相分布（图3）。
表 1 各海区 SSTA 前三个模态的方差贡献率及其累积贡献率Table1 The variance contribution rate and cumulative contribution rate of the first three SST modes |
| 空间模态 | 贡献率/% | |||
|---|---|---|---|---|
| TP | TI | SI | SP | |
| EOF1 | 72.5 | 52.1 | 27.5 | 40.1 |
| EOF2 | 12.8 | 12.4 | 21.4 | 12.5 |
| EOF3 | 4.0 | 9.8 | 11.8 | 9.4 |
| EOF1+2+3 | 89.3 | 74.3 | 60.7 | 62.0 |
将上述4个海区 12 个 EOF 空间模态作为各个海区的强迫场， 以 1978 -2014年青藏高原冬半年雪灾频数作为大气响应场， 根据广义平衡反馈方法的公式计算出高原雪灾频数对 4 个海区的响应振幅估计值， 并用蒙特卡洛方法检验响应的显著性。从表2可以看出， 除南印度洋第二空间模态（SI2）、 南太平洋第二空间模态（SP2）与青藏高原雪灾频数为负响应外， 其他关键海温模态与高原雪灾频数的响应均为正响应， 其中通过90%蒙特卡洛显著性检验的模态有热带太平洋第一空间模态（TP1）、 热带印度洋第二空间模态（TI2）、 南太平洋第三空间模态（SP3）、 南印度洋第一空间模态（SI1）及第二空间模态（SI2）， 对应的响应振幅分别为3.6， 3.9， 2.3， 4.5和-3.5， 即： 当TP1模的时间主分量增加 1 个标准差时， 雪灾频数增加3.6次； 当 T12模的时间主分量增加1个标准差时， 雪灾频数增加3.9次； 同理， 当SP3模、 SI1模及SI2模的时间主分量增加或下降 1 个标准差时， 雪灾频数分别增加2.3， 4.5和3.5次（图4）， 表明各关键海区异常海温对青藏高原雪灾频数有明显影响。
表2 各海区 SSTA 前三个模态响应振幅及蒙特卡洛检验统计表Table 2 Statistic table of response amplitude and Monte Carlo test of the first three SSTA modes |
| 统计项目 | TP1 | TP2 | TP3 | TI1 | TI2 | TI3 | SI1 | SI2 | SI3 | SP1 | SP2 | SP3 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 蒙特卡洛检验 | 90.9 | 29.0 | 25.6 | 64.5 | 92.2 | 20.0 | 98.4 | 90.1 | 79.1 | 70.5 | 71.8 | 94.6 |
| 响应振幅/次 | 3.6* | 0.5 | 0.3 | 1.8 | 3.9* | 0.3 | 4.5* | -3.5* | 2.6 | 1.8 | -1.8 | 2.3* |
*代表响应振幅通过90%蒙特卡洛显著性检验 |
图4 热带太平洋（TP1, TP2, TP3）、 热带印度洋（TI1, TI2, TI3）、 南印度洋（SI1, SI2, SI3）、 南太平洋（SP1, SP2, SP3）4个海区海温前三个模态的高原雪灾频数GEFA响应振幅空心柱代表模态响应振幅通过90%蒙特卡洛显著性检验 Fig.4 GEFA response amplitude of Plateau snow disaster frequency for the first three modalities of SST in tropical Pacific (TP1, TP2, TP3), tropical Indian Ocean (TI1, TI2, TI3), South Indian Ocean (SI1, SI2, SI3) and South Pacific (SP1, SP2, SP3).The hollow column represents the amplitude of modal response that have passed the 90% Monte Carlo significance test |
挑选上述显著模态， 计算其随时间变率的响应值序列（图5）可以看出， 5个显著模态的响应值合成序列与青藏高原频数序列有着较好的相关性， 相关系数0.56， 通过99%显著性检验， 能较好地体现青藏高原雪灾频数的变化特征。统计4个海区的前三个海温模态响应值占原始青藏高原雪灾频数序列的方差贡献， 得到12个海温模态的方差贡献在0.11%~23.81%， 其中TP1、 TI2两个海温模态对青藏高原雪灾频数的方差贡献最大， 分别为23.81%、 22.11%（两者之和达到 45.9%）， 可见各海区中热带太平洋 TP1、 热带印度洋TI2模态对青藏高原雪灾频数变率的强迫作用较为显著。
3.3 关键SSTA模影响雪灾生成的可能过程
雪灾对海温异常的响应必然通过大气环流的异常实现。以500 hPa位势高度场、 600 hPa相对湿度场及风场作为大气响应场， 将4个海区的前三个EOF模态组合在一起代表海洋对大气的主要强迫场进行 GEFA 分析， 探讨关键区热带太平洋TP1模态、 热带印度洋TI2模态对青藏高原雪灾频数异常的强迫作用及其影响雪灾发生的可能物理过程。
为了对比模式响应结果的合理性， 选取冬半年雪灾频数标准差大于1的典型年份作为多雪灾年， 从多雪灾年500 hPa高度场合成异常场［图6（a）］可以看出， 极地至亚洲中高纬度位势高度场整体偏低， 在65°N -85°N、 40°E -100°E区域内存在异常低值中心， 负异常高达-46 gpm以上； 而贝加尔湖至中国东北部一带则为位势高度正异常， 形成阻高， 亚欧大陆中纬度地区位势高度异常从西到东呈现“＋－＋”配置， 形成典型的两脊一槽型， 青藏高原一带同样为高度场异常低值区， 负异常高达-16 gpm以上， 致使东亚大槽相对于气候态偏西， 槽区主要位于我国西部地区， 来自高纬度冷空气不断东移下滑至槽区所在的青藏高原地区， 受东部阻塞系统影响， 大量南下冷空气在青藏高原地区堆积， 加之高原上低值系统影响， 降水天气过程频发， 导致高原冬半年降雪易偏多。可见， 亚欧大陆中纬度地区位势高度异常从西到东“＋－＋”配置有利于降雪形成。
图6 冬半年典型多雪灾年500 hPa高度场合成异常场(a, 单位: ×10 gpm)及高度场对赤道太平洋第一空间模态的GEFA响应(b, 单位: ×10 gpm·℃-1）阴影区为通过90%的显著性检验的区域 Fig.6 The 500 hPa height field of the composed anomaly field (a, unit: ×10 gpm) and the GEFA response of the height field to the first spatial mode of the equatorial Pacific Ocean (b, unit: ×10 gpm·℃-1) in a typical snowy disaster year in winter half year.The shadow area is the area that passes the 90% significance test |
图6（b）为冬半年500 hPa高度场对热带太平洋TP1的GEFA响应结果， 在TP1的强迫下， 500 hPa高度场上西北太平洋地区正响应明显， 中高纬乌拉尔山附近正响应同样很明显， 西伯利亚地区与我国范围内均为明显的负响应， 东亚大槽偏弱， 欧亚中高纬度自西向东呈现“＋－＋”响应， 这种环流配置与上述实际多雪灾年环流形势较为一致， 有利于高原降雪的形成。这与韦志刚等（2005）研究的高原冬季多雪年大气环流异常特征也基本吻合。
TI2对500 hPa高度场的GEFA响应同样在欧亚中高纬自西向东呈现“＋－＋”的环流响应模式， 但在印度洋上空的响应并不显著（图略）。在此， 进一步分析中低层600 hPa相对湿度场及风场的GEFA响应， 探讨TI2海温对高原积雪的影响。从冬半年多雪灾年600 hPa湿度场及风场合成异常场 ［图7（a）］ 可以看出， 欧亚大陆中纬度湿度场增加， 且伊朗高原湿度增加幅度大于青藏高原， 在伊朗高压影响下， 阿拉伯海暖湿气流经伊朗高原输送至青藏高原， 西太平洋湿度场为显著负响应， 以东风气流为主， 在西北太平洋异常反气旋影响下， 西北太平洋湿润气流进入高原北部， 为高原异常降雪提供了水汽条件。
图7 冬半年多雪灾年600 hPa湿度场（等值线, 单位: %）与风场（矢量, 单位: m·s-1）合成异常场(a)和湿度场对热带印度洋第二空间模态 (b, 单位： %·℃-1)的GEFA响应阴影区为通过90%的显著性检验的区域 Fig.7 The combination of 600 hPa humidity field (isoline, unit: %) and wind field (vector, unit: m·s-1) into anomaly field (a) and 600 hPa humidity field to the GEFA response of tropical Indian Ocean second spatial mode (b, unit： %·℃-1) in a multi-snow disaster year in winter half year.The shadow area is the area that passes the 90% significance test |
相较500 hPa高空场， TI2对中低层相对湿度场的GEFA响应更显著。从图7（b）可以看到， TI2对600 hPa湿度场的响应结果与实际湿度异常值的分布形态极其相似， 南印度洋海温偶极子对其湿度场的响应呈现东西“＋－”响应模式， 响应幅度均为实际观测的20倍以上。可见， 500 hPa高空“＋－＋”环流形式引导冷空气南下在高原堆积， 并与600 hPa中低层湿热气流相交汇， 产生降水， 加之低温影响， 易在高原形成多雪灾年。
综上分析， 赤道中东太平洋El Niño型海温异常主要引起高层大气环流异常， 使东亚大槽偏弱， 新地岛及乌拉尔山地区形成阻高， 偏北气流引导冷空气从西伯利亚通道南下， 在高原堆积； 而印度洋偶极子型海温异常主要引起中低层水汽异常， 西北太平洋湿润气流在东风气流作用下进入高原北部， 阿拉伯海暖湿气流经伊朗高原在中纬西风作用下进入高原南部， 为高原降雪提供了水汽条件。可见， 赤道中东太平洋及印度洋海温异常是青藏高原雪灾生成的关键因子。
4 结论与讨论
基于青藏高原积雪观测数据及Hadley 中心海温资料， 采用GEFA-EOF相结合的最优反馈模分析方法， 针对青藏高原雪灾频数对关键区海温异常模态的响应问题进行了分析， 并讨论了不同海域的贡献度及产生的可能机制， 得到如下结论：
（1） 青藏高原雪灾频数与热带海洋海温异常有着密切联系， 赤道中东太平洋海温异常的El Niño型（TP1）和赤道印度洋海温偶极子模态（TI2）对雪灾有显著的强迫作用， 两模态对雪灾频数的贡献达到45.9%， 其中TP1的贡献最大。
（2） 当热带太平洋或热带印度洋SSTA有TP1、 TI2型的正位相海温强迫时， 雪灾频数响应振幅分别达到3.6次和3.9次， 表明青藏高原雪灾对关键海区异常海温均有明显的响应。
（3） 高原雪灾频数对海温异常的响应主要通过与之相联系的大气环流异常实现， 赤道中东太平洋及印度洋海温异常是青藏高原雪灾生成的关键因子。在500 hPa高度场上， 赤道中东太平洋El Niño型海温异常主要引起高层大气环流异常， 使东亚大槽偏弱， 新地岛及乌拉尔山地区形成阻高， 偏北气流引导冷空气从西伯利亚通道南下， 在高原堆积； 而印度洋偶极子型海温异常主要引起中低层水汽异常， 西北太平洋湿润气流在东风气流作用下经进入高原北部， 阿拉伯海暖湿气流经伊朗高原在中纬西风作用下进入高原南部， 为高原异常降雪提供了水汽条件。
雪灾变化成因复杂， 本文针对雪灾频数对海温异常的响应分析一定程度上反映了海温的强迫作用， 但值得注意的是， 海表温度异常并不是高原雪灾变化的外部强迫作用的唯一形式， 如陆面热力状况的变化以及人类活动等都是影响积雪异常的因子。本文只讨论了海温的情况， 其他因子对高原雪灾的影响仍有待于通过大气环流模式进一步深入研究。