基于PEST的MODFLOW-HYDRUS耦合模型參數優化

胡 丹，伍靖偉，黃介生

(武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072)

0 引 言

土壤水和地下水運動是水文循環的重要環節，能夠準確模擬土壤水、地下水運動對水資源高效利用具有重要意義[1,2]。地下水與土壤水之間有著直接的水力聯系和水通量交換，為了更好地模擬區域尺度土壤水、地下水的運動，近年來發展了一些飽和-非飽和耦合模型，這些模型的主要區別在于對土壤水運動的模擬，其主要方式有3種：第一種是水均衡模型，如MODFLOW中的REC和ET模塊[3]，以土壤水質量平衡為原理，獨立考慮蒸發、根系吸水等水文過程，由于水均衡模型不涉及水頭，無法與地下水的計算直接聯系，過分簡化了土壤水運動，故會帶來較大誤差，但求解相對簡單；第二種是基于動力波方程，如MODFLOW-UZF1[4]，模擬土壤水運動時只考慮重力作用，而忽略毛管力作用；第三種是基于理查茲方程，如三維非飽和-三維飽和的耦合模型MODFLOW-VSF[5]、一維非飽和-三維飽和的耦合模型MODFLOW-HYDRUS[6]，理查茲方程是由達西定律和質量守恒定律推導而得，對土壤水運動的描述更加準確，但由于具有高度非線性，計算成本也更大。飽和-非飽和耦合模型能否做到模擬精度和計算成本的平衡直接關系到模型的運用推廣。Twarakavi等[7]對MODFLOW-UZF1、MODFLOW-VSF和MODFLOW-HYDRUS三種耦合模型進行了比較，發現進行區域模擬時，MODFLOW-HYDRUS和MODFLOW-VSF的模擬精度比MODFLOW-UZF1更高，在計算要求方面，MODFLOW-HYDRUS和MODFLOW-UZF1的計算量比MODFLOW-VSF更小，尤其對于淺地下水、蒸發強烈的地區，MODFLOW-HYDRUS進行土壤水-地下水區域模擬時，相比MODFLOW-UZF1和MODFLOW-VSF，在模擬精度和計算量方面能做到更好的平衡。

飽和-非飽和耦合模型進行區域水分運動模擬時，往往需要輸入較多的參數，且水文地質參數和土壤水力參數空間變異性較大，導致模擬結果存在較大不確定性。模型參數的獲得目前主要有實測法和參數率定法，其中實驗法測得的數據由于尺度效應的存在，很難直接運用到模型中。參數率定法包括試錯法、參數自動優化方法以及試錯法與參數自動優化法聯合優化。傳統的試錯法依賴使用者對模型的了解及經驗，主觀性太強，花費時間較長，且難以達到參數最優化。參數自動優化方法能夠較快得到優化解，但是對初值的依賴較大，不同的初值可能導致差別較大的率定結果[8]。試錯法與參數自動優化法聯合優化既能發揮模型使用者的知識和經驗，又能利用先進的優化技術[9]。

PEST[10]是獨立于計算模型的參數優化方法，已經分別在地下水、土壤水運動模型中得到廣泛運用[11-13]，但是在耦合模型中運用效果和效率如何尚未有研究。

本文以內蒙古河套灌區永聯試驗區地下水觀測井數據和土壤含水率實測數據為基礎，采用試錯法和PEST相結合的方法，對MODFLOW-HYDRUS模型的水文地質參數和土壤水力參數進行優化，以提升區域尺度上飽和-非飽和帶中水分運動的模擬效果。

1 材料與方法

1.1 研究區域及數據

研究區域為內蒙古河套灌區五原縣永聯村境內的永聯試驗區(見圖1)，東經108°00′～108°01′，北緯41°04′～41°05′，南北長約13 km，東西寬約5 km，總面積約2 875 hm2，其中灌溉面積約1 453 hm2，土地利用類型主要是耕地、裸地(包括荒地、渠道、排水溝等)和村莊。試驗區地勢自南向北稍有傾斜，中間較高，東西兩側略低，整體較為平坦，地面高程1 028.6～1 025.6 m。試驗區邊界條件較為清楚，北邊是六排干排水邊界，南邊是皂火渠，東西兩邊分別是永什分干溝和乃永分干溝(見圖1)。

試驗區內布置有11口地下水觀測井(見圖1)，每5 d觀測1次地下水位，每個地下水觀測井附近布置有剖面土壤含水率監測點，每10 d進行1次取樣測定，取樣深度為0～100 cm，分為0～10，10～30，30～50，50～70，70～100 cm 5個層次。地下水位觀測采用測鐘法，土壤含水率采用烘干法測定。

