基于Copula函數(shù)的HW一致性檢驗(yàn)修正方法

梁冀雨，林炳章

(南京信息工程大學(xué)應(yīng)用水文氣象研究院，南京 210044)

1 研究背景

相較于傳統(tǒng)的單站頻率分析，地區(qū)頻率分析法通過(guò)整合運(yùn)用某一特定地區(qū)內(nèi)所有站點(diǎn)歷史資料序列的總體平均特征，優(yōu)化對(duì)該區(qū)域內(nèi)各站點(diǎn)的頻率分析結(jié)果，使之更為準(zhǔn)確、穩(wěn)健[1]。在使用地區(qū)頻率分析法時(shí)，需滿足假設(shè)：研究區(qū)域內(nèi)所有站點(diǎn)除了有一個(gè)不同的尺度系數(shù)外，頻率分布的線型和參數(shù)應(yīng)保持一定程度的一致[2]。滿足上述假設(shè)時(shí)，稱該研究區(qū)域?yàn)橐恢聟^(qū)。一致區(qū)的劃分是地區(qū)頻率分析的基礎(chǔ)，區(qū)域的一致性將直接影響頻率估計(jì)值的準(zhǔn)確性[3]。因此,準(zhǔn)確有效的區(qū)域一致性判別方法具有重要的實(shí)踐意義和理論價(jià)值。

在地區(qū)頻率分析法的早期應(yīng)用中，Dalrymple[4]于1960年提出了一種基于樣本離差系數(shù)(Cv)或偏態(tài)系數(shù)(Cs)的區(qū)域一致性檢驗(yàn)方法，被大量應(yīng)用于早期的降雨以及徑流的地區(qū)頻率分析中。但是該法對(duì)經(jīng)驗(yàn)閾值的依賴較大，導(dǎo)致在應(yīng)用于某些特定類型的洪水資料序列時(shí)，對(duì)異質(zhì)站點(diǎn)的分辨能力較差。為了解決這一問(wèn)題，Wiltshire[5]在1984年提出兩種改進(jìn)方案，其一仍然基于樣本離差系數(shù)，其二基于區(qū)域的頻率分布曲線。前者因與分布無(wú)關(guān)，適用范圍更廣。之后，Lu等[6]于1992年利用歸一化的概化極值分布(GEV)的方差特性構(gòu)建了新的一致性檢驗(yàn)方法。在這些一致性檢驗(yàn)方法的基礎(chǔ)上，Hosking和Wallis[7]在1993年提出了至今應(yīng)用最為廣泛的Hosking-Wallis一致性檢驗(yàn)(后文簡(jiǎn)稱HW檢驗(yàn))。HW一致性檢驗(yàn)以樣本線性矩離差系數(shù)(L-Cv)和線性矩偏態(tài)系數(shù)(L-Cs)代替?zhèn)鹘y(tǒng)的離差系數(shù)和偏態(tài)系數(shù)，分析研究區(qū)域內(nèi)各站點(diǎn)資料序列線性矩系數(shù)的離散程度，能夠很好地與地區(qū)線性矩法(RLMA)相結(jié)合，對(duì)異質(zhì)性的鑒別能力較強(qiáng)。

雖然HW一致性檢驗(yàn)至今仍被大量應(yīng)用，但是卻存在理論上的缺陷。Hosking和Wallis在先前的研究中發(fā)現(xiàn)，站點(diǎn)資料間的相關(guān)性會(huì)降低地區(qū)頻率分析中一致性檢驗(yàn)的效力[3],其后2人提出的HW檢驗(yàn)也未針對(duì)該缺陷進(jìn)行調(diào)整。Castellarin等[8]發(fā)現(xiàn)HW檢驗(yàn)的效力也會(huì)隨著資料相關(guān)性的增大而減小，即當(dāng)資料存在相關(guān)性時(shí)，HW檢驗(yàn)可能將存在異質(zhì)站點(diǎn)的區(qū)域判定為一致。Castellarin等給出的解決方案是一個(gè)由站點(diǎn)平均相關(guān)系數(shù)確定的經(jīng)驗(yàn)修正系數(shù)，但該經(jīng)驗(yàn)修正系數(shù)在每次使用時(shí)需事先根據(jù)資料律定參數(shù)，無(wú)法直接應(yīng)用。

