Kendall檢驗在水文序列趨勢分析中的研究

申請免費試用、咨詢電話：400-8352-114

簡介： 本文通過蒙特卡洛模擬對Kendall非參數(shù)秩次相關(guān)檢驗法在水文時間序列趨勢檢驗中的識別能力做了研究。模擬結(jié)果表明，該方法的識別能力與之前給定的顯著性水平、樣本容量、趨勢度和變差系數(shù)有關(guān)。隨著趨勢度的絕對值、樣本容量和顯著性水平的增加，Kendall檢驗法的識別能力增強；隨著變差系數(shù)的增加，該方法的識別能力減?。划?dāng)趨勢存在時，該方法還與時間序列服從的分布類型和形狀參數(shù)有關(guān)。最后采用Kendall秩次檢驗法分析了普渡河流域?？?、蔡家村和三江口水文站的年徑流系列，發(fā)現(xiàn)蔡家村和三江口的年徑流有增加的趨勢，找出了增加的原因。
關(guān)鍵字：Mann-Kendall 非參數(shù)秩次相關(guān)檢驗法 蒙塔卡洛模擬 水文時間序列 趨勢分析

1、引言

　　水文序列是從工程所在地點或鄰近地點水文觀測的資料中選取表征水文過程特征值（如洪峰流量或水位、各種時段洪水總量等）的樣本。它們是進行頻率分析、估計設(shè)計水文過程的基礎(chǔ)資料。在水文資料的觀測期內(nèi)，如因流域上修建了蓄水、引水、水土保持等工程以及流域的氣候發(fā)生變化等等，這些人工或天然的原因使流域水文現(xiàn)象的形成條件發(fā)生了顯著的改變，因而水文變量的概率分布規(guī)律也發(fā)生了顯著的變異，我們把這一問題稱為水文資料的“非一致性”問題[1]。如果將“非一致性”的水文資料混雜在一起作為一個樣本進行水文頻率計算顯然違背了用于頻率分析的水文序列必須服從同一分布的要求。

　　因此，進行頻率分析之前，必須檢驗序列是否具有一致性，其中趨勢檢驗是一致性檢驗的內(nèi)容之一。在水文資料分析研究中，非參數(shù)秩次相關(guān)統(tǒng)計檢驗，即Kendall檢驗用于識別時間序列是否存在趨勢成分。本文首先簡單介紹了MK趨勢檢驗法；研究了該方法識別正態(tài)分布序列趨勢的能力與給定的顯著性水平、樣本容量、趨勢度和變差系數(shù)的關(guān)系；同時探討該方法識別非正態(tài)分布序列趨勢的能力與分布類型和形狀參數(shù)的關(guān)系；最后舉出該方法的一個應(yīng)用實例。

2、Kendall秩次相關(guān)檢驗

　　Kendall非參數(shù)秩次相關(guān)檢驗法已經(jīng)廣泛的用于檢驗水文氣象資料的趨勢成分，包括水質(zhì)、流量、氣溫和降雨序列等，究其原因主要是，與參數(shù)統(tǒng)計檢驗法相比，非參數(shù)檢驗法更適用于非正態(tài)分布或經(jīng)過刪檢（刪去低于或高于某水平的觀測值）的資料，而這些情況在時間序列分析中常常會遇到。過去20年里，國際上關(guān)于MK方法應(yīng)用研究的實例非常之多[2-4]。盡管該方法應(yīng)用如此之廣泛，但我們還是不清楚該方法是不是適用于各種情況下的時間序列的趨勢檢驗。

　　對序列，先確定所有對偶值中與的大小關(guān)系（設(shè)為）。趨勢檢驗的統(tǒng)計量[5]為：

　?。?）

　　式中：

　?。唬?）

　?。?）

　　當(dāng)大于10時，收斂于標(biāo)準(zhǔn)正態(tài)分布。

　　原假設(shè)為該序列無趨勢，采用雙邊趨勢檢驗，在給定顯著性水平下，在正態(tài)分布表中查得臨界值，當(dāng)時，接受原假設(shè)，即趨勢不顯著；若，則拒絕原假設(shè)，即認(rèn)為趨勢顯著。

