盧澤華,王君勤

(1.四川省農(nóng)田水利局,成都 610015;2.四川省水利科學(xué)院研究院,成都 610072;3.四川大學(xué) 水利水電學(xué)院,成都 610041)

隨著全球氣候變化和社會經(jīng)濟發(fā)展，水資源短缺問題日益嚴(yán)峻，已逐漸成為人類生存與發(fā)展的巨額的“生態(tài)資源赤字”[1]。2016年全國用水總量6 040.20億m3，其中農(nóng)業(yè)用水(3 768.00億m3)占比62.38%，灌溉用水占農(nóng)業(yè)用水的90%～95%，農(nóng)業(yè)是我國的用水大戶[2]。我國農(nóng)業(yè)水資源和農(nóng)業(yè)水管理面臨巨大挑戰(zhàn)，對水資源系統(tǒng)進行有效的風(fēng)險管理已成為水資源科學(xué)發(fā)展的必然趨勢[3]。農(nóng)業(yè)水資源短缺風(fēng)險評估作為農(nóng)業(yè)水資源風(fēng)險管理的基礎(chǔ)，可為區(qū)域水資源優(yōu)化配置提供基礎(chǔ)依據(jù)，對區(qū)域種植結(jié)構(gòu)調(diào)整，干旱風(fēng)險應(yīng)對具有重要意義[4]。

水資源短缺風(fēng)險是指在特定的時空環(huán)境下，由于來水和用水存在的不確定性，使區(qū)域水資源系統(tǒng)發(fā)生供水短缺的概率及由此產(chǎn)生的損失[5]。水資源系統(tǒng)是一個復(fù)雜的大系統(tǒng)，廣泛存在隨機性和模糊性[6]。目前，對水資源短缺風(fēng)險的研究大多從隨機模型或模糊模型的角度探討水資源短缺風(fēng)險問題。羅軍剛等[3]構(gòu)建了基于熵權(quán)的水資源短缺風(fēng)險評價模型；王紅瑞等[6]構(gòu)建了基于模糊概率的水資源短缺風(fēng)險評價模型；Feng等[7]構(gòu)建了基于信息擴散理論的水資源短缺風(fēng)險評價模型；Hsieh等[8]構(gòu)建了基于多站點徑流時空隨機模擬的灌溉用水短缺風(fēng)險評價模型；Moursi等[9]構(gòu)建了基于氣候響應(yīng)函數(shù)和大氣環(huán)流模式的未來氣候變化下農(nóng)業(yè)水資源短缺風(fēng)險評價模型。以上模型具有較好的數(shù)理基礎(chǔ)，評價指標(biāo)意義明確，但沒有充分考慮水資源系統(tǒng)組分間的非線性關(guān)系，無法對區(qū)域農(nóng)業(yè)供水量和作物需水量之間復(fù)雜的非線性關(guān)系定量描述。

基于變量間非線性相關(guān)關(guān)系構(gòu)建的Copula函數(shù)能獨立于隨機變量的邊緣分布反映變量間的相關(guān)性，可構(gòu)造任意邊緣分布的聯(lián)合分布函數(shù)，且在轉(zhuǎn)換過程中不會改變原始隨機變量，在干旱風(fēng)險分析、洪水遭遇風(fēng)險分析方面已得到廣泛應(yīng)用[10-11]，也有學(xué)者利用Copula函數(shù)進行水資源短缺風(fēng)險分析研究[12-13]，如丁志宏等[12]構(gòu)建了基于Copula函數(shù)的寧夏衛(wèi)寧灌區(qū)降水量(P)和參考作物騰發(fā)量(ET0)年際聯(lián)合分布模型；Zhang等[13]構(gòu)建了基于Copula函數(shù)的河南省陸渾灌區(qū)自然降水條件下的水資源短缺風(fēng)險模型。然而有效降水量(Pe)和作物需水量(ETc)是描述灌溉系統(tǒng)的基本變量[14]，二者的演變特征與匹配程度直接影響農(nóng)業(yè)水資源管理策略，比P和ET0更直觀反映自然降水條件下農(nóng)業(yè)水資源的供需平衡。關(guān)中平原地處內(nèi)陸、遠離海洋，其年降水量僅500～750 mm，對水資源短缺的響應(yīng)尤為敏感，然而目前對關(guān)中平原Pe和ETc相關(guān)關(guān)系的研究較少，且大多是闡釋其年內(nèi)、年際變化和簡單的相關(guān)關(guān)系[15-16]，對其非線性耦合關(guān)系還未報道。本文主要在現(xiàn)有研究基礎(chǔ)上，首次運用Copula函數(shù)構(gòu)建關(guān)中平原Pe和ETc的二維聯(lián)合概率分布模型，應(yīng)用該模型分析Pe和ETc的條件聯(lián)合概率和條件回歸周期，為關(guān)中平原灌溉農(nóng)業(yè)發(fā)展規(guī)劃制定、作物種植制度優(yōu)化及區(qū)域農(nóng)業(yè)水資源管理等提供理論依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)來源

