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

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

0 引言

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

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

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

1 研究區(qū)域

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

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

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

2 研究方法

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

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

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

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

在EnKF 中，定義第t個時刻的參數(shù)和狀態(tài)向量集合為：

式中：X為Nx維狀態(tài)變量；A為Na維參數(shù)向量；B為Nb維系統(tǒng)狀態(tài)向量；且

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

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

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

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

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

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

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

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

2.2.2 ESMDA

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

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

2.2.3 ESMDA-EnKF

ESMDA-EnKF 在每次循環(huán)中采用EnKF 對參數(shù)進行更新，其更新公式如下：

第i次循環(huán)：

其中，Yall為經(jīng)過αiRe擾動的觀測資料。

ESMDA-EnKF的具體實現(xiàn)如下偽代碼所示：

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

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

2.4 模型評價

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

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

3 結果與討論

3.1 算法參數(shù)分析

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

3.1.1 ESMDA執(zhí)行次數(shù)對不確定性分析的影響

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

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

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

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

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

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

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

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

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

3.2 地下水模型參數(shù)不確定性分析

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

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

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

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

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

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

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

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

表2 ESMDA-EnKF算法第2次迭代后地下水模型參數(shù)統(tǒng)計值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 結論

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

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

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

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