基于數據同化的地下水模型不確定性分析

陳 沖，張 偉，邢慶輝，豆沂宣

（1.中國石油大學（北京）信息科學與工程學院，北京 102249；2.中國科學院 西北生態環境資源研究院 冰凍圈科學國家重點實驗室/可可托海站，甘肅 蘭州 730000）

0 引言

地下水系統模擬已經成為地球科學的基本研究方法，在地球系統科學領域被譽為第二次哥白尼革命［1］，被證明是最有價值和最實用的工具之一［2-3］。地下水模型能夠刻畫水文過程的整體和局部行為；多次模擬，選擇最優的設計；通過情景分析進行情景預測，提出應對策略。圍繞黑河中游地表水、泉水和地下水，大量學者利用模型進行了有益探索：陳沖等［4］建立了黑河中游地區地下水模型，并初步分析了上游徑流與耕地面積對地下水資源的影響；王旭升等［5］總結了20 年來黑河流域的地下水模型研究，并指出地下水模型應加強與關聯過程（地表水、土壤水、水利工程）的集成；程國棟等［6］闡述了黑河流域生態—水文過程集成研究的進展及展望。目前的地下水模型一般基于確定性的參數、邊界條件，采用確定性機制來構建模型。然而，由于實際過程的復雜性及非線性因素，導致模型往往存在一定不確定性。正如George 所言，所有關于真實系統的模型都是“錯誤”的［7］，即模型永遠不能完美表達實際系統。地下水模型的不確定性來源一般可以分為：模型數據的不確定性、模型結構的不確定性以及模型參數的不確定性。模型數據的不確定性一般指由于監測誤差或者數據缺失導致的不確定性，只能通過提高監測技術和數據收集頻率來降低。而模型結構和參數不確定性則來源于建模過程、模型參數標定及驗證過程。這些不確定性將對進一步的科學研究以及決策提出等帶來諸多麻煩甚至風險。因此在使用地下水模型時，對其進行不確定性分析是十分必要的。

Oberkampf等［8］介紹了用來估計、分析模型中所有類型不確定性的多種方法。Beven［9］于1989 年在對物理模型的適用性進行討論的過程中，首次提出了對水文模型進行不確定性分析的概念。地下水模型不確定性分析方法主要分為蒙特卡羅（Monte Carlo）法、矩方程法與貝葉斯法三種［10］。雖然MC法是一種被廣泛采用的不確定性分析方法［11-15］，但是MC 法的缺點也不容忽視，其需對模型的參數進行大量的采樣才能保證算法的正確性和準確性，無法適用于計算消耗大的模型。矩方程法通過隨機偏微分方程直接求解模擬結果的各階統計矩，已在地下水模型不確定性研究中得到初步應用［16-18］。貝葉斯法利用觀測資料修正水文地質參數分布，其既可以用于模型參數識別反演［19］，也能夠用于對參數進行不確定性分析，其優點是通過更新參數分布，以更加準確的評估模型不確定性。利用貝葉斯法進行地下水模型不確定性分析的研究相對較少［20-21］。數據同化（DA，Data Assimilation）是基于貝葉斯理論的參數更新方法?？柭鼮V波算法（KF，Kalman Filter）是由Kalman 提出的順序數據同化算法［22］，其利用觀測資料自回歸地對模型的狀態變量進行更新，并在更新的整個過程中保證狀態變量估計值的誤差最小。自提出以來，KF算法已被廣泛應用于模型參數與狀態的估計中。黃春林等［23］基于集合卡爾曼濾波（EnKF，Ensemble Kalman Filter）利用土壤水分觀測數據同化了簡單生物圈模型（SiB2，Simple Biosphere Model 2），探討了簡單生物圈模型的單點土壤水分同化方案；褚楠等［24］進一步采用雙EnKF方法同時估計SiB2模型中土壤水分與土壤屬性參數，提高了模型對土壤水分的估計精度。Li 等［25］實現了采用EnKF 同時估計地下水模型的參數與狀態變量。Sly 等［26］基于EnKF算法采用SWOT（Surface Water and Ocean Topography）模型的 輸出對GHMs（Global Hydrological Models）進行了同化，提升了全球尺度水文模擬的精度。然而，僅有少量研究基于數據同化方法進行模型不確定性分析。例如，Li 等［25］采用了ESMDA（Ensemble Smoother with Multiple Data Assimilation）分析了焉耆盆地水文模型參數的不確定性。

