徐鵬程，李 帆，張昌盛，仇建春

（1.揚(yáng)州大學(xué)水利科學(xué)與工程學(xué)院，江蘇 揚(yáng)州 225009；2.江蘇省水利建設(shè)工程有限公司，江蘇 揚(yáng)州 225009）

0 引 言

準(zhǔn)確完整的水文氣象數(shù)據(jù)對(duì)于水利工程項(xiàng)目設(shè)計(jì)、規(guī)劃和施工是至關(guān)重要的。水文氣象站網(wǎng)可以為水資源管理、水庫(kù)運(yùn)行、洪水預(yù)報(bào)等提供必要和實(shí)時(shí)的基礎(chǔ)輸入信息?，F(xiàn)代水文氣象站網(wǎng)應(yīng)在滿足不同的用戶需求同時(shí)具備高效節(jié)約的特征。因此，水文氣象站網(wǎng)優(yōu)化常常從多個(gè)目標(biāo)函數(shù)的約束條件入手［1，2］綜合考慮站點(diǎn)的布設(shè)方式。

水文氣象站網(wǎng)的優(yōu)化設(shè)計(jì)方法大致可以歸納為：①基于地統(tǒng)計(jì)學(xué)的方法［4，5］；②基于信息熵的評(píng)價(jià)方法［6-8］；③基于兩種及以上理論耦合的方法［9，10］。其中，信息熵方法在國(guó)內(nèi)外水文站網(wǎng)優(yōu)化中得到了廣泛應(yīng)用，其基本原則是通過(guò)最大化站網(wǎng)信息總量，同時(shí)最小化信息冗余度從而獲得最優(yōu)站點(diǎn)布置方案。Mishra 和Coulibaly［11］提出了用基于信息傳遞指數(shù)（Transfer Index，TI）的多元回歸模型有效量化了徑流站網(wǎng)的顯著性。Markus等［12］將信息傳遞指標(biāo)法（Direct Transmission Index，DIT）和廣義最小二乘法（Generalized Least Square，GLS） 耦合并將站網(wǎng)信息量作為關(guān)鍵的懲罰函數(shù)對(duì)徑流站網(wǎng)進(jìn)行了優(yōu)化評(píng)估。Fuentes 等［13］提出了基于信息熵-效用（Entropy-Utility）標(biāo)準(zhǔn)對(duì)環(huán)境和空氣污染型樣本分析，并在非平穩(wěn)性的假設(shè)下建立了基于信息熵框架下的貝葉斯最優(yōu)子網(wǎng)絡(luò)和空間相關(guān)模型。Li等［14］采用了最大化信息量最小化冗余信息量的準(zhǔn)則對(duì)水文徑流站網(wǎng)進(jìn)行優(yōu)化設(shè)計(jì)并考慮站網(wǎng)優(yōu)化結(jié)果對(duì)序列離散化法的敏感性分析。

近年來(lái)，由于氣候變化和人類活動(dòng)雙重影響，未來(lái)氣候變化情景下的水文氣象情勢(shì)存在很大的不確定性，另一方面降雨、洪水極端事件呈現(xiàn)出顯著的趨勢(shì)性和突變性［15-17］，為此進(jìn)行站網(wǎng)優(yōu)化時(shí)有必要關(guān)注未來(lái)氣候變化條件下水文氣象序列的趨勢(shì)特性［18］。同時(shí)由于離散化方法對(duì)于站網(wǎng)優(yōu)化結(jié)果的影響較大［13，14］，為此Xu 等［19］采用阿基米德Copula 熵進(jìn)行多維互信息的估計(jì)，盡管克服了互信息估計(jì)的不確定性，但是由于Archimedean Copula 函數(shù)［20-22］對(duì)于高維聯(lián)合分布的擬合誤差使得其對(duì)于高維冗余信息量的估計(jì)精度不夠。為此本文擬從兩個(gè)方面進(jìn)行研究：①采用C-Vine Copula 函數(shù)擬合高維的聯(lián)合分布［］進(jìn)而求得Copula 熵以提高對(duì)總相關(guān)量的估計(jì)精度；②從不同流域由于站網(wǎng)內(nèi)部水文氣象序列的趨勢(shì)性出發(fā)，加強(qiáng)不同趨勢(shì)程度引發(fā)的序列非平穩(wěn)特性對(duì)站網(wǎng)優(yōu)化結(jié)果的敏感性分析。

1 基于C-Vine Copula 熵的站網(wǎng)優(yōu)化基礎(chǔ)理論

1.1 C-Vine Copula函數(shù)

