林 青，徐紹輝

（青島大學 環(huán)境科學與工程學院，山東 青島 266071）

1 研究背景

目前，常用的描述非飽和土壤中水分運動的基本方程是Richards方程，該方程的應用需要獲取土壤水分特征曲線及導水率與土壤水勢關(guān)系的相關(guān)信息，由于兩者關(guān)系的高度非線性，在實驗室以及野外直接測量這些屬性參數(shù)通常是耗時和昂貴的。另外，由于大孔隙、空間異質(zhì)性、測量尺度等問題的存在，針對點尺度樣本的室內(nèi)穩(wěn)態(tài)試驗所得到的土壤水分運動參數(shù)往往不能直接反映天然條件下的田間尺度土壤水分運動特征［1］，因此根據(jù)原位觀測的土壤含水率和土壤水勢等數(shù)據(jù)反演非飽和帶土壤水分運動參數(shù)的研究一直備受關(guān)注?；谧顑?yōu)化理論的反演方法可依據(jù)特定的初、邊值條件給出最能匹配觀測數(shù)據(jù)的模型參數(shù)，為參數(shù)的估計提供了一種簡單便捷的間接估算方法，很多研究者采用該方法估計了土壤水力特征［2-3］。但是，采用該方法進行反演計算時要求模型滿足適定性（可識別性、唯一性和穩(wěn)定性），其取決于測量變量的個數(shù)、邊界條件、測量誤差及土壤類型等因素［4-5］。雖然已有研究表明，額外變量的測量（如基質(zhì)勢）可以在一定程度上減少不確定性問題［6-7］，但額外變量的測定在實驗室內(nèi)容易實現(xiàn)，而在野外是較難進行的。另外，實際上模型可能有多組參數(shù)解，即“異參同效”現(xiàn)象，由于反演問題的不確定性獲得的最優(yōu)解只是這多組參數(shù)中的任意一個。因此，最優(yōu)化反演方法不能很好地處理觀測數(shù)據(jù)和參數(shù)的不確定性，存在非唯一解問題；而如果再考慮普遍存在的土壤空間異質(zhì)性，如層狀土壤，將使參數(shù)反演問題維數(shù)升高，不確定性更為復雜。目前，關(guān)于參數(shù)不確定性的定量描述及其對模型預測不確定性的影響評價，在國內(nèi)外已成為研究熱點，這有助于我們更為深刻地理解和認識現(xiàn)實世界和模型系統(tǒng)特征之間的本質(zhì)區(qū)別。

近十幾年來，貝葉斯方法在模型參數(shù)估計及預報中得到廣泛應用。但是，應用貝葉斯方法的最大障礙在于后驗統(tǒng)計推斷中分母積分的計算，尤其是對于具有高維參數(shù)向量的復雜總體分布，即便是用數(shù)值積分方法，后驗分布計算式中分母的計算也比較困難，而迭代逼近算法使這一難題得到了有效解決，其中的馬爾科夫鏈蒙特卡羅（MCMC）方法就是較好的方法。已有研究表明，MCMC法對參數(shù)后驗分布的搜索，無論是搜索性能還是搜索效率，均表現(xiàn)出了獨特的優(yōu)越性，其在貝葉斯推斷的參數(shù)識別過程中得到廣泛應用［8］。

貝葉斯方法可以最大限度地利用現(xiàn)有資料和已知信息，通過似然函數(shù)用觀測資料來修正輸入的參數(shù)分布，使得最終的參數(shù)分布與實際相吻合。因此，本文采用Hydrus_1d模擬自然條件下包氣帶土壤水分運動，基于貝葉斯推斷，選取基于DREAMZS（Differential Evolution Adaptive Metropolis）抽樣的MCMC方法，對大沽河流域農(nóng)田層狀分布土壤的水力特性參數(shù)進行識別，獲得參數(shù)后驗分布的抽樣樣本，并對其模擬預測的不確定性進行分析，探討造成模擬結(jié)果不確定性的因素，有助于提高模型應用的精度和預測結(jié)果的可靠性。

2 研究方法