本文根據(jù)HW檢驗(yàn)的理論基礎(chǔ)，重點(diǎn)分析資料相關(guān)性影響其效力的原因，并使用蒙特卡羅模擬法揭示資料相關(guān)性對(duì)HW檢驗(yàn)的影響。同時(shí)利用江西省年最大降雨量序列(AMS)研究實(shí)際資料的相關(guān)性特征，通過(guò)結(jié)合資料和理論，給出一種不受資料相關(guān)性影響的一致性檢驗(yàn)。最后再次通過(guò)蒙特卡羅模擬對(duì)新的一致性檢驗(yàn)方法進(jìn)行評(píng)估。結(jié)果顯示，該法可以為地區(qū)頻率分析法，特別是地區(qū)線性矩法提供更有效的理論支持。

2 Hosking-Wallis一致性檢驗(yàn)

為了將線性矩[9]應(yīng)用于地區(qū)頻率分析，Hosking和Wallis在1993年提出了基于資料線性矩系數(shù)的HW檢驗(yàn)。HW檢驗(yàn)將研究區(qū)域內(nèi)各站點(diǎn)資料的各階線性矩系數(shù)的離散程度，與符合一致區(qū)定義的、與研究區(qū)域具有相同站點(diǎn)數(shù)以及資料長(zhǎng)度的人工模擬資料相比較。2者的差別若在一定的范圍之內(nèi)，就認(rèn)定研究區(qū)域?yàn)橐恢聟^(qū)。

(1)

于是,區(qū)域內(nèi)各站點(diǎn)的平均加權(quán)標(biāo)準(zhǔn)差V為：

(2)

最后進(jìn)行蒙特卡羅模擬，生成500組模擬區(qū)域數(shù)據(jù)，根據(jù)式(2)計(jì)算得到各組模擬數(shù)據(jù)的平均加權(quán)標(biāo)準(zhǔn)差Vsim。利用Vsim的均值μV和標(biāo)準(zhǔn)差δV，對(duì)研究區(qū)域的平均加權(quán)標(biāo)準(zhǔn)差V進(jìn)行類正態(tài)化處理，得到HW檢驗(yàn)統(tǒng)計(jì)量H：

(3)

為了得到μV、σV，需進(jìn)行蒙特卡洛模擬。首先根據(jù)式(1)，使用研究區(qū)域各站點(diǎn)資料計(jì)算得到的前4階地區(qū)平均線性矩系數(shù)，擬合4參數(shù)的kappa分布[10]。接著模擬生成500個(gè)人工區(qū)域，每個(gè)人工區(qū)域都擁有和研究區(qū)域相同的站點(diǎn)數(shù)、資料長(zhǎng)度，同時(shí)符合先前得到的kappa分布。值得注意的是，由此得到的每個(gè)人工區(qū)域都符合一致區(qū)的假設(shè)，并且站點(diǎn)間不存在相關(guān)關(guān)系。最后，對(duì)每一個(gè)人工區(qū)域分別計(jì)算平均加權(quán)標(biāo)準(zhǔn)差Vsim，并得到500個(gè)Vsim的均值和標(biāo)準(zhǔn)差。

Hosking和Wallis認(rèn)為[9]，當(dāng)H小于1時(shí)，研究區(qū)域?yàn)榭山邮艿囊恢聟^(qū)；當(dāng)H大于1并小于2時(shí)，研究區(qū)為可能的異質(zhì)區(qū)域；而當(dāng)H大于2時(shí)，研究區(qū)域?yàn)榇_定的異質(zhì)區(qū)域。

另外，在考慮線性偏態(tài)系數(shù)和線性峰度系數(shù)的情況下，V具有2種變式，即：

(4)

式中：V2、V3分別代表同時(shí)考慮線性離差系數(shù)、線性偏態(tài)系數(shù)，以及線性偏態(tài)系數(shù)、線性峰度系數(shù)的情況。

