基于DEM模拟淹没区域随时间推演的实现研究

最近被逼着实现模拟闸口泄洪以后对下游淹没范围的影响程序，要求体现随着时间的推移，淹没区域发生变化的效果。于是自己琢磨了这个算法，因为不是水利水文相关专业出身，所以此处没有考虑到水汽蒸发、土壤、植被、岩石等要素，从地形一个因素上进行考虑，引入数字高程模型作为计算依据，该算法还存在很多疏漏之处，欢迎评论指正。

基本思路是将大坝地址作为初始点，不断向周围8领域扩散，搜索能够被淹没的像素点，标记被淹没。重复以上过程，直到整个dem区域计算完毕。

其中P(i,j)作为原始扩散点，搜索其周围八个像素，满足条件C则标记为淹没点，其中条件C定义如下:
1.其高程小于特定淹没高度
2.其高程小于P(i,j)的高程
以上两个条件满足其一即可。

被标记为淹没点的像素记录下来，作为下一轮的原始点继续进行八领域扩散计算。

为了体现时间推演的效果，该算法中添加了FloodLevel变量用于记录每个淹没点的淹没层级，比如起始点大坝地址的FloodLevel为1，其周围淹没点的FloodLevel为2，再次扩散的淹没点FloodLevel为3，依此类推，最后生成一张FloodLevel栅格图，其值大于0的则认为被淹没区域，而FloodLevel值不同代表不同时间点新增的淹没区：

从图中可以看出，该方法在平原地区扩散呈现方形边界，这并不符合实际水流的扩散方式，实际水流入平原应该呈现一种类似圆形的扩散图形。

这是比较理想的状态，即不考虑平原地区每个像素点之间高程的微小差别。

呈现方形扩散，而不是圆形扩散的原因是八领域扩散方法每次计算FloodLevel时都是在原始点的基础上加1，而没有考虑到每个点离中心点的距离不一样，而且每次计算时，中心点的位置发生了改变，而类似水波纹的同心圆状扩散的中心点位置是一样的。

算法改进的基本思路是每次计算FloodLevel时，不再是在上一个点的FloodLevel基础上加1，而是顺藤摸瓜找到最原始的点，然后计算距原始点的距离S，实际的FloodLevel为原始点的FloodLevel+S。

原始点的寻找

如下图所示，下游的所有点都应该把图中红色点位置作为原始点：

如下图所示，P1点的原始点应该是P2，而不是P3,原因是P1和P3之间隔着未被淹没的地区（黑色代表未被淹没区），水不可能越过未被淹没的点直线到达新的点，画面彩色部分的左上角为大坝位置（图中黄色直线标注）。

所以寻找原始点的方法为，搜索其父节点的原始点（父节点为扩充该像素的中心点），判断该像素和这个原始点中间是否隔着未被淹没的像素，如果没有则继续搜素这个原始点的原始点，如果有则把父节点作为原始点停止搜素，依此类推，直到找到大坝地址或者找到最终原始点。

判断地图上任意两个像素中间是否隔着未被淹没的像素

需要找到这两个像素点中间相隔的像素点，然后依此判定这些像素点是否被淹没，有一个未被淹没则返回false，退出循环。
此处认为八领域内所有像素都是和中心点相邻的，所以像如下图中红色的两个点认为相邻。

如下图所示，图中两个红色方格的中间相隔像素为黄色的方格：

其计算的方法为：
1.判断两个像素分别在x、y轴方向上相差的距离（即像素个数）
2.选择上一步中相差距离最大的那个方向作为基准方向
3.在基准方向上由起始点一个像素一个像素的移动，计算出对应的另一个方向的像素位置，两个方向的像素位置合起来定义了一个点的位置
4.重复步骤3得出所有中间点

代码请见：
https://blog.csdn.net/wqy248/article/details/87856722

计算结果体现了良好的效果：

上图中大坝起始点的第一处水流向拐弯处可以看出较上次计算有了明显的改观，水的前沿部分形状开始随着水流方向的改变而改变：