假定[X1，X2，…，Xd]是一個(gè)d個(gè)站點(diǎn)組成的站網(wǎng)，其聯(lián)合概率密度函數(shù)為p(x1，x2，…，xd)，變量X1，X2，…，和Xd對(duì)應(yīng)的邊緣密度函數(shù)分別為pX1(x1)，pX2(x2)，…，pXd(xd)。多維聯(lián)合概率密度函數(shù)可以用如下公式求得：

式中：P（·）函數(shù)是多維累計(jì)聯(lián)合分布函數(shù)（Joint Cumulative Distribution Function，JCDF）；PXi是單變量Xi的邊緣分布函數(shù)；θ是Copula參數(shù)值。

由于阿基米德Copula 在高維聯(lián)合分布模擬的不足，本文通過(guò)將多元概率密度分解為一系列二元Copula，基于C-Vine Copula 函數(shù)解決上述問(wèn)題［7］。按照變量間的相關(guān)性排序是確定Cvine Copula函數(shù)樹(shù)形結(jié)構(gòu)的關(guān)鍵步驟。n維C-Vine Copula 的概率密度函數(shù)可以表達(dá)如下：

式中：fk(xk；θk) (k=1，2，…，n)表達(dá)了第k個(gè)變量的邊緣概率密度函數(shù)；ci，i+j|1：(i-1)為二維Copula密度函數(shù)。

為了說(shuō)明C-Vine Copula函數(shù)的建立過(guò)程，同時(shí)由于篇幅的限制，本文以四維C-vine Copula 為例，圖1顯示了一個(gè)由4個(gè)變量和三棵樹(shù)組成的C-Vine結(jié)構(gòu)。為了便于說(shuō)明，數(shù)字1、2、3和4分別表示隨機(jī)變量x1，x2，x3和x4。這個(gè)圖示說(shuō)明了建立C-vine copula 的過(guò)程。在第一棵樹(shù)中，變量x1起著關(guān)鍵作用（基于Kendall秩的分析［23，24］），處于樹(shù)一的結(jié)點(diǎn)，其余變量處于樹(shù)一的邊。同理樹(shù)二中變量x2起著關(guān)鍵作用，處于樹(shù)二的結(jié)點(diǎn)，以此類推。

圖1 四維C-Vine Copula函數(shù)的樹(shù)形結(jié)構(gòu)Fig.1 The tree strcture of four-dimensional C-vine Copula

1.2 C-Vine Copula熵為核心的站網(wǎng)優(yōu)化模型

在獲得d個(gè)站點(diǎn)情形下的C-Vine Copula 函數(shù)，可以通過(guò)Copula 熵與總相關(guān)量（TC）之間互為相反數(shù)的數(shù)學(xué)關(guān)系，獲得所需站網(wǎng)優(yōu)化過(guò)程中的兩個(gè)關(guān)鍵站網(wǎng)信息目標(biāo)函數(shù)。

其中Copula熵可以表示為：

其中，c(u1，u2，…，ud) 是Copula 密度函數(shù)，主要可以采用兩種方法計(jì)算式（3）中描述的Copula 熵：①多重積分法；和②蒙特卡洛模擬法。由于多重積分法在維數(shù)較高時(shí)求解困難，本文采用了后一種方法計(jì)算Copula 熵。為此最終本文求得如下以Copula 熵為核心的兩個(gè)目標(biāo)函數(shù)（最大化聯(lián)合信息量-最小化總相關(guān)量）：

其中，邊緣熵（H(Xi)）本文擬通過(guò)極大熵原理（POME）進(jìn)行求解，而Copula 熵則是從C-Vine Copula 函數(shù)角度出發(fā)，首先估計(jì)C-Vine Copula參數(shù)，進(jìn)而采用蒙特卡洛模擬獲得隨機(jī)數(shù)的思路估計(jì)C-Vine Copula 熵；接著按照公式（4）獲得兩個(gè)目標(biāo)函數(shù)值；最終通過(guò)多目標(biāo)優(yōu)化求得帕累托最優(yōu)解集的方法獲得最優(yōu)站點(diǎn)組合。傳統(tǒng)的多維阿基米德Copula 函數(shù)只具備一個(gè)自由度的參數(shù)無(wú)法涵蓋高維情形下的所有可能的尾部依賴性，維數(shù)較高的阿基米德族Copula 函數(shù)更適合于擬合具有正相關(guān)的多變量相依性結(jié)構(gòu)，為此本文擬從C-Vine Copula角度進(jìn)行多站點(diǎn)間的相依性結(jié)構(gòu)的擬合從而獲得準(zhǔn)確的Copula熵值。

2 研究區(qū)域

