優(yōu)化法與貝葉斯估計(jì)法在非飽和水力參數(shù)反演中的比較

柯鳳喬，滿俊，曾令藻，吳勞生

（浙江大學(xué)環(huán)境與資源學(xué)院土水資源與環(huán)境研究所，杭州310058）

優(yōu)化法與貝葉斯估計(jì)法在非飽和水力參數(shù)反演中的比較

柯鳳喬，滿俊，曾令藻*，吳勞生

（浙江大學(xué)環(huán)境與資源學(xué)院土水資源與環(huán)境研究所，杭州310058）

非飽和土壤水分運(yùn)動(dòng)模型在農(nóng)業(yè)生產(chǎn)和環(huán)境保護(hù)等領(lǐng)域中具有重要的指導(dǎo)意義.非飽和土壤水力參數(shù)的準(zhǔn)確獲取是利用模型進(jìn)行可靠預(yù)測(cè)的前提.傳統(tǒng)的參數(shù)反演研究大多基于優(yōu)化方法，只能獲取一組最優(yōu)參數(shù)，不能量化其中的不確定性.最新發(fā)展的一種基于貝葉斯參數(shù)估計(jì)理論的馬爾科夫鏈蒙特卡羅算法（Markov Chain Monte Carlo，MCMC）—DREAM（ZS），可以有效進(jìn)行參數(shù)反演，且準(zhǔn)確量化不確定性.我們?cè)讷@取水頭觀測(cè)值的基礎(chǔ)上，利用MCMC和Levenberg-Marquardt（LM）非線性優(yōu)化算法分別對(duì)非飽和土壤水力參數(shù)進(jìn)行反演，通過數(shù)值模擬與一維沙柱入滲實(shí)驗(yàn)比較了2種方法對(duì)于分層異質(zhì)的水力參數(shù)估計(jì)與水頭預(yù)測(cè)的準(zhǔn)確性.結(jié)果表明:1）LM優(yōu)化方法使用廣泛，求解速度較快，但受參數(shù)初始值影響較大，預(yù)測(cè)結(jié)果與觀測(cè)值存在一定的偏差，同時(shí)由于只能給出一個(gè)單一的反演結(jié)果，無法量化結(jié)果的不確定性.2）與基于優(yōu)化的反演方法相比，MCMC反演方法不僅能夠更好地得出參數(shù)的單一估計(jì)值，同時(shí)預(yù)測(cè)結(jié)果與觀測(cè)值也具有更好的一致性；更重要的是可以給出未知參數(shù)的后驗(yàn)分布，從而準(zhǔn)確量化非飽和水力參數(shù)的不確定性，避免了基于單一參數(shù)反演結(jié)果進(jìn)行預(yù)測(cè)的風(fēng)險(xiǎn).但是相對(duì)于LM優(yōu)化方法，MCMC計(jì)算量大大增加.

貝葉斯算法；參數(shù)估計(jì)；非飽和流；不確定性

了解非飽和土壤中水分的運(yùn)動(dòng)在農(nóng)業(yè)生產(chǎn)、環(huán)境保護(hù)等方面具有十分重要的意義.土壤水分運(yùn)動(dòng)方程，即Richards方程［1］，可以描述土壤中水分的運(yùn)動(dòng)規(guī)律.目前，研究者已經(jīng)開發(fā)出了很多用于非飽和土壤中水分溶質(zhì)運(yùn)移的數(shù)值模型，其中影響力最大，使用最廣泛的是HYDRUS模型［2］.HYDRUS模型綜合考慮了水分運(yùn)動(dòng)、熱量遷移、溶質(zhì)運(yùn)移、氣體運(yùn)移、化學(xué)反應(yīng)、作物吸收等多個(gè)復(fù)雜的過程，而且可以處理各種復(fù)雜初始與邊界條件下的一維到三維的問題［34］.介于上述原因，本研究也采用HYDRUS來進(jìn)行非飽和土壤水分運(yùn)動(dòng)過程的模擬.

