基于Sobol法的寧夏固海揚水灌區(qū)ET0敏感性分析

葉 威，伍昊洋，鄒全程，趙 璐,2，崔寧博，梁 川

(1.四川大學(xué)水利水電學(xué)院，成都 610065 ；2 農(nóng)業(yè)部旱作節(jié)水農(nóng)業(yè)重點開放實驗室，北京 100081)

0 引 言

潛在蒸散發(fā)(ET0)是水文循環(huán)研究的重要參數(shù)，ET0、降水和徑流一起決定一個地區(qū)的水分平衡，在研究區(qū)域水循環(huán)中占有重要的地位。ET0也是農(nóng)業(yè)水資源評價、水量平衡模型運行和實際蒸散發(fā)估算中必不可少的重要參量，它獨立于作物類型、生長階段和作物管理方式，主要受氣象因子影響，因此，計算ET0對單個氣象因子或多個氣象因子變化的敏感系數(shù)[1-3]，可以定量的反應(yīng)主要氣象因子變化和ET0之間的關(guān)系。根據(jù)ET0對主要氣象因子的敏感性系數(shù)，可以快速估算當某個氣象因子(一階敏感性系數(shù))或多個氣象因子(總敏感性系數(shù))發(fā)生變化時ET0的相應(yīng)變化量。量化各種氣象因子的變化對ET0的影響，對于各個地區(qū)灌溉制度的制定和水循環(huán)的研究具有重要意義。

國內(nèi)外許多學(xué)者常采用敏感曲線法進行ET0的敏感性分析，即把因變量變化與自變量變化的比值繪成曲線來描述敏感系數(shù)的特征。侯蘭功[4]等人利用此方法對額濟納綠洲進行分析，得出太陽輻射為ET0最為敏感的氣象因子。也有學(xué)者采用敏感系數(shù)法，曹雯[5]等人利用此法對西北地區(qū)ET0的各因子敏感性進行分析，發(fā)現(xiàn)了太陽輻射是ET0最為敏感的氣象因子，這一結(jié)果與其他學(xué)者對我國長江流域、松嫩平原西部和東北地區(qū)進行研究所得的結(jié)果不同[6-8]，由此得出不同地區(qū)ET0的變化具有不同的控制因子。Sobol全局敏感性方法是進行敏感性分析的一種常用方法[9,10]，可以分析各氣象因子的獨立敏感性和總敏感性，這是其他敏感性分析方法無法實現(xiàn)的，在各個領(lǐng)域得到了廣泛的應(yīng)用，然而，此法卻極少應(yīng)用到ET0的敏感性分析，僅有泮蘇莉[11]等人利用Sobol全局敏感性系數(shù)分析方法對嵊泗列島進行分析，得到了海島地區(qū)ET0對各個氣象因子的敏感性系數(shù)。此外，關(guān)于寧夏固海揚水灌區(qū)ET0對各氣象因子的敏感性分析也未見報道?；谏?，本文擬采用Sobol[11-14]全局敏感性方法分析寧夏固海揚水灌區(qū)ET0對主要氣象因子的敏感性，得到ET0對每個氣象因子的一階敏感性系數(shù)和總敏感性系數(shù)。

1 研究區(qū)域概況

自20世紀70年代以來，為逐步解決寧夏回族自治區(qū)中部干旱帶日趨嚴重的干旱缺水狀況，興建了紅寺堡、固海、固海擴灌和鹽環(huán)定等揚水工程。該區(qū)域?qū)儆谥袦貛Ц珊祬^(qū)，太陽輻射強，大陸性氣候明顯，干旱少雨且強烈蒸發(fā)。

寧夏固海揚水灌區(qū)(東經(jīng)105°26′24″～106°14′24″，北緯36°26′24″～37°33′30″)位于清水河河谷平原下游，為研究固海揚水灌區(qū)ET0對各氣象因子的敏感性，本文擬采用同心站1961-2016年逐日的氣象數(shù)據(jù)(來自中國氣象數(shù)據(jù)共享網(wǎng)https:∥data.cma.cn/)，利用Penman-Monteith方法計算ET0，采用Sobol全局敏感性分析法定量分析ET0對各氣象因子的敏感性，對比興建揚水工程前后各氣象因子敏感系數(shù)的變化。

