一种基于多源观测数据的水质模型粒子滤波同化方法技术

本发明专利技术提供了一种基于多源观测数据的水质模型粒子滤波同化方法，构建二维水质模型；初始化粒子的状态变量和参数；生成粒子的边界条件；更新粒子的状态变量；判断是否有水质观测值；重采样获取新的粒子集合；计算二维水质模型的模拟状态变量和参数的最优估值；将粒子的参数从t时刻递推到t+1时刻；更新时刻，继续生成粒子的边界条件，直至所有时刻运行完成，实现对二维水质模型的粒子滤波同化。本发明专利技术采用基于自适应网格和并行计算的二维水质模型，考虑了二维水质模型模拟关键参数的时空变异性，利用粒子滤波算法将水质多源观测数据合理地融入二维水质模型，动态更新二维水质模型参数，提高了二维水质模型的模拟精度和预测能力。

A Particle Filter Assimilation Method for Water Quality Model Based on Multi-source Observation Data

【技术实现步骤摘要】
一种基于多源观测数据的水质模型粒子滤波同化方法
本专利技术属于水环境
，尤其涉及一种基于多源观测数据的水质模型粒子滤波同化方法。
技术介绍
面对日趋严重的水环境问题，观测分析和模型模拟是水环境管理中的重要手段。观测分析手段主要包括实地采样实验室分析手段和遥感监测手段，实地采样实验室分析手段通过布设断面采集水样，然后通过实验室分析获取水质状况，这种监测手段虽然监测精度高，但耗时耗力，且只能获取点上的水质状况，监测频次有限，难以满足水环境动态管理的需求；遥感监测能够快速获取湖泊水质的时空分布，越来越多地应用于湖泊水环境的监测和管理。水动力水质模型是水环境管理的重要工具，可以获取时空连续的水质分布，但受模型参数、模型输入和模型结构等不确定性因素的影响，模拟精度有待提高。数据同化可以实现观测和模型模拟优势互补，利用同化算法可以将多源观测数据合理地融入水质模型，校正模型模拟结果，同步更新模型参数，提升模型模拟精度和预测能力。粒子滤波在非线性非高斯模型数据同化中表现出色，逐渐应用于一维水质模型数据同化。现有水质模型粒子滤波数据同化研究中多采用空间均一的模型参数，忽略模型参数的空间变异性，难以同时得到不同空间(区域)水质模拟结果最优估计，无法适应二维水质模型数据同化；此外，水质模型粒子滤波数据同化中使用的多为水质原型观测数据，未能将水质原型观测和遥感观测的优势和水质模型模拟优势充分结合，限制了二维水质模型模拟和预测精度。
技术实现思路
针对现有技术中的上述不足，本专利技术提供的一种基于多源观测数据的水质模型粒子滤波同化方法解决了现有水质模型粒子滤波数据同化研究中多采用空间均一的模型参数，忽略模型参数的空间变异性，难以同时得到不同空间(区域)水质模拟结果最优估计，无法适应二维水质模型数据同化的问题，以及水质模型粒子滤波数据同化中不能将水质原型观测和遥感观测的优势和水质模型模拟优势充分结合，从而限制了二维水质模型模拟和预测精度的问题。为了达到以上目的，本专利技术采用的技术方案为：本方案提供一种基于多源观测数据的水质模型粒子滤波同化方法，包括如下步骤：S1、根据研究区的边界矢量、地形高程、水位、流量以及水质信息构建二维水质模型；S2、根据所述二维水质模型初始化t时刻粒子的水位、流量、水质和t+1时刻粒子的参数；S3、根据所述二维水质模型在t+1时刻的水位、流量和水质边界条件生成t+1时刻粒子的边界条件；S4、根据所述t时刻的粒子的水位、流量、水质和t+1时刻粒子的边界条件和参数驱动二维水质模型更新得到t+1时刻粒子状态变量；S5、判断t+1时刻是否有水质观测值，若有，则计算t+1时刻各粒子的似然函数值，更新粒子的权重，并进入步骤S6，反之，则维持粒子权重不变，得到的t+1时刻的原始粒子集合，并进入步骤S7；S6、根据所述粒子的权重对粒子进行多项式重采样，得到新的粒子集合；S7、根据所述新的粒子集合或t+1时刻的原始粒子集合分别计算所述二维水质模型的模拟水质状态变量和二维水质模型参数的最优估值；S8、令时刻t为t+1，并根据所述新的粒子集合的参数或原始粒子集合的参数利用核平滑方法将粒子的参数从t时刻递推到t+1时刻；S9、判断所有时刻是否运行完成，若否，则返回步骤S3，直至所有时刻运行完成，反之，则结束对粒子的滤波同化，从而实现对水质模型粒子的滤波同化。进一步地，所述步骤S1中构建二维水质模型的表达式如下：其中，U表示守恒向量，F和G分别x、y轴方向上的对流通量，分别表示U对时刻的偏导数、F对时刻的偏导数和G对时刻的偏导数，S表示源项，η表示水位，v和u分别表示x、y轴方向的水流流速，c表示物质的垂线平均浓度，g表示重力加速度，zb表示研究区的底部高程，qin和cin分别表示点源的流量强度和物质垂线平均浓度，τbx和τby分别表示x、y轴方向的床面摩擦应力，ρ表示水密度，分别表示zb对x、y轴的偏导数，Dx和Dy分别表示x、y轴方向的扩散系统，h表示水深，Sk表示水质浓度有关的生化反应项，和分别表示物质垂线浓度对x、y轴的偏导数，和分别表示对x、y轴的偏导数。再进一步地，所述步骤S2具体为：根据t时刻的二维水质模型在各个自适应网格处的状态变量和参数的先验分布采样生成N个等权重的粒子，并初始化各粒子的权重为从而实现对粒子的水位、流量、水质和t+1时刻粒子的参数的初始化。再进一步地，所述初始化粒子的水位、流量、水质和参数的表达式如下：其中，和分别表示t时刻第j个网格处第i个粒子、粒子模拟水位、流量、水质、模型参数和粒子权重，分别表示第j个网格处水位、流量、水质初始值和参数初始值，i表示粒子的编号，N表示粒子的总数，j表示网格单元的编号，ncell表示计算网格的总数，表示第j个网格处粒子参数的扰动误差，u表示状态变量和参数的均匀分布，n1,n2分别表示参数扰动误差取值的下限和上限，n表示水质模拟参数。再进一步地，所述步骤S3中生成粒子的边界条件的表达式如下：其中，表示t+1时刻边界条件观测值。再进一步地，所述步骤S4中更新粒子状态变量的表达式如下：其中，f(·)表示二维水质模型，Mt+1和Mt分别为t+1时刻二维水质模型水位和水质的状态变量，nMt+1和分别为t+1时刻二维水质模型的参数和边界条件，表示t时刻粒子的状态变量，表示t+1时刻粒子的参数，表示t+1时刻粒子的边界条件。再进一步地，所述步骤S5中计算t+1时刻各粒子的似然函数值的表达式如下：其中，表示t+1时刻第j个网格处第i个粒子权重，σo表示水质观测误差的标准差，exp(·)表示指数运算，表示t+1时刻第j个观测点处叶绿素a浓度观测值，表示t+1时刻第j个网格处第i个粒子的状态变量。再进一步地，所述步骤S6包括如下步骤：S601、利用多项式重采样算法将所述粒子的权重集合组成多项式分布其中，Mult(·)表示多项分布，表示表示t+1时刻第j个网格处第N个粒子权重，N表示粒子的总数；S602、随机生成Nr个服从(0,1]均匀分布的随机数rk～u(0,1]，其中，rk表示均匀分布的随机数，且k＝1,2...Nr，Nr＝N，N表示粒子的总数，u表示状态变量和参数的均匀分布；S603、计算各粒子的权重累积和序列，所述各粒子的权重累积和序列的表达式如下：其中，ci+1,j表示第j个网格处第i+1个粒子权重累加值，ci,j表示第j个网格处第i粒子个权重累加值，表示t+1时刻第j个网格处第i个粒子权重，i表示粒子的编号，j表示网格单元的编号，ncell表示计算网格的总数。S604、统计权重累积区间(ci,j,ci+1,j]中落入随机数rk的次数，并将所述落入随机数rk的次数记为Ni,j，得到粒子复制信息；S605、根据所述粒子复制信息依次对粒子的状态变量进行复制，从而得到重采样后权重均等新的粒子集合其中，分别表示重采样后得到新粒子的状态变量、参数和粒子权重。再进一步地，所述步骤S7中根据所述新的粒子集合计算所述二维水质模型的模拟水质状态变量的最优估值的表达式如下：其中，i表示粒子的编号，N表示粒子的总数，表示t+1时刻第j个网格处重采样得到的第i个粒子权重，表示t+1时刻第j个网格处重采样得到的第i个粒子的水质状态变量；根据所述新的粒子集合计算二维水质模型参数的最优估值的表达式本文档来自技高网...

