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

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

0 引言

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

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

1 研究區(qū)域

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

圖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ù)對(duì)研究區(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 口觀測(cè)井處的觀測(cè)水位，河流徑流量包括鶯落峽、高崖及正義峽處觀測(cè)的徑流數(shù)據(jù)。經(jīng)過(guò)校正（1986 年1 月—2008 年12 月）與驗(yàn)證（2009 年1 月—2010 年12 月），模型對(duì)研究區(qū)內(nèi)地下水位及徑流量的擬合程度較好，所構(gòu)建的地下水模型比較真實(shí)地反映了研究區(qū)地下水系統(tǒng)的情況［4］。

2 研究方法

本文采用ESMDA 與EnKF 結(jié)合的方法ESMDA-EnKF 對(duì)研究區(qū)地下水模型的參數(shù)進(jìn)行不確定性分析。圖2 給出了使用ESMDA 算法進(jìn)行不確定性分析的程序流程。在進(jìn)行不確定性分析之前，首先要對(duì)算法進(jìn)行初始化，并對(duì)模型的參數(shù)進(jìn)行設(shè)置，根據(jù)給定的概率分布（一般假定為正態(tài)分布或者對(duì)數(shù)正態(tài)分布）生成模型參數(shù)集合。將模型參數(shù)預(yù)測(cè)值分別輸入地下水模型中，計(jì)算得到模型輸出值；綜合地下水模型輸出、觀測(cè)值、觀測(cè)誤差協(xié)方差以及參數(shù)預(yù)測(cè)值集合輸入ESMDA-EnKF 算法中，從模型輸出值和觀測(cè)值集合計(jì)算得到增量場(chǎng)，根據(jù)集合中樣本的誤差統(tǒng)計(jì)計(jì)算出卡爾曼增益，由增量場(chǎng)與卡爾曼增益計(jì)算得到預(yù)測(cè)結(jié)果的更新量，將更新量疊加到初始場(chǎng)得到分析值集合；根據(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ù)集合，對(duì)狀態(tài)變量進(jìn)行預(yù)測(cè)，并根據(jù)獲取的觀測(cè)信息對(duì)狀態(tài)變量進(jìn)行更新，已經(jīng)有許多文獻(xiàn)對(duì)其理論和具體算法做了詳細(xì)的論述［38-39］。此處僅簡(jiǎn)要回顧一下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ù)測(cè)值（forecast）；a代表狀態(tài)變量的分析值（analysis）；為第i個(gè)狀態(tài)變量在t+1 時(shí)刻的預(yù)測(cè)值；Mt（?）為非線性模型算子；為t時(shí)刻的分析值；wk為期望為0，方差為Wk的高斯白噪聲。

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

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

對(duì)ESMDA 算法主要分為3 個(gè)方面進(jìn)行介紹，首先介紹集合平滑器（ES，Ensemble Smoother）；接下來(lái)在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ù)測(cè)t+1 時(shí)刻模型的狀態(tài)變量，在t+1 時(shí)刻利用觀測(cè)資料對(duì)狀態(tài)變量的預(yù)測(cè)值進(jìn)行加權(quán)更新，得到當(dāng)前時(shí)刻狀態(tài)變量的最優(yōu)估計(jì)值，多用于對(duì)狀態(tài)變量進(jìn)行實(shí)時(shí)標(biāo)定，以實(shí)現(xiàn)實(shí)時(shí)模型。然而，ES多用于參數(shù)反演以及參數(shù)的不確定性分析中，采用所有可用觀測(cè)資料對(duì)參數(shù)進(jìn)行更新。公式如下：

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

2.2.2 ESMDA

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

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

2.2.3 ESMDA-EnKF

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

第i次循環(huán)：

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

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

