简介崂山森林火灾扩散模拟分析与决策系统是一套面向森林防火应急响应与指挥决策的综合性软件工程适合GIS开发、应急管理及火灾建模方向的学习者与研究人员参考。系统融合地理信息系统、火灾扩散数学模型与决策支持技术涵盖数据采集预处理、火势蔓延预测、GIS空间分析、灭火方案评估、实时监控预警及后期效果评估等模块可帮助理解从火源定位到疏散路线规划的完整技术链路。资源包共1230个文件约215.19MB以177个C#源码文件与145个动态链接库构成核心程序122个gdbtable、122个gdbtablx及111个PNG、109个XML等地理数据库与配置资源支撑空间数据管理另有BMP、JPG、ICO等界面素材及SHP、MXD等GIS工程文件目录结构完整便于二次开发与模块拆解学习。目前已有221人学习下载适合作为森林火灾模拟与决策支持系统设计的实践参考。1. 崂山森林火灾扩散模拟分析与决策系统从火线到指挥屏的那条链路林火蔓延模拟这件事真正难的不是把火画出来而是让画出来的火在时间轴上站得住脚。崂山这类山地林区地形起伏大、植被类型交错、海陆风切换频繁一场火从山脊往下烧和从沟谷往上烧蔓延速度能差出好几倍。所谓「崂山森林火灾扩散模拟分析与决策系统」本质是把地形、植被、气象三类栅格数据喂给一个蔓延模型按时间步推演火线位置再把结果转成扑救力量该往哪投、隔离带该在哪开的决策建议。它适合三类人做林火研究的、做应急指挥系统开发的、以及需要把模拟结果接进大屏或预案演练的工程团队。这篇笔记不讲空概念只讲怎么把这条链路在本地跑通、参数怎么调、哪里最容易翻车。2. 蔓延模型选型Rothermel 为什么是山地场景的默认答案2.1 从三个候选模型说起林火蔓延模型不是只有一种。工程上常见的有三类元胞自动机类、Huygens 椭圆波扩展类、以及以 Rothermel 公式为核心的物理经验混合类。元胞自动机实现最简单一个栅格按邻居状态更新但它的蔓延方向是离散的八方向或十六方向山地坡度引起的上坡加速几乎没法自然表达结果往往是一个越来越圆的火团和真实火线形态差得远。Huygens 类把火线当成一圈不断向外扩展的椭圆波方向连续性好但椭圆长短轴比需要靠风速和坡度经验公式给参数标定工作量不小。Rothermel 模型是 1972 年提出的一套半经验公式核心思路是把火焰前锋单位面积的热量平衡拆成可燃物、风、坡度三部分贡献最后算出一个蔓延速度 R。它的优势在于坡度影响有显式项植被参数有相对成熟的查表体系而且被大量后续系统验证过。山地场景下坡度项直接决定上坡火加速、下坡火减速的行为这是元胞自动机给不了的。所以如果目标是崂山这种起伏地形Rothermel 是默认起点不是因为它最准而是因为它对地形和植被的响应最可解释、最容易调。2.2 Rothermel 的核心输入参数Rothermel 公式本身不复杂难的是输入参数怎么来。工程上一般把它拆成三组参数组代表量典型来源单位可燃物载量、表面积体积比、含水率、床层高度植被类型查表 实测kg/m²、1/m、%、m气象10m 风速、风向气象站或再分析数据m/s、度地形坡度、坡向DEM 计算度含水率是其中最敏感的一个。活可燃物含水率随季节变死可燃物含水率随日内湿度变同一个植被类型含水率从 8% 变到 12%蔓延速度可能掉一半。很多模拟结果离谱不是模型错是含水率随手填了个常数。2.3 用 Python 跑通最小蔓延计算下面这段代码是 Rothermel 蔓延速度的最小实现只算一个栅格在给定风、坡条件下的速度目的是让你先把公式跑通再谈栅格化。import math def rothermel_spread_rate(fuel_load, sav_ratio, moisture, bed_height, wind_speed, slope_deg): fuel_load: 可燃物载量 kg/m^2 sav_ratio: 表面积体积比 1/m moisture: 含水率小数如 0.08 bed_height: 床层高度 m wind_speed: 10m 风速 m/s slope_deg: 坡度 度 # 可燃物净载量扣除矿物质修正这里取简化系数 w0 fuel_load * 0.95 # 床层堆积比 rho_b w0 / bed_height # 最优堆积比经验式 rho_opt 0.0033 * (sav_ratio ** 1.15) # 相对堆积比 beta rho_b / 0.0012 beta_opt rho_opt / 0.0012 # 风因子简化形式 c 7.47 * math.exp(-0.133 * (sav_ratio ** 0.55)) b 0.02526 * (sav_ratio ** 0.54) e 0.715 * math.exp(-3.59e-4 * sav_ratio) phi_w c * (wind_speed ** b) * (beta / beta_opt) ** (-e) # 坡度因子 phi_s 5.275 * (beta ** -0.3) * (math.tan(math.radians(slope_deg)) ** 2) # 反应强度与传播通量这里用简化系数占位 gamma_max sav_ratio ** 1.5 / 495 beta_ratio beta / beta_opt a 133 * (sav_ratio ** -0.7913) gamma gamma_max * (beta_ratio ** a) * math.exp(a * (1 - beta_ratio)) # 含水率阻尼 eta_m 1 - 2.59 * (moisture / 0.3) 5.11 * (moisture / 0.3) ** 2 - 3.52 * (moisture / 0.3) ** 3 eta_m max(eta_m, 0.0) # 净反应强度 i_r gamma * w0 * 18000 * eta_m # 传播通量比 xi math.exp((0.792 0.681 * math.sqrt(sav_ratio)) * (beta 0.1)) / (192 0.2595 * sav_ratio) # 蔓延速度 m/min R (i_r * xi * (1 phi_w phi_s)) / (rho_b * 18000 * 0.5) return R # 示例中等载量、含水率 8%、风速 3m/s、上坡 20 度 rate rothermel_spread_rate(1.2, 2000, 0.08, 0.6, 3.0, 20) print(f蔓延速度: {rate:.2f} m/min)这段代码里phi_w是风因子phi_s是坡度因子两者相加再乘到基础传播项上这就是 Rothermel 对风和坡的响应方式。eta_m是含水率阻尼含水率越高净反应强度越低。参数说明sav_ratio对结果影响很大细燃料草、落叶能到 5000 以上粗燃料树枝只有几百bed_height是床层厚度不是树高填错会让堆积比整个偏掉。跑通之后你会发现坡度从 0 到 30 度速度可能翻三倍这就是山地场景必须用带坡度项模型的原因。3. 把模型装进栅格地形、植被、气象三套数据的对齐3.1 数据准备与坐标统一单点公式跑通只是第一步真实模拟要在栅格上逐像元算。崂山场景需要三套栅格DEM 高程、植被类型、以及随时间变化的气象场。三套数据必须对齐到同一套网格、同一分辨率、同一坐标系。常见做法是统一到 UTM 投影分辨率取 30m因为 DEM 公开数据大多是 30m植被分类图也常在这个尺度。对齐用 GDAL 或 rasterio 都行关键是重采样方式DEM 用双线性植被类型用最近邻气象场用双线性。植被类型如果用双线性会出现 0.5 类这种没有意义的中间值后面查表直接报错。import rasterio from rasterio.warp import reproject, Resampling def align_to_reference(src_path, ref_path, dst_path, method): with rasterio.open(ref_path) as ref: ref_crs ref.crs ref_transform ref.transform ref_shape (ref.height, ref.width) with rasterio.open(src_path) as src: data src.read(1) dst np.zeros(ref_shape, dtypedata.dtype) reproject( sourcedata, destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformref_transform, dst_crsref_crs, resamplingmethod ) profile src.profile.copy() profile.update(crsref_crs, transformref_transform, widthref_shape[1], heightref_shape[0]) with rasterio.open(dst_path, w, **profile) as out: out.write(dst, 1) # DEM 双线性植被最近邻 align_to_reference(dem.tif, ref.tif, dem_aligned.tif, Resampling.bilinear) align_to_reference(veg.tif, ref.tif, veg_aligned.tif, Resampling.nearest)reproject的resampling参数决定重采样方式Resampling.nearest对分类数据是必须的。ref.tif是基准网格一般选 DEM 或研究区边界裁剪后的模板。这一步做完三套数据每个像元一一对应后面循环才不会错位。3.2 逐像元蔓延与时间步推进栅格化蔓延有两种主流做法一种是每个时间步对全图算一遍速度然后按速度更新火线另一种是用最短路径或波前传播。工程上更稳的是前者配合一个火线状态栅格0 未燃、1 燃烧中、2 已燃尽。时间步不能太大。如果速度是 5 m/min分辨率 30m那一个像元烧完要 6 分钟时间步取 1 分钟比较稳。时间步太大火会「跳」过窄谷或窄脊形态失真。import numpy as np def step_fire(state, rate, dt, cell_size): state: 0未燃 1燃烧 2燃尽; rate: 蔓延速度 m/min new_state state.copy() burning np.argwhere(state 1) for i, j in burning: # 八邻域传播 for di in (-1, 0, 1): for dj in (-1, 0, 1): if di 0 and dj 0: continue ni, nj i di, j dj if 0 ni state.shape[0] and 0 nj state.shape[1]: if state[ni, nj] 0: # 距离按对角修正 dist cell_size * (1.414 if di and dj else 1.0) # 该方向速度取当前像元速度 t_burn dist / max(rate[i, j], 1e-6) if t_burn dt: new_state[ni, nj] 1 new_state[i, j] 2 return new_staterate是每个像元的速度栅格由上一章的公式逐像元算出来风向和坡度方向决定速度往哪个邻居传。这里简化成八邻域实际工程里会按风向做各向异性加权否则火会往各个方向等速跑。dt和cell_size的关系要盯住dt大于cell_size / rate时火会一步跨过多个像元形态就不可信了。3.3 气象场的时间插值气象数据一般是逐小时或逐三小时模拟时间步是分钟级中间要插值。风速风向不能直接线性插值风向角会出现 350 度和 10 度插出 180 度这种翻车。正确做法是把风向拆成 u、v 分量再插值最后反算角度。def interp_wind(wd1, ws1, wd2, ws2, t): wd 度ws m/st 0~1 u1 -ws1 * np.sin(np.radians(wd1)) v1 -ws1 * np.cos(np.radians(wd1)) u2 -ws2 * np.sin(np.radians(wd2)) v2 -ws2 * np.cos(np.radians(wd2)) u u1 (u2 - u1) * t v v1 (v2 - v1) * t ws np.hypot(u, v) wd (np.degrees(np.arctan2(-u, -v)) 360) % 360 return wd, ws风向角插值翻车是血泪经验直接对角度做线性插值在跨 0 度时结果完全错。拆分量是标准解法代价是多几行代码但省掉后面排查半天「为什么火突然往反方向烧」。4. 避坑与排查模拟结果不对劲时先看这五处4.1 火团越来越圆坡度像没起作用现象跑出来的火线是个近似圆形上坡下坡没区别。原因通常是坡度项没接进逐像元计算或者 DEM 对齐后坡度栅格全是 0。解决检查坡度栅格统计值确认不是常数再确认phi_s里用的坡度是当前像元朝火传播方向的坡度不是全局平均坡度。很多人把坡度算成标量场就完事忘了它是有方向的。4.2 火一步跳过整个山谷现象时间步设成 10 分钟火直接出现在对面山脊。原因dt大于cell_size / rate传播逻辑允许一步跨多格。解决把dt降到cell_size / max_rate以下或者改成子步循环每个子步只允许传播一个像元。代价是计算量上升但形态可信。4.3 含水率填常数导致季节间结果无差异现象同一块地春天和秋天模拟结果几乎一样。原因含水率写死成 0.08。解决按植被类型和季节建一张含水率查找表至少分活可燃物和死可燃物两列死可燃物再按日内湿度曲线调。这一步不做模型对季节的响应就是假的。4.4 风向插值跨 0 度后火往反方向烧现象气象数据从 350 度切到 10 度火突然掉头。原因直接对角度线性插值。解决拆 u、v 分量插值再反算见 3.3 的代码。这个坑几乎每个做气象驱动模拟的人都踩过。4.5 植被类型重采样出小数导致查表失败现象查植被参数时报 KeyError 或得到奇怪值。原因植被栅格用了双线性重采样出现 2.5 这种类型码。解决分类数据一律最近邻重采样重采样后检查唯一值集合是否还在原始类型码范围内。5. 从模拟结果到决策建议隔离带选址与力量投放的量化5.1 用到达时间场反推隔离带位置模拟跑完最有价值的不是某一帧火线图而是每个像元的到达时间场。到达时间场是一张和 DEM 同尺寸的栅格值是该像元被火到达的分钟数。有了它隔离带选址就变成一个约束优化问题在到达时间小于 T 的区域外找一条连续路径使得开挖代价最小、且能阻断火向目标保护区蔓延。常见做法是先按到达时间做等值线取 T60 分钟那条线作为候选带再叠加坡度太陡的机械上不去和道路可达性最后用最小成本路径算法在候选带里选一条。这一步不需要复杂优化库栅格上的 Dijkstra 就够。import heapq def min_cost_path(cost, start, end): cost: 2D 代价栅格; start/end: (row, col) h, w cost.shape dist np.full((h, w), np.inf) prev {} dist[start] 0 pq [(0, start)] while pq: d, (i, j) heapq.heappop(pq) if (i, j) end: break if d dist[i, j]: continue for di, dj in ((1,0),(-1,0),(0,1),(0,-1)): ni, nj idi, jdj if 0 ni h and 0 nj w: nd d cost[ni, nj] if nd dist[ni, nj]: dist[ni, nj] nd prev[(ni, nj)] (i, j) heapq.heappush(pq, (nd, (ni, nj))) # 回溯路径 path [] cur end while cur ! start: path.append(cur) cur prev[cur] path.append(start) return path[::-1]cost栅格由坡度、植被清除难度、道路距离加权得到。start和end是候选带两端。这条路径就是建议的隔离带走向。参数上坡度权重要给足超过 35 度的区域代价直接设成极大值因为机械上不去人工作业效率也极低。5.2 力量投放的优先级排序有了到达时间场力量投放就有了量化依据。把保护区、居民点、关键设施作为目标点每个目标点的「威胁时间」就是火到达它的最短时间。按威胁时间排序再结合道路通行时间就能排出先保谁、后保谁。这一步的输出是一张排序表不是一张图指挥屏上直接显示。目标点火到达时间(min)道路通行时间(min)优先级目标A4520高目标B9015中目标C15040低优先级不是简单按到达时间排而是按「到达时间减通行时间」的余量排余量越小越紧急。这个余量才是决策真正需要的数。5.3 结果验证用历史火场做回算模拟系统做完必须用历史火场回算验证。找一场有完整火线记录的火灾把当时的气象、植被、地形输入看模拟火线和实际火线的重合度。常用指标是交并比和面积误差。交并比低于 0.6 说明参数或模型有问题优先查含水率和风速。回算不是为了证明模型准而是为了标定参数——同一套参数在多个历史火场上都能到 0.7 以上这套参数才敢用于新火场预测。我自己的习惯是每换一个林区先拿两场历史火做标定标定完再跑预案。跳过这一步直接上指挥屏翻车是迟早的事。希望帮到你。本文还有配套的精品资源点击获取 崂山森林火灾扩散模拟分析与决策系统：Rothermel模型与栅格化实现
简介崂山森林火灾扩散模拟分析与决策系统是一套面向森林防火应急响应与指挥决策的综合性软件工程适合GIS开发、应急管理及火灾建模方向的学习者与研究人员参考。系统融合地理信息系统、火灾扩散数学模型与决策支持技术涵盖数据采集预处理、火势蔓延预测、GIS空间分析、灭火方案评估、实时监控预警及后期效果评估等模块可帮助理解从火源定位到疏散路线规划的完整技术链路。资源包共1230个文件约215.19MB以177个C#源码文件与145个动态链接库构成核心程序122个gdbtable、122个gdbtablx及111个PNG、109个XML等地理数据库与配置资源支撑空间数据管理另有BMP、JPG、ICO等界面素材及SHP、MXD等GIS工程文件目录结构完整便于二次开发与模块拆解学习。目前已有221人学习下载适合作为森林火灾模拟与决策支持系统设计的实践参考。1. 崂山森林火灾扩散模拟分析与决策系统从火线到指挥屏的那条链路林火蔓延模拟这件事真正难的不是把火画出来而是让画出来的火在时间轴上站得住脚。崂山这类山地林区地形起伏大、植被类型交错、海陆风切换频繁一场火从山脊往下烧和从沟谷往上烧蔓延速度能差出好几倍。所谓「崂山森林火灾扩散模拟分析与决策系统」本质是把地形、植被、气象三类栅格数据喂给一个蔓延模型按时间步推演火线位置再把结果转成扑救力量该往哪投、隔离带该在哪开的决策建议。它适合三类人做林火研究的、做应急指挥系统开发的、以及需要把模拟结果接进大屏或预案演练的工程团队。这篇笔记不讲空概念只讲怎么把这条链路在本地跑通、参数怎么调、哪里最容易翻车。2. 蔓延模型选型Rothermel 为什么是山地场景的默认答案2.1 从三个候选模型说起林火蔓延模型不是只有一种。工程上常见的有三类元胞自动机类、Huygens 椭圆波扩展类、以及以 Rothermel 公式为核心的物理经验混合类。元胞自动机实现最简单一个栅格按邻居状态更新但它的蔓延方向是离散的八方向或十六方向山地坡度引起的上坡加速几乎没法自然表达结果往往是一个越来越圆的火团和真实火线形态差得远。Huygens 类把火线当成一圈不断向外扩展的椭圆波方向连续性好但椭圆长短轴比需要靠风速和坡度经验公式给参数标定工作量不小。Rothermel 模型是 1972 年提出的一套半经验公式核心思路是把火焰前锋单位面积的热量平衡拆成可燃物、风、坡度三部分贡献最后算出一个蔓延速度 R。它的优势在于坡度影响有显式项植被参数有相对成熟的查表体系而且被大量后续系统验证过。山地场景下坡度项直接决定上坡火加速、下坡火减速的行为这是元胞自动机给不了的。所以如果目标是崂山这种起伏地形Rothermel 是默认起点不是因为它最准而是因为它对地形和植被的响应最可解释、最容易调。2.2 Rothermel 的核心输入参数Rothermel 公式本身不复杂难的是输入参数怎么来。工程上一般把它拆成三组参数组代表量典型来源单位可燃物载量、表面积体积比、含水率、床层高度植被类型查表 实测kg/m²、1/m、%、m气象10m 风速、风向气象站或再分析数据m/s、度地形坡度、坡向DEM 计算度含水率是其中最敏感的一个。活可燃物含水率随季节变死可燃物含水率随日内湿度变同一个植被类型含水率从 8% 变到 12%蔓延速度可能掉一半。很多模拟结果离谱不是模型错是含水率随手填了个常数。2.3 用 Python 跑通最小蔓延计算下面这段代码是 Rothermel 蔓延速度的最小实现只算一个栅格在给定风、坡条件下的速度目的是让你先把公式跑通再谈栅格化。import math def rothermel_spread_rate(fuel_load, sav_ratio, moisture, bed_height, wind_speed, slope_deg): fuel_load: 可燃物载量 kg/m^2 sav_ratio: 表面积体积比 1/m moisture: 含水率小数如 0.08 bed_height: 床层高度 m wind_speed: 10m 风速 m/s slope_deg: 坡度 度 # 可燃物净载量扣除矿物质修正这里取简化系数 w0 fuel_load * 0.95 # 床层堆积比 rho_b w0 / bed_height # 最优堆积比经验式 rho_opt 0.0033 * (sav_ratio ** 1.15) # 相对堆积比 beta rho_b / 0.0012 beta_opt rho_opt / 0.0012 # 风因子简化形式 c 7.47 * math.exp(-0.133 * (sav_ratio ** 0.55)) b 0.02526 * (sav_ratio ** 0.54) e 0.715 * math.exp(-3.59e-4 * sav_ratio) phi_w c * (wind_speed ** b) * (beta / beta_opt) ** (-e) # 坡度因子 phi_s 5.275 * (beta ** -0.3) * (math.tan(math.radians(slope_deg)) ** 2) # 反应强度与传播通量这里用简化系数占位 gamma_max sav_ratio ** 1.5 / 495 beta_ratio beta / beta_opt a 133 * (sav_ratio ** -0.7913) gamma gamma_max * (beta_ratio ** a) * math.exp(a * (1 - beta_ratio)) # 含水率阻尼 eta_m 1 - 2.59 * (moisture / 0.3) 5.11 * (moisture / 0.3) ** 2 - 3.52 * (moisture / 0.3) ** 3 eta_m max(eta_m, 0.0) # 净反应强度 i_r gamma * w0 * 18000 * eta_m # 传播通量比 xi math.exp((0.792 0.681 * math.sqrt(sav_ratio)) * (beta 0.1)) / (192 0.2595 * sav_ratio) # 蔓延速度 m/min R (i_r * xi * (1 phi_w phi_s)) / (rho_b * 18000 * 0.5) return R # 示例中等载量、含水率 8%、风速 3m/s、上坡 20 度 rate rothermel_spread_rate(1.2, 2000, 0.08, 0.6, 3.0, 20) print(f蔓延速度: {rate:.2f} m/min)这段代码里phi_w是风因子phi_s是坡度因子两者相加再乘到基础传播项上这就是 Rothermel 对风和坡的响应方式。eta_m是含水率阻尼含水率越高净反应强度越低。参数说明sav_ratio对结果影响很大细燃料草、落叶能到 5000 以上粗燃料树枝只有几百bed_height是床层厚度不是树高填错会让堆积比整个偏掉。跑通之后你会发现坡度从 0 到 30 度速度可能翻三倍这就是山地场景必须用带坡度项模型的原因。3. 把模型装进栅格地形、植被、气象三套数据的对齐3.1 数据准备与坐标统一单点公式跑通只是第一步真实模拟要在栅格上逐像元算。崂山场景需要三套栅格DEM 高程、植被类型、以及随时间变化的气象场。三套数据必须对齐到同一套网格、同一分辨率、同一坐标系。常见做法是统一到 UTM 投影分辨率取 30m因为 DEM 公开数据大多是 30m植被分类图也常在这个尺度。对齐用 GDAL 或 rasterio 都行关键是重采样方式DEM 用双线性植被类型用最近邻气象场用双线性。植被类型如果用双线性会出现 0.5 类这种没有意义的中间值后面查表直接报错。import rasterio from rasterio.warp import reproject, Resampling def align_to_reference(src_path, ref_path, dst_path, method): with rasterio.open(ref_path) as ref: ref_crs ref.crs ref_transform ref.transform ref_shape (ref.height, ref.width) with rasterio.open(src_path) as src: data src.read(1) dst np.zeros(ref_shape, dtypedata.dtype) reproject( sourcedata, destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformref_transform, dst_crsref_crs, resamplingmethod ) profile src.profile.copy() profile.update(crsref_crs, transformref_transform, widthref_shape[1], heightref_shape[0]) with rasterio.open(dst_path, w, **profile) as out: out.write(dst, 1) # DEM 双线性植被最近邻 align_to_reference(dem.tif, ref.tif, dem_aligned.tif, Resampling.bilinear) align_to_reference(veg.tif, ref.tif, veg_aligned.tif, Resampling.nearest)reproject的resampling参数决定重采样方式Resampling.nearest对分类数据是必须的。ref.tif是基准网格一般选 DEM 或研究区边界裁剪后的模板。这一步做完三套数据每个像元一一对应后面循环才不会错位。3.2 逐像元蔓延与时间步推进栅格化蔓延有两种主流做法一种是每个时间步对全图算一遍速度然后按速度更新火线另一种是用最短路径或波前传播。工程上更稳的是前者配合一个火线状态栅格0 未燃、1 燃烧中、2 已燃尽。时间步不能太大。如果速度是 5 m/min分辨率 30m那一个像元烧完要 6 分钟时间步取 1 分钟比较稳。时间步太大火会「跳」过窄谷或窄脊形态失真。import numpy as np def step_fire(state, rate, dt, cell_size): state: 0未燃 1燃烧 2燃尽; rate: 蔓延速度 m/min new_state state.copy() burning np.argwhere(state 1) for i, j in burning: # 八邻域传播 for di in (-1, 0, 1): for dj in (-1, 0, 1): if di 0 and dj 0: continue ni, nj i di, j dj if 0 ni state.shape[0] and 0 nj state.shape[1]: if state[ni, nj] 0: # 距离按对角修正 dist cell_size * (1.414 if di and dj else 1.0) # 该方向速度取当前像元速度 t_burn dist / max(rate[i, j], 1e-6) if t_burn dt: new_state[ni, nj] 1 new_state[i, j] 2 return new_staterate是每个像元的速度栅格由上一章的公式逐像元算出来风向和坡度方向决定速度往哪个邻居传。这里简化成八邻域实际工程里会按风向做各向异性加权否则火会往各个方向等速跑。dt和cell_size的关系要盯住dt大于cell_size / rate时火会一步跨过多个像元形态就不可信了。3.3 气象场的时间插值气象数据一般是逐小时或逐三小时模拟时间步是分钟级中间要插值。风速风向不能直接线性插值风向角会出现 350 度和 10 度插出 180 度这种翻车。正确做法是把风向拆成 u、v 分量再插值最后反算角度。def interp_wind(wd1, ws1, wd2, ws2, t): wd 度ws m/st 0~1 u1 -ws1 * np.sin(np.radians(wd1)) v1 -ws1 * np.cos(np.radians(wd1)) u2 -ws2 * np.sin(np.radians(wd2)) v2 -ws2 * np.cos(np.radians(wd2)) u u1 (u2 - u1) * t v v1 (v2 - v1) * t ws np.hypot(u, v) wd (np.degrees(np.arctan2(-u, -v)) 360) % 360 return wd, ws风向角插值翻车是血泪经验直接对角度做线性插值在跨 0 度时结果完全错。拆分量是标准解法代价是多几行代码但省掉后面排查半天「为什么火突然往反方向烧」。4. 避坑与排查模拟结果不对劲时先看这五处4.1 火团越来越圆坡度像没起作用现象跑出来的火线是个近似圆形上坡下坡没区别。原因通常是坡度项没接进逐像元计算或者 DEM 对齐后坡度栅格全是 0。解决检查坡度栅格统计值确认不是常数再确认phi_s里用的坡度是当前像元朝火传播方向的坡度不是全局平均坡度。很多人把坡度算成标量场就完事忘了它是有方向的。4.2 火一步跳过整个山谷现象时间步设成 10 分钟火直接出现在对面山脊。原因dt大于cell_size / rate传播逻辑允许一步跨多格。解决把dt降到cell_size / max_rate以下或者改成子步循环每个子步只允许传播一个像元。代价是计算量上升但形态可信。4.3 含水率填常数导致季节间结果无差异现象同一块地春天和秋天模拟结果几乎一样。原因含水率写死成 0.08。解决按植被类型和季节建一张含水率查找表至少分活可燃物和死可燃物两列死可燃物再按日内湿度曲线调。这一步不做模型对季节的响应就是假的。4.4 风向插值跨 0 度后火往反方向烧现象气象数据从 350 度切到 10 度火突然掉头。原因直接对角度线性插值。解决拆 u、v 分量插值再反算见 3.3 的代码。这个坑几乎每个做气象驱动模拟的人都踩过。4.5 植被类型重采样出小数导致查表失败现象查植被参数时报 KeyError 或得到奇怪值。原因植被栅格用了双线性重采样出现 2.5 这种类型码。解决分类数据一律最近邻重采样重采样后检查唯一值集合是否还在原始类型码范围内。5. 从模拟结果到决策建议隔离带选址与力量投放的量化5.1 用到达时间场反推隔离带位置模拟跑完最有价值的不是某一帧火线图而是每个像元的到达时间场。到达时间场是一张和 DEM 同尺寸的栅格值是该像元被火到达的分钟数。有了它隔离带选址就变成一个约束优化问题在到达时间小于 T 的区域外找一条连续路径使得开挖代价最小、且能阻断火向目标保护区蔓延。常见做法是先按到达时间做等值线取 T60 分钟那条线作为候选带再叠加坡度太陡的机械上不去和道路可达性最后用最小成本路径算法在候选带里选一条。这一步不需要复杂优化库栅格上的 Dijkstra 就够。import heapq def min_cost_path(cost, start, end): cost: 2D 代价栅格; start/end: (row, col) h, w cost.shape dist np.full((h, w), np.inf) prev {} dist[start] 0 pq [(0, start)] while pq: d, (i, j) heapq.heappop(pq) if (i, j) end: break if d dist[i, j]: continue for di, dj in ((1,0),(-1,0),(0,1),(0,-1)): ni, nj idi, jdj if 0 ni h and 0 nj w: nd d cost[ni, nj] if nd dist[ni, nj]: dist[ni, nj] nd prev[(ni, nj)] (i, j) heapq.heappush(pq, (nd, (ni, nj))) # 回溯路径 path [] cur end while cur ! start: path.append(cur) cur prev[cur] path.append(start) return path[::-1]cost栅格由坡度、植被清除难度、道路距离加权得到。start和end是候选带两端。这条路径就是建议的隔离带走向。参数上坡度权重要给足超过 35 度的区域代价直接设成极大值因为机械上不去人工作业效率也极低。5.2 力量投放的优先级排序有了到达时间场力量投放就有了量化依据。把保护区、居民点、关键设施作为目标点每个目标点的「威胁时间」就是火到达它的最短时间。按威胁时间排序再结合道路通行时间就能排出先保谁、后保谁。这一步的输出是一张排序表不是一张图指挥屏上直接显示。目标点火到达时间(min)道路通行时间(min)优先级目标A4520高目标B9015中目标C15040低优先级不是简单按到达时间排而是按「到达时间减通行时间」的余量排余量越小越紧急。这个余量才是决策真正需要的数。5.3 结果验证用历史火场做回算模拟系统做完必须用历史火场回算验证。找一场有完整火线记录的火灾把当时的气象、植被、地形输入看模拟火线和实际火线的重合度。常用指标是交并比和面积误差。交并比低于 0.6 说明参数或模型有问题优先查含水率和风速。回算不是为了证明模型准而是为了标定参数——同一套参数在多个历史火场上都能到 0.7 以上这套参数才敢用于新火场预测。我自己的习惯是每换一个林区先拿两场历史火做标定标定完再跑预案。跳过这一步直接上指挥屏翻车是迟早的事。希望帮到你。本文还有配套的精品资源点击获取