盧澤華,王君勤

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

隨著全球氣候變化和社會經(jīng)濟(jì)發(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)進(jìn)行有效的風(fēng)險(xiǎn)管理已成為水資源科學(xué)發(fā)展的必然趨勢[3]。農(nóng)業(yè)水資源短缺風(fēng)險(xiǎn)評估作為農(nóng)業(yè)水資源風(fēng)險(xiǎn)管理的基礎(chǔ)，可為區(qū)域水資源優(yōu)化配置提供基礎(chǔ)依據(jù)，對區(qū)域種植結(jié)構(gòu)調(diào)整，干旱風(fēng)險(xiǎn)應(yīng)對具有重要意義[4]。

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

基于變量間非線性相關(guān)關(guān)系構(gòu)建的Copula函數(shù)能獨(dú)立于隨機(jī)變量的邊緣分布反映變量間的相關(guān)性，可構(gòu)造任意邊緣分布的聯(lián)合分布函數(shù)，且在轉(zhuǎn)換過程中不會改變原始隨機(jī)變量，在干旱風(fēng)險(xiǎn)分析、洪水遭遇風(fēng)險(xiǎn)分析方面已得到廣泛應(yīng)用[10-11]，也有學(xué)者利用Copula函數(shù)進(jìn)行水資源短缺風(fēng)險(xiǎn)分析研究[12-13]，如丁志宏等[12]構(gòu)建了基于Copula函數(shù)的寧夏衛(wèi)寧灌區(qū)降水量(P)和參考作物騰發(fā)量(ET0)年際聯(lián)合分布模型；Zhang等[13]構(gòu)建了基于Copula函數(shù)的河南省陸渾灌區(qū)自然降水條件下的水資源短缺風(fēng)險(xiǎn)模型。然而有效降水量(Pe)和作物需水量(ETc)是描述灌溉系統(tǒng)的基本變量[14]，二者的演變特征與匹配程度直接影響農(nóng)業(yè)水資源管理策略，比P和ET0更直觀反映自然降水條件下農(nóng)業(yè)水資源的供需平衡。關(guān)中平原地處內(nèi)陸、遠(yuǎn)離海洋，其年降水量僅500～750 mm，對水資源短缺的響應(yīng)尤為敏感，然而目前對關(guān)中平原Pe和ETc相關(guān)關(guān)系的研究較少，且大多是闡釋其年內(nèi)、年際變化和簡單的相關(guān)關(guān)系[15-16]，對其非線性耦合關(guān)系還未報(bào)道。本文主要在現(xiàn)有研究基礎(chǔ)上，首次運(yùn)用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)斋@后即時(shí)播種，當(dāng)年9月底收獲[21]。本文以冬小麥播種—夏玉米收獲(10月—翌年9月底)作為一個單位年，計(jì)算單位年有效降水量(Pey)和作物需水量(ETcy)。

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

2 研究方法

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

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

(1)

(2)

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

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

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

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

采用Copula函數(shù)進(jìn)行Pey和ETcy2個特征變量的聯(lián)合之前，首先要確定其邊緣分布函數(shù)，同時(shí)要考慮變量間的相關(guān)性。設(shè)X和Y分別代表Pey和ETcy，其邊緣分布函數(shù)分別為FX(x)和FY(y)，選用常見的5種單變量分布函數(shù)(正態(tài)分布、伽馬分布、對數(shù)正態(tài)分布、廣義極值分布及韋爾布分布)[24-25]分別對其進(jìn)行擬合，利用K-S檢驗(yàn)確立最優(yōu)邊緣分布函數(shù)。同時(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ù)均采用極大似然法進(jìn)行估計(jì)[24]。

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

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

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

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

(4)

(5)

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

3 結(jié)果與分析

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

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

表2 5種分布函數(shù)對關(guān)中平原Pey和ETcy擬合的K-S檢驗(yàn)結(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)驗(yàn)累積概率與理論累積概率擬合較好，均呈極顯著相關(guān)(α<0.01)，進(jìn)一步表明邊緣分布及參數(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)驗(yàn)概率與理論概率的決定系數(shù)

3.1.2Pey和ETcy的二維聯(lián)合概率分布 采用Kendall秩相關(guān)系數(shù)τ分析度量關(guān)中東部、關(guān)中西部Pey和ETcy間的相關(guān)性，根據(jù)τ與θ的關(guān)系，計(jì)算Copula函數(shù)的參數(shù)θ，并對Copula函數(shù)進(jìn)行擬合優(yōu)度評價(jià)，結(jié)果見表5?？梢钥闯鲫P(guān)中東部、關(guān)中西部Pey和ETcy均呈負(fù)相關(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ù)進(jìn)行K-S檢驗(yàn)，取顯著性水平α=0.05，n=55時(shí)，查柯爾莫格洛夫檢驗(yàn)分位數(shù)表得出對應(yīng)分位點(diǎn)D0為0.179 8，表5中檢驗(yàn)統(tǒng)計(jì)量D小于分位點(diǎn)0.179 8，表明擬合優(yōu)度均達(dá)顯著水平(α<0.05)，Pey和ETcy聯(lián)合分布的經(jīng)驗(yàn)累積概率與理論累積概率擬合較好，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)險(xiǎn)模型應(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)度評價(jià)

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

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

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

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

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

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

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

因此關(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ū)的特點(diǎn)基礎(chǔ)之上，推進(jìn)農(nóng)村特色產(chǎn)業(yè)區(qū)域化布局、規(guī)?；l(fā)展、產(chǎn)業(yè)化經(jīng)營，形成產(chǎn)業(yè)聚集效應(yīng)和發(fā)展的比較優(yōu)勢，促進(jìn)高效節(jié)水產(chǎn)業(yè)的快速發(fā)展。

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

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