杜奕良 ，涂新軍 ，杜曉霞 ，陳曉宏 ，林凱榮 ，吳海鷗

（1.中山大學 水資源與環(huán)境系， 廣東 廣州　510275；2.廣東省華南地區(qū)水安全調控工程技術研究中心， 廣東 廣州　510275）

1　研究背景

農田灌溉需水常受到農作物類型、播種面積、土壤類型、耕作制度、灌溉水利用系數(shù)和氣候等因素的影響［1］，其中氣候特別是降水屬于隨機過程的不確定性因素，而其他因素可歸為確定性因素。因此在區(qū)域水資源規(guī)劃的農田灌溉需水分析中，通常需要推求從枯到豐多個降水頻率條件下的需水。當區(qū)域空間尺度較大時，農田灌溉需水需要考慮上述因素的空間差異［2］，一般是在各分區(qū)農田灌溉需水的基礎上進行空間疊加，其中如何得到給定區(qū)域降水頻率的分區(qū)組合是個難點。也就是說，空間尺度較大的區(qū)域農田灌溉需水分析，既要符合全區(qū)降水統(tǒng)計特征，還要反映各分區(qū)降水的空間關系。

以往的處理方法采用典型年法，即選擇與某一給定年降水量頻率相近的典型年，對各分區(qū)年降水量進行同倍比縮放［3-4］。但是降水空間組合關系并非一成不變的，某一區(qū)域降水頻率下對應多種分區(qū)組合，如何選擇分區(qū)組合屬于多變量統(tǒng)計模擬和設計問題。Copula作為多維隨機變量的連接函數(shù)，因其既可以呈現(xiàn)每一個隨機變量的分布特征，同時又能反映變量之間的相依性，近些年被廣泛應用于水文多維聯(lián)合分布模擬［5-7］。其中也包括水文變量的空間組合模擬，如不同水文區(qū)的降水或徑流豐枯遭遇［8-12］，而且當確定了水文多維聯(lián)合分布和邊緣分布的理論模型及參數(shù)后，可進一步根據(jù)同頻法和最大權函數(shù)法推求變量組合［13］，如雨潮遭遇分析［14］、推求洪水過程線［15］和洪水風險評估［16］等。

本文以廣東省農田灌溉需水為例，基于高維Gaussian Copula函數(shù)構建廣東省8個農業(yè)分區(qū)的降水聯(lián)合分布模型。考慮各分區(qū)降水和灌溉定額的空間差異，采用典型年法、同頻法和最大權函數(shù)法，得到給定全省降水頻率的分區(qū)及全省農田灌溉需水，并比較3種計算方法的區(qū)域需水、分區(qū)降水頻率及需水空間組合特征與差異，為空間尺度較大的區(qū)域農田灌溉需水分析提供新思路。

2　研究方法

2.1分區(qū)降水聯(lián)合分布模擬空間尺度較大區(qū)域的分區(qū)降水存在空間差異，但也存在相似降水成因。而分區(qū)降水的年系列一般認為是隨機的。因此，可引入適用于高維的Gaussian Copula函數(shù)模擬分區(qū)降水的聯(lián)合分布［17-18］。設為d個隨機變量，其邊緣分布函數(shù)分別為FX1(x1)，…，d維Gaussian Copula的分布函數(shù)C(u1，…，ud)和密度函數(shù)c(u1，…，ud)分別為：

其中

式中：ρij∈［-1，1］，i=1，…，d，j=1，…，d為任意兩個變量之間的相關系數(shù)。

邊緣分布模擬采用4種水文領域應用較為廣泛的3參數(shù)單變量分布模型，即廣義極值（Generalized Extreme Value，GEV）、廣義 Logistic（Generalized Logistic，GLO）、皮爾遜Ⅲ型（Pearson Ⅲ，PE3）和廣義正態(tài)（Generalized Normal，GNO）。擬合檢驗采用K-S方法，檢驗統(tǒng)計量D為［19］：

式中：n為樣本長度。

多維聯(lián)合分布擬合檢驗采用基于Rosenblatt變換的Cramervon Mises方法，檢驗統(tǒng)計量為［21］：

式中：E為基于Rosenblatt變換的偽隨機變量，令E1=u1，則

式中：i=2，…，d。

邊緣和聯(lián)合分布模型參數(shù)分別采用線性矩法［20］和極大似然法估計［21］，再基于Monte Carlo方法模擬分別計算統(tǒng)計量D和的P值大小。若P＞0.05，表明該模型在顯著性檢驗水平為0.05情況下是可接受的。若有多個模型通過了統(tǒng)計檢驗，則根據(jù)均方根差RMSE和赤遲信息準則AIC進一步優(yōu)選［21］。

