地下水位預(yù)測:集合卡爾曼濾波(EnKF)應(yīng)用概述

沈 曄，李海濤，黎 濤，李文鵬

（1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；2.中國地質(zhì)環(huán)境監(jiān)測院，北京 100081；3.中國地質(zhì)調(diào)查局水文地質(zhì)環(huán)境地質(zhì)調(diào)查中心，河北 保定 071051）

地下水資源是城市供水的重要來源之一，特別是在中國的北方城市，由于地表水資源不足，在很大程度上都依賴地下水資源進行城市供水，如北京市2011年地下水資源占城市總供水量的59%［1］。隨著城市化進程的不斷發(fā)展，未來對水資源的需求量將會越來越大。由于不合理地開采地下水資源，已經(jīng)造成了一系列的環(huán)境地質(zhì)問題，如局部地區(qū)地下水位持續(xù)下降，地面沉降等。而近年來，全球氣候不斷變化，極端氣候事件頻發(fā)(如2010年的西南旱情，2011年的華北旱情，2012年北京的強降雨等事件)，對地下水的影響也很大，人類對地下水資源的需求仍在不斷增加。在此背景之下，進行地下水位預(yù)報預(yù)警研究是一項十分重要和迫切的工作，已經(jīng)成為水資源安全保障工作中的重中之重［2］。其中，短、中期地下水預(yù)報預(yù)警工作顯得格外重要，對預(yù)報預(yù)警準(zhǔn)確性的要求也越來越高。

目前，地下水位預(yù)測的方法很多，諸如地下水均衡法、簡易類推法、數(shù)值法(有限差分法、有限元法)，周期分析法、地下水動力學(xué)法、相關(guān)分析法等［3］。隨著計算機技術(shù)的不斷發(fā)展，目前應(yīng)用較多的是數(shù)值法。本文將以數(shù)值法為例，分析討論集合卡爾曼濾波(Ensemble Kalman Filter，簡稱EnKF)在地下水位預(yù)報預(yù)警的應(yīng)用前景。

1 集合卡爾曼濾波的發(fā)展過程

針對離散數(shù)據(jù)線性濾波問題，Kalman在1960年發(fā)表一篇非常有影響力的文章［4］。在文章中，Kalman提出了一種新的離散數(shù)據(jù)線性濾波問題的遞歸算法。該算法就是最早的卡爾曼濾波，或稱為離散卡爾曼濾波。

該方法基本分為三個過程:①模型狀態(tài)預(yù)測，②獲取實際觀測數(shù)據(jù)，③利用實際觀測數(shù)據(jù)對模型狀態(tài)進行更新。在實際應(yīng)用中，可以不斷獲取觀測數(shù)據(jù)，不斷對模型狀態(tài)進行更新，從而增強模型預(yù)測的準(zhǔn)確性。自卡爾曼濾波的思想被提出來后，在線性優(yōu)化的研究中得到廣泛的應(yīng)用，特別是在信號處理、優(yōu)化控制、自主輔助導(dǎo)航等領(lǐng)域。

隨著推廣應(yīng)用的不斷深入，研究者們開始嘗試將其引入到非線性問題的研究中。先將非線性問題做Taylor展開，然后近似線性化，忽略高階導(dǎo)數(shù)項，再利用離散卡爾曼濾波進行求解。這種針對非線性問題的卡爾曼濾波后來被稱為擴展卡爾曼濾波。然而，對于一些非線性問題，往往被忽略的高階次項對計算結(jié)果具有不可忽略的影響，于是擴展卡爾曼濾波不再適用于這些研究中。