　　3、Kendall檢驗對正態(tài)分布序列趨勢的識別能力

　　3.1 兩類錯誤與識別能力

　　假設(shè)檢驗[7]中有兩類錯誤：第一類錯誤為“棄真”，即原假設(shè)本來是正確的，但檢驗結(jié)果卻拒絕了原假設(shè)，由數(shù)理統(tǒng)計知識可知，犯第一類錯誤的概率為顯著性水平；第二類錯誤為“納偽”，即原假設(shè)本來是錯誤的，但檢驗結(jié)果卻接受了原假設(shè)。犯第二類錯誤的概率為。

　　定義兩種檢驗方法的識別能力為，也就是當(dāng)原假設(shè)錯誤時，該檢驗方法正確拒絕原假設(shè)的能力。

3.2 蒙特卡洛模擬試驗步驟

　　因為水文時間序列可看成是多種成分組成。假定這些成分是線性疊加，可按下式表示[5]：

　　式中為確定性的非周期成分（包括趨勢，跳躍等暫態(tài)成分等）；為確定性的周期成分（包括簡單的或復(fù)合周期的成分等）；為隨機成分（包括平穩(wěn)的或非平穩(wěn)的隨機成分）。本次研究的時間序列只考慮確定性的非周期成分中的趨勢項和隨機成分。

　　蒙特卡洛模擬用于評價Kendall檢驗識別趨勢成分的能力，試驗步驟如下：

　?、倭?；

　?、谏烧龖B(tài)分布的隨機序列，樣本容量可取，樣本均值均取1.0，方差分別為，其中；樣本的標(biāo)準(zhǔn)偏差和變差系數(shù)為；

　?、凵哨厔莩煞?，，為趨勢度，可取-0.01，-0.005，0，0.005和0.01等；

　　④將②和③中生成的序列按序號對應(yīng)相加得到有趨勢的時間序列；

　?、菰僭O(shè)為序列不存在趨勢，采用Kendall檢驗法對序列進行趨勢檢驗。如果拒絕原假設(shè)，則，否則；

　?、拗貜?fù)步驟②~⑤次；

　?、吒鶕?jù)前面的定義，Kendall檢驗的識別能力為：

　　式中：為蒙特卡洛模擬的次數(shù)；為檢驗落在拒絕域里的次數(shù)。

　　3.3 結(jié)果分析

　　3.3.1 識別能力與顯著性水平、趨勢度的關(guān)系

　　樣本容量取50；均值取1.0；變差系數(shù)取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趨勢檢驗法的識別能力與顯著性水平和趨勢度的關(guān)系見表1和圖1。從蒙特卡洛模擬結(jié)果可以看出，對于同一顯著性水平，趨勢度絕對值越大，則檢驗法的識別能力越強；而對于同一趨勢度，顯著性水平越大，則檢驗法的識別能力越強；特別的，當(dāng)時，序列無趨勢存在，原假設(shè)是正確的，此時，檢驗法的識別能力均等于事先給定的顯著性水平。

　　表1Kendall的識別能力與顯著性水平和趨勢度的關(guān)系

趨勢度顯著性水平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的識別能力與顯著性水平和趨勢度的關(guān)系

3.3.2 識別能力與樣本容量、趨勢度的關(guān)系

　　顯著性水平取0.05；均值取1.0；變差系數(shù)取0.5；樣本容量取10~100；趨勢度取-0.01，-0.005，0，0.005和0.01。Kendall趨勢檢驗法的識別能力與樣本容量和趨勢度的關(guān)系見表2和圖2。從蒙特卡洛模擬結(jié)果可以看出，隨著樣本容量和趨勢度絕對值的增加，檢驗法的識別能力增強；特別的，當(dāng)時，序列無趨勢存在，原假設(shè)是正確的，此時，檢驗法的識別能力均等于事先給定的顯著性水平。

　　表2MK的識別能力與樣本容量和趨勢度的關(guān)系

趨勢度樣本容量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的識別能力與樣本容量和趨勢度的關(guān)系

