气候统计分析方法及其应用气候统计分析方法及其应用-1-1魏凤英魏凤英魏凤英魏凤英中国气象科学研究院中国气象科学研究院中国气象科学研究院中国气象科学研究院weifycams.cma.gov.cnweifycams.cma.gov.cnweifycams.cma.gov.cnweifycams.cma.gov.cn主要内容主要内容气候统计分析的目的和步骤气候统计分析的目的和步骤气候趋势的分析方法气候趋势的分析方法气候突变的检测方法气候突变的检测方法气候统计分析气候统计分析气候统计分析方法是利用气候系统的统计特性对气候变气候统计分析方法是利用气候系统的统计特性对气候变 化及其异常进行分析化及其异常进行分析. .主要内容包括：主要内容包括： （1 1）应用统计方法了解区域性或全球性气候变化的时）应用统计方法了解区域性或全球性气候变化的时 空分布特征、变化规律和气候异常的程度空分布特征、变化规律和气候异常的程度, , 检测检测 气候信号；气候信号； （2 2）通过统计方法探索气候变量之间及与其它物理因）通过统计方法探索气候变量之间及与其它物理因 素之间的联系，研究气候异常的原因；素之间的联系，研究气候异常的原因； （3 3）利用统计方法对气候数值模拟结果与观测结果之）利用统计方法对气候数值模拟结果与观测结果之 间的差异进行分析。间的差异进行分析。气候统计分析的一般步骤气候统计分析的一般步骤1.1.收集资料收集资料 准确,精确 均一性,代表性,比较性 研究对象与样本量长度,区域大小有关2.2.资料预处理资料预处理 标准化处理3.3.分析方法的选择分析方法的选择根据研究的目的选择合适的方法4.4.科学综合和分析科学综合和分析 利用气候学知识进行分析,判断,切忌 多种方法计算结果的简单罗列气候时间序列的特征气候时间序列的特征数据取值随时间坐标而变化数据取值随时间坐标而变化; ;每一时刻取值的随机性每一时刻取值的随机性; ;数据之间的相关性和持续性数据之间的相关性和持续性; ;序列整体的上升或下降趋势序列整体的上升或下降趋势; ;在某一时刻出现转折或突变在某一时刻出现转折或突变; ;序列存在周期性振荡序列存在周期性振荡. .赤道中东太平洋海温的滞后相关赤道中东太平洋海温的滞后相关我国夏季降水的显著周期我国夏季降水的显著周期2.5年年10.7年年气候变化趋势分析气候变化趋势分析气候变化趋势分析气候变化趋势分析任何一个气候序列都可以看作由以下几个分量构成任何一个气候序列都可以看作由以下几个分量构成: :方法方法: : 线性倾向估计, 滑动平均, 累积距平,多项 式拟合, Mann方法 趋势趋势固有固有周期周期 循环循环周期周期 平稳平稳滑动平均滑动平均功能功能: : 利用时间序列的平滑值显示趋势变化利用时间序列的平滑值显示趋势变化. . j=1,2,n-k+1 应用实例1950-2005年北京年降水量11年滑动平均累积距平累积距平功能功能: : 利用曲线直观判断变化趋势及发生转折利用曲线直观判断变化趋势及发生转折 或突变的大致时间或突变的大致时间. . t=1,2,n应用实例1951-2000年登陆台风累积距平1950-20071950-2007年中国年平均气温累积距平年中国年平均气温累积距平|Z|=0.4035Z0.05=0.178其他平滑方法其他平滑方法五五, ,七七, ,九点二次平滑九点二次平滑, ,五点三次平滑五点三次平滑, ,多项多项式拟合式拟合功能功能: : 起到低通滤波的作用起到低通滤波的作用, ,更适合短时期变更适合短时期变 化趋势的分析化趋势的分析, ,可以克服滑动平均削弱过多可以克服滑动平均削弱过多波幅的缺点波幅的缺点. .序列两端平滑值的处理方法序列两端平滑值的处理方法平滑造成缺少序列两端平滑值平滑造成缺少序列两端平滑值, ,很难反映两端的真实趋势很难反映两端的真实趋势; ;将平滑视为具有非唯一边界约束问题，这样至少有三种最将平滑视为具有非唯一边界约束问题，这样至少有三种最低阶边界约束方案可以应用到平滑过程中低阶边界约束方案可以应用到平滑过程中: : 方案方案1 1：滑动序列的零阶导数，它可以生成最小模的滑动序列的零阶导数，它可以生成最小模的 解，此方案有利于序列边界附近的平滑趋势接近于气候态，解，此方案有利于序列边界附近的平滑趋势接近于气候态，记为记为Norm(Norm(模模) )约束方案；约束方案； 方案方案2 2：滑动序列的一阶导数，它可以生成最小斜率的滑动序列的一阶导数，它可以生成最小斜率的约束，有利于序列边界附近的平滑趋势接近一个局部值，约束，有利于序列边界附近的平滑趋势接近一个局部值，记为记为Slope(Slope(斜率斜率) )约束方案；约束方案； 方案方案3 3：滑动序列的二阶导数，生成最小粗糙度的解，滑动序列的二阶导数，生成最小粗糙度的解，有利于边界平滑趋势由一个定常斜率来逼近，记为有利于边界平滑趋势由一个定常斜率来逼近，记为Roughness(Roughness(粗糙度粗糙度) )约束方案约束方案. . 序列两端平滑值的计算步骤序列两端平滑值的计算步骤 (1) (1) 首先使用低通滤波平滑器或其它滤波器，对首先使用低通滤波平滑器或其它滤波器，对气候序列进行平滑。气候序列进行平滑。 (2) (2) 分别用上述三种边界约束方案计算出序列两分别用上述三种边界约束方案计算出序列两端的平滑值。端的平滑值。 (3) (3) 分别计算利用上述三种方案得到的平滑序列分别计算利用上述三种方案得到的平滑序列的均方误差（的均方误差（mean-square error, MSEmean-square error, MSE）, ,可以证可以证明，最小明，最小MSEMSE的平滑序列就是最优的平滑方案。的平滑序列就是最优的平滑方案。 对对1900-2002年冬季（年冬季（12-2月）月）AO指数序列做低指数序列做低通滤波平滑，滑动尺度分别取通滤波平滑，滑动尺度分别取10年和年和20年，然后年，然后计算滑动序列的计算滑动序列的Norm、Slope和和Roughness方方案作为填补序列两端的平滑值。案作为填补序列两端的平滑值。 三种约束方案的三种约束方案的MSEMSE NormNorm方案方案 SlopeSlope方案方案 RoughnessRoughness方案方案1010年滑动长度年滑动长度 0.7242 0.7160 0.7242 0.7160 0.7026 0.7026 2020年滑动长度年滑动长度 0.8263 0.8263 0.8228 0.8228 0.8850 0.88501900-20021900-2002年冬季年冬季AOAO指数及其指数及其3 3种方案种方案的的2020年平滑年平滑使用三种边界约束方案计算北京使用三种边界约束方案计算北京1724-2005年年降水量的年年降水量的10年平滑序列。年平滑序列。三种方案平滑序列的三种方案平滑序列的MSE分别为分别为: 0.7368, 0.7329,0.7369。Slope方案的方案的MSE较小较小. 北京北京1724-20051724-2005年夏季降水量年夏季降水量及其及其1010年平滑年平滑 北京各气候阶段夏季降水量平均值北京各气候阶段夏季降水量平均值(mm) 时段时段 1724-1773 1774-1815 1816-1839 1840-1852 1853-1868 趋势趋势 少少 多多 少少 多多 少少平均值平均值 366 512 409 510 370 时段时段 1869-1899 1900-1947 1948-1964 1965-2005 趋势趋势 多多 少少 多多 少少平均值平均值 553 436 555 402变化趋势的显著性检验变化趋势的显著性检验-1非参数统计检验,对变化趋势序列计算秩统计量 +1 当i 时刻以后数值大于该数值时 0 否则 变化趋势的显著性检验变化趋势的显著性检验-2计算统计量Z值在1-1之间变化.给定显著性水平,假定=0.05,判据为:如果ZZ0.05,则认为变化趋势在0.05显著性水平下是显著的.应用实例 对1951-2000年登陆台风累积距平的变化趋势序列进行显著性检验: (1) 计算秩序列 (2) 计算统计量Z Z= 0 .5331 Z0.05= 0.1913 (3) Z Z0.05,因此认为登陆台风的变化趋势是显著的.线性倾向估计线性倾向估计功能功能: : 建立气候序列建立气候序列x x与时间与时间t t之间的一元线性回归之间的一元线性回归, ,用一条合理的直线表示用一条合理的直线表示x x与与t t之间的关系之间的关系, ,判断序列判断序列整体上升或下降趋势整体上升或下降趋势. . t=1,2,n t=1,2,n a a常数常数. . b b倾向值倾向值,b0,b0时说明序列随时间呈上升趋势时说明序列随时间呈上升趋势; ; b0 brr|r, ,表明序列随时间变化表明序列随时间变化趋势是显著的趋势是显著的. .应用实例11950-2005年北京年降水量线性倾向年北京年降水量线性倾向B=-4.95 r=-0.377 r0.05=0.261951-1960平平均为均为782mm1991-2005年平年平均均518mm应用实例2全国年平均气温线性趋势分布图（图中标全国年平均气温线性趋势分布图（图中标“1”1”为上升为上升趋势，趋势，“0”0”为下降趋势）为下降趋势） 应用实例3 全国夏季气温线性趋势分布图（图中标夏季气温线性趋势分布图（图中标“1”1”为上升为上升趋趋 势，势，“0”0”为下降趋势）为下降趋势） 气候突变检测气候突变检测气候突变定义气候突变定义序列不连续的突发性变化;突变理论的精髓是关于奇点的理论:系统或过程从一个稳定状态到另一个稳定状态的飞跃;从统计观点而言,突变现象定义为从一个统计特性到另一个统计特性的急剧变化;目前还没有成熟的突变分析方法,只能借助统计检验的手段,多种方法比较和依靠气候知识进行判断十分重要多种方法比较和依靠气候知识进行判断十分重要. .气候突变的统计特性气候突变的统计特性均值和方差的突变均值和方差的突变 均值突变均值突变 方差突变方差突变 崔建新崔建新,周尚哲周尚哲(2001)气候突变的统计特性气候突变的统计特性趋势趋势和频率的突变和频率的突变 趋势突变趋势突变 频率突变频率突变 崔建新崔建新,周尚哲周尚哲(2001)气候突变的统计特性气候突变的统计特性回归系数突变回归系数突变概率突变概率突变分布形式的突变分布形式的突变针对非线性系统针对非线性系统: 动力学结构动力学结构(控制方程控制方程)突变突变滑动滑动t-t-检验检验功能: 考察一气候序列两组样本平均值是否存在显著性差异来检验突变.步骤: 1)1) 设置基准点,确定子序列的长度. 2)2) 滑动连续设置基准点,计算统计量 3)3) 给定显著性水平,若|ti|t则认为在基准点发生了突变. 其中其中 应用实例用滑动用滑动t-检验检测检验检测1950-2005年北京年降水量突变点年北京年降水量突变点.Cramers法法功能: 与与t-t-检验类似检验类似, ,区别在于它是比较一子区别在于它是比较一子 序列与总序列平均值的显著性差异序列与总序列平均值的显著性差异. . Yamamoto法法功能功能: : 利用信噪比检测突变利用信噪比检测突变. 气候变化信号气候变化信号变变率率-噪噪音音应用实例用用YamamotoYamamoto检测检测1950-20051950-2005年北京年降水量年北京年降水量, ,无突变点无突变点用用YamamotoYamamoto检测检测1911-20001911-2000年中国年平均气温等级突变年中国年平均气温等级突变Mann-Mann-kendallkendall法法-1-1功能：功能： 利用非参数统计检验方法检测突变。利用非参数统计检验方法检测突变。步骤：步骤：1 1）计算计算顺序顺序气候序列的秩序列气候序列的秩序列, ,即第即第I I时刻数时刻数值大于值大于 第第j j时刻数值个数的累计数：时刻数值个数的累计数： 其中其中 2 2）计算统计量计算统计量UFUFK K Mann-Mann-kendallkendall法法-2-2 3 3）计算计算逆序逆序气候序列的秩序列，并按气候序列的秩序列，并按2 2）计）计算统计量算统计量UBUBK K。 4 4）给定显著性水平给定显著性水平a=0.05, UFa=0.05, UFK,K,大于大于0 0表明表明序列呈上升趋势，小于序列呈上升趋势，小于0 0表明呈下降趋势，超过临表明呈下降趋势，超过临界线表明趋势显著，两条线交点为突变点。界线表明趋势显著，两条线交点为突变点。应用实例应用实例用用M-K检测检测1861-1988年北半球气温突变点年北半球气温突变点PettittPettitt方法方法功能功能: : 与Mann-kendall方法相似的非参数检验检测突变.步骤步骤: : 1)1) 构造秩序列分三种情况,即 2)2) 用秩序列直接检测突变点 3)3) 计算统计量 P0.05, 则认为检测出的突变点在统计意义上是显著的.LepageLepage法法-1-1功能功能: : 无分布双样本的非参数检验无分布双样本的非参数检验, ,将序列中两个子序列看作两将序列中两个子序列看作两个独立总体个独立总体, ,检验两子序列有无显著差异检验两子序列有无显著差异. .步骤步骤: : 1)1) 基准点之前样本基准点之前样本n n1 1, ,基准点之后样本基准点之后样本n n2, 2, n n1212为为n n1 1和和n n2 2之之和和, ,在在n n1212范围内计算秩序列范围内计算秩序列: : 2) 2) 构造构造WilcoxonWilcoxon统计量统计量 3)3) 构造构造AnsariyAnsariy-Bradley-Bradley统计量统计量 LepageLepage法法-2-2 4)4) 构造两者的联合统计量构造两者的联合统计量(Lepage统计量统计量) 5) 给定显著性水平给定显著性水平, 查查2 2分布表得到自由度为分布表得到自由度为2 2的临界值的临界值, ,当当WAWAi i超过临界值时超过临界值时, ,表明表明i i时刻前与时刻前与i i时刻后序列直接存在显著性差异时刻后序列直接存在显著性差异,i,i时刻发生了突时刻发生了突变变. . 应用实例应用实例-1用用LapageLapage法检测天津法检测天津1890-20021890-2002年年8 8月降水量突变月降水量突变显著性水显著性水平为平为0.010.01显著性水平为0.001 Bernaola Galvan 分割算法分割算法-1-1 BernaolaBernaola Galvan Galvan 分割算法主要思想是将时间序分割算法主要思想是将时间序列的突变检测问题视为一个分割问题，即将序列列的突变检测问题视为一个分割问题，即将序列看作为由多个具有不同平均值的子序列构成，此看作为由多个具有不同平均值的子序列构成，此方法的目的就是要找出各子序列之间最大差值的方法的目的就是要找出各子序列之间最大差值的平均值的位置。平均值的位置。对非平稳时间序列更适用对非平稳时间序列更适用. . Bernaola Galvan 分割算法分割算法-2-2 (1) (1) 对于序列从左至右以滑动的方式，分别计算每一个点对于序列从左至右以滑动的方式，分别计算每一个点位置的左边子序列的平均值和右边子序列的平均值，分别位置的左边子序列的平均值和右边子序列的平均值，分别将它们记作为将它们记作为leftleft和和rightright。 (2) 检验检验left和和right之间差异是否显著，计算统计量之间差异是否显著，计算统计量 Bernaola Galvan 分割算法分割算法-3-3 (3) 确定确定T达到最大值达到最大值Tmax的位置，计算的位置，计算Tmax的统计显著的统计显著性。具有性。具有 的可能截取点的显著性水平的可能截取点的显著性水平P P（）可以定义为可以定义为的概的概率率 P（）可以近似地由下式求得可以近似地由下式求得 (4) 对新分割的两个子序列重复过程对新分割的两个子序列重复过程(1)-(3)(1)-(3)，如子序列左、，如子序列左、右的子序列间的均值差异又满足上述条件，则对子序列进行右的子序列间的均值差异又满足上述条件，则对子序列进行再分割。如此重复，直至分割的子序列不再可以分割为止。再分割。如此重复，直至分割的子序列不再可以分割为止。 北半球树木年轮距平宽度序列北半球树木年轮距平宽度序列(AD1-1980)BG(AD1-1980)BG突变突变 100年尺度滤波(a) 原序列 (b) 低频序列 (c) 高频序列 BG突变 封国林等封国林等(2006)(2006)突变检测遇到的问题突变检测遇到的问题样本量大小变量的变化趋势考察的时间尺度