隨著隨機理論的不斷發(fā)展，研究者們將其與卡爾曼濾波結(jié)合研究復(fù)雜的非線性問題。在此背景下，提出集合卡爾曼濾波(EnKF)。該方法通過隨機產(chǎn)生大量的參數(shù)集合，對模型狀態(tài)(如地下水位)進行預(yù)測，然后根據(jù)獲取的地下水位觀測數(shù)據(jù)對模型狀態(tài)進行更新。而且，隨著監(jiān)測技術(shù)和試驗手段的不斷發(fā)展，獲取的觀測數(shù)據(jù)呈現(xiàn)出多樣化的趨勢，如可以通過抽水試驗獲取水文地質(zhì)參數(shù)，通過監(jiān)測獲得高頻率的地下水位數(shù)據(jù)等。利用EnKF可以對模型中的多種參數(shù)進行聯(lián)合更新和數(shù)據(jù)同化，如在地下水模型中，可對模型中的狀態(tài)參數(shù)(如水位或水頭等)和模型參數(shù)(如滲透系數(shù)等)進行聯(lián)合更新。

目前，EnKF在很多領(lǐng)域都得到了應(yīng)用，如氣象學(xué)等［5］。圖1為集合卡爾曼濾波在氣象預(yù)報中的應(yīng)用實例:第一行為7月3日預(yù)測數(shù)據(jù)，第二行為7月4日預(yù)測數(shù)據(jù)，第三行為7月5日預(yù)測數(shù)據(jù)，每一天對后期預(yù)測的天氣情況都會有所變化。這就是通過預(yù)測氣象數(shù)據(jù)和實際監(jiān)測氣象數(shù)據(jù)相結(jié)合，利用EnKF及時將監(jiān)測數(shù)據(jù)與模型數(shù)據(jù)進行同化，對模型中的一些參數(shù)進行更新。經(jīng)過EnKF數(shù)據(jù)同化后，預(yù)測結(jié)果的準(zhǔn)確性更高。

圖1 EnKF在天氣預(yù)報中的應(yīng)用示意圖Fig.1 Schematic dirgram showing the application of EnKF to weather forecast

近5年來，國外研究者們才將EnKF引入到水文地質(zhì)相關(guān)領(lǐng)域的研究中來。如利用地下水水頭觀測數(shù)據(jù)和滲透系數(shù)試驗數(shù)據(jù)，對地下水模型中的狀態(tài)參數(shù)和模型參數(shù)進行同化［6］；利用實驗場地下水流和溶質(zhì)運移實驗的實際水頭與濃度等觀測資料對水流和運移模型參數(shù)進行有效估計［7］；利用觀測數(shù)據(jù)，并結(jié)合EnKF的特點，進行實時模擬［8］等。然而EnKF在地下水預(yù)測中的應(yīng)用在國內(nèi)并不多見。

2 集合卡爾曼濾波(EnKF)基本原理

本文將以地下水?dāng)?shù)值模型為例，闡述如何將EnKF應(yīng)用于地下水預(yù)測中。地下水?dāng)?shù)值模型中包括很多參數(shù)，既包括模型的水文地質(zhì)參數(shù)(如滲透系數(shù)、孔隙度等)，也包括了模型的狀態(tài)參數(shù)(如水頭、濃度、溫度等)［9］。

在利用EnKF時，首先需要產(chǎn)生一個參數(shù)向量集合，或稱為參數(shù)實現(xiàn)集合。隨機產(chǎn)生參數(shù)實現(xiàn)的集合具體表現(xiàn)形式:

式中:X——參數(shù)實現(xiàn)集合；

x——參數(shù)實現(xiàn)集合中的一個參數(shù)實現(xiàn)；

n——參數(shù)實現(xiàn)集合中實現(xiàn)的個數(shù)；

A——模型的水文地質(zhì)參數(shù)實現(xiàn)；

B——模型的狀態(tài)參數(shù)實現(xiàn)；

a——模型水文地質(zhì)參數(shù)實現(xiàn)的一個元素；

b——模型狀態(tài)參數(shù)實現(xiàn)中的一個元素；

p——模型水文地質(zhì)參數(shù)實現(xiàn)中參數(shù)的個數(shù)；

s——模型狀態(tài)參數(shù)實現(xiàn)中參數(shù)的個數(shù)。