此時(shí)，為區(qū)別3者不同，式(2)中的V改稱V1。而由3者得到的H對(duì)應(yīng)記作H1、H2和H3。事實(shí)上，3者結(jié)構(gòu)相似，計(jì)算得到的H值往往代表相同的結(jié)果，因此Hosking等[9]建議多數(shù)情況下可僅考慮H1的結(jié)果。其后，Viglione等[11]發(fā)現(xiàn)，H1相較于H2以及H3，具有更強(qiáng)的鑒別異質(zhì)性(heterogeneity)能力，且表現(xiàn)更為穩(wěn)定。綜上，本文僅對(duì)H1進(jìn)行分析，相應(yīng)的統(tǒng)計(jì)量，后文簡(jiǎn)稱H以及V。

3 實(shí)測(cè)與生成雨量資料的相關(guān)性

相關(guān)系數(shù)是表征2個(gè)資料序列間相關(guān)程度的統(tǒng)計(jì)量。其中最為常用的是統(tǒng)計(jì)學(xué)家K. Pearson于1895年首次提出的皮爾森積矩相關(guān)系數(shù)，簡(jiǎn)稱相關(guān)系數(shù)。對(duì)于2個(gè)成對(duì)的隨機(jī)變量序列X、Y，其皮爾森相關(guān)系數(shù)的表達(dá)式為：

(5)

式中：rX,Y、cov(X,Y)分別表示隨機(jī)序列X、Y的皮爾森相關(guān)系數(shù)和協(xié)方差。

由式(5)可知，皮爾森相關(guān)系數(shù)可以被看作2個(gè)隨機(jī)變量之間標(biāo)準(zhǔn)化的協(xié)方差，當(dāng)隨機(jī)序列X、Y為標(biāo)準(zhǔn)正態(tài)分布時(shí)，2者的標(biāo)準(zhǔn)差都為1，式(5)右邊可簡(jiǎn)化為cov(X,Y)。因此，選擇皮爾森相關(guān)系數(shù)作為相關(guān)性的度量能夠在進(jìn)行后文的蒙特卡羅模擬實(shí)驗(yàn)時(shí)，通過(guò)分解協(xié)方差矩陣的方法[12]，生成具有特定相關(guān)性的資料序列。

為了比較HW檢驗(yàn)中生成的模擬年最大降雨量資料與實(shí)測(cè)資料相關(guān)性的區(qū)別，本文以與HW檢驗(yàn)相同的設(shè)置[9]進(jìn)行蒙特卡羅模擬實(shí)驗(yàn)。從15 a到50 a，每一個(gè)序列長(zhǎng)度生成1 000個(gè)人造區(qū)域，每個(gè)人造區(qū)域都包含100個(gè)站點(diǎn)。為做到與實(shí)測(cè)資料一致，令每個(gè)站點(diǎn)的模擬資料都服從4參數(shù)的kappa分布[10]，并擁有江西全省平均的24 h年最大雨量資料線性矩系數(shù)：t=0.18，t3=0.20，t4=0.15。同時(shí)因?qū)Y料進(jìn)行歸一化處理，所以所有生成資料的均值為1。之后，求出每個(gè)人工區(qū)域的站點(diǎn)間平均相關(guān)系數(shù)，最后計(jì)算每個(gè)序列長(zhǎng)度1 000個(gè)人工區(qū)域平均相關(guān)系數(shù)的均值，繪于圖1。

圖1 實(shí)測(cè)與模擬雨量序列平均相關(guān)系數(shù)隨序列長(zhǎng)度的分布

圖1中，圓形標(biāo)記為實(shí)測(cè)年最大雨量平均相關(guān)系數(shù)隨資料長(zhǎng)度的分布，三角標(biāo)記為HW檢驗(yàn)中使用的模擬資料的情況。根據(jù)圖1，模擬相較實(shí)測(cè)資料折線圖下降趨勢(shì)更為明顯、平滑，且始終位于實(shí)測(cè)資料的下方，表明實(shí)測(cè)資料中普遍存在的相關(guān)性大于HW檢驗(yàn)中生成的與之比較的模擬數(shù)據(jù)，且相差達(dá)到25%以上。

4 資料相關(guān)性對(duì)Hosking-Wallis一致性檢驗(yàn)結(jié)果的影響

4.1 蒙特卡羅模擬中年最大雨量相關(guān)性對(duì)H統(tǒng)計(jì)量的影響

為了定量研究年最大雨量資料的相關(guān)性對(duì)H統(tǒng)計(jì)量的影

圖2 資料相關(guān)性對(duì)Hosking-Wallis一致性檢驗(yàn)的影響