本文基于前期工作中在黑河流域中游構建的地下水模型［4］，擬結合ESMDA 與EnKF 實現地下水模型參數的不確定性分析，探討ESMDA 與EnKF算法超參數對模型不確定性分析效果的影響，利用算法定量分析地下水模型參數的不確定性、對不同觀測數據的敏感性以及不同參數的不確定性范圍。

1 研究區域

黑河流域是我國第二大內陸河流域，位于我國西北部干旱區，面積約為14×104km2。流域上游的祁連山區是我國重要的冰凍圈分布區，氣溫較低、蒸發弱，年平均降水量大于300 mm，豐富的冰凍圈融水與降水是整個流域的主要水源。流域下游為額濟納旗盆地，為典型的大陸性氣候，具有降水少（年均降水量為43.7 mm）、蒸發強（年均蒸發量為2 248.8 mm）、日照時間長、晝夜溫差大等特性［27］。流域內支流共35 個，由于過度開采及蒸發，絕大部分支流在山前消失，中游主要河流只有黑河干流與梨園河（圖1）。黑河發源于祁連山北麓，由東南流向西北，途徑青海、甘肅及內蒙古三個省區，最終注入居延海，全長約821 km，水量主要由降水、冰雪融水以及地下水出流組成。鶯落峽水文站為黑河出山口徑流量觀測站，多年（1945—2012 年）平均年徑流量約為15.9×108m3，是中游張掖市、臨澤縣、高臺縣及下游金塔東部和額濟納旗綠洲等地城市工業、生活用水的主要水源［28］。黑河流域中游（38°38′～39°53′ N，98°53′～100°44′ E；圖1）為典型的干旱氣候，年降水量稀少（90～160 mm），潛在蒸發量大（2 000～2 500 mm），總面積約9 016 km2［29-30］。河流流出祁連山后，一部分被引入灌溉渠系和供水系統，其余則沿河床下泄，并于沿途滲入地下，補給地下水，黑河在山前沖積扇區損失33%的徑流量［31］。長期以來，黑河中游農業生產與下游濕地、尾閭湖等生態環境之間存在著用水矛盾，且由于社會、經濟發展和人口的增長，用水矛盾日益突出。

圖1 研究區地理位置Fig.1 Study area and observation sites

在前期工作中，我們構建了黑河流域中游的地下水模型，采用了SRTM（Shuttle Radar Topography Mission）的DEM（Digital Elevation Model）數據［32］、土地利用［33-34］、抽水［35］、灌溉［36］、地下水位［37］以及河流徑流量數據。將SRTM 90 m×90 m 的數據處理為1 km×1 km 的網格數據對研究區含水層系統進行空間離散化。時間上，將所獲取數據（1986 年1 月—2010 年12 月）離散為月時間間隔，采用以每個自然月為一個應力期的非穩定流（暫態流，Transient state）模擬。地下水位數據包括42 口觀測井處的觀測水位，河流徑流量包括鶯落峽、高崖及正義峽處觀測的徑流數據。經過校正（1986 年1 月—2008 年12 月）與驗證（2009 年1 月—2010 年12 月），模型對研究區內地下水位及徑流量的擬合程度較好，所構建的地下水模型比較真實地反映了研究區地下水系統的情況［4］。

2 研究方法

本文采用ESMDA 與EnKF 結合的方法ESMDA-EnKF 對研究區地下水模型的參數進行不確定性分析。圖2 給出了使用ESMDA 算法進行不確定性分析的程序流程。在進行不確定性分析之前，首先要對算法進行初始化，并對模型的參數進行設置，根據給定的概率分布（一般假定為正態分布或者對數正態分布）生成模型參數集合。將模型參數預測值分別輸入地下水模型中，計算得到模型輸出值；綜合地下水模型輸出、觀測值、觀測誤差協方差以及參數預測值集合輸入ESMDA-EnKF 算法中，從模型輸出值和觀測值集合計算得到增量場，根據集合中樣本的誤差統計計算出卡爾曼增益，由增量場與卡爾曼增益計算得到預測結果的更新量，將更新量疊加到初始場得到分析值集合；根據預先算法執行次數判定是否達到程序運行結束條件。

圖2 ESMDA-EnKF算法不確定性分析流程圖Fig.2 Diagram of uncertainty analysis using ESMDA-EnKF

2.1 集合卡爾曼濾波（Ensemble Kalman Filter）

EnKF 采用蒙特卡羅方法隨機產生參數集合，對狀態變量進行預測，并根據獲取的觀測信息對狀態變量進行更新，已經有許多文獻對其理論和具體算法做了詳細的論述［38-39］。此處僅簡要回顧一下EnKF的主要的工作原理。

