考慮預報偏差的迭代式集合卡爾曼濾波在地下水水流數(shù)據(jù)同化中的應用

楊 運，吳吉春，駱乾坤，錢家忠

（1.南京大學地球科學與工程學院水科學系，江蘇 南京 210023；2.水利部淮河水利委員會，安徽 蚌埠 233001；3.合肥工業(yè)大學資源與環(huán)境工程學院，安徽 合肥 230009）

集合卡爾曼濾波（Ensemble Kalman Filter，EnKF）是一種順序數(shù)據(jù)同化方法，可以融合多源觀測信息實現(xiàn)系統(tǒng)狀態(tài)動態(tài)更新[1-6]。其在大氣科學[7]、海洋科學[8]、石油工程[9]等諸多領域得到了廣泛的應用。在水文地質領域，研究者通過同化水位數(shù)據(jù)[10]、濃度數(shù)據(jù)[11]、地球物理數(shù)據(jù)[12]等觀測資料，開展了水文地質參數(shù)反演、地下水污染源識別、地下水水流和污染物運移預測等研究。同時，為使EnKF 方法適應強非線性、非高斯問題，一些學者對其進行了改進。南統(tǒng)超等[13]引入局域化修正方法，減小了小集合樣本抽樣帶來的偽相關的影響。張秋汝等[14]提出一種改進的全局迭代集合卡爾曼濾波方法，并通過非飽和土壤水數(shù)據(jù)同化問題研究了不同迭代式方法在降低參數(shù)和變量不一致性上的性能。周海燕[15]提出正態(tài)轉換—集合卡爾曼濾波方法，在正態(tài)分布空間里進行同化更新，成功用于非高斯含水層參數(shù)識別，并進一步用于識別非高斯含水層中地下水污染源[16]。

預報模型是集合卡爾曼濾波方法的重要組成部分，以往研究多建立在預報模型準確的基礎上（即預報模型不存在錯誤或預報模型僅存在服從零均值高斯分布的白噪聲）。許多實際問題如地下水水流數(shù)據(jù)同化，由于野外水文地質條件復雜及認知的局限性，地下水水流模型結構很難準確刻畫，模型概化存在不確定性，而模型概化不準確必將導致預報偏差[17-20]。研究表明忽視這種偏差會造成數(shù)據(jù)同化效果不佳，帶來錯誤的系統(tǒng)估計或低估系統(tǒng)的不確定性[21]。針對模型概化存在不確定性的數(shù)據(jù)同化問題，目前處理方式主要有2 種：（1）基于貝葉斯模型平均理論，采用幾種備選預報模型分別進行數(shù)據(jù)同化，再對同化結果進行模型平均從而降低單個備選模型預報偏差的影響。如Xue 等[22]將EnKF 集合到模型平均框架中，應用于地下水水流模型邊界水頭大小和滲透系數(shù)場變異函數(shù)模型不確定條件下的數(shù)據(jù)同化。實際問題中很難判斷備選模型是否已涵蓋所有合理的模型結構[23]，因此采用這種方式存在一定局限性。（2）給預報模型增加一個偏差項，偏差項也利用觀測數(shù)據(jù)同步更新，獲得的預報模型認為是無偏的模型參與數(shù)據(jù)同化。Kollat等[24]用這種方式求解了砂箱實驗中地下水污染物泄露過程不確定情況下的污染物運移預測問題，研究表明基于偏差的EnKF 方法能夠預報得到接近于真實情況的濃度穿透曲線。Erdal 等[21]研究了非飽和帶土壤分層結構不確定條件下的水分運移問題，基于偏差的EnKF 方法顯著優(yōu)于未考慮偏差的數(shù)據(jù)同化結果。然而，在基于偏差的EnKF 方法的應用中，對于引入偏差項可能增大參數(shù)和變量不一致性的問題并未進行深入分析，且其適用性也需進一步研究。

因此，本文提出考慮預報偏差的迭代式集合卡爾曼濾波（Biasaware Ensemble Kalman Filter with Confirming Option，Bias-CEnKF）方法，并將其應用于地下水水流模型中邊界條件、初始條件、源匯項等分別存在概化不準確的情形時非均質含水層滲透系數(shù)場識別和地下水水流預測問題。

1 研究方法

1.1 標準集合卡爾曼濾波（EnKF）

標準EnKF 通過模型預測和同化更新2 個步驟實現(xiàn)數(shù)據(jù)的順序同化。主要過程為：

（1）定義tk時刻增廣狀態(tài)矩陣Sk，包括狀態(tài)變量矩陣mk、模型參數(shù)矩陣θk，即：

（2）對tk時刻狀態(tài)變量進行預測得到mfk：

式中：F—預報算子；

—tk-1時刻同化更新后的狀態(tài)變量；

模型參數(shù)；

（3）定義tk時刻的觀測數(shù)據(jù)矩陣dobs,k：

式中：Hk—觀測算子；

Strue,k—tk時刻真實狀態(tài)向量；

?k—觀測誤差矩陣。

（4）利用觀測數(shù)據(jù)進行同化更新：

—預測的狀態(tài)矩陣；

Kk—卡爾曼增益矩陣；

Rk—觀測誤差協(xié)方差矩陣；

—狀態(tài)誤差協(xié)方差矩陣；

—預測的狀態(tài)均值向量；

N—集合樣本數(shù)目；

j—集合樣本編號。

1.2 考慮預報偏差的迭代式集合卡爾曼濾波（Bias-CEnKF）

由于EnKF 在觀測時間點對狀態(tài)變量和模型參數(shù)同時做一次線性校正，沒有考慮參數(shù)和變量之間的非線性關系，可能造成兩者的不一致性。當考慮預報偏差后，偏差項也同時進行線性校正，進一步增加了不一致性。Wen 等[25]提出的迭代式集合卡爾曼濾波（Ensemble Kalman Filter with Confirming Option，CEn-KF）是一種從上一時刻到當前時刻的局部迭代方法，在保證計算效率的前提下能夠有效降低不一致性影響，獲得了廣泛的應用[14，26-27]。本文采用該迭代式方法構建Bias-CEnKF 方法，相比于標準集合卡爾曼濾波，主要變化為：

（1）增加一個偏差項bk，將bk與狀態(tài)變量mk、模型參數(shù)θk一起組成增廣狀態(tài)矩陣Sk，即式（1）修改為：

（2）對偏差項進行預報，并將其直接反饋至狀態(tài)變量預報模型，對模型狀態(tài)變量預測值進行修正，即式（2）修改為：

—tk-1時刻同化更新后的偏差項；

Λk—偏差項的時間相關系數(shù)；

wk—偏差噪聲；

—tk-1時刻經(jīng)過同化、重新預測后的狀態(tài)變量。

（3）在式（4）中，偏差項與狀態(tài)變量、模型參數(shù)一起進行同化更新，然后將更新后的返回至tk-1時刻，重新進行tk時刻的模型狀態(tài)變量預測：

—tk時刻同化更新后的偏差項。

將f作為下一時段的初始狀態(tài)變量，順序開展數(shù)據(jù)同化。

2 算例

2.1 算例概況