圖1 研究區示意圖Fig.1 Sketch map of study area

1.2 飽和-非飽和耦合模型

非飽和帶水流運動考慮為一維垂向運動，用一維Richards方程來描述：

(1)

式中：θ為土壤的體積含水量，cm3/cm3；t為時間，d；h為壓力水頭，m；z為空間坐標，向上為正，m；S為匯項，m3/(m3·d)；K(h)為非飽和水力傳導度，m/d。

飽和帶水流運動為三維運動，其控制偏微分方程為：

(2)

式中：h為水頭，m；x，y，z為空間坐標，m；Kxx，Kyy，Kzz為滲透系數在x，y，z方向上的分量，這里假定滲透系數主軸方向與坐標軸方向一致，m/d；W為匯項，m3/(m3·d)；t為時間，d；SS為給水度，無量綱。

非飽和帶與飽和帶的耦合過程參考Seo等[5]的介紹。

本文采用2007年5月1日至11月30日(共214 d)期間實測的地下水位和土壤含水率數據對模型參數進行率定。

(1)網格劃分與土壤剖面設置。將研究區域劃分為25×50×2網格，其中有效網格為742個，水平方向上網格大小為160 m×264 m，垂直方向上，第一、二層底部高程分別為1 020.29 m和973.75 m。根據地面高程、地下水位、土地利用類型等情況，采用k均值聚類分析方法，在研究區內設置25個土壤剖面(見圖2)，其中1～7號土壤剖面為耕地，8～14號為排水溝，15～16號為人民支渠，17號為皂火渠，18～21號為荒地，22～25號為村莊。土壤剖面的高度為地面到1 020.29 m高程處。

圖2 網格與土壤剖面Fig.2 Grids and soil profiles

(2)初始輸入數據。以2007年5月1日實測的地下水位數據插值得到整個區域的地下水位作為初始地下水位[見圖3(b)]；土壤剖面的初始壓力水頭分布按照Seo等[5]的方法計算得到。

圖3 地面高程和初始地下水位Fig. 3 Ground surface elevation and initial groundwater table elevation

(3)邊界條件。上邊界條件為變化的大氣邊界和灌溉補給。從五原縣氣象站獲得該研究區域氣象資料，包括降雨量、氣溫、濕度、風速、日照時數等，采用FAO-56[14]介紹的方法計算每日的潛在蒸發量和蒸騰量，將人民支渠口部的凈引水量作為區域內的灌溉水量，而灌溉滲漏量概化為人民支渠上的線性補給。耕地(1～7號土壤剖面)的上邊界條件包括降雨、蒸發、蒸騰和灌溉[見圖4(a)]，荒地(18～21號土壤剖面)、排水溝(8～14號土壤剖面)和皂火渠(17號土壤剖面)的上邊界條件包括降雨和蒸發[見圖4(b)]，人民支渠(15～16號土壤剖面)的上邊界條件包括降雨、蒸發和灌溉滲漏量[見圖4(c)]，村莊(22～25號土壤剖面)的上邊界條件只有降雨[見圖4(d)]。第二含水層下是約1 m的黏土層，故將下邊界條件視為零通量邊界。將六排干、永什分干和乃永分干三條排水溝設置為排水邊界，皂火渠灌溉渠道設置為河流邊界。

圖4 不同土壤剖面的上邊界條件Fig.4 Upper boundary conditions of different soil profiles

(4)模型參數。研究區域兩個含水層分別是沙壤土和砂土，水文地質參數的初始值根據五原縣永聯村鉆孔水文地質報告取值(見表1)，van Genuchten[15]模型中土壤水力參數的初始值根據Carsel and Parrish[16]取值(見表1)。然后采用試錯法對模型各參數進行調整(見表1)，最后進行參數自動優化。

表1 土壤水力參數和水文地質參數取值Tab.1 Values of soil hydraulic parameters and hydrogeological parameters

1.3 PEST參數優化算法

PEST算法采用高斯牛頓法與梯度下降法結合快速收斂性，又具備梯度下降法的全局搜索性，它能夠使目標函數快速有效地收斂，進行參數優化時模型運行的次數比其他算法更少[17]。

PEST算法的核心是對目標函數求解最小值，目標函數是實驗觀測值與計算值的帶權重差異函數，根據參數優化的需求和實驗觀測值的種類，可將觀測值分組并賦予相應的權重，以達到模擬值與實測值的最佳匹配。本研究包含地下水位、土壤含水率2組目標函數，其公式如下：

(3)