【技术保护点】
1.一种基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，包括如下步骤：S1、根据研究区的边界矢量、地形高程、水位、流量以及水质信息构建二维水质模型；S2、根据所述二维水质模型初始化t时刻粒子的水位、流量、水质和t+1时刻粒子的参数；S3、根据所述二维水质模型在t+1时刻的水位、流量和水质边界条件生成t+1时刻粒子的边界条件；S4、根据所述t时刻的粒子的水位、流量、水质和t+1时刻粒子的边界条件和参数驱动二维水质模型更新得到t+1时刻粒子状态变量；S5、判断t+1时刻是否有水质观测值，若有，则计算t+1时刻各粒子的似然函数值，更新粒子的权重，并进入步骤S6，反之，则维持粒子权重不变，得到的t+1时刻的原始粒子集合，并进入步骤S7；S6、根据所述粒子的权重对粒子进行多项式重采样，得到新的粒子集合；S7、根据所述新的粒子集合或t+1时刻的原始粒子集合分别计算所述二维水质模型的模拟水质状态变量和二维水质模型参数的最优估值；S8、令时刻t为t+1，并根据所述新的粒子集合的参数或原始粒子集合的参数利用核平滑方法将粒子的参数从t时刻递推到t+1时刻；S9、判断所有时刻是否运行完成，若否，则返回步骤S3，直至所有时刻运行完成，反之，则结束对粒子的滤波同化，从而实现对水质模型粒子的滤波同化。...