運(yùn)用非飽和土壤水分運(yùn)動(dòng)方程進(jìn)行合理預(yù)測(cè)的前提是非飽和水力參數(shù)的準(zhǔn)確獲取.但在實(shí)際中，這些參數(shù)往往難以直接測(cè)量，需要從更容易獲取的水分狀態(tài)（例如水頭和含水量等）中反演得到.此外，這些參數(shù)往往具有空間變異性［5］，為參數(shù)的準(zhǔn)確反演帶來了難度.

傳統(tǒng)的參數(shù)反演方法大多基于優(yōu)化算法.以HYDRUS為例，其自帶的參數(shù)反演工具是基于Levenberg-Marquardt（LM）非線性優(yōu)化算法，已經(jīng)得到廣泛的應(yīng)用.但這種反演方法都只能給出1組最優(yōu)解，無法對(duì)內(nèi)在的不確定性進(jìn)行量化.此外，該組解也可能是局部最優(yōu)解，而非全局最優(yōu)解.為了對(duì)參數(shù)不確定性進(jìn)行準(zhǔn)確的量化，可以采用隨機(jī)的參數(shù)估計(jì)方法.近年來，在土壤水力特性參數(shù)估計(jì)領(lǐng)域，應(yīng)用最為廣泛的2種方法是集合卡爾曼濾波（Ensemble Kalman Filter，En KF）［6］和馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）［7］.前者基于最優(yōu)線性估計(jì)理論，后者基于貝葉斯參數(shù)估計(jì)理論.

在貝葉斯框架下，所有的待定參數(shù)都被視為隨機(jī)變量，都用一個(gè)相應(yīng)的分布函數(shù)（也稱為概率密度函數(shù)）來表征［8］.當(dāng)貝葉斯原理用于參數(shù)估計(jì)時(shí)，觀測(cè)數(shù)據(jù)對(duì)模型與參數(shù)的影響由后驗(yàn)分布來表示.在獲得觀測(cè)數(shù)據(jù)后，關(guān)于模型參數(shù)的統(tǒng)計(jì)信息，如常用的均值、方差、分位值、置信區(qū)間等都可以從后驗(yàn)分布得到.然而當(dāng)系統(tǒng)復(fù)雜且非線性時(shí)，往往不存在后驗(yàn)分布的解析式，于是較常用的一種方法是直接生成滿足后驗(yàn)分布的樣本［9］.最常用的方法是MCMC［10］.近年來比較流行的高效MCMC算法包括DRAM［11］，DREAM［12］.MCMC算法在水資源領(lǐng)域中獲得了非常廣泛的應(yīng)用，已有研究者利用MCMC算法來估計(jì)土壤水力特性的參數(shù).例如JASPER，等［13］分析了土壤水力特性參數(shù)的先驗(yàn)分布對(duì)MCMC反演效果的影響.

目前，已有的研究通常只針對(duì)均質(zhì)非飽和水流問題，即非飽和水力參數(shù)在空間上保持不變.而實(shí)際中的孔隙介質(zhì)則可能是分層甚至空間上連續(xù)性變化.同時(shí)，也缺乏在孔隙介質(zhì)異質(zhì)下貝葉斯參數(shù)估計(jì)方法與傳統(tǒng)優(yōu)化方法的比較研究.本文研究分層異質(zhì)時(shí)非飽和水流問題中的參數(shù)反演問題，通過數(shù)值模擬與沙柱實(shí)驗(yàn)，對(duì)比研究貝葉斯MCMC方法與傳統(tǒng)的LM優(yōu)化參數(shù)反演方法.

1 材料與方法

1.1 實(shí)驗(yàn)材料

＜0.25 mm，0.25～0.35 mm，0.35～0.5 mm 3種粒徑的沙子，樹脂玻璃制作的沙柱柱體，100 m L標(biāo)準(zhǔn)環(huán)刀，電子天平，秒表，溫度計(jì)，熱熔膠等.

1.2 實(shí)驗(yàn)方法

本實(shí)驗(yàn)通過圖1所示裝置，模擬水分在一維沙柱中的入滲過程.

沙柱的柱體為樹脂玻璃制作而成.其內(nèi)徑為10 cm，高度100 cm，柱子底部裝有陶瓷板（飽和導(dǎo)水率: 0.065 3 cm/d；厚度:0.714 cm）；沙柱由＜0.25 mm，0.25～0.35 mm，0.35～0.5 mm粒徑的沙子自下而上填裝，每層厚度為30 cm.沙柱的填裝密度見表1.

根據(jù)沙柱各層的填裝密度填裝環(huán)刀樣，并按以下方法測(cè)定環(huán)刀樣的飽和含水量、飽和導(dǎo)水率.

用電子天平稱取烘干后的填裝環(huán)刀質(zhì)量以及濾紙質(zhì)量（m1），再將其放置在水中浸泡24 h，此時(shí)環(huán)刀中土壤充滿水，用毛巾擦掉環(huán)刀外沾水，立即稱質(zhì)量（m2）.

飽和含水量/%=（m2-m1）/（ρ水×V）×100.

圖1 沙柱裝置圖Fig.1 Sketchofthesandcolumnsetup

表1 不同粒徑沙子的飽和水力參數(shù)Table1 Saturatedhydraulicpropertiesofsandwithdifferent particlesizes

式中:ρ水為水的密度，V為環(huán)刀容積.根據(jù)計(jì)算好的填裝密度填裝環(huán)刀，在環(huán)刀上再套1個(gè)空的環(huán)刀，接口處用熱熔膠封好，防止水分從接口處漏出，在環(huán)刀土表放1張大小適中的濾紙，再將其放置在水中浸泡24h，此時(shí)環(huán)刀中土壤充滿水.取出后將其放在漏斗上，漏斗下用50-mL燒杯盛接并向環(huán)刀內(nèi)加水，水層厚度保持為5cm.加水后待滲透速率接近穩(wěn)定開始計(jì)時(shí)接水，不同樣品設(shè)置不同的時(shí)間間隔，同時(shí)注意維持水頭高相同（5cm）.最后測(cè)定水溫.當(dāng)3次讀數(shù)一樣時(shí)，判定為飽和導(dǎo)水率.同一樣品做3次平行實(shí)驗(yàn).

飽和導(dǎo)水率=［10×穩(wěn)定水量×水柱高度（即環(huán)刀高）］/時(shí)間間隔×環(huán)刀橫截面積×［水頭高度（水層高）+土柱高度］.

測(cè)得環(huán)刀樣的飽和含水量和飽和導(dǎo)水率的數(shù)據(jù)見表1.

水分入滲試驗(yàn)從上部注水，在柱子側(cè)邊距沙柱表面10、30、45、70和85cm處設(shè)置水頭觀測(cè)點(diǎn)；下端開口，用于出水和連接真空泵.

水頭通過陶土頭和壓力傳感器（Honeywell 26PCBFA6G）組裝起來測(cè)得.張力計(jì)陶土頭部分進(jìn)氣值為30kPa，壓力傳感器的精度為0.4%.因此張力計(jì)的量程為-30～30kPa，測(cè)量精度為0.4%.壓力傳感器與數(shù)據(jù)采集器相連，每10min自動(dòng)采集1次各觀測(cè)點(diǎn)的水頭值.

沙柱下端的負(fù)壓設(shè)置分為恒壓和變壓2種處理.恒壓處理下，在下邊界處由真空泵不斷抽氣，提供1個(gè)恒定的負(fù)壓（-21.4kPa），整個(gè)過程時(shí)長(zhǎng)為76h；在變壓處理下，實(shí)驗(yàn)的總時(shí)長(zhǎng)為60h，其下邊界的負(fù)壓隨時(shí)間變化情況見圖2.

圖2 沙柱下邊界負(fù)壓隨時(shí)間變化Fig.2 Temporalevolutionofthevacuumpressureatthelower boundaryofthesandcolumn

2 模型與原理

2.1 模型描述

實(shí)驗(yàn)設(shè)計(jì)為一維沙柱，因此所選用的模型均為一維模型.控制方程（Richards方程）描述了土壤中水分的運(yùn)動(dòng)過程，形式如下:

