简介这是一套基于Matlab实现的澳大利亚山火模型AusFire源码与说明文档针对2020年澳洲山火场景进行建模适合防灾减灾、地理信息、气象环境等领域的科研人员、高校学生以及应急决策人员使用。资源包内共5个文件包含3个.m格式的Matlab源程序、1份PDF格式的模型说明文档及1个附加文件整体压缩包约12.92MB。源程序围绕火势蔓延过程进行数值仿真可结合风向、风速、温度、湿度、植被与地形等参数预测火线演化趋势为灭火策略制定与疏散方案评估提供科学依据。目前已有406人学习参考借助源码与说明文档读者可快速掌握AusFire模型的核心逻辑与参数调整方法并在此基础上开展二次开发或用于教学演示。1. 用MATLAB复现AusFire澳大利亚山火模型先要抓住这几点山火传播模型的价值不在“烧得准”而在“跑得快”。澳大利亚的AusFire模型正是围绕这个思路设计的它不追求单点燃烧的物理精度而是用半经验公式估计火线蔓延速度再在地理网格上推进边界。这个平衡点让AusFire能应用于数十公里尺度的火情推演也让它非常适合在MATLAB里做原型。常见做法是把地形高程、燃料类型、风速风向和湿度处理成矩阵然后逐时间步更新未燃单元运算逻辑和图像处理里的二值膨胀很像所以MATLAB的矩阵操作几乎是为这类模型量身定做的。这篇内容面向两类人一类是刚接触火灾模型、想把理论公式变成可跑代码的工程师另一类是已经有成熟模型、但希望用MATLAB做参数校准和成果展示的研究者。前者能从这里拿到一套最小可运行代码后者可以看到优化工具箱和PPT导出的衔接方式。注意AusFire本身并没有一套公开统一的官方实现不同项目里的版本差异很大所以这里以“能落地”为主用经典Rothermel方程和澳大利亚McArthur指数作为原型基础。2. AusFire模型的核心参数与MATLAB数学表达2.1 火线蔓延速度方程为什么选半经验形式AusFire类模型的底层逻辑是在某个地面单元火线蔓延速度由可燃物、气象和地形三个因子共同决定。最常用的原型是Rothermel速度方程虽然原始形式很长但工程实现时通常会把它分解成基准速度乘上若干无量纲修正系数R R0 * FW * FS * FM其中R0是基准蔓延速度零风、平坡、参考湿度下的速度FW是风速修正FS是坡度修正FM是燃料湿度修正。澳大利亚的McArthur森林火灾危险度尺FFDI也是类似思路只不过把温度、湿度、风速和干旱因子打包成一个指数。在MATLAB里做数值实现不需要纠结采用哪一套系数关键在于把这些修正项写成向量化函数能让矩阵上的每个单元并行计算。这里有一个常见误用把风速当常数。实际上风场在山地会随地形产生剧烈变化所以建议至少按坡向做一次风速遮蔽处理。MATLAB中可以用wind_shelter max(0, cos(slope_aspect - wind_direction))来近似把它乘到FW上。坡度修正则基于坡向与火线传播方向的夹角不是单纯的坡度值大小。2.2 关键参数表与取值范围参数符号典型范围单位MATLAB变量名基准蔓延速度R01 ~ 20m/minR0风速W0 ~ 50km/hWind风向WD0 ~ 360degWindDir坡度S0 ~ 45degSlope坡向Aspect0 ~ 360degAspect燃料湿度M2 ~ 30%Moisture燃料负载Load0.5 ~ 15kg/m²Load这些参数的敏感性并不相同。经验上风速对蔓延速度的非线性影响最强通常按(W / 30)^1.5作为FW燃料湿度超过15%后蔓延速度会急剧下降所以FM建议写成指数衰减形式FM exp(-0.15 * (Moisture - 5))上述公式在MATLAB里实现时要注意矩阵维度。所有修正系数必须保持和地形网格相同的尺寸否则矩阵运算会报维度错误。我一般会先读入数据再用meshgrid或直接依赖行列索引广播确保所有输入都是H x W的矩阵。2.3 用矩阵运算替代逐像元循环如果按照传统的C语言风格用双重循环遍历每个网格单元在MATLAB里会非常慢。正确做法是全部向量化。比如风速修正可以写成% Wind is a scalar in km/h, FW is a matrix of H x W speedRatio max(Wind / 30.0, 0); FW speedRatio .^ 1.5; % 注意点幂运算 FW FW .* (0.2 0.8 * cosd(WindDir - Aspect));这里cosd作用的矩阵MATLAB会逐元素运算。WindDir如果是一个常数它会被自动扩展成矩阵无需显式复制。关键在于.*和.^这两个点运算符缺一个就会变成矩阵乘法或矩阵乘方结果完全不对。这种坑几乎每个初做MATLAB的人都会踩所以调试时先检查变量尺寸再检查运算符。坡度修正FS需要坡度和坡向互动的方向效果。火线如果向下坡方向蔓延速度会减慢常用形式是FS (1 tan(degtorad(Slope)) .* cosd(Aspect - fireAdvanceDir)); FS max(FS, 0.05); % 防止出现负值或零值fireAdvanceDir是当前火线主传播方向在模型里通常初始由风决定后续根据局部火线法线更新。这样处理会让火线边界上不同位置的蔓延速度不同形成典型的舌头形状。3. 在MATLAB中写出最小可运行的AusFire传播代码3.1 网格初始化与燃料赋值我建议从一张简单的合成地型开始不要立刻读真实DEM。合成数据可以控制坡度、燃料和风场的行为方便验证模型是否按预期工作。下面代码生成一个 200x200 的网格中心有一块凸起地形燃料湿度左侧更高模拟一次从下风方向点燃的火灾。% 创建一个最小可运行的AusFire传播原型 nx 200; ny 200; [X, Y] ndgrid(1:nx, 1:ny); % 地形高程中心隆起其余平缓 h 20 * exp(-((X-100).^2 (Y-100).^2) / 2000); % 计算坡度简化直接用水平梯度 [dx, dy] gradient(h); Slope atand(sqrt(dx.^2 dy.^2)); Aspect atan2d(dy, dx); Aspect(Aspect 0) Aspect(Aspect 0) 360; % 燃料湿度左低右高模拟植被差异 Moisture 6 15 * (X / nx); % 燃料负载常数 Load 8 * ones(nx, ny); % 风参数 Wind 20; % km/h WindDir 90; % 从西向东吹 % 简化修正系数计算 FW (Wind / 30)^1.5; FS (1 tan(degtorad(Slope)) .* cosd(Aspect - WindDir)); FS max(FS, 0.05); FM exp(-0.15 * (Moisture - 5)); % 基准蔓延速度 R0 5.0; R R0 * FW * FS .* FM;参数说明ndgrid生成的X和Y矩阵尺寸是 200x200后面所有地形运算自动扩展到这个尺寸。gradient计算高程梯度时返回两个矩阵对应x和y方向坡度。atand和atan2d按 MATLAB 习惯使用角度制避免来回转换。这里的WindDir定义为风来的方向Aspect是坡向二者相减后取余弦表示“风顺着坡面吹”时的加速效果。3.2 时间推进与火线标记火线传播的本质是未燃单元被点燃的时间。这里采用一个简单的技术先给每个单元算出从初始火线传播到它所需的“最短旅行时间”。在MATLAB中可以用灰色距离变换或Dijkstra算法但作为原型元胞自动机更直观。每个时间步检查当前燃烧边界周围的未燃单元如果其累积传播时间达到阈值则点燃。% 初始化火线火线从左侧边缘点燃 ignited false(nx, ny); ignited(10:end-10, 1:3) true; % 燃烧状态矩阵0未燃, 1燃烧中, 2烧尽 state zeros(nx, ny); state(ignited) 1; burned state; % 模拟时间步每步代表1分钟 dt 1; time 0; tmax 60; % 邻域偏移8邻域 offsets [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1]; while time tmax % 找到所有燃烧中的单元 [row, col] find(state 1); growth false(nx, ny); for k 1:length(row) r0 row(k); c0 col(k); for oi 1:size(offsets,1) r r0 offsets(oi,1); c c0 offsets(oi,2); if r 1 r nx c 1 c ny state(r,c) 0 % 从当前燃烧单元到目标单元的速度 v R(r,c); % 目标单元的速度 travel_time sqrt((offsets(oi,1)^2 offsets(oi,2)^2)) * 1/v * 60; if travel_time dt growth(r,c) true; end end end end % 更新状态燃烧1分钟后变为烧尽 state(state 1) state(state 1) 1; % 新点燃 state(growth) 1; % 烧尽单元标记为2可选这里直接保留state0即已燃 time time dt; end这段循环性能不高但逻辑清晰。注意travel_time计算里乘了60因为速度R单位是 m/min而网格单元边长假设为 50 米常见分辨率。这里的对角邻域距离需要乘以sqrt(2)代码里用了sqrt((offsets(oi,1)^2 offsets(oi,2)^2))实现了这个权重。一个常用优化是直接用bwdist计算距离变换然后除以蔓延速度矩阵得到到达时间矩阵。这样可以将循环完全向量化但理解门槛更高。对于初学者先跑通循环再优化是合理路径。3.3 可视化与网格分辨率的影响用imagesc和contourf展示火线边界比直接画三维表面更直观figure; subplot(1,2,1); imagesc(state 0); axis off; colormap hot; title(最终过火区域); subplot(1,2,2); contourf(X, Y, R, 20); colorbar; axis equal; axis tight; title(蔓延速度分布);网格分辨率在这里是关键参数。如果单元边长从50米改成100米蔓延速度不变但每个时间步能穿越的单元数减半模型的时间步长也要相应调整。正确做法是固定“单步最大传播距离”即max_step dt * max(R(:))然后要求它小于网格单元边长。否则一个时间步内火线会跳过一整排网格出现不连续跳跃。一般建议保持dt * max(R) cell_size * 0.5。4. 用MATLAB优化工具箱校准AusFire参数并制作PPT成果4.1 定义目标函数与初始参数范围AusFire模型里的参数通常难以直接测量需要根据历史火情数据反演。常见做法是把风速修正指数、湿度衰减系数、坡度修正系数作为待优化变量用实际过火面积作为观测值。下面给出一个用lsqnonlin校准参数的最小框架。% 目标函数计算参数向量p下的模拟过火面积返回与观测值的残差 % 观测数据obs_burned_area (向量不同时刻) function resid calcResidual(p, meteo_data, terrain, obs) % 分解参数 windExp p(1); moistureCoef p(2); slopeCoef p(3); % 代入简化模型计算每时刻蔓延速度矩阵 R computeRos(meteo_data, terrain, windExp, moistureCoef, slopeCoef); % 用前面章节的传播代码得到每个时刻燃烧面积 sim_area runPropagation(R, meteo_data.time); % 残差模拟与观测的面积差 resid (sim_area - obs) ./ (obs 1); % 用百分比误差避免量纲差异 end % 初始参数猜测 p0 [1.5, 0.15, 1.0]; % 参数上下界 lb [0.5, 0.05, 0.2]; ub [3.0, 0.5, 3.0]; % 调用 lsqnonlin options optimoptions(lsqnonlin, Display, iter, MaxFunctionEvaluations, 200); p_opt lsqnonlin((p) calcResidual(p, meteo, terrain, obs_area), p0, lb, ub, options);参数说明p(1)是风速指数p(2)是湿度系数p(3)是坡度修正强度。lsqnonlin默认用信赖域反射算法对带边界约束的问题收敛较快。注意calcResidual里每次调用都会重新跑一遍传播模型因此最好把传播代码封装成一个独立函数并关闭所有图形输出否则优化会非常慢。这里的computeRos和runPropagation需要你根据前文思路自行封装。封装时可以把气象数据和地型变量统一放在一个 struct 里避免参数列表过长。经验上校准数据至少需要5个不同的风场条件否则会出现参数补偿现象比如风速指数增大、湿度系数减小两者互相抵消导致模型过拟合。4.2 用并行计算加速多组参数扫描校准过程中往往要跑几十次甚至上百次传播模拟。如果机器是多核MATLAB的parfor可以并行处理多组参数但lsqnonlin本身不支持并行评估目标函数。一个实用策略是先用patternsearch或遗传算法ga做全局粗搜再用parfor分散到多个worker% 使用parfor并行检查多个初始点 parpool(4); % 开启4个worker按需调整 paramSets [1.5 0.15 1.0; 1.2 0.2 0.8; 2.0 0.1 1.5]; residBest inf; pBest []; parfor i 1:size(paramSets,1) p_temp fminsearch((p) calcResidual(p, meteo, terrain, obs_area), paramSets(i,:)); resid_temp calcResidual(p_temp, meteo, terrain, obs_area); if sum(resid_temp.^2) sum(residBest.^2) pBest p_temp; residBest resid_temp; end end delete(gcp(nocreate));注意parfor循环里不能直接修改pBest需要在循环外用P parpool配合Parallel Pool的 Reduction 变量机制。上面代码里pBest和residBest作为 reduction 变量MATLAB 会为每个 worker 维护一份最后再合并。实际使用要注意parfor对变量分类的限制如果p_temp的定义依赖循环变量就允许但residBest这种累积量必须符合 reduction 规则。如果遇到错误优先检查是否所有 worker都能访问meteo和terrain等基础变量通常需要将它们放进基础工作区并确保不是匿名函数中的局部变量。4.3 把模型结果导出成PPT标题里“澳大利亚山火PPT”的常见做法有两种一是直接用copygraphics把MATLAB图形复制到剪贴板再手动粘贴二是用mlreportgen.ppt包自动生成整个PPT文档。后者适合需要输出十几页演示成果的场合。下面给出创建10页PPT并拆入图片和参数表的最小示例% 创建PPT生成器 import mlreportgen.ppt.*; slidesFile AusFire_Results.pptx; slides Presentation(slidesFile); open(slides); % 第1页标题页 titleSlide add(slides, Title Slide); replace(titleSlide, Title, 澳大利亚山火模型AusFire模拟结果); replace(titleSlide, Subtitle, sprintf(参数校准完成时间%s, datestr(now))); % 第2页插入模型参数表 tableSlide add(slides, Title and Table); tbl Table({参数名,优化值,初始猜测; ... 风速指数, num2str(p_opt(1)), 1.5; ... 湿度系数, num2str(p_opt(2)), 0.15; ... 坡度修正, num2str(p_opt(3)), 1.0}); tbl.Style {Border(solid), FontSize(12)}; replace(tableSlide, Table, tbl); % 第3页插入模拟火线图 figureSlide add(slides, Title and Picture); replace(figureSlide, Title, 模拟过火区域与速度分布); % 这里需要先导出当前图形为png exportgraphics(gcf, fire_spread.png, Resolution, 150); pic Picture(fire_spread.png); replace(figureSlide, Picture, pic); % 关闭并保存 close(slides);这段代码用到的mlreportgen.ppt在MATLAB R2019b之后的版本都内置无需额外工具箱。注意导出图片前要确保图形窗口已经绘制完成否则exportgraphics会抓到空白图。Picture对象可以指定宽度和高度比如pic.Width 5.5in避免图片撑破模板。如果PPT模板中某个占位符名称不对replace会报错可以先手动用PPT API查看add(slides, Title and Picture)返回的Block对象的Children列表来确认占位符名。5. AusFire模型验证技巧用图像分割结果对齐模拟火线验证AusFire模型时最容易被忽略的是空间精度。很多项目只对比过火总面积但两个面积相同的多边形可能形状完全不一致。一个实用方法是把模拟火线边界和从遥感影像分割出的真实火线图像做逐像素对比用Jaccard系数和轮廓距离评估。MATLAB图像处理工具箱提供了graythresh和imbinarize可以快速分割火线区域。真实火线通常用短波红外波段的亮温异常或归一化燃烧指数NBR来提取。这里假设你已经有了一个二值掩码observed_burn1表示烧过0表示未烧并且模拟结果sim_burn是同样尺寸的逻辑矩阵% 计算Jaccard交并比 intersection sum(sim_burn(:) observed_burn(:)); union sum(sim_burn(:) | observed_burn(:)); jaccard intersection / (union eps); fprintf(Jaccard系数: %.3f\n, jaccard); % 计算边界不匹配度对称距离 sim_edge edge(sim_burn, canny); obs_edge edge(observed_burn, canny); distSimToObs bwdist(obs_edge); distObsToSim bwdist(sim_edge); misalignment mean([distSimToObs(sim_edge(:)); distObsToSim(obs_edge(:))]); fprintf(平均边界距离: %.2f 像素\n, misalignment);Jaccard系数在0.7以上表示空间一致性较好0.5以下就说明模型参数还有明显偏差。平均边界距离则能反映系统性偏移如果真实火线总是比模拟火线偏北可能是风向场或坡度修正的方向定义反了。这里有一个很常见的坑MATLAB里的imagesc默认y轴向上而地理坐标通常y轴向北如果直接用矩阵行索引当作纬度模拟火线会比北方向偏移90度。解决方法是先把地理坐标翻转再计算距离sim_burn_geo flipud(sim_burn);另外bwdist计算的是欧几里得距离如果输入的是经纬度网格则需要先转换为投影坐标系。在澳大利亚范围内使用geoTrajectory或mapinterp可能更直接但最简单的方法是把经纬度格网用lla2utm需要Mapping Toolbox转成UTM投影然后再做像素对比。这样距离单位就是米而不是像素误差评估更符合实际。还有一种快速验证法模拟数百次随机参数组合把每次得到的Jaccard系数和对应参数绘制成散点图。这能直观看出模型对哪些参数不敏感。用scatter或boxchart展示参数敏感性比单独调参更高效。如果发现某个参数在很宽范围内Jaccard系数都差不多就说明该参数在当前数据下不可辨识应该从优化变量中剔除避免陷入局部极小。这个操作只需要几行代码对每个参数用linspace生成10个值固定其他参数在最优值上跑模型后记录Jaccard并绘图。放在整个AusFire模型开发的后期能大幅减少调参工作量。本文还有配套的精品资源点击获取 MATLAB实现AusFire山火模型：参数、代码与校准
简介这是一套基于Matlab实现的澳大利亚山火模型AusFire源码与说明文档针对2020年澳洲山火场景进行建模适合防灾减灾、地理信息、气象环境等领域的科研人员、高校学生以及应急决策人员使用。资源包内共5个文件包含3个.m格式的Matlab源程序、1份PDF格式的模型说明文档及1个附加文件整体压缩包约12.92MB。源程序围绕火势蔓延过程进行数值仿真可结合风向、风速、温度、湿度、植被与地形等参数预测火线演化趋势为灭火策略制定与疏散方案评估提供科学依据。目前已有406人学习参考借助源码与说明文档读者可快速掌握AusFire模型的核心逻辑与参数调整方法并在此基础上开展二次开发或用于教学演示。1. 用MATLAB复现AusFire澳大利亚山火模型先要抓住这几点山火传播模型的价值不在“烧得准”而在“跑得快”。澳大利亚的AusFire模型正是围绕这个思路设计的它不追求单点燃烧的物理精度而是用半经验公式估计火线蔓延速度再在地理网格上推进边界。这个平衡点让AusFire能应用于数十公里尺度的火情推演也让它非常适合在MATLAB里做原型。常见做法是把地形高程、燃料类型、风速风向和湿度处理成矩阵然后逐时间步更新未燃单元运算逻辑和图像处理里的二值膨胀很像所以MATLAB的矩阵操作几乎是为这类模型量身定做的。这篇内容面向两类人一类是刚接触火灾模型、想把理论公式变成可跑代码的工程师另一类是已经有成熟模型、但希望用MATLAB做参数校准和成果展示的研究者。前者能从这里拿到一套最小可运行代码后者可以看到优化工具箱和PPT导出的衔接方式。注意AusFire本身并没有一套公开统一的官方实现不同项目里的版本差异很大所以这里以“能落地”为主用经典Rothermel方程和澳大利亚McArthur指数作为原型基础。2. AusFire模型的核心参数与MATLAB数学表达2.1 火线蔓延速度方程为什么选半经验形式AusFire类模型的底层逻辑是在某个地面单元火线蔓延速度由可燃物、气象和地形三个因子共同决定。最常用的原型是Rothermel速度方程虽然原始形式很长但工程实现时通常会把它分解成基准速度乘上若干无量纲修正系数R R0 * FW * FS * FM其中R0是基准蔓延速度零风、平坡、参考湿度下的速度FW是风速修正FS是坡度修正FM是燃料湿度修正。澳大利亚的McArthur森林火灾危险度尺FFDI也是类似思路只不过把温度、湿度、风速和干旱因子打包成一个指数。在MATLAB里做数值实现不需要纠结采用哪一套系数关键在于把这些修正项写成向量化函数能让矩阵上的每个单元并行计算。这里有一个常见误用把风速当常数。实际上风场在山地会随地形产生剧烈变化所以建议至少按坡向做一次风速遮蔽处理。MATLAB中可以用wind_shelter max(0, cos(slope_aspect - wind_direction))来近似把它乘到FW上。坡度修正则基于坡向与火线传播方向的夹角不是单纯的坡度值大小。2.2 关键参数表与取值范围参数符号典型范围单位MATLAB变量名基准蔓延速度R01 ~ 20m/minR0风速W0 ~ 50km/hWind风向WD0 ~ 360degWindDir坡度S0 ~ 45degSlope坡向Aspect0 ~ 360degAspect燃料湿度M2 ~ 30%Moisture燃料负载Load0.5 ~ 15kg/m²Load这些参数的敏感性并不相同。经验上风速对蔓延速度的非线性影响最强通常按(W / 30)^1.5作为FW燃料湿度超过15%后蔓延速度会急剧下降所以FM建议写成指数衰减形式FM exp(-0.15 * (Moisture - 5))上述公式在MATLAB里实现时要注意矩阵维度。所有修正系数必须保持和地形网格相同的尺寸否则矩阵运算会报维度错误。我一般会先读入数据再用meshgrid或直接依赖行列索引广播确保所有输入都是H x W的矩阵。2.3 用矩阵运算替代逐像元循环如果按照传统的C语言风格用双重循环遍历每个网格单元在MATLAB里会非常慢。正确做法是全部向量化。比如风速修正可以写成% Wind is a scalar in km/h, FW is a matrix of H x W speedRatio max(Wind / 30.0, 0); FW speedRatio .^ 1.5; % 注意点幂运算 FW FW .* (0.2 0.8 * cosd(WindDir - Aspect));这里cosd作用的矩阵MATLAB会逐元素运算。WindDir如果是一个常数它会被自动扩展成矩阵无需显式复制。关键在于.*和.^这两个点运算符缺一个就会变成矩阵乘法或矩阵乘方结果完全不对。这种坑几乎每个初做MATLAB的人都会踩所以调试时先检查变量尺寸再检查运算符。坡度修正FS需要坡度和坡向互动的方向效果。火线如果向下坡方向蔓延速度会减慢常用形式是FS (1 tan(degtorad(Slope)) .* cosd(Aspect - fireAdvanceDir)); FS max(FS, 0.05); % 防止出现负值或零值fireAdvanceDir是当前火线主传播方向在模型里通常初始由风决定后续根据局部火线法线更新。这样处理会让火线边界上不同位置的蔓延速度不同形成典型的舌头形状。3. 在MATLAB中写出最小可运行的AusFire传播代码3.1 网格初始化与燃料赋值我建议从一张简单的合成地型开始不要立刻读真实DEM。合成数据可以控制坡度、燃料和风场的行为方便验证模型是否按预期工作。下面代码生成一个 200x200 的网格中心有一块凸起地形燃料湿度左侧更高模拟一次从下风方向点燃的火灾。% 创建一个最小可运行的AusFire传播原型 nx 200; ny 200; [X, Y] ndgrid(1:nx, 1:ny); % 地形高程中心隆起其余平缓 h 20 * exp(-((X-100).^2 (Y-100).^2) / 2000); % 计算坡度简化直接用水平梯度 [dx, dy] gradient(h); Slope atand(sqrt(dx.^2 dy.^2)); Aspect atan2d(dy, dx); Aspect(Aspect 0) Aspect(Aspect 0) 360; % 燃料湿度左低右高模拟植被差异 Moisture 6 15 * (X / nx); % 燃料负载常数 Load 8 * ones(nx, ny); % 风参数 Wind 20; % km/h WindDir 90; % 从西向东吹 % 简化修正系数计算 FW (Wind / 30)^1.5; FS (1 tan(degtorad(Slope)) .* cosd(Aspect - WindDir)); FS max(FS, 0.05); FM exp(-0.15 * (Moisture - 5)); % 基准蔓延速度 R0 5.0; R R0 * FW * FS .* FM;参数说明ndgrid生成的X和Y矩阵尺寸是 200x200后面所有地形运算自动扩展到这个尺寸。gradient计算高程梯度时返回两个矩阵对应x和y方向坡度。atand和atan2d按 MATLAB 习惯使用角度制避免来回转换。这里的WindDir定义为风来的方向Aspect是坡向二者相减后取余弦表示“风顺着坡面吹”时的加速效果。3.2 时间推进与火线标记火线传播的本质是未燃单元被点燃的时间。这里采用一个简单的技术先给每个单元算出从初始火线传播到它所需的“最短旅行时间”。在MATLAB中可以用灰色距离变换或Dijkstra算法但作为原型元胞自动机更直观。每个时间步检查当前燃烧边界周围的未燃单元如果其累积传播时间达到阈值则点燃。% 初始化火线火线从左侧边缘点燃 ignited false(nx, ny); ignited(10:end-10, 1:3) true; % 燃烧状态矩阵0未燃, 1燃烧中, 2烧尽 state zeros(nx, ny); state(ignited) 1; burned state; % 模拟时间步每步代表1分钟 dt 1; time 0; tmax 60; % 邻域偏移8邻域 offsets [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1]; while time tmax % 找到所有燃烧中的单元 [row, col] find(state 1); growth false(nx, ny); for k 1:length(row) r0 row(k); c0 col(k); for oi 1:size(offsets,1) r r0 offsets(oi,1); c c0 offsets(oi,2); if r 1 r nx c 1 c ny state(r,c) 0 % 从当前燃烧单元到目标单元的速度 v R(r,c); % 目标单元的速度 travel_time sqrt((offsets(oi,1)^2 offsets(oi,2)^2)) * 1/v * 60; if travel_time dt growth(r,c) true; end end end end % 更新状态燃烧1分钟后变为烧尽 state(state 1) state(state 1) 1; % 新点燃 state(growth) 1; % 烧尽单元标记为2可选这里直接保留state0即已燃 time time dt; end这段循环性能不高但逻辑清晰。注意travel_time计算里乘了60因为速度R单位是 m/min而网格单元边长假设为 50 米常见分辨率。这里的对角邻域距离需要乘以sqrt(2)代码里用了sqrt((offsets(oi,1)^2 offsets(oi,2)^2))实现了这个权重。一个常用优化是直接用bwdist计算距离变换然后除以蔓延速度矩阵得到到达时间矩阵。这样可以将循环完全向量化但理解门槛更高。对于初学者先跑通循环再优化是合理路径。3.3 可视化与网格分辨率的影响用imagesc和contourf展示火线边界比直接画三维表面更直观figure; subplot(1,2,1); imagesc(state 0); axis off; colormap hot; title(最终过火区域); subplot(1,2,2); contourf(X, Y, R, 20); colorbar; axis equal; axis tight; title(蔓延速度分布);网格分辨率在这里是关键参数。如果单元边长从50米改成100米蔓延速度不变但每个时间步能穿越的单元数减半模型的时间步长也要相应调整。正确做法是固定“单步最大传播距离”即max_step dt * max(R(:))然后要求它小于网格单元边长。否则一个时间步内火线会跳过一整排网格出现不连续跳跃。一般建议保持dt * max(R) cell_size * 0.5。4. 用MATLAB优化工具箱校准AusFire参数并制作PPT成果4.1 定义目标函数与初始参数范围AusFire模型里的参数通常难以直接测量需要根据历史火情数据反演。常见做法是把风速修正指数、湿度衰减系数、坡度修正系数作为待优化变量用实际过火面积作为观测值。下面给出一个用lsqnonlin校准参数的最小框架。% 目标函数计算参数向量p下的模拟过火面积返回与观测值的残差 % 观测数据obs_burned_area (向量不同时刻) function resid calcResidual(p, meteo_data, terrain, obs) % 分解参数 windExp p(1); moistureCoef p(2); slopeCoef p(3); % 代入简化模型计算每时刻蔓延速度矩阵 R computeRos(meteo_data, terrain, windExp, moistureCoef, slopeCoef); % 用前面章节的传播代码得到每个时刻燃烧面积 sim_area runPropagation(R, meteo_data.time); % 残差模拟与观测的面积差 resid (sim_area - obs) ./ (obs 1); % 用百分比误差避免量纲差异 end % 初始参数猜测 p0 [1.5, 0.15, 1.0]; % 参数上下界 lb [0.5, 0.05, 0.2]; ub [3.0, 0.5, 3.0]; % 调用 lsqnonlin options optimoptions(lsqnonlin, Display, iter, MaxFunctionEvaluations, 200); p_opt lsqnonlin((p) calcResidual(p, meteo, terrain, obs_area), p0, lb, ub, options);参数说明p(1)是风速指数p(2)是湿度系数p(3)是坡度修正强度。lsqnonlin默认用信赖域反射算法对带边界约束的问题收敛较快。注意calcResidual里每次调用都会重新跑一遍传播模型因此最好把传播代码封装成一个独立函数并关闭所有图形输出否则优化会非常慢。这里的computeRos和runPropagation需要你根据前文思路自行封装。封装时可以把气象数据和地型变量统一放在一个 struct 里避免参数列表过长。经验上校准数据至少需要5个不同的风场条件否则会出现参数补偿现象比如风速指数增大、湿度系数减小两者互相抵消导致模型过拟合。4.2 用并行计算加速多组参数扫描校准过程中往往要跑几十次甚至上百次传播模拟。如果机器是多核MATLAB的parfor可以并行处理多组参数但lsqnonlin本身不支持并行评估目标函数。一个实用策略是先用patternsearch或遗传算法ga做全局粗搜再用parfor分散到多个worker% 使用parfor并行检查多个初始点 parpool(4); % 开启4个worker按需调整 paramSets [1.5 0.15 1.0; 1.2 0.2 0.8; 2.0 0.1 1.5]; residBest inf; pBest []; parfor i 1:size(paramSets,1) p_temp fminsearch((p) calcResidual(p, meteo, terrain, obs_area), paramSets(i,:)); resid_temp calcResidual(p_temp, meteo, terrain, obs_area); if sum(resid_temp.^2) sum(residBest.^2) pBest p_temp; residBest resid_temp; end end delete(gcp(nocreate));注意parfor循环里不能直接修改pBest需要在循环外用P parpool配合Parallel Pool的 Reduction 变量机制。上面代码里pBest和residBest作为 reduction 变量MATLAB 会为每个 worker 维护一份最后再合并。实际使用要注意parfor对变量分类的限制如果p_temp的定义依赖循环变量就允许但residBest这种累积量必须符合 reduction 规则。如果遇到错误优先检查是否所有 worker都能访问meteo和terrain等基础变量通常需要将它们放进基础工作区并确保不是匿名函数中的局部变量。4.3 把模型结果导出成PPT标题里“澳大利亚山火PPT”的常见做法有两种一是直接用copygraphics把MATLAB图形复制到剪贴板再手动粘贴二是用mlreportgen.ppt包自动生成整个PPT文档。后者适合需要输出十几页演示成果的场合。下面给出创建10页PPT并拆入图片和参数表的最小示例% 创建PPT生成器 import mlreportgen.ppt.*; slidesFile AusFire_Results.pptx; slides Presentation(slidesFile); open(slides); % 第1页标题页 titleSlide add(slides, Title Slide); replace(titleSlide, Title, 澳大利亚山火模型AusFire模拟结果); replace(titleSlide, Subtitle, sprintf(参数校准完成时间%s, datestr(now))); % 第2页插入模型参数表 tableSlide add(slides, Title and Table); tbl Table({参数名,优化值,初始猜测; ... 风速指数, num2str(p_opt(1)), 1.5; ... 湿度系数, num2str(p_opt(2)), 0.15; ... 坡度修正, num2str(p_opt(3)), 1.0}); tbl.Style {Border(solid), FontSize(12)}; replace(tableSlide, Table, tbl); % 第3页插入模拟火线图 figureSlide add(slides, Title and Picture); replace(figureSlide, Title, 模拟过火区域与速度分布); % 这里需要先导出当前图形为png exportgraphics(gcf, fire_spread.png, Resolution, 150); pic Picture(fire_spread.png); replace(figureSlide, Picture, pic); % 关闭并保存 close(slides);这段代码用到的mlreportgen.ppt在MATLAB R2019b之后的版本都内置无需额外工具箱。注意导出图片前要确保图形窗口已经绘制完成否则exportgraphics会抓到空白图。Picture对象可以指定宽度和高度比如pic.Width 5.5in避免图片撑破模板。如果PPT模板中某个占位符名称不对replace会报错可以先手动用PPT API查看add(slides, Title and Picture)返回的Block对象的Children列表来确认占位符名。5. AusFire模型验证技巧用图像分割结果对齐模拟火线验证AusFire模型时最容易被忽略的是空间精度。很多项目只对比过火总面积但两个面积相同的多边形可能形状完全不一致。一个实用方法是把模拟火线边界和从遥感影像分割出的真实火线图像做逐像素对比用Jaccard系数和轮廓距离评估。MATLAB图像处理工具箱提供了graythresh和imbinarize可以快速分割火线区域。真实火线通常用短波红外波段的亮温异常或归一化燃烧指数NBR来提取。这里假设你已经有了一个二值掩码observed_burn1表示烧过0表示未烧并且模拟结果sim_burn是同样尺寸的逻辑矩阵% 计算Jaccard交并比 intersection sum(sim_burn(:) observed_burn(:)); union sum(sim_burn(:) | observed_burn(:)); jaccard intersection / (union eps); fprintf(Jaccard系数: %.3f\n, jaccard); % 计算边界不匹配度对称距离 sim_edge edge(sim_burn, canny); obs_edge edge(observed_burn, canny); distSimToObs bwdist(obs_edge); distObsToSim bwdist(sim_edge); misalignment mean([distSimToObs(sim_edge(:)); distObsToSim(obs_edge(:))]); fprintf(平均边界距离: %.2f 像素\n, misalignment);Jaccard系数在0.7以上表示空间一致性较好0.5以下就说明模型参数还有明显偏差。平均边界距离则能反映系统性偏移如果真实火线总是比模拟火线偏北可能是风向场或坡度修正的方向定义反了。这里有一个很常见的坑MATLAB里的imagesc默认y轴向上而地理坐标通常y轴向北如果直接用矩阵行索引当作纬度模拟火线会比北方向偏移90度。解决方法是先把地理坐标翻转再计算距离sim_burn_geo flipud(sim_burn);另外bwdist计算的是欧几里得距离如果输入的是经纬度网格则需要先转换为投影坐标系。在澳大利亚范围内使用geoTrajectory或mapinterp可能更直接但最简单的方法是把经纬度格网用lla2utm需要Mapping Toolbox转成UTM投影然后再做像素对比。这样距离单位就是米而不是像素误差评估更符合实际。还有一种快速验证法模拟数百次随机参数组合把每次得到的Jaccard系数和对应参数绘制成散点图。这能直观看出模型对哪些参数不敏感。用scatter或boxchart展示参数敏感性比单独调参更高效。如果发现某个参数在很宽范围内Jaccard系数都差不多就说明该参数在当前数据下不可辨识应该从优化变量中剔除避免陷入局部极小。这个操作只需要几行代码对每个参数用linspace生成10个值固定其他参数在最优值上跑模型后记录Jaccard并绘图。放在整个AusFire模型开发的后期能大幅减少调参工作量。本文还有配套的精品资源点击获取