基于Copula 函數(shù)的甌江設計潮位過程研究

林樂曼，周倩倩

（1.溫州市水利規(guī)劃發(fā)展研究中心，浙江 溫州 325000；2.溫州市水利局，浙江 溫州 325000）

1 問題的提出

設計潮位過程線是指相應于一定防御標準的潮位過程線，是河口海岸地區(qū)防洪排澇水動力模型計算的重要邊界條件。從防洪潮角度，高潮位毋庸置疑是決定設計潮位過程線至關重要的元素，但是潮差作為其重要組成特征，同樣也是排澇設計計算中的重要要素?？茖W地推求設計潮位過程線，對于確定工程規(guī)模以及識別分析工程防護風險程度而言，具有重要意義。

關于設計潮位過程線的計算方法，現(xiàn)行規(guī)范推薦的是同頻率設計法。該方法規(guī)避了同倍比法因基面不同而導致計算結(jié)果不同的問題，但是由于該方法是基于高潮位和潮差同頻率的假設，易出現(xiàn)設計潮差與實際潮差偏差較大乃至不合理的情況，缺乏對高潮位和潮差之間相依性的考慮，存在主觀性和隨意性。

近年來，得益于Copula 函數(shù)的發(fā)展，一些極端水文事件變量的相關性研究也有所進展。在設計潮位過程線的研究中，劉學等[1]利用Copula 函數(shù)建立年最高潮位和年最大潮差的二維聯(lián)合分布，并計算重現(xiàn)期，提出一種以高潮位為控制，結(jié)合同現(xiàn)重現(xiàn)期推求設計潮位過程線的方法。周月英等[2]將聯(lián)合分布應用于珠江口設計潮位過程線的計算，但是論文中采用的數(shù)據(jù)為年最高潮位與年最大潮差序列，事實上樣本選取時就破壞了潮位與潮差的相關性，結(jié)果并不能真實反映潮位的相關特性。薛曉鵬[3]等建立年最大潮位與同期潮差，保留了數(shù)據(jù)的相關性，采用最可能組合法計算設計高潮位條件下的最大潮差，并推求出設計潮位過程。為驗證該方法在甌江的應用效果，本文擬采用基于Copula函數(shù)建立甌江年最高潮位與同期潮差的聯(lián)合分布，并以高潮位為控制要素，采用條件最可能組合法，計算設計高潮位條件下概率最大的潮差，最后典型放大得到設計潮位過程。

2 基于Copula 函數(shù)聯(lián)合分布的設計潮位過程分析方法

根據(jù)Sklar 定理，存在一個Copula 函數(shù)C，能夠使多個邊緣分布聚合在一起構(gòu)建出多維聯(lián)合分布，若邊緣分布都是連續(xù)的，則C唯一。其數(shù)學表達式為[4]：

式（1）中：F（x1，x2，…，xd）表示聯(lián)合分布函數(shù)；ui=Fxi（x）是隨機變量的邊緣分布函數(shù)，ui∈[0，1] ；x1為設計高潮位，m；x2為同期潮差，m；d取2。

通過對C（u1，u2）函數(shù)求導，可得相應的Copula 聯(lián)合概率密度函數(shù)c（u1，u2），公式如下：

聯(lián)合分布的參數(shù)估計可采用邊際函數(shù)推斷法[4]，即求解參數(shù)通過2 個步驟完成，首先求解邊緣分布參數(shù)，然后根據(jù)極大似然原理由Kendall 秩相關系數(shù)推求出Copula 函數(shù)的參數(shù)。

2.1 邊緣分布建立

邊緣分布Fxi（x）即極端水文事件分析計算中經(jīng)常使用的單變量分布，本文采用標準規(guī)范推薦的P-Ⅲ型分布作為高潮位、同期潮差的邊緣分布[5]，概率密度函數(shù)公式為：

式（3）中：α=；β=EXCVCS/2；α0=EX（1-2CV/CS）；EX、CV、CS分別為樣本的均值、離差系數(shù)、偏態(tài)系數(shù)，一般可通過矩法結(jié)合目估適線法進行估計。為檢驗邊緣函數(shù)的擬合效果，采用卡方檢驗等方法[4]。

2.2 二維聯(lián)合分布建立

Copula 函數(shù)的類型和形式有很多種，Achimedean Copula 作為Copula 函數(shù)中重要的一簇函數(shù)，其結(jié)構(gòu)簡單，形式豐富，得到廣泛應用。本文采用Gumbel-Hougaard、Clayton、Frank、Ali-Mikhail-Haq Achimedean Copula 函數(shù)分別建立聯(lián)合分布，并通過多方案比選的方式來確定最優(yōu)Copula 聯(lián)合分布[4]，表達式見表1。表1 中，τ為Kendall 秩相關系數(shù)，公式為：

表1 3 種Achimedean Copula 函數(shù)表達式表

式（4）中：xi，xj，yi，yj為實測點據(jù)；n為樣本系列容量。

