楊 蕊，吳時(shí)強(qiáng)，高學(xué)平，張 晨

（1.天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350；2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029）

1 研究背景

21世紀(jì)以來(lái)，隨著我國(guó)工業(yè)化、城鎮(zhèn)化的快速發(fā)展，水資源短缺、水生態(tài)環(huán)境惡化問(wèn)題日益突出，為適應(yīng)新時(shí)期水利發(fā)展的需要，我國(guó)陸續(xù)興建了一批以水資源調(diào)配和水生態(tài)環(huán)境保護(hù)與修復(fù)為目的的河湖連通工程，如南水北調(diào)、引江濟(jì)太、引黃濟(jì)淀等河湖連通工程［1］。由于河湖水系連通是一個(gè)充斥著模糊性、隨機(jī)性以及未確知性信息的復(fù)雜系統(tǒng)工程［2］，這些工程在帶來(lái)社會(huì)經(jīng)濟(jì)效益和水生態(tài)環(huán)境效益的同時(shí)，也不可避免地帶來(lái)了一定的負(fù)面效應(yīng)［3-4］，如調(diào)水與洪水遭遇引起的防洪風(fēng)險(xiǎn)、調(diào)水引起的污染物轉(zhuǎn)移或滯留風(fēng)險(xiǎn)等。因此，如何識(shí)別與管控河湖水系連通潛在風(fēng)險(xiǎn)，實(shí)現(xiàn)科學(xué)、高效、安全的連通運(yùn)行管理，以減少連通過(guò)程各種潛在損失，是亟需解決的技術(shù)難題。

近年來(lái)，日益突出的環(huán)境污染問(wèn)題對(duì)生態(tài)系統(tǒng)及國(guó)民經(jīng)濟(jì)都造成了一定損失，使水環(huán)境風(fēng)險(xiǎn)研究得到廣泛的關(guān)注。借鑒調(diào)水工程水文風(fēng)險(xiǎn)的概念［5］，本文將河湖連通水環(huán)境風(fēng)險(xiǎn)定義為由連通工程引起的區(qū)域污染物遷移和滯留等風(fēng)險(xiǎn)事件造成的損失或不利后果的可能性。對(duì)于水環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)，學(xué)者們先后提出了健康評(píng)價(jià)模型［6］、隨機(jī)模型［7］、灰色模型［8］、模糊模型［9］、信息熵模型［10］等以及和其它理論耦合在一起的水環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)模型［11］。這些模型多用于對(duì)水體污染風(fēng)險(xiǎn)水平的衡量，在河湖連通工程環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)中也有部分應(yīng)用［12-15］。以上模型多應(yīng)用于單一水質(zhì)指標(biāo)的分析，并將多項(xiàng)關(guān)鍵指標(biāo)組合構(gòu)成指標(biāo)集對(duì)某一區(qū)域水環(huán)境狀態(tài)進(jìn)行整體評(píng)價(jià)，并未直接考慮各指標(biāo)間的相關(guān)性。而Copula函數(shù)作為構(gòu)造多變量的相關(guān)結(jié)構(gòu)和聯(lián)合分布的連接函數(shù)［16-18］，相較于其它模型，可以直接建立起各項(xiàng)關(guān)鍵指標(biāo)的聯(lián)合分布，進(jìn)而對(duì)由各項(xiàng)關(guān)鍵指標(biāo)超過(guò)標(biāo)準(zhǔn)限值形成的多因子聯(lián)合風(fēng)險(xiǎn)進(jìn)行分析［19］。但是，在使用Copula函數(shù)對(duì)多元相關(guān)結(jié)構(gòu)進(jìn)行建模時(shí)存在“維數(shù)災(zāi)難”的問(wèn)題，且有單一Copula函數(shù)擬合效果不理想的限制，如橢圓Copula函數(shù)只能刻畫(huà)對(duì)稱的相關(guān)結(jié)構(gòu)，三維及以上維度的阿基米德Copula函數(shù)只適合刻畫(huà)變量的正相關(guān)結(jié)構(gòu)［20-21］。當(dāng)前Copula函數(shù)在水環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)中多應(yīng)用于水體富營(yíng)養(yǎng)化評(píng)價(jià)，分析不同環(huán)境因子之間的相關(guān)性、單個(gè)富營(yíng)養(yǎng)化因子對(duì)風(fēng)險(xiǎn)的影響程度以及構(gòu)建富營(yíng)養(yǎng)化綜合評(píng)價(jià)指標(biāo)等［22-24］，多應(yīng)用于二維、三維。因此，如何突破維數(shù)限制，更好地評(píng)估高維變量聯(lián)合風(fēng)險(xiǎn)，反映風(fēng)險(xiǎn)事件多因子相互影響并共同作用的真實(shí)特征，也是當(dāng)前面臨的科學(xué)問(wèn)題。

本文旨在以下兩個(gè)方面對(duì)現(xiàn)有研究做出貢獻(xiàn)：首先，在方法上基于一種新的Copula 函數(shù)類型，即Vine Copula 函數(shù)對(duì)水環(huán)境風(fēng)險(xiǎn)多因子相關(guān)結(jié)構(gòu)進(jìn)行建模，并在建模過(guò)程中，通過(guò)比選15 種二元Copula（包含大部分可能的相關(guān)結(jié)構(gòu)［25］）優(yōu)化各變量（條件或無(wú)條件）間的相關(guān)結(jié)構(gòu)，以使模型能夠最恰當(dāng)?shù)乜坍?huà)多變量間相關(guān)結(jié)構(gòu)的真實(shí)特征；其次，在應(yīng)用上以當(dāng)下社會(huì)關(guān)注度較高的南水北調(diào)工程伴生水環(huán)境問(wèn)題為出發(fā)點(diǎn)，開(kāi)展河湖水系連通伴生水環(huán)境聯(lián)合風(fēng)險(xiǎn)的識(shí)別研究，為進(jìn)一步實(shí)現(xiàn)水系連通伴生水環(huán)境風(fēng)險(xiǎn)的管控提供科學(xué)依據(jù)。

2 研究方法

2.1 邊緣分布擬合構(gòu)建聯(lián)合分布之前，須先確定各風(fēng)險(xiǎn)因子的邊緣分布函數(shù)。隨機(jī)變量邊緣分布的經(jīng)驗(yàn)頻率采用水文中常用的Gringorten公式計(jì)算［26］：

選用6種常見(jiàn)的水文變量分布，即正態(tài)分布、指數(shù)分布、logistic分布、對(duì)數(shù)正態(tài)分布、伽馬分布和威布爾分布，分別對(duì)各變量觀測(cè)序列進(jìn)行擬合，采用極大似然法評(píng)估參數(shù)，并通過(guò)AIC 擬合優(yōu)度檢驗(yàn)比選出最適合的邊緣分布。AIC值計(jì)算公式［27］如下：

式中：lmax為極大似然函數(shù)值；k為參數(shù)個(gè)數(shù)。AIC值最小則擬合最優(yōu)。

2.2 Copula函數(shù)Copula函數(shù)是定義在[0，1]n區(qū)間上的多維聯(lián)合分布函數(shù)，它可以將多個(gè)隨機(jī)變量的邊際分布連接起來(lái)構(gòu)造聯(lián)合分布。

Sklar定理［28］：令H為n維聯(lián)合分布函數(shù)，F(xiàn)1，F(xiàn)2，…，F(xiàn)n為其邊緣分布，那么存在唯一的Copu?la函數(shù)C，使得對(duì)?x∈R n，有

若F1，F(xiàn)2，…，F(xiàn)n連續(xù)，則C是唯一的。相反地，如果C是一個(gè)Copula函數(shù)，F(xiàn)1，F(xiàn)2，…，F(xiàn)n是分布函數(shù)，則上式定義的H是邊緣分布為F1，F(xiàn)2，…，F(xiàn)n的一個(gè)聯(lián)合分布函數(shù)。

Copula 函數(shù)有很多種類型［29-31］，主要分為橢圓型和阿基米德型。表1列出了本文使用的各橢圓Copula和阿基米德Copula函數(shù)的二元分布函數(shù)表達(dá)式及其相關(guān)性質(zhì)。

2.3 Vine Copula 函數(shù)當(dāng)分析兩變量間的相關(guān)結(jié)構(gòu)時(shí)，有多種Copula 函數(shù)可供選擇，例如常用的Gaussian Copula［32］、Student-t-Copula［32］、Clayton Copula［29］、Frank Copula［32］、Gumbel Copula［32-34］和Joe Copula［35］等。但一旦需要對(duì)多元相關(guān)結(jié)構(gòu)進(jìn)行刻畫(huà)時(shí)，可供選擇的Copula函數(shù)就變得十分有限，尤其對(duì)于阿基米德類Copula函數(shù)，隨著維度的增加其估計(jì)和推斷難度會(huì)迅速增加［21］，并且當(dāng)維度達(dá)到三維及以上時(shí)只能刻畫(huà)正相關(guān)結(jié)構(gòu)［20］?，F(xiàn)有研究中常用的多元Copula函數(shù)有Gaussian，Student t和Clay?ton Copula，但以上多元Copula 都具有一定局限性，例如Gaussian 和Student t Copula 只能刻畫(huà)對(duì)稱的相關(guān)結(jié)構(gòu)，Clayton Copula 雖能刻畫(huà)非對(duì)稱相關(guān)結(jié)構(gòu)，但卻只存在下尾相關(guān)性［21］。因此，本文引入一種更為靈活和直觀的多元相關(guān)結(jié)構(gòu)建模工具Vine Copula。

表1 二元Copula函數(shù)表達(dá)式及其相關(guān)性質(zhì)

Vine Copula 是由Bedford 和Cooke（2001）引入的一種具有圖形結(jié)構(gòu)的多變量Copula 建模函數(shù)［36］。它基于條件Copula函數(shù)和藤式（Vine）圖形建模工具，利用Pair Copula將多元聯(lián)合分布進(jìn)行分解，從而建立起Vine Copula 模型。一般來(lái)說(shuō)，n維藤結(jié)構(gòu)可以由n（n-1）個(gè)二元Copula 函數(shù)以（n-1）層次樹(shù)構(gòu)建。每個(gè)樹(shù)上有多個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)之間的連線為邊，而每一個(gè)節(jié)點(diǎn)代表的就是變量或條件變量，每一條邊代表著相應(yīng)的Pair Copula，最后每棵樹(shù)組合在一起構(gòu)成一個(gè)完整的藤結(jié)構(gòu)。常用的藤結(jié)構(gòu)有Cvine 和D-vine，其中C-vine 多用于擬合多變量中具有一個(gè)關(guān)鍵變量的情形，本文選用C-vine 進(jìn)行模型構(gòu)建。C-vine Copula的密度函數(shù)表達(dá)式如下：

式中，下標(biāo)j表示樹(shù)層級(jí)的遍歷，下標(biāo)i表示每棵樹(shù)中節(jié)點(diǎn)的遍歷。四維變量C-vine Copula結(jié)構(gòu)如圖1所示。

圖1 C-vine Copula四維變量結(jié)構(gòu)圖

該方法的優(yōu)勢(shì)在于，相較于傳統(tǒng)多元Copula函數(shù)，其在構(gòu)建聯(lián)合分布函數(shù)時(shí)不僅可以將各變量的邊緣分布函數(shù)與相關(guān)結(jié)構(gòu)分開(kāi)建模，并且允許多個(gè)變量?jī)蓛芍g的相關(guān)結(jié)構(gòu)有所不同［21］。由于不同的Copula 函數(shù)代表著特定的相關(guān)結(jié)構(gòu)，在選取Pair Copula 函數(shù)時(shí)，本文從包含大部分相關(guān)結(jié)構(gòu)的15種Copula函數(shù)中，選擇最為恰當(dāng)?shù)腃opula函數(shù)來(lái)刻畫(huà)變量（條件或無(wú)條件）間的真實(shí)相關(guān)結(jié)構(gòu)，所構(gòu)建的聯(lián)合分布函數(shù)比傳統(tǒng)的多元Copula函數(shù)更具靈活性［21，25］。

2.4 C-vine Copula模型構(gòu)建及驗(yàn)證計(jì)算多個(gè)變量?jī)蓛芍g的Kendall系數(shù)，將與其它變量Kendall系數(shù)絕對(duì)值之和最大的變量作為中樞變量［21］。根據(jù)AIC準(zhǔn)則，從表1所列出的二元橢圓Copula、阿基米德Copula以及Clayton、Gumbel、Joe Copula 旋轉(zhuǎn)90°、180°、270°共15個(gè)Copula函數(shù)中，分別比選出最佳的Copula函數(shù)作為第一層樹(shù)（T1）的中樞變量與其它變量之間的Pair Copula，并通過(guò)貝葉斯估計(jì)方法估計(jì)出T1中各Pair Copula密度函數(shù)的參數(shù)值。利用T1所估參數(shù)生成T2各變量的觀察值，并重復(fù)建立T1的過(guò)程，可得到T2最優(yōu)Pair Copula類型及密度函數(shù)的參數(shù)值。重復(fù)以上過(guò)程直至最后一層樹(shù)建立完成［25］。最后利用極大似然估計(jì)對(duì)藤結(jié)構(gòu)中每個(gè)Pair Copula函數(shù)的參數(shù)估計(jì)進(jìn)行修正［37］。

