王現(xiàn)領，劉嘉男

(1.天津市水利科學研究院，天津 300061；2.山東省水利勘測設計院，山東 濟南 250013)

通過現(xiàn)場勘察，遴選出封閉的農(nóng)田排水小區(qū)并進行了水文數(shù)據(jù)量測，采集近兩年的降水及對應降水時段的排水量、土壤含水量、地下水位等實測數(shù)據(jù)，實測數(shù)據(jù)以1 h為計算時段進行整理，排水小區(qū)水文數(shù)據(jù)匯總見表1。

表1 排水小區(qū)水文數(shù)據(jù)匯總表Tab.1 Summary of hydrological data in drainage area

結(jié)合北方平原區(qū)農(nóng)田產(chǎn)流特性，充分考慮在北方平原區(qū)可能同時發(fā)生超滲和蓄滿產(chǎn)流現(xiàn)象，建立超滲-蓄滿垂向耦合產(chǎn)流模型，采用有更好魯棒性的Sobol'法進行敏感性分析，通過模型參數(shù)的一階敏感度和總敏感度分析，確定產(chǎn)流模型的敏感性參數(shù)。模型參數(shù)率定采用參數(shù)優(yōu)選中最有效的SCE-UA算法，通過參數(shù)自動率定與人工率定結(jié)合法，得到優(yōu)選參數(shù)。選取徑流深相對誤差、地下水位相對變化量與土壤含水量變化量的相對誤差和Nash效率系數(shù)為指標，評價模型的適用性。并通過與超滲產(chǎn)流模型、蓄滿產(chǎn)流模型進行對比檢驗，進一步驗證超滲-蓄滿垂向耦合產(chǎn)流模型的合理性。

1 產(chǎn)流模型研究現(xiàn)狀

目前，國內(nèi)外有關產(chǎn)流的研究已涌現(xiàn)出很多水文模型，主要分為集總式水文模型和分布式水文模型，集總式水文模型是以流域整體為研究對象，而分布式水文模型則是將流域劃分為很多較小單元，并假定各單元內(nèi)部屬性的相對一致性，從而進行產(chǎn)匯流模擬。我國水文研究者分析改進降水-徑流相關圖，創(chuàng)建了以新安江模型為典型的蓄滿產(chǎn)流模式和以Horton模型為代表的超滲產(chǎn)流模式[1]，該產(chǎn)流理論對于流域的面積大小無特殊要求，所需數(shù)據(jù)較易獲得，可用于小尺度平原地區(qū)產(chǎn)流計算。但對于多數(shù)流域，并不單單只有一種產(chǎn)流模式，通常兩種產(chǎn)流方式交叉存在，混合式產(chǎn)流模型因其模型機制更加完善，更符合區(qū)域產(chǎn)流實況而逐漸得到推廣，并成為今后產(chǎn)流模型研究的重點。

2 產(chǎn)流計算

采用超滲-蓄滿垂向耦合產(chǎn)流模型，即空間分布的下滲曲線與蓄水容量分布曲線的耦合，首先判斷是否發(fā)生超滲產(chǎn)流，當發(fā)生超滲產(chǎn)流時，采用具有流域分布特征的格林-安普特下滲曲線計算超滲產(chǎn)流量和下滲量[2]，超滲產(chǎn)流的下滲量即為蓄滿產(chǎn)流中的降水輸入量。在土壤缺水量大的面積上，下滲的雨量首先補充土壤含水量不產(chǎn)流，而在缺水量小的面積上，補足土壤缺水量后，發(fā)生蓄滿產(chǎn)流。

超滲產(chǎn)流量計算過程為：

PE=P-E-I

(1)

(2)

FMM=FM×(1+B1)

(3)

(4)

R1=PE-FA

(5)

式中：P為降水量，mm；E為蒸發(fā)量，mm；I為植物截留量，mm；PE為凈雨量，mm；FM為農(nóng)田小區(qū)平均下滲能力，mm/h；FC1為穩(wěn)定下滲率，mm/h；K2為土壤缺水量對下滲率影響的靈敏度系數(shù)；FMM為平均下滲率為FM時上層最大點的下滲率，mm/h；FA為實際下滲量，mm；WM為土壤平均蓄水容量，mm；W為土壤含水量，mm；B1為下滲能力分配曲線的指數(shù)，反映了下滲能力分布的不均勻性；R1為超滲產(chǎn)流量，mm。

地面以下的徑流，采用蓄滿產(chǎn)流，計算公式為：

(6)

WMM=WM×(1+B2)

(7)

(8)

