专利名称：一种水文地质参数的估计方法
技术领域：
本发明涉及水文学及水资源领域，具体是一种水文地质参数估计方法。
背景技术：
含水层性质因表现出明显的空间变异性，目前无法进行完全测量或描述。然而， 表征含水层性质的主要水文地质参数又是地下水流和运移数值模型的关键因素和主要的 不确定性来源。随着数量和种类都不断增加的连续原位观测资料的出现，人们迫切需要高 效实用的方法挖掘这些数据的价值以获得对模型参数更客观详细的认识。集合卡尔曼滤 波(Ensemble Kalman Filter, EnKF)作为一种有效的连续数据同化工具，由于其对非线 性和大规模问题的适用性、执行的简便性以及低计算消耗，得到了广泛的认可和应用。自 EVENSEN于1994年首次提出EnKF至今，EnKF已在海洋学、气象学、石油工程、气候学、水文 学等领域发挥重要作用。过去三年中，EnKF方法开始应用于水文地质学领域的数据同化和 含水层参数估计。CHEN和ZHANG在二维和三维理想模型中应用EnKF同化地下水水头和渗 透系数观测资料，获得了整个渗透系数场的合理估计。LIU等借助EnKF融合MADE实验场地 下水流和溶质运移调查实验的实际水头与浓度观测资料，对渗透系数和运移模型参数进行 了有效估计。FRANSSEN和KINZELBACH在二维理想模型中研究了不同"导水系数_补给量" 设置下EnKF的执行效果，并详细讨论了处理滤波近交的方法。 实际上，随着实际地下水问题复杂程度和研究尺度的迅速增加，使得标准集合卡 尔曼滤波需要的集合规模也迅速增加。另一方面，模拟复杂的地下水系统往往需要更多的 计算量。单个实现复杂度的增加和实现总数的增加双重因素造成的计算负担使当前的计算
能力很难承受。若减小实现数即集合大小，会使所获得的状态误差协方差受到的取样噪声 影响更加严重，从而使协方差精度降低，出现不合实际的伪相关，最终影响EnKF的执行效 果。这种干扰对于强烈非线性或高维问题尤其明显。地下水系统中两点间空间距离增大时 内在联系通常也会减弱。距离增加而造成的真实相关性减弱，以及集合减小而造成的绝对 取样误差增加，都会造成相对误差的增大，使估计的相关性受到更严重的干扰[8]。取样噪 声使本应趋于零的相关性不为零，相当于虚构了变量之间的相关性；伪相关会导致集合的 过度修正和过早收敛，其后果是集合变异性极易崩溃，EnKF无法有效吸收观测资料和同化 失败。
发明内容
本发明所要解决的技术问题是提供一种改善集合卡尔曼滤波进行水文地质参数
估计效果的方法——局域化集合卡尔曼滤波，该方法能有效过滤有限集合估计中的协方差
伪相关，大大提高协方差和水文地质参数的估计精度，改善滤波特性。 本发明所述的一种水文地质反求参数的数据同化方法，其包括以下步骤 1)针对实地水文地质条件建立概念模型和数值模型，充分考虑地质条件和边界条
件，确认研究区地质统计学的适用性，即满足二阶平稳假设和本征假设；
2)对研究区应用地质统计学方法估计参数场的协方差结构C(h)和协方差参数 入x， A y及、作为地质统计先验信息； 3)根据参数的地质统计先验信息，生成一组Ne个参数场的随机实现，作为取样集 合(即参数场的随机模拟)； 4)将取样集合中的Ne个实现分别用于先前建立的数值模型，进行短时间上的演化 计算，即获得k个模型预报；用公式表达如下 s: W = f[s;; W 其中S是状态向量，包含模型参数，变量及其他需要的观测数据；t表示时间；F是
预报算子，代表所建立的模型的作用；上标'f'和'u'分别代表预报和更新过程；下标'r' 是集合成员序数； 5)统计这些实现的预报结果，根据集合预报原理估计模型中所有参数及变量之间 的均值和协方差，即一阶矩和二阶矩，分别由下面两式计算
1
Sf|
c/ (,) ^it|:{[s〗(0-〈sf w〉][g 〈sf w〉]T}+c"o 上标'T
Cs是估计的状态误差协方差矩阵；CME是模型误差矩阵；上标'f'代表预报新过程； 表示矩阵转置(下同)；
6)考虑参数和变量的空间位置关系，并将这种关系以舒尔积权重的形式添加到集
合估计的协方差中得到修正的协方差
小-xi)ftii
■SP
p。cs = 7)将修改后的协方差用于卡尔曼滤波计算，输入当前可用的新观测数据，基于最 小二乘原理获得状态，包括参数和变量，的优化估计，更新集合状态； 8)以更新的集合状态为初始状态，进行下一时间段的预报，即重复4) 7)，直至同 化时间全部完成。 本发明考虑了集合估计中忽视的变量间的空间相关性，在协方差估计结果中通过 人为修正的方式强调这种基于空间位置的内在联系，有效地减弱了有限集合中获得的估计 结果所包含的取样误差。该方法有机集成了多个理论和方法的长处，且使用方便。本方法 的每一步都有充分的理论保证其合理性、必要性，有着同类参数估计问题下的普遍适用性。
图1是一个理想二维平面矩形承压含水层模拟计算区； 图2是Y参考场； 图3是第一时段真实等水头线和点(80m， 720m)处水头值与整个水头场的协方差， 其中(a)是等水头线；(b) (c) (d)分别为ensemble A、B和C的水头协方差。空心方形表示
4点(80m，720m)的位置，实心圆和三角分别表示抽、注水井的位置； 图4是ensemble A、 B和C中的均方根误差(RMSE)和展形(SPREAD)变化过程;
图5是ensemble A和C中求得的参数均值场估计； 图6是五组初始Y场在无/有Y观测的情况下的十次无局域化集合同化中RMSE 和SPREAD变化过程； 图7是第一时段点(80m，720m)处Y值与整个水头场的交叉协方差(第一行)及 与整个Y场的协方差(第二行)。(a) (d)来自ensemble A ; (b) (e)来自ensemble D ; (c) (f)来自ensemble C ; 图8是五组初始Y场在无/有Y观测的情况下的十次有局域化集合同化中RMSE 和SPREAD变化过程。
具体实施例方式1、本发明的具体步骤如下 步骤1针对实地水文地质条件建立概念模型和数值模型，充分考虑边界条件和研 究区地质条件，同时确认地质统计学对研究区内待估参数的适用性，即满足二阶平稳假设 和本征假设。 步骤2对研究区应用地质统计学方法估计参数场的协方差结构C(h)和协方差参 数Ax，、及、作为地质统计先验信息。 常用的协方差类型主要有球形模型、指数模型、高斯模型、幂函数模型和迪维生模 型。目前随着野外水文地质实验效率的提高、新的地球物理探测等方法的引入，参数场的协 方差先验估计正趋于简便和准确。 步骤3根据参数的地质统计先验信息，生成一组Ne个参数场的随机实现，作为取 样集合(即参数场的随机模拟)。生成随机实现的途径也比较多，如GSLIB工具包、卡亨 南-刘维展开、傅里叶展开等等。由于模型复杂度和计算能力的限制，通常Ne的取值非常 有限( 一般不超过103)，而状态估计涉及的变量维数有远大于此(105 107)，这种差异是 造成估计矩阵病态的根本原因，源于此的取样噪声正是本方法要处理的问题。
步骤4将取样集合中的Ne个实现分别用于先前建立的数值模型，进行短时间上的 演化计算，即获得Ne个模型预报。用公式表达如下
<formula>formula see original document page 5</formula> 其中S是状态向量，包含模型参数，变量及其他需要的观测数据；t表示时间；F是
预报算子，代表所建立的模型的作用；上标'f'和'u'分别代表预报和更新过程；下标'r' 是集合成员序数。 步骤5统计这些实现的预报结果，根据集合预报原理估计模型中状态量(包括所 有参数及变量)之间的均值和协方差，即一阶矩和二阶矩，分别由下面两式计算
<formula>formula see original document page 5</formula>
<formula>formula see original document page 5</formula>
Cs是估计的状态误差协方差矩阵；CME是模型误差矩阵；上标'f'代表预报新过程； 上标'T'表示矩阵转置(下同)。
步骤6考虑参数和变量的空间位置关系，并将这种关系以权重P的舒尔积形式添
加到集合估计的协方差中，得到修正的协方差
<formula>formula see original document page 6</formula> 由于这种修正使协方差也具有了空间局部性，因此称为局域化修正。 步骤7将修改后的协方差用于卡尔曼滤波计算，输入当前可用的新观测数据d。bs,
r(t)，基于最小二乘原理获得状态(包括参数和变量)的优化估计，以此作为集合的更新后状态。
<formula>formula see original document page 6</formula>
d。bs,r是观测向量；H是观测算子；C。E是观测误差矩阵； 步骤8以更新的集合状态为初始状态进行下一时段计算，重复步骤4) 7)，直至 同化时段全部完成。 2、以下是采用本发明对一个理想算例的检验实验，据以说明本方法的有效性。
为了说明修正的EnKF估计渗透系数场的能力和研究局域化的影B向，这里构造 了一个理想二维平面矩形承压含水层模拟计算区长X宽为1440mX800m，均匀离散为 72X40个正方形单元(见图l)，其中Y观测位于实心方形所示的15个点；水头观测位于实 心和空心方形所示的47个点；实心圆为抽水井，实心三角为注水井。上下边界为隔水边界， 左右边界分别是水头为203m和197m的定水头边界。含水层顶板标高设为Om，厚度为lm。 流量均为150m3/day的抽水井和注水井分别位于点(240m， 560m)和点(1200m， 240m)。两井 从t = Oday时刻起开始运行并分别保持抽、注水流量恒定。 假设计算区初始水头为197m ;贮水系数恒为0. 0001 ;对数渗透系数场Y(x)= lnKs(x)为零均值、单位方差的高斯随机函数，二阶平稳并可以用以下可分离指数型协方差
函数刻画<formula>formula see original document page 6</formula>
其中<formula>formula see original document page 6</formula>T是x-y平面内的二维分离向量；、和、是沿x和y方向 的积分步长。这里A x = 200m， A y = 100m。因此饱和渗透系数的几何均值为Ke = 1. Om/ day。用Karhunen-Loeve (K-L)展开产生一组符合上述统计特性的对数渗透系数场，从中任 取一个作为真实参考场(见图2)。按图1所示位置选取区域内15个点作为渗透系数测量 点，认为这些位置上的Y值已知但带有均值为零、标准差为0. 001的正态分布的随机观测误 差。此外，每个时段有47个水头值作为水头观测数据(位置如图l所示)。假设水头观测 误差也是正态分布的，均值为零，标准差为0. 05m。同时假设所有观测数据相互无关，则C。E 是对角矩阵。本研究中所用的初始随机Y场同样使用K-L展开法产生，但取Ax = 220m， 入y二120m及oY=1.2。模拟总时间为9.0day，平均分为30个时段，步长0.3day。 47个 水头观测点的观测值每时段都向EnKF输送，而观测数据只在第一时段进行同化。同化过程中所有水流模型均使用相同代码(M0DFL0W-2000)进行模拟，故同化过程中无模型误差。
用一个二维可分离指数型协方差函数作为局域化函数(以下简称"L函数")，
p (h' ) = exp(_3| Ax|/9 x_3| Ay|/9 y) ， 9 x > 0 9y>0 其中e x和e y是表示滤波长度的参数，建议取为 《(A^)《A(V9 + 87^-5)/4
《(Ag A(V^^—5)/4 其中13 x和13 y是x和y方向的执行范围；Ne同前，表示集合的大小。
考虑到水头观测的影响半径比Y观测更大[3]，假设
pw(h' ) = [phY(h' )]2=[ptt(h' )]4 其中p YY， P hY和P tt分别是Y与Y之间、水头与Y之间及水头与水头之间的L函 数值。同时取水头间的Px = 200m， |3y = 100m并根据式(8)的假设计算其他的执行范围。 集合实现数为100，则有Ne = 100， e x = 1172m， e y = 586m。
2. 1结果与讨论 2. 1. 1局域化对EnKF参数估计能力的提高 为了更好地理解取样噪声的影响及局域化的机理，在图3中列出了第一时段末真 实等水头线(图3(a))以及从ensemble A、B和C中估算出的第一时段上点(80m， 720m)处 水头值与整个水头场的协方差(图3(b)、(c)、(d))。 ensemble A是大小为100的局域化集 合；ensemble B是大小为100的无局域化集合；ensemble C是大小为2000的无局域化集 合。实际上ensemble A与ensemble B中的初始Y场(即初始对数渗透系数场)取样完全 相同，局域化是它们的唯一差别。尽管ensemble C无法完全避免取样噪声，但在此认为它 提供的估计协方差相对准确，将其作为评价ensemble A和B的参考。
从ensemble C中可以看出点(80m， 720m)处(空心方块所示)的水头对附近一窄 带状分布的区域非常敏感。该窄带分布的区域沿定水头边界延伸且在抽水井附近有明显突 出。考虑到初始水头和边界条件的影响，除了右下角注水井附近有极小的协方差异常外，整 个结果是合理的。以此为标准，在ensemble B给出的协方差中(图3 (c))，敏感区域看起 来分布更广，特别是多出了注水井附近的大片高敏感区域。这正是伪相关造成的，因为在抽 水开始仅O. 3day的时候，降落漏斗仍局限在计算区域的左半区(见图3(a)等水位线)，注 水井附近的水头不可能跳过抽水井与更远距离外点(80m，720m)处的水头发生如此明显的 相关。更为重要的是，集合大小从100增至2000可以使这部分敏感区域基本消失(如图 3(d))，说明ensemble C增加实现数很大程度地去除了远距离相关性的噪声。
另一方面，ensemble A中的协方差敏感区域分布(图3(b))与ensemble C的非 常相似，只是前者比后者略小。这充分说明L函数也除去了远距离协方差的噪声。
应该注意的是，滤波长度太短的话可能会造成ensemble A估计的敏感区域过小， 制约水头观测数据的使用。合理的敏感区域(或相关性)估计意味着集合修正比较恰当， 以及滤波近交不严重，这是滤波过程成功的基础。 为了评价估计结果，我们引入常用的RMSE和SPREAD作为度量标准。RMSE(全称 Root Mean Square Error)代表估计量与参考场的接近程度 RMSE = ^!i(S' -Su)2//V
7
由于只对参数估计进行评价，因此只考虑Y均值场。Su为集合更新后的Y均值场； St为真实Y场；N是每个场中Y参数个数。而SPREAD (全称Ensemble Spread)可以表示为
SPREAD 一^iVAR五"(/)/W 其中VAREn (i)代表每个Y参数的集合方差。SPREAD是集合包含的变化性(或不 确定性)的直接度量。若集合对状态进行了恰当的估计，则SPREAD应该与RMSE相吻合。 ensemble A、 B禾P C同化过程的RMSE与SPREAD如图4所示。 从图4中可以明显看出ensemble B中的RMSE基本上随时间增加最后稳定于1. 6 左右，说明集合实现未能向参考场靠拢，完全偏离了实际，滤波发散了 ；而相应的SPREAD则 快速降到很低的水平，表明ensemble B未能保持变化性(滤波近交问题严重)，也未能从观 测中有效吸收有用信息，最终导致滤波的发散。相反，ensemble C中的RMSE与SPREAD稳定 且持续地下降。在第30时段，ensemble C中的RMSE与SPREAD分别降为0. 558和0. 555并 且相互吻合得极好，说明该集合运行正常，原因在于大集合取样时取样误差更小，滤波近交 问题也弱得多。ensemble A虽然只有100个实现，但其结果也相当令人满意RMSE和SPREAD 都与ensemble C的结果接近，而且彼此也吻合得很好(0. 605和0. 580) 。 ensemble A和C 的参数估计结果相差不大，均值场估计也与参考场相似(见图5)。值得注意的是，ensemble A的两个指标下降速度都比ensemble C的慢，而且从始至终高于后者，原因在于L函数过滤 取样噪声的同时也造成少量有用信息的丢失。 简言之，ensemble A依赖远小于ensemble C的实现数量(100 : 2000)却获得了
几乎同样准确的估计结果；ensemble B与ensemble A使用相同的初始Y场，但由于取样噪
声的影响完全偏离了真实场(即滤波发散)。说明对协方差估计进行局域化可以大大提高
EnKF的参数估计能力。 2. 1. 2Y观测数据对局域化的依赖 EnKF的一大优势是仅仅考虑数据间的相关性，因此可以同化来自不同来源的不同 类型信息。但也正是由于EnKF仅仅考虑相关性，相关性估计中的误差会对相关长度较小的 信息(如渗透系数相对于水头)造成更加明显的影响。换言之，一种观测数据，即使精度很 高，若其相关长度很小，则稍远距离的相关性估计就容易带有很大的相对误差，因此只在很 短的距离内才有同化的价值。正如Chen和Zhang (CHEN Y，ZHANG D. Dataassimilation for transient flow in geologic formations via ensemble KalmanfiIter[J]. Advances in Water Resources, 2006， 29 (8) :1107-1122.)所指出的，Y观测由于本身的相关长度较小， 只在一定半径范围内有贡献。他们在一个41X41X1的网格模型中添加了 9个Y观测数 据，使参数估计结果有了轻微改善。模拟计算实验表明，这种"改善"由于与取样噪声和取 样代表性有关，因此并不稳定。本节用五组各500个实现的不同初始Y场集合(记为setA、 B、 C、 D和E)分别在有Y观测数据和没有Y观测数据的情况下进行EnKF同化，共进行十次 参数估计过程，RMSE和SPREAD结果如图6所示。不同的估计过程SPREAD差别微弱而RMSE 则显著不同。比较各组集合有/无Y观测数据的过程，可以发现添加额外的Y观测数据使 RMSE在同化初期快速下降；SPREAD也有明显降低。与Chen和Zhang的结论不同点在于最 终的RMSE :set C和E中添加Y观测数据后RMSE降低；但set A和B有Y观测数据的RMSE 反而比没有Y观测数据的RMSE更高；而set D则几乎不变。这说明增加加Y观测数据给估 计过程造成的影响是不稳定的，RMSE可能更高也可能更低。另一方面，增加的Y观测数据
8对SPREAD的影响十分有限，只是使其下降速度在同化过程中略有加快；这种影响也随时间逐渐减弱。 从set A、 B、 C、 D禾P E五个Y集合中任取一个作为ensemble D (操作时选取的是set A)。为了直观地比较ensemble A、 D和C中与Y观测数据相关的取样噪声，在此将点(80m，720m)处的Y值与状态向量(包括水头和Y)的协方差表示于图7。从其中ensembleD和C的结果可以看出，点(80m，720m)处的Y观测数据与整个Y场的估计协方差在很大的区域内远不为零。虽然Y与水头的相关性可因水流环境而异，但是Y参考场设定的相关长度最大不过200m，超过此距离的Y与Y的相关性应趋于零。换句话说，ensembleD估计出的Y数据相关性比实际存在的相关性强得多(甚至ensemble C也是如此)，在Y数据之间尤其明显；正是由于Y数据相关长度小，随距离增加真实相关性快速减弱，所以更容易被取样噪声干扰。相关性估计中存在的这种显著偏差，会造成输入的Y观测数据修改原本几乎无关的变量，导致估计结果的歪曲。而一般情况下同化初期的集合方差明显大于后期的集合方差，前期观测数据获得的权重会远大于后期观测获得的权重，因此初期的同化过程(尤其是第一、二步)常有决定性的影响，被不恰当修改的实现在后续过程中也很难得到修正，常会造成意料之外的问题(例如Y值超出计算程序的合理存储范围导致的水流模拟程序运行失败)。Y观测数据对噪声的相对易感性使得将Y观测值的影响限制在观测点的周围变得更加必要。 为了进一步进行说明，这里用五组100个实现的带局域化集合(记为set a、b、c、d和e)在无Y观测数据和有Y观测数据的情况下进行了总共十个同化过程实验，RMSE和SREAD如图8所示。与无局域化的EnKF不同，带局域化的EnKF可以很好地整合增加的Y观测数据，在5组不同集合下均能保证相应的RMSE和SPREAD低于未增加Y观测数据的同化过程。可以看出RMSE的改善不明显，与Chen和Zhang的结论是一致的，不过局域化使得这种改善更加稳定(如图8(a)所示)；SPREAD也只有微小的降低。最终RMSE和SPREAD都吻合得非常好。 Y数据的相关长度比水头数据的更小，即前者的实际物理影响更局限。向系统输入
Y观测数据，对该观测点附近的变量是有积极影响的，但稍远的变量却可能因高估的相关性
在状态更新步骤中不恰当地被这个观测数据修改，反而偏离真实场，给全局带来消极影响。
因此这种观测数据的添加更依赖局域化对取样噪声的过滤作用。 2. 2结论 本发明利用协方差局域化方法过滤与小集合取样有关的伪相关，大大改善了EnKF估计渗透系数场的执行效果。 (1)带局域化的EnKF可以从同样大小的集合中获得更为准确的参数估计，小集合中尤其明显。同时局域化极大地减轻了对集合方差的系统性低估，从而有效避免滤波近交和滤波发散。带局域化的EnKF可以显著提高水文地质参数估计能力。
(2)进一步研究表明，相关长度更小的Y数据更容易受取样噪声干扰；由于噪声使协方差估计中带有伪相关，Y观测数据在改善观测点附近变量估计的同时有可能歪曲距离较远处的变量估计，给同化系统带来负面影响。局域化的EnKF中，Y观测数据的影响被限制在相关长度的范围内，因此可以很大程度上避免对全局的歪曲。这种处理使Y观测数据能够稳定地改善渗透系数场的估计。
9
权利要求
一种水文地质反求参数的数据同化方法，其特征在于包括以下步骤1)针对实地水文地质条件建立概念模型和数值模型，充分考虑地质条件和边界条件，确认研究区地质统计学的适用性，即满足二阶平稳假设和本征假设；2)对研究区应用地质统计学方法估计参数场的协方差结构C(h)和协方差参数λx，λy及λz作为地质统计先验信息；3)根据参数的地质统计先验信息，生成一组Ne个参数场的随机实现，作为取样集合，即参数场的随机模拟；4)将取样集合中的Ne个实现分别用于先前建立的数值模型，进行短时间上的演化计算，即获得Ne个模型预报；用公式表达如下 <mrow><msubsup> <mi>S</mi> <mi>r</mi> <mi>f</mi></msubsup><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow><mo>=</mo><mi>F</mi><mo>[</mo><msubsup> <mi>S</mi> <mi>r</mi> <mi>u</mi></msubsup><mrow> <mo>(</mo> <mi>t</mi> <mo>-</mo> <mi>Δt</mi> <mo>)</mo></mrow><mo>]</mo><mo>+</mo><msub> <mi>e</mi> <mrow><mi>ME</mi><mo>,</mo><mi>r</mi> </mrow></msub><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow> </mrow>其中S是状态向量，包含模型参数，变量及其他需要的观测数据；t表示时间；F是预报算子，代表所建立的模型的作用；上标‘f’和‘u’分别代表预报和更新过程；下标‘r’是集合成员序数；5)统计这些实现的预报结果，根据集合预报原理估计模型中所有参数及变量的均值和协方差，即一阶矩和二阶矩，分别由下面两式计算 <mrow><mo>⟨</mo><msup> <mi>S</mi> <mi>f</mi></msup><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow><mo>⟩</mo><mo>≈</mo><mfrac> <mn>1</mn> <msub><mi>N</mi><mi>e</mi> </msub></mfrac><munderover> <mi>Σ</mi> <mrow><mi>r</mi><mo>=</mo><mn>1</mn> </mrow> <msub><mi>N</mi><mi>e</mi> </msub></munderover><msubsup> <mi>S</mi> <mi>r</mi> <mi>f</mi></msubsup><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow> </mrow> <mrow><mrow> <msup><msub> <mi>C</mi> <mi>S</mi></msub><mi>f</mi> </msup> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow> <mo>≈</mo> <mfrac><mn>1</mn><mrow> <msub><mi>N</mi><mi>e</mi> </msub> <mo>-</mo> <mn>1</mn></mrow> </mfrac> <munderover><mi>Σ</mi><mrow> <mi>r</mi> <mo>=</mo> <mn>1</mn></mrow><msub> <mi>N</mi> <mi>e</mi></msub> </munderover> <mo>{</mo> <mo>[</mo> <msubsup><mi>S</mi><mi>r</mi><mi>f</mi> </msubsup> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow> <mo>-</mo> <mo>⟨</mo> <msup><mi>S</mi><mi>f</mi> </msup> <mrow><mo>(</mo><mi>t</mi><mo>)</mo> </mrow> <mo>⟩</mo> <mo>]</mo></mrow><msup> <mrow><mo>[</mo><msubsup> <mi>S</mi> <mi>r</mi> <mi>f</mi></msubsup><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow><mo>-</mo><mo>⟨</mo><msup> <mi>S</mi> <mi>f</mi></msup><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow><mo>⟩</mo><mo>]</mo> </mrow> <mi>T</mi></msup><mo>}</mo><mo>+</mo><msub> <mi>C</mi> <mi>ME</mi></msub><mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo></mrow> </mrow>CS是估计的状态误差协方差矩阵；CME是模型误差矩阵；上标‘f’代表预报新过程；上标‘T’表示矩阵转置；6)考虑参数和变量的空间位置关系，并将这种关系以舒尔积权重的形式添加到集合估计的协方差中得到修正的协方差7)将修改后的协方差用于卡尔曼滤波计算，输入当前可用的新观测数据，基于最小二乘原理获得状态，包括参数和变量，的优化估计，更新集合状态；8)以更新的集合状态为初始状态，进行下一时间段的预报，即重复4)~7)，直至同化时间全部完成。F2010100183606C00014.tif
全文摘要
本发明公开了一种水文地质参数的估计方法，其将经典卡尔曼滤波在已有的集合估计的基础上添加协方差局域化系统形成了一个新方法。本方法考虑了复杂水文地质条件下多类型、多来源测量数据之间的相关性；针对集合估计理论的不足，采用舒尔积权重最大程度地降低有限集合中的取样噪声；具有连续同化的特点，可以在观测的同时进行求参，随着获得的观测数据不断的输入，求参效果也会不断改善；获得参数估计的同时也得到了对应估计所包含的不确定性，更加科学、准确。本方法有机集成了多个理论和方法的长处，且使用方便，有着同类复杂水文地质求参问题下的普遍适用性。
文档编号G01V9/00GK101793977SQ201010018360
公开日2010年8月4日 申请日期2010年1月14日 优先权日2010年1月14日
发明者南统超, 吴吉春, 施小清 申请人:南京大学