CN111191372B - 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 - Google Patents
一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 Download PDFInfo
- Publication number
- CN111191372B CN111191372B CN202010000075.5A CN202010000075A CN111191372B CN 111191372 B CN111191372 B CN 111191372B CN 202010000075 A CN202010000075 A CN 202010000075A CN 111191372 B CN111191372 B CN 111191372B
- Authority
- CN
- China
- Prior art keywords
- flood
- field
- grid
- water
- flow
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 70
- 238000004088 simulation Methods 0.000 title claims abstract description 56
- 230000003993 interaction Effects 0.000 title claims abstract description 33
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 100
- 238000005516 engineering process Methods 0.000 claims abstract description 16
- 230000000694 effects Effects 0.000 claims abstract description 14
- 230000010354 integration Effects 0.000 claims abstract description 13
- 239000003086 colorant Substances 0.000 claims abstract description 7
- 230000000007 visual effect Effects 0.000 claims description 14
- 238000011156 evaluation Methods 0.000 claims description 12
- 238000013507 mapping Methods 0.000 claims description 11
- 230000008859 change Effects 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 7
- 238000009792 diffusion process Methods 0.000 claims description 7
- 238000009499 grossing Methods 0.000 claims description 6
- 239000002245 particle Substances 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 5
- 238000009877 rendering Methods 0.000 claims description 5
- 230000001133 acceleration Effects 0.000 claims description 4
- 239000006185 dispersion Substances 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000010606 normalization Methods 0.000 claims description 4
- 238000012545 processing Methods 0.000 claims description 3
- 238000012886 linear function Methods 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 10
- 238000004891 communication Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000002452 interceptive effect Effects 0.000 description 5
- 239000012530 fluid Substances 0.000 description 4
- 238000006424 Flood reaction Methods 0.000 description 3
- 238000007726 management method Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000011218 segmentation Effects 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 239000003344 environmental pollutant Substances 0.000 description 2
- 231100000719 pollutant Toxicity 0.000 description 2
- 230000002265 prevention Effects 0.000 description 2
- 230000000452 restraining effect Effects 0.000 description 2
- 230000007480 spreading Effects 0.000 description 2
- 239000002352 surface water Substances 0.000 description 2
- 238000012800 visualization Methods 0.000 description 2
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 229920000642 polymer Polymers 0.000 description 1
- 239000013049 sediment Substances 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/001—Texturing; Colouring; Generation of texture or colour
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于二维水动力学方法的大规模洪水场景模拟预警交互方法和系统。1)对大规模洪水这一灾害现象进行了建模与绘制。2)将半拉格朗日隐式积分法求解浅水方程方法应用于洪水模拟，达到了大时间步长并进行了体积、动量守恒修正。3)应用GPGPU编程基于Computer Shader框架，实现洪水模拟的大规模实时并行。4)引入自适应网格曲面细分技术和纹理平铺技术，实现了真实感水面，且有效传达了洪水的高度、流速相关信息。5)应用三次B样条插值方式解决水体边界锯齿问题。6)提出了综合洪水淹没水深和建筑物类别的灾害评估模型，根据评估结果实时变换建筑物颜色实现分级预警。7)支持洪水演进过程中的实时交互，展示设置挡板阻挡洪水的效果。
Description
技术领域
本发明涉及洪水场景模拟方法，尤其涉及一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统。
背景技术
洪水是指河流、湖泊、海洋等水文环境上涨超过常规水位的一种自然现象。我国洪水灾害在空间上既有普遍性，又具有区域性；在时间上既表现出无序的非稳定性，又存在有序的韵律性、周期性。大量研究表明：洪水灾害具有不均匀性与差异性、多样性、随机性与可预测性、突变性与规律性等复杂性的特点。对于洪水的研究和仿真有着十分重要的意义，一个成熟完善的模拟系统给普通人洪水预警，提高对于洪水的防范意识，也能够给决策者做出决策提供可靠的参考，可以大大减少洪水带来的人员伤亡和财产损失。
洪水模拟的特殊性在于其不仅仅是流体模拟应用于视觉效果或演示，它对于水文专家或是可能会受到洪水影响的普通民众都具有很强的现实意义，一个可以交互的洪水自然灾害模拟系统可以提高普通民众的防洪意识和预警作用，也可以辅助水文专家或是防洪决策者预见到洪水的灾害做出正确的决策。
欧洲方面，奥地利的维也纳VRVis研究中心与维也纳技术大学合作已经发展出一套比较成熟的理论体系，并且在此基础上已经有成熟的、投入现实使用的洪水的决策系统Visdom。目前该系统很好的模拟了洪水和暴雨事件，集成了真实地形建模、洪水模拟和决策管理模块，对于业界洪水的模拟仿真系统的建设有着重要的参考意义。国内方面，西安理工大学的水模拟及灾害管理研究团队主要从事地表水及其附随过程数值模型的理论推导和实际应用、城市及流域洪涝管理、城市水利和水利遥感技术的研究。对于洪水研究，提出了GAST(GPU Accelerated Surface Water Flow and Transport)模型，即利用GPU加速的模拟水体流动和传输的模型，该模型整合了水体的模拟、泥沙和污染物的转移、地下输水管网以及海绵城市技术，该模型对比欧洲团队，不仅注重洪水本身的模拟计算，对于洪水泛滥携带的污染物和土质颗粒也进行了模拟计算，系统更加注重洪水整体的长远决策，例如海绵城市相关的技术探讨。不足点在于缺少比较完善的可视化和交互模块，无法给决策者的即时决策提供一些指导。
发明内容
本发明的目的在于开发一个大规模洪水场景模拟预警交互系统，拥有完善的可视化模块和交互模块，为此提出了一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统，其利用半拉格朗日隐式积分法求解的浅水方程为框架，采用GPU加速的方法达到实时性，并且综合了多学科的优点，生成真实有效的洪水模拟预警场景。
一种基于二维水动力学的大规模洪水场景模拟预警交互方法，包括以下步骤：
1)将模拟场景在水平面上细分为均匀的网格，将水体高度和基础地形高度存储在网格中心，将速度存储在网格边界；采用基于浅水方程的二维欧拉网格法模拟洪水演进，将浅水方程转化为拉格朗日表达形式，得到初步建立的物理场；所述物理场包括速度场和高度场；
2)在交错网格上对物理场进行对流更新，并在更新对流项时采用大时间步长并引入体积守恒与动量守恒，得到基于体积守恒和动量守恒变形后的浅水方程；根据变形后的浅水方程，采用隐式积分法和雅克比迭代法求解更新后的高度场，根据更新后的高度场梯度进一步求解得到更新后的速度场；
3)根据步骤2)更新后的物理场，采用动态自适应网格细分技术抬升渲染水面，采用纹理平铺技术实现速度场的水流绘制，并采用双三次B样条曲线插值法对绘制的洪水边缘做光滑处理；所述水流绘制包括单块平板单一流速的绘制、固定流场的绘制以及时变流场的绘制；
4)基于Hydrotec模型，根据淹没水深和建筑物类别设计洪灾损失评估函数，其中每一类建筑物对应一个洪灾损失评估函数；对洪灾损失评估函数值做归一化处理，并将归一化结果设为不同等级，每一个等级使用不同的颜色来表示建筑的受损程度，用于洪水场景绘制时实时改变建筑的外观颜色，实现可视化预警；
5)当收到预警信号时，选择在建筑物周围模拟修建挡板操作，调整挡板所设位置对应网格水底地形的高度，然后重复步骤2)-步骤4)完成洪水场景的实时更新，实现洪水模拟场景下的实时交互。
本发明应用隐式积分法结合半拉格朗日法达到大时间步长。积分方法有显式、隐式之分，显式积分法时间步长需满足Courant-Friedrichs-Lewy(CFL)条件：Δt·vi+1/2≤min(Δxi,Δxi+1),vi+1/2表示i与i+1相邻网格间的流体速度，ΔxiΔxi+1分别为两个网格的间隔。隐式积分法除了结合上一时间步长的状态，还结合了当前时间步长的状态，因此允许的时间步长相对于显示积分法较大，结合半拉格朗日法可将时间步长突破CFL限制，达到大时间步长。
本发明还公开了一种用于实现上述基于二维水动力学的大规模洪水场景模拟预警交互方法的预警交互系统，包括：预处理模块、求解模块、绘制模块、可视化预警模块和交互模块。所述的可视化预警模块可以根据设计的洪灾损失评估函数，对评估函数结果归一化处理和分级，并根据分级结果使用不同的颜色表示建筑的受损程度；所述的交互模块可以实现用户和系统的实时交互，当用户收到预警信号时，能够选择在建筑物周围模拟修建挡板操作，并实时显示修建挡板后的可视化效果。
本发明的优点在于：
洪水模拟问题涵盖了计算流体力学、数学、计算机等多门学科，往往包含大量的物理性质和数学运算，需要进行合理准确的抽象和建模，并通过计算机高效准确地模拟仿真，基于二维水动力学、计算机图形学等多学科的知识，在计算机图形学领域对大规模洪涝这一灾害现象进行了建模与实时绘制，这在我国的灾难防治领域有较大的应用前景。
(1)本发明在传统浅水方程显式积分求解方法的基础上应用半拉格朗日隐式积分法，引入体积守恒和动量守恒，达到了大时间步长，并采用GPU对这一过程进行加速，从而能实时模拟洪水演进；
(2)本发明引入了自适应网格曲面细分技术和纹理平铺技术，解决了水体边界锯齿问题，实现了真实感水面的同时，有效传达了洪水的高度、流速相关信息；
(3)本发明提出了一种可视化预警模块及其实现可视化的方法，构造了基于淹没水深和建筑物类别的洪水灾害评估预警模型，使得洪水灾害的损失得以量化，利用热力图进行预警，提高了预警有效性；
(4)本发明考虑到了洪水模拟的特殊性，其不仅仅是流体模拟应用于视觉效果或演示，它对于水文专家或是可能会受到洪水影响的普通民众都具有很强的现实意义，一个可以交互的洪水自然灾害模拟系统可以提高普通民众的防洪意识和预警作用，也可以辅助水文专家或是防洪决策者预见到洪水的灾害做出正确的决策；本发明提出的交互模块，可以在洪水泛滥时添加一些挡板，从而改变洪水泛滥的状态，这对于遏制洪水漫延、保护重要建筑物有着重要的意义。对于本系统和发明而言，修建挡板带来的结果在于相关地形的高度抬升，由于洪水模拟中的地形抬升是在GPU中实现的，而交互操作则与CPU相关，为减少CPU与GPU的通信代价，因此选择直线构建的实现方式尽可能的减少通信数据，从而达到实时交互的目的。
最后，本发明将上述方法整合，构造了一个基于二维水动力学的大规模洪水场景模拟预警交互系统，成功解决了计算机图形学领域对此涉及较少的问题，说明书最后对实际城市场景进行了洪水模拟，展示出了本方法的场景真实感。
附图说明
图1本发明选取的实际上海浦东江水场景；
图2本发明实际洪水场景模拟不同阶段演示；
图3本发明的仿真系统整体流程展示；
图4本发明曲面细分技术展示；
图5本发明真实感水面展示；
图6本发明时变流场水面展示；
图7本发明光滑边缘提取；
图8本发明预警等级颜色划分；
图9本发明洪水漫延阶段预警；
图10添加挡板前后效果图对比。
具体实施方式
以下通过具体实施例对本发明进行详细描述。
如图3所示为本发明的流程示意图，具体包括以下步骤：
步骤一、预处理、初始化：
1.1)将模拟场景在水平面上细分为m*n均匀网格，x方向上m等分，y方向上n等分，基础网格取边长为Δx的正方形，时间步长取Δt。为保证稳定性，使用交错网格表示物理量，将水体高度和基础地形高度存储在网格中心，将速度存储在网格边界；图1为本发明选取的上海浦东江水场景。
1.2)将浅水方程写成拉格朗日形式：
其中，u代表x方向的速度，v表示y方向的速度，h表示水面总高度，b表示水底地形的高度，η＝h-b为水的深度，g表示重力加速度，等式左边为各物理量与速度场相关的对流项，等式右边为扩散项。
步骤二、洪水模拟：
洪水模拟过程就是对步骤1.2中的浅水方程进行求解的过程，求解步骤为：在交错网格上1.将物理量场(高度场、速度场)根据当前速度场分别进行对流更新。2.对高度场进行扩散项更新。3.利用更新后的高度场计算速度场的扩散项。
2.1)基于体积守恒和动量守恒，在交错网格上对物理场进行对流更新：
在每个网格位置执行虚拟粒子的向后回溯来计算网格上的对流，对于需要更新的物理场，以物理位置p处的网格单元(i,j)为例，该粒子在上一个时间步长的位置为pbackward＝p-Δt·ν(p)，其中p为当前网格中心点的位置，Δt为时间步长，ν(p)为当前网格中心点通过差值获得的速度；由于pbackward计算结果可能为网格任何位置，因此使用双线性插值估计中心值。其可保证估算数值不超过周围网格的最大、最小值，保证稳定性。
对物理场更新对流项时，先更新高度场，再更新速度场；将当前物理场F中的网格的物理量，按照一定比例加权到下一时间步物理场F′的网格中，更新方式为：F′＝Interpolate(F,pbackward)，其中Interpolate(·)表示插值，确保分配时F′中每个网格的物理量比例都为1，保证整体的体积守恒，具体步骤为：
S1：构建一个数组W，用于存储各个网格实际被分配的比例；
S2：执行对流计算，得到对应网格的回溯点，即该粒子在上一个时间步长的位置pbackward，根据双线性插值原则得到需要从周围四个网格提取的比例；
S3：针对比例权重大于1的部分做归一化处理，针对比例小于1的网格进行分散操作；所述的分散操作具体为：
将多余的物理量分散到对流后F′的网格中，计算映射关系qforward＝q+Δt·ν(q),其中y表示比例小于1的网格，ν(q)为网格q中心点通过差值获得的速度，qforward表示q映射后的网格，对qforward周围四个网格做双线性插值估算比例，按照估算后的比例将网格q中的多余物理量分散到四个网格中；上述操作可同时用于高度场和速度场，来保证体积、动量的守恒。值得注意的是，必须首先进行高度场的对流更新，然后才能完成速度场对流更新的补正。
2.2)更新浅水方程形式：
更新对流项引入体积守恒和动量守恒后，浅水方程变形为：
其中，上标n+1表示当前时间步长的物理量，n表示前一时间步长的物理量，
分别表示更新过对流项后的物理量；
2.3)更新高度场：
计算物理场中的扩散项，使用隐式积分法，根据更新后的浅水方程得到：
其中，
然后使用雅克比迭代法并行求解大型稀疏线性方程组hn+1 i,j＝f(hn+1 i-1,j,hn +1 i+1,j,hn+1 i,j-1,hn+1 i,j+1)，其中f为线性函数，下标(i,j)表示网格位置，从而得到更新后的高度场；
2.4)更新速度场：
根据步骤2.3)得到的高度场梯度，得到：
进一步求解得到更新后的速度场。
采用Computer Shader框架进行网格计算，使用两张分辨率与网格粒度相同分辨率的RGBAFloat形式的纹理贴图M、M’来代表高度场网格，R通道存储地形高度b,G通道存储水深η,B通道存储速度u,A通道存储速度v，使用上述半拉格朗日隐式积分法求解浅水方程，由M更新对流项后存储到M’，而后针对M’计算扩散项后得到对应于下一时间步长的M，不断迭代求解得到洪水演进过程。
步骤三、洪水绘制：
3.1)根据更新后的物理场，采用动态自适应网格细分技术抬升渲染水面：
对于水面高度变化剧烈的地方会相应的在原有的三角网格顶点的基础上增加一些新的顶点从而提高采样水体高度的精确性，而在水面变化平缓的地方则分布少一些顶点，具体而言：计算水面高度变化程度
其中hi和hj表示相邻的两个网格对应的水面总高度，定义阈值1.5和3，将水面高度变化划分为三个等级：当
时，无需细分网格，当
时，中度细分网格，当
时，高度细分网格；图4为本发明使用自适应曲面细分技术绘制水面高度的效果展示，其中(a)为未细分效果，(b)为中度细分效果，(c)为高度细分效果。
3.2)采用纹理平铺技术实现速度场的水流绘制：
a.单块平板单一流速的绘制：
水面流动由流速大小和方向决定，流速大小决定采样法线贴图的偏移幅度，流速方向决定偏移的角度；以(1,0)方向为默认方向，默认的偏移后的UV纹理坐标为(速度*时间步长,0)，根据水流方向生成对应的旋转矩阵：
其中，flowdir为水流方向矢量，flowdir.x表示flowdir在X方向上的分量大小，flowdir.y表示flowdir在Y方向上的分量大小；
根据旋转矩阵对纹理坐标进行旋转变化，再根据水流的速度计算出UV偏移：
uv＝mul(rotateMatrix,uv)-(flowSpeed·Δt)
其中，flowSpeed表示流速大小，uv代表采样坐标，mul(·)表示乘法运算；
所述的采样法线贴图具体为：为每个三角网格光栅化为的多个片段提供一个法线，利用一张2D纹理存储法线数据，进行采样和映射后可得到对应的法线向量，使得水面表现为由很多微小平面组成，得到波动效果；图5为本发明使用纹理平铺技术绘制得到的波动效果真实感水面展示。
b.固定流场的绘制：
根据所述的单块平板单一流速的绘制方法，利用多个平板平铺来表达整个流场，每块平板按照对应流场的流速分别绘制；为解决简单拼凑带来的间隙问题，引入更多的纹理块，并且将引入的纹理块在原有块的基础上以固定的偏移平铺展开，结果块的法线贴图由展开块的双线性插值得到；
c.时变流场的绘制：
为解决前后时间步长间速度方向大小轻微改变后引起的水面抖动和混乱问题，以10°一个刻度划分将本来连续变化的流速方向离散成36个固定方向的流场，将采样到的流速进行流速方向上的分类过滤后，再利用分类过滤后的方向进行法线贴图的偏移；将本来连续变化的速度大小，按照0.1m/s的刻度分割，将采样到的流速进行流速大小的分类过滤后，再利用分类过滤后的方向进行法线贴图的偏移；图6为本发明对时变流场江面和溃堤局部绘制展示，有效解决了前后时间步长间流速方向大小变化带来的剧烈抖动现象。
3.3)采用双三次B样条曲线插值法对绘制的洪水边缘做光滑处理：
图7为本发明使用双三次B样条曲线插值法平滑洪水边缘效果(c)和使用最近邻域插值法(a)、双线性插值法(b)获得的洪水边缘细节对比。(a)图水流边缘有明显的锯齿，(b)图中洪水边缘锯齿相对(a)得到一定改善，(c)图洪水边缘相较于(a)(b)明显更为光滑自然。采用双三次B样条曲线作为三次样条插值的基函数，
对于待估算点p，其坐标为(x,y),提取出周围最靠近的16个网格节点，X方向坐标为x1,x2,x3,x4；Y方向坐标为y1,y2,y3,y4,先在X方向进行插值求出距待估算点p最近的a(x,y1),b(x,y2),c(x,y3),d(x,y4)四个网格点的物理场量值：
其中，m＝1,2,3,4分别对应a,b,c,d，x表示X方向上坐标数值，xt表示下标t对应的X方向上的坐标，ym表示下标m对应的Y方向上的坐标；f表示物理场的数值，f(x,ym)为坐标为(x,ym)的点的物理场数值，k为基函数；
然后在Y方向上对a,b,c,d四个点再次插值求出p点的数值：
其中，yt表示下标t对应的Y方向上的坐标，f(x,y)为p点物理场的数值；
通过步骤3.3)，扩展的16个邻域网格允许将插值样条线的起点和终点的导数与相邻的插值面片匹配，从而生成导数平滑连续的表面。图2为本发明在图1场景模拟洪水溃堤演进的不同阶段的展示，其中(a)为平静江面图；(b)为堤坝决口，江水开始从决口出漫延；(c)为洪水持续漫延，开始淹没建筑物；(d)较多区域被洪水淹没。
步骤四、建筑物预警：
4.1)设计基于淹没水深和建筑物类别的洪灾损失评估函数：
基于Hydrotec模型，水深与灾害程度呈现平方根级别的正相关，其相关系数与建筑物的4个类别：办公楼、商业用地、住宅用地、基础设施相关，共得到四种建筑物的洪灾损失评估函数；
4.2)设计可视化预警：
根据步骤4.1)得到的四种建筑物的洪灾损失评估函数，将洪灾损失评估函数值归一化到[0,1]之间，每隔0.2取一个等级，每一个等级使用不同的颜色来表示建筑的受损程度，图8为本发明预警等级颜色划分，使得洪水场景实时绘制时改变建筑的外观颜色。
步骤五、交互操作：
对于已经发生的洪水灾害，修建挡板、沙袋这一决策算是最常见也是最容易实施的，可以在洪水泛滥时添加一些挡板，从而改变洪水泛滥的状态。这对于遏制洪水漫延、保护重要建筑物有着重要的意义。
对于本系统而言，修建挡板带来的结果在于相关地形的高度抬升，由于洪水模拟中的地形抬升是在GPU中实现的，而交互操作则与CPU相关，为减少CPU与GPU的通信代价，因此选择直线构建的实现方式尽可能的减少通信数据，从而达到实时交互的目的。当收到预警信号时，选择在建筑物周围模拟修建挡板操作，调整挡板所设位置对应网格水底地形的高度，由于洪水模拟中的地形抬升是在GPU中实现的，而交互操作则与CPU相关，为减少CPU与GPU的通信代价，因此选择直线构建的实现方式尽可能的减少通信数据。地形抬升后重复步骤2)-步骤4)完成洪水场景的实时更新，实现洪水模拟场景下的实时交互。
图9为本发明模拟过程中不同阶段预警：(a)为第一栋建筑物受灾，其颜色发生变化；(b)为建筑物受灾程度加剧，显示颜色加深；(c)为洪水持续漫延，更多建筑物处于警急状态，更多建筑物变色；(d)为集成了谷歌地球卫星纹理的效果图。在模拟过程中实时使用热力图预警，预处理阶段，对场景中的建筑物编号、分类、提取出建筑物的边缘图。对于待模拟的城市模型，首先对其中的各个建筑物进行编号(buildingIndex)，为每个房屋绑定有关洪水预警的材质并且设定与建筑物有关的参数，将提取出的建筑物边缘轮廓图传递给GPU端中的洪水模拟模块。维持一张代表各个建筑物周围最大水位的纹理图BuildingWaterDepth(按照建筑物的标号顺序存储)，实时更新洪水状态图的同时，若一个线程计算的网格属于某一建筑物的外部轮廓(通过建筑物轮廓图来判断)，则将该网格的水位深度填入BuildingWaterDepth纹理的相应位置，这其中涉及到互斥访问，所以需要通过原子操作InterLockedMax来执行。直接访问BuildingWaterDepth纹理，结合自身建筑物有关的参数和上文提到的损失函数计算公式，输出相应预警等级对应的色彩。
图10添加挡板前后效果图对比：(a)为添加前的洪水淹没建筑预警图，可以看到中心建筑物被洪水淹没；(b)为添加后的阻洪效果，可以看到中心建筑物未被洪水淹没。在模拟过程的实时交互中修建挡板的操作的具体实现：CPU端通过用户界面交互传输给GPU起点、终点、宽度、高度四个参数，GPU内部每一个线程会计算自己所代表的网格是否落在相应的矩形内，若在内部则相应的代表地形高度的纹理通道增加相应的高度。
以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。
Claims (5)
1.一种基于二维水动力学的大规模洪水场景模拟预警交互方法，其特征在于包括以下步骤：
1)将模拟场景在水平面上细分为均匀的网格，将水体高度和基础地形高度存储在网格中心，将速度存储在网格边界；采用基于浅水方程的二维欧拉网格法模拟洪水演进，将浅水方程转化为拉格朗日表达形式，得到初步建立的物理场；所述物理场包括速度场和高度场；
2)在交错网格上对物理场进行对流更新，并在更新对流项时采用大时间步长并引入体积守恒与动量守恒，得到基于体积守恒和动量守恒变形后的浅水方程；根据变形后的浅水方程，采用隐式积分法和雅克比迭代法求解更新后的高度场，根据更新后的高度场梯度进一步求解得到更新后的速度场；所述的步骤2)具体为：
2.1)基于体积守恒和动量守恒，在交错网格上对物理场进行对流更新：
在每个网格位置执行虚拟粒子的向后回溯来计算网格上的对流，对于需要更新的物理场，以物理位置p处的网格单元(i,j)为例，该粒子在上一个时间步长的位置为pbackward＝p-Δt·ν(p)，其中p为当前网格中心点的位置，Δt为时间步长，ν(p)为当前网格中心点通过差值获得的速度；
对物理场更新对流项时，先更新高度场，再更新速度场；将当前物理场F中的网格的物理量，按照一定比例加权到下一时间步物理场F′的网格中，更新方式为：F′＝Interpolate(F,pbackward)，其中Interpolate(·)表示插值，确保分配时F′中每个网格的物理量比例都为1，保证整体的体积守恒，具体步骤为：
S1：构建一个数组W，用于存储各个网格实际被分配的比例；
S2：执行对流计算，得到对应网格的回溯点，即该粒子在上一个时间步长的位置pbackward，根据双线性插值原则得到需要从周围四个网格提取的比例；
S3：针对比例权重大于1的部分做归一化处理，针对比例小于1的网格进行分散操作；所述的分散操作具体为：
将多余的物理量分散到对流后F′的网格中，计算映射关系qforward＝q+Δt·ν(q),其中q表示比例小于1的网格，ν(q)为网格q中心点通过差值获得的速度，qforward表示q映射后的网格，对qforward周围四个网格做双线性插值估算比例，按照估算后的比例将网格q中的多余物理量分散到四个网格中；
2.2)更新浅水方程形式：
更新对流项引入体积守恒和动量守恒后，浅水方程变形为：
其中，u代表x方向的速度，v表示y方向的速度，h表示水面总高度，b表示水底地形的高度，η＝h-b为水的深度，g表示重力加速度；上标n+1表示当前时间步长的物理量，n表示前一时间步长的物理量，
分别表示更新过对流项后的物理量；
2.3)更新高度场：
计算物理场中的扩散项，使用隐式积分法，根据更新后的浅水方程得到：
其中，
然后使用雅克比迭代法并行求解大型稀疏线性方程组hn+1 i,j＝f(hn+1 i-1,j,hn+1 i+1,j,hn +1 i,j-1,hn+1 i,j+1)，其中f为线性函数，下标(i,j)表示网格位置，从而得到更新后的高度场；
2.4)更新速度场：
根据步骤2.3)得到的高度场梯度，得到：
进一步求解得到更新后的速度场；
3)根据步骤2)更新后的物理场，采用动态自适应网格细分技术抬升渲染水面，采用纹理平铺技术实现速度场的水流绘制，并采用双三次B样条曲线插值法对绘制的洪水边缘做光滑处理；所述水流绘制包括单块平板单一流速的绘制、固定流场的绘制以及时变流场的绘制；
4)基于Hydrotec模型，根据淹没水深和建筑物类别设计洪灾损失评估函数，其中每一类建筑物对应一个洪灾损失评估函数；对洪灾损失评估函数值做归一化处理，并将归一化结果设为不同等级，每一个等级使用不同的颜色来表示建筑的受损程度，用于洪水场景绘制时实时改变建筑的外观颜色，实现可视化预警；
5)当收到预警信号时，选择在建筑物周围模拟修建挡板操作，调整挡板所设位置对应网格水底地形的高度，然后重复步骤2)-步骤4)完成洪水场景的实时更新，实现洪水模拟场景下的实时交互。
2.根据权利要求1所述的一种基于二维水动力学的大规模洪水场景模拟预警交互方法，其特征在于，所述的步骤1)具体为：
1.1)将模拟场景在水平面上细分为m*n均匀网格，为保证稳定性，使用交错网格表示物理量，将水体高度和基础地形高度存储在网格中心，将速度存储在网格边界；
1.2)将浅水方程写成拉格朗日形式：
其中，u代表x方向的速度，v表示y方向的速度，h表示水面总高度，b表示水底地形的高度，η＝h-b为水的深度，g表示重力加速度，等式左边为各物理量与速度场相关的对流项，等式右边为扩散项。
3.根据权利要求1所述的一种基于二维水动力学的大规模洪水场景模拟预警交互方法，其特征在于，所述的步骤3)具体为：
3.1)根据更新后的物理场，采用动态自适应网格细分技术抬升渲染水面：
计算水面高度变化程度
其中hi和hj表示相邻的两个网格对应的水面总高度，定义阈值Scap1和Scap2，将水面高度变化划分为三个等级：当
时，无需细分网格，当
时，中度细分网格，当
时，高度细分网格；
3.2)采用纹理平铺技术实现速度场的水流绘制：
a.单块平板单一流速的绘制：
水面流动由流速大小和方向决定，流速大小决定采样法线贴图的偏移幅度，流速方向决定偏移的角度；以(1,0)方向为默认方向，默认的偏移后的UV纹理坐标为(速度*时间步长,0)，根据水流方向生成对应的旋转矩阵：
其中，flowdir为水流方向矢量，flowdir.x表示flowdir在X方向上的分量大小，flowdir.y表示flowdir在Y方向上的分量大小；
根据旋转矩阵对纹理坐标进行旋转变化，再根据水流的速度计算出UV偏移：
uv＝mul(rotateMatrix,uv)-(flowSpeed·Δt)
其中，flowSpeed表示流速大小，uv代表采样坐标，mul(·)表示乘法运算；
所述的采样法线贴图具体为：为每个三角网格光栅化为的多个片段提供一个法线，利用一张2D纹理存储法线数据，进行采样和映射后可得到对应的法线向量，使得水面表现为由很多微小平面组成，得到波动效果；
b.固定流场的绘制：
根据所述的单块平板单一流速的绘制方法，利用多个平板平铺来表达整个流场，每块平板按照对应流场的流速分别绘制；为解决简单拼凑带来的间隙问题，引入更多的纹理块，并且将引入的纹理块在原有块的基础上以固定的偏移平铺展开，结果块的法线贴图由展开块的双线性插值得到；
c.时变流场的绘制：
为解决前后时间步长间速度方向大小轻微改变后引起的水面抖动和混乱问题，以10°一个刻度划分将本来连续变化的流速方向离散成36个固定方向的流场，将采样到的流速进行流速方向上的分类过滤后，再利用分类过滤后的方向进行法线贴图的偏移；将本来连续变化的速度大小，按照0.1m/s的刻度分割，将采样到的流速进行流速大小的分类过滤后，再利用分类过滤后的方向进行法线贴图的偏移；
3.3)采用双三次B样条曲线插值法对绘制的洪水边缘做光滑处理：
采用双三次B样条曲线作为三次样条插值的基函数：
对于待估算点p，其坐标为(x,y),提取出周围最靠近的16个网格节点，X方向坐标为x1,x2,x3,x4，Y方向坐标为y1,y2,y3,y4,先在X方向进行插值求出距待估算点p最近的a(x,y1),b(x,y2),c(x,y3),d(x,y4)四个网格点的物理场量值：
其中，m＝1,2,3,4分别对应a,b,c,d，x表示X方向上坐标数值，xt表示下标t对应的X方向上的坐标，ym表示下标m对应的Y方向上的坐标；f表示物理场的数值，f(x,ym)为坐标为(x,ym)的点的物理场数值，k为基函数；
然后在Y方向上对a,b,c,d四个点再次插值求出p点的数值：
其中，yt表示下标t对应的Y方向上的坐标，f(x,y)为p点物理场的数值；
通过步骤3.3)，扩展的16个邻域网格允许将插值样条线的起点和终点的导数与相邻的插值面片匹配，从而生成导数平滑连续的表面。
4.根据权利要求1所述的一种基于二维水动力学的大规模洪水场景模拟预警交互方法，其特征在于，所述的步骤4)具体为：
4.1)设计基于淹没水深和建筑物类别的洪灾损失评估函数：
基于Hydrotec模型，水深与灾害程度呈现平方根级别的正相关，其相关系数与建筑物的4个类别：办公楼、商业用地、住宅用地、基础设施相关，共得到四种建筑物的洪灾损失评估函数；
4.2)设计可视化预警：
根据步骤4.1)得到的四种建筑物的洪灾损失评估函数，将洪灾损失评估函数值归一化到[0,1]之间，每隔0.2取一个等级，每一个等级使用不同的颜色来表示建筑的受损程度，使得洪水场景实时绘制时改变建筑的外观颜色。
5.一种用于实现权利要求1所述的基于二维水动力学的大规模洪水场景模拟预警交互方法的预警交互系统，其特征在于，所述预警交互系统包括：
预处理模块：用于将模拟场景在水平面上细分为均匀的网格，构建初始的基于拉格朗日表达的浅水方程；
求解模块：用于对预处理模块构建的浅水方程进行求解，不断更新物理场；
绘制模块：根据求解模块得到的更新物理场参数，实现渲染水面、水流绘制、以及洪水边缘的光滑处理；
可视化预警模块：用于设计洪灾损失评估函数，并对评估函数结果归一化处理和分级，根据分级结果使用不同的颜色表示建筑的受损程度，实现可视化预警；
交互模块：用于用户和系统的交互，当用户收到预警信号时，选择在建筑物周围模拟修建挡板操作，并实时显示修建挡板后的可视化效果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010000075.5A CN111191372B (zh) | 2020-01-01 | 2020-01-01 | 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010000075.5A CN111191372B (zh) | 2020-01-01 | 2020-01-01 | 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111191372A CN111191372A (zh) | 2020-05-22 |
| CN111191372B true CN111191372B (zh) | 2021-12-21 |
Family
ID=70708037
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010000075.5A Active CN111191372B (zh) | 2020-01-01 | 2020-01-01 | 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111191372B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111915158A (zh) * | 2020-07-15 | 2020-11-10 | 云南电网有限责任公司带电作业分公司 | 一种基于Flood Area模型的暴雨灾害天气风险评估方法、装置及设备 |
| CN111931275B (zh) * | 2020-07-27 | 2022-03-29 | 南昌大学 | 一种尾矿库坝身渗透破坏引起的溃坝过程模拟方法 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4431545B2 (ja) * | 2006-03-24 | 2010-03-17 | 株式会社日立エンジニアリング・アンド・サービス | 氾濫シミュレーション装置及びプログラム |
| EP2390801A1 (en) * | 2010-05-26 | 2011-11-30 | Livermore Software Technology Corporation | Hybrid element enabling finite element/smoothed particle hydrodynamics coupling |
| US8831916B2 (en) * | 2011-05-05 | 2014-09-09 | Siemens Aktiengesellschaft | Simplified smoothed particle hydrodynamics |
| CN103258085B (zh) * | 2013-04-25 | 2015-10-28 | 交通运输部天津水运工程科学研究所 | 一种污水深海排放扩散器上升管数量确定的方法 |
| CN108021780B (zh) * | 2018-01-24 | 2021-08-13 | 国电南瑞科技股份有限公司 | 一种基于无规则非结构网格模型的山洪动态仿真方法 |
| CN110287590B (zh) * | 2019-06-24 | 2023-08-18 | 天津大学 | 基于算子分裂及改进半拉格朗日求解污染物传播的方法 |
-
2020
- 2020-01-01 CN CN202010000075.5A patent/CN111191372B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN111191372A (zh) | 2020-05-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Št'ava et al. | Interactive terrain modeling using hydraulic erosion | |
| Cordonnier et al. | Large scale terrain generation from tectonic uplift and fluvial erosion | |
| Zyserman et al. | Modelling morphological processes in the vicinity of shore-parallel breakwaters | |
| van Maanen et al. | Modelling the effects of tidal range and initial bathymetry on the morphological evolution of tidal embayments | |
| Szmytkiewicz et al. | Coastline changes nearby harbour structures: comparative analysis of one-line models versus field data | |
| CN108021780B (zh) | 一种基于无规则非结构网格模型的山洪动态仿真方法 | |
| CN110362925A (zh) | 一种包含库区的土石坝漫顶溃决洪水数值模拟方法 | |
| CN111191372B (zh) | 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统 | |
| Roudier et al. | Landscapes synthesis achieved through erosion and deposition process simulation | |
| CN114511995B (zh) | 一种基于cesium模型的洪水分级预警方法 | |
| Liu et al. | Dynamic visualisation of storm surge flood routing based on three‐dimensional numerical simulation | |
| Zhang et al. | Synthetic modeling method for large scale terrain based on hydrology | |
| CN114020943B (zh) | 流域水面混合绘制方法及系统、电子设备、存储介质 | |
| Guisado-Pintado | Shallow water wave modelling in the nearshore (SWAN) | |
| Woodruff et al. | Subgrid corrections in finite-element modeling of storm-driven coastal flooding | |
| Borthwick et al. | The shallow flow equations solved on adaptive quadtree grids | |
| Begmohammadi et al. | Subgrid surface connectivity for storm surge modeling | |
| El Rahi et al. | Numerical investigation of wave-induced flexible vegetation dynamics in 3D using a coupling between DualSPHysics and the FEA module of Project Chrono | |
| Gainza et al. | A process based shape equation for a static equilibrium beach planform | |
| Hagen et al. | An unstructured mesh generation algorithm for shallow water modeling | |
| Favaretto et al. | A model of coastal flooding using linearized bottom friction and its application to a case study in Caorle, Venice, Italy | |
| Do Carmo et al. | Near-shore sediment dynamics computation under the combined effects of waves and currents | |
| CN104517299A (zh) | 视频流体物理驱动模型恢复及重新仿真的方法 | |
| Pytel et al. | Self-organized approach to modeling hydraulic erosion features | |
| Krámer et al. | An adaptively refined, finite-volume model of wind-induced currents in lake neusiedl |
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 |