2 研究方法

2.1 FAO Penman-Monteith方法[15-17]

(1)

式中：ET0為潛在蒸散發(fā)，mm/d；Rn為凈輻射，MJ/(m2·d)；G為土壤熱通量, MJ/(m2·d)；G相對于Rn取值很小，在計算步長大于等于1 d的情況下，可以忽略為0；T為2 m處日平均氣溫，℃；U2為2 m處的風速，m/s；es為飽和水汽壓，kPa；ea為實際水汽壓，kPa；es-ea為飽和水氣壓差，kPa；Δ為水汽壓曲線斜率，kPa/℃；γ為濕度計算常數(shù)，kPa/℃。

2.2 Mann-Kendall趨勢檢驗[18-20]

Mann-Kendall方法是基于秩序的非參數(shù)方法，被國內(nèi)外學(xué)者廣泛應(yīng)用于分析降水、氣溫和ET0等因子時間序列的趨勢變化。本文選取Mann-Kendall方法對寧夏揚水灌區(qū)1961-2016年各氣象因子和ET0時間序列進行趨勢分析，并檢驗ET0在這56年間的突變。

2.3 Sobol全局敏感分析法

Sobol方法是基于方差分配的全局敏感性分析方法，原理大致如下：模型函數(shù)可以分解為單個變量與單個變量之間相互結(jié)合的函數(shù)，公式(2)中的方差由近似蒙特卡洛數(shù)值積分進行估算。

(2)

式中：D為模型的總方差；Di為變量Xi通過變量Xi作用所貢獻的方差；Di,j為變量Xi通過變量Xi，Xj所貢獻的方差；Di,j,…,n為變量Xi通過變量Xi,j,…,n所貢獻的方差。

定義變量及變量相互作用的方差與總方差的比值為敏感性指數(shù)，反映變量Xi對模型輸出總方差的直接貢獻率,即變量Xi一階敏感性指數(shù)Si可表示如下：

(3)

同理，變量Xi的二階敏感性系數(shù)可表示為：

(4)

變量Xi的總敏感性指數(shù)即為各階敏感性指數(shù)之和，表示如下：

STi=Si+Si,j+…+Si,j,…,n

(5)

Si體現(xiàn)的是某一輸入?yún)?shù)在輸出方差上的影響，而STi體現(xiàn)的是輸入?yún)?shù)i的方差的總貢獻，包括它的一階和更高階方差。公式(3)用于計算一階敏感性系數(shù)，代表ET0對各氣象因子的敏感性，公式(5)用于計算總敏感性系數(shù)，代表ET0對各氣象因子的總敏感性。

3 結(jié)果分析

3.1 Mann-Kendall檢驗結(jié)果

根據(jù)Penman-Monteith公式計算ET0并利用Mann-Kendall趨勢檢驗法檢測其在1961-2016年的突變，結(jié)果如圖1中所示。

圖1 ET0突變檢驗圖Fig.1 ET0 mutation test chart

其中，實線為順序變化曲線UFk、虛線表示逆序變化曲線UBk，2條臨界線±1.96(顯著性水平為0.05)。由圖1可知，兩條序列曲線的交點，大致位于1980年，因此，1980年為ET0在1961-2016年這56年間的突變年。 寧夏固海揚水灌區(qū)是20世紀70年代以來發(fā)展的大型灌區(qū)，通過該灌區(qū)的修建，引入大量的黃河水，局部的水循環(huán)有一定的變化，對ET0有一定的影響，故1980年作為ET0的突變年具有一定的合理性。

根據(jù)突變點檢驗的結(jié)果，本文將平均氣溫、最高氣溫、最低氣溫、相對濕度、日照時數(shù)、風速以及ET0分為1961-1980年、1981-2016年和1961-2016年三個時間段分別進行時間趨勢分析，分析結(jié)果見表1。

表1 氣象因子趨勢分析結(jié)果Tab.1 Meteorological factor trend analysis results

注：Z>0時各因素變化趨勢上升，Z<0時為下降， *表示顯著，**表示極顯著，下同。