3.3.3 識別能力與趨勢度、變差系數(shù)的關(guān)系

　　樣本容量取50；顯著性水平取0.05；均值為1.0；變差系數(shù)取0.01~1.0；趨勢度取-0.01，-0.005，0，0.005和0.01。Kendall趨勢檢驗法的識別能力與變差系數(shù)、趨勢度的關(guān)系見表3和圖3。從蒙特卡洛模擬結(jié)果可以看出，對于同一變差系數(shù)，趨勢度絕對值越大，則檢驗法的識別能力越強；而對于同一趨勢度，變差系數(shù)越大，則檢驗法越難識別出序列存在趨勢成分；特別的，當(dāng)時，序列無趨勢存在，原假設(shè)是正確的，此時，檢驗法的識別能力均等于事先給定的顯著性水平。

　　表3Kendall的識別能力與變差系數(shù)和趨勢度的關(guān)系

趨勢度變差系數(shù)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的識別能力與變差系數(shù)和趨勢度的關(guān)系

3.3.4 識別能力與樣本容量、變差系數(shù)的關(guān)系

　　顯著性水平取0.05；均值為1.0；趨勢度取0.005；樣本容量取10~100；變差系數(shù)取0.01~1.0。兩種趨勢檢驗法的識別能力與樣本容量和變差系數(shù)的關(guān)系見表4和圖4。從蒙特卡洛模擬結(jié)果可以看出，對于同一變差系數(shù)，檢驗法的識別能力隨著樣本容量的增加而增強；而對于同一樣本容量，檢驗法的識別能力則隨著變差系數(shù)的增加而減小；特別的，當(dāng)時，序列無趨勢存在，原假設(shè)是正確的，此時，檢驗法的識別能力均等于事先給定的顯著性水平。

　　表4Kendall的識別能力與樣本容量和變差系數(shù)的關(guān)系

樣本容量方差系數(shù)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的識別能力與樣本容量和變差系數(shù)的關(guān)系

4、Kendall檢驗對非正態(tài)分布序列趨勢的識別能力

　　前面已經(jīng)討論了Kendall檢驗對正態(tài)分布序列趨勢的識別能力，然而，水文時間序列通常是有偏的，很少有服從正態(tài)分布的，因此，下面進一步探討Kendall檢驗法識別非正態(tài)分布序列趨勢的能力。這里比較的都是水文中常用的幾種分布類型，包括三類極值分布（EV1、EV2和EV3）、對數(shù)正態(tài)分布（LN）和皮爾遜Ⅲ型分布（PE3）。

4.1 識別能力與分布類型的關(guān)系

　　為了方便比較研究，樣本容量取為50，均值和方差分別取為1.0和0.5，形狀參數(shù)取值見表5；趨勢度仍取-0.01，-0.005，0，0.005和0.01等；蒙特卡洛模擬次數(shù)取2000；顯著性水平取0.05。根據(jù)3.2的步驟評價兩種檢驗法的識別能力，檢驗結(jié)果列于表5中，并以圖5來表示。

　　表5 Kendall的識別能力與分布類型的關(guān)系

分布類型趨勢度-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

　　從模擬結(jié)果可以看出，當(dāng)，即沒有趨勢存在時，檢驗法對所有分布的識別能力等于之前給定的顯著性水平0.05，說明檢驗統(tǒng)計量對序列服從的分布類型并不敏感。然而，當(dāng)時，即有趨勢存在時，檢驗法對不同分布趨勢的識別能力大不一樣。Kendall檢驗對EV3型極值分布趨勢的識別能力最高，而對LN趨勢的識別能力最低。由此說明，當(dāng)序列存在趨勢時，Kendall趨勢檢驗法與序列的分布類型有關(guān)，而通常我們錯誤地認(rèn)為秩次檢驗法與分布類型無關(guān)。

　　圖5Kendall檢驗法的識別能力與分布類型的關(guān)系

4.2 識別能力與形狀參數(shù)的關(guān)系