2.2分區(qū)降水頻率組合的推求方法區(qū)域降水頻率是指分區(qū)降水經(jīng)過面積加權平均后，擬合得到的理論分布函數(shù)。理論上給定一個區(qū)域降水頻率，存在無數(shù)種分區(qū)降水頻率組合。為此，可分別根據(jù)典型年法、同頻法和最大權函數(shù)法推求分區(qū)降水頻率組合。典型年法是指在實測樣本中選擇區(qū)域降水量較接近的年份，對降水空間組合進行同倍比縮放，進而得到分區(qū)降水頻率組合的方法。同頻法是指假定各分區(qū)降水頻率相同情況下的區(qū)域降水空間組合方法。最大權函數(shù)法是在多維聯(lián)合分布模擬基礎上提出的［22-23］，其定義為［13］：

具體分析步驟如下：

（1）根據(jù)樣本確定的理論聯(lián)合分布模型，采用Monte Carlo法模擬生成l組分區(qū)降水頻率組合［u1(j)，…，ui(j)，…，ud(j)］，j=1，…，l，計算分區(qū)降水量組合為［X1(j)，…，Xi(j)，…，Xd(j)］，j=1，…，l及區(qū)域降水量X0(j)，j=1，…，l與其降水頻率u0(j)，j=1，…，l。

（4）考慮到Monte Carlo法生成分區(qū)降水頻率組合的隨機性，對于模擬優(yōu)選的分區(qū)降水頻率組合進一步給出置信區(qū)間。即重復分析步驟（1）—（3）B次，每個給定區(qū)域降水頻率得到B組分區(qū)降水頻率組合B和相應的聯(lián)合分布概率密度則1-α置信區(qū)間的組合S()fα應滿足：

式中：fα為（1-α）分位數(shù)，

在B個組合中，聯(lián)合概率密度值最大的組合即為最終的選擇。

另外為了盡可能擴大模擬數(shù)和提高誤差控制精度，本文研究l、Ru0、B和檢驗水平α分別取106、10-4、103和0.05。

2.3農田灌溉需水計算灌溉農業(yè)受氣候、土壤、植被等自然條件的制約，土壤和植被等條件是相對穩(wěn)定的，但降水卻時刻在變化，是個隨機因素。因此，農田灌溉需水是與降水有關的隨機過程。水資源規(guī)劃中分析農田灌溉需水，通常是明確灌溉農田面積A和灌溉水利用系數(shù)η，并給出不同降水頻率下單位面積的農田灌溉凈定額q(u)。因此，降水頻率u的農田灌溉需水W(u)計算公式為：

對于空間尺度較大的區(qū)域來說（如省級行政區(qū)），認為區(qū)域內不同分區(qū)（如水資源分區(qū)、行政分區(qū)或農業(yè)分區(qū)等）的降水和灌溉定額是有差異的。因此，整個區(qū)域降水頻率u0的農田灌溉需水是各分區(qū)的農田灌溉需水之和。

3　實例研究

3.1區(qū)域概括與數(shù)據(jù)廣東省位于中國大陸最南端，全省面積17.96萬km2。全省屬于熱帶、南亞熱帶季風氣候區(qū)，雨量充沛，但時空分布不均，旱澇常見。自1970年代末中國改革開放以來，隨著當?shù)厣鐣?jīng)濟快速發(fā)展，區(qū)域用水量不斷增加。至2015年，全省總用水量443.1億m3，其中農業(yè)用水227.0億m3，且80%以上的農業(yè)用水為農田灌溉用水。

依據(jù)氣候條件、土壤類型、耕作制度及作物布局等，廣東省分為8個農業(yè)分區(qū)，見圖1。目前，廣東省農田灌溉基本采用一年兩熟的耕作制度。根據(jù)灌溉試驗和歷史資料［24］，各農業(yè)分區(qū)農田灌溉凈定額，其理論統(tǒng)計分布（P-Ⅲ分布）參數(shù)特征值如表1所示。廣東省在實行最嚴格水資源管理制度時，要求2020年各地市農田灌溉水有效利用系數(shù)介于0.51～0.55，因此本文灌溉水有效利用系數(shù)采用全省均值0.511。各分區(qū)農田灌溉面積則參照廣東省水資源規(guī)劃的2020年水田有效灌溉面積目標（見表1）。年降水量采用1953—2013年國家地面氣候站點的逐日降水量過程，其中廣東省內共25個水文站點，基于站點算術平均法得出各分區(qū)年降水量系列，其基本特征如表2。