由表1可知：寧夏固海揚水灌區(qū)ET0在1961-1980年呈下降趨勢，1981-2016年呈顯著的上升趨勢， 1961-2016年呈上升趨勢；平均氣溫、最高氣溫和最低氣溫三個時間段均呈上升趨勢，1961-2016年呈極顯著的上升趨勢；相對濕度在1961-1980年呈下降趨勢，1981-2016年呈顯著上升趨勢，1961-2016年呈上升趨勢；風速在1961-1980年呈極顯著的下降趨勢，1981-2016年呈極顯著的上升趨勢，1961-2016年呈上升趨勢；僅日照時數(shù)在1961-2016年呈顯著的下降趨勢。綜上所述，1981-2016年各氣象因子和ET0較1961-1980年都有上升趨勢；1961-2016年除日照時數(shù)呈下降趨勢，其余氣象因子和ET0均呈上升趨勢。

3.2 一階和總敏感系數(shù)年內(nèi)變化結(jié)果

利用Sobol全局敏感性分析方法分別計算1961-1980年和1981-2016年兩個時間段年內(nèi)的一階和總敏感性系數(shù)，圖2為1961-1980年一階和總敏感性系數(shù)年內(nèi)變化，圖3為1981-2016年一階和總敏感性系數(shù)年內(nèi)變化。

圖3 1981-2016年各氣象因子一階和總敏感性系數(shù)年內(nèi)變化Fig.3 Changes in the first-order and total sensitivity coefficients of various meteorological factors from 1981 to 2016

由圖2、圖3可知，各氣象因子的總敏感性系數(shù)大體上小于一階敏感性系數(shù)，其中，也有個別月份出現(xiàn)總敏感性系數(shù)大于一階敏感性系數(shù)的情況，如：1961-1980年風速的總敏感性系數(shù)在1-3月和9-12月均大于一階敏感性系數(shù)。1961-1980年和1981-2016年兩個時間段各氣象因子的總敏感性系數(shù)的年內(nèi)變化大致呈相同的變化趨勢。平均氣溫、最低氣溫和相對濕度的總敏感性系數(shù)年內(nèi)近似呈一條水平線，前兩者敏感性系數(shù)較小，約為0.02，相對濕度的敏感性系數(shù)較大，約為0.2。最高氣溫和風速在年內(nèi)的變化呈波谷型，年內(nèi)敏感性系數(shù)的最小值一般出現(xiàn)在8月左右，而最大值一般出現(xiàn)在12月左右。1961-1980年和1981-2016年兩個時間段，最高氣溫的總敏感性系數(shù)最大值分別為0.26和0.28，最小值分別為0.04和0.039，風速的總敏感性系數(shù)最大值分別為0.43和0.40，最小值分別為0.093和0.092。1961-1980年和1981-2016年日照時數(shù)總敏感性系數(shù)在年內(nèi)的變化呈波峰型，其中最大值一般出現(xiàn)在7月左右，分別為0.20和0.25，最小值一般出現(xiàn)在12月左右，近似為0?？偟恼f來，總敏感性系數(shù)的大小排序為風速>相對濕度>最高氣溫>日照時數(shù)>最低氣溫>平均氣溫。

各氣象因子的一階敏感性系數(shù)在年內(nèi)呈多樣性和不確定性變化的特點，平均氣溫和日照時數(shù)的一階敏感性系數(shù)在年內(nèi)的變化大致為波峰型，最小值均大致出現(xiàn)在1月左右，平均氣溫的最大值出現(xiàn)在4月左右，日照時數(shù)出現(xiàn)在7月左右。兩個時間段，平均氣溫的一階敏感性系數(shù)最大值分別為0.26和0.21，最小值分別為0.01和0.08；日照時數(shù)的一階敏感性系數(shù)最大值分別為0.55和0.63，最小值分別為0.006和0.008。最低氣溫、最高氣溫、相對濕度和風速的一階敏感性系數(shù)年內(nèi)呈峰谷交替變化的特點，以1961-1980年最高氣溫的一階敏感性系數(shù)為例，其極大值分別在4月和11月出現(xiàn)，分別為0.32和0.37，極小值分別在3月和6月出現(xiàn)，分別為0.15和0.004?？偟恼f來，一階敏感性系數(shù)的大小排序為相對濕度>日照時數(shù)>最高氣溫>風速>平均氣溫>最低氣溫。

3.3 1961-2016年敏感性年際變化