如圖1 所示，研究區(qū)為一個二維承壓含水層。含水層東西方向長500 m，南北方向長300 m，垂向厚2 m，平面均勻離散成50×30 個網(wǎng)格。含水層東部和西部邊界為定水頭邊界，水頭值分別為100，103 m，南部和北部為隔水邊界。含水層內有2 口抽水井和2 口注水井，自初始時刻開始運行，流量均恒定為100 m3/d。研究區(qū)共設置64 個水位（H）觀測點（含抽、注水井所在位置），12 個對數(shù)滲透系數(shù)（Y）觀測點。設定含水層滲透系數(shù)場（Y=lnK）服從對數(shù)正態(tài)分布，二階平穩(wěn)并可用二維可分離指數(shù)型協(xié)方差函數(shù)描述。真實Y場均值為0.5 m/d，標準差為1.2，x、y方向相關長度分別為120，60 m。采用Karhunen-Loeve(K-L)展開[28]生成一個對數(shù)滲透系數(shù)場，作為Y參考場，見圖1（b）。該滲透系數(shù)場下含水層的初始流場為無其他干擾情況下水流達到穩(wěn)定狀態(tài)時的流場，見圖1（a）。模型模擬期總時長10.0 d，分為20 個時段，時間步長0.5 d。采用MODFLOW 代碼模擬真實的水流模型，獲取水位觀測點每個時段的水位觀測數(shù)據(jù)[29]。模擬期前15 個時段的水位觀測值用于數(shù)據(jù)同化，后5 個時段的水位觀測值用于預測檢驗。

圖1 研究區(qū)初始流場與觀測位置和對數(shù)滲透系數(shù)Y 參考場（文獻[26]有修改）Fig.1 Initial flow field and observation locations (a) and reference field of Y (b) (modified from Ref.[26])

2.2 情景設置

本次研究將對數(shù)滲透系數(shù)場作為未知參數(shù)，通過數(shù)據(jù)同化對其進行識別。同化系統(tǒng)中地下水水流模型概化存在不確定性，對模型初始條件、邊界條件、源匯項的概化可能不準確，設置4 種情景：

情景1：真實地下水水流模型東西邊界為定水頭邊界，同化系統(tǒng)中地下水水流模型東西邊界為隔水邊界。

情景2：真實地下水水流模型東部和西部邊界水頭值分別為100，103 m，同化系統(tǒng)中地下水水流模型的東部和西部邊界分別為99.5，103.5 m。

情景3：真實地下水水流模型的初始流場見圖1（a），同化系統(tǒng)中地下水水流模型的初始水頭均設定為100 m。

情景4：真實地下水水流模型包含降雨入滲補給過程，入滲補給量0.001 m/d，同化系統(tǒng)中地下水水流模型未考慮降雨入滲補給過程。

2.3 參數(shù)設置

集合樣本數(shù)目設置為500，經(jīng)驗證能夠保證本文模擬結果的統(tǒng)計穩(wěn)定性。初始Y場統(tǒng)計特征與真實場有一定差距，均值為0，標準差為1.1，代表對Y場的先驗認知。水位觀測點觀測誤差服從均值為0、方差為2.5×10-5的高斯分布。Y觀測點觀測誤差服從均值為0、方差為1×10-6的高斯分布[10]。偏差噪聲服從均值為0、方差為0.01，x和y方向相關長度分別為300，180 m 的空間相關高斯分布。偏差項的時間相關系數(shù)設置為0.99[24]。

2.4 數(shù)據(jù)同化實現(xiàn)過程

采用Bias-CEnKF 同化水位觀測數(shù)據(jù)，識別對數(shù)滲透系數(shù)場的主要過程為：

①首先生成初始對數(shù)滲透系數(shù)場樣本。

②針對每個滲透系數(shù)場，采用MODFLOW 計算得到時段末各網(wǎng)格點水位預報值。

③預報偏差項，并修正時段末各網(wǎng)格點水位預報結果。

④根據(jù)修正后的各網(wǎng)格點水位預報值集合，計算狀態(tài)誤差協(xié)方差矩陣、卡爾曼增益矩陣。

⑤利用本時段水位觀測數(shù)據(jù)，同化更新各網(wǎng)格點對數(shù)滲透系數(shù)、水位、偏差項。

⑥重新返回到時段初，將同化更新后的對數(shù)滲透系數(shù)場代入MODFLOW 計算得到時段末各網(wǎng)格點水位值，利用更新后的偏差項予以修正。

⑦將修正后的各網(wǎng)格點水位值作為下一時段的初始狀態(tài)變量，重復步驟②～⑦，至同化結束。