全省年降水均值約為1819 mm。其中，高值區(qū)為分區(qū)V，均值為2789 mm，低值區(qū)為分區(qū)Ⅷ和Ⅲ，均值不足1600 mm；其他分區(qū)均值介于1600～2000 mm之間。分區(qū)最大年降水可達4071 mm（V，1973年）最小年份不足800 mm（Ⅵ，1963年；Ⅷ，1955）。各分區(qū)年降水變差系數(shù)介于0.146～0.221，全省則為0.153。

圖1　廣東省雨量站點分布及農業(yè)區(qū)區(qū)劃圖

表1　分區(qū)土地面積、灌溉面積及凈灌溉定額參數(shù)

表2　年降水量基本特征

3.2降水分布特征分區(qū)及全省降水分布擬合檢驗結果如表3。所有單變量模型擬合P值均大于顯著性檢驗水平0.05，表明均可作為年降水分布的備選模型。從RMSE和AIC值來看，模型之間差別非常小?；谧钚≈禍蕜t，得出分區(qū)Ⅰ—Ⅷ及全省的最優(yōu)模型分別為GLO、GEV、GEV、GNO、GNO、GNO、GEV、GEV、GNO，分布參數(shù)見表4。

分區(qū)降水的相依性結構如圖2所示，對角線為分區(qū)編號，左下三角的數(shù)值為兩兩之間的相關系數(shù)ρ，右上三角區(qū)域為兩兩之間的平滑擬合線和置信橢圓［25］。橢圓越圓相依性越弱，越扁相依性越強?？傮w上來看，分區(qū)Ⅰ—Ⅵ兩兩之間、分區(qū)Ⅶ與Ⅱ、Ⅲ、Ⅴ、Ⅷ之間的置信橢圓較扁，相關系數(shù)均大于0.5，降水相依性較強；分區(qū)Ⅷ與Ⅰ—Ⅵ的置信橢圓較圓，相依性相對較弱。表明廣東省東部、北部及珠江三角洲地區(qū)之間的年降水量相依性強，而西部與東部、珠江三角洲地區(qū)的相依性較弱。

如圖3所示，8個分區(qū)的年降水聯(lián)合分布擬合P值為0.258，大于顯著性檢驗水平0.05，經(jīng)驗和理論聯(lián)合分布的Q-Q點據(jù)基本位于等值線附近，表明Gaussian Copula能夠很好地模擬廣東省分區(qū)年降水聯(lián)合分布。

表3　年降水分布擬合檢驗

3.3分區(qū)降水頻率組合在水資源規(guī)劃中，為了反映特枯、枯水、平水、豐水、特豐的降水條件，通常關注降水頻率分別為0.10、0.25、0.50、0.75和0.90的農田灌溉需水。因此將以上頻率作為設定的全省降水頻率，對應的典型年份分別為1958、1989、1953、1994和1961年。3種方法推求的分區(qū)降水頻率組合如圖4。

設定的全省頻率由小到大依次對應的分區(qū)同頻頻率，分別為0.15、0.30、 0.51、0.71、 0.85。在平水時，同頻分區(qū)頻率與全省頻率差別甚小，而在特枯、枯水、豐水、特豐時，兩者相差均約為0.05。其中偏枯時（特枯、枯水）分區(qū)頻率大于全省頻率，偏豐時（豐水、特豐）分區(qū)頻率小于全省頻率。

表4　年降水理論分布模型及參數(shù)

圖2　分區(qū)降水相依性結構

典型年法中分區(qū)之間的頻率差異明顯。當設定的全省降水頻率為0.25、0.5和0.75時，約一半地區(qū)低于同頻頻率，而另一半地區(qū)則高于同頻頻率。全省降水頻率為0.1的特枯年時，分區(qū)Ⅰ、Ⅱ、Ⅴ和Ⅶ的降水頻率更小，但分區(qū)Ⅷ降水頻率較大，屬于豐水年。全省降水頻率為0.9的特豐年時，分區(qū)Ⅰ—Ⅵ的降水頻率更大，但分區(qū)Ⅶ和Ⅷ降水頻率較小，屬于枯水年。

圖3　分區(qū)降水聯(lián)合分布擬合Q-Q圖