式中：A1為前期影響雨量，mm；WMM為土壤內(nèi)最大點蓄水容量，mm；B2為蓄水容量～面積分配曲線的指數(shù)，反映了蓄水容量分布的不均勻性；R2為蓄滿產(chǎn)流量，mm；式中其他符號含義同前。

總產(chǎn)流量R為超滲產(chǎn)流量R1與蓄滿產(chǎn)流量R2之和，即：

R=R1+R2

(9)

蓄滿產(chǎn)流計算出的產(chǎn)流量R2包括地面徑流R2S和地下徑流R2G兩部分，可以用穩(wěn)定下滲率進一步劃分。計算過程如下：

(10)

式中：FC2為下層土壤穩(wěn)定下滲能力，mm/h。

蒸發(fā)量與土壤含水量的計算采用三層蒸發(fā)模式，將植物的截留作用合并于降水的蒸發(fā)損失中，構(gòu)成了完整的超滲-蓄滿垂向耦合產(chǎn)流模型。

3 模型參數(shù)敏感性分析

3.1 分析方法

敏感性分析方法采用Sobol'法，與傳統(tǒng)的敏感性分析方法相比具有更好的魯棒性，能更好地解決非線性問題[3]。模型表示為：

y=f(X)=f(X1,…,Xm) 0≤X1≤1,i=1,2,…,m

(11)

式中：y為模型輸出的目標函數(shù)；X=X(Xi,…,Xp)為模型的參數(shù)。那么y的方差D(y)可以作如下分解：

(12)

式中：Di為參數(shù)i產(chǎn)生的方差；Dij為兩個參數(shù)i和j相互作用產(chǎn)生的方差；Dijk為3個參數(shù)i、j、k相互作用產(chǎn)生的方差；D1,2,…,m為m個參數(shù)共同作用產(chǎn)生的方差。將上式歸一化后得到各參數(shù)和參數(shù)相互作用的敏感性：

(13)

產(chǎn)流模型參數(shù)i對目標函數(shù)的敏感度值可表示為：

一階敏感度：

(14)

二階敏感度：

(15)

總敏感度：

(16)

式中：Si為參數(shù)i作用的敏感度；二階敏感度Sij為兩個參數(shù)i和j相互作用的敏感度；總敏感度STi為參數(shù)i主要作用和相互作用的敏感度；D～j為除了參數(shù) 之外的參數(shù)的方差。

計算參數(shù)的各階敏感度實際上是計算參數(shù)的各階方差，由于水文模型具有很大的非線性和復雜性，方差的計算參照前人的建議[4]選用拉丁超立方體抽樣對產(chǎn)流模型中待分析的變量進行兩次獨立抽樣，得到ON×m和PN×m兩個矩陣，其中 為抽樣次數(shù)，m為參數(shù)個數(shù)，通過下列公式進行計算，即：

(17)

(18)

(19)

(20)

(21)

(22)

結(jié)合式(11)、式(12)、式(13)就可以得到產(chǎn)流模型參數(shù)各階敏感度系數(shù)值。

因為沒有關于參數(shù)的先驗信息，進行敏感性分析輸入的參數(shù)值是從均勻分布中抽取的，不同的參數(shù)范圍在0～1之間進行線性變換[5]。產(chǎn)流模型中參數(shù)的取值上下限見表2。利用拉丁超立方體抽樣，在參數(shù)的取值范圍內(nèi)隨機生成包含2 000組參數(shù)的兩個樣本。

表2 模型參數(shù)意義及取值范圍Tab.2 Model parameter meaning and range of value

3.2 分析結(jié)果

根據(jù)趙人俊[6]的研究，結(jié)合垂向耦合模型結(jié)構(gòu)，待分析參數(shù)可以劃分為4個層次，見表3。

表3 模型參數(shù)分層Tab.3 Model parameter stratification

其中，在該模型中IMP、α、β均為定值，IMP為不透水面積比例，該農(nóng)田小區(qū)為0.41，不透水面上的產(chǎn)流系數(shù)α取0.95，透水面積流入溝渠的徑流分配系數(shù)β取0.3。

五場降水中，通過對模型中十個參數(shù)對徑流深目標函數(shù)的一階敏感度及總敏感度進行分析，確定該產(chǎn)流模型的敏感性參數(shù)為：WUM、WLM、WM(WDM)、FC1、K2、FC2，在進行模型參數(shù)率定過程中，綜合考慮這些參數(shù)對模型結(jié)果的影響。

4 模型參數(shù)率定

4.1 自動率定方法