指定的分布模型能否很好的擬合變量的實(shí)際分布，對(duì)Copula函數(shù)能否正確的描述變量間的相關(guān)性結(jié)構(gòu)至關(guān)重要，因此要建立分布的檢驗(yàn)和擬合度評(píng)價(jià)。常見(jiàn)的評(píng)價(jià)方法主要分為圖解法和解析法，圖形法相對(duì)于解析法比較直觀，解析法可以量化擬合效果。本文以箱型圖的方式，對(duì)變量間的原始Kendall系數(shù)和模型所得變量間的Kendall系數(shù)進(jìn)行比對(duì)，以驗(yàn)證模型的可靠性［38］。

3 案例應(yīng)用

圖2 南四湖水系及連通工程示意圖

3.1 研究區(qū)域概況及數(shù)據(jù)來(lái)源南四湖屬沂沭泗的泗河水系，是南陽(yáng)湖、獨(dú)山湖、昭陽(yáng)湖和微山湖4個(gè)相連湖泊的總稱，位于山東省西南部，自北向南流動(dòng)，全湖面積1266 km2，湖泊容積16.06億m3。作為南水北調(diào)工程?hào)|線上重要的輸水通道和調(diào)蓄湖泊，南水北調(diào)東線工程的運(yùn)行意味著其連通性的改變，必定引起其水環(huán)境因子的變化，區(qū)域概況如圖2所示?；谥笜?biāo)的針對(duì)性、實(shí)用性和可行性，以及對(duì)南四湖水環(huán)境風(fēng)險(xiǎn)的綜合反映能力，本文從水質(zhì)指標(biāo)和富營(yíng)養(yǎng)化指標(biāo)中選取代表性指標(biāo)總磷（TP）、總氮（TN）、氨氮（NH3-N）和葉綠素a（Chl-a）2008—2014年月均觀測(cè)數(shù)據(jù)（數(shù)據(jù)來(lái)源于南陽(yáng)監(jiān)測(cè)點(diǎn)［39］，各指標(biāo)的統(tǒng)計(jì)特征列于表2），采用C-vine Copula函數(shù)構(gòu)建四維聯(lián)合分布，分別構(gòu)建了南水北調(diào)東線一期工程運(yùn)行前（2008.01—2013.10）和運(yùn)行后（2013.11—2014.12）南四湖水環(huán)境多因子聯(lián)合風(fēng)險(xiǎn)模型。結(jié)合南水北調(diào)東線工程Ⅲ類水水質(zhì)要求，定義TP、TN、NH3-N濃度值劣于Ⅲ類水（Chl-a對(duì)應(yīng)中營(yíng)養(yǎng)）標(biāo)準(zhǔn)限值時(shí)為存在水環(huán)境聯(lián)合風(fēng)險(xiǎn)。通過(guò)分析南四湖在工程運(yùn)行前后水環(huán)境聯(lián)合風(fēng)險(xiǎn)對(duì)TP、TN、NH3-N 和Chl-a 各因子敏感性的變化，提取出工程運(yùn)行后聯(lián)合風(fēng)險(xiǎn)的關(guān)鍵風(fēng)險(xiǎn)因子。之后，設(shè)置運(yùn)行后關(guān)鍵風(fēng)險(xiǎn)因子不同風(fēng)險(xiǎn)狀態(tài)組合情景，計(jì)算并分析風(fēng)險(xiǎn)發(fā)生的不利情況，識(shí)別出運(yùn)行后潛在水環(huán)境聯(lián)合風(fēng)險(xiǎn)。

表2 運(yùn)行前/后南陽(yáng)監(jiān)測(cè)站水質(zhì)指標(biāo)統(tǒng)計(jì)特征表

3.2 C-vine Copula模型構(gòu)建

3.2.1 建立邊緣分布 表3所示為南水北調(diào)東線工程運(yùn)行前、后南四湖TP、TN、NH3-N 和Chl-a 環(huán)境因子分別擬合6種理論分布所對(duì)應(yīng)的AIC值，比選可得運(yùn)行前TP、TN、NH3-N和Chl-a的最優(yōu)擬合分布為對(duì)數(shù)正態(tài)分布；運(yùn)行后TP、TN的最優(yōu)擬合分布為對(duì)數(shù)正態(tài)分布，NH3-N和Chl-a的最優(yōu)擬合分布為伽馬分布和威布爾分布。