考慮預報偏差的集合卡爾曼濾波（Bias-EnKF）與Bias-CEnKF 相比，少了迭代過程，即無上述步驟⑥，直接采用同化更新后的水位作為下一時段的初始狀態(tài)變量。標準EnKF 與Bias-CEnKF 相比，無偏差項，也無迭代過程。

2.5 同化效果評價標準

采用常用的均方根誤差（RMSE）作為評價標準，RMSE越小，代表估計量與真實值越接近。RMSE計算公式為：

式中：M—網(wǎng)格點數(shù)目；

—網(wǎng)格點i處的集合均值；

—網(wǎng)格點i處的真實值。

3 結果

3.1 同化結果

采用標準EnKF、Bias-EnKF 和Bias-CEnKF 分別對設置的4 種情景進行數(shù)據(jù)同化。

圖2 表示含水層Y場的同化結果與參考值的RMSE隨同化時間的變化情況。EnKF 計算得到的RMSE均大于Bias-EnKF 和Bias-CEnKF 的計算結果，且對情景1、情景3 和情景4 而言，隨著同化的進行，基于EnKF的RMSE逐漸增大或趨于一個較大值，說明濾波發(fā)散，集合未能收斂到一個合理的狀態(tài)。而在4 種情景下，Bias-EnKF 和Bias-CEnKF 計算得到的RMSE均隨著同化的進行逐步減小至穩(wěn)定。在情景1 和情景2 下，隨著同化步數(shù)的增加，Bias-CEnKF 的RMSE均小于Bias-EnKF。說明Bias-CEnKF 考慮了參數(shù)和變量的非線性關系，降低了兩者的不一致性，同化效果優(yōu)于Bias-EnKF。如對情景1 而言，Bias-EnKF 和Bias-CEnKF 計算得到的RMSE分別收斂到0.79，0.75。圖3列舉了情景1 下3 種方法同化得到的Y均值場?？梢钥闯?，在初始Y均值場比較均勻，與參考場（圖1b），差距較大的情況下，EnKF 同化得到的Y均值場與參考場顯著不同，而Bias-EnKF 和Bias-CEnKF 得到的Y均值場與參考場相似，基本反映了參考場的非均質分布特征。

圖2 不同情景下含水層對數(shù)滲透系數(shù)場RMSE 變化情況Fig.2 RMSE of Y under different scenarios

圖3 情景1 不同方法同化得到的對數(shù)滲透系數(shù)均值場Fig.3 Mean field of Y identified with different methods under scenario 1

圖4 為含水層水頭場的RMSE隨同化時間的變化情況。與Y場同化效果類似，EnKF 同化得到的水頭場RMSE最大。這表明存在模型預報偏差時，隨著同化的進行，觀測數(shù)據(jù)提供的信息不但未改善EnKF 的預測性能，而且對情景1 和情景4 而言，其對模型狀態(tài)的刻畫反而變得更差。相比之下，Bias-EnKF 和Bias-CEnKF 同化得到的水頭場RMSE在同化的中后期趨于穩(wěn)定，同化水頭場與真實水頭場接近。圖5 列舉了情景1 下含水層中幾個代表性點的水頭擬合結果。代表性點分別為靠近含水層西邊界、中間和東邊界處的點A（30 m，150 m）、B（250 m，150 m）和C（470 m，150 m）?？梢钥闯觯瑢拷吔缥恢玫狞cA 和點C，在同化初期，Bias-EnKF 和Bias-CEnKF 同化值與實際水頭有一定差距，隨著同化的進行，同化值與實際水頭擬合效果較好，對中間位置的點B，Bias-EnKF 和Bias-CEnKF 同化過程中同化值與實際水頭一直保持較好的一致性；而對EnKF 而言，點B 和點C 處的同化值與實際水頭差距隨著同化的進行并沒有減小，反而有變大趨勢。綜合說明，考慮預報偏差后，含水層水頭場的識別效果得到了有效提升。

圖4 不同情景下含水層水頭場RMSE 變化情況Fig.4 RMSE of H under different scenarios

