地下水位预测:集合卡尔曼滤波(EnKF)应用概述

沈 晔，李海涛，黎 涛，李文鹏

（1.中国地质大学(北京)水资源与环境学院，北京 100083；2.中国地质环境监测院，北京 100081；3.中国地质调查局水文地质环境地质调查中心，河北 保定 071051）

地下水资源是城市供水的重要来源之一，特别是在中国的北方城市，由于地表水资源不足，在很大程度上都依赖地下水资源进行城市供水，如北京市2011年地下水资源占城市总供水量的59%［1］。随着城市化进程的不断发展，未来对水资源的需求量将会越来越大。由于不合理地开采地下水资源，已经造成了一系列的环境地质问题，如局部地区地下水位持续下降，地面沉降等。而近年来，全球气候不断变化，极端气候事件频发(如2010年的西南旱情，2011年的华北旱情，2012年北京的强降雨等事件)，对地下水的影响也很大，人类对地下水资源的需求仍在不断增加。在此背景之下，进行地下水位预报预警研究是一项十分重要和迫切的工作，已经成为水资源安全保障工作中的重中之重［2］。其中，短、中期地下水预报预警工作显得格外重要，对预报预警准确性的要求也越来越高。

目前，地下水位预测的方法很多，诸如地下水均衡法、简易类推法、数值法(有限差分法、有限元法)，周期分析法、地下水动力学法、相关分析法等［3］。随着计算机技术的不断发展，目前应用较多的是数值法。本文将以数值法为例，分析讨论集合卡尔曼滤波(Ensemble Kalman Filter，简称EnKF)在地下水位预报预警的应用前景。

1 集合卡尔曼滤波的发展过程

针对离散数据线性滤波问题，Kalman在1960年发表一篇非常有影响力的文章［4］。在文章中，Kalman提出了一种新的离散数据线性滤波问题的递归算法。该算法就是最早的卡尔曼滤波，或称为离散卡尔曼滤波。

该方法基本分为三个过程:①模型状态预测，②获取实际观测数据，③利用实际观测数据对模型状态进行更新。在实际应用中，可以不断获取观测数据，不断对模型状态进行更新，从而增强模型预测的准确性。自卡尔曼滤波的思想被提出来后，在线性优化的研究中得到广泛的应用，特别是在信号处理、优化控制、自主辅助导航等领域。

随着推广应用的不断深入，研究者们开始尝试将其引入到非线性问题的研究中。先将非线性问题做Taylor展开，然后近似线性化，忽略高阶导数项，再利用离散卡尔曼滤波进行求解。这种针对非线性问题的卡尔曼滤波后来被称为扩展卡尔曼滤波。然而，对于一些非线性问题，往往被忽略的高阶次项对计算结果具有不可忽略的影响，于是扩展卡尔曼滤波不再适用于这些研究中。

随着随机理论的不断发展，研究者们将其与卡尔曼滤波结合研究复杂的非线性问题。在此背景下，提出集合卡尔曼滤波(EnKF)。该方法通过随机产生大量的参数集合，对模型状态(如地下水位)进行预测，然后根据获取的地下水位观测数据对模型状态进行更新。而且，随着监测技术和试验手段的不断发展，获取的观测数据呈现出多样化的趋势，如可以通过抽水试验获取水文地质参数，通过监测获得高频率的地下水位数据等。利用EnKF可以对模型中的多种参数进行联合更新和数据同化，如在地下水模型中，可对模型中的状态参数(如水位或水头等)和模型参数(如渗透系数等)进行联合更新。

目前，EnKF在很多领域都得到了应用，如气象学等［5］。图1为集合卡尔曼滤波在气象预报中的应用实例:第一行为7月3日预测数据，第二行为7月4日预测数据，第三行为7月5日预测数据，每一天对后期预测的天气情况都会有所变化。这就是通过预测气象数据和实际监测气象数据相结合，利用EnKF及时将监测数据与模型数据进行同化，对模型中的一些参数进行更新。经过EnKF数据同化后，预测结果的准确性更高。

图1 EnKF在天气预报中的应用示意图Fig.1 Schematic dirgram showing the application of EnKF to weather forecast

近5年来，国外研究者们才将EnKF引入到水文地质相关领域的研究中来。如利用地下水水头观测数据和渗透系数试验数据，对地下水模型中的状态参数和模型参数进行同化［6］；利用实验场地下水流和溶质运移实验的实际水头与浓度等观测资料对水流和运移模型参数进行有效估计［7］；利用观测数据，并结合EnKF的特点，进行实时模拟［8］等。然而EnKF在地下水预测中的应用在国内并不多见。

2 集合卡尔曼滤波(EnKF)基本原理

本文将以地下水数值模型为例，阐述如何将EnKF应用于地下水预测中。地下水数值模型中包括很多参数，既包括模型的水文地质参数(如渗透系数、孔隙度等)，也包括了模型的状态参数(如水头、浓度、温度等)［9］。

在利用EnKF时，首先需要产生一个参数向量集合，或称为参数实现集合。随机产生参数实现的集合具体表现形式:

式中:X——参数实现集合；

x——参数实现集合中的一个参数实现；

n——参数实现集合中实现的个数；

A——模型的水文地质参数实现；

B——模型的状态参数实现；

a——模型水文地质参数实现的一个元素；

b——模型状态参数实现中的一个元素；

p——模型水文地质参数实现中参数的个数；

s——模型状态参数实现中参数的个数。

EnKF计算方法的核心步骤有两个:预测和更新。在数据同化的过程中，这两个步骤不断的交替进行。在整个参数实现集合中，对于每一个参数实现都需要进行预测和更新的步骤。为方便起见，在此以模型的一个状态参数——水头h为例进行说明。

第一步:模型状态参数预测。

模型状态参数预测是指根据t－1时刻的h值，通过模型计算，得到t时刻的h值:

式中:F——模型预测函数。

从t－1时刻到t时刻，一方面，模型中的水文地质参数基本保持不变，通过模型获取t时刻模型的水头预测值；另一方面，通过一定的监测手段(如人工监测或自动监测)，获得从t－1到t时刻的水头观测值。在得到同一时间段的水头预测值和观测值后，就可以进行EnKF的第二步——模型状态参数更新。

第二步:模型状态参数更新。

模型状态参数的更新是指通过一定的计算方法，将模型状态参数向量中的数值进行更新:

ε——观测误差，一般服从正态分布；

H——观测点矩阵(根据模型中观测点坐标获得)；

Gt——卡尔曼增益矩阵。

卡尔曼增益矩阵是假设后验误差方差达到最小而获得的，见式(5)。Evensen在2003年对集合卡尔曼滤波中卡尔曼增益矩阵的具体计算方法和步骤进行了详细描述［10］。

R——误差协方差矩阵。

完成了t时刻模型状态参数水头h的更新之后，就可以用更新的t时刻的h值来预测t+1时刻的h值，到了t+1时刻又获得一些观测数据，然后再进行t+1时刻的h值的更新，…。以此类推，通过不断的获得观测数据，对模型的状态参数h值进行不断的更新，始终使模型的预测值不至于与实际观测值之间偏差太大，从而增强了模型预测的准确性。

3 地下水位预测

地下水位预测是对地下水位的未来状态进行预测的一项工作:通过分析以往与地下水位相关的数据，如地下水位监测数据，地下水源汇项监测数据等，分析其动态变化规律，对未来的地下水位进行预测。地下水位预测数据的准确性在地下水位预警中占有很重要的地位。地下水位预警是通过一定的技术方法，结合本区的水文地质条件，根据地下水位数据，对区域的地下水水位进行预警，划分不同的预警等级，以提醒水资源管理部门和相关政府部门以及社会公众予以重视，防止地下水位持续下降，从而达到水资源的可持续开发利用的目的［11］。其预警的结果就是在预测的地下水位基础上分析而来的。因此，地下水位预测的准确性直接影响地下水预警结果的准确性。