在EnKF 中，定義第t個時刻的參數和狀態向量集合為：

式中：X為Nx維狀態變量；A為Na維參數向量；B為Nb維系統狀態向量；且

則狀態變量在第t+1時刻的預報值為

式中：f代表狀態變量的預測值（forecast）；a代表狀態變量的分析值（analysis）；為第i個狀態變量在t+1 時刻的預測值；Mt（?）為非線性模型算子；為t時刻的分析值；wk為期望為0，方差為Wk的高斯白噪聲。

在分析步驟（analysis step）中，采用公式（4）對觀測數據進行擾動

2.2 多重數據同化集合平滑器（Ensemble Smoother with Multiple Data Assimilation）

對ESMDA 算法主要分為3 個方面進行介紹，首先介紹集合平滑器（ES，Ensemble Smoother）；接下來在ES 算法中引入膨脹系數構成ESMDA 算法框架；最后，介紹ESMDA與EnKF算法的結合，并給出偽代碼的實現。

2.2.1 集合平滑器（Ensemble Smoother，ES）

EnKF 屬于順序數據同化算法，其根據t時刻狀態變量值初始化模型，預測t+1 時刻模型的狀態變量，在t+1 時刻利用觀測資料對狀態變量的預測值進行加權更新，得到當前時刻狀態變量的最優估計值，多用于對狀態變量進行實時標定，以實現實時模型。然而，ES多用于參數反演以及參數的不確定性分析中，采用所有可用觀測資料對參數進行更新。公式如下：

式中：各個符號定義與EnKF中類似。

2.2.2 ESMDA

ES 的本質是將所有時刻的觀測資料輸入同化算法中進行一次參數更新［40］。然而，正如算法名字所示，ESMDA 在ES 的基礎上設定一個算法執行次數，進行多次數據同化，在每個循環中，ESMDA 并不是簡單重復了ES過程，而是在每個循環中對觀測誤差添加了一個膨脹系數αi，并令

式中：Na為循環次數。顯然地，αi有許多種取值方式，然而文獻［41］指出，隨著循環次數的增加而逐漸減小的膨脹系數對于同化效果并無明顯提升。因此，在每次循環中，采用相同的膨脹系數。

2.2.3 ESMDA-EnKF

ESMDA-EnKF 在每次循環中采用EnKF 對參數進行更新，其更新公式如下：

第i次循環：

其中，Yall為經過αiRe擾動的觀測資料。

ESMDA-EnKF的具體實現如下偽代碼所示：

輸入：Na，觀測資料Oall，Xf，Re

輸出：Xa，K，Pa

read（Oall）

fori=1 to Na

Xf’=Xf-mean（Xf）

Yall=Oall+sqrt（αi* Re）

HXf’=HXf-mean（HXf）

PfHT=（Xf’* HXf’）/（N-1）

HPfHT=（HXf’* transpose（HXf’）/（N-1）

K=PfHT/（HPfHT+αi* Re）

Xa=Xf+K *（Yall-HXf）

Pa=Pf-K * PfHT

update（Xf，Pf）；

return（K，Xa，Pa）

2.3 地下水模型

黑河流域中游地下水模型基于MODFLOW（MODular three-dimension finite-difference groundwater FLOW model）構建，MODFLOW 采用有限差分法將時間與空間離散化以解決地下水在三維空間中的流動問題。研究區平面面積約為9 016 km2，采用1 km×1 km 的正方形網格將研究區含水層系統在水平方向上進行空間離散化，離散化之后，研究區在平面上剖分成為132 行×165 列網格（如圖3）。研究區內定義為活動單元，研究區外定義為非活動單元。時間上，以每月為一個應力期，每天為一個時間步。研究區邊界條件參考文獻［42］以及自然邊界確定（圖3）。由于地下水分水嶺的存在，A～E 設置為無水流邊界；E 處為正義峽水文站，黑河從此處流出研究區；據調查資料分析，北山地區地下水含水介質與南部祁連山相似，但由于降水稀少，無常年地表徑流，地下水含量無法與祁連山區相比較［27］，因此，將D～E 邊界假設為無水流邊界；由于多種側向流從南部祁連山流入研究區，導致A～D邊界最為復雜，因此，將此邊界分為3 段進行定義：如圖3 所示，A～B 與C～D 定義為定流量邊界條件（第二類邊界條件或Neumam 條件），B～C 之間定義為無流量邊界。垂直方向上，模型頂部邊界為空氣—土壤接觸面邊界，模型底板定義為不透水邊界。黑河干支流采用MODFLOW 內的河流（STR，STreamflow-Routing）程序包［43］進行模擬。農業灌溉通過采用地下水補給程序包（RCH，Re-CHarge）［44］在網格上添加補給率實現。蒸散發采用了MODFLOW 內置的蒸散發程序包（EVT，EVapo-Transpiration）進行模擬。機井采用機井程序包（Well）進行模擬。