由圖2可知，相關(guān)性對(duì)不同的地區(qū)分布的影響程度是大致相同的，資料間的平均相關(guān)系數(shù)越大，H統(tǒng)計(jì)量曲線的位置越靠下，其值越小。當(dāng)區(qū)域平均相關(guān)系數(shù)從0.2增加到0.8時(shí)，H值減小量均達(dá)到1.5以上，已經(jīng)超過(guò)了判定有效一致區(qū)的H值區(qū)間[0,1]的長(zhǎng)度，這大大降低了HW檢驗(yàn)對(duì)資料相關(guān)區(qū)域異質(zhì)性的辨別能力。

4.2 不同站點(diǎn)間相關(guān)性空間分布情況對(duì)H統(tǒng)計(jì)量的影響

(5)

式中：R為資料的相關(guān)系數(shù)矩陣，其形式?jīng)Q定了生成的資料間相關(guān)系數(shù)的均一的。

(6)

圖3展示了地區(qū)分布為GEV，區(qū)域站點(diǎn)數(shù)為66時(shí)，不同相關(guān)系數(shù)分布情況對(duì)H值的影響。圖3的橫縱坐標(biāo)與圖2相同，改進(jìn)后的模擬結(jié)果以圓點(diǎn)虛線分別繪于2張圖上。圖3(a)中黑色實(shí)線和三角形點(diǎn)虛線分別表示以式(5)的形式生成的平均相關(guān)系數(shù)為0.4、0.8的10 000個(gè)模擬區(qū)域的HW檢驗(yàn)結(jié)果的均值，與圖2(e)相同；而圓點(diǎn)虛線代表的人造區(qū)域中，僅有1/2的站點(diǎn)間具有0.8的相關(guān)系數(shù)，根據(jù)式(6)，其區(qū)域平均相關(guān)系數(shù)為0.4。圖3(b)與圖3(a)類似，其圓點(diǎn)虛線代表的人造區(qū)域中 的站點(diǎn)間擁有0.6的相關(guān)系數(shù)，區(qū)域平均相關(guān)系數(shù)為0.2。圖3(a)、(b)中，圓點(diǎn)虛線幾乎與黑色實(shí)線重合，而與三角形點(diǎn)虛線相距較遠(yuǎn)。由此不難發(fā)現(xiàn)，在區(qū)域平均相關(guān)系數(shù)相同時(shí)，相關(guān)性對(duì)H值的影響是近似相同的，而與站點(diǎn)間相關(guān)性空間分布無(wú)關(guān)。

基本指令是可編程控制編程中的重要組成部分，有 LD、LD NOT、AND、AND NOT、OR 等?；局噶疃嘤糜诤?jiǎn)單控制，編程較為靈活；對(duì)于繁瑣的控制，若僅使用基本指令，編程較復(fù)雜，可用功能指令配合編程使其程序簡(jiǎn)化，且所編的程序較易閱讀。

圖3 不同相關(guān)系數(shù)分布對(duì)H值的影響

5 對(duì)Hosking-Wallis一致性檢驗(yàn)的修正

5.1 修正Hosking-Wallis檢驗(yàn)的算法及步驟

由第4節(jié)知,HW檢驗(yàn)的偏差主要來(lái)自研究區(qū)域資料間固有的相關(guān)性與檢驗(yàn)中生成的人工資料序列間相關(guān)性的不同。而即使站點(diǎn)資料序列間相關(guān)系數(shù)分布情況不同，只要整個(gè)區(qū)域的平均相關(guān)系數(shù)相同，其對(duì)HW檢驗(yàn)結(jié)果的影響就基本相同?；谝陨辖Y(jié)論，對(duì)HW檢驗(yàn)進(jìn)行修正，使檢驗(yàn)中生成的人工資料擁有與研究區(qū)域相同的平均相關(guān)系數(shù)，改進(jìn)后的HW檢驗(yàn)步驟如下。

(7)

式中：ρij為研究區(qū)域相關(guān)系數(shù)矩陣第i行第j列的元素；N為區(qū)域站點(diǎn)數(shù)。

(2)由式(5)得到N維平均相關(guān)系數(shù)矩陣R。