2.1 試驗數(shù)據(jù) 試驗區(qū)域位于青島即墨移風店鎮(zhèn)上泊村（36°33'N，120°12'E），地處大沽河中游河谷沖積平原，面積約8 hm2，東西長275.0 m，南北長290.0 m。試驗區(qū)農(nóng)作物主要為小麥-玉米輪作，土壤類型為棕壤，土壤結(jié)構(gòu)分層明顯，結(jié)合地質(zhì)工程勘察結(jié)果，將土壤剖面劃分為5層，對應的深度分別為0～30、30～60、60～100、100～200和200～450 cm，其中0～100 cm為耕土層，100～200 cm為粉質(zhì)黏土層，200～450 cm為中細砂。對耕層土壤用圓盤滲透儀（Disc permeameter）進行飽和導水率Ks的原位測定，人工開挖剖面，用環(huán)刀取原狀土壤進行水分特征曲線和容重的測定，并取擾動土樣進行土壤粒徑分析。采用土鉆法對土壤含水量進行長期監(jiān)測，監(jiān)測點分別設在20、40、80和150 cm處，監(jiān)測時間從2013年4月1日至2014年10月5日，2013年每7 d監(jiān)測1次，生育期、雨后或灌溉后加密監(jiān)測1次，2014年每10 d監(jiān)測一次。為了降低土壤含水量觀測數(shù)據(jù)的誤差，共設置8個觀測點（南北各3個，中部2個），取其平均值進行參數(shù)反演。地下水位每5 d觀測一次。2013年為模型的率定期，利用這一階段的監(jiān)測數(shù)據(jù)進行參數(shù)識別和不確定性分析，2014年為驗證期，用于預測計算和模型驗證。

2.2 作物生長條件下土壤水分運動模型 本文采用經(jīng)典的Richards方程模擬變飽和土壤中垂向一維水流運動，如式（1）所示，利用Hydrus_1d模型對方程進行求解。

式中：θ為土壤體積含水量，cm3/cm3；h為壓力水頭，cm；K（h）為非飽和土壤導水率，cm/d；S（z，t）為t時刻z深度處根系吸水速率，表示單位時間單位土體的根系吸水量，cm3/(cm3·d)，采用van Genu?chten宏觀根系吸水模型計算［9］。其中，t為時間，d；z為土壤深度，坐標向下為正，cm。土壤水分脅迫函數(shù)選取S-Shape方程，方程中只有兩個參數(shù)p0和p50，對冬小麥和夏玉米而言，p0≈3；p50為作物潛在蒸騰率減少50%時相應的土壤基質(zhì)勢，cm。p50與作物生理特性有關(guān)，其絕對值越大，作物耐旱吸水能力越強，這里取為-1600 cm。

土壤的持水和導水特征由VGM（van Genuchten-Mualem）模型來計算，

式中：θr、θs分別為殘余含水量和飽和含水量，cm3/cm3；α（1/cm）、n和m為曲線形狀參數(shù)，m=1-1/n；為有效飽和度。在模擬過程中l(wèi)取經(jīng)驗值0.5。因此，在模型中僅對各層土壤的θr，θs，α，n和Ks進行率定。

模擬期內(nèi)試驗區(qū)地下水位變化情況如圖1所示，地形標高為17.8 m，地下水埋深為2.7～4.2 m，故取模擬土層厚度為4.5 m，剖分為不等間隔的100個單元。

圖1 2013—2014年試驗區(qū)地下水位

模型初始條件為：

式中：h0（z）為初始壓力值，設為-100 cm，目的是為了避免由于層狀非均質(zhì)土壤含水量不連續(xù)造成模型不收斂，而且為了降低模擬結(jié)果對初始條件的依賴性，在模型校準期設置30 d預熱期（spin-up）；L為土體深度，cm。

邊界條件為：

式中：hL（t）為下邊界處（L=4.5m）的壓力水頭，模擬期間地下水最大埋深為4.2 m，模擬土層的底部始終位于地下水面以下，故以變水頭作為模型的下邊界條件；上邊界采用通量已知的第二類邊界條件，在模擬時段內(nèi)逐日輸入通過上界面的變量值，包括降水量、灌溉量、作物潛在蒸散發(fā)量，E（t）表示土壤水分最大蒸發(fā)或最大入滲強度，cm/d；ha為地表允許的最小壓力水頭，設定為-16 000 cm。降雨量數(shù)據(jù)來自移風鎮(zhèn)雨量站，冬小麥和夏玉米的潛在蒸散發(fā)率參照FAO推薦的Hargreaves公式來求得，如圖2所示，其中作物系數(shù)、葉面積指數(shù)、根系深度等參照廖凱華在大沽河流域土壤水資源評價研究中的數(shù)據(jù)［10］。