式中:θ為容積含水量，cm3/cm3；h為壓力水頭，cm；K（h）為非飽和導(dǎo)水率，cm/min；z表示垂向深度，cm；t為時(shí)間，min；S是源匯項(xiàng)，1/min.Richards方程描述土壤中含水量與水頭隨時(shí)間的變化關(guān)系.為求解該方程，還需要知道含水量與水頭之間的關(guān)系，以及非飽和導(dǎo)水率與水頭之間的關(guān)系，即VAN GENUCHTEN-Mualem（VGM）模型.

m=1-1/n.

式中:α為和平均粒徑大小有關(guān)的參數(shù)；n為和粒徑均勻性有關(guān)的參數(shù)；Se為有效飽和度［14］.在給定的邊界條件和初始條件下，HYDRUS利用有限元法數(shù)值求解Richards方程.本文中，需要反演的參數(shù)為θs（飽和含水量）、θr（殘余含水量）、α、n、Ks（飽和導(dǎo)水率）.

2.2 反演方法

在基于優(yōu)化的反演算法中，首要一步是定義一個(gè)目標(biāo)函數(shù)Φ，可表示為:

式中:Φ為目標(biāo)函數(shù)；β為待估計(jì)參數(shù)；h為模型預(yù)測(cè)值；w為權(quán)重因子；h*為實(shí)測(cè)值；z為垂向深度坐標(biāo)；t為時(shí)間.當(dāng)選取合適的權(quán)重后，利用Levenberg-Marquardt非線性優(yōu)化算法調(diào)整參數(shù)，使得實(shí)測(cè)值與模型預(yù)測(cè)值的差值最小，可獲得目標(biāo)函數(shù)的最小值，此時(shí)對(duì)應(yīng)的參數(shù)即為反演結(jié)果.該方法只能給出1個(gè)結(jié)果，無法量化內(nèi)在的不確定性.在貝葉斯框架下，所有待定的模型與參數(shù)都被視為隨機(jī)變量，都用1個(gè)相應(yīng)的分布函數(shù)（也稱為概率密度函數(shù)）來表征.當(dāng)貝葉斯原理用于反演問題時(shí)，觀測(cè)數(shù)據(jù)對(duì)模型與參數(shù)的影響用后驗(yàn)分布來表示.在獲得觀測(cè)數(shù)據(jù)后，關(guān)于模型參數(shù)的統(tǒng)計(jì)信息，如常用的均值、方差、分位值、置信區(qū)間等都可以從后驗(yàn)分布得到.貝葉斯原理可表示為

式中:m為待估計(jì)的非飽和水力參數(shù)，d為觀測(cè)到的水頭數(shù)據(jù)，P（m）為獲得觀測(cè)數(shù)據(jù)前的參數(shù)分布函數(shù)，即先驗(yàn)分布（prior），反映了我們?cè)诘玫接^測(cè)值之前對(duì)目標(biāo)參數(shù)的認(rèn)識(shí)，可通過查詢文獻(xiàn)、地質(zhì)調(diào)查等手段獲取，具有一定的主觀性.在信息不夠齊全的時(shí)候，可假設(shè)為在某個(gè)區(qū)間內(nèi)的均勻分布.P（d|m）稱為似然比（likelihood），表示在參數(shù)取m時(shí)觀測(cè)值d的分布函數(shù)，m與d由模型（這里即Richards方程）聯(lián)系起來.P（m|d）即為后驗(yàn)分布，反映了從觀測(cè)數(shù)據(jù)中獲取信息后對(duì)目標(biāo)參數(shù)的更新認(rèn)識(shí).后驗(yàn)分布是貝葉斯參數(shù)估計(jì)法里面的核心目標(biāo).由于Richards方程具有高度非線性，即使在未知參數(shù)具有簡(jiǎn)單均勻分布時(shí)，后驗(yàn)分布也無法得到解析解.本文使用最新的MCMC算法，DREAM（ZS），來生成滿足后驗(yàn)分布的參數(shù)樣本.

3 結(jié)果與分析

3.1 數(shù)值模擬

