pettitt检验r语言_多种均一性检验方法比较研究

引言

随着全球气候变暖，气候要素序列的均一性检验愈来愈受到科学家的关注[。所谓均一性序列，是指测站得到的气象记录序列仅仅是气候实际变化的反映[。然而大多数气候要素序列都受到非气候因素的影响，包括台站迁移、仪器变更和台站环境改变等[。利用这些非均一性气候序列分析气候变化，可能导致非真实结论，进而影响气候研究以及业务工作的开展。

国外均一性研究起步较早。自Mitchell在1966年首次提出均一性检验后，Alexandersson等[、Solow[、Lund等[相继发展了SNHT、TPR、Potter等均一性检验方法，并用于北美和欧洲气候数据处理。Slonosky等[和Wijingard等[利用Buishand、Pettitt、von Neumann方法检验ECA日平均气温、降水和气压数据的均一性。国内，气象资料的均一性检验和订正工作也逐渐受到学者的关注，刘小宁[、吴利红[和宋超辉等[总结了气候资料序列的均一性检验方法，并应用于我国部分台站气温和降水序列。李庆祥等[、李庆祥[和高理[利用多种方法对中国东南部地区气候序列进行了均一性检验，力求减少对台站沿革数据的依赖。

以往均一性检验主要集中在我国东部地区，对西南区域气候资料研究较少。四川省位于中国西南部，地形复杂多样。根据地形特点，本文将四川省气象台站分为盆地底部地区(99个)，盆地边缘和川西南山地区(14个)，川西北高原区(47个)(见

图 1 四川省台站分布图 Fig. 1 Station distribution in Sichuan Province

1 资料与检验方法

1.1 均一化思路

考虑到四川省1960年以来气温资料的相对完整性，本文从160个台站1960—2009年气温资料中，选取具有连续观测记录50年以上的站点105个。以周边台站为参考站，选相关系数最高的5个台站为待检站，利用比值法构造参考序列Yi。依据Wijingaard等[所设定均一性检验思路：通过不同统计方法判断序列的均一性：4种及4种以上方法检验均通过的台站资料为可用，2种不通过的为可疑，至少3种不通过的为不可信。

1.2 方法介绍

本文选取SNHT、TPR、Buishand、Pettitt、MASH[和von Neumann法对年平均气温序列进行均一性检验。各方法原理在文献[中已有详细介绍。本文总结各方法具体差异，见

表 1

Table 1

表 1 各种检验方法对比

Table 1 Comparison of various test methods

检验方法基本原理检验变量的正态性是否建立参考序列主要原理公式断点判断依据断点年断点敏感度

SNHT最大似然法假设正态选取相关性高的邻近台站$begin{align} &{{Z}_{i}}=({{Y}_{i}}-bar{Y})/{{sigma }_{Q}} \ &T_{max }^{s}=maxlimits_{1le ale n-1}{{T}^{s}}=maxlimits_{1le ale n-1}{a{{{bar{z}}}_{1}}+(n-a){{{bar{z}}}_{2}}} end{align}$$T_{max }^{s}$大于某临界值水平提供待检序列前端和靠近结尾部分

Buishand距平累加法假设正态选取相关性高的邻近台站$begin{align} &S_{0}^{*}=0,S_{k}^{*}=sumlimits_{i=1}^{k}{({{Y}_{i}}-bar{Y})} \ &R=(maxlimits_{0le kle n}S_{k}^{*}-minlimits_{0le kle n}S_{k}^{*})/s end{align}$$R/sqrt{n}$大于某临界值水平提供待检序列中部

Pettitt非参数秩检验秩检验选取相关性高的邻近台站r1，…，rn是检验序列Y1，…，Yn的秩序列

${{X}_{k}}=2sumlimits_{i=1}^{k}{{{r}_{i}}-k(n+1)}$Xk大于某临界值水平提供待检序列中部

von Neumann最大似然法假设随机分布选取相关性高的邻近台站$N={{sumlimits_{i=1}^{n-1}{{{({{Y}_{i}}-{{Y}_{i+1}})}^{2}}/sumlimits_{i=1}^{n}{({{Y}_{i}}-bar{Y})}}}^{2}}$N值大于某临界值水平不提供不适用

MASH序列均一性的多元分析检验假设正态选取相关性高的邻近台站$begin{align}

&Z_{c}^{(m)}(t)={{X}_{c}}(t)-sumlimits_{iin R_{c}^{(m)}}{{{omega }_{i}}cdot {{X}_{i}}(t)}, \

&sum{{{omega }_{i}}}=1,{{omega }_{i}}ge 0, \

&V(Z_{c}^{(m)})=Variance(Z_{c}^{(m)})=underset{omega }{mathop{min imum}}

end{align}$最优差异序列的统计量大于某临界值提供待检序列全部

TPR线性回归拟合假设正态不需要建立对任一时间点前后的序列进行线性回归模拟$X=left{ begin{align}

&{{mu }_{1}}+{{alpha }_{1}}t+{{varepsilon }_{t}},1le tle c \

&{{mu }_{2}}+{{alpha }_{2}}t+{{varepsilon }_{t}},cle tle n \

end{align} right},$将最小二乘估计$left{ begin{align}

&{{{hat{mu }}}_{1}}={{{bar{x}}}_{1}}-{{{hat{alpha }}}_{1}}{{{bar{t}}}_{1}} \

&{{{hat{mu }}}_{2}}={{{bar{x}}}_{2}}-{{{hat{alpha }}}_{2}}{{{bar{t}}}_{2}} \

end{align} right}$调整为$begin{align}

&F=frac{(SS{{E}_{text{red}}}-SS{{E}_{text{full}}})/2}{SS{{E}_{text{full}}}/(n-4)} \

&SS{{E}_{text{full}}}=sumlimits_{t=1}^{c}{{{({{x}_{t}}-{{{hat{mu }}}_{1}}-{{{hat{alpha }}}_{1}}t)}^{2}}}+ \

&sumlimits_{t=c+1}^{n}{{{({{x}_{t}}-{{{hat{mu }}}_{2}}-{{{hat{alpha }}}_{2}}t)}^{2}}} \

&SS{{E}_{text{red}}}=sumlimits_{t=1}^{n}{{{({{x}_{t}}-{{{hat{mu }}}_{text{red}}}-{{{hat{alpha }}}_{text{red}}}t)}^{2}}} \

end{align}$F统计值>Fmax提供待检序列中部

表 1 各种检验方法对比

Table 1 Comparison of various test methods

2 四川省年平均气温序列个例均一性检验分析

以汶川站(56183) 为例，6种方法均以α=0.01显著性水平下的临界值检验。

图 2 汶川站SNHT(a)、Buishand检验(b)和Pettitt检验(c)方法年平均气温的均一性检验图(虚线为α=0.01显著性水平检验，点线为α=0.05显著性水平检验) Fig. 2 The average temperature of Wenchuan Station in the homogeneity test chart by using SNHT (a), Buishand range test (b), and Pettitt test (c)

(Dashed line is the 1% significance level, dotted lines for the 5% significance level)

图 3 汶川站TPR方法平均气温的均一性检验图(虚线为α=0.01显著性水平检验) Fig. 3 The average temperature of Wenchuan Station in the homogeneity test chart by using TPR test (Dashed line is the 1% significance level)

图 4 MASH检验汶川站平均气温差异序列及其转折值 Fig. 4 The difference series and the shift values identified by MASH for the average temperature of Wenchuan Station

由α=0.01显著性水平检验。由

根据台站沿革资料，1977年7月1日，该站由威州镇姜维城“山腰”迁至汶川县城南部边缘处。

表 2

Table 2

表 2 12个台站不同方法间断序列的统计

Table 2 The statistics of discontinuous series under different methods in 12 stations

地形站名各方法检验检验不均一方法个数检验结果原因

SNHT断点年份Buishand断点年份Pettitt断点年份TPR断点年份MASH断点年份Von Neumann结果

盆地剑阁19871987198519861986不均一6不均一1987迁站

高坪无无无19801971均一2怀疑1980迁站

珙县19791979无20021978均一3不均一1979迁站，2002仪器变更

遂宁20012001无19992001均一4不均一2001迁站

山地都江堰1992199219861989无均一3不均一原因不明

宝兴2007无无1996无均一1均一无迁站，无仪器变更

米易无无无无1965、1971均一1均一无迁站，无仪器变更

石棉1963无无19631964均一3不均一1963迁站

高原金阳20062006无20061978、2006不均一5不均一2006仪器变更

盐源19771977无19771977不均一5不均一1977迁站

木里19901990198919891989不均一6不均一1990迁站

乡城无无无无无均一0均一无迁站，无仪器变更

表 2 12个台站不同方法间断序列的统计

Table 2 The statistics of discontinuous series under different methods in 12 stations

3 气温均一性综合检验分析

3.1 断点检验的敏感度

文中所选6种方法均要求原假设为独立、同分布，但拒绝假设存在差异：SNHT、Buishand、Pettitt和TPR法假设在平均值上存在逐步转变(即断点)，MASH检验中最优差异序列的统计量通过Monte Carlo法与临界值相比来确定断点位置，这5种方法检验时可定位断点所在年，为特定位置检验法。von Neumann法在拒绝假设中设定该序列非随机分布，因此无法给出断点所在年份的信息。此外，SNHT和Buishand法均假设Yi值为正态分布，但Pettitt法无此限制，它基于序列中要素的排列方式，并非值本身。Pettitt法的排列方式使其对“离群点”的敏感度减弱。TPR法将断点时间作为不可知参数处理，避免data-dredging[。为进一步说明特定位置检验法对断点敏感度的区别，现将检验结果去掉突变点位于序列前后两年的站点(虚假不连续点)后分年代统计(见[的研究结论相符。

图 5 不同方法检验温度序列断点敏感度 Fig. 5 The breakpoint sensitivity testing of temperature series in different methods

3.2 检验结果相符率

为了表征检验结果的相似性，借鉴吴利红[的“相符率”概念，相符率=检验结果相同站数/被检验站总数。对温度检验结果相符率分析(见

表 3

Table 3

表 3 温度均一性检验结果相符率(单位：%)

Table 3 Matching rate of temperature homogeneity test results (unit: %)

2种方法之间的相符率SBSPSTBPBTPTSMBMPMTM

68.644.819.069.522.918.122.915.212.424.8

3种方法之间的相符率SBPSBTSBMSPTSPMSTMBPTBPMBTMPTM

48.619.014.315.29.59.516.28.68.65.7

4种方法之间的相符率SBPTSBPMSBTMSPTMBPTM

14.38.97.63.84.8

5种方法之间的相符率SBPTM

5.7

表中S为SNHT法，B为Buishand法，P为Pettitt法，T为TPR法，M为MASH法。

表 3 温度均一性检验结果相符率(单位：%)

Table 3 Matching rate of temperature homogeneity test results (unit: %)

3.3 年平均气温序列均一性检验结果

由四川省105个站点气温序列的检验结果(见

表 4

Table 4

表 4 四川省年平均气温序列均一性检验结果

Table 4 Homogeneity test results of temperature sequence in Sichuan Province

拒绝原假设的个数是否通过均一性检验站数所占比率

3~6不通过3028.6%

2可疑4341.0%

0~1通过3230.5%

表 4 四川省年平均气温序列均一性检验结果

Table 4 Homogeneity test results of temperature sequence in Sichuan Province

图 6 四川省气温序列通过均一性检验(实心点)和未通过均一性检验(空心点)台站分布图 Fig. 6 Station distribution of the temperature series in Sichuan Province which passed the homogeneity test (dark dots) and not (circles)

分测站类型进行温度均一性检验时发现，基准站有5个，无不均一台站；基本站24个，均一占20.8%；一般站76个，均一台站占28.9%。由于四川地形变化明显，气候资料代表的区域范围较为有限。因此，大范围的气候变化研究宜选择代表性范围尽可能大的站点。

借鉴吴利红[对断点原因的分类，本文将非均一性原因作如下定义：间断点与测站迁站时间基本吻合(允许误差5年)为迁站原因；由于仪器重新安置而导致资料序列的非均一性为仪器变化原因；其他无明显环境破坏记载，定义为不明原因。联系台站沿革数据，将非均一性台站分类，见

表 5

Table 5

表 5 四川省年平均气温序列非均一性原因统计

Table 5 The causes for inhomogenization on the temperature sequence

断点原因台站迁移仪器变更不明原因

台站个数281836

表 5 四川省年平均气温序列非均一性原因统计

Table 5 The causes for inhomogenization on the temperature sequence

由于地理条件的差异，各台站对检验方法的适用性也不同，以降低误判断点为前提，利用台站沿革数据验证各方法所得断点的真实性(见

表 6

Table 6

表 6 不同地形各方法判断气温间断的统计

Table 6 The statistics of different methods to determine the breakpoint of temperature in different regions

SNHTBuishandPettittvon NenmannTPRMASH

N1N2N1N2N1N2N1N2N1N2N1N2

盆地底部区域342114614330不可查57123813

盆地边缘和川西南山地区域7431527不可查6163

川西高原区域2111181120615不可查145147

注：N1为非均一台站数，N2为沿革数据验证的非均一台站数。

Note: N1 is number of inhomogeneous stations, N2 is number of inhomogeneous stations for data valodation in historical evolution.

表 6 不同地形各方法判断气温间断的统计

Table 6 The statistics of different methods to determine the breakpoint of temperature in different regions

表 7

Table 7

表 7 不同地形下SNHT法对温度均一性检验

Table 7 SNHT to determine the breakpoint of temperature in different regions

无缺测台站数显露数台站迁移台站数仪器变更台站数误判台站数漏检台站数

盆地底部区域63341211144

盆地边缘和川

西南山地区域1073140

川西高原区域322184100

表 7 不同地形下SNHT法对温度均一性检验

Table 7 SNHT to determine the breakpoint of temperature in different regions

同理对SNHT、Buishand、Pettitt、TPR及MASH法进行比较，并参考各台站具体数据(见

表 8

Table 8

表 8 温度均一性检验各方法漏检和误判率(单位：%)

Table 8 The missed rate and false positives of temperature in different regions with different methods (unit: %)

SNHT法Buishand法Pettitt法TPR法MASH法

漏测率误判率漏测率误判率漏测率误判率漏测率误判率漏测率误判率

盆地底部区域6.1521.525.412.725.415.94.873.017.525.4

盆地边缘和川西南山地区域0402020103010703040

川西高原区域031.3028.13.143.89.421.92518.75

表 8 温度均一性检验各方法漏检和误判率(单位：%)

Table 8 The missed rate and false positives of temperature in different regions with different methods (unit: %)

4 结论与讨论

4.1 资料的可用性分析

(1) 在显著性水平α=0.01的检验中，四川省105个待检台站年平均气温序列通过均一性检验的有32个，占30.5%；其中均一的基准站有5个，基本站有5个，一般站有22个。均一的基本站和一般站的使用将扩展气候资料代表的区域范围，提高气候资料的使用率，有助于四川省的气候变化研究。

(2) 可疑台站有31个，占29.5%；造成这类台站无法确定断点的原因，一方面可能由其他因素产生，如气候变化、环境变化、观测方法改变等原因，需要进一步分析；另一方面可能由于台站沿革资料不完备造成，致使断点信息无法查阅。因此对于不明原因的断点序列，必须通过更加详细的历史沿革资料判断其真实性。

(3) 不均一台站42个，占40%，利用这部分温度资料序列进行气候气候变化分析研究时应注意资料的订正。

4.2 非均一性原因

序列不均一原因中，因迁站导致的有28个，更换仪器引起的有18个。分年代统计，20世纪60—70、80—90年代以及2000年以后不均一台站有增加的趋势，其中60—70和80—90年代不均一性均由台站迁移所致。再次证实，迁站是引起序列间断的主要原因之一。

4.3 检验方法评估分析

对不同方法的敏感度和适用性作评估分析发现，SNHT和TPR法易在序列前端和后端检测断点，Pettitt法对序列中部断点较敏感，MASH对序列各部断点检测均较敏感，这与其差异序列的数学原理及检验效率有关。SNHT与Buishand和Pettitt法检验结果相符率为48.6%，TPR法与其他方法相符率偏低，但该法不涉及邻近站是否均一，因此检验中可做补充判断。根据四川省特殊地形分析，发现对盆地及山地台站而言，SNHT法显露率最高，分别为67.6%和57.1%，漏检和误判较少；高原地区SNHT法和Buishand法显露率较高，其中Buishand法漏检和误判率最低。

因此，鉴于减小对元数据依赖度的均一化检验思路，在缺乏元数据时，应采用多种检验方法综合分析，有助于气候资料在区域气候变化研究中的有效利用。

致谢：上海台风研究所陈葆德老师、占瑞芬老师，中科院大气所李珍老师以及四川省气候中心孙昭萱、钟燕川同志对本工作给予很多指点和帮助，谨致谢忱。