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

（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）

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

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

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

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

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

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

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

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

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

近5年來，國外研究者們才將EnKF引入到水文地質(zhì)相關(guān)領(lǐng)域的研究中來。如利用地下水水頭觀測數(shù)據(jù)和滲透系數(shù)試驗(yàn)數(shù)據(jù)，對地下水模型中的狀態(tài)參數(shù)和模型參數(shù)進(jìn)行同化［6］；利用實(shí)驗(yàn)場地下水流和溶質(zhì)運(yùn)移實(shí)驗(yàn)的實(shí)際水頭與濃度等觀測資料對水流和運(yùn)移模型參數(shù)進(jìn)行有效估計(jì)［7］；利用觀測數(shù)據(jù)，并結(jié)合EnKF的特點(diǎn)，進(jìn)行實(shí)時(shí)模擬［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時(shí)，首先需要產(chǎn)生一個(gè)參數(shù)向量集合，或稱為參數(shù)實(shí)現(xiàn)集合。隨機(jī)產(chǎn)生參數(shù)實(shí)現(xiàn)的集合具體表現(xiàn)形式:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Gt——卡爾曼增益矩陣。

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

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

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

3 地下水位預(yù)測

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

隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展，現(xiàn)階段，區(qū)域地下水位預(yù)測一般采用數(shù)值模擬的手段。在進(jìn)行數(shù)值模型模擬時(shí)，需要大量的數(shù)據(jù)作為支撐，如地下水位監(jiān)測數(shù)據(jù)、地下水開采量，降雨入滲補(bǔ)給量以及水文地質(zhì)參數(shù)等。利用數(shù)值模擬手段進(jìn)行地下水位預(yù)測的基本步驟是利用某一時(shí)刻以前的地下水監(jiān)測數(shù)據(jù)(如地下水位等)對模型的參數(shù)進(jìn)行識(shí)別、校正，然后再利用校正后的模型對該時(shí)刻后的地下水位進(jìn)行預(yù)測。這個(gè)過程中存在兩個(gè)方面的問題:①在實(shí)際工作中，由于種種原因，這些數(shù)據(jù)并不能完全獲取或者統(tǒng)計(jì)數(shù)據(jù)不夠準(zhǔn)確，給模型預(yù)測的準(zhǔn)確性帶來很大的影響，從而降低了地下水位預(yù)測的精度。②在模型預(yù)測的應(yīng)用過程中，隨著時(shí)間的推移，又獲取了新的地下水位監(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)測自動(dòng)化技術(shù)的不斷發(fā)展，地下水自動(dòng)監(jiān)測傳輸技術(shù)已經(jīng)在國內(nèi)外得到了廣泛應(yīng)用［12］，如荷蘭的DIVER，瑞士的KELLER，加拿大的LEVEL等，中國也出現(xiàn)了一系列的地下水自動(dòng)監(jiān)測傳輸設(shè)備，如中科光大產(chǎn)品，精誠華通產(chǎn)品等，大大提高了地下水監(jiān)測的頻率。從2003年以來，通過試點(diǎn)以及推廣運(yùn)行，地下水位自動(dòng)監(jiān)測網(wǎng)已經(jīng)形成相當(dāng)了規(guī)模，同時(shí)，也獲取了大量的地下水監(jiān)測數(shù)據(jù)，亟待二次開發(fā)利用!

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

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

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

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

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

［2］李麗，任永泰.地下水預(yù)警研究的發(fā)展與現(xiàn)狀［J］.邊疆經(jīng)濟(jì)與文化，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ù)報(bào)方法的探討［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é)進(jìn)展，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)］