表3 運(yùn)行前/后各變量擬合不同分布類型對(duì)應(yīng)的AIC值

3.2.2 C-vine Copula模型的構(gòu)建和模型檢驗(yàn)

（1）構(gòu)建C-vine Copula模型。運(yùn)行前后TN與其他3種變量的Kendall'sτ絕對(duì)值之和Si最大，即表明TN 與其他3 個(gè)變量的相關(guān)性最強(qiáng)，將其作為中樞變量。按照上文所述模型構(gòu)建方法，構(gòu)建運(yùn)行前、后兩個(gè)C-vine Copula 四維水環(huán)境風(fēng)險(xiǎn)因子聯(lián)合分布模型，模型的Pair Copula 函數(shù)選擇結(jié)果及參數(shù)列于表4，模型圖如圖3所示。

表4 運(yùn)行前/后C-vine Pair Copula的函數(shù)選擇結(jié)果及參數(shù)估計(jì)

圖3 運(yùn)行前/后C-vine Copula四維水環(huán)境風(fēng)險(xiǎn)因子聯(lián)合分布模型

（2）模型驗(yàn)證。將構(gòu)建的C-vine Copula四維水環(huán)境因子聯(lián)合分布模型中的Pair Copula類型和參數(shù)估計(jì)的結(jié)果帶入仿真函數(shù)產(chǎn)生100組仿真數(shù)據(jù)，計(jì)算TP、TN、NH3-N和Chl-a兩兩因子之間的Kend?all 系數(shù)，循環(huán)以上步驟100 次，將100 組仿真數(shù)據(jù)所得Kendall 系數(shù)整理為箱型圖，并將TP、TN、NH3-N和Chl-a各環(huán)境因子之間的原始Kendall系數(shù)標(biāo)于圖中。如圖4所示，運(yùn)行前后各環(huán)境因子之間的原始Kendall系數(shù)均落在箱型圖上、下四分位數(shù)之間，并且大部分在中位數(shù)附近，說(shuō)明所建立的兩個(gè)C-vine Copula模型對(duì)變量的相關(guān)結(jié)構(gòu)刻畫(huà)較為可靠。

圖4 運(yùn)行前/后Kendall系數(shù)模型驗(yàn)證

3.3 風(fēng)險(xiǎn)情景設(shè)計(jì)

3.3.1 各因子敏感性分析情景 為分析南水北調(diào)東線工程運(yùn)行前后南四湖多因子聯(lián)合水環(huán)境風(fēng)險(xiǎn)中TP、TN、NH3-N和Chl-a 各因子的敏感性，識(shí)別出關(guān)鍵風(fēng)險(xiǎn)因子進(jìn)而實(shí)現(xiàn)其風(fēng)險(xiǎn)組合分析，在完成對(duì)C-vine Copula模型的構(gòu)建后，結(jié)合我國(guó)地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)（GB 3838-2002）及OECD（Organization for Economic Co-operation and Development）提出的富營(yíng)養(yǎng)化湖泊（溫帶湖泊）營(yíng)養(yǎng)水平劃分標(biāo)準(zhǔn)，以工程運(yùn)行前后TP、TN、NH3-N和Chl-a原始數(shù)據(jù)序列的均值作為基準(zhǔn)值，分別按I—V類水對(duì)應(yīng)的標(biāo)準(zhǔn)限值（Chl-a采用貧營(yíng)養(yǎng)—超富營(yíng)養(yǎng)標(biāo)準(zhǔn)限值）改變其中一個(gè)變量的濃度值，分析工程運(yùn)行前后各環(huán)境因子對(duì)水環(huán)境聯(lián)合風(fēng)險(xiǎn)的影響程度，并識(shí)別出關(guān)鍵風(fēng)險(xiǎn)因子。工程運(yùn)行前后各因子敏感性分析計(jì)算情景如表5所示。