2.3 基于Bayes理論的參數(shù)反演

式中：y為所有觀測點的土壤含水量數(shù)據(jù)；p（y|θ）為給定參數(shù)θ時數(shù)據(jù)y被觀測到的概率，即似然函數(shù)；p（y）為觀測數(shù)據(jù)的產(chǎn)生概率，該值的計算比較困難，通常的做法是把其從方程中去掉，則：

圖2 2013—2014年試驗區(qū)冬小麥-夏玉米逐日潛在蒸散發(fā)率及降雨量

假設各觀測點的數(shù)據(jù)相互獨立，則所有數(shù)據(jù)對應的似然函數(shù)是各觀測點數(shù)據(jù)對應的似然函數(shù)的乘積，則似然函數(shù)可表示為：

式中：n為土壤剖面觀測點的個數(shù)。為簡化計算和達到數(shù)值的穩(wěn)定性，文中似然函數(shù)取其對數(shù)值。

2.3.2 基于DREAMZS抽樣的MCMC模擬 MCMC方法是在貝葉斯理論框架下，建立一個平穩(wěn)分布為所求后驗分布p（θ|y）的隨機樣本，應用中若產(chǎn)生轉(zhuǎn)移軌跡的馬爾科夫鏈的建議分布選取不適當，會造成該方法收斂速度較慢。DREAM算法作為一種自適應抽樣方法，結(jié)合了差分進化算法和自適應Me?tropolis算法的優(yōu)點，可同時運行多條馬爾科夫鏈進行全局搜索，并可以自動調(diào)整建議分布的范圍和方向，在復雜高維、非線性和多峰目標分布問題中表現(xiàn)出極高的效率［13-14］。抽樣過程中，假設第i條鏈的當前狀態(tài)用d維（待估參數(shù)個數(shù)）向量θi（i=1，…，N）表示，則其候選樣本可根據(jù)下式計算：

式中：δ為用于產(chǎn)生候選樣本的平行鏈對數(shù)；其中d′為更新的參數(shù)維數(shù)；r1（j）和r2（h）∈{1，…，N}，且r1（j）≠r2（h）≠i；e和ε為隨機產(chǎn)生的很小的數(shù)，且e服從均勻分布Ud（-b，b），ε服從于正態(tài)分布Nd（0，b*），b、b*為定義的極小值。

2.3.3 模型參數(shù)先驗分布的確定 本課題組曾對試驗區(qū)所在大沽河流域上100個樣點的土壤飽和導水率Ks進行了測定，統(tǒng)計結(jié)果顯示其服從對數(shù)正態(tài)分布。已有研究表明，VGM模型中的形狀特性參數(shù)α和n分別服從對數(shù)正態(tài)分布和正態(tài)分布，而其中的兩個水力特性參數(shù)（θs和θr）都遵從正態(tài)分布［17-19］，因此，在這里我們設定參數(shù)α服從對數(shù)正態(tài)分布，n、θr和θs服從正態(tài)分布，如表1所示。對于0～200 cm深度內(nèi)的土壤水分運動參數(shù)，其先驗分布由試驗區(qū)8個樣點32個土樣的顆粒組成和容重通過土壤轉(zhuǎn)換函數(shù)（Rosetta）獲得，由于200～450 cm缺少先驗信息，參數(shù)設為均勻分布，取值范圍根據(jù)經(jīng)驗并參考相關(guān)文獻給出［9，20］。

表1 不同深度處土壤水分運動參數(shù)的先驗分布和極限值

2.4 不確定性預測區(qū)間的評價指標 通常評價預測置信區(qū)間優(yōu)良性指標為預測區(qū)間包含的實測數(shù)據(jù)點比例（PCI）和平均相對區(qū)間寬度（ARIL）。PCI值越大，預測區(qū)間覆蓋實測數(shù)據(jù)點的比例越高，模型可靠性越大。ARIL值越小，預測的平均相對區(qū)間寬度越窄，精度越高。研究表明，隨著PCI的增高，ARIL也在增大，即兩者往往是相互矛盾的。因此，在前兩個指標的基礎上，本文還采用單位平均相對寬度所包含的實測點數(shù)據(jù)比例（PUCI）這個綜合指標［21］，來評價模型的可信性，PUCI值越大，置信區(qū)間的性能越優(yōu)良。ARIL和PUCI表達式如下：