采用BIC 法（Bayesian information criterial）比選最優(yōu)Copula，BIC 統(tǒng)計量是描述理論聯(lián)合分布值與經(jīng)驗聯(lián)合分布值之間差距的統(tǒng)計值，也就是說統(tǒng)計量越小，擬合效果越好，公式如下：

例如，在講解《賣火柴的小女孩》一課的過程中，當學生對課文內(nèi)容有一定了解之后，教師就可以將角色扮演這一形式引入進來，讓學生以小組為單位，對課文進行改編，并將改編好的課本劇以角色扮演的方式進行情景再現(xiàn)，在這其中不僅可以彰顯學生的主體地位，也能使學生站在不同的視角下對文章內(nèi)容進行重新審視，使其的理解程度大大加深，使教學效果得以進一步的深化。

式（5）中：m為Copula 參數(shù)估計數(shù)；MSE=；Pc、P0分別為Copula 多元聯(lián)合分布計算值、聯(lián)合分布經(jīng)驗值。

最后采用CPI Rosenblatt 轉(zhuǎn)換法對建立的最優(yōu)模型進行擬合優(yōu)度檢驗，統(tǒng)計量選擇A—D，若統(tǒng)計觀測值小于置信水平為α的臨界值，則認為擬合效果較好，否則分布不合理[4]。

2.3 設計組合風險分析

在一般的防洪擋潮工程規(guī)劃建設中，高潮位是影響安全的主要因素。實際工程中，發(fā)生一定設計標準的高潮位時，潮差往往存在多種可能性，而它們出現(xiàn)的概率也不盡相同，相應概率的風險情況也不同。研究和掌握不同設計標準高潮位下，高潮位和潮差的組合概率風險值對于工程防汛風險分析評估具有重要意義。在滿足設計高潮位x1條件下，將潮差超過x2發(fā)生的概率作為其組合風險率[3]，數(shù)學表達式為：

2.4 條件最可能組合

理論上潮位和潮差的組合情況有無數(shù)種，且發(fā)生概率不同，工程中往往比較關心的是大概率事件。對于設計潮位過程而言，特別在防洪排澇設計中，我們往往以設計高潮位做為控制要素，基于此，本文重點研究特定重現(xiàn)期的設計高潮位和條件概率最有可能發(fā)生潮差的組合。條件概率最可能組合定義如下：當出現(xiàn)設計高潮位時，最可能出現(xiàn)的潮差，用數(shù)學公式可表達為：

式（7）顯然是一個非線性最優(yōu)化問題。由于本文中u2為x2的P-Ⅲ型概率分布函數(shù)，其值一般可用數(shù)值積分近似求解，無顯式函數(shù)表達式，故無法通過函數(shù)微分方式直接求解fx2 |x1的極值點，而遺傳算法在求解非線性最優(yōu)化問題中表現(xiàn)出色[6]，可通過生成一定數(shù)量的個體，以數(shù)值積分的方法求得分布函數(shù)值u2，繼而計算出由fx2 |x1構(gòu)建的適應度函數(shù)，再經(jīng)過一代代的交叉、選擇、變異、淘汰、遺傳等操作，最終找到最優(yōu)解。因此，本文嘗試將遺傳算法應用于條件最可能模型的求解。

遺傳算法求解的關鍵問題設置如下：①種群和種群規(guī)模。種群即所有個體的集合，x1k，x2k也就是個體，又稱為決策向量，個體數(shù)量又稱為種群規(guī)模（POP）。②個體編碼方法。采用實數(shù)編碼，考慮潮位實際情況，將上下限設定為各變量在TP=1.000 001 a、100 000 a 時相應的邊緣分布反函數(shù)值。xi=ximin+rand*（ximax-ximin），其中，ximax，ximin分別為變量xi的取值上下限，rand為[0，1]的隨機數(shù)。③適應度函數(shù)。極大值優(yōu)化問題，選用負構(gòu)造法構(gòu)造適應度函數(shù)g（x），即g（x）=B-fx2 |x1，其中，B為目標函數(shù)fx2 |x1界限的保守估計值。④停止準則。采用最大迭代次數(shù)停止準則。

3 甌江設計潮位過程計算

3.1 基于Copula 函數(shù)的聯(lián)合分布

3.1.1 邊緣分布

選取甌江溫州站1970—2018 年潮位數(shù)據(jù)，根據(jù)年最大法（AM 法）篩選年最大高潮位和同期相應潮差，共獲得49 組樣本數(shù)據(jù)，假定年最高潮位和同期潮差服從P-Ⅲ型分布，采用矩法結(jié)合目估適線法估計參數(shù)，并用K-S 檢驗法檢驗假設的合理性，計算結(jié)果見表2。

表2 邊緣分布函數(shù)參數(shù)計算結(jié)果表

由表2 可見，樣本統(tǒng)計檢驗值小于K-S 檢驗臨界值，即上述建立的邊緣分布通過了顯著性水平α=0.05 的K-S 檢驗，即假設的邊緣分布合理。