式中：Φ為目標函數；Oi和Si分別表示第i個實驗觀測值和第i個模擬值；p表示實驗觀測值的個數；ω表示權重系數；下標wt和sm分別對應地下水位和土壤含水率。

1.4 評價方法

采用決定系數(R2)、納什系數(NSE)、均方根誤差(RMSE)、標準化均方根誤差(NRMSE)4種指標來評價模型的模擬效果：

(7)

決定系數(R2)用來衡量模擬值和實測值的相關密切程度，R2越接近1表示相關性越好；納什系數(NSE)可用來評價模型的模擬質量，NSE越接近1說明模擬質量越好；均方根誤差(RMSE)用來衡量模擬值與實測值之間的偏差，RMSE越接近0，表示模擬值與觀測值的差異越小，模擬效果越好；標準化均方根誤差(NRMSE)也是用來衡量模擬誤差的大小，NRMSE值越小，說明模擬效果越好。

2 結果與分析

參數自動優化共選取19個參數作為優化對象，包括15個土壤水力參數和4個水文地質參數，參數優化初值為試錯法調整后的結果，各個參數的取值范圍和優化結果見表1。

將試錯法和PEST參數優化后的地下水位模擬值和觀測井實測值進行比較(見圖5)，同時將模擬開始后第46、86、128、168天土壤剖面含水率的模擬值和實測值進行比較(見圖6)。將所有的地下水位觀測值與對應時刻的模擬值繪制在一張圖中(見圖7)，將所有的土壤含水率實測值與對應時刻的模擬值繪制在一張圖中(見圖7)，并計算評價模型效果的4個指標(見表2)。

采用2008年5月1日至11月31日期間的地下水位和土壤含水率實測值對模型進行驗證，將地下水位觀測值與對應時刻的模擬值繪制在一張圖中(見圖8)，將土壤含水率實測值與對應時刻的模擬值繪制在一張圖中(見圖8)，并計算評價模型效果的4個指標(見表2)。

圖5 地下水位模擬值與實測值對比Fig.5 Computed vs. observed groundwater table

圖6 不同時刻土壤剖面含水率模擬值與實測值對比Fig.6 Computed vs. observed soil moisture of different time

模擬期觀測變量試錯法R2NSERMSENRMSEPESTR2NSERMSENRMSE率定期地下水位0．7650．7310．5210．1040．8080．7860．4640．093土壤含水率0．4160．4080．0640．1370．5290．5280．0570．122驗證期地下水位0．8040．6030．5800．1290．8240．6410．5520．123土壤含水率0．3250．0320．0660．1740．4380．3950．0530．140

圖7 率定期地下水位與土壤含水率散點圖Fig. 7 Scatter plots of groundwater table and soil moisture for calibration period

圖8 驗證期地下水位與土壤含水率散點圖Fig. 8 Scatter plots of groundwater table and soil moisture for validation period

從圖5可以看出，地下水位模擬值與觀測值的變化趨勢基本一致；從圖6可以看出，模擬開始后不同時刻土壤剖面含水率模擬值與實測值較為接近；從圖7、圖8可以發現，地下水位和土壤含水率散點都較好地集中在1:1直線附近。

表2的模型效果評價指標顯示，PEST算法進行參數優化，提高了模型的模擬效果。

在率定期，地下水位的決定系數由0.765提高到0.808，提高了6%；納什系數由0.731提高到0.786，提高了8%；均方根誤差由0.521減小到0.464，減小了11%；標準化均方根誤差由0.104減小到0.093，減小了11%。

在率定期，土壤含水率的決定系數由0.416提高到0.529，提高了27%；納什系數由0.408提高到0.528，提高了29%；均方根誤差由0.064減小到0.057，減小了11%；標準化均方根誤差由0.137減小到0.122，減小了11%。

在驗證期，地下水位的決定系數由0.804提高到0.824，提高了2%；納什系數由0.603提高到0.641，提高了6%；均方根誤差由0.580減小到0.552，減小了5%；標準化均方根誤差由0.129減小到0.123，減小了5%。

在驗證期，土壤含水率的決定系數由0.325提高到0.438，提高了35%；納什系數由0.032提高到0.395，提高了11.3倍；均方根誤差由0.066減小到0.053，減小了20%；標準化均方根誤差由0.174減小到0.140，減小了20%。

