林樂(lè)曼，周倩倩

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

1 問(wèn)題的提出

設(shè)計(jì)潮位過(guò)程線是指相應(yīng)于一定防御標(biāo)準(zhǔn)的潮位過(guò)程線，是河口海岸地區(qū)防洪排澇水動(dòng)力模型計(jì)算的重要邊界條件。從防洪潮角度，高潮位毋庸置疑是決定設(shè)計(jì)潮位過(guò)程線至關(guān)重要的元素，但是潮差作為其重要組成特征，同樣也是排澇設(shè)計(jì)計(jì)算中的重要要素。科學(xué)地推求設(shè)計(jì)潮位過(guò)程線，對(duì)于確定工程規(guī)模以及識(shí)別分析工程防護(hù)風(fēng)險(xiǎn)程度而言，具有重要意義。

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

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

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

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

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

通過(guò)對(duì)C（u1，u2）函數(shù)求導(dǎo)，可得相應(yīng)的Copula 聯(lián)合概率密度函數(shù)c（u1，u2），公式如下：

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

2.1 邊緣分布建立

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

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

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

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

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

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

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

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

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

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

2.3 設(shè)計(jì)組合風(fēng)險(xiǎn)分析

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

2.4 條件最可能組合

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

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

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

3 甌江設(shè)計(jì)潮位過(guò)程計(jì)算

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

3.1.1 邊緣分布

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

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

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

3.1.2 聯(lián)合分布

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

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

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

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

3.2 設(shè)計(jì)組合風(fēng)險(xiǎn)分析

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

表4 設(shè)計(jì)高潮位與同期潮差組合風(fēng)險(xiǎn)率表

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

3.3 條件最可能組合及設(shè)計(jì)潮位過(guò)程線

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

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

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

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

圖2 不同方法推求的設(shè)計(jì)潮位過(guò)程線圖

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

4 結(jié) 論

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

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

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