圖3 黑河中游模型概化及邊界設置Fig.3 The conceptualization and boundary conditions of the groundwater model

2.4 模型評價

本文采用均方根誤差（RMSE，Root Mean Square Error）對模擬結果進行量化評價，方程如下：

式中：N為觀測值總數（觀測井個數×應力期個數）；Mi和Oi分別是模型模擬值和觀測值。

3 結果與討論

3.1 算法參數分析

在使用ESMDA-EnKF 算法對地下水模型進行不確定性分析的過程中，算法存在一些超參數直接影響到同化算法對模型參數的更新效果（例如：算法執行次數（即膨脹系數）Na、觀測資料數n、參數集合大小N等），因此，本文首先為ESMDA-EnKF確定最優超參數，之后采用算法分析模型參數的不確定性以及模型參數對觀測數據的敏感性。

3.1.1 ESMDA執行次數對不確定性分析的影響

ESMDA 采用EnKF 算法對地下水模型參數進行更新，因此，EnKF 算法的執行次數將直接影響模型參數的更新效果。然而，目前尚沒有關于執行次數的理論研究，而在文獻［41］給出的例子中，分別執行了2次、4次算法以對比同化效果?？紤]到計算消耗問題，本節評價了執行1 次、2 次、3 次、4 次En-KF算法后的同化效果（圖4）。

圖4 不同執行次數對ESMDA-EnKF算法效果的影響Fig.4 Results of ESMDA-EnKF for different iterations

圖4 中展示了EnKF 算法執行不同次數后的同化效果，Last iteration 代表采用不同的執行次數時，最后一次執行結束之后的同化效果。由圖中可以看出隨著EnKF 算法對模型參數的更新，地下水位的RMSE 值逐漸降低，表明地下水位模擬值逐漸接近地下水位的觀測值。經過隨機采樣之后的第一次、第二次參數更新對RMSE 影響較大；第三次及第四次參數更新之后，與第二次參數更新之后的RMSE值相比，變化不大。這表明，使用觀測數據對地下水模型參數進行的前2 次更新最為有效，可能是因為：（1）地下水模型的參數采用隨機初始化，因此，ESMDA 在隨機模型參數的基礎上進行的第一次更新力度相對較大；（2）ESMDA 采用了所有時刻的觀測數據進行更新，最大程度上利用了觀測數據中的信息，因此，能夠更加有效且大幅度的更新參數。

3.1.2 觀測資料數對不確定性分析的影響

數據同化的目的即是最大限度的融合不同來源、不同時間、空間分辨率的直接、間接觀測數據以用于提高模型的估計精度，以獲得更加準確的狀態變量的時空分布［45-46］。然而由于成本原因，對系統內相關因素進行無限制觀測是不現實的。因此，研究觀測資料對同化效果的影響是必要的。選擇42個水位觀測井處的觀測水位以及高崖和正義峽水文觀測站處的黑河流量觀測數據，分別研究觀測數為7、14、21、28、35、42 以及44 時對同化效果的影響。圖5 顯示了經過2 次參數更新過程之后，不同觀測資料數的RMSE 值。由于涉及到不同類型的觀測值、輸出值（地下水位、流量），所以，圖中對地下水水位、流量的觀測值及輸出值進行了歸一化處理，將觀測值、輸出值限定于［0，1］之間。由圖5 的總體趨勢上可以看出，隨著觀測資料的增加，參數逐漸接近真實值，地下水位的RMSE 值逐漸降低。然而，觀測數從7 個增加至14 個對同化效果影響不明顯；觀測井數增加至21 個之后，同化效果基本保持不變；增加流量觀測（42 至44）對同化效果影響較為顯著。由此可見，一種類型的觀測數據確實會提升同化效果，但是當一種類型的觀測數據增加到一定程度時（本實驗中為21 個），觀測數據對同化效果的影響有限；然而，繼續增加不同類型的觀測數據（河流流量觀測數據）將進一步提升同化效果。