3.1.2 聯(lián)合分布

利用Gumbel-Hougaard、Frank、Clayton、Ali-Mikhail-Haq Copula 分別建立二元聯(lián)合模型，按照邊際函數(shù)推斷法計算得到聯(lián)合分布模型參數(shù)，并相應計算BIC 統(tǒng)計量，結(jié)果見表3。按照指標最小值對應的即為最優(yōu)模型準則，F(xiàn)rank Copula 即為最佳聯(lián)合模型。

表3 二元Copula 聯(lián)合分布模型參數(shù)計算結(jié)果表

通過CPI Rosenblatt 轉(zhuǎn)換檢驗法對Fank Copula擬合效果進行檢驗，統(tǒng)計量取A—D，樣本的檢驗統(tǒng)計值An2 為0.446 4，小于顯著性水平α=0.05 時的臨界值=0.604 2，滿足An檢驗要求，即二元Frank Copula 聯(lián)合分布是合理的。擬合的聯(lián)合分布見圖1。

圖1 建立的Frank-Copula 聯(lián)合分布圖

3.2 設計組合風險分析

由式（6）分別計算出重現(xiàn)期為2～100 a 的高潮位和同期不同重現(xiàn)期的潮差組合風險率，見表4。

表4 設計高潮位與同期潮差組合風險率表

由表4 可知，溫州站高重現(xiàn)期高潮位，遭遇高重現(xiàn)期潮差的組合風險值較低，遭遇可能性較小，隨著遭遇潮差的重現(xiàn)期減小，遭遇組合風險值逐漸變大，遭遇可能性變大；同一重現(xiàn)期的潮差，遭遇高潮位的重現(xiàn)期變小，組合風險率變小，但變化不大?？傮w上，溫州站設計高潮位與潮差不具備同頻性，不同重現(xiàn)期的高潮位遭遇低重現(xiàn)期潮差的概率更高，也就是說高潮位與低潮差的組合出現(xiàn)概率更高。

3.3 條件最可能組合及設計潮位過程線

實際防洪（潮）排澇設計中，對安全影響更大的因素是高潮位，所以潮位過程線的設計考慮以高潮位設計為主，同時考慮在給定高潮位條件下，出現(xiàn)幾率最大的潮差。根據(jù)式（5），利用遺傳算法計算給定重現(xiàn)期的設計高潮位下的條件最可能潮差，計算結(jié)果見表5。

表5 設計高潮位及條件最可能潮差計算結(jié)果表

由表5 可知，溫州站重現(xiàn)期為2～200 a 的設計高潮位下，條件最可能發(fā)生的潮差為1.90～2.34 a一遇，也就是說不同重現(xiàn)期的設計高潮位最可能遭遇的是低重現(xiàn)期的潮差。分析原因，主要是發(fā)生高重現(xiàn)期的潮水過程中，高潮位變大，事實上此時整個潮位抬升，風浪爬高，所以低潮位也會相應抬升，因此潮差不會加劇反而相對維持一個較小值。

由求出的設計高潮位與條件最可能潮差組合，按照雙重控制，縮放典型潮型即可得到設計潮位過程線。為比較不同方法的差異，將本文所用方法與僅設計高潮位控制、同頻率法等方法所推求的設計潮位過程線繪于同一圖（見圖2）。

圖2 不同方法推求的設計潮位過程線圖

由圖2 可見，幾種方法的設計高潮位相同，但是同頻率設計方法推求的設計潮位過程線的低潮位更低，不一定能反映真實的風暴增水過程。而采用本文方法推求的設計潮差是在考慮年最高潮位和同期潮位相關性的基礎上，根據(jù)實測資料擬合確定，數(shù)據(jù)統(tǒng)計基礎更強，在防洪排澇規(guī)劃與設計中，更不利于排澇，用于計算安全裕度更大，對于以考慮高潮位為主的防潮排澇工程設計而言更具合理性，有一定參考意義。

4 結(jié) 論

（1）采用Frank Copula 函數(shù)構(gòu)建甌江溫州站年最高潮位和同期潮差二元聯(lián)合分布合理可行。

（2）從概率統(tǒng)計上分析，甌江高潮位與潮差不具備同頻率的特性；高重現(xiàn)期的設計高潮位與低重現(xiàn)期的潮差組合概率風險更大，遭遇可能性更高。

（3）基于聯(lián)合分布，構(gòu)造給定重現(xiàn)期的設計高潮位下最可能的潮差模型，首次嘗試用遺傳算法求解該模型，結(jié)果證明遺傳算法有效，利用設計高潮位與條件最可能潮差組合推求潮位過程線的方法可行。該方法統(tǒng)計理論依據(jù)充分，較同頻率法更合理，且該方法設計低潮位較同頻率法更高，潮位過程偏不利，用于防洪排澇工程設計計算安全裕度更大。