式中：LimitUpper，t、LimitLower，t分別為t時刻模擬值95%置信區(qū)間的上限和下限；K為觀測值的個數(shù)；Qobs，t為t時刻的觀測值。

3 結(jié)果與討論

3.1 參數(shù)的后驗分布和統(tǒng)計特征 采用DREAMZS抽樣方法，基于計算時間（約5 h）與樣本收斂性（R＜1.2）的考慮，設定3條馬爾科夫鏈，每條鏈最大樣本數(shù)為50 000個，按1∶10進行稀釋。圖3為收斂準則R隨進化代數(shù)的變化曲線，不同顏色的曲線代表了不同參數(shù)采樣結(jié)果的收斂性變化過程?？梢钥闯?，當平行馬爾科夫鏈進化代數(shù)達到45 000代時，各參數(shù)的R值均已小于1.2的基準值，說明各參數(shù)都收斂于穩(wěn)定的后驗分布。本文取采樣結(jié)果最后25%的樣本（3750個）作為參數(shù)后驗分布樣本進行統(tǒng)計分析。

圖3 DREAMzs算法采樣過程的收斂性

表2 參數(shù)后驗分布的均值、標準差、變異系數(shù)、最大后驗概率參數(shù)（MAP）及95%的置信區(qū)間

表2列出了水分運動參數(shù)的后驗分布特征值及95%的置信區(qū)間。從后驗均值來看，各層土壤的參數(shù)之間差別較明顯，尤其是Ks，不同土層之間相差接近10倍，反映出土壤剖面持水與導水特征的變異性。從置信區(qū)間和變異系數(shù)來看，Ks最不敏感，置信區(qū)間較寬，變異系數(shù)較大，θs最為敏感，可識別性較高，置信區(qū)間較窄，變異系數(shù)較小。研究區(qū)8個點位耕層土壤Ks的原位分層測定結(jié)果分別為5.69～15.12、8.47～23.71和5.91～19.74 cm/d，實測值全部位于95%置信區(qū)間內(nèi)，說明參數(shù)Ks后驗分布的可靠性較高；但是無論是實測值還是后驗分布的標準差都較大，說明Ks的不確定性較大。對比分析室內(nèi)水分特征曲線反演獲得的水分特征參數(shù)（如表2試驗值所示），發(fā)現(xiàn)反演得到的θs與田間獲得的參數(shù)相差不大，認為可把室內(nèi)獲得的θs用于田間水分的模擬，這與Kumar等的研究結(jié)果相似［22］；其它參數(shù)則差別較大，這主要是由于田間土壤的大孔隙、空間變異、尺度效應等因素的影響，因此，在使用這些參數(shù)對田間土壤水分進行模擬預測時要慎重考慮其適用性。

3.2 最大后驗概率參數(shù)（MAP）模擬分析 采用后驗分布中的最大概率參數(shù)組對2013—2014年（2013年為率定期和2014年為驗證期）土壤水分動態(tài)變化進行模擬分析，不同深度處土壤含水量的模擬值和實測值如圖4所示。從圖中可以看出，模擬結(jié)果基本能反映出各深度處土壤水分的動態(tài)變化趨勢，相較于率定期而言，驗證期的水分變化相對平緩，但是模型仍能較好的模擬預測驗證期水分的變化。另外，雖然率定期缺少220 cm深度處的實測含水量數(shù)據(jù)，但仍能獲得該深度處土壤的水力學參數(shù)，且驗證期采用該參數(shù)的模擬值與實測值吻合程度較好。然而，由于土壤含水量受到降水、蒸發(fā)等外界因素干擾影響較大，而且土壤水分狀況實測值還受到一些不確定因素的影響，加上一些儀器和人為的誤差，模擬值和實測值存在一定的偏差，尤其是降雨后土壤水分的實測值和模擬值差別較大。就土壤含水量隨土層變化而言，受降雨、灌水及蒸散發(fā)的影響，上層土壤含水量變化劇烈，而深層土壤含水量變化相對平緩且接近飽和。