圖5 不同的觀測數對ESMDA-EnKF算法效果的影響Fig.5 Results of ESMDA-EnKF for different number of observations

3.1.3 集合大小對不確定性分析的影響

EnKF 采用樣本集合的方式表示模型狀態變量的先驗概率分布，估計誤差協方差［47］。Hamill 等［48］研究發現，EnKF 中的背景誤差協方差估計和濾波函數的最優相關尺度等均與樣本集合大小相關；當樣本集合比較小時，協方差的估計噪音較大、最優相關尺度較小，出現濾波發散現象，且集合大小直接影響算法運行時間。因此，采用不同的集合大小運行算法，以探索集合大小對算法同化效果的影響。圖6 顯示了模型參數集合大小分別為30、60、90、120、150 情況下，經過兩次參數更新之后的RMSE 值。由圖中可以看出，當參數集合大小為30時，并沒有出現濾波發散現象，且增加采樣對同化效果的影響并不明顯。因此，本文將集合大小設置為30以平衡算法效果與計算消耗。

圖6 不同集合大小對ESMDA算法效果的影響Fig.6 Results of ESMDA for different ensemble size

3.2 地下水模型參數不確定性分析

基于以上分析確定了算法參數之后，本文利用ESMDA-EnKF 對地下水模型參數進行了不確定性分析，模型參數分布見圖7，相關設置見表1。從模型參數先驗概率分布中進行30 次隨機采樣，使用42 個水位觀測井、2 個徑流觀測站處的觀測數據對模型參數進行2 次（Na=2）更新。模型參數隨ESMDA-EnKF 運行所得的不確定性如圖8（注：由于參數與觀測點數較多，圖8 中只顯示了部分參數與觀測點處的更新情況）。由圖8 可以看出，無論參數的初始分布如何，每次執行ESMDA-EnKF 算法都會對模型參數進行更新，使參數值更接近其標定值。第一次執行算法對模型參數的更新效果顯著，大多模型參數在進行第一次更新之后，參數值已經較為接近參數標定值，不確定性明顯減??；在進行第二次更新之后，模型參數均分布于標定值附近，模型參數的不確定性進一步減小。圖8 第三列顯示了各個參數在ESMDA-EnKF 算法執行完成之后的分布情況，由此分布情況可以看出，各個分區的滲透系數（P1～P8）的不確定性顯著減小，基本收斂于最優值附近。對比各個分區的滲透系數可以看出，經過ESMDA-EnKF 算法更新之后，第Ⅷ個分區的滲透系數（P8）的分布最分散，不確定性最大，這可能是由于在此分區內觀測數據較少，導致算法無法對參數進行有效更新。在算法執行完成之后，黑河子分區Ⅰ與子分區Ⅱ的河床水力傳導系數（P9與P10）雖然也收斂到了最優值附近，但是其分布比黑河子分區Ⅲ的河床水力傳導系數（P11）分散（即其標準差比較大）。這可能是由于在黑河子分區Ⅰ與子分區Ⅱ河段，黑河與地下水相互作用頻繁（子分區Ⅰ河段河流補給地下水，子分區Ⅱ河段地下水出流轉化為河流），導致河床水力傳導系數不確定性較大。由圖8 可見，ESMDA-EnKF 算法能夠對模型參數進行有效更新，從而減小模型參數的不確定性。

圖7 地下水模型參數分區Fig.7 Zonation of the parameters：hydraulic conductivities of the aquifer（a）；hydraulic conductance of the streambed（b）

圖8 經過ESMDA-EnKF更新之后的模型參數分布（參數編號解釋見表1）Fig.8 Parameter distributions after the updating by ESMDA-EnKF（The numbers of parameters are shown in Table 1）

表1 進行不確定性分析的參數及其設置Table 1 Parameters involved in uncertainty analysis

地下水位隨ESMDA-EnKF 運行所得的不確定性如圖9 所示。由圖9 第一列可以看出，受模型參數不確定性的影響，研究區內地下水位［圖9（a）、（b）］以及河流徑流量［圖9（c）、（d）］均存在較大的不確定性，且隨著時間變化，不確定性區間有變大的趨勢，這可能是受預報前期模型運行結果不確定性的累積影響，模型運行后期地下水位與河流徑流量運行結果產生了更大的不確定性。圖9 第二、三列分別為ESMDA-EnKF 算法運行1、2 次之后的模型輸出值?？梢钥闯?，ESMDA-EnKF算法運行一次之后，模型輸出的地下水位與河流徑流量的不確定性范圍明顯縮小，無論是數據起伏較小的地下水位還是數據起伏較大的河流徑流量，都收斂到了最優值附近。第二次執行ESMDA-EnKF 算法則基本消除了模型輸出結果的不確定性（對應于圖9 的第三列）。

