简介元胞自动机CA结合Matlab实现的地区城市增长模拟项目面向城市地理、规划及GIS方向的学生与研究者用于理解土地利用变化和城市扩张过程。资源共12个文件包含7个m脚本、1个TXT运行说明、1张PNG结果图及3张测试图片压缩包仅105KB。核心主程序Main_code.m实现状态转移与迭代逻辑其余脚本覆盖影像读取、邻域筛选、道路距离计算、建成区提取等关键环节可直接在Matlab中复现。示例图像以不同颜色直观展示各时间步的土地利用格局运行说明则逐步交代参数设置与调用流程便于二次改造。已有265人学习适合希望借助完整CA框架快速搭建城市增长模型、开展预测实验的Matlab使用者。1. 把城市增长装进格子世界一份能跑的Matlab元胞自动机模拟方案拿到一份基于Matlab元胞自动机的城市增长模拟资源第一反应不该是“又一个跑图的Demo”而是先问一句它能不能解释“这块地为什么在第五年突然被开发”。城市扩张本质是一个自下而上的涌现过程每一块地的状态变化只取决于自身条件和邻居状态这正是元胞自动机的用武之地。这套资源的核心价值在于它把复杂的土地利用演变拆成“状态、邻居、规则、时间”四个可调变量让你在Matlab里跑一遍就能看到未来十几年的城市轮廓变化。适合做城乡规划课题、土地覆盖变化研究、或者单纯想理解CA模型怎么落地的从业者。下面从模型设定开始逐层拆解这套代码的真实用法。2. 元胞自动机的模型设定从栅格数据到状态转换2.1 状态空间与邻域类型Moore邻域为什么是默认选择元胞自动机里每一块地就是一个“元胞”它的状态通常是整数编码0代表不可开发区域水体、保护区1代表农业用地2代表已开发的城市用地3代表未来可能被开发的自然用地。这套编码不是随便定的后续的转换规则、概率计算、结果统计全部依赖这套数值体系。你拿到手的Matlab代码第一步就是确认输入矩阵的取值是否规范否则后面所有逻辑都会因为“出现了状态4”这种意外而翻车。邻域类型直接决定模拟结果的空间形态。常见的两种是Von Neumann邻域上下左右四个方向和Moore邻域周围8格。城市增长模拟里Moore邻域几乎是不二之选原因是城市开发倾向于沿已有建成区边缘蔓延8个方向的邻居能更好地捕捉“边缘效应”。代码里如果看到邻域半径参数默认值是1也就是只看紧邻的一圈这个参数不是不能改但改动后模拟结果会明显偏碎或者偏聚合属于典型的“微调出大偏差”。% 定义邻域偏移量Moore邻域半径1共8个方向 offsets [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1];这段代码的核心是偏移量矩阵每一行代表一个邻居相对当前元胞的坐标差。计算邻居时只要把当前坐标加上偏移量就能取出对应的邻居状态值。为什么不用全窗口扫描因为逐格扫描整个矩阵再判断邻居复杂度是O(NM8)而用偏移量配合边界判断写起来更直观也更容易和后面的并行化改造对接。2.2 转换规则的确定性与随机性开发概率如何驱动增长纯粹确定性的CA规则会得出机械化的棋盘格图案和现实城市相去甚远。实际项目几乎都采用“确定性约束 随机概率阈值”的组合方式。每一轮迭代元胞从“可开发”状态转为“已开发”状态的概率通常由三部分乘起来基础开发概率、邻居开发密度的加权影响、地形与政策约束的系数。% 开发概率计算伪代码逻辑真实实现见主程序 dev_prob base_prob * neighbor_factor * constraint_factor; if rand dev_prob current_state 1 grid_new(row, col) 2; % 状态2表示已开发 end转换规则是整套模型的灵魂也是最容易让新手产生“玄学”感的地方。开发概率的取值区间是0到1rand函数的输出也落在0到1之间当随机数小于开发概率时该格发生状态转换。这种随机性处理方式模拟的是现实中“开发商何时决定动工”的不确定性。参数标定时如果直接用经验值往往会出现模拟过快饱和或完全不长城的极端情况后文会专门讲怎么校准。2.3 初始数据准备把遥感分类图转成CA输入矩阵项目里能直接跑通的前提是你手上的数据格式得对得上。最常见的做法是用遥感分类结果比如把landsat影像在ENVI或ArcGIS里做监督分类得到土地利用类型栅格再在ArcGIS里把栅格转成ASCII文本导出。Matlab里用readmatrix读进来注意头部6行是投影信息要跳过。% 读取ArcGIS导出的ASCII栅格跳过元数据行 raw readmatrix(landuse.asc, NumHeaderLines, 6); landuse raw(:, 2:end); % 首列是列号标识实际栅格从第2列开始这里有个容易踩的细节ASCII栅格的第二行写的是ncols和nrows真正的栅格数据在头部信息之后而且不同软件导出的格式不完全一样——有的带逗号分隔符有的是空格分隔有的头部是6行有的是5行。在写读取代码之前先打开文件看三眼“头部几行”“分隔符是什么”“值域是什么”。宁可多花两分钟做数据探查也不要等跑出全零矩阵再回头排查。3. Matlab实现核心模块邻居遍历、状态更新与规则骨架3.1 读入数据与参数初始化把一切可变项集中到文件头部一份可维护的CA模型代码最重要的不是算法有多炫而是参数是否有“集中管理”的意识。拿到手的这份资源头部参数区写得很清晰包括模拟年份总数、行列数、初始状态矩阵路径、输出路径等。初始化时要把所有可变项拉到文件最前面避免在几百行代码里到处翻魔数。% 参数区 years_total 15; % 模拟15年每年迭代一次 rows size(landuse, 1); cols size(landuse, 2); base_prob 0.05; % 基础开发概率标定后确定 neighbor_weight 0.8; % 邻居开发密度权重 slope_limit 25; % 坡度超过25度禁止开发单位度 % 参数区写清楚每一行的含义是这份代码被反复复用的关键。某开发者曾经拿这份代码做课题时直接把neighbor_weight从0.8改成0.3模拟结果立刻从“沿主干道蔓延”变成了“四处飞地”后来对着参数区一条条检查才发现改错了地方。养成习惯每次调参前先复制一份原始参数区备份或者在代码里写清参数含义和推荐范围省掉后面的后悔药。3.2 并行迭代的核心循环为什么用双缓冲区CA模型的迭代有个经典陷阱如果直接在原矩阵上逐个更新那么前一个格子的新状态会污染后一个格子的邻居计算导致转换规则实际变成了“更新过的状态再去影响邻居”。正确做法是维护两个矩阵一个保存当前状态一个保存下一时间步的候选状态。每一轮迭代结束后把新矩阵赋给旧矩阵。grid landuse; % 当前状态 for t 1:years_total grid_next grid; % 复制一份作为缓冲 for i 2:rows-1 for j 2:cols-1 if grid(i,j) 1 % 只对可开发用地做状态判断 dev_prob compute_dev_prob(grid, i, j, offsets, params); if rand dev_prob grid_next(i,j) 2; % 更新到缓冲矩阵 end end end end grid grid_next; % 正式切换状态 % 每轮保存一张快照用于后续动图制作 save_snapshot(grid, t); end双缓冲区的意义在于“时间步的隔离”。如果不做这一步第一行第一列更新完后第二行第一列看到的邻居状态里就已经混入了“未来状态”模拟结果会偏向某个方向形态失真。另一个细节是循环边界从2到rows-1把边缘一圈留给边界处理避免索引越界报错。内层循环对每个格子的访问要做到O(1)取邻居值而不是每次重新从文件读入。3.3 约束条件叠加河流、保护区、坡度如何写入规则实际城市增长模拟很少只靠邻域概率基本都要叠加空间约束。比如河流水系两侧一定范围内不允许开发生态红线区域直接锁死状态坡度超过阈值的区域禁止转换。这些约束在代码里就是几个布尔索引矩阵。% 约束叠加禁建区河流、保护区、高坡度统一处理 forbidden is_water | is_protected | slope slope_limit; % 在迭代中forbidden格子的开发概率强制置0 dev_prob(forbidden_current) 0;约束叠加的处理时机有讲究是在计算开发概率之前就排除还是在计算之后强行归零常规做法是前者因为如果先算概率再归零那些高强度开发倾向的区域会反复被“尝试”浪费算力不说逻辑上也说不通——禁建区根本不在开发决策的考虑范围内。还有个容易忽略的细节约束矩阵如果和主栅格尺寸不一致Matlab不会报错而是直接自动扩展产生的模拟结果完全不可信。每次载入新约束数据后都要加一句矩阵尺寸断言。assert(isequal(size(forbidden), size(grid)), 约束矩阵尺寸与主栅格不一致);3.4 模型参数表一份可以直接抄的初始化配置参数是CA模型最说不清也最需要经验的部分。整理一份参数表有助于参照着调拿到资源后先按下表跑一遍确认结果符合常识再开始慢慢改参数名推荐初值作用调参方向base_prob0.05基础开发概率偏大则城市扩张过快偏小则基本不长neighbor_weight0.8邻居开发密度的权重偏大则聚集效应强偏小则产生飞地iteration_years15模拟时间跨度按研究需求设置过长需要验证规则稳定性slope_limit25度坡度约束阈值山区城市适当放宽neighborhood_radius1邻域半径大于1会显著增加计算量参数之间不是独立的。某开发者经常看到初学者只调base_prob结果模拟结果始终不对因为base_prob和neighbor_weight是乘积关系两个一起调到合理区间才能出现“沿路蔓延、边缘扩张”的自然形态。类似的模拟起始年份的数据和验证年份的数据要保持相同的分类标准否则状态编码含义不一致对不上历史趋势。4. 真实案卷演示模拟某区域未来15年的城市增长4.1 校准阶段历史数据反推开发概率直接拍脑袋定参数属于“自嗨式模拟”真正的项目流程是先校准再预测。常规校准方法是拿两期历史数据比如某地区2010年和2015年的土地利用图用2010年的数据作为初始状态模拟到2015年把模拟结果和2015年的真实数据做对比。通过多次调整参数让模拟结果的吻合度达到某个阈值比如整体精度80%以上此时的参数组合就是最终要用到未来预测的配置。% 校准流程跑一遍模拟对比真实结果 sim_2015 run_ca_simulation(landuse_2010, params); accuracy calc_overall_accuracy(sim_2015, landuse_2015); fprintf(整体精度 %.2f%%\n, accuracy * 100);校准阶段有个常见误区只盯着整体精度看。整体精度高不代表形态模拟得好因为有些区域压根不开发全预测为“保持原状”也能拿到高分。更稳妥的做法是同时检查Kappa系数和“新增城市像元的位置准确率”后者重点看新开发的格子是否落在真实增长区域附近。某开发者做过一次摸底发现只调整体精度到85%的模型在新增城市用地上实际命中率只有61%调整权重后才提到74%。4.2 执行模拟从2025推演到2040参数定下来后预测未来的操作本身很简单就是把初始状态换成最新一期建设用地数据把迭代年份设成目标年份差。但有一个关键点历史校准阶段用的转换规则固定不变等于假设“未来的城市开发行为和过去一致”。这个假设在短期内成立放到20年以上就靠不住了所以在执行推演之前要明确在文档里写下这个假设否则决策者拿你去预测一个25年后的结果本身就超出了模型的能力边界。% 预测模拟以2025年为起点推演至2040年 landuse_2025 readmatrix(landuse_2025.asc, NumHeaderLines, 6); landuse_2025 landuse_2025(:, 2:end); params.iter_years 15; future_result run_ca_simulation(landuse_2025, params);执行模拟的时候建议把每一年迭代的快照都保存下来不要只存最终结果。这样做有两个好处一是做动图时可以直接调用序列帧不用重新跑一遍二是如果发现某一轮迭代结果异常比如大面积同步开发能快速定位是哪个时间段出了非线性变化。存储格式上保存为逐年的TIF或者MAT文件都行但每轮都要在文件名里带上年份避免后面的可视化阶段去猜哪张图对应哪一年。4.3 结果验证Kappa系数与形态指标怎么看模拟跑了15年之后总得回答一个问题模拟得到的新城轮廓空间上到底靠不靠谱Kappa系数是常用的整体一致性指标能反映分类结果与真实数据的总体吻合程度但它是个“全局平均值”城市增长模拟里新增建设用地的位置准不准比整体分类准不准更重要。这时候要看两个补充指标新增元胞的命中率以及景观格局指数中的斑块数量和平均斑块面积。% 计算新增城市像元的命中率 sim_new (sim_result 2) (landuse_start ~ 2); real_new (landuse_real ~ 2) (landuse_end 2); hit_rate sum(sim_new(:) real_new(:)) / sum(real_new(:));如果命中率偏低先别急着调参数看一眼新增城市像元都跑到哪去了。常见情况是模拟结果把开发放在了大片连片区域而真实增长其实是沿道路呈枝状分布。这说明邻居权重和约束矩阵里缺少“交通廊道”这一层信息。交通可达性对城市扩张的导向作用非常明显纯粹的邻域模型很难复现枝状蔓延这时候要做的不是硬调参数而是在约束和概率层加入距离主干道的衰减函数属于模型本身的升级。5. 避坑手册Matlab城市增长模拟的五个翻车现场5.1 矩阵边缘出现了大片不开发区域模拟结果往往在栅格的四边出现一条整齐的“未开发带”看起来像是被人为裁掉了一圈。原因是邻居遍历时为了省事直接从2循环到rows-1边缘一圈元胞的邻居数量不足导致开发概率被系统性低估。更隐蔽的是某些角落位置只有3个邻居和中心区域的8个邻居相比开发概率天然偏低。解决方法是把边缘行和边缘列单独处理或者干脆在初始状态矩阵外围扩充一圈“值为-1”的虚拟元胞循环范围改到0到rows2碰到虚拟元胞直接跳过。% 外围填充虚拟元胞避免边界截断 grid_padded padarray(grid, [1 1], -1, both); for i 2:size(grid_padded,1)-1 for j 2:size(grid_padded,2)-1 if grid_padded(i,j) 1 % 处理邻居遇到-1视为不可开发 end end end填充虚拟元胞的方法还有一个好处它让所有内部元胞的邻居数目都严格等于8规则在空间上完全一致。如果你坚持用“跳过边缘”的方式后面做邻域半径为2的参数实验时边缘效应会被放大得更明显到时候再去补救就很被动了。5.2 模拟刚跑到第十年城市面积突然饱和这个现象多半是同步更新和随机比较之间出现了“正反馈失控”。场景是这样的某块地从农田转成城市后下一轮它成了邻居把周围一圈农地的开发概率全部顶高这些农地又变成城市再下一轮影响范围继续扩大。如果neighbor_weight设置过高这个扩散过程就是指数级的模拟结果在几年内飞速饱和整张图几乎全红。要判断是不是这个原因把neighbor_weight降到0.5以下重新跑一轮如果曲线突然平缓那说明反馈强度爆表了。另一个常见原因是每轮迭代里rand重估次数过多——有些实现会把转换判定拆成两三次随机判断导致累积概率远超名义值。检查代码里是否同一个格子在一轮里做了多次rand比较有就合并成一次。5.3 模拟结果每次都完全不一样无法复现CA模型的随机性来自rand函数而Matlab的rand是伪随机数依靠种子seed启动序列。如果不设定随机种子每次运行的结果固然有差异这个差异在研究场景下没问题但要做参数对比实验、重复实验或者给导师复现时同一套参数得出完全不同的两张图会让人直接质疑你的结果。解决办法很简单在代码最开头固定随机种子。% 固定随机种子保证模拟结果可复现 rng(2024001);固定种子之后同一套参数永远跑出同一套结果这不是为了“作弊”而是为了控制变量。做灵敏度分析时每次只改一个参数其他条件包括随机种子保持不变才能确定结果差异确实来自于被改的参数而不是随机波动。另一个相关的小技巧是如果你怀疑某一组参数的结果是“运气好才好看”换五个不同的种子跑五遍看形态特征是否稳定这比单次结果更有说服力。5.4 小型数据跑得动换成大范围数据直接内存溢出城市增长模拟的研究范围常常是几十公里乘几十公里栅格分辨率如果做到30米矩阵规模直奔百万级。Matlab里对double类型的矩阵做zeros(rows, cols)8字节一个元素几千万个元素就是几百MB。而这份代码在迭代中同时存在2到3个同尺寸矩阵内存压力翻倍。先看数据能不能精简状态矩阵用uint8类型存储值域0到255完全够用内存占用直接降为原来的八分之一。再看循环写法如果内层用了逐格索引三次以上同一个大矩阵Matlab的访存开销会非常高。保守做法是先把变量类型都检查一遍再考虑是否重构成矩阵化操作。% 使用uint8压缩状态矩阵 landuse uint8(landuse); grid landuse;类型转换这个动作在代码里只占一行但能解决大部分中小规模数据的内存瓶颈。改类型后原来那些和double类型的常量比较是否受影响Matlab会自动做类型提升比较和赋值都没问题。真正要注意的是算术运算的结果类型uint8加uint8有可能溢出涉及数学计算时先转成double。5.5 转换规则没生效约束区域仍然出现了开发遇到过好几次类似反馈“我把保护区在代码里置为不可开发了跑完发现保护区里全是新建城区”。检查代码后找到的根因是读取约束数据时执行了find(forbidden 1)但在迭代循环里条件判断写成了if grid(i,j) 1 dev_prob 0并没有把约束矩阵挂进去。说得直接点就是规则绕过了约束约束只在初始化时“看了一眼”后面再没参与过判定。正确做法是每轮迭代都实时计算约束掩膜并且把约束判断放在开发概率计算之前作为一票否决项。代码逻辑上要写成先查表判断该格是否允许开发再算概率再掷骰子三步各自独立顺序不能乱。% 正确的约束判定顺序不可开发则直接跳出 if forbidden(i,j) continue; end dev_prob compute_dev_prob(...); if rand dev_prob grid_next(i,j) 2; end6. 一个压箱底技巧把逐年快照做成带图例的动图让结果看得见也讲得清模拟跑完之后光有一张终态图不够团队汇报时最容易被追问的问题是“这个过程是怎么演变的”。动态可视化不只是好看它能帮你发现空间演变中的突变时刻比如第五年到第七年之间是不是出现了某个方向的加速蔓延。这里给出一个我常用的制作时间序列动图的方法用Matlab的imwrite把每年的结果写成PNG帧再用ffmpeg合成视频。% 逐年保存为PNG帧颜色固定图例在最后统一制作 cmap [0.2 0.6 0.2; 0.9 0.8 0.2; 0.6 0.2 0.2]; % 绿、黄、红 for t 1:years_total img uint8(255 * reshape(cmap(snapshot(:,:,t)1, :), rows, cols, 3)); imwrite(img, sprintf(frame_%02d.png, t)); end这个代码块的核心有三个要点。一是颜色映射表每一行对应一种状态状态值加1作为索引这样状态0是绿色状态1是黄色状态2是红色视觉上“从绿到黄再到红”正好对应未开发到已开发的演变。二是imwrite直接写PNG比saveas绕一圈figure快很多而且不占屏幕资源。三是文件名用%02d补零保证后续按字典序拼接视频时帧顺序不会跳乱。合成视频在Windows和Mac下略有不同Windows下用system(ffmpeg -framerate 4 -i frame_%02d.png -c:v libx264 output.mp4)Mac下先确认ffmpeg版本和路径。如果不想依赖ffmpeg也可以在Matlab里直接用VideoWriter类逐帧写入但代码会多几行。常用的配置是每帧停留0.25秒4帧率15年的模拟结果正好是4秒钟的动图汇报时反复播放三遍都不会觉得冗长。从这份资源的实际使用经历来看最稳的工作流是第一遍用默认参数跑通第二遍用历史数据校准第三遍带着约束条件做未来推演最后再做动图。我曾经因为跳过校准直接拿默认参数模拟未来20年结果输出的建成区面积比真实面积还大被同组的人点名质疑。从那以后我每次拿到CA模型代码第一件事就是在参数区写清“来源”和“标定日期”任何参数必须能往回查到出处没有标定的参数宁可先用保守值也不要拍脑袋填一个好看的数。希望这份拆解能帮你顺利跑出自己的城市增长模拟也希望你在调参的路上少踩几个和曾经的我一模一样的坑。本文还有配套的精品资源点击获取