EnKF計算方法的核心步驟有兩個:預(yù)測和更新。在數(shù)據(jù)同化的過程中，這兩個步驟不斷的交替進行。在整個參數(shù)實現(xiàn)集合中，對于每一個參數(shù)實現(xiàn)都需要進行預(yù)測和更新的步驟。為方便起見，在此以模型的一個狀態(tài)參數(shù)——水頭h為例進行說明。

第一步:模型狀態(tài)參數(shù)預(yù)測。

模型狀態(tài)參數(shù)預(yù)測是指根據(jù)t－1時刻的h值，通過模型計算，得到t時刻的h值:

式中:F——模型預(yù)測函數(shù)。

從t－1時刻到t時刻，一方面，模型中的水文地質(zhì)參數(shù)基本保持不變，通過模型獲取t時刻模型的水頭預(yù)測值；另一方面，通過一定的監(jiān)測手段(如人工監(jiān)測或自動監(jiān)測)，獲得從t－1到t時刻的水頭觀測值。在得到同一時間段的水頭預(yù)測值和觀測值后，就可以進行EnKF的第二步——模型狀態(tài)參數(shù)更新。

第二步:模型狀態(tài)參數(shù)更新。

模型狀態(tài)參數(shù)的更新是指通過一定的計算方法，將模型狀態(tài)參數(shù)向量中的數(shù)值進行更新:

ε——觀測誤差，一般服從正態(tài)分布；

H——觀測點矩陣(根據(jù)模型中觀測點坐標(biāo)獲得)；

Gt——卡爾曼增益矩陣。

卡爾曼增益矩陣是假設(shè)后驗誤差方差達到最小而獲得的，見式(5)。Evensen在2003年對集合卡爾曼濾波中卡爾曼增益矩陣的具體計算方法和步驟進行了詳細描述［10］。

R——誤差協(xié)方差矩陣。

完成了t時刻模型狀態(tài)參數(shù)水頭h的更新之后，就可以用更新的t時刻的h值來預(yù)測t+1時刻的h值，到了t+1時刻又獲得一些觀測數(shù)據(jù)，然后再進行t+1時刻的h值的更新，…。以此類推，通過不斷的獲得觀測數(shù)據(jù)，對模型的狀態(tài)參數(shù)h值進行不斷的更新，始終使模型的預(yù)測值不至于與實際觀測值之間偏差太大，從而增強了模型預(yù)測的準(zhǔn)確性。

3 地下水位預(yù)測

地下水位預(yù)測是對地下水位的未來狀態(tài)進行預(yù)測的一項工作:通過分析以往與地下水位相關(guān)的數(shù)據(jù)，如地下水位監(jiān)測數(shù)據(jù)，地下水源匯項監(jiān)測數(shù)據(jù)等，分析其動態(tài)變化規(guī)律，對未來的地下水位進行預(yù)測。地下水位預(yù)測數(shù)據(jù)的準(zhǔn)確性在地下水位預(yù)警中占有很重要的地位。地下水位預(yù)警是通過一定的技術(shù)方法，結(jié)合本區(qū)的水文地質(zhì)條件，根據(jù)地下水位數(shù)據(jù)，對區(qū)域的地下水水位進行預(yù)警，劃分不同的預(yù)警等級，以提醒水資源管理部門和相關(guān)政府部門以及社會公眾予以重視，防止地下水位持續(xù)下降，從而達到水資源的可持續(xù)開發(fā)利用的目的［11］。其預(yù)警的結(jié)果就是在預(yù)測的地下水位基礎(chǔ)上分析而來的。因此，地下水位預(yù)測的準(zhǔn)確性直接影響地下水預(yù)警結(jié)果的準(zhǔn)確性。