【技术特征摘要】
1.一种基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，包括如下步骤：S1、根据研究区的边界矢量、地形高程、水位、流量以及水质信息构建二维水质模型；S2、根据所述二维水质模型初始化t时刻粒子的水位、流量、水质和t+1时刻粒子的参数；S3、根据所述二维水质模型在t+1时刻的水位、流量和水质边界条件生成t+1时刻粒子的边界条件；S4、根据所述t时刻的粒子的水位、流量、水质和t+1时刻粒子的边界条件和参数驱动二维水质模型更新得到t+1时刻粒子状态变量；S5、判断t+1时刻是否有水质观测值，若有，则计算t+1时刻各粒子的似然函数值，更新粒子的权重，并进入步骤S6，反之，则维持粒子权重不变，得到的t+1时刻的原始粒子集合，并进入步骤S7；S6、根据所述粒子的权重对粒子进行多项式重采样，得到新的粒子集合；S7、根据所述新的粒子集合或t+1时刻的原始粒子集合分别计算所述二维水质模型的模拟水质状态变量和二维水质模型参数的最优估值；S8、令时刻t为t+1，并根据所述新的粒子集合的参数或原始粒子集合的参数利用核平滑方法将粒子的参数从t时刻递推到t+1时刻；S9、判断所有时刻是否运行完成，若否，则返回步骤S3，直至所有时刻运行完成，反之，则结束对粒子的滤波同化，从而实现对水质模型粒子的滤波同化。2.根据权利要求1所述的基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，所述步骤S1中构建二维水质模型的表达式如下：其中，U表示守恒向量，F和G分别x、y轴方向上的对流通量，分别表示U对时刻的偏导数、F对时刻的偏导数和G对时刻的偏导数，S表示源项，η表示水位，v和u分别表示x、y轴方向的水流流速，c表示物质的垂线平均浓度，g表示重力加速度，zb表示研究区的底部高程，qin和cin分别表示点源的流量强度和物质垂线平均浓度，τbx和τby分别表示x、y轴方向的床面摩擦应力，ρ表示水密度，分别表示zb对x、y轴的偏导数，Dx和Dy分别表示x、y轴方向的扩散系统，h表示水深，Sk表示水质浓度有关的生化反应项，和分别表示物质垂线浓度对x、y轴的偏导数，和分别表示对x、y轴的偏导数。3.根据权利要求1所述的基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，所述步骤S2具体为：根据t时刻的二维水质模型在各个自适应网格处的状态变量和参数的先验分布采样生成N个等权重的粒子，并初始化各粒子的权重为从而实现对粒子的水位、流量、水质和t+1时刻粒子的参数的初始化。4.根据权利要求3所述的基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，所述初始化粒子的水位、流量、水质和参数的表达式如下：其中，和分别表示t时刻第j个网格处第i个粒子、粒子模拟水位、流量、水质、模型参数和粒子权重，分别表示第j个网格处水位、流量、水质初始值和参数初始值，i表示粒子的编号，N表示粒子的总数，j表示网格单元的编号，ncell表示计算网格的总数，表示第j个网格处粒子参数的扰动误差，u表示状态变量和参数的均匀分布，n1,n2分别表示参数扰动误差取值的下限和上限，n表示水质模拟参数。5.根据权利要求1所述的基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，所述步骤S3中生成粒子的边界条件的表达式如下：其中，表示t+1时刻边界条件观测值。6.根据权利要求1所述的基于多源观测数据的水质模型粒子滤波同化方法，其特征在于，所述步骤S4中更新粒子状态变量的表达式如下：其中，f(·)表示二维水质模型，Mt+1和Mt分别为t+1时刻二维水质模型水位和水质的状态变量，nMt+1和分别为t+1时刻二维水质模型的参数和边界条件，表示t时刻粒子的状态变量，表示t+1时刻粒子的参数，表示t+1...

【专利技术属性】
技术研发人员：冶运涛，曹引，蒋云钟，梁犁丽，顾晶晶，赵红莉，尚毅梓，龚家国，张双虎，张海涛，杜军凯，段浩，
申请(专利权)人：中国水利水电科学研究院，
类型：发明
国别省市：北京,11