准备工作1——洼地识别

       高程数据是淹没深度计算的基础，我最初的设想是基于Arcgis Pro的流向算法（D8和MFD）先将研究区域内的流向计算出来，然后根据输入的水量数据（如雨量、溢流体积或地表径流分布）进行深度模拟。但是要运行流向算法要先进行填挖操作，对于只流向坡度下降最陡方向的D8算法，如果存在洼地，则该像元会被赋予一个编号，该编号意味着该洼地没有流向。对于MFD算法，洼地会被强行赋予一个流向，这又导致了回流问题。因此我决定采用洼地识别+洼地淹没和溢流+动态流向计算的框架，先采用Arcgis Pro的填挖工具，得到填挖后的高程，再进一步通过栅格计算器提取出洼地所在的空间像元。备注：如果填挖工具输出的填挖后的高程的无法显示，那大概率是输入的高程数据有问题，先检查高程数据是否有负值，如果有可以通过栅格计算器将最小值变为零，再导出为栅格数据，注意定义数据类型和无效值（我定义为：0）。如果填挖失败或者实在找不到原因也没有关系，因为代码中提供了一个填挖算法，只需指定高程数据和对应的无效值，便可以生成填挖后的高程数据。

inundation = Simulation(100, 141)
# 读取栅格
dem, profile = inundation.load(r'em.tif')
# 填挖
fill_dem = inundation.fill_in_surface_raster(dem, 0)
# 保存
inundation.save_tif(r'fill_dem.tif', fill_dem, profile)

# 输出
填洼完成!
总处理像素: 11708788
填洼像素数: 5657723
填洼后高程范围: 0.00 - 1707.00
已保存: fill_dem.tif

准备工作2——水量数据

       我计划将SWAT模型的地表径流模拟模块与流向驱动的淹没深度快速评估方法进行耦合，但是这需要理清一个前提，即SWAT模型输出的模拟结果是什么？实际上SWAT模型的地表径流模拟结果（地表径流深度）反映了特定降雨事件后，单位面积上产生的净地表产流量，通常以毫米（mm）为单位，其计算基于总降水量扣除植被截留、洼蓄、下渗及蒸发等损失后的剩余部分。这一过程主要受降雨强度、土壤入渗能力、土地利用/覆被以及前期土壤湿度等因素控制，其生成机制可分为超渗产流（当降雨强度超过土壤入渗速率时）和蓄满产流（当土壤剖面完全饱和后）两大类。相比之下，淹没深度则是指洪水或地表积水在某一具体位置的实际垂直水层厚度，同样以毫米或米为单位，它是地表径流经过坡面汇流、沟道传输、地形调蓄、排水系统干预及可能的外部边界条件（如潮汐、河流水位顶托）共同作用后的最终空间分布结果。在理想化的无约束平地且忽略所有损失的情况下，若全部径流完全蓄积于原地，则二者数值可能趋近。

       所以SWAT模型的地表径流模拟结果（地表径流深度）在一定情况下（不考虑地表径流会沿地形坡度迅速汇集，通过地表漫流、排水管网或天然河道向下游输送）可以当作淹没水深去进行后续的洪涝危险性评估。那为什么要研究流向驱动的淹没深度快速评估方法呢？因为从地表径流到最终淹没深度的转化，是一个受多重物理过程和边界条件调控的复杂水动力学问题。准确评估洪水风险，必须超越对单一产流指标的关注，转而构建能够整合产流、汇流与淹没全过程的综合性模拟框架。因此，我打算先从淹没深度的快速评估方法入手，对现有的基于高程的淹没深度快速评估方法进行改进，后续再转向对水动力原理的研究。

       我的水量输入数据是SWAT模型的地表径流模拟结果，我需要确保它和高程数据的栅格数目一致（空间精度对齐），因此我需要先进行重采样，再进行按掩膜提取，最后导出数据时记得定义无效值。

流向驱动的淹没深度快速评估方法

       评估过程基于优先队列，会对地表漫流和洼地淹没进行分开模拟，其核心框架都储存在类Simulation中，包括：

• 洼地的连通域识别（连通的洼地），为洼地模拟做准备，输入是diff.tif(填挖后与填挖前的高程数据之差，无效数据被定义为0)，输出为labeled.tif(连通域结果，相互连通的洼地都被标记为唯一的颜色)

# Simulation需要两个参数
像元大小CELL_SIZE（代表像元的大小，一般假设像元是正方形，即X=Y）
斜距slantDistance(代表对角像元中心点之间的距离)
inundation = Simulation(100, 141)
dem, profile = inundation.load(r'dem.tif')
diff, _ = inundation.load(r'data\diff.tif')
labeled_array, num_labels = inundation.label_filled_pixels(diff, 0)
inundation.save_tif(r'data\labeled.tif', labeled_array, profile)
# 输出
dem.tif加载成功，形状:(1697, 1359)
data\diff.tif加载成功，形状:(1697, 1359)
共找到 19019 个填挖连通域
计算连通域: 100%|██████████| 19019/19019 [00:33<00:00, 568.13it/s]
已保存: data\labeled.tif

