本发明涉及水利工程技术领域，特别是涉及一种非统一高精度曲面网格水流水质模拟及可视化方法和系统。
背景技术：
在对水流水质进行数值模拟过程中，计算域内地形的概化精度以及计算网格单元对水流水质状态变化的适应性对模拟结果有较大影响。实际工程中计算域内往往包含地形变化显著的狭窄河道区域、堤防部分，相对平坦的河道滩地及面积广大的洪泛区，采用统一网格建模，会无法平衡模型计算效率和精度。目前水流水质模拟中较少涉及复杂地形区域下垫面的处理，通常采用已有的前处理软件建模，这些软件采用的建模方法以反距离加权法为主，建模精度无法得到保证。多数的水流水质模拟计算多采用对计算完成后才进行可视化，而且以二维可视化居多，缺乏对计算过程的实时干预方法以及可视化缺乏真实环境衬托，无法进行精确空间定位，沉浸感不足，实时性不强。
技术实现要素：
针对上述问题，本发明提供了非统一高精度曲面网格水流水质模拟及可视化方法和系统，用于为洪水预报调度、水环境保护和水资源管理提供技术支撑。
为了解决上述问题，本发明公开了一种非统一高精度曲面网格水流水质模拟及可视化方法，包括步骤：
利用预先获取的河道地形的散点数据及研究区域边界和预先选取的预设空间网格尺度建立河道地形高精度曲面模型，对所述河道地形高精度曲面模型进行离散求解，生成初始的地形高精度曲面网格；
根据预先设定的地形网格加密层数、对所述地形高精度曲面网格的局部区域进行剖分处理，生成初始的地形高精度自适应曲面网格，并采用四叉树对所述地形高精度自适应曲面网格进行组织；
在所述地形高精度自适应曲面网格的基础上，离散求解二维水流水质模型，在求解所述二维水流水质模型计算过程中，确定网格调整的判定指标和阈值，根据水流水质实时状态，对满足条件的网格单元进行合并和剖分，动态调整所述地形高精度自适应曲面网格，并采用动态时间步长，使所述地形高精度只适应曲面网格适应模型高精度计算；
利用通过三维可视化软件和三维建模软件开发的三维虚拟仿真平台实现对所述二维水流水质模型的计算的实时交互，将模型计算结果与三维虚拟环境相结合，实现流速、水位、水深、水质的三维可视化。
优选的，所述根据预先设定的地形网格加密层数、对所述地形高精度曲面网格的局部区域进行剖分处理，包括：
对高程变化梯度较大的所述局部区域，利用局部曲率熵方法确定加密网格，生成n层嵌套网格。
优选的，所述河道地形高精度曲面模型中相邻网格之间的尺寸满足预设整数倍的关系。
优选的，所述离散求解二维水流水质模型，包括：
利用有限体积法对控制方程离散，并采用hllc格式计算数值通量，进一步利用muscl求解所述河道的水流水质方程。
优选的，所述实时交互过程中：
采用调整显示时间间隔和插值处理的方法来提高计算速度，采用管道传输与事件触发机制实现数据传输；
所述的流速三维可视化采用粒子系统的方式进行表达；
所述的水位、水深、水质的三维可视化，采用颜色映射的方式进行表达。
另外，还提供了一种非统一高精度曲面网格水流水质模拟及可视化系统，包括：
数据获取模块，用于获取河道地形的散点数据，所述散点数据包括地形数据、边界条件数据、卫星遥感数据、矢量数据、工程设计数据和纹理数据中的部分或全部；
曲面建模模块，用于利用所述散点数据及研究区域边界和预先选取的预设空间网格尺度建立河道地形高精度曲面模型，对所述河道地形高精度曲面模型进行离散求解，生成初始的地形高精度曲面网格；
分层建模模块，用于根据预先设定的地形网格加密层数、对所述地形高精度曲面网格的局部区域进行剖分处理，生成初始的地形高精度自适应曲面网格，并采用四叉树对所述地形高精度自适应曲面网格进行组织；
模型计算模块，用于在所述地形高精度自适应曲面网格的基础上，离散求解二维水流水质模型，在求解所述二维水流水质模型计算过程中，确定网格调整的判定指标和阈值，根据水流水质实时状态，对满足条件的网格单元进行合并和剖分，动态调整所述地形高精度自适应曲面网格，并采用动态时间步长，使所述地形高精度自适应曲面网格适应模型高精度计算；
虚拟仿真模块，用于利用通过三维可视化软件和三维建模软件开发的三维虚拟仿真平台实现对所述二维水流水质模型的计算的实时交互，将模型计算结果与三维虚拟环境相结合，实现流速、水位、水深、水质的三维可视化。
可选的，所述分层建模模块具体用于对高程变化梯度较大的所述局部区域，利用局部曲率熵方法确定加密网格，生成n层嵌套网格。
可选的，所述河道地形高精度曲面模型中相邻网格之间的尺寸满足预设整数倍的关系。
可选的，所述模型计算模块具体用于利用有限体积法对控制方程离散，并采用hllc格式计算数值通量，进一步利用muscl求解所述河道的水流水质方程。
可选的，所述虚拟仿真模块在实时交互过程中：
采用调整显示时间间隔和插值处理的方法来提高计算速度，采用管道传输与事件触发机制实现数据传输；
所述的流速三维可视化采用粒子系统的方式进行表达；
所述的水位、水深、水质的三维可视化，采用颜色映射的方式进行表达。
上述技术方案可以看出，本发明实施例提供了一种非统一高精度曲面网格水流水质模拟及可视化方法和系统，该方法和系统利用高精度曲面建模方法生成河道地形，生成地形的精度高于传统方法，有效避免了传统方法生成地形精度存在缺陷的不足；对高精度曲面建模生成网格进行局部区域分层加密，分层网格不仅能够对地形梯度变化大的区域进行加密，而且对重点区域所在的网格单元进行加密，能够提高局部区域的模拟精度；利用四叉树数据结构组织河道地形网格，在此基础上离散求解二维水流水质模型，可以根据地形梯度、水流水质时空变化特征，动态调整网格，并采用动态时间步长，提高了模拟精度和效率；通过三维仿真平台与水流水质模型的实时交互，实现了基于三维虚拟环境的水流水质的实时可视化仿真，能够直观地展现计算结果的动态过程，从而能够为洪水预报调度、水环境保护和水资源管理提供技术支撑。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的非统一高精度曲面网格水流水质模拟及可视化方法的步骤流程图；
图2为本发明实施例提供的分层曲面网格示意图；
图3为本发明实施例提供的四叉树结构组织示意图；
图4为本发明实施例提供的侧边、边和角点的示意图；
图5为本发明实施例提供的四叉树的等级结构；
图6为本发明实施例提供的子网格在父网格中所处位置的整数表示示意图；
图7为本发明实施例提供的网格位置及其各个相邻方向的数字表示；
图8为本发明实施例提供的模型计算与三维虚拟仿真平台实时交互流程图；
图9为本发明实施例提供的非统一高精度曲面网格水流水质模拟及可视化系统的结构框图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例一
图1为本发明实施例提供的非统一高精度曲面网格水流水质模拟及可视化方法的步骤流程图。
参照图1所示，具体包括步骤：
s101：根据河道地形数据生成高精度曲面网格。
利用预先获取的所述河道地形的散点数据及研究区域边界，选取预设空间网格尺度，并根据所述空间网格尺度建立河道地形高精度曲面模型；对所述河道地形高精度曲面模型进行离散求解，生成初始的地形高精度曲面网格m1。
高精度曲面模型如下：
根据曲面论基本定义(somasundaram,2005)，当曲面的第一基本量e,f,g，第二基本量l,m,n满足对称性，且e,f,g为正定时，e,f,g,l,m,n满足gauss-codazii方程组，全微分方程组在初始条件f(x,y)＝f(x0,y0)(x＝x0,y＝y0)下存在着唯一的解z＝f(x,y)。
曲面的第一类基本量可表达为：
第二类基本量可表达为：
其中，
选用中心差分格式，分别得到第一类基本量和第二类基本量的有限差分格式
第一类基本量的有限差分形式为：
第二类基本量的有限差分形式为：
式中，hx和hy为差分步长，
对方程中的fx,fy采用中心差分离散格式，对fxx,fyy,fxy采用高阶差分离散格式：
式中：h是计算格网尺寸，将以上两式带入式，可以得到hasm模型的最终求解方程组。
在采样点位置，hasm模型要求满足：
式中：a,b,c分别为方程组(1)左边的系数矩阵，可以由式(19)计算；d,q,p为方程组(1)右边的常数向量；z为待求解向量，n为迭代次数。
基于最小二乘优化理论，hasm最终转化为求解下述最小二乘问题：
式中，s和k分别为采样点系数矩阵和采样点的值，与传统hasm相同。通过引入权重系数λ，上述约束最小二乘问题转化为：
最优化问题，即式(22)等价于：
wz＝v(23)
w＝ata+btb+ctc+λ2sts(24)
v＝atd+btq+ctp+λ2stk(25)
式中，w为对称正定大型系数矩阵。
hasm计算最终转化为大型线性方程组式(23)的求解。一般需要采用迭代法进行线性方程组的求解，例如高斯-塞尔德迭代法、预处理共轭梯度迭代法、多重网格算法、自适应算法等。大型线性方程组计算出满足精度的解后，将解向量的元素值对应相应的空间位置，即得到河道高程值的空间分布曲面，将其记为河道高精度曲面网格m1。
s102：剖分生成四叉树组织的分层曲面网格。
在得到上述的河道高精度曲面网格后，设定地形网格剖分层数n，在所述曲面网格m1基础上，对局部区域进行剖分处理，生成初始的地形高精度自适应曲面网格m2。采用四叉树对分层自适应网格进行组织。
对所述曲面网格m1的局部区域进行剖分处理，包括：
对高程变化梯度较大局部区域，利用局部曲率熵方法确定剖分网格，生成n层嵌套网格。设定局部曲率熵的阈值rc，当某网格单元的曲率熵大于rc时，将该网格单元剖分n层，形成嵌套网格。
p点的局部曲率熵r利用下式计算:
式中，ri为任意点i处的局部曲率熵，k为p点临近范围内点的个数，临近范围为预设的范围，pi为i处的曲率概率分布，其计算公式为：
pj为j处的曲率概率分布，其计算公式为：
式中，ki为点i处的曲率，kj为i邻域内点j处的曲率，曲率为已知值。
对重点关注局部区域(如研究区边界、河道弯道处等)，确定剖分网格，生成n层嵌套网格。
为使数值计算保持稳定和降低模型的复杂度，通过局部剖分生成的初始自适应曲面网格m2中相邻网格之间的尺寸满足1倍、2倍或1/2倍的关系。如图2所示。
对于生成的初始自适应曲面网格m2，采用四叉树结构进行组织存储，以便后面在求解水流水质模型时实现网格的快速剖分和合并。四叉树是一种树状数结构，在每一个节点上会有四个子区块。
图2所示的自适应网格的网格形状为正方形，属于非结构网格，它是依据一定的剖分规则，在父网格基础上递归剖分为四个相等的正方形而形成的，符合四叉树结构特征。同时满足相邻单元的尺寸倍数为1/2、1或2的关系。
四叉树网格属于结构网格，网格形状为正方形。每个正方形网格对应四叉树的一个结点，可以有上一层父网格和下一层子网格(如图2和3所示，ne、nw、sw和se分别表示根网格的4个象限)，尺寸最大的网格所在层为0层，对应的节点为根结点，没有子网格的为末端的叶网格。每个0层网格所包含的全部网格的集合称为一个四叉树结构，若干个四叉树形成的网格结构群称为四叉树森林。每个四叉树网格有4个角点，相邻角点的连线称为网格的侧边，而侧边上相邻网格节点的连线称为网格的边(图4)。只包含一条边的侧边称为最小公共边。包含相同公共边的2个网格互为邻居。为叙述方便，引入指针来说明建立四叉树网格结构的基本思路。每个网格有4个指针指向其4个子网格，1个指针指向其父网格。同时，每个网格还有分别指向其4个角点和4个侧边的指针。为了方便寻找邻居，每个网格节点有4个分别指向其4个象限处叶网格的指针，很明显，t型节点会有2个指针指向同一个叶网格。每个最小公共边有2个指针指向相邻的叶网格。如图5所示，网格ⅱ、ⅲ、ⅳ和ⅴ有共同的父网格ⅰ，网格ⅳ有4个子网格ⅵ、ⅶ、ⅷ和ⅸ。网格ⅰ和子网格ⅳ有共同角点c。若网格ⅱ、ⅲ、ⅴ和ⅶ是叶网格，网格节点e的4个指针将分别指向这4个网格。若网格ⅸ和ⅴ是叶网格，网格ⅸ的右侧边l是最小公共边，则最小公共边l的2个指针分别指向ⅸ和ⅴ。在遍历四叉树网结构时，采用链表连接同层的叶网格十分方便。
为提高自适应网格调整效率，利用邻居查询算法来减少数据的存储量，提高网格生成效率；但是已有的算法只是沿着四叉树上寻和下寻的过程中确定邻居，而没有考虑任一网格在其父网格中的具体位置的差异，没有对12个邻居方向作出规律性区分而同一对待，从而造成网格查寻时判断较多，效率不高。
本发明利用“方向数字化”邻居查寻算法。由于任一网格都是由最初的根网格划分而来，且它们都具有在各自父网格中的特定位置，如图4所示，因此在网格邻居查寻时应将这种特定性考虑进去；在任一网格的8个邻居法向(依据本发明的规定0)中，其规律性可总结如图5所示。
其网格本身所处其父网格的位置及其各个相邻方向以数字表示，用(1,1)，(1,0)，(0,1)和(0,0)表示指定网格所处其父网格中的位置，这样用两个数组存储这些信息以确定该网格相对于根网格的位置。边相邻的特征定义为x或y，取定值为0或1；角相邻的方向以数字代替，如se为(1,1)，nw为为(0,0)。
针对任一既定网格，只需比较此网格的位置数组和所需查寻邻居的数字方向表示即可完成查寻。
边邻居查寻时，通过比较给定网格的位置数组和查寻方向的数字特征表示上寻，找到能同时包含此网格和所需查找邻居的最小父网格，再按照存储的上寻路径以一定的特征变化反向下寻，直至达到规定的边邻居。
角邻居查寻相对复杂，但有了边邻居查寻算法和方向的数字表示，查寻变得直观、方便。考虑到任一网格在其父网格中的位置差异，并从此点出发，将四个方向的角邻居查寻分为四种不同位置网格的四种不同情况；其中，针对任一种位置网格，记录一个自然生成的角邻，另两种情况考虑如比表1，其基本步骤为：首先查寻该网格父网格的既定边邻(显然，如果此边邻居不存在，则所需查找的角邻居亦不复存在。)，再沿着此边邻居下寻一步即可到达，如果此边邻居没有子网格，则此边邻居即为所寻。
表1角邻居查寻的步骤(考虑两种情况)
针对任一种位置网格只剩下一种情况，即1网格的(0，0)方向角邻居查寻；2网格的(0，1)方向角邻居查寻；3网格的(1，0)方向角邻居查询；4网格的(1，1)方向角邻居查寻(如图4和图5所示)。考虑到这种情况的特殊性，以1网格的(0,0)方向角邻居查寻为例，步骤如下：
(1)沿着该网格相对于根网格的二进制信息数组上寻，并与方向数(0,0)比较，同时记入上寻次数；
(2)如果该网格的父网格的位置数未发生变化，则继续上寻(如果一直上寻到根网格而位置数仍未变化，则同样要作适当位置坐标判断)；否则，记入此时父网格的自然角邻居，以备下寻之用；
(3)下寻次数取执行完第1步所记入的次数，下寻得子网格位置数为(1,1)；同样，如果下寻的当前网格的子网格不存在，则当前网格即为所寻。
同理，另外三种位置网格亦可按以上步骤做适当修改，完成此最后一种情况的角网格邻居查寻。角邻居查寻以边邻居查寻算法为基础，将给定网格分为两种情况考虑(即在同一父网格中的四个子网格的位置特征)；在任一网格的四个方向的角邻居查寻过程中，判断次数减少，所需查寻的最多步骤数分别为1步、2步、3步和网格最大划分层数。
s103：基于自适应分层曲面网格的水流水质模型计算。
在自适应曲面网格m2的基础上，离散求解二维水流水质模型。在模型计算过程中，确定网格调整的判定指标和阈值，根据水流水质实时状态，对满足条件的网格单元进行合并和剖分，动态调整自适应曲面网格m2，并采用动态时间步长，使其适应模型高精度计算。
所述离散求解二维水流水质模型，包括：
利用有限体积法对控制方程离散，并采用hllc格式计算数值通量，进一步利用muscl求解所述河道的水流水质方程。
本发明利用的二维水流水质模型的控制方程以向量表示为：
式中，u为守恒向量；e为对流项的通量向量，e＝(f,g)；f、g分别为x、y方向的对流通量；s为源项向量，包括床面底坡项和摩阻项等。
向量u、f、g和s具体表示为：
式中，η定义为基准面上的水位，水深h＝η-zb，zb为底部高程，u和v分别为沿x和y方向的深度平均速度分量，g为重力加速度，ρ为水密度，τbx和τby为床面摩擦应力，
所述的网格调整的判定指标和阈值，包括：
水流判定指标θw是由水位梯度、水深梯度、弗劳德数组成的综合指标,用下式表示：
式中，h和η分别为水深和水位，θη为水位梯度、θh为水深梯度，fr为流态弗劳德数，限制在[0,1]区间内，即fr＝max{0,min{1,v/(gh)1/2}},v为网格单元的绝对速度，g为重力加速度。
水质判定指标θc是由水质浓度变化梯度指标，用下式表示：
水流判定指标阈值上、下限分别用
水质判定指标阈值上、下限分别用
所述的动态调整自适应网格m2，包括
当某网格单元水流判定指标θw大于给定阈值
所述的动态时间步长，根据cfl条件来确定。cfl条件表达形式为：
在模型计算中，一般ncfl可以取值为0.8。
s104：基于三维虚拟环境的水流水质可视化仿真。
利用三维可视化软件和三维建模软件，开发三维虚拟环境平台。通过二维水流水质模型计算与三维虚拟环境平台实时交互，将模型计算结果与三维虚拟环境相结合，实现流速、水位、水深、水质的三维可视化。
地形地物建模，本发明实施例选用terravista进行地形建模，multigencreator进行地物建模。terravista软件能集成数字高程模型、遥感影像数据和矢量数据一次性生成大范围流域地形。对于其中的关键地物，如大坝、闸门等水利工程，利用multigencreator在获取地物模型的工程设计图和纹理图片基础上建模生成。
三维虚拟环境平台，本发明实施例利用开发平台visualstudio2010和三维可视化软件包openscenegraph开发三维可视化平台，将建好的地形地物三维模型导入到三维可视化平台中，进行视点设置、光照参数设置、纹理映射参数设置、三维可视化模型调度参数设置等相关步骤，生成三维虚拟环境平台。
所述的三维可视化软件，选择openscenegraph、opengvs、skyline等其中一个。本实施例选择openscenegraph。
所述的三维建模软件，选择multigencreator、terravista、terrabuilder、3dmax等其中一种或两种。本实施例选择multigencreator和terravista。
所述的实时交互方式，采用调整显示时间间隔和插值处理的方法来提高计算速度，采用共享内存方式实现计算进程与显示进程的数据传输。具体流程如图所示。
对计算二维流场的水流水质模型而言，两次计算的时间tc大约以数10s计，而三维场景的刷新显示频率要在25帧/s以上时，才具有良好的视觉效果。假定系统循环一次的时间为t，将两次计算结果按m步等间距内插，内插后数组的时间步长δt＝t/m，选取适当的速度比尺，将每一步的速度矢量场按空间位置插值，使每一条速度矢量线按比尺前进，并以(1/25)s的速度刷屏，直至累计时间接近于δt，转入下一步的速度矢量场并插值和显示。完成m步循环，这时的显示总时间t已大于一次计算所需的时间tc，即后台的数学模型已完成一次计算，等待显示平台调入下一次的计算结果，将新的流场数据分步内插、空间插值和实时显示，继续循环，完成大计算量的模拟计算与三维显示同时并存下的实时显示功能。
以系统运行时间为自变量，在三维可视化程序中用三个数组previous，current，next依次存储3次计算结果。在previous和current数组之间插值得出的流场过程作为当前显示的数据。当系统运行一个周期t后，数据进行更新传递，previous＝current，current＝next，此时向计算程序发出下一轮计算的指令，等计算完毕后将结果存入next数组。在此期间显示仍然以previous和current插值所得的数据继续运行。如此反复，就可以实现计算与显示的同步，系统的运行和数据、指令传递的关系如图所示。当然这种同步并非真正意义上的同步，当计算程序进行第三次计算时，显示程序显示的却是第一、二次计算的插值结果，不过这种滞后对流域模拟是可以接受的。
本发明实施例采用事件触发机制完成计算与显示的指令传递与控制，在三维可视化程序中创建内核对象事件(event)，在数模计算程序中打开该事件，这样就获得了该事件的使用权，从而使事件成为计算和显示程序的公共变量。在计算程序的每个计算步长起始处，加入等待事件触发的语句waitforsingleobject，这样当无事件触发时，计算程序挂起，不进行计算；当系统程序需要进行下一步计算，发出事件触发时，计算程序中由于事件处于触发状态，程序开始执行。同样，系统程序接收数据也用同样的方法，当计算完成并将结果输出后，发送事件给系统程序，系统程序收到触发后读取数据，存入next数组。
数据传输，三维可视化程序采用visualstudio开发。对数学模型研究，大量的计算程序主要是采用fortran编写。需要实现计算与显示进程间的数据传输。两个进程间数据通讯方式包括共享内存、管道、消息队列等方式。本发明实施例采用共享内存的方法。
共享内存允许两个或多个进程共享一个给定的存储区，这一段存储区可以被两个或两个以上的进程映射至自身的地址空间中，一个进程写入共享内存的信息，可以被其他使用这个共享内存的进程，通过一个简单的内存读取操作读出，从而实现进程间的通信。
采用共享内存进行通信的一个主要好处是效率高，因为进程可以直接读写内存，而不需要任何数据的拷贝，对于像管道和消息队列等通信方式，则需要在内核和用户空间进行四次的数据拷贝，而共享内存则只拷贝两次：一次从输入文件到共享内存区，另一次从共享内存到输出文件。如图所示。
一般而言，进程之间在共享内存时，并不总是读写少量数据后就解除映射，有新的通信时，在重新建立共享内存区域；而是保持共享区域，直到通信完毕为止，这样，数据内容一直保存在共享内存中，并没有写回文件。共享内存中的内容往往是在解除映射时才写回文件，因此，采用共享内存的通信方式效率非常高。
所述的流速三维可视化，采用粒子系统的方式进行表达。
粒子系统利用openscenegraph提供的粒子系统类来通过粒子的产生、运动、消亡和绘制的全生命周期的管理来实现对流场的可视化模拟。osg提供了专门的粒子系统工具，如osgparticle::particle类是粒子模板类，生成粒子模板，决定粒子的大小、颜色、生命周期等，osgparticle::particlesystem决定粒子的总体属性，包括粒子总数、粒子纹理等，那么由osgparticle::program类实现对粒子的操作，即进行粒子运动的控制。本发明实施例通过设计一个继承osgparticle::operator类的flowfield-operator类，重写其operate函数，实现对粒子运动的控制。在流场模拟中，粒子需要表现流动特性，涉及空间位置、矢量大小、矢量方向随时间的变化过程等，通常用移动的箭头表示，箭头的长短代表流速大小，箭头指向代表流速的方向。要通过插值得出粒子响应位置的流速，还需确定粒子所处计算网格的编号信息，另外为判断粒子的存亡，需要加入粒子是否跑出模拟区域的标识变量。
粒子的产生主要是基于水流水质模型计算网格单元来限制，每个网格单元产生一个粒子。
粒子产生之后，粒子的运动由水流水质模型计算控制。由于水流水质模型计算的流速场为欧拉场，为了更加直观可视化表现流体运动特征，需要将其转化拉格朗日场，模拟流体质点随时间的运动过程。
设粒子在t时刻的位置为x(t)，运动速度为v(x)＝dx(t)/dt。经过δt时段后，粒子位移到新的位置：x(t+δt)＝x(t)+∫v(x(t))dt。利用改进欧拉算法计算上述积分方程，可以达到二阶精度。先用欧拉法求粒子的预测位置x*(t+δt)，即x*(t+δt)＝x(t)+v(x(t))·δt，然后利用x*估计t+δt时刻粒子的速度值，与原位置速度值求平均，得到t+δt粒子位置x(t+δt)：x(t+δt)＝x(t)+[v(x(t))+v(x*(t+δt))]·δt/2。判断粒子运动后所在网格单元的节点流速值，可以插值得出粒子的流速值和流速方向。空间插值方法有反距离权重法、最小曲率法等方法，本发明实施例选用反距离权重方法。所有当前时刻的粒子的属性值求出之后，就可以绘制研究区域的粒子系统，实时显示在屏幕上。
为了平衡粒子分布，当粒子运动出研究区域外，这个粒子就死亡了。如果一个网格单元内粒子数大于1，则保留1个粒子，其余粒子置于死亡。由于粒子运动及调整，可能导致部分网格单元内没有粒子，在没有粒子的网格单元内重新生成粒子，粒子属性可通过空间插值获取。
水位、水深、水质的三维可视化，采用颜色映射的方式进行表达。本发明实施例以水质三维可视化为例说明，绘制水面需要将离散的四边形单元转化为三角形单元，颜色映射就需要给三角形顶点的着色，颜色值的确定需要依据水质浓度值的大小。计算出各顶点的颜色值后，openscenegraph内置函数就会自动计算出三角形表面各点的颜色值。本发明实施例确定颜色的方式如下：
首先建立一个颜色表，颜色表可由上百种渐变颜色组成(或者设计成离散色)，并给这些颜色依次编号为0，1，2，3，…，对于第i个顶点，显示时先取出它的水质浓度值ci代入下式计算：
式中，cmin为水质浓度最小值，cmax为最大高程值，num为颜色表上的颜色总数，int表示取整，求出k值后，以颜色表中第k种颜色为第i个顶点的颜色。建立颜色表时，先确定最大浓度值、中间浓度值和最小浓度值对应的颜色，然后用分段函数内插出其他水质浓度值对应的颜色。
假设颜色表上的颜色总数为num，最小浓度值cmin对应的颜色为clr1＝(r1,g1,b1)，中间浓度值对应的颜色为clr2＝(r2,g2,b2)，最大浓度值cmax对应的颜色为clr3＝(r3,g3,b3)，从c1到c2有m种颜色，c2到c3有num-m种颜色，则颜色表上第i种颜色值clri的计算公式为：
从上述技术方案可以看出，本实施例提供了一种非统一高精度曲面网格水流水质模拟及可视化方法，该方法利用高精度曲面建模方法生成河道地形，生成地形的精度高于传统方法，有效避免了传统方法生成地形精度存在缺陷的不足；对高精度曲面建模生成网格进行局部区域分层加密，分层网格不仅能够对地形梯度变化大的区域进行加密，而且对重点区域所在的网格单元进行加密，能够提高局部区域的模拟精度；利用四叉树数据结构组织河道地形网格，在此基础上离散求解二维水流水质模型，可以根据地形梯度、水流水质时空变化特征，动态调整网格，并采用动态时间步长，提高了模拟精度和效率；通过三维仿真平台与水流水质模型的实时交互，实现了基于三维虚拟环境的水流水质的实时可视化仿真，能够直观地展现计算结果的动态过程，从而能够为洪水预报调度、水环境保护和水资源管理提供技术支撑。
需要说明的是，对于方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明实施例并不受所描述的动作顺序的限制，因为依据本发明实施例，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明实施例所必须的。
实施例二
图2为本发明实施例提供的非统一高精度曲面网格水流水质模拟及可视化系统的结构框图。
参照图2所示，本实施例提供的非统一高精度曲面网格水流水质模拟及可视系统具体包括数据获取模块10、曲面建模模块20、分层建模模块30、模型计算模块40和虚拟仿真模块50。
数据获取模块用于获取河道地形的散点数据，该散点数据包括地形数据、边界条件数据、卫星遥感数据、矢量数据、工程设计数据和纹理数据中的部分或全部。
曲面建模模块用于根据河道地形数据生成高精度曲面网格。
利用预先获取的所述河道地形的散点数据及研究区域边界，选取预设空间网格尺度，并根据所述空间网格尺度建立河道地形高精度曲面模型；对所述河道地形高精度曲面模型进行离散求解，生成初始的地形高精度曲面网格m1。
选用中心差分格式，分别得到第一类基本量和第二类基本量的有限差分格式
一般需要采用迭代法进行线性方程组的求解，例如高斯-塞尔德迭代法、预处理共轭梯度迭代法、多重网格算法、自适应算法等。大型线性方程组计算出满足精度的解后，将解向量的元素值对应相应的空间位置，即得到河道高程值的空间分布曲面，将其记为河道高精度曲面网格m1。
分层建模模块30用于剖分生成四叉树组织的分层曲面网格。
在得到上述的河道高精度曲面网格后，设定地形网格剖分层数n，在所述曲面网格m1基础上，对局部区域进行剖分处理，生成初始的地形高精度自适应曲面网格m2。采用四叉树对分层自适应网格进行组织。
对所述曲面网格m1的局部区域进行剖分处理，包括：
对高程变化梯度较大局部区域，利用局部曲率熵方法确定剖分网格，生成n层嵌套网格。设定局部曲率熵的阈值rc，当某网格单元的曲率熵大于rc时，将该网格单元剖分n层，形成嵌套网格。
对重点关注局部区域(如研究区边界、河道弯道处等)，确定剖分网格，生成n层嵌套网格。
为使数值计算保持稳定和降低模型的复杂度，通过局部剖分生成的初始自适应曲面网格m2中相邻网格之间的尺寸满足1倍、2倍或1/2倍的关系。如图2所示。
对于生成的初始自适应曲面网格m2，采用四叉树结构进行组织存储，以便后面在求解水流水质模型时实现网格的快速剖分和合并。四叉树是一种树状数结构，在每一个节点上会有四个子区块。
图2所示的自适应网格的网格形状为正方形，属于非结构网格，它是依据一定的剖分规则，在父网格基础上递归剖分为四个相等的正方形而形成的，符合四叉树结构特征。同时满足相邻单元的尺寸倍数为1/2、1或2的关系。
本发明利用“方向数字化”邻居查寻算法。由于任一网格都是由最初的根网格划分而来，且它们都具有在各自父网格中的特定位置，如图4所示，因此在网格邻居查寻时应将这种特定性考虑进去；在任一网格的8个邻居法向(依据本发明的规定0)中，其规律性可总结如图5所示。
同理，另外三种位置网格亦可按以上步骤做适当修改，完成此最后一种情况的角网格邻居查寻。角邻居查寻以边邻居查寻算法为基础，将给定网格分为两种情况考虑(即在同一父网格中的四个子网格的位置特征)；在任一网格的四个方向的角邻居查寻过程中，判断次数减少，所需查寻的最多步骤数分别为1步、2步、2步和网格最大划分层数。
模型计算模块用于基于自适应分层曲面网格的水流水质模型计算。
在自适应曲面网格m2的基础上，离散求解二维水流水质模型。在模型计算过程中，确定网格调整的判定指标和阈值，根据水流水质实时状态，对满足条件的网格单元进行合并和剖分，动态调整自适应曲面网格m2，并采用动态时间步长，使其适应模型高精度计算。
虚拟仿真模块用于基于三维虚拟环境的水流水质可视化仿真。
利用三维可视化软件和三维建模软件，开发三维虚拟环境平台。通过二维水流水质模型计算与三维虚拟环境平台实时交互，将模型计算结果与三维虚拟环境相结合，并输出到相应的显示设备上实现流速、水位、水深、水质的三维可视化。
从上述技术方案可以看出，本实施例提供了一种非统一高精度曲面网格水流水质模拟及可视化系统，该系统利用高精度曲面建模方法生成河道地形，生成地形的精度高于传统方法，有效避免了传统方法生成地形精度存在缺陷的不足；对高精度曲面建模生成网格进行局部区域分层加密，分层网格不仅能够对地形梯度变化大的区域进行加密，而且对重点区域所在的网格单元进行加密，能够提高局部区域的模拟精度；利用四叉树数据结构组织河道地形网格，在此基础上离散求解二维水流水质模型，可以根据地形梯度、水流水质时空变化特征，动态调整网格，并采用动态时间步长，提高了模拟精度和效率；通过三维仿真平台与水流水质模型的实时交互，实现了基于三维虚拟环境的水流水质的实时可视化仿真，能够直观地展现计算结果的动态过程，从而能够为洪水预报调度、水环境保护和水资源管理提供技术支撑。
对于装置实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
本说明书中的各个实施例均采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似的部分互相参见即可。
尽管已描述了本发明实施例的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明实施例范围的所有变更和修改。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者终端设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者终端设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者终端设备中还存在另外的相同要素。
以上对本发明所提供的技术方案进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。