基于Copula函數的HW一致性檢驗修正方法

梁冀雨，林炳章

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

1 研究背景

相較于傳統的單站頻率分析，地區頻率分析法通過整合運用某一特定地區內所有站點歷史資料序列的總體平均特征，優化對該區域內各站點的頻率分析結果，使之更為準確、穩健[1]。在使用地區頻率分析法時，需滿足假設：研究區域內所有站點除了有一個不同的尺度系數外，頻率分布的線型和參數應保持一定程度的一致[2]。滿足上述假設時，稱該研究區域為一致區。一致區的劃分是地區頻率分析的基礎，區域的一致性將直接影響頻率估計值的準確性[3]。因此,準確有效的區域一致性判別方法具有重要的實踐意義和理論價值。

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

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

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

2 Hosking-Wallis一致性檢驗

為了將線性矩[9]應用于地區頻率分析，Hosking和Wallis在1993年提出了基于資料線性矩系數的HW檢驗。HW檢驗將研究區域內各站點資料的各階線性矩系數的離散程度，與符合一致區定義的、與研究區域具有相同站點數以及資料長度的人工模擬資料相比較。2者的差別若在一定的范圍之內，就認定研究區域為一致區。

(1)

于是,區域內各站點的平均加權標準差V為：

(2)

最后進行蒙特卡羅模擬，生成500組模擬區域數據，根據式(2)計算得到各組模擬數據的平均加權標準差Vsim。利用Vsim的均值μV和標準差δV，對研究區域的平均加權標準差V進行類正態化處理，得到HW檢驗統計量H：

(3)

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

Hosking和Wallis認為[9]，當H小于1時，研究區域為可接受的一致區；當H大于1并小于2時，研究區為可能的異質區域；而當H大于2時，研究區域為確定的異質區域。

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

(4)

式中：V2、V3分別代表同時考慮線性離差系數、線性偏態系數，以及線性偏態系數、線性峰度系數的情況。

此時，為區別3者不同，式(2)中的V改稱V1。而由3者得到的H對應記作H1、H2和H3。事實上，3者結構相似，計算得到的H值往往代表相同的結果，因此Hosking等[9]建議多數情況下可僅考慮H1的結果。其后，Viglione等[11]發現，H1相較于H2以及H3，具有更強的鑒別異質性(heterogeneity)能力，且表現更為穩定。綜上，本文僅對H1進行分析，相應的統計量，后文簡稱H以及V。

3 實測與生成雨量資料的相關性

相關系數是表征2個資料序列間相關程度的統計量。其中最為常用的是統計學家K. Pearson于1895年首次提出的皮爾森積矩相關系數，簡稱相關系數。對于2個成對的隨機變量序列X、Y，其皮爾森相關系數的表達式為：

(5)

式中：rX,Y、cov(X,Y)分別表示隨機序列X、Y的皮爾森相關系數和協方差。

由式(5)可知，皮爾森相關系數可以被看作2個隨機變量之間標準化的協方差，當隨機序列X、Y為標準正態分布時，2者的標準差都為1，式(5)右邊可簡化為cov(X,Y)。因此，選擇皮爾森相關系數作為相關性的度量能夠在進行后文的蒙特卡羅模擬實驗時，通過分解協方差矩陣的方法[12]，生成具有特定相關性的資料序列。

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

圖1 實測與模擬雨量序列平均相關系數隨序列長度的分布

圖1中，圓形標記為實測年最大雨量平均相關系數隨資料長度的分布，三角標記為HW檢驗中使用的模擬資料的情況。根據圖1，模擬相較實測資料折線圖下降趨勢更為明顯、平滑，且始終位于實測資料的下方，表明實測資料中普遍存在的相關性大于HW檢驗中生成的與之比較的模擬數據，且相差達到25%以上。

4 資料相關性對Hosking-Wallis一致性檢驗結果的影響

4.1 蒙特卡羅模擬中年最大雨量相關性對H統計量的影響

為了定量研究年最大雨量資料的相關性對H統計量的影

圖2 資料相關性對Hosking-Wallis一致性檢驗的影響

由圖2可知，相關性對不同的地區分布的影響程度是大致相同的，資料間的平均相關系數越大，H統計量曲線的位置越靠下，其值越小。當區域平均相關系數從0.2增加到0.8時，H值減小量均達到1.5以上，已經超過了判定有效一致區的H值區間[0,1]的長度，這大大降低了HW檢驗對資料相關區域異質性的辨別能力。

4.2 不同站點間相關性空間分布情況對H統計量的影響

(5)

式中：R為資料的相關系數矩陣，其形式決定了生成的資料間相關系數的均一的。

(6)

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

基本指令是可編程控制編程中的重要組成部分，有 LD、LD NOT、AND、AND NOT、OR 等。基本指令多用于簡單控制，編程較為靈活；對于繁瑣的控制，若僅使用基本指令，編程較復雜，可用功能指令配合編程使其程序簡化，且所編的程序較易閱讀。

圖3 不同相關系數分布對H值的影響

5 對Hosking-Wallis一致性檢驗的修正

5.1 修正Hosking-Wallis檢驗的算法及步驟

由第4節知,HW檢驗的偏差主要來自研究區域資料間固有的相關性與檢驗中生成的人工資料序列間相關性的不同。而即使站點資料序列間相關系數分布情況不同，只要整個區域的平均相關系數相同，其對HW檢驗結果的影響就基本相同。基于以上結論，對HW檢驗進行修正，使檢驗中生成的人工資料擁有與研究區域相同的平均相關系數，改進后的HW檢驗步驟如下。

(7)

式中：ρij為研究區域相關系數矩陣第i行第j列的元素；N為區域站點數。

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

(3)計算研究區域各站點的各階線性系數和區域的各階平均線性離差系數。

(4)根據式(2)計算統計量V。

(5)生成與研究區域相同站點數、資料長度的人工區域資料yik，i=1,2,…,N，k=1,2,…,ni。ni為第i個站點的資料長度。運用正態Copula函數，使人工區域中站點資料服從多元正態分布，擁有均值零，協方差矩陣R。

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

Qik=KAP[φ(yik)]

(8)

式中：KAP為擬合得到的kappa函數的分位數函數；Φ為標準正態分布的累積分布函數。

(7)對人工區域資料序列Qik計算Vsim。

(8)重復1 000次步驟(5)～(7)，計算得到的1 000個Vsim的均值μV和標準差σV，并根據式(3)計算得到修正后的統計量H*。

5.2 資料相關性對修正前后Hosking-Wallis檢驗的影響

圖4 改進前后HW檢驗結果對比

6 結 語

(1)利用蒙特卡洛模擬實驗，得出HW一致性檢驗中生成的人工資料序列間的平均相關系數隨資料長度的分布，并將其與實測資料進行對比。結果表明人工資料間的相關性均大于實測資料25%以上。

(2)從HW檢驗的定義出發，定性分析資料相關性影響H統計量值的原因，并通過蒙特卡羅模擬實驗，定量分析資料相關性以及不同的相關性分布對H值大小的影響。定量分析發現當資料間存在較大相關性時，HW檢驗對異質性的檢驗能力大幅下降，而當區域平均相關系數相同時，相關性導致H值的減小量相同。

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

(4)修正后的一致性檢驗通過優化水文一致區劃分，提升了針對各水文時間序列資料，特別是降雨、徑流序列頻率分析的效率及其結果的準確性。將改進后的降雨、徑流頻率分析結果與氣象降水量預報或流域水文模型徑流量預報相比較，能夠更為準確地劃定預報降水或洪水的量級、重現期，為洪水預警預報提供更加堅實的理論基礎。

□