最大權函數(shù)法推求的結果顯示，大部分分區(qū)降水頻率處于同頻位置附近，只有分區(qū)Ⅷ的降水頻率在特枯和偏豐時明顯偏離同頻位置。在全省降水頻率為0.5時（平水年），最大權函數(shù)法分區(qū)頻率置信區(qū)間的上下邊界，相對于同頻位置基本對稱。在全省降水頻率較小時（如枯水年、特枯年），置信區(qū)間偏向于大于同頻頻率。在全省降水頻率較大時（如豐水年、特豐年），置信區(qū)間偏向于小于同頻頻率。分區(qū)Ⅷ的置信區(qū)間范圍明顯大于其他分區(qū)，表明組合推求過程中分區(qū)Ⅷ不確定性更大。其主要原因是由于分區(qū)Ⅷ和其他分區(qū)的降水相依性相對較小。

圖4　分區(qū)降水頻率組合

表5　廣東省農田灌溉需水及分區(qū)組合

3.4農田灌溉需水全省農田灌溉需水結果（見表5）顯示，3種方法之間的全省灌溉需水相差甚小。全省降水頻率分別為0.1、0.25、0.5、0.75和0.9時，基于同頻法的灌溉需水分別為198.4億m3、177.3億m3、156.5億m3、138.3億m3和123.9億m3。與同頻法比較，典型年法的需水相對變化分別為-2.5%、1.1%、-1.5%、0.8%和4.0%，最大權函數(shù)法的需水相對變化分別為-0.7%、0.1%、0.0%、0.5%和0.8%，其95%置信區(qū)間邊界的需水相對變化分別為-2.0%～0.9%、-1.5%～1.3%、-1.0%～1.5%、-0.7%～2.0%和-0.6%～2.3%。

與同頻法比較，典型年法的部分分區(qū)需水變化較大，如在全省降水頻率為0.1（特枯年）和0.9（特豐年）時，分區(qū)Ⅷ需水相對變化分別達到-26.6%和37.5%，分區(qū)Ⅶ需水相對變化在全省特豐年時則達到了44.9%。相比較而言，最大權函數(shù)法的分區(qū)Ⅰ—Ⅶ需水變化較小，分區(qū)Ⅷ需水變化相對大些，前者的相對變化介于-5.1%～5.9%之間，后者約為-6.3%～14.0%。分區(qū)Ⅰ—Ⅷ的95%置信區(qū)間兩端邊界的需水變化差值，分別為19.4%～20.8%、18.6%～20.7%、18.1%～21.0%、25.2%～29.7%、18.5%～23.5%、30.1%～39.2%、25.9%～32.4%、30.7%～44.0%，顯示分區(qū)Ⅷ需水在置信區(qū)間的取值范圍明顯大于其他分區(qū)。

4　結論

廣東省東部、北部及珠江三角洲地區(qū)之間的年降水量相依性強，而西部與東部、珠江三角洲地區(qū)的相依性較弱。Gaussian Copula函數(shù)能夠很好地模擬廣東省8個農業(yè)分區(qū)年降水聯(lián)合分布，分區(qū)及全省年降水最優(yōu)分布主要為廣義極值分布和廣義正態(tài)分布。

典型年法推求的分區(qū)降水頻率差異明顯。最大權函數(shù)法推求的大部分分區(qū)降水頻率處于同頻位置附近，其中分區(qū)Ⅷ的降水頻率在特枯和特豐時明顯偏離同頻位置，而且分區(qū)Ⅷ的置信區(qū)間范圍明顯大于其他分區(qū)，表明分區(qū)Ⅷ不確定性更大。

3種推求方法之間的全省農田灌溉需水相差甚小，但是分區(qū)需水差別相對較大。與同頻法比較，典型年法和最大權函數(shù)法的全省需水相對變化絕對值分別在4.0%和0.8%之內，最大權函數(shù)法95%置信區(qū)間邊界的全省需水相對變化絕對值不超過2.3%；典型年法和最大權函數(shù)的分區(qū)（如分區(qū)Ⅶ和分區(qū)Ⅷ）需水相對變化絕對值分別可達44.9%和14.0%，最大權函數(shù)法95%置信區(qū)間邊界的分區(qū)需水變化差值，最高可達44.0%。

總之，若只要計算大尺度區(qū)域農田灌溉需水，3種方法均可以采用。但是如果需要進一步給出灌溉需水分區(qū)組合，基于高維Copula函數(shù)的聯(lián)合分布模擬分析，既考慮分區(qū)降水的獨立分布，又考慮了他們之間的相依性，而且能夠給出分區(qū)需水的置信范圍。因此可以認為，最大權函數(shù)法更為合理。

參考文獻：

［1］宋悅，粟曉玲.變化環(huán)境下涇惠渠灌區(qū)凈灌溉需水的響應［J］.節(jié)水灌溉，2015（7）：90-94.