圖4 率定期（2013）和驗證期（2014）不同深度處土壤含水量模擬值和實測值的對比

采用納什系數(shù)（NSE）、均方根誤差（RMSE）、平均絕對誤差（MAE）和偏差百分比（PBIAS）4個評價指標來綜合衡量模擬的效果。NSE越接近于1，MAE、RMSE和PBIAS越接近于0，說明模擬效果越好，評價指標值如表3所示。

就納什系數(shù)NSE而言，NSE＞0.5時為可接受模擬值。通過計算，不同深度處土壤水分模擬值和實測值的NSE均在可接受的范圍內(nèi)（由于率定期缺乏220 cm深度處的實測數(shù)據(jù)，故驗證期該處的NSE較低為0.498，但是其模擬結(jié)果仍能反應土壤水分的變化趨勢）。另外，計算得到的RMSE、MAE和PBIAS指標值都較小，其中RMSE小于0.02 cm3/cm3，MAE小于0.015 cm3/cm3，PBIAS處于-1.57%～4.32%之間，接近于0，皆說明模擬值和實測值吻合良好。所以，基于DREAMZS方法的Bayes估計得到的參數(shù)可實現(xiàn)對田間層狀分布土壤水分動態(tài)的良好模擬，尤其是在實測數(shù)據(jù)不足（缺乏底層土壤含水量實測數(shù)據(jù)）的情況下，其仍能實現(xiàn)參數(shù)的識別及土壤水分的模擬預測。從不同土層來看，20、40、80 cm剖面處模擬值和實測值的NSE相差不大（0.613～0.661），而150 cm深度處的NSE則較小，主要是因為該深度正處于粉質(zhì)黏土層，黏土水分運動具有明顯的滯后性［23］，而模型沒有考慮這種滯后性。RMSE、MAE和|PBIAS|3個評價指標都表現(xiàn)出隨著土層深度增大而減少的趨勢，所以，總體而言模型對深層土壤含水量的模擬精度高于表層，主要原因在于深層土壤含水量受外界條件影響較小，變化相對表層比較平緩。

表3 不同深度處模擬值與實測值吻合度的評價指標值

3.3 模擬結(jié)果的不確定性分析 不確定性理論認為，對存在不確定性的事件進行數(shù)值模擬時，一條最優(yōu)的擬合曲線并不能真實反映客觀事件本身的規(guī)律，科學的模擬應該是確定一組合理的取值區(qū)間。因此，針對率定期（2013）模型預測結(jié)果的不確定性，用采樣的后25%的參數(shù)集合運行土壤水分運動模型，根據(jù)模擬的結(jié)果分析模型預測的參數(shù)不確定性及總不確定性。這里參數(shù)不確定性指由參數(shù)非唯一性所造成的模擬結(jié)果的不確定性，主要源自于模型的復雜性、參數(shù)間的相互作用等；總不確定性是指由參數(shù)的不確定性、模型結(jié)構(gòu)的不確定性及輸入數(shù)據(jù)的不確定性等所引起的模擬預測的不確定性，其由模擬結(jié)果加上一個正態(tài)分布的隨機誤差項獲得，隨機誤差的標準差為最大后驗概率參數(shù)（MAP）經(jīng)Hydrus_1d模擬得到的均方根誤差。圖5為95%置信區(qū)間下模型預測值的上下限及土壤含水量的實測值，黑色區(qū)間為參數(shù)不確定性所引起的模型預測結(jié)果的不確定性，灰色區(qū)間為由參數(shù)、模型結(jié)構(gòu)、邊界條件等所引起的模型預測的總不確定性，點為實測值。

圖5 不同深度處土壤含水量模擬值的95%的置信區(qū)間和實測值

從圖5可以看出，除了220 cm處的預測區(qū)間較寬外，其它深度處的預測區(qū)間相對較窄，參數(shù)不確定性區(qū)域（黑色區(qū)間帶）沒有包含所有觀測值，總的不確定性區(qū)域（灰色區(qū)間帶）包含了所有的觀測值，但是部分時間段內(nèi)預測的不確定性范圍較大，說明在目前進行土壤水分預測時模型結(jié)構(gòu)、邊界條件及輸入數(shù)據(jù)等造成的不確定性是值得考慮的。對比降雨數(shù)據(jù)（圖1），可以看出，模型低估了由于降雨所引起的土壤水分變化預測的不確定性，而高估了無雨期間土壤水分變化的不確定性，表明在進行田間土壤水分模擬的不確定性分析時，采用標準可加誤差模型可能是不適當?shù)?，這在Jiang等的研究中也有提到［21］。