3.3.2 關(guān)鍵因子風(fēng)險(xiǎn)組合情景 由于南水北調(diào)東線工程南四湖供水水質(zhì)要求為Ⅲ類水，故關(guān)鍵因子風(fēng)險(xiǎn)組合情景設(shè)置主要考慮工程運(yùn)行后關(guān)鍵風(fēng)險(xiǎn)因子TN、NH3-N、Chl-a濃度值劣于Ⅲ類水（Chl-a對(duì)應(yīng)中營(yíng)養(yǎng)）標(biāo)準(zhǔn)限值時(shí)，不同風(fēng)險(xiǎn)狀態(tài)（即Ⅳ類和Ⅴ類，Chl-a對(duì)應(yīng)富營(yíng)養(yǎng)和超營(yíng)養(yǎng)）組合下的系列風(fēng)險(xiǎn)情景，從而識(shí)別出工程運(yùn)行后伴生水環(huán)境多因子聯(lián)合風(fēng)險(xiǎn)發(fā)生的最不利情況。具體設(shè)置如下：當(dāng)TP不是關(guān)鍵風(fēng)險(xiǎn)因子時(shí)，其濃度值分別取工程運(yùn)行前、運(yùn)行后數(shù)據(jù)序列的均值（0.06/0.04 mg·L-1），TN、NH3-N分別按照Ⅳ（1.5 mg·L-1）、Ⅴ（2 mg·L-1）類水限值，Chl-a按照富營(yíng)養(yǎng)（25 μg·L-1）、超營(yíng)養(yǎng)（35μg·L-1）狀態(tài)限值進(jìn)行組合，可得到工程運(yùn)行前后各6組計(jì)算工況，如表6所示。

表5 各因子敏感性分析情景（單位：TP/TN/NH3-N（mg·L-1），Chl-a（μg·L-1））

表6 關(guān)鍵因子風(fēng)險(xiǎn)組合情景（單位：TP/TN/NH3-N（mg·L-1），Chl-a（μg·L-1））

3.4 風(fēng)險(xiǎn)識(shí)別結(jié)果分析

3.4.1 各因子敏感性分析結(jié)果 由表2可知，工程運(yùn)行后TP、TN、Chl-a濃度值均有所降低，TP由Ⅳ類水變?yōu)棰箢愃?；NH3-N濃度略有上升，由Ⅱ類水變?yōu)棰箢愃??？傮w上，南水北調(diào)東線工程運(yùn)行后南四湖南陽(yáng)站處水質(zhì)有所改善，與已有研究結(jié)果一致［40］。

采用四維水環(huán)境因子聯(lián)合分布的概率密度值分析水環(huán)境聯(lián)合風(fēng)險(xiǎn)對(duì)各因子的敏感性，在各因子處于同一分類標(biāo)準(zhǔn)限值時(shí)，其所對(duì)應(yīng)的聯(lián)合分布概率密度值越大，說(shuō)明該因子對(duì)聯(lián)合分布的變化率影響越大，即聯(lián)合風(fēng)險(xiǎn)對(duì)其越敏感。由圖5（a）可見(jiàn)，當(dāng)TP、TN、NH3-N和Chl-a分別為Ⅳ、Ⅴ類水（Chl-a對(duì)應(yīng)富營(yíng)養(yǎng)和超營(yíng)養(yǎng)）限值時(shí)，各因子對(duì)應(yīng)的聯(lián)合分布概率密度值相對(duì)其他工況均較大，但各因子之間的差別較小，說(shuō)明各因子對(duì)于聯(lián)合風(fēng)險(xiǎn)均較為敏感，其中最敏感的因子是TP，NH3-N 次之。對(duì)比工程運(yùn)行前（圖5（a））和運(yùn)行后（圖5（b）），在Ⅳ、Ⅴ類水（Chl-a 對(duì)應(yīng)富營(yíng)養(yǎng)和超營(yíng)養(yǎng)）工況下，運(yùn)行后各因子聯(lián)合概率密度值均有所降低，說(shuō)明南水北調(diào)東線工程運(yùn)行后，水環(huán)境聯(lián)合風(fēng)險(xiǎn)對(duì)TP、TN、NH3-N和Chl-a的敏感性均有所降低。其中TP敏感性變換最大，從最敏感的因子變?yōu)槊舾行宰钊醯囊蜃?，敏感性減小幅度較大，可將其從關(guān)鍵風(fēng)險(xiǎn)因子中排除；而TN、NH3-N、Chl-a 概率密度值相近且較運(yùn)行前減小幅度不大，因此，認(rèn)為這3個(gè)環(huán)境因子是工程運(yùn)行后水環(huán)境聯(lián)合風(fēng)險(xiǎn)的關(guān)鍵因子。