自動率定方法采用SCE-UA(ShuffledComplexEvolutionAlgorithm)算法[7]，該方法綜合了隨機搜索算法、聚類分析法及生物競爭演化等方法的優(yōu)點，收斂速度較快，魯棒性強，穩(wěn)定性好，為高維非線性參數(shù)的全局優(yōu)化問題提供了高效的解決辦法，且無需構(gòu)建顯式的目標函數(shù)或者對目標函數(shù)進行求導計算，被認為是流域水文模型參數(shù)優(yōu)選中最有效的方法。

SCE-UA算法詳細步驟如下所示：

(1)算法參數(shù)初始化。若待優(yōu)化的問題是n維，參與進化的復合形的個數(shù)p≥2，則每個復合形所包含的頂點數(shù)目m≥n+1，計算樣本點的數(shù)目s=pm。

(2)隨機產(chǎn)生樣本點。由于缺少先驗信息，采用均勻抽樣在可行空間隨機抽取s個樣本點，分別記為x1，x2，…，xs，并計算每一個點si的函數(shù)值，記為fi=f(xi)，其中i=1,2,…,s。

(3)樣本點排序。把得到的s個樣本點按目標函數(shù)值升序排列，以數(shù)組形式(xi,fi)記錄，其中i=1,2,…,s；f1≤f2≤…≤fs，并存入數(shù)組空間：D={(xi,fi),i=1,2,…,s}，i=1代表了目標函數(shù)值最小的數(shù)組。

(5)復合形演化。依據(jù)競爭的復合形演化算法分別對各個復合形進行演化。

(6)復合形混合。把進化后的每個復合形A1,A2,…,Ap的所有頂點組合成新的點集，存入數(shù)組空間D中，并對D按目標函數(shù)的升序進行排序。

(7)收斂性判斷。若達到收斂條件則停止，否則返回第(4)步。

在進行參數(shù)優(yōu)化時，將率定期五場降水綜合考慮，以徑流深相對誤差為核心，采用表3所示目標函數(shù)，通過SCE-UA算法，自動尋找到使目標函數(shù)值達到最小，且滿足單場次降水徑流深相對誤差在20%范圍內(nèi)的最優(yōu)參數(shù)組。

4.2 自動率定結(jié)果

SCE-UA算法自身包含多個參數(shù)，算法參數(shù)取值[8,9]見表4。參數(shù)優(yōu)化的可行域與表4-1所示相同。

通過拉丁超立方體抽樣方法抽取42組參數(shù)，對率定期五場降水綜合考慮，最大迭代次數(shù)T分別設置為100、200、500、1 000、2 000、5 000、10 000，并分別連續(xù)運行五次，并初步率定出一組較優(yōu)參數(shù)，結(jié)果見表5。此時，總目標函數(shù)值以及單場次徑流深相對誤差值見表6，此時SCE-UA算法總目標函數(shù)值最小，場次徑流深相對誤差均在20%以內(nèi)，滿足要求。

表4 SCE-UA算法參數(shù)設置Tab.4 SCE-UA algorithm parameter setting

表5 模型參數(shù)初步率定結(jié)果Tab.5 Initial rate determination of model parameters

表6 初步率定目標函數(shù)值及單場次徑流深相對誤差Tab.6 Preliminarily rate determination objective function value and relative error of single runoff depth

4.3 參數(shù)人工率定與自動率定相結(jié)合

模型參數(shù)優(yōu)化的結(jié)果不僅依賴于參數(shù)的取值范圍，參數(shù)的敏感性對其也有較大影響，敏感參數(shù)通過自動尋優(yōu)較容易得到最優(yōu)值，但是不敏感參數(shù)容易在“不敏感區(qū)間”中陷入次優(yōu)，導致結(jié)果的不確定性。

本文在SCE-UA參數(shù)優(yōu)選過程中，發(fā)現(xiàn)蒸發(fā)計算與兩水源劃分參數(shù)，即K1、C、WUM、WLM、FC2波動比較大；產(chǎn)流計算參數(shù)FC1、WDM相對比較收斂；產(chǎn)流計算參數(shù)K2、B1、B2較穩(wěn)定。其中K1、C、B1、B2為不敏感參數(shù)，在率定過程中，可以對K1取多次尋優(yōu)結(jié)果的平均值0.75；C根據(jù)經(jīng)驗可以取0.17；B1、B2取自動率定結(jié)果值，分別為0.3和0.2。

結(jié)合趙人俊提出的客觀優(yōu)化理論[10]，不僅以徑流深相對誤差為評價指標，還需考慮地下水位相對變化與土壤含水量變化過程線，提高模型的模擬精度。地下水位主要受地下徑流的影響，著重調(diào)節(jié)FC1與FC2使地下水位變化過程線滿足要求；K2、WUM、WLM、WDM與土壤含水量有密切關系，調(diào)節(jié)這四個參數(shù)使土壤含水量評價指標滿足要求。自動率定與人工率定結(jié)合得到的參數(shù)數(shù)值見表7。