模型預測置信區(qū)間的評價指標值如表4所示。從表中可以看出，由模型參數(shù)的不確定性所造成的模擬結(jié)果的不確定性較小，預測帶平均相對寬度為0.026～0.135，預測帶對實測值的覆蓋度為0.545～0.773。綜合評價指標值PUCI為4.407～26.385，且隨著土層深度的增加而增大，表明隨著土壤含水量預測深度的增加，模型預測的性能（可靠性和精度的綜合表現(xiàn)）越高。模型預測結(jié)果的總不確定性明顯高于參數(shù)預測結(jié)果的不確定性，ARIL較大為0.086～0.472，PUCI較小為2.013～11.047，說明模型結(jié)構(gòu)等不確定性造成了很大的模型預測結(jié)果的不確定性。相對于模型參數(shù)和觀測數(shù)據(jù)而言，基于現(xiàn)有科學認知體系構(gòu)建的模型結(jié)構(gòu)往往是模擬過程中不確定性的根本來源［24］，所以，未來要提高模型對土壤水分動態(tài)模擬的精度，需要對模型結(jié)構(gòu)做進一步的改進。

表4 土壤含水量模擬值95%置信區(qū)間評價指標值

4 結(jié)論

本文基于貝葉斯理論，采用基于自適應差分演化抽樣的MCMC方法，對大沽河流域農(nóng)田層狀分布土壤的水力特性參數(shù)進行識別，并分析了土壤水分模擬預測的不確定性，得到了以下幾點結(jié)論：

（1）采用Bayes方法對模型參數(shù)進行識別，不僅得到了確定的模型輸出（最大后驗概率參數(shù)組的模擬結(jié)果），而且還獲得由參數(shù)不確定性和總不確定性所造成的模型預測的區(qū)間范圍，相比最優(yōu)化方法獲得的一組結(jié)果，Bayes方法的模擬結(jié)果給出了更多的信息，提高了預測的可靠性。

（2）土壤水分運動參數(shù)識別結(jié)果顯示，Ks最不敏感，置信區(qū)間較寬，θs最為敏感，置信區(qū)間較窄，可識別性較高。通過與室內(nèi)試驗反演獲得的參數(shù)相比，發(fā)現(xiàn)室內(nèi)反演得到的θs可用于田間水分的模擬，而其他參數(shù)差別較大，在進行田間水分模擬時要謹慎使用。

（3）最大后驗概率參數(shù)模擬結(jié)果基本能反應出各深度處土壤水分的動態(tài)變化趨勢，即使是在缺失率定期（220 cm深度處）實測數(shù)據(jù)的情況下，其驗證期的模擬值和實測值的吻合程度也較好。所以，基于DREAMZS方法的Bayes估計得到的參數(shù)可實現(xiàn)對田間層狀分布土壤水分動態(tài)的良好模擬。

（4）隨著土壤含水量預測深度的增加，PUCI值越大，表明模型預測的性能（可靠性和精度的綜合表現(xiàn)）越高，未來要提高模型預測的精度，需要對模型結(jié)構(gòu)做適當?shù)男薷耐晟啤?/p>

參 考 文 獻：

［1］W?HLING T，VRUGT J A.Combining multiobjective optimization and Bayesian model averaging to calibrate forecast ensembles of soil hydraulic models［J］.Water Resources Research，2008，44（12）：181-198.

［2］查元源，周發(fā)超，楊金忠.一種由土壤剖面含水率估算土壤水力參數(shù)的方法［J］.水利學報，2011，42（8）：883-891.

［3］李春友，任理，李保國.利用優(yōu)化方法求算Van Genuchten方程參數(shù)［J］.水科學進展，2001，12（4）：473-478.

［4］CARRERA J，NEUMAN S P.Estimation of aquifer parameters under transient and steady state conditions：1.maximum likelihood method incorporating prior information［J］.Water Resources Research，1986，22（2）：199-210.