為進(jìn)一步分析運(yùn)行后模型中各變量對(duì)概率密度的影響，比對(duì)運(yùn)行后模型在運(yùn)行前20組工況和運(yùn)行后20組工況下的概率密度結(jié)果。運(yùn)行前和運(yùn)行后20組工況的區(qū)別在于各變量的均值不同，運(yùn)行后工況中各變量均值除NH3-N升高外，其它3個(gè)變量均有所降低。如圖5（c）、（d）所示，TP在Ⅰ—Ⅴ級(jí)時(shí)，運(yùn)行后工況概率密度較運(yùn)行前各工況呈現(xiàn)增大后減小趨勢(shì)，并在Ⅲ-Ⅳ級(jí)時(shí)出現(xiàn)極值點(diǎn)，說(shuō)明Ⅰ-Ⅲ級(jí)時(shí)TN、NH3-N、Chl-a相較于TP發(fā)揮主導(dǎo)作用，后者在Ⅳ-Ⅴ級(jí)時(shí)發(fā)揮主導(dǎo)作用。TN在Ⅰ—Ⅴ級(jí)時(shí)，概率密度始終呈現(xiàn)上升趨勢(shì)，可認(rèn)為T(mén)N 始終對(duì)概率密度的影響較大。NH3-N 處于I 級(jí)時(shí)，當(dāng)其它3個(gè)變量濃度值降低，概率密度降低，但從Ⅱ級(jí)便出現(xiàn)上升趨勢(shì)，說(shuō)明NH3-N對(duì)概率密度的影響跨度最大。Chl-a 在Ⅰ—Ⅴ級(jí)時(shí)，概率密度呈現(xiàn)降低后升高走勢(shì)，并在Ⅲ-Ⅳ級(jí)之間出現(xiàn)極值點(diǎn)，說(shuō)明Ⅰ-Ⅲ級(jí)時(shí)TP、TN、NH3-N 相較于Chl-a 發(fā)揮主導(dǎo)作用，而隨著Chl-a 濃度的增加，其在Ⅳ-Ⅴ級(jí)時(shí)發(fā)揮主導(dǎo)作用。

圖5 運(yùn)行前/后模型（聯(lián)合分布）不同工況下的概率密度值

同時(shí)，應(yīng)用傳統(tǒng)多元Copula函數(shù)構(gòu)建四維水環(huán)境因子聯(lián)合分布，即橢圓Copula函數(shù)：（1）Gauss?ian Copula，（2）Student-t-Copula；阿基米德Copula 函數(shù)：（3）Clayton Copula，（4）Frank Copula，（5）Gumbel Copula。通過(guò)聯(lián)合分布概率密度值分析水環(huán)境聯(lián)合風(fēng)險(xiǎn)對(duì)各因子的敏感性，并與Vine Copula函數(shù)計(jì)算結(jié)果進(jìn)行比對(duì)，得到橢圓族Copula 函數(shù)計(jì)算結(jié)果（1）（2）相似，阿基米德族Copula 函數(shù)（3）（4）（5）計(jì)算結(jié)果相似。對(duì)比工程運(yùn)行前和運(yùn)行后的計(jì)算結(jié)果，在IV、V類水（Chl-a 對(duì)應(yīng)富營(yíng)養(yǎng)和超營(yíng)養(yǎng)）工況下，二者均能排除TP為關(guān)鍵風(fēng)險(xiǎn)因子，與Vine Copula 函數(shù)計(jì)算結(jié)果一致；但是，運(yùn)行后隨著NH3-N、Chl-a 濃度值增大，傳統(tǒng)多元Copula 函數(shù)的聯(lián)合分布概率密度值始終為零，無(wú)法反映NH3-N、Chl-a為關(guān)鍵風(fēng)險(xiǎn)因子。具體分析可能存在兩方面原因，一是由于橢圓Copula函數(shù)只能刻畫(huà)對(duì)稱的相關(guān)結(jié)構(gòu)，而不同變量間的相關(guān)關(guān)系并非都是對(duì)稱結(jié)構(gòu)；二是在計(jì)算過(guò)程中由于阿基米德Copula函數(shù)只能刻畫(huà)變量間的正相關(guān)結(jié)構(gòu)，在擬合負(fù)相關(guān)變量時(shí)會(huì)以零值覆蓋負(fù)值Kendall系數(shù)，進(jìn)而影響Copula函數(shù)參數(shù)的計(jì)算，致使模型擬合效果失真。受文章篇幅限制，未將對(duì)比結(jié)果圖列于正文中。