作為中國(guó)最為典型的兩個(gè)代表性流域，在過(guò)去幾十年以來(lái)，其流域的水文氣象情勢(shì)和氣候變化、人類活動(dòng)存在顯著的響應(yīng)關(guān)系。本文為了對(duì)比不同趨勢(shì)程度對(duì)站網(wǎng)優(yōu)化結(jié)果的影響，分別選取了黃河流域和淮河流域43個(gè)雨量站組成的初始站網(wǎng)（圖2）并都選用了1992-2018年的日降水量觀測(cè)序列作為研究對(duì)象，考慮到Copula 函數(shù)的獨(dú)立同分布假設(shè)并對(duì)降雨數(shù)據(jù)進(jìn)行了前處理以消除零降雨值對(duì)序列自相關(guān)的影響。首先分析了全序列的優(yōu)化結(jié)果；為了加強(qiáng)優(yōu)化結(jié)果對(duì)趨勢(shì)性的敏感性分析，分別采用25年，20年，10年，5年和2年的滑動(dòng)窗口法進(jìn)行子序列的二次優(yōu)化分析。

圖2 黃河和淮河流域雨量站分布Fig.2 Distribution of rainfall stations in the Yellow River and Huaihe River Basin

3 站網(wǎng)優(yōu)化結(jié)果

3.1 基于全序列的站網(wǎng)優(yōu)化結(jié)果

由于帕累托解集的不唯一性且考慮到站點(diǎn)組合數(shù)的眾多，為了提高優(yōu)化效率，本文通過(guò)隨機(jī)生成10 000 個(gè)站點(diǎn)組合，從中選取滿足式（4）中兩個(gè)目標(biāo)函數(shù)的帕累托解集，并采用統(tǒng)計(jì)學(xué)的方法統(tǒng)計(jì)帕累托解中的各站點(diǎn)出現(xiàn)累計(jì)頻率，最終站點(diǎn)優(yōu)化結(jié)果見(jiàn)圖3。其中紅色圓圈的直徑大小代表了帕累托解中站點(diǎn)出現(xiàn)頻率的多少。由圖3可知，對(duì)于淮河流域而言，高頻率站點(diǎn)主要出現(xiàn)海拔較高的山丘區(qū)域低頻率站點(diǎn)主要存在于平原地區(qū)；而對(duì)于黃河流域而言，高頻率和低頻率的站點(diǎn)分布混合分布在流域內(nèi)并沒(méi)有呈現(xiàn)顯著的海拔分布效應(yīng)。

圖3 黃河和淮河流域全序列優(yōu)化結(jié)果Fig.3 The optimization results based on whole sequence of Yellow River and Huaihe River Basins

3.2 總相關(guān)量估計(jì)精度分析

為了對(duì)比分析阿基米德Copula 函數(shù)和C-Vine Copula 函數(shù)估計(jì)總相關(guān)量時(shí)的差異性，本節(jié)分別計(jì)算了k站點(diǎn)組成（k=1，2，3，…，43）最優(yōu)站點(diǎn)組合中的總相關(guān)量；并將離散化方法（直方圖法）計(jì)算得到總相關(guān)量作為參照樣本。具體結(jié)果見(jiàn)圖4，其中維度為10 代表10站點(diǎn)組成的最優(yōu)站點(diǎn)組合所包含的冗余信息量（即總相關(guān)量）。由圖4可知，由于C-Vine Copula 函數(shù)在高維聯(lián)合擬合的優(yōu)勢(shì)，其按照蒙特卡洛隨機(jī)模擬求得總相關(guān)量和直方圖法求得總相關(guān)量之間是比較接近的，而阿基米德Copula 函數(shù)由于在高維聯(lián)合分布擬合的不足，其總相關(guān)的估計(jì)值隨著維數(shù)越來(lái)越高越發(fā)的偏離直方圖法計(jì)算得到冗余信息量。這在一定程度上表明了采用C-Vine Copula 函數(shù)去估計(jì)總相關(guān)量這一思路是可行的且具有比阿基米德Copula 函數(shù)精度更優(yōu)的優(yōu)勢(shì)。

圖4 總相關(guān)量估計(jì)對(duì)比Fig.4 Comparison of total correlation estimation

3.3 站點(diǎn)優(yōu)化結(jié)果對(duì)序列趨勢(shì)性的敏感性分析