• 读取水量数据，并提取其中的有效数据及其空间像元位置，作为迭代的对象，get_valid_data_mask(输入数据, 无效值)定义了一个简单方法，通过大于无效值提取有效值空间位置。水量数据（input_list）需采用列表形式，列表的元素均为三元元组，分别代表地表径流深度、所在栅格的行和列

suf_runoff, _ = inundation.load(r'hrus.tif')
vaule, row, col = inundation.get_valid_data_mask(suf_runoff, 0)
input_list = (list(zip(vaule, row, col))

• 计算区域内的最小坡度和最大坡度，用于随坡度变化的流量分配系数的计算

inundation.min_max_slope(dem,0)
# 输出
坡度最大值：1.8700
坡度最小值：0.0071

• 最后是进行淹没深度的快速评估，输入数据和参数说明如下：

inundation.water_flow(list(zip(vaule, row, col)), dem, labeled_array, 0, (0.007092198356986046, 1.8700000047683716), volume_loss=0.001, min_volume=1)
inundation.simulate_water_flow()
# 保存最大淹没水深和最终淹没水深
water_max_depth = inundation.water_max_depth
inundation.save_tif(r'water_max_depth_test2604202320.tif', water_max_depth, profile)
water_final_depth = inundation.water_final_depth
inundation.save_tif(r'water_max_depth_final2604210020.tif', water_final_depth, profile)

• 完整运行代码如下：分为初次运行和二次运行

# 第一次运行
# 初始化
inundation = Simulation(100, 141)
# 读取高程数据
dem, profile = inundation.load(r'dem.tif')
# 读取填挖前后的高程差
diff, _ = inundation.load(r'data\diff.tif')
# 计算洼地的连通域
labeled_array, num_labels = inundation.label_filled_pixels(diff, 0)
# 保存下来，下次直接通过labeled_array, profile_labeled = inundation.load(r'data\labeled.tif')读取
inundation.save_tif(r'data\labeled.tif', labeled_array, profile)
# 读取地表径流栅格
suf_runoff, _ = inundation.load(r'hrus.tif')
# 提取有效数据
vaule, row, col = inundation.get_valid_data_mask(suf_runoff, 0)
# 计算区域的坡度最值，运行较慢，记住值后下次无需运行
min_slope, max_slope = inundation.min_max_slope(dem,0)
# 输入参数
inundation.water_flow(list(zip(vaule, row, col)), dem, labeled_array, 0, (min_slope, max_slope), volume_loss=0.001, min_volume=1)
#开始运行
inundation.simulate_water_flow()
# 保存最大淹没水深和最终淹没水深
water_max_depth = inundation.water_max_depth
inundation.save_tif(r'water_max_depth_test2604202320.tif', water_max_depth, profile)
water_final_depth = inundation.water_final_depth
inundation.save_tif(r'water_max_depth_final2604210020.tif', water_final_depth, profile)

# 第二次运行
inundation = Simulation(100, 141)
# 读取高程数据
dem, profile = inundation.load(r'dem.tif')
# 读取填挖前后的高程差
diff, _ = inundation.load(r'data\diff.tif')
# 读取洼地连通域
labeled_array, profile_labeled = inundation.load(r'data\labeled.tif')
# 读取地表径流栅格
suf_runoff, _ = inundation.load(r'hrus.tif')
# 提取有效数据
vaule, row, col = inundation.get_valid_data_mask(suf_runoff, 0)
# 运行
inundation.water_flow(list(zip(vaule, row, col)), dem, labeled_array, 0, (0.007092198356986046, 1.8700000047683716), volume_loss=0.001, min_volume=1)
inundation.simulate_water_flow()
# 保存最大淹没水深和最终淹没水深
water_max_depth = inundation.water_max_depth
inundation.save_tif(r'water_max_depth_test2604202320.tif', water_max_depth, profile)
water_final_depth = inundation.water_final_depth
inundation.save_tif(r'water_max_depth_final2604210020.tif', water_final_depth, profile)

参考文献：
[1]秦承志,李宝林,朱阿兴,杨琳,裴韬,周成虎.水流分配策略随下坡坡度变化的多流向算法[J].水科学进展,2006,17(4):450-456

[2]于淼, 任立良. 基于DEM模型的新填洼算法[J]. 地球信息科学学报, 2009, 11(1): 50-55

[3]周倩倩,许小文,苏炯恒. GIS几何分析法在易涝区识别和空域构建中的应用[J]. 中国给水排水,2018,34(19):114-117,123.