随着计算机技术的不断发展，现阶段，区域地下水位预测一般采用数值模拟的手段。在进行数值模型模拟时，需要大量的数据作为支撑，如地下水位监测数据、地下水开采量，降雨入渗补给量以及水文地质参数等。利用数值模拟手段进行地下水位预测的基本步骤是利用某一时刻以前的地下水监测数据(如地下水位等)对模型的参数进行识别、校正，然后再利用校正后的模型对该时刻后的地下水位进行预测。这个过程中存在两个方面的问题:①在实际工作中，由于种种原因，这些数据并不能完全获取或者统计数据不够准确，给模型预测的准确性带来很大的影响，从而降低了地下水位预测的精度。②在模型预测的应用过程中，随着时间的推移，又获取了新的地下水位监测数据。但是，这些地下水位监测数据并没有被好好的在模型中利用。

4 地下水监测

中国传统的地下水监测主要是以人工监测为主，地下水监测频率较低，有的地方6次/月，有的地方3次/月，有的地方更低。而从区域上控制地下水位状态主要依靠统测数据(一般2次/年，丰水期1次，枯水期1次)。随着地下水监测自动化技术的不断发展，地下水自动监测传输技术已经在国内外得到了广泛应用［12］，如荷兰的DIVER，瑞士的KELLER，加拿大的LEVEL等，中国也出现了一系列的地下水自动监测传输设备，如中科光大产品，精诚华通产品等，大大提高了地下水监测的频率。从2003年以来，通过试点以及推广运行，地下水位自动监测网已经形成相当了规模，同时，也获取了大量的地下水监测数据，亟待二次开发利用!

5 集合卡尔曼滤波在地下水位预测中的应用分析

集合卡尔曼滤波是一个可以将连续监测数据与模型预测数据进行较好融合的工具。正如在第二节中阐述的，对于模型预测可以分为很多个预测时段。在地下水模型中，参数较多，模型中的不确定性(如概念模型的不确定性，模型参数的不确定性等)会给模型预测的结果带来一定的不确定性。某一时段预测地下水位的不确性会给后期预测结果带来更大的不确定，最终影响地下水位预测的准确性。随着时间推移，模型预测时段成为实际发生的时段。在这个过程中，产生新的地下水位监测数据。在获得某一时段的模型预测数据和实际监测数据后，利用EnKF对该时段模型的状态参数进行更新，进而利用更新后的模型状态参数进行下一时段的模型预测，使模型预测数据和实际监测数据不至于偏差太大。通过不断的获取监测数据，不断将其与模型预测数据进行融合，交替进行，以提高或保证模型预测的准确度。

6 EnKF的应用前景与展望

综合分析集合卡尔曼滤波的实际功能和地下水监测的实际情况，集合卡尔曼滤波在地下水位预测中的应用前景十分光明，不仅可以使获取的海量地下水位监测数据得到充分的利用，而且可以提高或保证地下水位预测的准确性，进而提高地下水位预报预警工作的准确性。因此，相信在不久的将来，随着集合卡尔曼滤波技术应用研究的不断深入，其在地下水位预测以至水文地质相关领域的研究中都将得到广泛的应用，它将对地下水位预报预警的实际工作产生长远的影响。因此，建议在对地下水位进行预报的工作中，引进集合卡尔曼滤波技术，将地下水位监测数据与模型预测数据进行较好的融合，提高模型预测的准确性。

［1］范庆莲.北京市水资源公报［R］.北京:北京市水务局，2011.［Beijing Water Resources Bulletin［R］.Beijing:Beijing Water Authority，2011.(in Chinese)］

［2］李丽，任永泰.地下水预警研究的发展与现状［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］吴林娜.地下水位预报方法的探讨［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滤波在气象数据同化中的发展与应用［J］.地球科学进展，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］薛禹群，谢春红.地下水数值模拟［M］.北京:科学出版社，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］李文鹏，郑跃军，郝爱兵.北京平原区地下水位预警初步研究［J］.地学前缘，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］周仰效，李文鹏.地下水监测、信息系统、模型及可持续发展［M］.北京:科学出版社，2011.［ZHOU Y X，LI W P.Groundwater Monitoring，Information System，Modeling and Sustainable Development［M］.Science Publishing House，2011.(in Chinese)］