由于兩個(gè)流域研究降雨序列的趨勢(shì)程度的不同（圖5），可以發(fā)現(xiàn)淮河流域43 個(gè)站點(diǎn)中有41 個(gè)站點(diǎn)都通過(guò)了5%顯著水平的Mann-Kendall（MK）檢驗(yàn)，而黃河流域只有6個(gè)站點(diǎn)呈現(xiàn)顯著的上升趨勢(shì)性。為此有必要研究趨勢(shì)性程度的不同對(duì)站網(wǎng)優(yōu)化結(jié)果的影響。

圖5 黃河與淮河流域序列趨勢(shì)分析結(jié)果（MK檢驗(yàn)）Fig.5 Trend analysis results of the Yellow River and Huaihe River Basin Series （MK test）

為了量化站網(wǎng)優(yōu)化結(jié)果對(duì)趨勢(shì)性的敏感性分析，本節(jié)基于不同窗寬（Sliding Windows，SW）條件下的降雨子序列進(jìn)行站網(wǎng)的二次優(yōu)化。由圖3中結(jié)果可知，站點(diǎn)入選帕累托解集呈現(xiàn)頻率差異性，后續(xù)的敏感性分析按照入選頻率對(duì)各站點(diǎn)進(jìn)行了降序排序（秩次越小代表其越重要）。由圖6可以發(fā)現(xiàn)淮河流域在不同窗寬條件下的優(yōu)化結(jié)果存在顯著的差異性，只有在站點(diǎn)S8，S31 和S43 有一定的相似性（圖6中黃河流域3 個(gè)站點(diǎn)所在行處于藍(lán)色系，其余站點(diǎn)所在行在藍(lán)色、黃色和紅色色系間不斷波動(dòng)）。而黃河流域不同窗寬條件下站點(diǎn)的排序結(jié)果相似性較為顯著。

圖6 不同窗寬條件下黃河和淮河流域站點(diǎn)秩次圖Fig.6 Rank map of Yellow River and Huaihe River Basin station under different window width conditions

為了能夠直觀顯示站點(diǎn)排序結(jié)果的差異性，又對(duì)各窗寬條件下的子序列（秩次序列）進(jìn)行了相關(guān)性分析；如窗寬等于5年時(shí)，有23 個(gè)子序列，兩兩之間配對(duì)可以產(chǎn)生C223=253 個(gè)相關(guān)系數(shù)，并對(duì)這些子序列相關(guān)系數(shù)的經(jīng)驗(yàn)分布形態(tài)進(jìn)行了計(jì)算分析（見(jiàn)圖7）。圖7中橫向比較可以發(fā)現(xiàn)，各窗寬下淮河流域排序子序列的相關(guān)程度明顯弱于黃河流域的排序子序列?？v向比較可以發(fā)現(xiàn)，隨著窗寬的不斷減小，排序子序列的相關(guān)程度在不斷減弱。

圖7 不同窗寬條件下排序子序列的相關(guān)性分析Fig.7 Correlation analysis of rank subsequences with different window widths

總體來(lái)看，由于淮河流域的顯著趨勢(shì)性使得其站網(wǎng)的優(yōu)化結(jié)果對(duì)于窗寬的敏感性較強(qiáng)而趨勢(shì)性程度較弱的黃河流域其站網(wǎng)優(yōu)化結(jié)果對(duì)于窗寬的敏感性較弱；趨勢(shì)程度較大的淮河流域站網(wǎng)秩次的年際變化相比于趨勢(shì)程度小的黃河流域站網(wǎng)更加顯著。

4 分析與討論

基于信息熵的站網(wǎng)優(yōu)化方法常常是水文氣象站網(wǎng)優(yōu)化領(lǐng)域常用的方法，但是受限于離散化手段的高度敏感性，使得其計(jì)算過(guò)程有很大的不確定性［6，7］。前人研究表明動(dòng)態(tài)化箱寬和固定箱寬的離散化方法進(jìn)行了比較，發(fā)現(xiàn)對(duì)最佳站網(wǎng)方案的影響不容忽視［9］。受到其計(jì)算能力的限制，早期大多數(shù)研究使用擬合分布或核密度估計(jì)來(lái)估計(jì)信息熵值，但其計(jì)算精度，特別對(duì)于多變量的冗余信息量估計(jì)值效果不佳［9，10］。圖4中的直方圖法較寬的TC 估計(jì)置信區(qū)間很好地應(yīng)證了這一點(diǎn)。本研究采用C-Vine Copula 函數(shù)較為精確地模擬多站點(diǎn)間高維相依性結(jié)構(gòu)，一方面可以克服信息熵離散化方法引發(fā)的優(yōu)化結(jié)果的不確定性問(wèn)題，另一方面相比于傳統(tǒng)阿基米德Copula 函數(shù)提高了高維依賴性結(jié)構(gòu)的模擬精度，從而實(shí)現(xiàn)對(duì)站網(wǎng)優(yōu)化的兩個(gè)目標(biāo)函數(shù)：信息總量和總相關(guān)量的精準(zhǔn)估計(jì)。