隨著計算機技術(shù)的不斷發(fā)展，現(xiàn)階段，區(qū)域地下水位預(yù)測一般采用數(shù)值模擬的手段。在進行數(shù)值模型模擬時，需要大量的數(shù)據(jù)作為支撐，如地下水位監(jiān)測數(shù)據(jù)、地下水開采量，降雨入滲補給量以及水文地質(zhì)參數(shù)等。利用數(shù)值模擬手段進行地下水位預(yù)測的基本步驟是利用某一時刻以前的地下水監(jiān)測數(shù)據(jù)(如地下水位等)對模型的參數(shù)進行識別、校正，然后再利用校正后的模型對該時刻后的地下水位進行預(yù)測。這個過程中存在兩個方面的問題:①在實際工作中，由于種種原因，這些數(shù)據(jù)并不能完全獲取或者統(tǒng)計數(shù)據(jù)不夠準(zhǔn)確，給模型預(yù)測的準(zhǔn)確性帶來很大的影響，從而降低了地下水位預(yù)測的精度。②在模型預(yù)測的應(yīng)用過程中，隨著時間的推移，又獲取了新的地下水位監(jiān)測數(shù)據(jù)。但是，這些地下水位監(jiān)測數(shù)據(jù)并沒有被好好的在模型中利用。

4 地下水監(jiān)測

中國傳統(tǒng)的地下水監(jiān)測主要是以人工監(jiān)測為主，地下水監(jiān)測頻率較低，有的地方6次/月，有的地方3次/月，有的地方更低。而從區(qū)域上控制地下水位狀態(tài)主要依靠統(tǒng)測數(shù)據(jù)(一般2次/年，豐水期1次，枯水期1次)。隨著地下水監(jiān)測自動化技術(shù)的不斷發(fā)展，地下水自動監(jiān)測傳輸技術(shù)已經(jīng)在國內(nèi)外得到了廣泛應(yīng)用［12］，如荷蘭的DIVER，瑞士的KELLER，加拿大的LEVEL等，中國也出現(xiàn)了一系列的地下水自動監(jiān)測傳輸設(shè)備，如中科光大產(chǎn)品，精誠華通產(chǎn)品等，大大提高了地下水監(jiān)測的頻率。從2003年以來，通過試點以及推廣運行，地下水位自動監(jiān)測網(wǎng)已經(jīng)形成相當(dāng)了規(guī)模，同時，也獲取了大量的地下水監(jiān)測數(shù)據(jù)，亟待二次開發(fā)利用!

5 集合卡爾曼濾波在地下水位預(yù)測中的應(yīng)用分析

集合卡爾曼濾波是一個可以將連續(xù)監(jiān)測數(shù)據(jù)與模型預(yù)測數(shù)據(jù)進行較好融合的工具。正如在第二節(jié)中闡述的，對于模型預(yù)測可以分為很多個預(yù)測時段。在地下水模型中，參數(shù)較多，模型中的不確定性(如概念模型的不確定性，模型參數(shù)的不確定性等)會給模型預(yù)測的結(jié)果帶來一定的不確定性。某一時段預(yù)測地下水位的不確性會給后期預(yù)測結(jié)果帶來更大的不確定，最終影響地下水位預(yù)測的準(zhǔn)確性。隨著時間推移，模型預(yù)測時段成為實際發(fā)生的時段。在這個過程中，產(chǎn)生新的地下水位監(jiān)測數(shù)據(jù)。在獲得某一時段的模型預(yù)測數(shù)據(jù)和實際監(jiān)測數(shù)據(jù)后，利用EnKF對該時段模型的狀態(tài)參數(shù)進行更新，進而利用更新后的模型狀態(tài)參數(shù)進行下一時段的模型預(yù)測，使模型預(yù)測數(shù)據(jù)和實際監(jiān)測數(shù)據(jù)不至于偏差太大。通過不斷的獲取監(jiān)測數(shù)據(jù)，不斷將其與模型預(yù)測數(shù)據(jù)進行融合，交替進行，以提高或保證模型預(yù)測的準(zhǔn)確度。

6 EnKF的應(yīng)用前景與展望