圖5 情景1 代表性點水頭擬合結果對比Fig.5 Comparison of head-fitting results of the representative points under scenario 1

表1 為EnKF、Bias-EnKF 和Bias-CEnKF 同化末期的RMSE統(tǒng)計結果。與圖3 和圖4 綜合分析，在情景1 下Y場和水頭場的RMSE均最大，而情景4 下Y場和水頭場的RMSE均最小。說明對本文研究問題而言，含水層邊界性質的不準確刻畫帶來的預報偏差對同化結果的影響最大，降雨補給量的不準確刻畫帶來的預報偏差對同化結果的影響相對最小。就3 種方法同化效果而言，Bias-CEnKF 同化得到的RMSE均是最小的，表明其同化性能最好。

表1 EnKF、Bias-EnKF 和Bias-CEnKF 同化結果對比Table 1 Assimilation results of EnKF, Bias-EnKF and Bias-CEnKF

3.2 預測結果

因Bias-CEnKF 的同化性能最好，利用其同化結果，對后5 個時段的水頭場進行預測，并綜合分析同化階段及預測階段水頭預測值與實際值的接近程度。同樣，分別選取靠近含水層西邊界、中間和東邊界處的點A（30 m，150 m）、B（250 m，150 m）和C（470 m，150 m）作為觀測點進行分析。如圖6（a）—（c）所示，在情景1 下除邊界處觀測點A 和C 在同化前期水頭預測值與實際值有一定差距外，其余時段水頭預測值和實際值均擬合較好。同化前期觀測點A 和C 處水頭值與實際值差距較大是由于情景1 將模型東西邊界誤認為隔水邊界，而邊界性質對靠近邊界位置處的水頭影響較大所致。但是隨著同化步數(shù)的增加，邊界概化不準確的影響逐漸消失。在預測階段不同觀測點處水頭預測結果均與實際水頭值接近。圖6（d）—（f）為情景2 下的水頭擬合結果。與情景1 相比，不同階段3 個觀測點處的水頭預測值均與實際值擬合較好。這說明相對邊界性質而言，邊界水頭數(shù)值大小對同化效果的影響較小，同化前期可以通過偏差項彌補邊界水頭值誤差帶來的影響。圖6（g）—（i）為情景3下的水頭擬合結果。由于初始流場水頭值統(tǒng)一設置為100 m，在同化初期對水頭預測值也有較大影響。但是隨著同化的進行，初始流場概化不準確帶來的影響逐步減弱。圖6（j）—（l）為情景4 下的水頭擬合結果。受降雨補給影響，點B 處水頭呈現(xiàn)逐漸增大趨勢，邊界處點A 和C 因為靠近定水頭邊界，受降雨補給的影響相對較小。但是不同觀測點處同化和預測階段的水頭估計值與真實值均擬合得較好。

圖6 不同觀測點處水頭預測值與實際值對比Fig.6 Predicted and observed head values at different observation locations

圖7 為預測末期y=150 m 剖面處各點的水頭預測值與實際值對比結果。4 種情景下水頭實際值基本位于95%置信區(qū)間范圍內，且除在高水頭和低水頭位置有一點差距外，水頭預測均值與真實值均擬合較好。高水頭和低水頭位置處水頭擬合結果稍差的原因主要是通過數(shù)據(jù)同化并未完全復制真實的Y場，而滲透系數(shù)空間變異性對水頭的影響非常關鍵。

圖7 預測末期y=150 m 剖面上水頭預測值與實際值對比Fig.7 Predicted and observed head values in cross section y=150 m at the end of the prediction period

Bias-CEnKF 在同化初期的預測水頭值可能和真實值有一定差距，但是同化中后期和預測階段，水頭真實值均在預測的置信區(qū)間范圍內，且預測水頭值均與真實值擬合較好。在同一剖面位置處，由于受到含水層滲透系數(shù)空間變異性的影響，在高水頭及低水頭點處的水位擬合效果稍微變差，但整體擬合結果較好。因此，Bias-CEnKF 方法可以很好地彌補模型概化不準確對水位預報帶來的影響，具有較好的同化性能。