　　為了說明趨勢檢驗法的識別能力與分布類型的形狀參數(shù)的關(guān)系，我們以PE3和EV3分布為例做了研究。除分布的形狀參數(shù)取值（見表6）不同以外，其它參數(shù)取值與4.1中所取值相同。兩種檢驗法的識別能力與形狀參數(shù)取值的關(guān)系見表6，并以圖6來表示。

　　表6Kendall的識別能力與分布的形狀參數(shù)的關(guān)系

分布類型形狀參數(shù)γ趨勢度-0.010 -0.008 -0.006 -0.004 0.000 0.004 0.006 0.008 0.010 PE32.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 EV31.50.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

　　從模擬結(jié)果中我們可以找出檢驗法識別能力與形狀參數(shù)取值之間的相關(guān)規(guī)律，那就是同一趨勢度下，形狀參數(shù)值越大，檢驗法的識別能力越大。因此，即使同一個站點的不同時間序列（例如降雨和徑流）有相類似的趨勢成份存在，但由于他們服從的分布類型的形狀參數(shù)不同，采用Kendall去進行趨勢檢驗，也可能會出現(xiàn)不一樣的結(jié)果。

　　圖6Kendall檢驗法的識別能力與PE3型分布形狀參數(shù)的關(guān)系

5、應(yīng)用實例

　　以云南省昆明市普渡河流域滇池以下的?？?、蔡家村和三江口水文站為例，采用Kendall非參數(shù)秩次檢驗法對三站的年徑流系列進行趨勢檢驗，徑流系列情況及檢驗結(jié)果見表7，其中顯著性水平均取0.002。從表7可以看出，在該顯著性水平下，蔡家村和三江口的徑流系列存在明顯的趨勢成分，且這兩個站的年徑流在實測期內(nèi)有上升的趨勢，這一點從圖7中也可以看出，從1997年開始，三江口站的年徑流就存在明顯的增加。經(jīng)調(diào)查發(fā)現(xiàn)，該地區(qū)的人類活動頻繁，1997年在?？趡蔡家村區(qū)間興建一調(diào)水工程，將上游滇池水引入該區(qū)間，導(dǎo)致蔡家村和三江口站年徑流增加。因此，這兩站的實測年徑流系列不滿足一致性條件，在對這兩站進行頻率分析之前，必須根據(jù)前面分析的原因，對其年徑流系列進行還原計算，把全部系列建立在同一基礎(chǔ)上。

　　表7?？凇⒉碳掖搴腿谡灸陱搅髭厔莘治龀晒?/p>

站名系列系列參數(shù)（矩法計算）Kendall檢驗成果均值CVCS是否存在趨勢趨勢類型?？?957-200312.30.620.92否/蔡家村1953-200328.00.520.71是　上升三江口1954-200388.50.571.87是上升6、結(jié)語

　　本文主要研究了Kendall非參數(shù)秩次檢驗法識別水文時間序列趨勢成分的能力，大量的模擬研究表明：

　?。?）Kendall檢驗法的識別能力與之前給定的顯著性水平、樣本容量、趨勢度和變差系數(shù)有關(guān)。其識別能力是趨勢度的絕對值、樣本容量和顯著性水平的增函數(shù)，是變差系數(shù)的減函數(shù)。

　?。?）當(dāng)時間序列存在趨勢時，Kendall檢驗法與時間序列服從的分布類型和形狀參數(shù)有關(guān)，但我們還無法判斷識別能力與形狀參數(shù)的關(guān)系。

　?。?）通過采用Kendall非參數(shù)秩次檢驗法分析了普渡河流域滇池以下的海口、蔡家村和三江口水文站的年徑流序列，發(fā)現(xiàn)蔡家村和三江口站年徑流存在明顯的上升趨勢，造成徑流增加的原因是人類活動（上游調(diào)水工程）的影響。

參考文獻

　　[1] 葉守澤, 詹道江主編. 工程水文學(xué)[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] 丁晶, 鄧育仁. 隨機水文學(xué) [M]. 成都: 成都科技大學(xué)出版社, 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] 朱勇華, 邰淑彩, 孫韞玉. 應(yīng)用數(shù)理統(tǒng)計 [M]. 武漢: 武漢水利電力大學(xué)出版社, 2000年.