綜合分析集合卡爾曼濾波的實際功能和地下水監(jiān)測的實際情況，集合卡爾曼濾波在地下水位預(yù)測中的應(yīng)用前景十分光明，不僅可以使獲取的海量地下水位監(jiān)測數(shù)據(jù)得到充分的利用，而且可以提高或保證地下水位預(yù)測的準(zhǔn)確性，進而提高地下水位預(yù)報預(yù)警工作的準(zhǔn)確性。因此，相信在不久的將來，隨著集合卡爾曼濾波技術(shù)應(yīng)用研究的不斷深入，其在地下水位預(yù)測以至水文地質(zhì)相關(guān)領(lǐng)域的研究中都將得到廣泛的應(yīng)用，它將對地下水位預(yù)報預(yù)警的實際工作產(chǎn)生長遠的影響。因此，建議在對地下水位進行預(yù)報的工作中，引進集合卡爾曼濾波技術(shù)，將地下水位監(jiān)測數(shù)據(jù)與模型預(yù)測數(shù)據(jù)進行較好的融合，提高模型預(yù)測的準(zhǔn)確性。

［1］范慶蓮.北京市水資源公報［R］.北京:北京市水務(wù)局，2011.［Beijing Water Resources Bulletin［R］.Beijing:Beijing Water Authority，2011.(in Chinese)］

［2］李麗，任永泰.地下水預(yù)警研究的發(fā)展與現(xiàn)狀［J］.邊疆經(jīng)濟與文化，2010，80(8):24 －25.［LI L，REN Y T.Development and status of groundwater warning Research［J］.The Board Economy and Culture，2010，80(8):24 －25.(in Chinese)］

［3］吳林娜.地下水位預(yù)報方法的探討［J］.陜西水利科技，1996，113(4):63 － 65.［WU L N.Discussion on methods for groundwater level prediction［J］.Shanxi Hydro-techniques，1996，113(4):63 － 65.(in Chinese)］

［4］KalmanR W.A new approach to linear filtering and prediction problems［J］.Transaction of the ASME Journal of basic engineering，1960，82(Series D):35－45.

［5］高山紅，吳增茂，謝紅琴.Kalman濾波在氣象數(shù)據(jù)同化中的發(fā)展與應(yīng)用［J］.地球科學(xué)進展，2000，15(5):571 －575.［GAO S H，WU Z M，XIE H Q.The development and application of Kalman filters in meteorological data assimilation［J］.Advances in earth sciences，2000，15(5):571 －575.(in Chinese)］

［6］CHEN Y，ZHANG D.Data assimilation for transient flow in geologic formations via ensemble Kalman filter［J］.Advances in Water Resources，2006，29(8):1107－122.

［7］LIU G S，CH EN Y，ZHANG D.Investigation of flow and transport processes at the MADE site using ensemble Kalman filter［J］.Advances in Water Resources，2008，31(7):975 －986.

［8］FRANSSEN H J H，K INZELBACH W.Real time groundwater flow modeling with the Ensemble Kalman Filter:Joint estimation of states and parameters and the filter inbreeding problem ［J］.Water Resources Research，2008，44(9):21.

［9］薛禹群，謝春紅.地下水?dāng)?shù)值模擬［M］.北京:科學(xué)出版社，2007.［XUE Y Q，XIE C H.Numerical Simulation for Groundwater［M］.Science Publishing House，2007.(in Chinese)］

［10］EvensenG.The Ensemble Kalman Filter:theoretical formation and practical implementation［J］.Ocean Dynamics，2003，53:343 －367.

［11］李文鵬，鄭躍軍，郝愛兵.北京平原區(qū)地下水位預(yù)警初步研究［J］.地學(xué)前緣，2010，17(6):166 －173.［LI W P，ZHENG Y J，HAO A B.A preliminary study of groundwater level pre-warning in Beijing Plain［J］.Earth Science Frontiers，2010，17(6):166 － 173.(in Chinese)］

［12］周仰效，李文鵬.地下水監(jiān)測、信息系統(tǒng)、模型及可持續(xù)發(fā)展［M］.北京:科學(xué)出版社，2011.［ZHOU Y X，LI W P.Groundwater Monitoring，Information System，Modeling and Sustainable Development［M］.Science Publishing House，2011.(in Chinese)］