表7 模型參數(shù)率定結(jié)果Tab.7 Model parameter rate determination result

5 模型適用性評價

5.1 評價指標

產(chǎn)流計算的準確度與可信度是進行排澇模數(shù)計算的保障，選取徑流深相對誤差作為水量平衡評價指標，選取地下水位相對變化量與土壤含水量變化量的相對誤差Re和Nash效率系數(shù)Ens為指標，共同評價模型模擬精度。一般認為，誤差控制在±20%以內(nèi)，Ens>0.6時模擬結(jié)果較好，模擬精度越高，誤差越小。計算公式為：

(23)

(24)

式中：Re表示水文要素的相對誤差；Sobj表示總徑流深(地下水位總相對變化量/土壤含水量總變化量)的實測值，mm；Ssim表示總徑流深(地下水位總相對變化量/土壤含水量總變化量)的模擬值，mm；Ens表示Nash效率系數(shù)；Si,obj為實測地下水位相對變化量(土壤含水量變化量)序列，mm；Si,sim為模擬地下水位相對變化量(土壤含水量變化量)序列，mm；n為序列長度。

5.2 評價結(jié)果

研究中對地下水位相對變化量和土壤含水量使用Nash系數(shù)、相對誤差來評價，水量平衡校核即模擬徑流深使用相對誤差來評價。模型評價指標分析見表8。

由表8可知，農(nóng)田小區(qū)的率定期與驗證期地下水相對變化量與土壤含水量的 效率系數(shù)均大于0.6，各項相對誤差均控制在20%之內(nèi)，模型適用性較好。

表8 模型評價指標分析Tab.8 Model evaluation index analysis

6 模型檢驗

6.1 模型檢驗方法

為了進一步驗證構(gòu)建的產(chǎn)流模型，分別采用單一的超滲產(chǎn)流模型、蓄滿產(chǎn)流模型以及已構(gòu)建的超滲-蓄滿垂向耦合模型進行計算，計算模型相對有效性。計算公式如下：

(25)

(26)

式中：CEFc代表超滲-蓄滿垂向耦合產(chǎn)流相對于超滲產(chǎn)流的有效性；CEFx代表超滲-蓄滿垂向耦合產(chǎn)流相對于蓄滿產(chǎn)流的有效性；R0為實測降水徑流深，mm；Rc、Rx、RH分別為用超滲、蓄滿和超滲-蓄滿垂向耦合模型計算的降水徑流深，mm。

6.2 模型檢驗結(jié)果

按照公式(25)、(26)選取八場降水計算垂向耦合模型相較于單一超滲模型和單一蓄滿模型的有效性，結(jié)果見表9。

由表8可知，超滲-蓄滿垂向耦合模型相對于超滲產(chǎn)流的有效性為0.869，相對于蓄滿產(chǎn)流的有效性為0.912，進一步證明了在農(nóng)田小區(qū)采用超滲-蓄滿垂向耦合產(chǎn)流模型是合理的。

表9 三種產(chǎn)流模型下產(chǎn)流結(jié)果計算成果表Tab.9 Calculation results of runoff yield under three kinds of runoff models

7 結(jié) 論

通過農(nóng)田產(chǎn)流模型構(gòu)建及檢驗，證明超滲-蓄滿垂向耦合產(chǎn)流模型應用于華北平原地區(qū)是合理的，主要結(jié)論如下：

(1)產(chǎn)流模型的敏感性參數(shù)為：上層土壤蓄水容量、下層土壤蓄水容量、深層土壤蓄水容量、上層穩(wěn)定下滲率、土壤缺水量靈敏度系數(shù)、下層土壤穩(wěn)定下滲率。

(2)通過自動參數(shù)率定與人工率定結(jié)合法，得到優(yōu)選參數(shù)，模擬結(jié)果的有效性系數(shù)均在0.6以上，相對誤差均在±20%之內(nèi)，說明模型精度較高，模型具有較好的適用性。

(3)通過與單一模型對比檢驗，超滲-蓄滿垂向耦合模型相對于超滲產(chǎn)流模型的有效性為0.869，相對于蓄滿產(chǎn)流模型的有效性為0.912，說明超滲-蓄滿垂向耦合模型相對單一模型具有更高的模擬精度，也進一步證明了超滲-蓄滿垂向耦合產(chǎn)流模型的合理性。

□