輸入：Na，觀測(cè)資料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)問(wèn)題。研究區(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è)置為無(wú)水流邊界；E 處為正義峽水文站，黑河從此處流出研究區(qū)；據(jù)調(diào)查資料分析，北山地區(qū)地下水含水介質(zhì)與南部祁連山相似，但由于降水稀少，無(wú)常年地表徑流，地下水含量無(wú)法與祁連山區(qū)相比較［27］，因此，將D～E 邊界假設(shè)為無(wú)水流邊界；由于多種側(cè)向流從南部祁連山流入研究區(qū)，導(dǎo)致A～D邊界最為復(fù)雜，因此，將此邊界分為3 段進(jìn)行定義：如圖3 所示，A～B 與C～D 定義為定流量邊界條件（第二類邊界條件或Neumam 條件），B～C 之間定義為無(wú)流量邊界。垂直方向上，模型頂部邊界為空氣—土壤接觸面邊界，模型底板定義為不透水邊界。黑河干支流采用MODFLOW 內(nèi)的河流（STR，STreamflow-Routing）程序包［43］進(jìn)行模擬。農(nóng)業(yè)灌溉通過(guò)采用地下水補(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 模型評(píng)價(jià)

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

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

3 結(jié)果與討論

3.1 算法參數(shù)分析

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

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

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

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

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

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

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

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

3.1.3 集合大小對(duì)不確定性分析的影響

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)行算法，以探索集合大小對(duì)算法同化效果的影響。圖6 顯示了模型參數(shù)集合大小分別為30、60、90、120、150 情況下，經(jīng)過(guò)兩次參數(shù)更新之后的RMSE 值。由圖中可以看出，當(dāng)參數(shù)集合大小為30時(shí)，并沒有出現(xiàn)濾波發(fā)散現(xiàn)象，且增加采樣對(duì)同化效果的影響并不明顯。因此，本文將集合大小設(shè)置為30以平衡算法效果與計(jì)算消耗。

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

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

基于以上分析確定了算法參數(shù)之后，本文利用ESMDA-EnKF 對(duì)地下水模型參數(shù)進(jìn)行了不確定性分析，模型參數(shù)分布見圖7，相關(guān)設(shè)置見表1。從模型參數(shù)先驗(yàn)概率分布中進(jìn)行30 次隨機(jī)采樣，使用42 個(gè)水位觀測(cè)井、2 個(gè)徑流觀測(cè)站處的觀測(cè)數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行2 次（Na=2）更新。模型參數(shù)隨ESMDA-EnKF 運(yùn)行所得的不確定性如圖8（注：由于參數(shù)與觀測(cè)點(diǎn)數(shù)較多，圖8 中只顯示了部分參數(shù)與觀測(cè)點(diǎn)處的更新情況）。由圖8 可以看出，無(wú)論參數(shù)的初始分布如何，每次執(zhí)行ESMDA-EnKF 算法都會(huì)對(duì)模型參數(shù)進(jìn)行更新，使參數(shù)值更接近其標(biāo)定值。第一次執(zhí)行算法對(duì)模型參數(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)值附近。對(duì)比各個(gè)分區(qū)的滲透系數(shù)可以看出，經(jīng)過(guò)ESMDA-EnKF 算法更新之后，第Ⅷ個(gè)分區(qū)的滲透系數(shù)（P8）的分布最分散，不確定性最大，這可能是由于在此分區(qū)內(nèi)觀測(cè)數(shù)據(jù)較少，導(dǎo)致算法無(wú)法對(duì)參數(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)過(guò)ESMDA-EnKF更新之后的模型參數(shù)分布（參數(shù)編號(hào)解釋見表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ū)間有變大的趨勢(shì)，這可能是受預(yù)報(bào)前期模型運(yùn)行結(jié)果不確定性的累積影響，模型運(yùn)行后期地下水位與河流徑流量運(yùn)行結(jié)果產(chǎn)生了更大的不確定性。圖9 第二、三列分別為ESMDA-EnKF 算法運(yùn)行1、2 次之后的模型輸出值。可以看出，ESMDA-EnKF算法運(yùn)行一次之后，模型輸出的地下水位與河流徑流量的不確定性范圍明顯縮小，無(wú)論是數(shù)據(jù)起伏較小的地下水位還是數(shù)據(jù)起伏較大的河流徑流量，都收斂到了最優(yōu)值附近。第二次執(zhí)行ESMDA-EnKF 算法則基本消除了模型輸出結(jié)果的不確定性（對(duì)應(yīng)于圖9 的第三列）。

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

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

為了定量的分析地下水模型參數(shù)的不確定性，同時(shí)進(jìn)一步說(shuō)明ESMDA-EnKF 算法對(duì)地下水模型參數(shù)及輸出不確定性的更新效果，本文運(yùn)用Origin 2018 軟件對(duì)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)過(guò)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ù)同化的不確定性分析方法，通過(guò)包含觀測(cè)資料信息減小模型參數(shù)的不確定性?；谝褬?gòu)建的黑河流域中游地下水模型，利用提出的算法對(duì)模型參數(shù)進(jìn)行不確定性分析，探索了不同超參數(shù)對(duì)算法效果的影響，分析了含水層滲透系數(shù)、河床水力傳導(dǎo)系數(shù)、含水層給水度、灌溉回流系數(shù)的不確定性，評(píng)估了不同系數(shù)對(duì)觀測(cè)信息的敏感程度，主要結(jié)論如下：

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

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

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