简介面向自然灾害研究与应急管理领域这份资源包以2020年澳大利亚山火为背景提供了基于MATLAB的AusFire山火蔓延模型。包内共5个文件包含3个M源码脚本、1份PDF技术说明及配套演示文稿M脚本分别负责火线扩散主算法、气象与地形参数解析、模拟结果可视化PDF文档详细介绍了模型假设、差分方程求解及参数校核方法PPT则汇总了不同场景下的火势演进案例。压缩包整体体积仅12.92MB目录结构紧凑方便科研人员、高校师生及防灾部门快速借鉴。目前已有406人学习使用显示其在实际研究中的适用性。借助这份资源可完整复现山火蔓延动态过程理解风速、湿度、植被与地形如何影响火线走向为制定疏散方案和优化灭火资源调度提供量化分析工具具有较高的学术与实践价值。1. 澳大利亚山火模型AusFire用Matlab把林火蔓延从公式推到PPT先下一个判断澳大利亚山火模型AusFire并没有全行业统一的官方案例更常见的是把它当成一套“McArthur燃烧公式 火线几何推进 结果可视化”的组合实现。决定模型质量的关键不在界面而在对澳大利亚干燥桉树林环境参数适配。我按类似林火模拟项目里反复走过的流程从FFDI火险指数、燃速函数、Huygens子波火线算法讲到Matlab画图与PPT汇报让气象、林业、应急和GIS方向的人在一台装了Matlab的机器上把一场山火从气象输入推到可以放进汇报里的结论。新手能照着命令跑通老手可以直接拿走参数边界和坑位清单。2. 燃速计算的物理量把FFDI与McArthur公式变成Matlab可复用函数2.1 为什么先算FFDIAusFire类模型绕不开火险指数林火蔓延的核心输入不是温度或风速单个变量而是它们的组合结果。澳大利亚气象部门每天发布的FFDIForest Fire Danger Index森林火险指数背后是McArthur Mark 5公式温度、相对湿度、风速、干旱因子四个量共同决定可燃物燃烧的剧烈程度。在AusFire类模型的参数文件里常见做法是把一天之内的气象站观测插值到模拟网格逐栅格算FFDI再套一层燃速公式得到每个位置的理论蔓延速度。气象站的CSV数据直接用readmatrix读进来下面这个函数是FFDI的最小实现。注意Matlab的log是自然对数不需要额外数学库function FFDI calcFFDI(T, RH, V, DF) % calcFFDI: McArthur Mark 5 火险指数 % T 气温(°C) | RH 相对湿度(%) | V 10m平均风速(km/h) | DF 干旱因子(0~10) FFDI 2 * exp(-0.45 0.987 * log(DF) ... - 0.0345 * RH 0.0338 * T 0.0234 * V); end这里有个数据单位陷阱相对湿度用百分比风速用km/h而不是m/s。拿到m/s的风速先乘3.6再传入否则FFDI会被系统性压低。DF干旱因子通常由KBDI干旱指数换算少数站点没有KBDI时用当地消防部门公开的每日FFDI产品反推DF这比从土壤湿度自己凑一个数稳妥。四个参数里DF的跨度最大对结果影响也最大它是模型的“底噪”。参数符号单位典型范围量级感气温T°C10~45每升10°CFFDI约放大1.4倍相对湿度RH%5~60每降10%约放大1.4倍风速Vkm/h0~60每增10km/h约放大1.26倍干旱因子DF无量纲0~10十倍跨度决定整体量级用一组夏季极端条件做个单点验证FFDI calcFFDI(35, 25, 30, 9); % 35°C / 25% / 30km/h / DF9这个算例对应干热大风天FFDI落在70~90区间已经进入澳大利亚火险等级里的最高档。拿你所在时区当天官方发布的FFDI对拍偏差在5以内说明气象数据单位处理没有问题偏差到了两位数先回头查风速和湿度单位再查DF是否来自同一套干旱指数体系。提示公式中的风速取10米高度平均风不是阵风。静态推演阶段用阵风直接代公式会把FFDI系统性抬高阵风只能放在逐分钟重算的动态场景里。2.2 用Matlab写燃速函数McArthur公式与坡度修正FFDI本身不等于蔓延速度还要过一道燃速公式。McArthur森林模型在文献里常见的分段形式是FFDI小于5走轻量级大于等于5走线性高段两段在5处连续。我习惯把坡度修正和可燃物载荷修正一起封装进同一个函数function ROS mcarthurROS(FFDI, slopeDeg, fuelLoad) % mcarthurROS: 林火蔓延速度, 单位 km/h % FFDI 火险指数 | slopeDeg 上坡坡度(°) | fuelLoad 可燃物载荷(t/ha, 可选) if FFDI 5 ROS 0.0012 * FFDI; else ROS 0.0027 * FFDI - 0.0075; end % 分段连续, FFDI5时两段都等于0.006 if slopeDeg 0 ROS ROS * exp(0.069 * slopeDeg); end if nargin 2 fuelLoad 0 ROS ROS * sqrt(fuelLoad / 2.0); % 以2 t/ha为基准 end ROS max(ROS, 1e-3); % 数值保护, 防止极低火险除零 end两个细节值得展开。坡度因子用指数形式20°坡时燃速放大exp(0.069×20)≈4倍这是山地火蔓延速度的主要推力可燃物载荷修正用平方根意思是载荷翻倍燃速只涨约41%不要高估载荷的作用。跑气候情景时fuelLoad从燃料图直接读DF和fuelLoad是两个独立变量前者表示燃料干湿程度后者表示单位面积可燃物总量混用是最常见的建模错误。2.3 风与坡不一致时有效风向量合成与反演标定坡度和风的方向经常不一致不能把风速和坡度分别乘进去完事。我一般把风向量和坡度向量做一次矢量合成坡度向量指向上坡方向、模长由坡度值决定风向量由风向风速决定合成后的方向作为椭圆长轴方向合成后的模长再映射回燃速。这里给出合成代码的骨架% 注意: 下面角度约定为x轴正方向0°, 逆时针为正 wx V * cosd(windDirCart); wy V * sind(windDirCart); sx 20 * cosd(slopeDir); sy 20 * sind(slopeDir); effX wx sx; effY wy sy; effDirCart atan2d(effY, effX); % 有效风向(笛卡尔角) effSpeed hypot(effX, effY); % 有效风速坡度向量模长取20理由是McArthur体系里20°坡爬升的燃速放大和20km/h风速贡献量级相当具体权重应该用所在区域的历史火行为反演确定。手里有几次已知火场蔓延记录时用优化工具箱的fminsearch以火头位置误差为目标函数反求这个权重比手调参数靠谱。气象站给的风向通常是“正北0°顺时针”的来向角要转成上面代码的笛卡尔去向角换算关系是取270−dir再mod 360或者读取时直接做90−dir两种写法等价关键在于全项目统一。3. 火线怎么长出来Huygens子波法与元胞自动机的Matlab实现3.1 两条路线的选型先跑通CA再升级到Huygens火线推进算法主流是两条路线Huygens子波法和元胞自动机CA。选型结论先给只想要燃烧面积曲线或教学演示用CA做应急业务、火头位置要精确到几百米用Huygens。AusFire这个方向公开的演示实现大多是CA但业务系统最终都会往子波法靠。两个维度的差别可以看表维度Huygens子波法元胞自动机火线形态连续光滑折线网格锯齿边界方向分辨率子波长轴连续可设受8邻域限制时间步长5~30s20~120s实现成本中上低典型定位业务推演快速预案、教学我一般建议先做CA拿到面积曲线和大致火场形态再写Huygens验证火头位置。这个顺序对调试也有好处两类算法的报错特征完全不同第5章的边界泄漏和时间步长问题分不清你在跑哪种算法时会浪费很多时间。3.2 Huygens子波法火线点集外推与凸包重采样Huygens子波法的核心假设是火线上每个点独立向外辐射一个椭圆子波椭圆长轴与有效风向一致长轴长度等于ROS×dt短轴按燃料类型取长轴的0.3~0.7倍所有子波的外包络就是下一时刻的火线。Matlab里外包络恰好有现成函数convhullfunction newFront huygensAdvance(front, effDir, ROSmax, dt) % front: Nx2火线点集(米) | effDir: 风去向(气象角, 北0°顺时) % ROSmax: 火头燃速(km/h) | dt: 时间步长(s) a ROSmax * 1000 / 3600 * dt; % 长半轴(m), 单位换算 b max(0.5 * a, 10); % 短半轴, 按燃料类型调整 theta deg2rad(90 - effDir); % 气象去向角转笛卡尔角 phi linspace(0, 2*pi, 36); % 每个子波36个采样点 R [cos(theta) -sin(theta); sin(theta) cos(theta)]; allPts zeros(0, 2); for i 1:size(front, 1) xy [a*cos(phi), b*sin(phi)] * R; allPts [allPts; xy front(i, :)]; end k convhull(allPts(:, 1), allPts(:, 2)); newFront allPts(k, :); newFront resampleFront(newFront, 200); end代码逻辑分四步换算长半轴并构造旋转矩阵每个火线点生成36个椭圆采样点并平移到该点全部采样点做凸包得到下一条火线最后重采样恢复均匀点距。重采样这步不能省凸包顶点在平直段会很稀、在尖角处很密点距不均匀会让下一步的子波发生重叠或空洞火线慢慢变形。重采样我按周长坐标插值function frontOut resampleFront(front, nPts) % 先把火线顶点闭合成环 front [front; front(1, :)]; d [0; cumsum(hypot(diff(front(:,1)), diff(front(:,2))))]; s linspace(0, d(end), nPts); frontOut [interp1(d, front(:,1), s), interp1(d, front(:,2), s)]; frontOut(end, :) []; % 去掉与起点重复的闭合点 end风向约定是整个函数最容易错的地方。气象站给的是“来向”即风从哪个方向吹来这里effDir要求的是“去向”指向火头移动方向两者差180°。算错的话整条火线会朝上风向长而且形态上看不出异常只能通过和卫星过火区对比才能发现。3.3 元胞自动机版本网格点火与方向权重CA实现更接近直觉每个栅格三态未燃、燃烧、燃尽。燃烧单元以一定概率点燃8邻域概率由局部燃速、时间步长、格网尺寸和风向夹角共同决定。下面是一个可直接循环调用的单步函数function [gridNew, areaBurned] caStep(grid, ROSmap, cellSize, dt, windTo) % grid: 1未燃 2燃烧 3燃尽 | ROSmap: 各格最大燃速(m/s) % cellSize: 格网边长(m) | dt: 时间步长(s) % windTo: 风去向(笛卡尔角, x轴正方向0°, 逆时针为正) [ny, nx] size(grid); gridNew grid; [rB, cB] find(grid 2); for k 1:numel(rB) r rB(k); c cB(k); for dr -1:1 for dc -1:1 nr rdr; nc cdc; if nr1 || nc1 || nrny || ncnx, continue; end if grid(nr, nc) ~ 1, continue; end ang atan2d(dr, dc); % 邻居方向角 cosw max(0, cosd(ang - windTo)); % 方向权重 p min(1, ROSmap(nr,nc) * dt / cellSize ... * (0.3 0.7 * cosw)); if rand p gridNew(nr, nc) 2; end end end end gridNew(grid 2) 3; % 燃烧格燃尽 areaBurned nnz(gridNew 2) * cellSize^2; % 过火面积(m^2) end点火概率里括号内的0.3是逆风下限用于近似飞火和辐射热的零星引燃这个系数在不同燃料图上要重新标定。方向权重只用了cos的一次方严格来说应该按椭圆方程投影但对30~60m的格网CA本身的离散误差远大于这个方向精度差别。如果风的数据是气象来向进函数前先加180°再换算成笛卡尔角换算规则和3.2节一致。3.4 两种算法共同的时间步长约束不管哪条路线都要守一条类似CFL条件的约束单步推进距离不能超过一个格网尺度量级。CA版本直接体现在概率上要求ROS×dt/cellSize远小于1取0.3以下安全Huygens版本要求a大于子波采样点间隔、同时小于火线局部曲率半径否则凸包会出现虚假尖角火线看起来“带刺”。按cellSize30m、ROSmax0.5m/s估算dt取20~30s跑6小时约720~1080步。每步对500×500栅格全扫描纯循环在Matlab里大约几秒总时间在可接受范围做参数敏感性分析时建议改自适应步长全局最大燃速低于阈值时放宽到60s火头接近高燃速区自动收紧能省一半以上的计算时间代价只是多写几行判定逻辑。4. 把AusFire结果变成汇报素材Matlab画图与PPT导出链路4.1 三张必出图热力图、火线序列、燃烧面积曲线模拟跑完先出三张图第一张用imagesc画燃烧状态热力图第二张把每15分钟的火线叠成序列图第三张画燃烧面积随时间曲线。这三张是澳大利亚山火PPT里最常用的标准配置一张交代空间、一张交代过程、一张交代量级。% 燃烧状态热力图 figure(Position, [100 100 800 500]); imagesc(xVec, yVec, burnGrid); axis xy; colormap(flipud(hot)); colorbar; xlabel(x / m); ylabel(y / m); title(AusFire 模拟末期燃烧状态); % 火线时间序列(每隔15分钟一条) figure; for i 1:numel(Fronts) plot(Fronts{i}(:,1)/1000, Fronts{i}(:,2)/1000, LineWidth, 1.2); hold on; end legend(cellstr(string(times, HH:mm))); % times为datetime数组 % 燃烧面积曲线 figure; plot(minutesElapsed/60, areaKm2, b-o, MarkerSize, 4); xlabel(时间 / h); ylabel(燃烧面积 / km^2); grid on;第二张图的图例用了datetime转string的写法string(times,HH:mm)在R2021b之后都支持。三张图的保存统一用exportgraphics分辨率给300这套函数对中文标题和坐标轴标签的字形支持比旧版print稳定。输出文件名里带上模拟时刻后面对PPT和复盘都方便fname AusFire_frame_ string(simClock, yyyyMMdd_HHmm) .png; exportgraphics(fig, fname, Resolution, 300);4.2 用Matlab直接生成PPT的三个办法把结果图拼成PPT常见的路线有三条mlreportgen.pp官方API、Windows下的COM自动化、先导出PNG再用python-pptx汇总。我一般优先mlreportgen.pp因为跨平台Linux和macOS也能批处理适合每天自动出一版报告。三者的取舍列成表办法跨平台批量能力主要坑mlreportgen.pp是强占位符名称依赖模板COM自动化仅Windows弱慢、依赖Office安装先PNG再python-pptx是强多一道文件交接mlreportgen.pp最小可用代码import mlreportgen.pp.* ppt Presentation(AusFire_Result.pptx); open(ppt); s add(ppt, Title and Content); replace(s, Title, AusFire模拟2020-01-06 14:00); img Picture(AusFire_frame_20200106_1400.png); img.Width 6in; img.Height 4in; replace(s, Content, img); close(ppt);两个高频坑第一占位符名称依赖所选模板默认模板图片位叫Content换成自定义模板后要先读模板的Placeholder列表再replace第二Picture的宽高最好不要同时写死先算好图片宽高比只设一个值否则图片会被拉伸变形。COM自动化在Windows上可以actxserver(PowerPoint.Application)交互性强但慢且依赖Office安装生产线环境很少用。4.3 澳大利亚山火PPT的页面结构建议面对非建模背景的汇报对象页面顺序建议按“结果倒推参数”安排先放最终燃烧范围与卫星影像叠加图再放燃烧面积曲线最后才放参数表。模拟动画放中间不要放第一页。参数表只列对结果影响最大的四项FFDI、风速风向、燃料载荷、坡度因子。燃料类型标签缺失时需要用遥感波段自己分类Matlab自带的kmeans函数可以直接用。把几个波段的反射率排成N×B矩阵kmeans聚类出3~5类分类结果转成燃料类型图再喂给模型比手工勾画边界快一个量级。聚类数量用silhouette轮廓图选不要拍脑袋定K值。5. 调试AusFire的三个高频问题边界泄漏、时间步长与燃料图分辨率5.1 边界泄漏与虚假尖角最典型的现象是火线接近模拟区边缘后趋向平直或者Huygens凸包在角点甩出很长的尖刺。前者是CA在边界外直接continue表现为“火到边界就停”后者是子波在角点的包络被人为截断。我一般做法是在网格四周补10圈不可燃单元模拟区外不做任何特殊处理出图时裁掉补边。这样CA不会在边界外静默停止Huygens的子波也不会在角点生成异常长轴。5.2 时间步长与库朗数确认逻辑没错之后火头推进距离明显超过ROS×dt或者火线出现锯齿优先缩短时间步长。经验值是库朗数不超过0.3dtSafe 0.3 * cellSize / max(ROS(:)); dt min(dtUser, dtSafe);ROS(:)里的最大值来自之前插值好的燃速场注意单位必须是m/s。自适应时间步长实现也不复杂每跑完一步重新求max(ROS)平稳段放宽到0.5高燃速段自动收紧。5.3 燃料图分辨率与投影对齐第三种高频问题不是算法而是数据。模拟火场与卫星过火区整体错位、且错位方向固定八成是投影坐标系不一致。GDA94/MGA分区和WGS84经纬度混用是这个方向的日常。Matlab里用readgeoraster读GeoTIFF检查RasterReference的坐标系再用projcrs按EPSG码创建投影对象做转换GDA2020/MGA zone 55对应EPSG 7855。不要手工加减经纬度偏移不同椭球体的差值不是常数。燃料图分辨率比模拟网格粗时插值前先想清楚变量性质。燃料载荷是连续变量griddata双线性插值没问题燃料类型是类别变量直接插值会插出不存在的中间类型。我通常把类别转成one-hot后再逐类插值最后取argmax决定类型。做完这一套火线形状和卫星观测的差异才真正落在模型参数上而不是数据管道上。最后补一个验证技巧不要只用燃烧面积曲线评判模型。面积可能因为网格粗而“看起来接近”火头位置和火线形状可能差很远。把模拟火线与卫星过火矢量叠加后计算双向Hausdorff距离这个指标对火头偏差远比面积误差敏感也是模型交付时最常被追问的一项。本文还有配套的精品资源点击获取 用Matlab实现AusFire山火模型：从FFDI计算到Huygens火线模拟
简介面向自然灾害研究与应急管理领域这份资源包以2020年澳大利亚山火为背景提供了基于MATLAB的AusFire山火蔓延模型。包内共5个文件包含3个M源码脚本、1份PDF技术说明及配套演示文稿M脚本分别负责火线扩散主算法、气象与地形参数解析、模拟结果可视化PDF文档详细介绍了模型假设、差分方程求解及参数校核方法PPT则汇总了不同场景下的火势演进案例。压缩包整体体积仅12.92MB目录结构紧凑方便科研人员、高校师生及防灾部门快速借鉴。目前已有406人学习使用显示其在实际研究中的适用性。借助这份资源可完整复现山火蔓延动态过程理解风速、湿度、植被与地形如何影响火线走向为制定疏散方案和优化灭火资源调度提供量化分析工具具有较高的学术与实践价值。1. 澳大利亚山火模型AusFire用Matlab把林火蔓延从公式推到PPT先下一个判断澳大利亚山火模型AusFire并没有全行业统一的官方案例更常见的是把它当成一套“McArthur燃烧公式 火线几何推进 结果可视化”的组合实现。决定模型质量的关键不在界面而在对澳大利亚干燥桉树林环境参数适配。我按类似林火模拟项目里反复走过的流程从FFDI火险指数、燃速函数、Huygens子波火线算法讲到Matlab画图与PPT汇报让气象、林业、应急和GIS方向的人在一台装了Matlab的机器上把一场山火从气象输入推到可以放进汇报里的结论。新手能照着命令跑通老手可以直接拿走参数边界和坑位清单。2. 燃速计算的物理量把FFDI与McArthur公式变成Matlab可复用函数2.1 为什么先算FFDIAusFire类模型绕不开火险指数林火蔓延的核心输入不是温度或风速单个变量而是它们的组合结果。澳大利亚气象部门每天发布的FFDIForest Fire Danger Index森林火险指数背后是McArthur Mark 5公式温度、相对湿度、风速、干旱因子四个量共同决定可燃物燃烧的剧烈程度。在AusFire类模型的参数文件里常见做法是把一天之内的气象站观测插值到模拟网格逐栅格算FFDI再套一层燃速公式得到每个位置的理论蔓延速度。气象站的CSV数据直接用readmatrix读进来下面这个函数是FFDI的最小实现。注意Matlab的log是自然对数不需要额外数学库function FFDI calcFFDI(T, RH, V, DF) % calcFFDI: McArthur Mark 5 火险指数 % T 气温(°C) | RH 相对湿度(%) | V 10m平均风速(km/h) | DF 干旱因子(0~10) FFDI 2 * exp(-0.45 0.987 * log(DF) ... - 0.0345 * RH 0.0338 * T 0.0234 * V); end这里有个数据单位陷阱相对湿度用百分比风速用km/h而不是m/s。拿到m/s的风速先乘3.6再传入否则FFDI会被系统性压低。DF干旱因子通常由KBDI干旱指数换算少数站点没有KBDI时用当地消防部门公开的每日FFDI产品反推DF这比从土壤湿度自己凑一个数稳妥。四个参数里DF的跨度最大对结果影响也最大它是模型的“底噪”。参数符号单位典型范围量级感气温T°C10~45每升10°CFFDI约放大1.4倍相对湿度RH%5~60每降10%约放大1.4倍风速Vkm/h0~60每增10km/h约放大1.26倍干旱因子DF无量纲0~10十倍跨度决定整体量级用一组夏季极端条件做个单点验证FFDI calcFFDI(35, 25, 30, 9); % 35°C / 25% / 30km/h / DF9这个算例对应干热大风天FFDI落在70~90区间已经进入澳大利亚火险等级里的最高档。拿你所在时区当天官方发布的FFDI对拍偏差在5以内说明气象数据单位处理没有问题偏差到了两位数先回头查风速和湿度单位再查DF是否来自同一套干旱指数体系。提示公式中的风速取10米高度平均风不是阵风。静态推演阶段用阵风直接代公式会把FFDI系统性抬高阵风只能放在逐分钟重算的动态场景里。2.2 用Matlab写燃速函数McArthur公式与坡度修正FFDI本身不等于蔓延速度还要过一道燃速公式。McArthur森林模型在文献里常见的分段形式是FFDI小于5走轻量级大于等于5走线性高段两段在5处连续。我习惯把坡度修正和可燃物载荷修正一起封装进同一个函数function ROS mcarthurROS(FFDI, slopeDeg, fuelLoad) % mcarthurROS: 林火蔓延速度, 单位 km/h % FFDI 火险指数 | slopeDeg 上坡坡度(°) | fuelLoad 可燃物载荷(t/ha, 可选) if FFDI 5 ROS 0.0012 * FFDI; else ROS 0.0027 * FFDI - 0.0075; end % 分段连续, FFDI5时两段都等于0.006 if slopeDeg 0 ROS ROS * exp(0.069 * slopeDeg); end if nargin 2 fuelLoad 0 ROS ROS * sqrt(fuelLoad / 2.0); % 以2 t/ha为基准 end ROS max(ROS, 1e-3); % 数值保护, 防止极低火险除零 end两个细节值得展开。坡度因子用指数形式20°坡时燃速放大exp(0.069×20)≈4倍这是山地火蔓延速度的主要推力可燃物载荷修正用平方根意思是载荷翻倍燃速只涨约41%不要高估载荷的作用。跑气候情景时fuelLoad从燃料图直接读DF和fuelLoad是两个独立变量前者表示燃料干湿程度后者表示单位面积可燃物总量混用是最常见的建模错误。2.3 风与坡不一致时有效风向量合成与反演标定坡度和风的方向经常不一致不能把风速和坡度分别乘进去完事。我一般把风向量和坡度向量做一次矢量合成坡度向量指向上坡方向、模长由坡度值决定风向量由风向风速决定合成后的方向作为椭圆长轴方向合成后的模长再映射回燃速。这里给出合成代码的骨架% 注意: 下面角度约定为x轴正方向0°, 逆时针为正 wx V * cosd(windDirCart); wy V * sind(windDirCart); sx 20 * cosd(slopeDir); sy 20 * sind(slopeDir); effX wx sx; effY wy sy; effDirCart atan2d(effY, effX); % 有效风向(笛卡尔角) effSpeed hypot(effX, effY); % 有效风速坡度向量模长取20理由是McArthur体系里20°坡爬升的燃速放大和20km/h风速贡献量级相当具体权重应该用所在区域的历史火行为反演确定。手里有几次已知火场蔓延记录时用优化工具箱的fminsearch以火头位置误差为目标函数反求这个权重比手调参数靠谱。气象站给的风向通常是“正北0°顺时针”的来向角要转成上面代码的笛卡尔去向角换算关系是取270−dir再mod 360或者读取时直接做90−dir两种写法等价关键在于全项目统一。3. 火线怎么长出来Huygens子波法与元胞自动机的Matlab实现3.1 两条路线的选型先跑通CA再升级到Huygens火线推进算法主流是两条路线Huygens子波法和元胞自动机CA。选型结论先给只想要燃烧面积曲线或教学演示用CA做应急业务、火头位置要精确到几百米用Huygens。AusFire这个方向公开的演示实现大多是CA但业务系统最终都会往子波法靠。两个维度的差别可以看表维度Huygens子波法元胞自动机火线形态连续光滑折线网格锯齿边界方向分辨率子波长轴连续可设受8邻域限制时间步长5~30s20~120s实现成本中上低典型定位业务推演快速预案、教学我一般建议先做CA拿到面积曲线和大致火场形态再写Huygens验证火头位置。这个顺序对调试也有好处两类算法的报错特征完全不同第5章的边界泄漏和时间步长问题分不清你在跑哪种算法时会浪费很多时间。3.2 Huygens子波法火线点集外推与凸包重采样Huygens子波法的核心假设是火线上每个点独立向外辐射一个椭圆子波椭圆长轴与有效风向一致长轴长度等于ROS×dt短轴按燃料类型取长轴的0.3~0.7倍所有子波的外包络就是下一时刻的火线。Matlab里外包络恰好有现成函数convhullfunction newFront huygensAdvance(front, effDir, ROSmax, dt) % front: Nx2火线点集(米) | effDir: 风去向(气象角, 北0°顺时) % ROSmax: 火头燃速(km/h) | dt: 时间步长(s) a ROSmax * 1000 / 3600 * dt; % 长半轴(m), 单位换算 b max(0.5 * a, 10); % 短半轴, 按燃料类型调整 theta deg2rad(90 - effDir); % 气象去向角转笛卡尔角 phi linspace(0, 2*pi, 36); % 每个子波36个采样点 R [cos(theta) -sin(theta); sin(theta) cos(theta)]; allPts zeros(0, 2); for i 1:size(front, 1) xy [a*cos(phi), b*sin(phi)] * R; allPts [allPts; xy front(i, :)]; end k convhull(allPts(:, 1), allPts(:, 2)); newFront allPts(k, :); newFront resampleFront(newFront, 200); end代码逻辑分四步换算长半轴并构造旋转矩阵每个火线点生成36个椭圆采样点并平移到该点全部采样点做凸包得到下一条火线最后重采样恢复均匀点距。重采样这步不能省凸包顶点在平直段会很稀、在尖角处很密点距不均匀会让下一步的子波发生重叠或空洞火线慢慢变形。重采样我按周长坐标插值function frontOut resampleFront(front, nPts) % 先把火线顶点闭合成环 front [front; front(1, :)]; d [0; cumsum(hypot(diff(front(:,1)), diff(front(:,2))))]; s linspace(0, d(end), nPts); frontOut [interp1(d, front(:,1), s), interp1(d, front(:,2), s)]; frontOut(end, :) []; % 去掉与起点重复的闭合点 end风向约定是整个函数最容易错的地方。气象站给的是“来向”即风从哪个方向吹来这里effDir要求的是“去向”指向火头移动方向两者差180°。算错的话整条火线会朝上风向长而且形态上看不出异常只能通过和卫星过火区对比才能发现。3.3 元胞自动机版本网格点火与方向权重CA实现更接近直觉每个栅格三态未燃、燃烧、燃尽。燃烧单元以一定概率点燃8邻域概率由局部燃速、时间步长、格网尺寸和风向夹角共同决定。下面是一个可直接循环调用的单步函数function [gridNew, areaBurned] caStep(grid, ROSmap, cellSize, dt, windTo) % grid: 1未燃 2燃烧 3燃尽 | ROSmap: 各格最大燃速(m/s) % cellSize: 格网边长(m) | dt: 时间步长(s) % windTo: 风去向(笛卡尔角, x轴正方向0°, 逆时针为正) [ny, nx] size(grid); gridNew grid; [rB, cB] find(grid 2); for k 1:numel(rB) r rB(k); c cB(k); for dr -1:1 for dc -1:1 nr rdr; nc cdc; if nr1 || nc1 || nrny || ncnx, continue; end if grid(nr, nc) ~ 1, continue; end ang atan2d(dr, dc); % 邻居方向角 cosw max(0, cosd(ang - windTo)); % 方向权重 p min(1, ROSmap(nr,nc) * dt / cellSize ... * (0.3 0.7 * cosw)); if rand p gridNew(nr, nc) 2; end end end end gridNew(grid 2) 3; % 燃烧格燃尽 areaBurned nnz(gridNew 2) * cellSize^2; % 过火面积(m^2) end点火概率里括号内的0.3是逆风下限用于近似飞火和辐射热的零星引燃这个系数在不同燃料图上要重新标定。方向权重只用了cos的一次方严格来说应该按椭圆方程投影但对30~60m的格网CA本身的离散误差远大于这个方向精度差别。如果风的数据是气象来向进函数前先加180°再换算成笛卡尔角换算规则和3.2节一致。3.4 两种算法共同的时间步长约束不管哪条路线都要守一条类似CFL条件的约束单步推进距离不能超过一个格网尺度量级。CA版本直接体现在概率上要求ROS×dt/cellSize远小于1取0.3以下安全Huygens版本要求a大于子波采样点间隔、同时小于火线局部曲率半径否则凸包会出现虚假尖角火线看起来“带刺”。按cellSize30m、ROSmax0.5m/s估算dt取20~30s跑6小时约720~1080步。每步对500×500栅格全扫描纯循环在Matlab里大约几秒总时间在可接受范围做参数敏感性分析时建议改自适应步长全局最大燃速低于阈值时放宽到60s火头接近高燃速区自动收紧能省一半以上的计算时间代价只是多写几行判定逻辑。4. 把AusFire结果变成汇报素材Matlab画图与PPT导出链路4.1 三张必出图热力图、火线序列、燃烧面积曲线模拟跑完先出三张图第一张用imagesc画燃烧状态热力图第二张把每15分钟的火线叠成序列图第三张画燃烧面积随时间曲线。这三张是澳大利亚山火PPT里最常用的标准配置一张交代空间、一张交代过程、一张交代量级。% 燃烧状态热力图 figure(Position, [100 100 800 500]); imagesc(xVec, yVec, burnGrid); axis xy; colormap(flipud(hot)); colorbar; xlabel(x / m); ylabel(y / m); title(AusFire 模拟末期燃烧状态); % 火线时间序列(每隔15分钟一条) figure; for i 1:numel(Fronts) plot(Fronts{i}(:,1)/1000, Fronts{i}(:,2)/1000, LineWidth, 1.2); hold on; end legend(cellstr(string(times, HH:mm))); % times为datetime数组 % 燃烧面积曲线 figure; plot(minutesElapsed/60, areaKm2, b-o, MarkerSize, 4); xlabel(时间 / h); ylabel(燃烧面积 / km^2); grid on;第二张图的图例用了datetime转string的写法string(times,HH:mm)在R2021b之后都支持。三张图的保存统一用exportgraphics分辨率给300这套函数对中文标题和坐标轴标签的字形支持比旧版print稳定。输出文件名里带上模拟时刻后面对PPT和复盘都方便fname AusFire_frame_ string(simClock, yyyyMMdd_HHmm) .png; exportgraphics(fig, fname, Resolution, 300);4.2 用Matlab直接生成PPT的三个办法把结果图拼成PPT常见的路线有三条mlreportgen.pp官方API、Windows下的COM自动化、先导出PNG再用python-pptx汇总。我一般优先mlreportgen.pp因为跨平台Linux和macOS也能批处理适合每天自动出一版报告。三者的取舍列成表办法跨平台批量能力主要坑mlreportgen.pp是强占位符名称依赖模板COM自动化仅Windows弱慢、依赖Office安装先PNG再python-pptx是强多一道文件交接mlreportgen.pp最小可用代码import mlreportgen.pp.* ppt Presentation(AusFire_Result.pptx); open(ppt); s add(ppt, Title and Content); replace(s, Title, AusFire模拟2020-01-06 14:00); img Picture(AusFire_frame_20200106_1400.png); img.Width 6in; img.Height 4in; replace(s, Content, img); close(ppt);两个高频坑第一占位符名称依赖所选模板默认模板图片位叫Content换成自定义模板后要先读模板的Placeholder列表再replace第二Picture的宽高最好不要同时写死先算好图片宽高比只设一个值否则图片会被拉伸变形。COM自动化在Windows上可以actxserver(PowerPoint.Application)交互性强但慢且依赖Office安装生产线环境很少用。4.3 澳大利亚山火PPT的页面结构建议面对非建模背景的汇报对象页面顺序建议按“结果倒推参数”安排先放最终燃烧范围与卫星影像叠加图再放燃烧面积曲线最后才放参数表。模拟动画放中间不要放第一页。参数表只列对结果影响最大的四项FFDI、风速风向、燃料载荷、坡度因子。燃料类型标签缺失时需要用遥感波段自己分类Matlab自带的kmeans函数可以直接用。把几个波段的反射率排成N×B矩阵kmeans聚类出3~5类分类结果转成燃料类型图再喂给模型比手工勾画边界快一个量级。聚类数量用silhouette轮廓图选不要拍脑袋定K值。5. 调试AusFire的三个高频问题边界泄漏、时间步长与燃料图分辨率5.1 边界泄漏与虚假尖角最典型的现象是火线接近模拟区边缘后趋向平直或者Huygens凸包在角点甩出很长的尖刺。前者是CA在边界外直接continue表现为“火到边界就停”后者是子波在角点的包络被人为截断。我一般做法是在网格四周补10圈不可燃单元模拟区外不做任何特殊处理出图时裁掉补边。这样CA不会在边界外静默停止Huygens的子波也不会在角点生成异常长轴。5.2 时间步长与库朗数确认逻辑没错之后火头推进距离明显超过ROS×dt或者火线出现锯齿优先缩短时间步长。经验值是库朗数不超过0.3dtSafe 0.3 * cellSize / max(ROS(:)); dt min(dtUser, dtSafe);ROS(:)里的最大值来自之前插值好的燃速场注意单位必须是m/s。自适应时间步长实现也不复杂每跑完一步重新求max(ROS)平稳段放宽到0.5高燃速段自动收紧。5.3 燃料图分辨率与投影对齐第三种高频问题不是算法而是数据。模拟火场与卫星过火区整体错位、且错位方向固定八成是投影坐标系不一致。GDA94/MGA分区和WGS84经纬度混用是这个方向的日常。Matlab里用readgeoraster读GeoTIFF检查RasterReference的坐标系再用projcrs按EPSG码创建投影对象做转换GDA2020/MGA zone 55对应EPSG 7855。不要手工加减经纬度偏移不同椭球体的差值不是常数。燃料图分辨率比模拟网格粗时插值前先想清楚变量性质。燃料载荷是连续变量griddata双线性插值没问题燃料类型是类别变量直接插值会插出不存在的中间类型。我通常把类别转成one-hot后再逐类插值最后取argmax决定类型。做完这一套火线形状和卫星观测的差异才真正落在模型参数上而不是数据管道上。最后补一个验证技巧不要只用燃烧面积曲线评判模型。面积可能因为网格粗而“看起来接近”火头位置和火线形状可能差很远。把模拟火线与卫星过火矢量叠加后计算双向Hausdorff距离这个指标对火头偏差远比面积误差敏感也是模型交付时最常被追问的一项。本文还有配套的精品资源点击获取