4 討論

含水層Y場和水頭場的同化結果表明，當模型概化不準確時，Bias-CEnKF 的同化性能顯著優(yōu)于EnKF，分析原因主要為前者利用偏差項對模型狀態(tài)變量預測值進行了修正。因此，為說明不同情景下偏差項如何減小模型概化不準確帶來的影響，對Bias-CEnKF 同化得到的偏差項結果作進一步分析。如圖8 所示，情景1 的偏差整體呈現(xiàn)西北部和東南部為負值、西南部和東北部為正值的特征，且越靠近邊界位置偏差越大。這是由于研究區(qū)西北部和東南部為抽水井，真實的定水頭邊界條件下，會引起邊界水流補給含水層，而情景1 假定東西邊界為隔水邊界，無水流補給含水層，在系統(tǒng)同化過程中需要通過偏差項增加源項（負值）以彌補模型概化不準確的影響。同理，研究區(qū)西南部和東北部為注水井，真實的定水頭邊界條件下，會有水流通過定水頭邊界流出含水層，而情景1 假定邊界為隔水邊界，同化系統(tǒng)通過偏差項增加匯項（正值）以彌補模型概化不準確的影響。

圖8 不同情景下Bias-CEnKF 同化得到的偏差項的均值場Fig.8 Mean field of bias corrections based on the Bias-CEnKF under different scenarios

情景2 的偏差整體呈現(xiàn)西部為正值、東部為負值、向中間位置趨于0 的特征。這是由于真實的西部邊界水頭為103 m，情景2 假定西部邊界水頭（103.5 m）高于真實水頭值，系統(tǒng)同化過程中需要通過匯項（正值）降低高水頭的影響。同理，真實的東部邊界水頭為100 m，情景2 假定東部邊界水頭（99.5 m）低于真實水頭值，在系統(tǒng)同化過程中需要通過源項（負值）降低低水頭的影響。此外，邊界水頭對含水層水頭分布的影響通常表現(xiàn)為由邊界向中間逐漸減弱的趨勢。因此，偏差項也由邊界向中間逐漸趨于0。

情景3 的偏差項無明顯規(guī)律性。這是由于理論上初始水頭對同化系統(tǒng)的影響在順序數(shù)據(jù)同化過程中會逐漸減弱甚至消失，且初始水頭不像邊界或源匯項等給予含水層以持續(xù)的流入或流出影響。

情景4 的偏差基本為負值并由中間向邊界逐漸趨于0。這是由于情景4 忽視了真實的降雨補給，同化系統(tǒng)需要通過源項（負值）體現(xiàn)降雨補給的影響。同時，越往邊界受定水頭邊界影響越大，降雨補給影響相對減弱，因此，偏差也趨于0。

5 結論

（1）模型概化不確定條件下，模型概化不準確時，使用標準EnKF 進行數(shù)據(jù)同化，可能會導致濾波發(fā)散，造成同化失敗，而Bias-CEnKF 能夠識別得到接近于真實場的Y場和水頭場，且其同化性能優(yōu)于Bias-EnKF。

（2）Bias-CEnKF 同化系統(tǒng)中的偏差項在物理意義上相當于用源匯項彌補模型概化不準確帶來的影響。對本文4 種情景而言，Bias-CEnKF 能夠識別得到合理的偏差項分布特征。

（3）利用Bias-CEnKF 同化得到的狀態(tài)變量、模型參數(shù)以及偏差項進行預測，均獲得了可靠的預測結果。對于實際復雜問題，模型概化不準確對模擬結果的影響可能會隨著時間增加而累積。要進行長期的合理預報，一方面應盡可能對地下水水流模型進行準確刻畫，另一方面應及時采用最新觀測數(shù)據(jù)進行同化更新，實時修正偏差項。