楊 蕊，吳時強，高學平，張 晨

（1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300350；2.南京水利科學研究院 水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210029）

1 研究背景

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

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

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

2 研究方法

2.1 邊緣分布擬合構建聯合分布之前，須先確定各風險因子的邊緣分布函數。隨機變量邊緣分布的經驗頻率采用水文中常用的Gringorten公式計算［26］：

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

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

2.2 Copula函數Copula函數是定義在[0，1]n區(qū)間上的多維聯合分布函數，它可以將多個隨機變量的邊際分布連接起來構造聯合分布。

Sklar定理［28］：令H為n維聯合分布函數，F1，F2，…，Fn為其邊緣分布，那么存在唯一的Copu?la函數C，使得對?x∈R n，有

若F1，F2，…，Fn連續(xù)，則C是唯一的。相反地，如果C是一個Copula函數，F1，F2，…，Fn是分布函數，則上式定義的H是邊緣分布為F1，F2，…，Fn的一個聯合分布函數。

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

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

表1 二元Copula函數表達式及其相關性質

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

式中，下標j表示樹層級的遍歷，下標i表示每棵樹中節(jié)點的遍歷。四維變量C-vine Copula結構如圖1所示。

圖1 C-vine Copula四維變量結構圖

該方法的優(yōu)勢在于，相較于傳統多元Copula函數，其在構建聯合分布函數時不僅可以將各變量的邊緣分布函數與相關結構分開建模，并且允許多個變量兩兩之間的相關結構有所不同［21］。由于不同的Copula 函數代表著特定的相關結構，在選取Pair Copula 函數時，本文從包含大部分相關結構的15種Copula函數中，選擇最為恰當的Copula函數來刻畫變量（條件或無條件）間的真實相關結構，所構建的聯合分布函數比傳統的多元Copula函數更具靈活性［21，25］。

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

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

3 案例應用

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

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

表2 運行前/后南陽監(jiān)測站水質指標統計特征表

3.2 C-vine Copula模型構建

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

表3 運行前/后各變量擬合不同分布類型對應的AIC值

3.2.2 C-vine Copula模型的構建和模型檢驗

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

表4 運行前/后C-vine Pair Copula的函數選擇結果及參數估計

圖3 運行前/后C-vine Copula四維水環(huán)境風險因子聯合分布模型

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

圖4 運行前/后Kendall系數模型驗證

3.3 風險情景設計

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

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

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

表6 關鍵因子風險組合情景（單位：TP/TN/NH3-N（mg·L-1），Chl-a（μg·L-1））

3.4 風險識別結果分析

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

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

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

圖5 運行前/后模型（聯合分布）不同工況下的概率密度值

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

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

圖6 運行前/后關鍵因子風險組合情景聯合分布概率密度值

4 結論

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

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

（3）在本文工作基礎上，進一步研究可在以下方面拓展：增加考慮水質指標項目及數據序列長度，以提高各因子之間的相關性；探討其他Vine Copula函數的適用性；考慮實際變量間的結構動態(tài)性與Vine Copula模型的耦合關系。