本算例參照沙柱恒壓的實(shí)驗(yàn)條件，從先驗(yàn)范圍（表2）中產(chǎn)生1組參數(shù)作為“真實(shí)”值輸入非飽和水流模型，在觀測(cè)時(shí)刻點(diǎn)產(chǎn)生相應(yīng)的水頭輸出值，通過對(duì)其加入標(biāo)準(zhǔn)差為5 cm的白噪聲擾動(dòng)，得到1組模擬觀測(cè)值.分別用LM優(yōu)化算法和MCMC算法進(jìn)行參數(shù)反演.對(duì)于LM算法，選擇在先驗(yàn)范圍內(nèi)隨機(jī)選取1組參數(shù)作為優(yōu)化搜索的初始值；對(duì)于MCMC算法，設(shè)置最長(zhǎng)鏈條長(zhǎng)度為50 000，選取后面5 000長(zhǎng)度的鏈條進(jìn)行后驗(yàn)分析.因此，以計(jì)算量衡量，LM約需要70次HYDRUS求解器調(diào)用，計(jì)算時(shí)間約30 s；MCMC需要50 000次求解器調(diào)用，計(jì)算時(shí)間約4 h.

表2 各層沙子非飽和水力參數(shù)先驗(yàn)范圍Table2 Priori ranges of the unsaturated hydraulic parameters for different sand layers

圖3 各層沙體非飽和水力參數(shù)反演結(jié)果Fig.3!Inversionresultsofhydraulicparametersfordifferentsandlayers

MCMC算法和LM優(yōu)化法反演參數(shù)的結(jié)果見圖3.紅色的實(shí)線表示各水力參數(shù)的真實(shí)值，藍(lán)色的線表示MCMC得到后驗(yàn)概率分布，紅色的虛線表示用LM優(yōu)化算法得到的結(jié)果，是一個(gè)確定的值. LM優(yōu)化法得到的反演參數(shù)與真實(shí)值偏差不穩(wěn)定，同時(shí)容易限于局部最優(yōu).MCMC也同樣可以給出具體的估計(jì)值.例如，若選用后驗(yàn)概率分布的最大值作為參數(shù)估計(jì)值（即最大后驗(yàn)估計(jì)，Maximum-a-Posterior MAP估計(jì)），從圖3中可見，大部分的MAP估計(jì)值也較靠近真實(shí)參數(shù)值.更重要的是，MCMC可以充分量化參數(shù)的不確定性，即使最大后驗(yàn)估計(jì)值偏離了真實(shí)值，后驗(yàn)分布的主要分布區(qū)間（也即置信區(qū)間）基本上都能將真實(shí)值涵蓋其中，因此可以避免只使用單一反演結(jié)果進(jìn)行預(yù)測(cè)帶來的風(fēng)險(xiǎn).

在實(shí)際問題中，由于真實(shí)參數(shù)都是未知的.另一種判定反演結(jié)果準(zhǔn)確性的方法是查看使用反演參數(shù)值得到水頭的預(yù)測(cè)值和實(shí)測(cè)值之間的擬合程度.圖4表示不同觀測(cè)位置處不同反演方法得到的壓力水頭的預(yù)測(cè)值隨時(shí)間的變化關(guān)系.紅色表示壓力水頭真實(shí)值，藍(lán)線表示LM優(yōu)化法對(duì)應(yīng)的預(yù)測(cè)結(jié)果，而灰色的表示用MCMC反演得到的水頭預(yù)測(cè)值80%的置信區(qū)間.可以看出MCMC的反演結(jié)果更能將真實(shí)值覆蓋住，而LM優(yōu)化法的結(jié)果較MCMC的反演結(jié)果，與實(shí)測(cè)值偏差較大.因此，MCMC的效果要優(yōu)于優(yōu)化算法.

3.2 沙柱試驗(yàn)

在沙柱實(shí)驗(yàn)中，由于真實(shí)參數(shù)值未知，因此無法從參數(shù)的后驗(yàn)分布中判定MCMC方法與LM優(yōu)化法的優(yōu)劣.利用下邊界恒壓處理獲取的水頭觀測(cè)數(shù)據(jù)與用MCMC算法對(duì)沙柱各層沙的水力參數(shù)反演結(jié)果見圖5.可以看出，有些參數(shù)的后驗(yàn)分布呈現(xiàn)出明顯的多峰分布，意味著不止一個(gè)最優(yōu)解.因此用MCMC的方法可以有效規(guī)避得到一個(gè)局部最優(yōu)解的問題.