由圖2和圖3可知，1961-1980年和1981-2016年兩個時間段，總敏感性系數(shù)年內(nèi)呈現(xiàn)基本相同的變化特點，而一階敏感性系數(shù)在兩個時間段呈不同的變化特點，其大小有明顯差異，為探究這種差別，對各氣象因子的一階敏感性系數(shù)在1961-2016年的變化進行了趨勢分析，計算結(jié)果如表2和表3。由表2可知，一階敏感性系數(shù)相對濕度最大，最低氣溫最小，三個時間段，最高氣溫和相對濕度的變化量較大，日照時數(shù)和風速次之，平均氣溫和最低氣溫較小。表3根據(jù)Mann-Kendall趨勢檢驗法的結(jié)果，更準確的反映了各氣象因子在三個時間段的變化趨勢。1961-1980年平均氣溫和相對濕度的一階敏感性系數(shù)呈上升趨勢，最高氣溫、最低氣溫、日照時數(shù)和風速的一階敏感性系數(shù)呈下降趨勢。1981-2016年平均氣溫基本無變化，風速的一階敏感性系數(shù)呈上升趨勢，最高氣溫、最低氣溫、相對濕度和日照時數(shù)呈下降趨勢。1961-2016年，平均氣溫、最高氣溫、最低氣溫和相對濕度的一階敏感性系數(shù)呈現(xiàn)下降的趨勢，日照時數(shù)的一階敏感性系數(shù)呈顯著的下降趨勢，風速的一階敏感性系數(shù)呈上升趨勢。

表2 一階敏感性系數(shù)平均值Tab.2 Average value of first-order sensitivity coefficient

表3 一階敏感性系數(shù)趨勢分析結(jié)果Tab.3 First-order sensitivity coefficient trend analysis results

4 結(jié) 語

根據(jù)寧夏固海揚水灌區(qū)1961-2016年同心站的逐日氣象數(shù)據(jù)，計算該地區(qū)的ET0，利用Sobol全局敏感性分析方法計算ET0對各氣象因子的一階和總敏感性系數(shù)，并利用Mann- Kendall趨勢檢驗探究ET0和各氣象因子的敏感性趨勢變化，得出主要結(jié)論如下。

(1)各氣象因子的總敏感性系數(shù)在年內(nèi)基本上小于一階敏感性系數(shù)，個別月份例外。各因子的總敏感性系數(shù)年內(nèi)呈直線型(平均氣溫、最低氣溫和相對濕度)、波谷型(最高氣溫和風速)和波峰型(日照時數(shù))三種變化特點，而一階敏感性系數(shù)呈不確定性和多樣性變化的特點。

(2)1981-2016年各氣象因子和ET0較1961-2016年都有上升趨勢，1961-2016年除日照時數(shù)外的氣象因子和ET0均呈上升趨勢，日照時數(shù)呈下降趨勢。1961-2016年，平均氣溫、最高氣溫、最低氣溫和相對濕度的一階敏感性系數(shù)呈現(xiàn)下降的趨勢，日照時數(shù)的一階敏感性系數(shù)呈顯著的下降趨勢，風速的一階敏感性系數(shù)呈上升趨勢。

(3)本文利用Mann-Kendall檢驗出1980年為ET0突變年，并將1961-1980年和1981-2016年進行對比，探究寧夏固海揚水灌區(qū)的修建對ET0、各氣象因子以及ET0對各氣象因子的敏感性的變化趨勢的影響。根據(jù)計算結(jié)果可以看出， 1981-2016年各氣象因子和ET0均有上升的趨勢，ET0對各氣象因子的一階敏感性系數(shù)中，風速的一階敏感性系數(shù)呈上升趨勢，而平均氣溫、最高氣溫、最低氣溫、相對濕度和日照時數(shù)的一階敏感性系數(shù)呈下降趨勢。這些變化特點表明大型灌區(qū)的修建在一定程度上可以影響ET0、各氣象因子以及ET0對各氣象因子的敏感性，ET0呈上升趨勢說明人工引黃河水對灌區(qū)的水分平衡起到一定的調(diào)節(jié)作用。此外，本文僅基于有限的氣象數(shù)據(jù)分析了ET0、各氣象因子和ET0對各氣象因子的敏感性，具有一定局限性。

□