關(guān)中平原(106°48′—110°36′E，33°35′—35°51′N)地處陜西省中部的渭河流域，介于秦嶺和渭北北山之間，西起寶雞，東至潼關(guān)，長約350 km，平均海拔約500 m，面積約3.6萬km2，屬暖溫帶半干旱氣候，是陜西省糧棉油主要產(chǎn)區(qū)[17-18]。根據(jù)降水量與氣候特征，關(guān)中平原可分為關(guān)中東部和關(guān)中西部[19-20]，其中關(guān)中東部年降水量500～650 mm，關(guān)中西部年降水量550～700 mm。關(guān)中平原糧食生產(chǎn)以一年兩熟的冬小麥和夏玉米輪作為主，冬小麥一般于10月初播種，翌年6月上旬收獲，夏玉米在冬小麥?zhǔn)斋@后即時播種，當(dāng)年9月底收獲[21]。本文以冬小麥播種—夏玉米收獲(10月—翌年9月底)作為一個單位年，計算單位年有效降水量(Pey)和作物需水量(ETcy)。

結(jié)合資料系列的完整性，選取關(guān)中平原6個代表性站點1962—2016年56 a逐日氣象數(shù)據(jù)作為基本資料，各站點氣象概況見表1。氣象數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http:∥data.cma.cn/)，包括日最高、最低和平均氣溫(Tmax，Tmin，Tmean)、相對濕度(RH)、日照時數(shù)(n)和距地面2 m高處的風(fēng)速u2(計算時采用FAO推薦方法由10 m風(fēng)速換算得出風(fēng)速u2)，以及站點經(jīng)緯度和海拔高度。少數(shù)的缺失數(shù)據(jù)(占全部數(shù)據(jù)的0.56%)采用線性內(nèi)插法和多年平均值法補全。

2 研究方法

2.1 反距離權(quán)重法

反距離權(quán)重法IDW(inverse distance weighted)是一種充分考慮各因素之間地域性聯(lián)系的空間展布方法，因其原理簡便，結(jié)果精確，已被廣泛應(yīng)用[22]。關(guān)中平原地形地貌復(fù)雜，不同站點海拔差異明顯，所以直接將區(qū)域內(nèi)各站點Pey和ETcy取平均或選擇典型站點來代表區(qū)域Pey和ETcy具有明顯局限性。為增強選擇站點的代表性，本文將關(guān)中東部和關(guān)中西部各站點經(jīng)度、緯度、海拔做平均，得到2個虛擬站點，根據(jù)反距離權(quán)重法分別計算各站點對應(yīng)虛擬站點的權(quán)重，將各站點氣象因子乘以權(quán)重后相加，得到2個虛擬站點56 a逐日氣象資料，再計算2個虛擬站點的Pey和ETcy，分別代表關(guān)中東部和關(guān)中西部的Pey和ETcy，權(quán)重計算公式如下：

(1)

(2)

式中:xi，yi，zi分別為各站點的經(jīng)度(°)、緯度(°)和海拔(m);xm,ym,zm分別為虛擬站點的經(jīng)度、緯度、海拔;ri為第i個站點到虛擬站點的距離;n為區(qū)域站點個數(shù);wi為第i個站點的權(quán)重。

2.2 有效降水量(Pey)和作物需水量(ETcy)