圖9 使用ESMDA-EnKF對模型參數進行更新之后的模型輸出值Fig.9 Model outputs after the update of model parameters by ESMDA-EnKF

對比圖8～9 可以看出，隨著ESMDA-EnKF 算法的執行，圖8中的模型參數逐漸收斂于參數最優值，相應地，通過模型參數計算得到的模型輸出也逐漸接近觀測值（圖9）。由算法對參數的更新可以看出，對觀測數據最為敏感的模型參數為第二個參數分區以及第六個參數分區的滲透系數P2 和P6。最不敏感的參數為第八個參數分區的滲透系數，這可能是因為此區域緊鄰不透水邊界，且區域內沒有河流，其與其他區域的水力聯系較弱，而在此分區內缺乏有效的觀測數據，無法對模型參數進行有效更新。圖9顯示了在不同算法運行階段的模型輸出值變化，由灰色逐漸加深為黑色，分別為不同迭代次數的模型輸出值，由圖中可以看出含水層滲透系數（P1～P8）、含水層給水度（P12）以及灌溉回流系數（P13）對模型地下水位影響顯著。由圖8～9 中河流流量與河床水力傳導系數的同化效果可以看出，水力傳導系數對河流流量觀測較為敏感。河床傳導系數會影響河流與地下水之間的水量交換，然而，河流流量還受模型模擬區域內整體地下水水位以及上游來水量的影響，所以，僅使用高崖和正義峽水文站處流量觀測值對河床水力傳導系數進行更新，無法完全消除其不確定性。

為了定量的分析地下水模型參數的不確定性，同時進一步說明ESMDA-EnKF 算法對地下水模型參數及輸出不確定性的更新效果，本文運用Origin 2018 軟件對ESMDA-EnKF 第2 次迭代后的地下水模型參數以及第2 次迭代后地下水模型最后時刻（2008年12月）的輸出結果（地下水位、河流流量）進行了統計分析。地下水模型參數的統計指標見表2，地下水模型輸出結果的統計指標見表3。由表2中可以看出，經過ESMDA-EnKF 算法更新后，地下水模型參數的標準差比更新之前明顯變小，且參數均值與表1中參數均值基本相同。雖然更新之后的參數均值與參數的標定值相似，但是仍有細微差別，這可能是由于算法中僅按照相應的分布隨機采樣了30 個樣本，樣本數較小導致抽樣誤差（sampling error）較大，但是根據中心極限定理（central limit theorem）可知，采樣之后的均值服從正態分布，因此，均值誤差在可接受的范圍之內。

表2 ESMDA-EnKF算法第2次迭代后地下水模型參數統計值Table 2 Statistical values of groundwater model parameters after two iterations of ESMDA-EnKF

表3 ESMDA-EnKF算法第2次迭代后地下水模型最后時刻的輸出結果Table 3 Statistical values of groundwater model outputs in the last time step after two iterations of ESMDA-EnKF

4 結論

本文提出了基于數據同化的不確定性分析方法，通過包含觀測資料信息減小模型參數的不確定性?；谝褬嫿ǖ暮诤恿饔蛑杏蔚叵滤Ｐ?，利用提出的算法對模型參數進行不確定性分析，探索了不同超參數對算法效果的影響，分析了含水層滲透系數、河床水力傳導系數、含水層給水度、灌溉回流系數的不確定性，評估了不同系數對觀測信息的敏感程度，主要結論如下：

（1）基于貝葉斯理論的數據同化方法是不確定性分析的有效手段，ESMDA-EnKF算法通過包含觀測資料更新參數，能夠有效的減小模型參數的不確定性；對不確定性分析結果的統計分析結果表明，經過ESMDA-EnKF 算法更新后的參數均值收斂到最優均值附近，不確定性范圍明顯減小。

（2）不同超參數對算法效果影響不一，觀測資料的增加將提升算法對模型參數的更新程度，減小參數的不確定性，且加入新類型的觀測資料會進一步減小參數的不確定性。

（3）不同模型參數的不確定性不同，地表水與地下水相互作用頻繁的區域參數不確定性較大；含水層滲透系數、含水層給水度以及灌溉回流系數對模型輸出的地下水位影響顯著，水力傳導系數對模型輸出的河流流量影響較大。