［5］RUSSO D，BRESLER E，SHANI U，et al.Analysis of infiltration events in relation to determining soil hydraulic properties by inverse problems methodology［J］.Water Resources Research，1991，27（6）：1361-1373.

［6］KOOL J B，PARKER J C.Analysis of the inverse problem for transient unsaturated flow［J］.Water Resources Re?search，1988，24：817-830.

［7］PARKER J C，KOOL J B，van GENUCHTEN M T.Determining soil hydraulic properties from one-step outflow experiments by parameter estimation.Part 2：experimental studies［J］.Soil Science Society of America Journal，1985，49：1354-1359.

［8］朱嵩，毛根海，劉國華，等.改進的MCMC方法及其應用［J］.水利學報，2009，40（8）：1019-1023.

［9］?IM?NEK J，?EJNA M，SAITO H，et al.The HYDRUS-1D Software Package for Simulating the One-Dimen?sional Movement of Water，Heat，and Multiple Solutes in Variably-Saturated Media：version 4.17［EB/OL］.［2017-09-25］.https：//www.pc-progress.com/Downloads/Pgm_hydrus1D/HYDRUS1D-4.17.pdf.

［10］廖凱華.大沽河流域土壤水資源評價及農(nóng)業(yè)節(jié)水灌溉模式研究［D］.青島：青島大學，2009.

［11］W?HLING T，VRUGT J A.Multiresponse multilayer vadose zone model calibration using Markov chain Monte Carlo simulation and field water retention data［J］.Water Resources Research，2011，47（4）：W04510，doi：10.1029/2010WR009265.

［12］KAVETSKI D，KUCZERA G，F(xiàn)RANKS S W.Bayesian analysis of input uncertainty in hydrological modeling：1.Theory［J］.Water Resources Research，2006，42（3）：W03407，doi：10.1029/2005WR004368

［13］SCHARNAGL B，VRUGT J A，VEREECKEN H，et al.Information content of incubation experiments for inverse estimation of pools in the Rothamsted carbon model：a Bayesian perspective［J］.Biogeosciences Discussions，2009，6（5）：1736-1751.

［14］SCHARNAGL B，VRUGT J A，VEREECKEN H，et al.Inverse modelling of in situ soil water dynamics：investi?gating the effect of different prior distributions of the soil hydraulic parameters［J］.Hydrology&Earth System Sci?ences，2011，15（10）：3043-3059.

［15］LALOY E，VRUGT J A.High-dimensional posterior exploration of hydrologic models using multiple-try DREAM（ZS）and high performance computing［J］.Water Resources Research，2012，48（1）：182-205.

［16］GELMAN A，RUBIN D B.Inference from iterative simulation using multiple sequences［J］.Statistical Science，1992，7（4）：457-472.

［17］MOREIRA P H S，van GENUCHTEN M T，ORLANDE H R B，et al.Bayesian estimation of the hydraulic and solute transport properties of a small-scale unsaturated soil column［J］.Journal of Hydrology&Hydromechan?ics，2016，64（1）：30-44.

［18］RUSSO D，BOUTON M.Statistical analysis of spatial variability in unsaturated flow parameters［J］.Water Re?sources Research，1992，28（7）：1911-1925.

［19］楊金忠.多孔介質(zhì)中水分及溶質(zhì)運移的隨機理論［M］.北京：科學出版社，2000.

［20］山東省環(huán)境水文地質(zhì)總站.山東省青島市大沽河地下水庫勘察報告［R］.1987.

［21］JIANG S，JOMAA S，BüTTNER O，et al.Multi-site identification of a distributed hydrological nitrogen model using Bayesian uncertainty analysis［J］.Journal of Hydrology，2015，529（20）：940-950.

［22］KUMAR S，SEKHAR M，REDDY D V，et al.Estimation of soil hydraulic properties and their uncertainty：com?parison between laboratory and field experiment［J］.Hydrological Processes，2010，24（23）：3426-3435.

［23］王金生，楊志峰，陳家軍，等.包氣帶土壤水分滯留特征研究［J］.水利學報，2000（2）：1-7.

［24］BECK M B.Water quality modeling：A review of the analysis of uncertainty［J］.Water Resources Research，1987，23（8）：1393-1442.