简介： 本文通过蒙特卡洛模拟对Kendall非参数秩次相关检验法在水文时间序列趋势检验中的识别能力做了研究。模拟结果表明，该方法的识别能力与之前给定的显著性水平、容量、趋势度和变差系数有关。随着趋势度的绝对值、容量和显著性水平的增加，Kendall检验法的识别能力增强；随着变差系数的增加，该方法的识别能力减小；当趋势存在时，该方法还与时间序列服从的分布类型和形状参数有关。最后采用Kendall秩次检验法分析了普渡河流域海口、蔡家村和三江口水文站的年径流系列，发现蔡家村和三江口的年径流有增加的趋势，找出了增加的原因。
关键字：Mann-Kendall 非参数秩次相关检验法 蒙塔卡洛模拟 水文时间序列 趋势分析
1、引言
水文序列是从工程所在地点或邻近地点水文观测的资料中选取表征水文过程特征值（如洪峰流量或水位、各种时段洪水总量等）的。它们是进行频率分析、估计设计水文过程的基础资料。在水文资料的观测期内，如因流域上修建了蓄水、引水、水土保持等工程以及流域的气候发生变化等等，这些人工或天然的原因使流域水文现象的形成条件发生了显著的改变，因而水文变量的概率分布规律也发生了显著的变异，我们把这一问题称为水文资料的“非一致性”问题[1]。如果将“非一致性”的水文资料混杂在一起作为一个进行水文频率计算显然违背了用于频率分析的水文序列必须服从同一分布的要求。
因此，进行频率分析之前，必须检验序列是否具有一致性，其中趋势检验是一致性检验的内容之一。在水文资料分析研究中，非参数秩次相关统计检验，即Kendall检验用于识别时间序列是否存在趋势成分。本文首先简单介绍了MK趋势检验法；研究了该方法识别正态分布序列趋势的能力与给定的显著性水平、容量、趋势度和变差系数的关系；同时探讨该方法识别非正态分布序列趋势的能力与分布类型和形状参数的关系；最后举出该方法的一个应用实例。
2、Kendall秩次相关检验
Kendall非参数秩次相关检验法已经广泛的用于检验水文气象资料的趋势成分，包括水质、流量、气温和降雨序列等，究其原因主要是，与参数统计检验法相比，非参数检验法更适用于非正态分布或经过删检（删去低于或高于某水平的观测值）的资料，而这些情况在时间序列分析中常常会遇到。过去20年里，国际上关于MK方法应用研究的实例非常之多[2-4]。尽管该方法应用如此之广泛，但我们还是不清楚该方法是不是适用于各种情况下的时间序列的趋势检验。
对序列
，先确定所有对偶值
中
与
的大小关系（设为
）。趋势检验的统计量[5]为：
（1）
式中：
；
（2）
（3）
当
大于10时，
收敛于标准正态分布。
原假设为该序列无趋势，采用双边趋势检验，在给定显著性水平
下，在正态分布表中查得临界值
，当
时，接受原假设，即趋势不显著；若
，则拒绝原假设，即认为趋势显著。
3、Kendall检验对正态分布序列趋势的识别能力
3.1 两类错误与识别能力
假设检验[7]中有两类错误：第一类错误为“弃真”，即原假设本来是正确的，但检验结果却拒绝了原假设，由数理统计知识可知，犯第一类错误的概率为显著性水平
；第二类错误为“纳伪”，即原假设本来是错误的，但检验结果却接受了原假设。犯第二类错误的概率为
。
定义两种检验方法的识别能力为
，也就是当原假设错误时，该检验方法正确拒绝原假设的能力。
3.2 蒙特卡洛模拟试验步骤
因为水文时间序列
可看成是多种成分组成。假定这些成分是线性叠加，
可按下式表示[5]：
式中
为确定性的非周期成分（包括趋势
，跳跃
等暂态成分等）；
为确定性的周期成分（包括简单的或复合周期的成分等）；
为随机成分（包括平稳的或非平稳的随机成分）。本次研究的时间序列只考虑确定性的非周期成分中的趋势项和随机成分。
蒙特卡洛模拟用于评价Kendall检验识别趋势成分的能力，试验步骤如下：
①令
；
②生成正态分布的随机序列
，容量可取
，均值
均取1.0，方差
分别为
，其中
；的标准偏差
和变差系数
为
；
③生成趋势成分
，
，
为趋势度，可取-0.01，-0.005，0，0.005和0.01等；
④将②和③中生成的序列按序号对应相加得到有趋势的时间序列
；
⑤原假设为序列不存在趋势，采用Kendall检验法对序列
进行趋势检验。如果拒绝原假设，则
，否则
；
⑥重复步骤②~⑤
次；
⑦根据前面的定义，Kendall检验的识别能力为：
式中：
为蒙特卡洛模拟的次数；
为检验落在拒绝域里的次数。
3.3 结果分析
3.3.1 识别能力与显著性水平、趋势度的关系
容量
取50；均值取1.0；变差系数
取0.5；显著性水平
取0.002，0.005，0.01，0.025，0.05，0.10，0.15，0.20；趋势度取-0.01，-0.005，0，0.005和0.01。Kendall趋势检验法的识别能力与显著性水平和趋势度的关系见表1和图1。从蒙特卡洛模拟结果可以看出，对于同一显著性水平，趋势度绝对值越大，则检验法的识别能力越强；而对于同一趋势度，显著性水平越大，则检验法的识别能力越强；特别的，当
时，序列无趋势存在，原假设是正确的，此时，检验法的识别能力均等于事先给定的显著性水平。
表1Kendall的识别能力与显著性水平和趋势度的关系
|
趋势度
|
显著性水平
|
|||||||
|
0.002
|
0.005
|
0.010
|
0.025
|
0.050
|
0.100
|
0.150
|
0.200
|
|
|
-0.010
|
0.114
|
0.180
|
0.252
|
0.380
|
0.498
|
0.629
|
0.702
|
0.757
|
|
-0.005
|
0.009
|
0.026
|
0.050
|
0.101
|
0.161
|
0.269
|
0.335
|
0.391
|
|
0.000
|
0.002
|
0.005
|
0.010
|
0.025
|
0.050
|
0.100
|
0.150
|
0.200
|
|
0.005
|
0.015
|
0.029
|
0.051
|
0.098
|
0.162
|
0.252
|
0.320
|
0.378
|
|
0.010
|
0.103
|
0.170
|
0.244
|
0.361
|
0.488
|
0.615
|
0.694
|
0.754
|
图1Kendall的识别能力与显著性水平和趋势度的关系
3.3.2 识别能力与容量、趋势度的关系
显著性水平
取0.05；均值取1.0；变差系数
取0.5；容量
取10~100；趋势度取-0.01，-0.005，0，0.005和0.01。Kendall趋势检验法的识别能力与容量和趋势度的关系见表2和图2。从蒙特卡洛模拟结果可以看出，随着容量和趋势度绝对值的增加，检验法的识别能力增强；特别的，当
时，序列无趋势存在，原假设是正确的，此时，检验法的识别能力均等于事先给定的显著性水平。
表2MK的识别能力与容量和趋势度的关系
|
趋势度
|
容量
|
||||||||
|
20
|
30
|
40
|
50
|
60
|
70
|
80
|
90
|
100
|
|
|
-0.010
|
0.069
|
0.126
|
0.284
|
0.498
|
0.729
|
0.897
|
0.975
|
0.996
|
0.999
|
|
-0.005
|
0.049
|
0.063
|
0.110
|
0.161
|
0.240
|
0.370
|
0.527
|
0.666
|
0.790
|
|
0.000
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
|
0.005
|
0.045
|
0.079
|
0.094
|
0.162
|
0.251
|
0.367
|
0.493
|
0.646
|
0.781
|
|
0.010
|
0.065
|
0.141
|
0.252
|
0.488
|
0.727
|
0.894
|
0.975
|
0.997
|
0.999
|
图2Kendall的识别能力与容量和趋势度的关系
3.3.3 识别能力与趋势度、变差系数的关系
容量
取50；显著性水平
取0.05；均值为1.0；变差系数
取0.01~1.0；趋势度取-0.01，-0.005，0，0.005和0.01。Kendall趋势检验法的识别能力与变差系数、趋势度的关系见表3和图3。从蒙特卡洛模拟结果可以看出，对于同一变差系数，趋势度绝对值越大，则检验法的识别能力越强；而对于同一趋势度，变差系数越大，则检验法越难识别出序列存在趋势成分；特别的，当
时，序列无趋势存在，原假设是正确的，此时，检验法的识别能力均等于事先给定的显著性水平。
表3Kendall的识别能力与变差系数和趋势度的关系
|
趋势度
|
变差系数
|
||||||||||
|
0.01
|
0.10
|
0.20
|
0.30
|
0.40
|
0.50
|
0.60
|
0.70
|
0.80
|
0.90
|
1.00
|
|
|
-0.010
|
1.000
|
1.000
|
0.999
|
0.904
|
0.686
|
0.498
|
0.368
|
0.292
|
0.238
|
0.193
|
0.161
|
|
-0.005
|
1.000
|
0.999
|
0.686
|
0.368
|
0.238
|
0.161
|
0.129
|
0.107
|
0.096
|
0.085
|
0.076
|
|
0.000
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
0.050
|
|
0.005
|
1.000
|
0.999
|
0.669
|
0.354
|
0.223
|
0.162
|
0.123
|
0.104
|
0.088
|
0.078
|
0.072
|
|
0.010
|
1.000
|
1.000
|
0.999
|
0.893
|
0.669
|
0.488
|
0.354
|
0.279
|
0.223
|
0.185
|
0.162
|
图3Kendall的识别能力与变差系数和趋势度的关系
3.3.4 识别能力与容量、变差系数的关系
显著性水平
取0.05；均值为1.0；趋势度取0.005；容量
取10~100；变差系数
取0.01~1.0。两种趋势检验法的识别能力与容量和变差系数的关系见表4和图4。从蒙特卡洛模拟结果可以看出，对于同一变差系数，检验法的识别能力随着容量的增加而增强；而对于同一容量，检验法的识别能力则随着变差系数的增加而减小；特别的，当
时，序列无趋势存在，原假设是正确的，此时，检验法的识别能力均等于事先给定的显著性水平。
表4Kendall的识别能力与容量和变差系数的关系
|
容量
|
方差系数
|
||||||||||
|
0.01
|
0.10
|
0.20
|
0.30
|
0.40
|
0.50
|
0.60
|
0.70
|
0.80
|
0.90
|
1.00
|
|
|
20
|
1.000
|
0.200
|
0.080
|
0.056
|
0.050
|
0.045
|
0.041
|
0.041
|
0.040
|
0.041
|
0.041
|
|
30
|
1.000
|
0.596
|
0.199
|
0.112
|
0.088
|
0.079
|
0.070
|
0.063
|
0.060
|
0.059
|
0.057
|
|
40
|
1.000
|
0.920
|
0.387
|
0.210
|
0.134
|
0.106
|
0.084
|
0.073
|
0.066
|
0.061
|
0.055
|
|
50
|
1.000
|
0.998
|
0.684
|
0.358
|
0.249
|
0.166
|
0.138
|
0.106
|
0.090
|
0.089
|
0.083
|
|
60
|
1.000
|
1.000
|
0.878
|
0.559
|
0.361
|
0.238
|
0.182
|
0.148
|
0.134
|
0.117
|
0.098
|
|
70
|
1.000
|
1.000
|
0.984
|
0.765
|
0.528
|
0.364
|
0.259
|
0.202
|
0.159
|
0.142
|
0.128
|
|
80
|
1.000
|
1.000
|
0.999
|
0.908
|
0.691
|
0.483
|
0.378
|
0.310
|
0.237
|
0.199
|
0.169
|
|
90
|
1.000
|
1.000
|
1.000
|
0.979
|
0.856
|
0.663
|
0.516
|
0.397
|
0.316
|
0.281
|
0.218
|
|
100
|
1.000
|
1.000
|
1.000
|
0.996
|
0.937
|
0.817
|
0.627
|
0.518
|
0.428
|
0.329
|
0.278
|
图4Kendall的识别能力与容量和变差系数的关系
4、Kendall检验对非正态分布序列趋势的识别能力
前面已经讨论了Kendall检验对正态分布序列趋势的识别能力，然而，水文时间序列通常是有偏的，很少有服从正态分布的，因此，下面进一步探讨Kendall检验法识别非正态分布序列趋势的能力。这里比较的都是水文中常用的几种分布类型，包括三类极值分布（EV1、EV2和EV3）、对数正态分布（LN）和皮尔逊Ⅲ型分布（PE3）。
4.1 识别能力与分布类型的关系
为了方便比较研究，容量取为50，均值和方差分别取为1.0和0.5，形状参数取值见表5；趋势度
仍取-0.01，-0.005，0，0.005和0.01等；蒙特卡洛模拟次数取2000；显著性水平取0.05。根据3.2的步骤评价两种检验法的识别能力，检验结果列于表5中，并以图5来表示。
表5 Kendall的识别能力与分布类型的关系
|
分布类型
|
趋势度
|
||||||||||
|
-0.010
|
-0.008
|
-0.006
|
-0.004
|
-0.002
|
0.000
|
0.002
|
0.004
|
0.006
|
0.008
|
0.010
|
|
|
PE3（γ=1.5）
|
0.685
|
0.495
|
0.318
|
0.173
|
0.079
|
0.050
|
0.086
|
0.162
|
0.317
|
0.499
|
0.656
|
|
EV1（γ=0.0）
|
0.400
|
0.273
|
0.180
|
0.111
|
0.064
|
0.050
|
0.076
|
0.100
|
0.182
|
0.282
|
0.399
|
|
EV2（γ=-1.5）
|
0.363
|
0.263
|
0.187
|
0.117
|
0.069
|
0.050
|
0.080
|
0.109
|
0.188
|
0.277
|
0.357
|
|
EV3（γ= 1.5）
|
0.913
|
0.831
|
0.656
|
0.447
|
0.188
|
0.050
|
0.172
|
0.424
|
0.691
|
0.841
|
0.917
|
|
LN（γ=1.5）
|
0.692
|
0.548
|
0.373
|
0.236
|
0.100
|
0.050
|
0.101
|
0.222
|
0.391
|
0.545
|
0.673
|
从模拟结果可以看出，当
，即没有趋势存在时，检验法对所有分布的识别能力
等于之前给定的显著性水平0.05，说明检验统计量对序列服从的分布类型并不敏感。然而，当
时，即有趋势存在时，检验法对不同分布趋势的识别能力大不一样。Kendall检验对EV3型极值分布趋势的识别能力最高，而对LN趋势的识别能力最低。由此说明，当序列存在趋势时，Kendall趋势检验法与序列的分布类型有关，而通常我们错误地认为秩次检验法与分布类型无关。
图5Kendall检验法的识别能力与分布类型的关系
4.2 识别能力与形状参数的关系
为了说明趋势检验法的识别能力与分布类型的形状参数的关系，我们以PE3和EV3分布为例做了研究。除分布的形状参数取值（见表6）不同以外，其它参数取值与4.1中所取值相同。两种检验法的识别能力与形状参数取值的关系见表6，并以图6来表示。
表6Kendall的识别能力与分布的形状参数的关系
|
分布类型
|
形状参数γ
|
趋势度
|
||||||||
|
-0.010
|
-0.008
|
-0.006
|
-0.004
|
0.000
|
0.004
|
0.006
|
0.008
|
0.010
|
||
|
PE3
|
2.0
|
0.820
|
0.648
|
0.440
|
0.251
|
0.045
|
0.245
|
0.445
|
0.652
|
0.799
|
|
1.5
|
0.685
|
0.495
|
0.318
|
0.173
|
0.045
|
0.162
|
0.317
|
0.499
|
0.656
|
|
|
1.2
|
0.613
|
0.434
|
0.274
|
0.154
|
0.045
|
0.142
|
0.273
|
0.433
|
0.586
|
|
|
1.0
|
0.571
|
0.400
|
0.253
|
0.144
|
0.045
|
0.128
|
0.255
|
0.396
|
0.556
|
|
|
EV3
|
1.5
|
0.913
|
0.831
|
0.656
|
0.447
|
0.045
|
0.424
|
0.691
|
0.841
|
0.917
|
|
1.2
|
0.862
|
0.727
|
0.531
|
0.321
|
0.045
|
0.305
|
0.547
|
0.729
|
0.865
|
|
|
0.9
|
0.777
|
0.597
|
0.388
|
0.225
|
0.045
|
0.208
|
0.405
|
0.603
|
0.755
|
|
|
0.6
|
0.631
|
0.446
|
0.273
|
0.158
|
0.045
|
0.146
|
0.285
|
0.447
|
0.610
|
|
|
0.3
|
0.489
|
0.340
|
0.215
|
0.126
|
0.045
|
0.113
|
0.209
|
0.339
|
0.475
|
|
从模拟结果中我们可以找出检验法识别能力与形状参数取值之间的相关规律，那就是同一趋势度下，形状参数值越大，检验法的识别能力越大。因此，即使同一个站点的不同时间序列（例如降雨和径流）有相类似的趋势成份存在，但由于他们服从的分布类型的形状参数不同，采用Kendall去进行趋势检验，也可能会出现不一样的结果。
图6Kendall检验法的识别能力与PE3型分布形状参数的关系
5、应用实例
以云南省昆明市普渡河流域滇池以下的海口、蔡家村和三江口水文站为例，采用Kendall非参数秩次检验法对三站的年径流系列进行趋势检验，径流系列情况及检验结果见表7，其中显著性水平
均取0.002。从表7可以看出，在该显著性水平下，蔡家村和三江口的径流系列存在明显的趋势成分，且这两个站的年径流在实测期内有上升的趋势，这一点从图7中也可以看出，从1997年开始，三江口站的年径流就存在明显的增加。经调查发现，该地区的人类活动频繁，1997年在海口~蔡家村区间兴建一调水工程，将上游滇池水引入该区间，导致蔡家村和三江口站年径流增加。因此，这两站的实测年径流系列不满足一致性条件，在对这两站进行频率分析之前，必须根据前面分析的原因，对其年径流系列进行还原计算，把全部系列建立在同一基础上。
表7海口、蔡家村和三江口站年径流趋势分析成果
|
站名
|
系列
|
系列参数（矩法计算）
|
Kendall检验成果
|
|||
|
均值
|
CV
|
CS
|
是否存在趋势
|
趋势类型
|
||
|
海口
|
1957-2003
|
12.3
|
0.62
|
0.92
|
否
|
/
|
|
蔡家村
|
1953-2003
|
28.0
|
0.52
|
0.71
|
是
|
上升
|
|
三江口
|
1954-2003
|
88.5
|
0.57
|
1.87
|
是
|
上升
|
本文主要研究了Kendall非参数秩次检验法识别水文时间序列趋势成分的能力，大量的模拟研究表明：
（1）Kendall检验法的识别能力与之前给定的显著性水平、容量、趋势度和变差系数有关。其识别能力是趋势度的绝对值、容量和显著性水平的增函数，是变差系数的减函数。
（2）当时间序列存在趋势时，Kendall检验法与时间序列服从的分布类型和形状参数有关，但我们还无法判断识别能力与形状参数的关系。
（3）通过采用Kendall非参数秩次检验法分析了普渡河流域滇池以下的海口、蔡家村和三江口水文站的年径流序列，发现蔡家村和三江口站年径流存在明显的上升趋势，造成径流增加的原因是人类活动（上游调水工程）的影响。
参考文献
[1] 叶守泽, 詹道江主编. 工程水文学[M]. 北京: 中国水利电力出版社, 2000年.
[2] Hirsch, R.M., Slack, J.R. Non-parametric trend test for seasonal data with serial dependence [J]. Water Resource Research, 1984, 20 (6): 727-732.
[3] Gan, T.Y. Hydroclimatic trends and possible climatic warming in the Canadian Prairies [J]. Water Resource Research, 1998, 34 (11): 3009-3015.
[4] Douglas, E.M., Vogel, R.M., Kroll, C.N. Trends in floods and low flows in the United States: impact of spatial correlation [J]. Journal of Hydrology, 2000 (240): 90-105.
[5] 丁晶, 邓育仁. 随机水文学 [M]. 成都: 成都科技大学出版社, 1988年.
[6] Hirsch, R.M., Slack, J.R., Smith, R.A. Techniques of trend analysis for monthly water quality data [J]. Water Resource Research, 1982, 18 (1): 107-121.
[7] 朱勇华, 邰淑彩, 孙韫玉. 应用数理统计 [M]. 武汉: 武汉水利电力大学出版社, 2000年.