［2］顧世祥，何大明，崔遠來，等.近50多年來瀾滄江流域農業(yè)灌溉需水的時空變化［J］.地理學報，2010，65（11）：1355-1362.

［3］劉冬梅，張京恩.淺論灌溉水庫設計訴典型年法和長系列法［J］.泥沙研究，2000（2）：72-75.

［4］邵孝候，戴琳，鐘華，等.不同水平年烤煙灌溉制度優(yōu)化設計［J］.中國農村水利水電，2008（5）：22-24.

［5］涂新軍，陳曉宏，趙勇，等.變化環(huán)境下東江流域水文干旱特征及缺水響應［J］.水科學進展，2016，26（6）：810-821.

［6］涂新軍，陳曉宏，刁振舉，等.珠江三角洲Copula徑流模型及西水東調缺水風險分析［J］.農業(yè)工程學報，2016，32（18）：162-168.

［7］高超，梅亞東，涂新軍.基于Copula函數(shù)的區(qū)域降水聯(lián)合分布與特征分析［J］.水電能源科學，2013（6）：1-5.

［8］謝華，羅強，黃介生 .基于三維Copula函數(shù)的多水文區(qū)豐枯遭遇分析［J］.水科學進展，2012，23（2）：186-193.

［9］莫淑紅，沈冰，張曉偉，等.基于Copula函數(shù)的河川徑流豐枯遭遇分析［J］.西北農林科技大學學報（自然科學版），2009，37（6）：131-136.

［10］馮平，牛軍宜，張永，等.南水北調西線工程水源區(qū)河流與黃河的豐枯遭遇分析［J］.水利學報，2010，41（8）：900-907.

［11］閆寶偉，郭生練，肖義.南水北調中線水源區(qū)與受水區(qū)降水豐枯遭遇研究［J］.水利學報，2007，38（10）：1178-1185.

［12］TU X J，DU X X，SINGH V P，et al.Joint risk of interbasin water transfer and impact of the window size of sampling low flows under environmental change［J］.Journal of Hydrology，2017，554.doi：10.1016/j.jhydrol.2017.08.037.

［13］SALVADORI G，MICHELE C D，DURANTE F.On the return period and design in a multivariate framework［J］.Hydrology and Earth System Sciences，2011，15（11）：3293-3305.

［14］涂新軍，杜奕良，陳曉宏，等 .濱海城市雨潮遭遇聯(lián)合分布模擬與設計［J］.水科學進展，2017，27（1）：49-58.

［15］李天元，郭生練，閆寶偉，等.基于多變量聯(lián)合分布推求設計洪水過程線的新方法［J］.水力發(fā)電學報，2013，32（3）：10-14.

［16］黃強，陳子燊.基于二次重現(xiàn)期的多變量洪水風險評估［J］.湖泊科學，2015，27（2）：352-360.

［17］RENARD B，LANG M.Use of a Gaussian copula for multivariate extreme value analysis：Some case studies in hydrology［J］.Advances in Water Resources，2007，30（4）：897-912.

［18］GENEST C，F(xiàn)AVRE A，BéLIVEAU J，et al.Metaelliptical copulas and their use in frequency analysis of multi-variate hydrological data［J］.Water Resources Research，2007，43（9）：223-236.

［19］TU X J，SINGH V P，CHEN X H，et al.Uncertainty and variability in bivariate modeling of hydrological droughts［J］.Stochastic Environmental Research Risk Assessment，2016，30（5）：1317-1334.

［20］HOSKING J R M.L-Moments：analysis and estimation of distributions using linear combinations of order statistics［J］.Journal of the Royal Statistical Society，1990，52（1）：105-124.

［21］GENEST C，RéMILLARD B，BEAUDOIN D.Goodness-of-fit tests for copulas：A review and a power study［J］.Insurance Mathematics&Economics，2009，44（2）：199-213.

［22］VOLPI E，F(xiàn)IORI A.Design event selection in bivariate hydrological frequency analysis［J］.Hydrological Sciences Journal，2012，57（8）：1506-1515.

［23］李天元，郭生練，劉章君，等.基于峰量聯(lián)合分布推求設計洪水［J］.水利學報，2014，45（3）：269-276.

［24］廣東省水利水電科學研究所.廣東省一年三熟灌溉定額［M］.廣州：暨南大學出版社，1999.

［25］MURDOCH D J，CHOW E D.A graphical display of large correlation matrices［J］.American Statistician，1996，50（2）：178-180.