☰
崂山森林火灾扩散模拟分析与决策系统:Rothermel模型与栅格化实现
2026/10/9 16:53:51 网站建设 项目流程

简介:崂山森林火灾扩散模拟分析与决策系统是一套面向森林防火应急响应与指挥决策的综合性软件工程,适合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, dtype=data.dtype) reproject( source=data, destination=dst, src_transform=src.transform, src_crs=src.crs, dst_transform=ref_transform, dst_crs=ref_crs, resampling=method ) profile = src.profile.copy() profile.update(crs=ref_crs, transform=ref_transform, width=ref_shape[1], height=ref_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_state

rate是每个像元的速度栅格,由上一章的公式逐像元算出来,风向和坡度方向决定速度往哪个邻居传。这里简化成八邻域,实际工程里会按风向做各向异性加权,否则火会往各个方向等速跑。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/s,t 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 的区域外,找一条连续路径,使得开挖代价最小、且能阻断火向目标保护区蔓延。

常见做法是:先按到达时间做等值线,取 T=60 分钟那条线作为候选带,再叠加坡度(太陡的机械上不去)和道路可达性,最后用最小成本路径算法在候选带里选一条。这一步不需要复杂优化库,栅格上的 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 = i+di, j+dj 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 以上,这套参数才敢用于新火场预测。

我自己的习惯是:每换一个林区,先拿两场历史火做标定,标定完再跑预案。跳过这一步直接上指挥屏,翻车是迟早的事。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询