用MCMC反演出來的參數(shù)樣本帶入模型，可以得到沙柱恒壓實(shí)驗(yàn)中壓力水頭的模擬值，即圖6中的灰色區(qū)域；紅線表示沙柱恒壓實(shí)驗(yàn)中壓力水頭的觀測(cè)值.二者擬合效果很好，可以認(rèn)為MCMC反演方法比較可靠.

為進(jìn)一步驗(yàn)證MCMC反演得到的參數(shù)在系統(tǒng)外界條件變化時(shí)的預(yù)測(cè)能力，我們利用下邊界變壓處理下的壓力水頭觀測(cè)值作為校驗(yàn)，將其與恒壓處理下反演參數(shù)結(jié)果對(duì)應(yīng)的預(yù)測(cè)值進(jìn)行對(duì)比.如圖7，紅線為壓力水頭觀測(cè)值，灰色區(qū)域表示壓力水頭模擬值80%的置信區(qū)間.除10 cm觀測(cè)點(diǎn)處二者擬合效果不是很好，其他觀測(cè)點(diǎn)處擬合性很高.

圖4 下邊界恒定負(fù)壓下，沙柱各觀測(cè)點(diǎn)壓力水頭預(yù)測(cè)值隨時(shí)間的變化（數(shù)值算例）Fig.4 Temporal evolution of predicted pressure head at measurement locations（numerical case）with the constant vacuum pressure at bottom boundary

圖5 各參數(shù)的后驗(yàn)分布Fig.5!Posteriordistributionsoftheparameters

圖6 下邊界恒定負(fù)壓下，沙柱各觀測(cè)點(diǎn)壓力水頭預(yù)測(cè)值隨時(shí)間的變化（沙柱試驗(yàn)）Fig.6 Temporal evolution of predicted pressure head at measurement locations（sand experiment）with the constant vacuum pressure at bottom boundary

圖7 下邊界變化負(fù)壓下，沙柱各觀測(cè)點(diǎn)壓力水頭預(yù)測(cè)值隨時(shí)間的變化（沙柱試驗(yàn)）Fig.7 Temporal evolution of predicted pressure head at measurement locations（sand experiment）with the changing vacuum pressure at bottom boundary

4 結(jié)論

4.1 對(duì)于本研究的多層異質(zhì)下非飽和水流的反演問題，LM優(yōu)化方法和MCMC方法都能給出較為合理的結(jié)果.

4.2 LM優(yōu)化方法使用廣泛，求解速度較快，但受參數(shù)初始值影響較大，預(yù)測(cè)結(jié)果與觀測(cè)值存在一定的偏差，同時(shí)由于只能給出一個(gè)單一的反演結(jié)果，無法量化結(jié)果的不確定性.

4.3 與基于優(yōu)化的反演方法相比，MCMC反演方法不僅能夠更好地得出參數(shù)的單一估計(jì)值，同時(shí)預(yù)測(cè)結(jié)果與觀測(cè)值也具有更好的一致性；更重要的是可以給出未知參數(shù)的后驗(yàn)分布，從而準(zhǔn)確量化非飽和水力參數(shù)的不確定性，避免了基于單一參數(shù)反演結(jié)果進(jìn)行預(yù)測(cè)的風(fēng)險(xiǎn).但是，相對(duì)于LM優(yōu)化方法，MCMC計(jì)算量大大增加.

（References）:

［1］ RICH ARDS L A.Capillary conduction of liquid sin porous mediums.Physics，1931（1）:30-33.

［2］ RADCLIFFE D E，SIMUNEK J.Soil Physics with HYDRUS:Modeling and Applications.Boca Raton，F(xiàn)L: CRC Press，2010.

［3］ SIMUNEK J，HUANG K，VAN GENUCHTEN M T.The HYDRUS-ET Software Package for Simulating the One-Dimentional Movement of Water，Heat and Multiple Solutes in Variably-Saturated Media，Version 1.1. Bratislava:Inst.Hydrology Slovak Acad.Sci，1997.