關(guān)中平原Pey的估算參考段愛旺[23]的研究成果，ETcy采用單作物系數(shù)法計算[17]。其中，關(guān)中平原冬小麥和夏玉米生育期的劃分及Kc的選取參考康紹忠[17]和時學(xué)雙等[19]的研究成果。

2.3 構(gòu)建Pey和ETcy的邊緣分布函數(shù)

采用Copula函數(shù)進行Pey和ETcy2個特征變量的聯(lián)合之前，首先要確定其邊緣分布函數(shù)，同時要考慮變量間的相關(guān)性。設(shè)X和Y分別代表Pey和ETcy，其邊緣分布函數(shù)分別為FX(x)和FY(y)，選用常見的5種單變量分布函數(shù)(正態(tài)分布、伽馬分布、對數(shù)正態(tài)分布、廣義極值分布及韋爾布分布)[24-25]分別對其進行擬合，利用K-S檢驗確立最優(yōu)邊緣分布函數(shù)。同時采用Kendall秩相關(guān)系數(shù)(τ)、Spearman秩相關(guān)系數(shù)(ρ)和Pearson相關(guān)系數(shù)(γ)度量關(guān)中東部和關(guān)中西部1962—2016年P(guān)ey和ETcy的相關(guān)性。本文中單變量分布函數(shù)及聯(lián)合分布函數(shù)中所含參數(shù)均采用極大似然法進行估計[24]。

2.4 構(gòu)建Pey和ETcy的Copula聯(lián)合概率分布

Copula函數(shù)主要基于變量間相關(guān)性進行邊緣變量聯(lián)合，根據(jù)其定義可知Pey和ETcy的聯(lián)合分布函數(shù)F(x，y)為：

F(x,y)=P(X≤x,Y≤y)=C[FX(x),FY(y)]

(3)

常見的Copula函數(shù)有Archimedean Copula，橢圓Copula，Plackett Copula以及經(jīng)驗Copula函數(shù)。其中，Archimedean Copula函數(shù)結(jié)構(gòu)簡單，計算簡便，可以構(gòu)造出形式多樣、適應(yīng)性強的多變量聯(lián)合分布函數(shù)，能滿足多領(lǐng)域應(yīng)用要求，在實際應(yīng)用中占有重要地位。目前常用的Archimedean Copula函數(shù)有Frank，Ali-Mikhail-Haq，Clayton和Gumbel-Hougaard Copula函數(shù)，常見的相關(guān)性測度τ與Archimedean Copula函數(shù)的參數(shù)θ之間保持著一種對應(yīng)關(guān)系，具體形式見參考文獻[11]。其中τ是可以描述變量之間非線性相關(guān)性的Kendall相關(guān)系數(shù)。由參考文獻中Copula函數(shù)參數(shù)θ和Kendall秩相關(guān)系數(shù)τ間關(guān)系式可知：Clayton Copula和Gumbel-Hougaard Copula函數(shù)秩相關(guān)系數(shù)τ的范圍為τ>0，Ali-Mikhail-Haq Copula函數(shù)秩相關(guān)系數(shù)τ的范圍為-0.1817≤τ≤0.3333，F(xiàn)rank Copula函數(shù)秩相關(guān)系數(shù)τ的范圍為τ∈R。本文基于Pey和ETcy兩個特征變量，分別采用參考文獻中Copula函數(shù)擬合這兩個變量的聯(lián)合分布函數(shù)，并采用K-S檢驗對Copula函數(shù)進行擬合檢驗[10]。

2.5 基于Pey和ETcy的條件概率和回歸周期