此外，時(shí)間變異性和空間差異對(duì)水文氣象站網(wǎng)優(yōu)化設(shè)計(jì)的影響得到了廣泛關(guān)注［14，16］。在水文氣象站網(wǎng)優(yōu)化設(shè)計(jì)過(guò)程中，氣候變化、人類活動(dòng)和水文過(guò)程的非平穩(wěn)性不容忽視，尤其是對(duì)于高度異質(zhì)性、局部性且受地理、地形和氣候因素強(qiáng)烈影響的降雨事件。不同的時(shí)間域范圍可能會(huì)對(duì)站網(wǎng)設(shè)計(jì)結(jié)果產(chǎn)生不同的影響，這意味著最佳的站網(wǎng)布設(shè)可能只適用于特定的觀測(cè)時(shí)段。盡管前人對(duì)于站網(wǎng)優(yōu)化開(kāi)展了大量工作，但是從水文氣象時(shí)間序列的非平穩(wěn)性入手應(yīng)對(duì)氣候變化背景下水文氣象站網(wǎng)優(yōu)化的動(dòng)態(tài)分析工作是有限的。本研究著重討論了中國(guó)兩個(gè)典型的流域（淮河流域和黃河流域）由于各自對(duì)于氣候變化響應(yīng)程度的不同（黃河流域趨勢(shì)型非平穩(wěn)性較弱，淮河流域較強(qiáng)）其最終站點(diǎn)優(yōu)化秩次呈現(xiàn)不同的變化形態(tài)。因此，本研究表明在處理水文氣象站網(wǎng)優(yōu)化時(shí)應(yīng)非常謹(jǐn)慎，因?yàn)樵谄渌麠l件相同情形下，基于固定的全序列剔除或者增加某些臺(tái)站會(huì)導(dǎo)致整個(gè)站網(wǎng)系統(tǒng)的水文氣象信息損失或信息冗余。根據(jù)本文的研究結(jié)果，站網(wǎng)的優(yōu)化設(shè)計(jì)應(yīng)當(dāng)以使其更適應(yīng)不斷變化的水文氣象條件可能更為可取。

5 結(jié) 論

（1）Archimedean Copula 函數(shù)由于僅僅適用于正相關(guān)的多變量相依性結(jié)構(gòu)，特別是對(duì)于高維聯(lián)合分布的擬合精度無(wú)法準(zhǔn)確刻畫潛在的高維相依性結(jié)構(gòu)的尾部特征，最終使得其對(duì)于高維冗余信息量的估計(jì)精度明顯不夠。為此本研究擬從C-Vine Copula 函數(shù)出發(fā)開(kāi)發(fā)基于C-Vine Copula 熵的多目標(biāo)站網(wǎng)優(yōu)化模型，一方面有效提高了阿基米德Copula熵總相關(guān)量的估計(jì)精度，另一方面有效實(shí)現(xiàn)了對(duì)于水文氣象站網(wǎng)過(guò)程中最大化站網(wǎng)信息總量和最小化信息冗余量的兩個(gè)目標(biāo)函數(shù)的求解，基于兩個(gè)典型研究區(qū)域的計(jì)算結(jié)果可知，本文提出的C-Vine Copula熵具有較好的適用性和可行性。

（2）在水文氣象站網(wǎng)優(yōu)化設(shè)計(jì)過(guò)程中，氣候變化、人類活動(dòng)和水文過(guò)程的非平穩(wěn)性不容忽視。一方面趨勢(shì)性引發(fā)的序列非平穩(wěn)特征增加了未來(lái)氣候變化背景下水文氣象站網(wǎng)優(yōu)化的不確定性，另一方面不同的時(shí)間域范圍可能會(huì)對(duì)站網(wǎng)設(shè)計(jì)結(jié)果產(chǎn)生不同的影響?；诨春雍忘S河流域的計(jì)算結(jié)果可知：趨勢(shì)性較強(qiáng)的淮河流域隨著窗口期的不同，子序列的秩次年際變化明顯強(qiáng)于趨勢(shì)性較小的黃河流域。為此，在氣候變化的背景下進(jìn)行水文氣象站網(wǎng)優(yōu)化分析時(shí)，有必要考慮站點(diǎn)序列的趨勢(shì)性從而引發(fā)的非平穩(wěn)特性對(duì)站網(wǎng)優(yōu)化結(jié)果的影響。