(3)計(jì)算研究區(qū)域各站點(diǎn)的各階線性系數(shù)和區(qū)域的各階平均線性離差系數(shù)。

(4)根據(jù)式(2)計(jì)算統(tǒng)計(jì)量V。

(5)生成與研究區(qū)域相同站點(diǎn)數(shù)、資料長(zhǎng)度的人工區(qū)域資料yik，i=1,2,…,N，k=1,2,…,ni。ni為第i個(gè)站點(diǎn)的資料長(zhǎng)度。運(yùn)用正態(tài)Copula函數(shù)，使人工區(qū)域中站點(diǎn)資料服從多元正態(tài)分布，擁有均值零，協(xié)方差矩陣R。

(6)由(3)中求得的地區(qū)平均線性矩系數(shù)擬合4參數(shù)kappa分布。將資料yik轉(zhuǎn)化為符合擬合的kappa分布的資料Qik：

Qik=KAP[φ(yik)]

(8)

式中：KAP為擬合得到的kappa函數(shù)的分位數(shù)函數(shù)；Φ為標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)。

(7)對(duì)人工區(qū)域資料序列Qik計(jì)算Vsim。

(8)重復(fù)1 000次步驟(5)～(7)，計(jì)算得到的1 000個(gè)Vsim的均值μV和標(biāo)準(zhǔn)差σV，并根據(jù)式(3)計(jì)算得到修正后的統(tǒng)計(jì)量H*。

5.2 資料相關(guān)性對(duì)修正前后Hosking-Wallis檢驗(yàn)的影響

圖4 改進(jìn)前后HW檢驗(yàn)結(jié)果對(duì)比

6 結(jié) 語(yǔ)

(1)利用蒙特卡洛模擬實(shí)驗(yàn)，得出HW一致性檢驗(yàn)中生成的人工資料序列間的平均相關(guān)系數(shù)隨資料長(zhǎng)度的分布，并將其與實(shí)測(cè)資料進(jìn)行對(duì)比。結(jié)果表明人工資料間的相關(guān)性均大于實(shí)測(cè)資料25%以上。

(2)從HW檢驗(yàn)的定義出發(fā)，定性分析資料相關(guān)性影響H統(tǒng)計(jì)量值的原因，并通過(guò)蒙特卡羅模擬實(shí)驗(yàn)，定量分析資料相關(guān)性以及不同的相關(guān)性分布對(duì)H值大小的影響。定量分析發(fā)現(xiàn)當(dāng)資料間存在較大相關(guān)性時(shí)，HW檢驗(yàn)對(duì)異質(zhì)性的檢驗(yàn)?zāi)芰Υ蠓陆担?dāng)區(qū)域平均相關(guān)系數(shù)相同時(shí)，相關(guān)性導(dǎo)致H值的減小量相同。

(3)針對(duì)雨量資料間的平均相關(guān)系數(shù)，對(duì)HW檢驗(yàn)進(jìn)行修正。實(shí)踐證實(shí)，資料相關(guān)性對(duì)修正后HW檢驗(yàn)結(jié)果的影響大大降低，同時(shí)修正后的HW檢驗(yàn)依舊保留了其對(duì)區(qū)域一致性的鑒別能力。這種新的一致性檢驗(yàn)方法適用條件等同HW檢驗(yàn)，無(wú)需額外調(diào)整，計(jì)算速度快，可直接運(yùn)用于地區(qū)線性矩頻率法，幫助劃分水文一致區(qū)。相較傳統(tǒng)HW檢驗(yàn)更為符合實(shí)測(cè)雨量資料情況、檢驗(yàn)結(jié)果更加可靠。

(4)修正后的一致性檢驗(yàn)通過(guò)優(yōu)化水文一致區(qū)劃分，提升了針對(duì)各水文時(shí)間序列資料，特別是降雨、徑流序列頻率分析的效率及其結(jié)果的準(zhǔn)確性。將改進(jìn)后的降雨、徑流頻率分析結(jié)果與氣象降水量預(yù)報(bào)或流域水文模型徑流量預(yù)報(bào)相比較，能夠更為準(zhǔn)確地劃定預(yù)報(bào)降水或洪水的量級(jí)、重現(xiàn)期，為洪水預(yù)警預(yù)報(bào)提供更加堅(jiān)實(shí)的理論基礎(chǔ)。

□