本文主要考慮在Pey分別處于高(p≤37.5%)、中(37.5%

(4)

(5)

同理，可以計算ETcy分別處于高(p≤37.5%)、中(37.5%

3 結(jié)果與分析

3.1 基于Copula函數(shù)的農(nóng)業(yè)水資源短缺風(fēng)險模型

3.1.1Pey和ETcy的邊緣分布函數(shù)確定 選用5種單變量分布函數(shù)分別對關(guān)中東部、關(guān)中西部Pey和ETcy進行擬合，利用K-S檢驗確立最優(yōu)邊緣分布函數(shù)，K-S檢驗結(jié)果見表2，可以看出5種單變量分布函數(shù)對關(guān)中東部、關(guān)中西部Pey和ETcy的擬合優(yōu)度均達顯著水平(α<0.05)。其中，關(guān)中東部擬合最優(yōu)的分布函數(shù)分別是對數(shù)正態(tài)分布和廣義極值分布，關(guān)中西部擬合最優(yōu)的分布函數(shù)分別是伽馬分布和廣義極值分布。

表2 5種分布函數(shù)對關(guān)中平原Pey和ETcy擬合的K-S檢驗結(jié)果

關(guān)中東部、關(guān)中西部Pey和ETcy邊緣分布函數(shù)的具體參數(shù)值見表3，可以看出關(guān)中東部Pey多年平均值小于關(guān)中西部，而ETcy多年平均值大于關(guān)中東部；關(guān)中東部、關(guān)中西部均表現(xiàn)為Pey的變異系數(shù)大于ETcy，即Pey相比ETcy變化更不穩(wěn)定；關(guān)中東部、關(guān)中西部Pey和ETcy的經(jīng)驗累積概率與理論累積概率擬合較好，均呈極顯著相關(guān)(α<0.01)，進一步表明邊緣分布及參數(shù)的選擇較為合理。劉俊民等[20]研究發(fā)現(xiàn)關(guān)中平原P的多年變幅較大，薛璐[26]研究發(fā)現(xiàn)關(guān)中平原ET0年際變化具有較好的穩(wěn)定性，均與本文結(jié)果一致。

表3 關(guān)中平原Pey和ETcy邊緣分布的參數(shù)值、經(jīng)驗概率與理論概率的決定系數(shù)

3.1.2Pey和ETcy的二維聯(lián)合概率分布 采用Kendall秩相關(guān)系數(shù)τ分析度量關(guān)中東部、關(guān)中西部Pey和ETcy間的相關(guān)性，根據(jù)τ與θ的關(guān)系，計算Copula函數(shù)的參數(shù)θ，并對Copula函數(shù)進行擬合優(yōu)度評價，結(jié)果見表5?？梢钥闯鲫P(guān)中東部、關(guān)中西部Pey和ETcy均呈負相關(guān)，τ分別為-0.403和-0.228，因此僅Frank Copula函數(shù)適合描述Pey和ETcy間的相關(guān)關(guān)系。對關(guān)中東部、關(guān)中西部Pey和ETcy聯(lián)合概率分布的Frank Copula函數(shù)進行K-S檢驗，取顯著性水平α=0.05，n=55時，查柯爾莫格洛夫檢驗分位數(shù)表得出對應(yīng)分位點D0為0.179 8，表5中檢驗統(tǒng)計量D小于分位點0.179 8，表明擬合優(yōu)度均達顯著水平(α<0.05)，Pey和ETcy聯(lián)合分布的經(jīng)驗累積概率與理論累積概率擬合較好，R2分別為0.976，0.987，且均呈極顯著相關(guān)(α<0.01)，因此關(guān)中東部、關(guān)中西部Pey和ETcy的聯(lián)合分布擬合效果較好，即選用Frank Copula函數(shù)是合理的。

3.2 農(nóng)業(yè)水資源短缺風(fēng)險模型應(yīng)用

關(guān)中平原不同頻率的Pey和ETcy閾值見表5，可以看出關(guān)中東部的Pey小于關(guān)中西部的Pey，而關(guān)中東部的ETcy高于關(guān)中西部ETcy，因此關(guān)中東部的農(nóng)業(yè)水資源短缺程度明顯高于關(guān)中西部。

表4 關(guān)中平原Pey和ETcy間的相關(guān)性度量及Copula函數(shù)的參數(shù)和擬合優(yōu)度評價

表5 關(guān)中平原不同頻率的Pey和ETcy閾值

圖1給出了ETcy分別處于高(p≤37.5%)、中(37.5%

圖1 關(guān)中平原年作物需水量(ETcy)處于不同水平時年有效降水量(Pey)小于等于某特定值的條件概率和條件回歸周期

圖2 關(guān)中平原年有效降水量(Pey)處于不同水平時年作物需水量(ETcy)大于等于某特定值的條件概率和條件回歸周期

Pey分別為頻率37.5%，62.5%，87.5%的閾值時，由圖1可知，關(guān)中東部的ETcy在處于高(p≤37.5%)、中(37.5%

ETcy分別為頻率62.5%，37.5%，12.5%的閾值時，由圖2可知，關(guān)中東部Pey分別處于高(p≤37.5%)、中(37.5%

因此，關(guān)中東部、關(guān)中西部均表現(xiàn)為當(dāng)ETcy處于高水平(p≤37.5%)時，Pey不超過某一特定值的條件概率最大，重現(xiàn)期最短；當(dāng)Pey處于低水平(p≥62.5%)時，ETcy超過某一特定值的條件概率最大，重現(xiàn)期最短，即當(dāng)ETcy處于高水平(p≤37.5%)或Pey處于低水平(p≥62.5%)時，關(guān)中平原自然降水和作物需水處于不協(xié)調(diào)狀況的可能性較高，供水不能滿足需水要求的概率較大，重現(xiàn)期較短(1～4 a)，農(nóng)業(yè)水資源短缺風(fēng)險較高。

因此關(guān)中東部和關(guān)中西部根據(jù)該區(qū)域的自然降水情況，有效地調(diào)節(jié)作物種植結(jié)構(gòu)，充分利用區(qū)域自然水資源，以水資源總量控制作物布局，以水布局產(chǎn)業(yè)發(fā)展，建立與水資源承載能力相適應(yīng)、與節(jié)水增收目標(biāo)相配套的種植業(yè)結(jié)構(gòu)。此外應(yīng)該在考慮關(guān)中東部和西部灌區(qū)的特點基礎(chǔ)之上，推進農(nóng)村特色產(chǎn)業(yè)區(qū)域化布局、規(guī)模化發(fā)展、產(chǎn)業(yè)化經(jīng)營，形成產(chǎn)業(yè)聚集效應(yīng)和發(fā)展的比較優(yōu)勢，促進高效節(jié)水產(chǎn)業(yè)的快速發(fā)展。

本文對Pey和ETcy兩個特征變量進行分析只適用于描述自然降水條件下的農(nóng)業(yè)水資源短缺風(fēng)險，增加灌溉用水量、作物產(chǎn)量等特征變量，可進一步建立針對作物不同生育期不同環(huán)境下的農(nóng)業(yè)水資源短缺風(fēng)險，自然降水和灌溉水量組合條件下的農(nóng)業(yè)水資源短缺風(fēng)險等也是今后研究的主要內(nèi)容。

4 結(jié) 論

(1)關(guān)中東部和關(guān)中西部Pey和ETcy擬合最優(yōu)的邊緣分布函數(shù)為正態(tài)分布和廣義極值分布，關(guān)中西部Pey和ETcy擬合最優(yōu)的邊緣分布函數(shù)為伽馬分布和廣義極值分布，F(xiàn)rank Copula函數(shù)可較好反映Pey和ETcy間的聯(lián)合分布特性。

(2)關(guān)中東部的農(nóng)業(yè)水資源短缺程度明顯高于關(guān)中西部，當(dāng)Pe處于低水平(p≥62.5%)或ETc處于高水平(p≤37.5%)條件下，關(guān)中平原自然降水不能滿足作物需水的概率較大，重現(xiàn)期較短(1～4 a)，農(nóng)業(yè)水資源短缺風(fēng)險較高。

(3)本文基于Pey和ETcy兩個變量來進行關(guān)中平原農(nóng)業(yè)水資源短缺風(fēng)險分析，因為歷史數(shù)據(jù)資料局限性，且隨機變量Pey和ETcy本身可能具有自相關(guān)性，本文所建立的聯(lián)合分布模型初步表征了Pey和ETcy在普遍意義下的統(tǒng)計特征，但二者同期空間相依關(guān)系仍需進一步探究。對隨機變量Pey和ETcy的自相關(guān)性進行研究，并基于此考慮單變量的短期時序相依關(guān)系，將是下一階段的主要研究內(nèi)容。