［4］ SIMUNEK J，SEJNA M，VAN GENUCHTEN M T. HYDRUS-2D:Simulating Water Flow and Solute Transport in Two-dimensional Variably Saturated Media. Golden，Colorado:International Groundwater Modeling Center，Colorado School of Mines，1996.

［5］ RUSSO D，BRESLER E.Soil hydraulic properties as stochastic processes:Ⅰ.an analysis of field spatial variability.Soil Science Society of America Journal.1981，45（4）:682-687.

［6］ EVENSEN G.The ensemble Kalman filter:theoretical formulation and practical implementation.Ocean Dynamics，2003，53（4）:343-367.

［7］ GILKS W R.Markov Chain Monte Carlo.Wiley Online Library，2005.

［8］ BERNARDO J E M，SMITH A F.Bayesian Theory.John Wiley&Sons，2009.

［9］ ZHOU H，GOMEZ-HERNANDEZ J J，LI L.Inverse methods in hydrogeology:evolution and recent trends. Advances in Water Resources，2014，63:22-37.

［10］ GAMERMAN D，LOPES H F.Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference.CRC Press，2006.

［11］ HAARIO H，LAINE M，MIRA A，et al.DRAM:efficient adaptive MCMC.Statistics and Computing，2006，16（4）: 339-354.

［12］ VRUGT J A，TER BRAAK C，DIKS C，et al.Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling. International Journal of Nonlinear Sciences and Numerical Simulation，2009，10（3）:273-290.

［13］ SCHARNAGL B，VRUGT J A，VEREECKEN H，et al. Inverse modelling of in situ soil water dynamics:investigating the effect of different prior distributions of the soil hydraulic parameters.Hy drology and Earth System Sciences，2011，15（10）:3043-3059.

［14］ VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils.Soil Science Society of America Journal，1980，44（5）:892-898.

Comparative study on inversion of the unsaturated hydraulic parameters using optimization and
Bayesian estimation methods.Journal of Zhejiang University（Agric.&Life Sci.），2016，42（5）:598- 606

KE Fengqiao，MAN Jun，ZENG Lingzao*，WU Laosheng
（Institute of Soil Water Resources and Environment，College of Environmental&Resource Sciences，Zhejiang University，Hangzhou 310058，China）

SummaryThe model of water movement in variably saturated flow is of guiding significance in agricultural production and environmental protection.The accurate acquisition of soil hydraulic parameters is the precondition of reliable prediction.Based on searching for one set of parameters that best fit the measurements，traditional optimization methods can not quantify the uncertainty of parameters.Now，a newly developed Markov Chain Monte Carlo（MCMC）algorithm，i.e.，DREAM（ZS）was adopted for efficient estimation and accurate uncertainty quantification of soil hydraulic parameters.With the measurements obtained from a one dimensional sand column infiltration experiment，the soil hydraulic parameters were estimated by two approaches，i.e.，the MCMC algorithm and the Levenberg-Marquardt（LM）nonlinear optimization algorithm.Then the two approaches were compared in terms of parameter estimation and state prediction.It can be concluded that:

1）The LM algorithm can provide a single set of model parameter estimations with only a few model runs. However，this method is sensitive to the initial guess of parameters，and the obtained predictions occasionally deviate from the measurements.2）The MCMC algorithm can provide state predictions which better fit measurements.More importantly，it accurately quantifies the uncertainty of parameters，which can avoid the potential risk introduced by making predictions via a single estimated value.

Bayesian estimation；parameter inversion；unsaturated flow；uncertainty

S 121

A

10.3785/j.issn.1008-9209.2016.11.091

國(guó)家自然科學(xué)基金（41371237；41271470）.

*通信作者（Corresponding author）:曾令藻（http://orcid.org/0000-0002-4094-1310），E-mail:lingzao@zju.edu.cn

聯(lián)系方式:柯鳳喬（http://orcid.org/0000-0002-6313-8790），E-mail:kfqiao.cool@163.com

（Received）:2015 11 09；接受日期（Accepted）:2016 05 13；

日期（Published online）:2016 09 18 URL:http://www.cnki.net/kcms/detail/33.1247.S.20160918.1533.010.html