3.4.2 關(guān)鍵因子風(fēng)險(xiǎn)組合情景結(jié)果分析 基于各因子敏感性分析結(jié)果，計(jì)算了關(guān)鍵風(fēng)險(xiǎn)因子TN、NH3-N、Chl-a 不同風(fēng)險(xiǎn)狀態(tài)組合下6 種工況在南水北調(diào)東線工程運(yùn)行前后的聯(lián)合分布概率密度值，結(jié)果如圖6所示。由圖6坐標(biāo)定義可知，若計(jì)算工況所得值處于對(duì)角線上方則說(shuō)明運(yùn)行后聯(lián)合分布的概率密度值大于運(yùn)行前概率密度值，且距離對(duì)角線越遠(yuǎn)說(shuō)明二者差異越大，即該工況下運(yùn)行后風(fēng)險(xiǎn)變化速率越大，響應(yīng)時(shí)間越短。對(duì)比工況（1，2）（3，4）（5，6）（括號(hào)中第二工況相對(duì)于前一工況僅Chl-a濃度值由富營(yíng)養(yǎng)變?yōu)槌瑺I(yíng)養(yǎng)）計(jì)算結(jié)果可知，隨著TN、NH3-N、Chl-a濃度的增加，相應(yīng)工況的計(jì)算結(jié)果距離對(duì)角線越遠(yuǎn)，說(shuō)明各因子濃度越大時(shí)，運(yùn)行后該濃度值下的聯(lián)合風(fēng)險(xiǎn)變化速率越大，即不確定性越大。其中，工況1，即TN和NH3-N為Ⅳ類水標(biāo)準(zhǔn)限值、Chl-a為富營(yíng)養(yǎng)狀態(tài)限值時(shí)，運(yùn)行前后四維水環(huán)境聯(lián)合風(fēng)險(xiǎn)變化不大；而工況6是運(yùn)行后四維水環(huán)境聯(lián)合風(fēng)險(xiǎn)最不利的狀況，即TN和NH3-N達(dá)到Ⅴ類水標(biāo)準(zhǔn)限值、Chl-a為超營(yíng)養(yǎng)狀態(tài)（35 μg·L-1）時(shí)，運(yùn)行后水環(huán)境聯(lián)合風(fēng)險(xiǎn)變化速率增幅最大，不確定性提高，需要及時(shí)響應(yīng)。同時(shí)，由于運(yùn)行后模型中TN、NH3-N對(duì)概率密度的影響范圍較其它變量更大，建議重點(diǎn)關(guān)注TN和NH3-N指標(biāo)。

圖6 運(yùn)行前/后關(guān)鍵因子風(fēng)險(xiǎn)組合情景聯(lián)合分布概率密度值

4 結(jié)論

（1）本文構(gòu)建了一種基于Vine Copula函數(shù)的水環(huán)境多因子聯(lián)合風(fēng)險(xiǎn)模型，該模型避免了構(gòu)建多變量聯(lián)合分布函數(shù)時(shí)各邊緣分布不同的問(wèn)題、二元Copula 函數(shù)“維數(shù)災(zāi)難”的問(wèn)題以及單一多元Copula函數(shù)對(duì)兩兩變量間相關(guān)結(jié)構(gòu)刻畫(huà)失真等弊端，能夠較好地綜合多風(fēng)險(xiǎn)因子，為水系連通伴生水環(huán)境風(fēng)險(xiǎn)識(shí)別提供可行的技術(shù)方法。

（2）將所構(gòu)建的多因子聯(lián)合風(fēng)險(xiǎn)模型應(yīng)用于南水北調(diào)東線工程運(yùn)行后南四湖水環(huán)境聯(lián)合風(fēng)險(xiǎn)的識(shí)別，得出以下結(jié)論：南水北調(diào)東線工程運(yùn)行后，當(dāng)TN和NH3-N達(dá)到V類水標(biāo)準(zhǔn)限值、Chl-a為超營(yíng)養(yǎng)狀態(tài)（35 μg·L-1）時(shí)，南四湖水環(huán)境聯(lián)合風(fēng)險(xiǎn)最不利，發(fā)生聯(lián)合風(fēng)險(xiǎn)的可能性和速率較工程運(yùn)行前大幅增加，需及時(shí)采取措施。建議重點(diǎn)關(guān)注TN和NH3-N指標(biāo)。

（3）在本文工作基礎(chǔ)上，進(jìn)一步研究可在以下方面拓展：增加考慮水質(zhì)指標(biāo)項(xiàng)目及數(shù)據(jù)序列長(zhǎng)度，以提高各因子之間的相關(guān)性；探討其他Vine Copula函數(shù)的適用性；考慮實(shí)際變量間的結(jié)構(gòu)動(dòng)態(tài)性與Vine Copula模型的耦合關(guān)系。