同時從表2可以看出，耦合模型對地下水位的模擬效果要優于土壤含水率。從試錯法的結果看，地下水位在率定期和驗證期決定系數分別為0.765和0.804，納什系數分別為0.731和0.603；而土壤含水率在率定期和驗證期決定系數分別為0.416和0.325，納什系數分別為0.408和0.032。通過PEST參數優化后，地下水位在率定期和驗證期決定系數分別為0.808和0.824，納什系數分別為0.786和0.641；而土壤含水率在率定期和驗證期決定系數分別為0.529和0.438，納什系數分別為0.528和0.395。土壤含水率的各項評價指標都低于地下水位的相應指標。對土壤水運動模擬效果較差，可能的原因是本次模擬將整個土壤剖面視為均質土壤，與實際土壤剖面為非均質土壤的情況不符。

3 結 語

(1)本文采用試錯法和PEST優化算法相結合的方法對MODFLOW-HYDRUS耦合模型的水文地質參數和土壤水力參數進行優化，提高了模擬精度。與試錯法的結果相比，PEST參數優化后決定系數提高了2%～35%，納什系數提高了6～11.3倍，均方根誤差減小了5%～20%，標準化均方根誤差減小了5%～20%，其中驗證期土壤含水率的各項指標提高或減小的幅度均為最大，決定系數、納什系數、均方根誤差和標準化均方根誤差變化幅度分別為35%、11.3倍、20%和20%。

(2)MODFLOW-HYDRUS耦合模型模擬區域尺度上地下水、土壤水運動時，在地下水運動模擬上的表現要優于土壤水運動的模擬。從試錯法與PEST參數優化后結果來看，地下水位的決定系數在0.765～0.824之間，納什系數在0.603～0.786之間；而土壤含水率的決定系數在0.325～0.529之間，納什系數在0.032～0.528之間，明顯低于地下水位相對應的指標。導致土壤水運動模擬效果不佳的原因可能是將整個土壤剖面視為均質土壤，與實際土壤剖面非均質土壤的情況不符，在今后運用模型時可以考慮按照實際的非均質土壤剖面來設置。

□

[1] 楊欣坤, 王 宇, 趙蘭坡, 等. 土壤水動力學參數及其影響因素研究進展[J]. 中國農學通報, 2014,30(3):38-43.

[2] 薛禹群. 中國地下水數值模擬的現狀與展望[J]. 高校地質學報, 2010,16(1):1-6.

[3] Harbaugh A W, E R Banta, M C Hill, et al. MODFLOW-2000, the U.S. Geological Survey modular ground-water model user guide to modularization concepts and the ground-water fl ow process.USGS, Denver, CO,2000.

[4] Niswonger R G, D E Prudic, R S Regan. Documentation of the Unsaturated-Zone Flow (UZF1) package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods 6-A19. Available at http:∥pubs.usgs.gov/tm/2006/tm6a19/ (verified 14 Apr. 2008). USGS, Denver, CO,2006.

[5] Thoms R B, R L Johnson, R W Healy. User’s guide to the Variably Saturated Flow (VSF) process for MODFLOW. Techniques and Methods 6-A18. Available at http://pubs.usgs.gov/tm/2006/ tm6a18/ (verified 15 Apr. 2008). USGS, Reston, VA, 2006.

[8] 舒曉娟, 陳洋波, 黃鋒華, 等. PEST在WetSpa分布式水文模型參數率定中的應用[J]. 水文, 2009,29(5):45-49.

[9] 章四龍, 劉九夫. 通用模型參數率定技術研究[J]. 水文, 2005,25(1):9-12.

[10] Doherty J. PEST Model-Independent Parameter Estimation User Manual[M]. 5th Edition. Brisbane, Australia: Watermark Numerical Computing, 2004.

[11] Zyvoloski G, E Kwicklis, A A. Eddebbarh, et al. The site-scale saturated zone flow model for Yucca Mountain: calibration of different conceptual models and their impact on flow paths[J]. Journal of Contaminant Hydrology, 2003,62-63(2):731-750.

[12] 董艷輝, 李國敏, 郭永海, 等. 應用并行PEST 算法優化地下水模型參數[J]. 工程地質學報, 2010,18(1):140-145.

[13] Fang Q X, Green T R, Ma L, et al. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions[J]. Soil Science Society of America Journal, 2010,74(6):1 897-1 913.

[14] Allen R G, Pereira L S, Raes D, et al. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements[D]. Rome, Italy: United Nations Food and Agriculture Organization, 1998.

[15] van Genuchten, M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(44):892-898.

[16] Carsel R F, Parrish R S. Developing joint probability distributions of soil water retention characteristics[J]. Water Resources Research, 1988,24(5):755-769.

[17] Bahremand A, De Smedt F. Distributed hydrological modeling and sensitivity analysis in Torysa Watershed, Slovakia[J]. Water Resources Management, 2008,22(3):393-408.