CN108334660B - 一种基于数据同化的强冲积河流的水沙预测方法和系统 - Google Patents

一种基于数据同化的强冲积河流的水沙预测方法和系统 Download PDF

Info

Publication number
CN108334660B
CN108334660B CN201711461256.2A CN201711461256A CN108334660B CN 108334660 B CN108334660 B CN 108334660B CN 201711461256 A CN201711461256 A CN 201711461256A CN 108334660 B CN108334660 B CN 108334660B
Authority
CN
China
Prior art keywords
water
sand
grid
river
strong
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711461256.2A
Other languages
English (en)
Other versions
CN108334660A (zh
Inventor
冶运涛
曹引
尚毅梓
梁犁丽
蒋云钟
龚家国
仇亚琴
杜龙江
王浩
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Institute of Water Resources and Hydropower Research
Original Assignee
China Institute of Water Resources and Hydropower Research
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Institute of Water Resources and Hydropower Research filed Critical China Institute of Water Resources and Hydropower Research
Priority to CN201711461256.2A priority Critical patent/CN108334660B/zh
Publication of CN108334660A publication Critical patent/CN108334660A/zh
Application granted granted Critical
Publication of CN108334660B publication Critical patent/CN108334660B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Strategic Management (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Development Economics (AREA)
  • Geometry (AREA)
  • Game Theory and Decision Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明实施例提供了一种基于数据同化的强冲积河流的水沙预测方法和系统，具体为利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元；获取强冲积河流进出口边界的水沙信息和初始场信息；以地形网格的三角形单元为控制体，将强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到控制体的计算值；实时采集强冲积河流的河道水沙状态信息，根据水沙状态信息建立基于数据同化的实时预测模型，得到同化状态变量和参数变量；最后利用二维水沙数据同化模型对强冲积河流进行水沙预测。通过上述处理可以对强冲积河流的水沙实时状态进行有效预测，从而可以支持水沙调控，并为防洪减灾、水环境保护和水资源管理提供依据。

Description

一种基于数据同化的强冲积河流的水沙预测方法和系统
技术领域
本发明涉及水利工程技术领域，特别是涉及一种基于数据同化的强冲积河流的水沙预测方法和系统。
背景技术
随着全球气候变化，流域下垫面随之发生相应的变化，且随着大坝修建等人类活动也改变了天然的来水来沙情势和水沙输移边界条件，由此可能会引发溃坝等极端事件的发生，这会进一步导致活跃的水沙运动和快速的河床变形，从而给科学认知冲积河流过程水沙相互作用机理带来了困难，给预测水沙运动变化带来了更大的挑战。
水沙模型参数受人类对泥沙输移认知的局限和模型次网格效应的影响，在实际操作过程中，模型参数确定不易，难以通过理论推导获得，一般基于原型观测和实验等直接测量数据确定，或者是通过历史数据调试或率定以适应特定的研究区域，这就意味着泥沙模型通常不是基于物理机制而是通过经验建立，从而导致水沙模型参数不易被直接测量。因此，水沙模型发展中的关键问题是如何估计水沙模型参数和状态变量，降低模型计算的不确定性，以便实现对水沙模型状态变量的准确模拟预测。
传统的预测系统主要采用历史回归的方法，如神经网络、支持向量机等对河道的水沙现状进行预测。这些方法忽略了河道水流演进中动力学特性和守恒规律，从而导致利用模型进行实时水位、流量、水质等方面预报时出现了误差大、精度低的问题，因此此类系统只适合于瞬时预测，无法辅助决策，达不到水资源实时管理的目的。
发明内容
有鉴于此，本发明提供了一种基于数据同化的强冲积河流的水沙预测方法和系统，以预测强冲积河流实时的水沙状态，以支持水沙调控，并为防洪减灾、水环境保护和水资源管理提供依据。
为了解决上述问题，本发明公开了一种基于数据同化的强冲积河流的水沙预测方法，包括步骤：
利用采样的地形高程散点数据，生成不规则三角形地形网格单元；
获取强冲积河流的进出口边界的水沙信息，并设置模型的初始场信息；
以所述地形网格的三角形单元为控制体，将所述强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到所述控制体的计算值，作为模型预测场信息；
实时采集所述强冲积河流的河道水沙状态信息，根据所述水沙状态信息建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量，作为模型分析场信息；
根据预设的强冲积河流水沙模型的计算间隔，将步骤所述同化状态变量和所述参数变量的分析场作为各所述三角形单元的状态变量初始场和参数变量初始场，利用所述水沙信息和所述初始场信息以及所述二维水沙模型对所述强冲积河流进行水沙预测。
可选的，所述利用采样的地形高程散点数据，生成不规则三角形地形网格单元，包括：
利用采样得到的地形高程散点数据生成不规则三角形地形网格单元；
标识所述地形网格的三角形单元及其边、节点之间的拓扑关系；
计算所述三角形单元的几何特征；
赋予所述三角形单元、所述边和所述节点相对应的水沙状态变量和参数变量。
可选的，所述水沙信息包括河道进出口边界的流量边界、水位边界、含沙量边界、河流悬沙级配组成、河流床沙级配组成、泥沙干容重和泥沙湿容重中的部分或全部。
可选的，所述计算值包括水位计算值、流速计算值、含沙量计算值和河床高程计算值中的部分或全部。
可选的，所述水沙状态信息包括水位测量值、流量或流速测量值和含沙量测量值的部分或全部。
相应的，为了保证上述方法的实施，本发明还提供了一种基于数据同化的强冲积河流的水沙预测系统，包括：
网格单元生成模块，用于利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元；
信息获取模块，用于获取强冲积河流的进出口边界水沙信息和初始场信息；
模型计算模块，用于以所述地形网格的三角形单元为控制体，将所述强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到所述控制体的计算值；
信息采集模块，用于实时采集所述强冲积河流的河道水沙状态信息，根据所述水沙状态信息建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量；
预测输出模块，用于根据预设的强冲积河流水沙模型的计算间隔，将步骤所述同化状态变量和所述参数变量作为各所述三角形单元的状态变量初始场和参数变量初始场，利用所述水沙信息和所述初始场信息以及所述二维水沙模型对所述强冲积河流进行水沙预测。
可选的，所述网格单元生成模块包括：
生成处理单元，用于利用采样的地形高程散点数据生成不规则三角形地形网格单元；
标识处理单元，用于标识所述地形网格的三角形单元及其边、节点之间的拓扑关系；
计算处理单元，用于计算所述三角形单元的几何特征；
赋值处理单元，用于赋予所述三角形单元、所述边和所述节点相对应的水沙状态变量和参数变量。
可选的，所述水沙信息包括河道进出口边界的流量边界、水位边界、含沙量边界、河流悬沙级配组成、河流床沙级配组成、泥沙干容重和泥沙湿容重中的部分或全部。
可选的，所述计算值包括水位计算值、流速计算值、含沙量计算值和河床高程计算值中的部分或全部。
可选的，所述水沙状态信息包括水位测量值、流量或流速测量值和含沙量测量值的部分或全部。
从上述技术方案可以看出，本发明提供了一种基于数据同化的强冲积河流的水沙预测方法和系统，具体为利用采样得到的地形高程散点数据，生成不规则三角形网格单元；获取强冲积河流的进出口边界水沙信息和初始场信息；以三角形地形网格单元为控制体，将强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到控制体单元的计算值；实时采集强冲积河流的河道水沙状态信息，根据水沙状态信息建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量；根据预设的强冲积河流水沙模型的计算间隔，将同化状态变量和参数变量作为各三角形单元的状态变量初始场和参数变量初始场，利用水沙信息和初始场信息以及二维水沙模型对强冲积河流进行水沙预测。通过上述处理可以对强冲积河流水沙的实时状态进行有效预测，从而可以支持水沙调控，并为防洪减灾、水环境保护和水资源管理提供依据。
并且，还能针对实际工程需要，结合已有的实时水沙观测手段，构建水沙实时数据接收系统，实时获取水位、流速或流量和含沙量等数据，将实测的水位、流速或流量和含沙量等数据考虑在模型中，使得原有的水沙数学模型的应用范围从工程设计和规划领域拓展到河道水沙的实时预报领域。
另外，本发明采用先进的数据同化方法，同时考虑实测数据的测量误差和水沙模型的计算误差，对实测数据融入水沙模型的过程进行优化，使得模型预报的初始场取得最优值，从而有效提高水沙模型实时预报河道水沙的预报精度。
再有，本发明方法系统地提出了从实时水沙数据接收、水沙模型演算、数据同化初始场、水沙状态量预报等全部模块，完善了水沙模型实时预报系统框架，可实现水位流量预测、泥沙预测和洪水预报等功能，具有针对性强、功能齐全、方便实用等特点，可应用在大江大河的河道洪水预报中，为实际防汛指挥工作提供决策支持。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种基于数据同化的强冲积河流的水沙预测方法的步骤流程图；
图2为本发明实施例提供的一种不规则三角形网格单元的示意图；
图3为本发明实施例提供的一种基于数据同化的强冲积河流的水沙预测系统的结构框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例一
图1为本发明实施例提供的一种基于数据同化的强冲积河流的水沙预测方法的步骤流程图。
参照图1所示，本实施例提供的水沙预测方法用于基于数据同化的方法对强冲积河流的水沙情况进行预测，具体的水沙预测方法包括如下步骤：
S101：利用采样得到的地形高程散点数据生成不规则三角形地形网格单元。
地形网格中会包括相应的三角形单元，该三角形单元具体可以作为强冲积河流水沙模型离散的控制体。具体的生成步骤为：首先利用采样的地形高程散点数据，生成不规则三角形地形网格单元；然后，对于不规则三角形地形网格，标识三角形单元、边、节点之间的拓扑关系；再后，计算三角形单元的几何特征；最后赋予三角形单元、边和节点相应的水沙状态变量和参数变量数据，并设计数据结构进行存储。
对于每一三角形、边和节点都对应一个记录，三角形的记录包括指向3个边的记录的指针。边的记录有4个指针字段，包括两个指向相邻三角形的指针和它的两个顶点指针。每个节点包括三个坐标值的字段，分别存储X、Y、Z坐标，这种拓扑网络结构的特点对于三角形的每一条边都记录了与三角形的拓扑邻接关系。
一个完整的构网原始数据，即用来构网的离散点可用如下结构来描述：
Figure BDA0001530303620000061
每个构网过程中新生成的三角形可用如下结果来描述:
Figure BDA0001530303620000062
TRIANGLE对应着三角网中的每一个三角形，数组vertex描述了构成三角形的三个离散点的序号和顺序，而数组neighbor记录了三角形三条边依次的邻接三角形序号。
建模优化算法中，首先定义一个结点队列来描述原始数据：CArray<CTinPoint*,CTinPoint*>m_aPoints。每个离散点都用CTinPoint数据结构来表示。
定义一个三角形队列来记录构建的三角形：CArray<CTriangle*,CTriangl-e*>m_aTriangles。每个三角形都用CTriangle数据结构来表示。
参见图2所示，变量L表示优化算法生成的三角形总数量，每扩展生成一个三角形，L就增加一个数值。变量k表示当前扩展三角形在整个三角形队列中的序号，即当前正在处理第几号三角形。每完成一个三角形的扩展，变量k就增加一个数值，即开始扩展下一个三角形。变量i表示扩展三角形正在处理的边。函数FindFirstTriangle用来确定第一个三角形，并以此三角形为基础进行三角网扩展。函数FindNeighborTriangle用来确定与扩展三角形，最终建立三角形的离散点，并记录新生成的三角形信息。因为在扩展新三角形时要以数据预处理阶段得到的山谷线、山脊线等地性线作为建立三角网的控制条件，即在构建三角网时，所建立的三角网不能跨越地性线，也就是要保证地性线无条件地作为三角形的边。
因此，地形的三维重建问题可以看作一个如何将离散点构建成约束Dela-unay三角网的问题。该优化算法的主要实现步骤如下。
(1)根据本优化算法对数据的要求进行数据预处理，生成离散点集合；
(2)在离散点中找出符合条件的3个点生成第一个三角形，编号为1，并加入到三角形队列中；
(3)从该三角形的边向外扩展，确定扩展边；
(4)从离散点中挑选扩展候选点，使这些点与扩展三角形第三顶点位于扩展边的异侧；
(5)在候选点中确定一个点，以扩展边的两端点和该点建立新三角形，并保证新建立的三角形在地性线的同侧，然后加入到三角形队列。
(6)判断扩展三角形的三边是否都已扩展，若是，则转到步骤(7)；否则转向步骤(3)，重复步骤(3)和(4)，直到三角形三边全部扩展完毕，转到步骤(7)。
(7)判断三角形队列中是否还有要扩展的三角形，若是，则转到步骤(3)，再次重复步骤(3)～(6)；若否，结束扩展，程序结束。当程序结束时，也就是三角形队列中所有的三角形都被扩展了，这时，所有的离散点都参与到构网过程中，并最终成为约束三角网的结点。由于每个离散点都带有高程信息，这样约束三角网中的每个三角形实际上就是空间三角形对象，地形的三维重建工作就完成了。
S102：获取强冲积河流的进出口边界水沙信息，并设置模型初始场信息。
水沙信息包括河道进出口边界的流量Q边界、水位Z边界、泥沙量S边界、河流床沙组成、泥沙干容重和湿容重等信息；初始场信息为以上各个计算单元的初始水位场Z初始、初始流速场u初始和v初始、初始含沙量场S初始以及初始河底高程z0初始，将上述信息按时间和空间进行组织得到一个数据文件。
S103：以地形网格的三角形单元为控制体系进行离散计算。
具体为以三角形单元为控制体，将强冲积河流二维水沙模型控制方程进行离散，选择求解方法进行计算，获得冲积河流各控制单元的水位计算值Z预测、流速计算值u预测和v预测、泥沙量计算值S预测和河床高程值z0预测。
根据S101和S102步信息，求解以下强冲积河流二维水沙模型：
Figure BDA0001530303620000081
Figure BDA0001530303620000082
Figure BDA0001530303620000083
Figure BDA0001530303620000084
Figure BDA0001530303620000085
式中，h为水深；u、v分别为x、y方向上的平均流速；z0为河床高程；S为垂线平均含沙量；ρm为混合物的容量；ρ0为泥沙的饱和湿容重；ρ′s为泥沙干容重；εx、εy分别为x、y方向上的泥沙紊动扩散系数；Cz为谢才系数；υt为水流紊动粘滞系数；E为近底泥沙上扬通量；D为近底泥沙沉降通量。
本实施例选用三角形网格单元来说明数值计算，便于控制方程数值离散及建立时间显式的求解数值格式。通过消除方程(1)、(2)和(3)中左端的浑水密度，得到新的浑水连续方程和浑水动量方程，并将不平衡泥沙输移方程统一写成散度形式的控制方程。
Figure BDA0001530303620000091
式中，
Figure BDA0001530303620000092
Figure BDA0001530303620000093
其中，
Figure BDA0001530303620000094
Figure BDA0001530303620000095
数值求解的控制体选用不规则三角形单元，以准确拟合实际河道的不规则岸边界和复杂的实际地形，且便于局部加密网格。计算变量(h,hu,hv,hS)置于三角形形心，控制体采用如图2的CC型网格。对方程(6)进行空间积分，并用格林公式可以离散得到：
Figure BDA0001530303620000096
令：
Figure BDA0001530303620000097
其中，nx、ny为边lij的外法线方向的矢量，nx＝Δy/Δlij；ny＝Δx/Δlij，可以得到控制方程的离散格式：
Figure BDA0001530303620000098
河床变形方程的数值离散格式为：
Figure BDA0001530303620000101
溃坝洪水计算涉及间断水流问题，间断前后，水位、流速和含沙量变化梯度大，数值格式要具有捕捉激波的性能。因此，控制体界面通量Fn的计算显得尤为重要，这是有限体积法的核心问题之一，已有许多方法可以采用，如TVD格式、MacCormack格式、BGK格式、KFVS格式、HLLC格式、Roe格式等。本实例采用基于近似黎曼解的Roe格式计算法向数值通量。格式的空间精度取决于单元界面两侧物理量的插值精度。如取界面两侧单元形心的值，并采用分片常数的重构方法，则只有一阶精度。为了获得空间二阶精度，通常采用分片线性重构的方法对计算变量进行重构，同时为了得到一个高阶且稳定的计算格式，通常对变量坡度进行限制，以调节和控制数值耗散和频散效应。本实施例采用类似MUSCL方法将空间一阶精度提高到二阶。控制体界面通量及空间二阶精度的数值计算格式参阅相关文献，如(朱宝土,史英标,穆锦斌,等.平面二维高精度水质输移数学模型及其应用[J].中国农村水利水电,2010(1):47-50)。
在溃坝洪水计算中，由于水位变化使实际计算域不断变化，为准确模拟这种动边界变化过程，通常进行单元界面的干、湿处理，有多种方法可采用，常用的有冻结法、窄缝法及最小水深等。本发明实施例采用限制水深的方法处理动边界问题，并按单元水深将网格分为干、湿及半干单元等三类。设定水深临界阈值为Epsel。
干网格：n时刻，如网格水深h<Epse1，如果相邻单元的水深也小于Epse1，则没有流量和动量通量通过公共边，如果所有相邻的水深小于Epse1，则该网格为干网格，单元的流速取为0；
半干网格：n时刻，如网格水深Epse1<h<Epse2，如果相邻单元的水深小于Epse2，则相邻边界只有流量通量，无动量通量。
湿网格：n时刻，网格水深h>Epse2。
离散式(8)一般仅适用于湿单元，对于半干半湿的非平底单元水深或水位的计算按文献(Begnudelli L,Sanders B F.Unstructured grid finite-volume algorithm forshallow-water flow and scalar transport with wetting and drying[J].Journal ofhydraulic engineering,2006,132(4):371-384)的方法进行。由于单元各边的数值通量计算在时间上为显式，对即将由湿变干的单元通量计算可能使流出该单元的水量或沙量过多导致水深(h)或沙量(hS)为负的不合理现象，需要对单元周围的单元水沙通量按流入的水沙量的比例进行校正，当周围单元的水沙通量校正后，该单元的(h,hS)赋为0，这样整个计算域的水沙量守恒且(h,hS)≥0。
S104：根据实时采集的水沙状态信息建立实时预测模型。
具体为利用在线监测设备实时采集河道水沙状态的信息，包括水位Z测 量、流速u测量和v测量、含沙量S测量，根据这些实时采集的信息，采用残差重采样粒子滤波算法，建立基于数据同化的强冲积河流水沙实时预测模型。
粒子滤波算法是对贝叶斯理论的一种蒙特卡洛算法实现。它的基本思想是用一系列加权粒子来近似状态变量的后验概率分布，随着粒子数目的增加，粒子的概率密度函数逐渐逼近状态的真实概率密度函数。
对于一个特定的水流状态过程，把网格单元状态变量(水位、流速、泥沙量、河床高程)和参数变量(糙率系数、紊动粘滞系数、扩散系数)作为滤波的基本粒子考虑，每一个滤波粒子代表一种可能的水沙状态。第k个时间步序号为i的粒子
Figure BDA0001530303620000111
Figure BDA0001530303620000112
Figure BDA0001530303620000113
Figure BDA0001530303620000114
Figure BDA0001530303620000115
Figure BDA0001530303620000116
Figure BDA0001530303620000117
Figure BDA0001530303620000118
Figure BDA0001530303620000119
Figure BDA00015303036200001110
Figure BDA00015303036200001111
式中，
Figure BDA00015303036200001112
和
Figure BDA00015303036200001113
分别为k时刻j单元处的水位、水平流速、垂向流速、含沙量、河床高程、糙率系数、紊动粘滞系数、水平扩散系数和垂向扩散系数，M为单元数目。
假设能够独立从状态的后验概率分布p(xk|z1:k)中抽取N个粒子
Figure BDA0001530303620000121
其中
Figure BDA0001530303620000122
表示粒子的状态值，
Figure BDA0001530303620000123
表示粒子的权重，则状态后验概率密度分布可通过下式近似得到：
Figure BDA0001530303620000124
式中，δ为Dirac函数。由于状态变量的后验概率分布通常未知，很难直接从p(xk|z1:k)中采样，因此采用容易采样的重要性分布函数q(xk|z1:k)进行采样。在实际应用中，多采用序贯重要性采样方法(SIS)，它采用递推的形式计算重要性权值如下式：
Figure BDA0001530303620000125
式中，
Figure BDA0001530303620000126
为重要性采样函数。一般选择先验概率分布函数为重要性采样函数，即令
Figure BDA0001530303620000127
则权值计算公式变为:
Figure BDA0001530303620000128
最终状态变量的估计值即为所有粒子状态值的加权平均：
Figure BDA0001530303620000129
在粒子滤波算法存在着粒子退化问题，即经过一定时间的迭代后，只有少数几个粒子具有较大权值，大量的计算负担用于更新对后验概率的计算贡献几乎为0的粒子。
本发明实施例采用改进的残差重采样方法来解决粒子退化问题。该方法采用哈尔顿序列和指数函数生成新的粒子，避免了原始残差重采样方法中直接通过复制生成新粒子的方法，因此能在解决粒子退化问题的同时保持粒子的多样性。
本发明实施例采用的参数优化方法，将模型参数和状态变量一起构成增广的状态变量，在估算状态变量的同时对模型参数进行同步优化。
利用残差重采样粒子滤波算法对模型参数进行同步优化时，状态变量随时间的递推靠水沙模型，但没有显式的模型对参数进行递推，因此必须构建参数递推模型。
本发明实施例采用核平滑方法，它既能实现参数递推，又能保证参数分布的方差不会随时间不断增大。假设用一系列基于蒙特卡洛采样得到的粒子
Figure BDA0001530303620000131
来近似参数向量的后验概率分布函数p(θk|z1:k-1)，
Figure BDA0001530303620000132
为粒子的状态值，
Figure BDA0001530303620000133
为粒子的权重，粒子的均值为
Figure BDA0001530303620000134
方差为Vk-1。可利用一系列高斯函数的加权和来近似参数向量的后验概率分布，即
Figure BDA0001530303620000135
式中，
Figure BDA0001530303620000136
表示均值；h2Vk-1表示方差；h是平滑参数，用来控制参数变化快慢。为了防止参数的分布过于发散，可利用下式计算
Figure BDA0001530303620000137
Figure BDA0001530303620000138
则利用核平滑方法，参数的递推过程为：
Figure BDA0001530303620000139
(1)通过上述同化水沙模型，得到同化的状态变量，包括水位Z同化、流速u同化和v同化、含沙量S同化、河底高程z0同化，以及同化的参数变量，包括糙率系数n同化、同化紊动粘滞系数υt同化、同化扩散系数εx同化和εy同化。
假设状态变量为x，代表水位、水平流速、垂向流速、含沙量、河床高程，参数变量为θ，代表糙率系数、紊动粘滞系数、水平扩散系数和垂向扩散系数，初始时刻k＝0，从初始的参数变量和状态变量的分布函数中采样，得到N个权值均等的参数变量和状态变量粒子；
(2)利用核平滑方法实现参数粒子从上一时刻到当前时刻的状态递推，即：
(3)准备水沙模型所需的文件，运行水沙模型，实现状态变量从上一时刻到当前时刻的状态递推。
(4)判断当前是否有观测数据，若有水位、泥沙含量、流速或流量数据，计算似然函数，最后对参数变量和状态变量粒子进行权重更新，并进行归一化。否则，直接进入步骤(6)，计算得到参数变量和状态变量的状态估计值。
(5)判断粒子是否发生退化，若退化，则根据权重对参数变量和状态变量粒子进行残差重采样，得到权值均等的新的N个粒子。否则，执行步骤(6)，计算得到参数变量和状态变量的状态估计值。
(6)计算参数变量和状态变量的状态估计值。
(7)令k＝k+1，返回步骤(2)继续进行迭代运算，直至所有时刻运行结束。
S105：利用上述参数和状态变量对强冲积河流进行水沙预测。
具体为根据设定的强冲积河流水沙模型计算间隔，将得到的同化状态变量和参数变量作为各单元的状态变量初始场和参数变量初始场，即初始水位场Z初始＝Z同化、初始流速场u初始＝u同化和v初始＝v同化、含沙量初始场S初始＝S同化和河底高程初始场z0初始＝z0同化，以及糙率系数初始场n初始＝n同化、紊动粘滞系数初始场υt初始＝υt同化、扩散系数初始场εx初始＝εx同化和εy初始＝εy同化。并利用前面得到的开边界的流量或流速、水位、含沙量，并利用上面的强冲积河流二维水沙模型模型，进行河道水沙预测计算，得到河道未来时刻各单元的未来水位Z预测、未来流速u预测和v预测，未来含沙量S预测和未来河底高程z0预测。
从上述技术方案可以看出，本实施例提供了一种基于数据同化的强冲积河流的水沙预测方法，具体为利用采样的地形高程散点数据，生成不规则三角形地形网格单元；获取强冲积河流的进出口边界水沙信息和初始场信息；以地形网格的三角形单元为控制体，将强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到控制体的计算值；实时采集强冲积河流的河道水沙状态信息，根据水沙状态信息建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量；根据预设的强冲积河流水沙模型的计算间隔，将同化状态变量和参数变量作为各三角形单元的状态变量初始场和参数变量初始场，并利用水沙信息和初始场信息，运行二维水沙模型对强冲积河流进行水沙预测。通过上述处理可以对强冲积河流的实时状态进行有效预测，从而可以支持水沙调控，并为防洪减灾、水环境保护和水资源管理提供依据。
并且，还能针对实际工程需要，结合已有的实时水沙观测手段，构建水沙实时数据接收系统，实时获取水位、流速或流量和含沙量等数据，将实测的水位、流速或流量和含沙量等数据考虑在模型中，使得原有的水沙数学模型的应用范围从工程设计和规划领域拓展到河道水沙实时预报领域。
另外，本发明采用先进的数据同化方法，同时考虑实测数据的测量误差和水沙模型的计算误差，对实测数据融入水沙模型的过程进行优化，使得模型预报的初始场取得最优值，从而有效提高水沙模型实时预报河道水沙的预报精度。
再有，本发明方法系统地提出了从实时水沙数据接收、水沙模型演算、数据同化初始场、水沙状态量预报等全部模块，完善了水沙模型实时预报系统框架，可实现水位流量预测、洪水预报等功能，具有针对性强、功能齐全、方便实用等特点，可应用在大江大河的河道洪水预报中，为实际防汛指挥工作提供决策支持。
需要说明的是，对于方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明实施例并不受所描述的动作顺序的限制，因为依据本发明实施例，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明实施例所必须的。
实施例二
图3为本发明实施例提供的一种基于数据同化的强冲积河流的水沙预测系统的结构框图。
参照图3所示，本实施例提供的水沙预测系统用于基于数据同化的方法对强冲积河流的水沙情况进行预测，该水沙预测系统包括网格单元生成模块10、信息获取模块20、模型计算模块30、信息采集模块40和预测输出模块50。
网格单元生成模块用于利用采样的地形高程散点数据生成不规则三角形地形网格单元。
该地形网格中会包括相应的三角形单元，该三角形单元具体可以作为该强冲积河流水沙模型离散的控制体。该模块包括生成处理单元、标识处理单元、计算处理单元和赋值处理单元，生成处理单元用于利用采样的地形高程散点数据，生成不规则三角形地形网格单元；标识处理单元用于对于不规则三角形网格，标识三角形单元、边、节点之间的拓扑关系；计算处理单元用于计算三角形单元的几何特征；赋值处理单元用于赋予三角形单元、边和节点相应的水沙状态变量和参数变量数据，并设计数据结构进行存储。
对于每一三角形、边和节点都对应一个记录，三角形的记录包括指向3条边的记录的指针。边的记录有4个指针字段，包括两个指向相邻三角形的指针和它的两个顶点指针。每个节点包括三个坐标值的字段，分别存储X、Y、Z坐标，这种拓扑网络结构的特点对于三角形的每一条边都记录了与三角形的拓扑邻接关系。
一个完整的构网原始数据，即用来构网的离散点可用如下结构来描述：
Figure BDA0001530303620000161
每个构网过程中新生成的三角形可用如下结果来描述:
Figure BDA0001530303620000162
TRIANGLE对应着三角网中的每一个三角形，数组vertex描述了构成三角形的三个离散点的序号和顺序，而数组neighbor记录了三角形三条边依次的邻接三角形序号。
建模优化算法中，首先定义一个结点队列来描述原始数据：CArray<CTinPoint*,CTinPoint*>m_aPoints。每个离散点都用CTinPoint数据结构来表示。
定义一个三角形队列来记录构建的三角形：CArray<CTriangle*,CTriangl-e*>m_aTriangles。每个三角形都用CTriangle数据结构来表示。
参见图2所示，变量L表示优化算法生成的三角形总数量，每扩展生成一个三角形，L就增加一个数值。变量k表示当前扩展三角形在整个三角形队列中的序号，即当前正在处理第几号三角形。每完成一个三角形的扩展，变量k就增加一个数值，即开始扩展下一个三角形。变量i表示扩展三角形正在处理的边。函数FindFirstTriangle用来确定第一个三角形，并以此三角形为基础进行三角网扩展。函数FindNeighborTriangle用来确定与扩展三角形边最终建立三角形的离散点，并记录新生成的三角形信息。因为在扩展新三角形时要以数据预处理阶段得到的山谷线、山脊线等地性线作为建立三角网的控制条件，即在构建三角网时，所建立的三角网不能跨越地性线，也就是要保证地性线无条件地作为三角形的边。
因此，地形的三维重建问题可以看作一个如何将离散点构建成约束Dela-unay三角网的问题。该优化算法的主要实现步骤如下。
(1)根据本优化算法对数据的要求进行数据预处理，生成离散点集合；
(2)在离散点中找出符合条件的3个点生成第一个三角形，编号为1，并加入到三角形队列中；
(3)从该三角形的边向外扩展，确定扩展边；
(4)从离散点中挑选扩展候选点，使这些点与扩展三角形第三顶点位于扩展边的异侧；
(5)在候选点中确定一个点，以扩展边的两端点和该点建立新三角形，并保证新建立的三角形在地性线的同侧，然后加入到三角形队列。
(6)判断扩展三角形的三边是否都已扩展，若是，则转到步骤(7)；否则转向步骤(3)，重复步骤(3)和(4)，直到三角形三边全部扩展完毕，转到步骤(7)。
(7)判断三角形队列中是否还有要扩展的三角形，若是，则转到步骤(3)，再次重复步骤(3)～(6)；若否，结束扩展，程序结束。当程序结束时，也就是三角形队列中所有的三角形都被扩展了，这时，所有的离散点都参与到构网过程中，并最终成为约束三角网的结点。由于每个离散点都带有高程信息，这样约束三角网中的每个三角形实际上就是空间三角形对象，地形的三维重建工作就完成了。
信息获取模块用于获取强冲积河流的进出口边界水沙信息和初始场信息。
水沙信息包括河道进出口边界的流量Q边界、水位Z边界、泥沙量S边界、河流悬沙组成、河流床沙组成、泥沙干容重和湿容重等信息；初始场信息为以上各计算单元的初始水位场Z初始、初始流速场u初始和v初始、初始泥沙场S初始以及初始河底高程z0初始，将上述信息按时间和空间进行组织得到一个数据文件。
模型计算模块用于以三角形单元为控制体进行离散计算。
具体为以三角形单元为控制体，将强冲积河流二维水沙模型控制方程进行离散，选择求解方法进行计算，获得冲积河流各控制单元的水位计算值Z预测、流速计算值u预测和v预测、泥沙量计算值S预测和河床高程值z0预测。
根据S101和S102步信息，求解以下强冲积河流二维水沙模型：
Figure BDA0001530303620000181
Figure BDA0001530303620000182
Figure BDA0001530303620000183
Figure BDA0001530303620000191
Figure BDA0001530303620000192
式中，h为水深；u、v分别为x、y方向上的平均流速；z0为河床高程；S为垂线平均含沙量；ρm为混合物的容量；ρ0为泥沙的饱和湿容重；ρs′为泥沙干容重；εx、εy分别为x、y方向上的泥沙紊动扩散系数；Cz为谢才系数；υt为水流紊动粘滞系数；E为近底泥沙上扬通量；D为近底泥沙沉降通量。
本实施例选用三角形网格单元来说明数值计算，便于控制方程数值离散及建立时间显式的求解数值格式。通过消除方程(1)、(2)和(3)中左端的浑水密度，得到新的浑水连续方程和浑水动量方程，并将不平衡泥沙输移方程统一写成散度形式的控制方程。
Figure BDA0001530303620000193
式中，
Figure BDA0001530303620000194
Figure BDA0001530303620000195
其中，
Figure BDA0001530303620000196
Figure BDA0001530303620000197
数值求解的控制体选用不规则三角形单元，以准确拟合实际河道的不规则岸边界和复杂的实际地形，且便于局部加密网格。计算变量(h,hu,hv,hS)置于三角形形心，控制体采用如图2的CC型网格。对方程(6)进行空间积分，并用格林公式可以离散得到：
Figure BDA0001530303620000198
令：
Figure BDA0001530303620000201
其中，nx、ny为边lij的外法线方向的矢量，nx＝Δy/Δlij；ny＝Δx/Δlij，可以得到控制方程的离散格式：
Figure BDA0001530303620000202
河床变形方程的数值离散格式为：
Figure BDA0001530303620000203
溃坝洪水计算涉及间断水流问题，间断前后，水位、流速和含沙量变化梯度大，数值格式要具有捕捉激波的性能。因此，控制体界面通量Fn的计算显得尤为重要，这是有限体积法的核心问题之一，已有许多方法可以采用，如TVD格式、MacCormack格式、BGK格式、KFVS格式、HLLC格式、Roe格式等。本实例采用基于近似黎曼解的Roe格式计算法向数值通量。格式的空间精度取决于单元界面两侧物理量的插值精度。如取界面两侧单元形心的值，并采用分片常数的重构方法，则只有一阶精度。为了获得空间二阶精度，通常采用分片线性重构的方法对计算变量进行重构，同时为了得到一个高阶且稳定的计算格式，通常对变量坡度进行限制，以调节和控制数值耗散和频散效应。本实施例采用类似MUSCL方法将空间一阶精度提高到二阶。控制体界面通量及空间二阶精度的数值计算格式参阅相关文献，如(朱宝土,史英标,穆锦斌,等.平面二维高精度水质输移数学模型及其应用[J].中国农村水利水电,2010(1):47-50)。
在溃坝洪水计算中，由于水位变化使实际计算域不断变化，为准确模拟这种动边界变化过程，通常进行单元界面的干、湿处理，有多种方法可采用，常用的有冻结法、窄缝法及最小水深等。本发明实施例采用限制水深的方法处理动边界问题，并按单元水深将网格分为干、湿及半干单元等三类。设定水深临界阈值为Epsel。
干网格：n时刻，如网格水深h<Epse1，如果相邻单元的水深也小于Epse1，则没有流量和动量通量通过公共边，如果所有相邻的水深小于Epse1，则该网格为干网格，单元的流速取为0；
半干网格：n时刻，如网格水深Epse1<h<Epse2，如果相邻单元的水深小于Epse2，则相邻边界只有流量通量，无动量通量。
湿网格：n时刻，网格水深h>Epse2。
离散式(8)一般仅适用于湿单元，对于半干半湿的非平底单元水深或水位的计算按文献(Begnudelli L,Sanders B F.Unstructured grid finite-volume algorithm forshallow-water flow and scalar transport with wetting and drying[J].Journal ofhydraulic engineering,2006,132(4):371-384)的方法进行。由于单元各边的数值通量计算在时间上为显式，对即将由湿变干的单元通量计算可能使流出该单元的水量或沙量过多导致水深(h)或沙量(hS)为负的不合理现象，需要对单元周围的单元水沙通量按流入的水沙量的比例进行校正，当周围单元的水沙通量校正后，该单元的(h,hS)赋为0，这样整个计算域的水沙量守恒且(h,hS)≥0。
信息采集模块用于根据实时采集的水沙状态信息建立实时预测模型。
具体为利用在线监测设备实时采集河道水沙状态的信息，包括水位Z测 量、流速u测量和v测量、含沙量S测量，根据这些实时采集的信息，采用残差重采样粒子滤波算法，建立基于数据同化的强冲积河流水沙实时预测模型。
粒子滤波算法是对贝叶斯理论的一种蒙特卡洛算法实现。它的基本思想是用一系列加权粒子来近似状态变量的后验概率分布，随着粒子数目的增加，粒子的概率密度函数逐渐逼近状态的真实概率密度函数。
对于一个特定的水流状态过程，把网格单元状态变量(水位、流速、泥沙量、河床高程)和参数变量(糙率系数、紊动粘滞系数、扩散系数)作为滤波的基本粒子考虑，每一个滤波粒子代表一种可能的水沙状态。第k个时间步序号为i的粒子
Figure BDA0001530303620000211
Figure BDA0001530303620000212
Figure BDA0001530303620000213
Figure BDA0001530303620000214
Figure BDA0001530303620000215
Figure BDA0001530303620000216
Figure BDA0001530303620000217
Figure BDA0001530303620000221
Figure BDA0001530303620000222
Figure BDA0001530303620000223
Figure BDA0001530303620000224
式中，
Figure BDA0001530303620000225
和
Figure BDA0001530303620000226
分别为k时刻j单元处的水位、水平流速、垂向流速、含沙量、河床高程、糙率系数、紊动粘滞系数、水平扩散系数和垂向扩散系数，M为单元数目。
假设能够独立从状态的后验概率分布p(xk|z1:k)中抽取N个粒子
Figure BDA0001530303620000227
其中
Figure BDA0001530303620000228
表示粒子的状态值，
Figure BDA0001530303620000229
表示粒子的权重，则状态后验概率密度分布可通过下式近似得到：
Figure BDA00015303036200002210
式中，δ为Dirac函数。由于状态变量的后验概率分布通常未知，很难直接从p(xk|z1:k)中采样，因此采用容易采样的重要性分布函数q(xk|z1:k)进行采样。在实际应用中，多采用序贯重要性采样方法(SIS)，它采用递推的形式计算重要性权值如下式：
Figure BDA00015303036200002211
式中，
Figure BDA00015303036200002212
为重要性采样函数。一般选择先验概率分布函数为重要性采样函数，即令
Figure BDA00015303036200002213
则权值计算公式变为:
Figure BDA00015303036200002214
最终状态变量的估计值即为所有粒子状态值的加权平均：
Figure BDA00015303036200002215
在粒子滤波算法存在着粒子退化问题，即经过一定时间的迭代后，只有少数几个粒子具有较大权值，大量的计算负担用于更新对后验概率的计算贡献几乎为0的粒子。
本发明实施例采用改进的残差重采样方法来解决粒子退化问题。该方法采用哈尔顿序列和指数函数生成新的粒子，避免了原始残差重采样方法中直接通过复制生成新粒子的方法，因此能在解决粒子退化问题的同时保持粒子的多样性。
本发明实施例采用的参数优化方法，将模型参数和状态变量一起构成增广的状态变量，在估算状态变量的同时对模型参数进行同步优化。
利用残差重采样粒子滤波算法对模型参数进行同步优化时，状态变量随时间的递推靠水沙模型，但没有显式的模型对参数进行递推，因此必须构建参数递推模型。
本发明实施例采用核平滑方法，它既能实现参数递推，又能保证参数分布的方差不会随时间不断增大。假设用一系列基于蒙特卡洛采样得到的粒子
Figure BDA0001530303620000231
来近似参数向量的后验概率分布函数p(θk|z1:k-1)，
Figure BDA0001530303620000232
为粒子的状态值，
Figure BDA0001530303620000233
为粒子的权重，粒子的均值为
Figure BDA0001530303620000234
方差为Vk-1。可利用一系列高斯函数的加权和来近似参数向量的后验概率分布，即
Figure BDA0001530303620000235
式中，
Figure BDA0001530303620000236
表示均值；h2Vk-1表示方差；h是平滑参数，用来控制参数变化快慢。为了防止参数的分布过于发散，可利用下式计算
Figure BDA0001530303620000237
Figure BDA0001530303620000238
则利用核平滑方法，参数的递推过程为：
Figure BDA0001530303620000239
(1)通过上述同化水沙模型，得到同化的状态变量，包括水位Z同化、流速u同化和v同化、含沙量S同化、河底高程z0同化，以及同化的参数变量，包括糙率系数n同化、同化紊动粘滞系数υt同化、同化扩散系数εx同化和εy同化。
假设状态变量为x，代表水位、水平流速、垂向流速、含沙量、河床高程，参数变量为θ，代表糙率系数、紊动粘滞系数、水平扩散系数和垂向扩散系数，初始时刻k＝0，从初始的参数变量和状态变量的分布函数中采样，得到N个权值均等的参数变量和状态变量粒子；
(2)利用核平滑方法实现参数粒子从上一时刻到当前时刻的状态递推，即：
(3)准备水沙模型所需的文件，运行水沙模型，实现状态变量从上一时刻到当前时刻的状态递推。
(4)判断当前是否有观测数据，若有水位、泥沙含量、流速或流量数据，计算似然函数，最后对参数变量和状态变量粒子进行权重更新，并进行归一化。否则，直接进入步骤(6)，计算得到参数变量和状态变量的状态估计值。
(5)判断粒子是否发生退化，若退化，则根据权重对参数变量和状态变量粒子进行残差重采样，得到权值均等的新的N个粒子。否则，执行步骤(6)，计算得到参数变量和状态变量的状态估计值。
(6)计算参数变量和状态变量的状态估计值。
(7)令k＝k+1，返回步骤(2)继续进行迭代运算，直至所有时刻运行结束。
预测输出模块用于利用上述参数对强冲积河流进行水沙预测。
具体为根据设定的强冲积河流水沙模型计算间隔，将得到的同化状态变量和参数变量作为各单元的状态变量初始场和参数变量初始场，即初始水位场Z初始＝Z同化、初始流速场u初始＝u同化和v初始＝v同化、含沙量初始场S初始＝S同化和河底高程初始场z0初始＝z0同化，以及糙率系数初始场n初始＝n同化、紊动粘滞系数初始场υt初始＝υt同化、扩散系数初始场εx初始＝εx同化和εy初始＝εy同化。并利用前面得到的开边界的流量或流速、水位、含沙量，并运行上面的强冲积河流二维水沙模型，进行河道水沙预测计算，得到河道未来时刻各单元的未来水位Z预测、未来流速u预测和v预测，未来含沙量S预测和未来河底高程z0预测。
从上述技术方案可以看出，本实施例提供了一种基于数据同化的强冲积河流的水沙预测系统，具体为利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元；获取强冲积河流的进出口边界水沙信息和初始场信息；以三角形地形网格单元为控制体，将强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到控制体的计算值；实时采集强冲积河流的河道水沙状态信息，根据水沙状态信息建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量；根据预设的强冲积河流水沙模型的计算间隔，将同化状态变量和参数变量作为各三角形单元的状态变量初始场和参数变量初始场，并利用水沙信息和初始场信息，运行二维水沙模型对强冲积河流进行水沙预测。通过上述处理可以对强冲积河流的实时状态进行有效预测，从而可以支持水沙调控，并为防洪减灾、水环境保护和水资源管理提供依据。
对于装置实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。
本领域内的技术人员应明白，本发明实施例的实施例可提供为方法、装置、或计算机程序产品。因此，本发明实施例可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明实施例可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明实施例是参照根据本发明实施例的方法、终端设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理终端设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理终端设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理终端设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理终端设备上，使得在计算机或其他可编程终端设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程终端设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
尽管已描述了本发明实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
以上对本发明所提供的技术方案进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于数据同化的强冲积河流的水沙预测方法，其特征在于，包括步骤：
利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元；
对组成所述三角形地形网格单元的三角形单元，所述三角形单元的边和所述三角形单元的节点，赋予各自对应的水沙状态变量和参数变量数据；所述状态变量包括：水位、流速、泥沙量、河床高程；所述参数变量包括：糙率系数、紊动粘滞系数、扩散系数；
获取强冲积河流进出口边界的水沙信息，并设置模型的初始场信息；以所述三角形地形网格单元为控制体，将所述强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到所述控制体水沙状态信息的计算值，作为模型的预测场信息；
实时采集所述强冲积河流的河道水沙状态信息，根据所述水沙状态信建立基于数据同化的实时预测模型，并得到同化状态变量和参数变量，作为模型的分析场信息；
根据预设的强冲积河流水沙模型的计算间隔，将步骤所述同化状态变量和所述参数变量的分析场作为各所述三角形单元的状态变量初始场和参数变量初始场，利用所述水沙信息和所述初始场信息以及所述二维水沙模型对所述强冲积河流进行水沙预测；
根据所述水沙状态信息建立基于数据同化的实时预测模型，包括：
采用残差重采样粒子滤波算法，建立基于数据同化的强冲积河流水沙实时预测模型；其中，所述残差重采样粒子滤波算法中的基本粒子是所述状态变量以及所述参数变量；每一个基本粒子表示一种水沙状态；
所述方法还包括：
按单元水深将网格分为干、湿及半干单元，设定水深临界阈值为Epsel；其中，在网格水深h<Epse1并且其相邻单元的水深小于Epse1，该网格为干网格，单元的流速取为0；在网格水深Epse1<h<Epse2并且其相邻单元的水深小于Epse2，该网格为半干网格，相邻边界只有流量通量，无动量通量；在网格水深h>Epse2，该网格为湿网格。
2.如权利要求1所述的水沙预测方法，其特征在于，所述利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元，包括：
利用采样的地形高程散点数据生成不规则三角形地形网格单元；
标识所述地形网格的三角形单元及其边、节点之间的拓扑关系；
计算所述三角形单元的几何特征。
3.如权利要求1所述的水沙预测方法，其特征在于，所述水沙信息包括河道进出口边界的流量边界、水位边界、含沙量边界、河流悬沙级配组成、河流床沙级配组成、泥沙干容重和泥沙湿容重中的部分或全部。
4.如权利要求1所述的水沙预测方法，其特征在于，所述计算值包括水位计算值、流速计算值、含沙量计算值和河床高程计算值中的部分或全部。
5.如权利要求1所述的水沙预测方法，其特征在于，所述实时采集的水沙状态信息包括水位测量值、流量或流速测量值和含沙量测量值的部分或全部。
6.一种基于数据同化的强冲积河流的水沙预测系统，其特征在于，包括：
网格单元生成模块，用于利用采样得到的地形高程散点数据，生成不规则三角形地形网格单元；
所述网格单元生成模块包括赋值处理单元，所述赋值处理单元用于对组成所述三角形地形网格单元的三角形单元，所述三角形单元的边和所述三角形单元的节点，赋予各自对应的水沙状态变量和参数变量数据；
信息获取模块，用于获取强冲积河流的进出口边界水沙信息和初始场信息；
模型计算模块，用于以所述地形网格的三角形单元为控制体，将所述强冲积河流的二维水沙模型进行离散，选择求解方法进行计算，得到所述控制体的计算值；
信息采集模块，用于实时采集所述强冲积河流的河道水沙状态信息，根据所述水沙状态信息建立基于数据同化的实时预测模型，并得到同化的状态变量和参数变量；
预测输出模块，用于根据预设的强冲积河流水沙模型的计算间隔，将步骤所述同化状态变量和所述参数变量作为各所述三角形单元的状态变量初始场和参数变量初始场，利用所述水沙信息和所述初始场信息以及所述二维水沙模型对所述强冲积河流进行水沙预测；
所述信息采集模块具体用于采用残差重采样粒子滤波算法，建立基于数据同化的强冲积河流水沙实时预测模型；其中，所述残差重采样粒子滤波算法中的基本粒子是所述状态变量以及所述参数变量；每一个基本粒子表示一种水沙状态；
还包括用于按单元水深将网格分为干、湿及半干单元，设定水深临界阈值为Epsel的模块；其中，在网格水深h<Epse1并且其相邻单元的水深小于Epse1，该网格为干网格，单元的流速取为0；在网格水深Epse1<h<Epse2并且其相邻单元的水深小于Epse2，该网格为半干网格，相邻边界只有流量通量，无动量通量；在网格水深h>Epse2，该网格为湿网格。
7.如权利要求6所述的水沙预测系统，其特征在于，所述网格单元生成模块包括：
生成处理单元，用于利用采样的地形高程散点数据生成不规则三角形地形网格单元；
标识处理单元，用于标识所述地形网格的三角形单元及其边、节点之间的拓扑关系；
计算处理单元，用于计算所述三角形单元的几何特征。
8.如权利要求6所述的水沙预测系统，其特征在于，所述水沙信息包括河道进出口边界的流量边界、水位边界、含沙量边界、河流悬沙级配组成、河流床沙级配组成、泥沙干容重和泥沙湿容重中的部分或全部。
9.如权利要求6所述的水沙预测系统，其特征在于，所述计算值包括水位计算值、流速计算值、含沙量计算值和河床高程计算值中的部分或全部。
10.如权利要求6所述的水沙预测系统，其特征在于，所述实时采集的水沙状态信息包括水位测量值、流量或流速测量值和含沙量测量值的部分或全部。
CN201711461256.2A 2017-12-28 2017-12-28 一种基于数据同化的强冲积河流的水沙预测方法和系统 Active CN108334660B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711461256.2A CN108334660B (zh) 2017-12-28 2017-12-28 一种基于数据同化的强冲积河流的水沙预测方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711461256.2A CN108334660B (zh) 2017-12-28 2017-12-28 一种基于数据同化的强冲积河流的水沙预测方法和系统

Publications (2)

Publication Number Publication Date
CN108334660A CN108334660A (zh) 2018-07-27
CN108334660B true CN108334660B (zh) 2021-07-20

Family

ID=62924726

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711461256.2A Active CN108334660B (zh) 2017-12-28 2017-12-28 一种基于数据同化的强冲积河流的水沙预测方法和系统

Country Status (1)

Country Link
CN (1) CN108334660B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110059870B (zh) * 2019-04-04 2020-08-25 长江航道规划设计研究院 基于bim和gis的航道整治建筑物维护分析方法
CN110046469B (zh) * 2019-05-13 2020-02-11 水利部交通运输部国家能源局南京水利科学研究院 多约束条件下水电站坝前河床冲淤变形的计算方法
CN110119590B (zh) * 2019-05-22 2020-08-11 中国水利水电科学研究院 一种基于多源观测数据的水质模型粒子滤波同化方法
CN110795792B (zh) * 2019-11-13 2020-06-23 水利部交通运输部国家能源局南京水利科学研究院 一种由工程建设引发的河道强紊动区的河床变形预测方法
CN112381402B (zh) * 2020-11-13 2023-06-27 中国科学院自动化研究所 平行智能风沙防护治理决策支持方法及系统
CN112945859B (zh) * 2021-01-19 2023-03-24 生态环境部南京环境科学研究所 一种用于游荡性河流整治的生态护岸系统及控制方法
CN112857505B (zh) * 2021-02-23 2022-05-27 长江水利委员会水文局 一种急速涨落水位全过程应急递测方法
CN113327323B (zh) * 2021-06-09 2022-11-11 四川大学 基于散点数据的水体环境地形构建方法
CN114091359B (zh) * 2022-01-21 2022-05-17 中国长江三峡集团有限公司 一种水库水沙预测模型训练、水库水沙预测方法及装置
CN114547951B (zh) * 2022-04-24 2022-07-22 浙江远算科技有限公司 一种基于数据同化的大坝状态预测方法及系统
CN116415508B (zh) * 2023-06-12 2023-10-13 珠江水利委员会珠江水利科学研究院 一种河口二维泥沙模型生成方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004106974A2 (en) * 2003-05-23 2004-12-09 Exxonmobil Upstream Research Company Method for predicting grain size distribution from the shape of a sedimentary body
CN101482516A (zh) * 2009-02-12 2009-07-15 彭宣戈 一种图像在线测量河流泥沙含沙量的方法
CN103886187A (zh) * 2014-03-06 2014-06-25 清华大学 一种基于数据同化的河道水沙实时预测方法
CN107451383A (zh) * 2017-09-29 2017-12-08 中国水利水电科学研究院 一种平面二维水沙数学模型初始床沙级配的率定方法
CN107480384A (zh) * 2017-08-24 2017-12-15 北方民族大学 河道水流泥沙平面二维数值模拟方法和系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2849211B1 (fr) * 2002-12-20 2005-03-11 Inst Francais Du Petrole Methode de modelisation pour constituer un modele simulant le remplissage multilithologique d'un bassin sedimentaire

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004106974A2 (en) * 2003-05-23 2004-12-09 Exxonmobil Upstream Research Company Method for predicting grain size distribution from the shape of a sedimentary body
CN101482516A (zh) * 2009-02-12 2009-07-15 彭宣戈 一种图像在线测量河流泥沙含沙量的方法
CN103886187A (zh) * 2014-03-06 2014-06-25 清华大学 一种基于数据同化的河道水沙实时预测方法
CN107480384A (zh) * 2017-08-24 2017-12-15 北方民族大学 河道水流泥沙平面二维数值模拟方法和系统
CN107451383A (zh) * 2017-09-29 2017-12-08 中国水利水电科学研究院 一种平面二维水沙数学模型初始床沙级配的率定方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Harmonize input selection for sediment transport prediction;Haitham Abdulmohsin Afan 等;《Journal of Hydrology》;ELSEVIER;20170608;第366-375页 *
Water and sediment transport modeling of a large temporary river basin in Greece;C. Gamvroudis 等;《Science of the Total Environment》;ELSEVIER;20141210;第354-365页 *
水沙实时预测数学模型研究;赖瑞勋 等;《水利学报》;20140831;第45卷(第8期);第930-937页 *

Also Published As

Publication number Publication date
CN108334660A (zh) 2018-07-27

Similar Documents

Publication Publication Date Title
CN108334660B (zh) 一种基于数据同化的强冲积河流的水沙预测方法和系统
CN108256177B (zh) 一种河流水沙模型的参数优化方法和系统
CN109657418B (zh) 一种基于mike21的湖泊水环境容量计算方法
Rodriguez et al. High-resolution numerical simulation of flow through a highly sinuous river reach
Hung et al. An artificial neural network model for rainfall forecasting in Bangkok, Thailand
Kim et al. Mesh type tradeoffs in 2D hydrodynamic modeling of flooding with a Godunov-based flow solver
JP6736731B2 (ja) エネルギー過程の類似性に基づく水文モデル流域スケール決定方法
CN110232471B (zh) 一种降水传感网节点布局优化方法及装置
Roushangar et al. Evaluation of genetic programming-based models for simulating friction factor in alluvial channels
CN109726433B (zh) 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
CN102968529A (zh) 一种供水管网模型计算结果不确定性区间的量化方法
CN112464584A (zh) 自由表面流的水位和流量推求方法
Allen et al. Characterizing the impact of geometric simplification on large woody debris using CFD
CN108205713B (zh) 一种区域风电功率预测误差分布确定方法和装置
Üneş et al. Determination of groundwater level fluctuations by artificial neural networks
CN115688246A (zh) 一种局部坐标系下的水库库容模拟方法及装置
CN117195603B (zh) 基于高分辨率遥感要素的洪涝灾害推演方法、设备及介质
CN108920737A (zh) 一种水动力模型的粒子滤波同化方法、装置和计算设备
Zhang et al. Integrating 1D and 2D hydrodynamic, sediment transport model for dam-break flow using finite volume method
CN115952753B (zh) 一种ca与lbm结合的碎屑流模拟方法、系统及设备
Meng et al. Variable infiltration capacity model with BGSA-based wavelet neural network
Vosoughi et al. The application of Bayesian model averaging based on artificial intelligent models in estimating multiphase shock flood waves
KR101487846B1 (ko) 임계마름수심 기법을 이용하여 마름/젖음 조건의 부여를 통해 2차원 천수흐름을 해석하는 방법
Thakur et al. Exploring CCHE2D and its sediment modelling capabilities
CN115563908A (zh) 一种水库洪峰预测方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant