基于數(shù)據(jù)同化的地下水模型不確定性分析

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

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

0 引言

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

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

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

1 研究區(qū)域

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

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

在前期工作中，我們構(gòu)建了黑河流域中游的地下水模型，采用了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)進(jìn)行空間離散化。時(shí)間上，將所獲取數(shù)據(jù)（1986 年1 月—2010 年12 月）離散為月時(shí)間間隔，采用以每個(gè)自然月為一個(gè)應(yīng)力期的非穩(wěn)定流（暫態(tài)流，Transient state）模擬。地下水位數(shù)據(jù)包括42 口觀測井處的觀測水位，河流徑流量包括鶯落峽、高崖及正義峽處觀測的徑流數(shù)據(jù)。經(jīng)過校正（1986 年1 月—2008 年12 月）與驗(yàn)證（2009 年1 月—2010 年12 月），模型對研究區(qū)內(nèi)地下水位及徑流量的擬合程度較好，所構(gòu)建的地下水模型比較真實(shí)地反映了研究區(qū)地下水系統(tǒng)的情況［4］。

2 研究方法

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

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

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

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

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

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

則狀態(tài)變量在第t+1時(shí)刻的預(yù)報(bào)值為

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

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

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

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

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

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

式中：各個(gè)符號定義與EnKF中類似。

2.2.2 ESMDA

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

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

2.2.3 ESMDA-EnKF

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

第i次循環(huán)：

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

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

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

2.4 模型評價(jià)

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

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

3 結(jié)果與討論

3.1 算法參數(shù)分析

在使用ESMDA-EnKF 算法對地下水模型進(jìn)行不確定性分析的過程中，算法存在一些超參數(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ù)進(jìn)行更新，因此，EnKF 算法的執(zhí)行次數(shù)將直接影響模型參數(shù)的更新效果。然而，目前尚沒有關(guān)于執(zhí)行次數(shù)的理論研究，而在文獻(xiàn)［41］給出的例子中，分別執(zhí)行了2次、4次算法以對比同化效果?？紤]到計(jì)算消耗問題，本節(jié)評價(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ù)時(shí)，最后一次執(zhí)行結(jié)束之后的同化效果。由圖中可以看出隨著EnKF 算法對模型參數(shù)的更新，地下水位的RMSE 值逐漸降低，表明地下水位模擬值逐漸接近地下水位的觀測值。經(jīng)過隨機(jī)采樣之后的第一次、第二次參數(shù)更新對RMSE 影響較大；第三次及第四次參數(shù)更新之后，與第二次參數(shù)更新之后的RMSE值相比，變化不大。這表明，使用觀測數(shù)據(jù)對地下水模型參數(shù)進(jìn)行的前2 次更新最為有效，可能是因?yàn)椋海?）地下水模型的參數(shù)采用隨機(jī)初始化，因此，ESMDA 在隨機(jī)模型參數(shù)的基礎(chǔ)上進(jìn)行的第一次更新力度相對較大；（2）ESMDA 采用了所有時(shí)刻的觀測數(shù)據(jù)進(jìn)行更新，最大程度上利用了觀測數(shù)據(jù)中的信息，因此，能夠更加有效且大幅度的更新參數(shù)。

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

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

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

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

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

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

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

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

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

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

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

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

表2 ESMDA-EnKF算法第2次迭代后地下水模型參數(shù)統(tǒng)計(jì)值Table 2 Statistical values of groundwater model parameters after two iterations of ESMDA-EnKF

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

4 結(jié)論

本文提出了基于數(shù)據(jù)同化的不確定性分析方法，通過包含觀測資料信息減小模型參數(shù)的不確定性?；谝褬?gòu)建的黑河流域中游地下水模型，利用提出的算法對模型參數(shù)進(jìn)行不確定性分析，探索了不同超參數(shù)對算法效果的影響，分析了含水層滲透系數(shù)、河床水力傳導(dǎo)系數(shù)、含水層給水度、灌溉回流系數(shù)的不確定性，評估了不同系數(shù)對觀測信息的敏感程度，主要結(jié)論如下：

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

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

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