简介面向城市规划与地理信息建模场景基于元胞自动机的Matlab项目围绕艾哈迈达巴德地区构建城市增长预测模型利用离散网格状态和邻域规则迭代模拟住宅、商业等土地利用类型的演化趋势。适合希望将CA理论应用到实际区域模拟的研究者、地理信息专业学生也适合初步接触Matlab建模的入门者。资源共12个文件压缩包仅105KB主体为7个m脚本覆盖主程序、图像读取、邻域过滤等核心功能另配jpg/png图像作为输入样例与结果展示txt文件说明运行方法整体结构清晰。目前已有265人学习。通过该项目可掌握用Matlab搭建元胞自动机状态转移逻辑的方法理解交通邻近度、邻域范围等因素对城市扩张的影响同时获得模型调参和扩展的基础思路适合作为相关课程设计或小型科研实验的起点。1. 为什么要用元胞自动机做城市增长模拟城市蔓延是个典型的空间过程新增建设用地往往贴着已有建成区向外扩张邻域里的开发强度越高下一时段这里被开发的可能性就越大。元胞自动机Cellular Automata简称 CA恰好能把这种“局部相互作用决定全局演变”的规律讲清楚模型简单、机理直观尤其适合在栅格数据上逐像元推演未来几十年的城市形态变化。用 Matlab 做这件事优势在于矩阵运算天然贴合栅格结构代码量比用通用 GIS 脚本少很多还能灵活接进坡度、水域、生态保护区、规划管控这些约束条件。这套方法在土地利用模拟、规划方案比选、生态风险预警里用得非常多。读这篇文章的人多半是刚接触 SD、FLUS、CA-Markov 这类模型的学生或规划方向的从业者手头可能已经有一份土地利用数据想在未来情景下看看城市会往哪边长、长多少但又不想一上来就投入复杂商业软件。我们就用纯 Matlab 把整个流程走通从数据准备到参数率定再到未来情景模拟和结果验证每一段都给可直接抄的代码和参数说明。城市增长 CA 建模的核心不是代码而是三个问题怎样定义元胞状态和邻域怎样表达开发适宜性以及怎样让模拟结果经得起检验。下面按建模、实现、率定、避坑、进阶的顺序一层层拆开讲。2. 核心概念与模型框架状态、邻域、转换规则与约束城市增长 CA 在原理上并不复杂研究区域划分成一个个规则格网每个格网就是一个元胞。元胞的状态通常定义为“非城市”或“城市”时间步长一般取 1 年每个时间步内所有元胞根据一套转换规则决定是否从非城市转为城市。设计这套规则需要想清楚元胞邻域、开发概率计算、约束条件、随机扰动四个部分。下面逐个说。2.1 元胞状态与邻域结构模型里每个格网像元都有一个状态值。我一般用 0 表示非城市用地1 表示城市用地水域和永久性保护区用特殊标记例如 -1在计算时直接锁定为不可开发。邻域结构最常用的是 3×3 摩尔邻域也就是中心像元周围一圈的 8 个像元也可以是 5×5 或 7×7。理论上邻域越大越能反映区域尺度的城市辐射效应但计算量和内存也会随之增长。对 30 米分辨率的栅格我一般用 3×3 邻域因为城市扩张的主要驱动力来自紧邻的已开发区域再远的像元影响已经很弱设置过大的邻域反而会带来空间平滑过度、城市形态失真。邻域影响因子的计算是核心定义为中心像元邻域内城市像元数占有效像元总数的比例范围从 0 到 1。这个值越高说明周围建成环境越成熟中心像元被开发的概率也越大。它反映的是“集聚效应”——城市开发有很强的路径依赖跟着已有设施走。2.2 开发概率怎么算自然因子、约束因子与随机扰动每个非城市元胞在每个时间步都有一个开发概率 P它是自然适宜性、邻域影响和约束条件共同作用的结果。我采用一种加权乘积的形式P (1 邻域发展密度 × 邻域权重) × 自然适宜性 × 各约束条件系数 × (1 随机扰动参数)自然适宜性来自坡度、高程、距道路距离、距建成区距离等静态因子。多个因子归一化后乘以各自的权重相加得到一个 0 到 1 的适宜性得分。约束系数则是二值乘数若像元位于水源保护地、生态红线、基本农田等不可开发区域该系数为 0若在允许开发区内则为 1。随机扰动参数是一个 0 到 0.05 的均匀随机数用来体现规划之外偶然因素带来的开发行为。2.3 迭代过程每年重算一次直到跑完整个时段模拟从初始年份的土地利用图开始按时间步逐年更新每一年先重算每个非城市像元的开发概率再通过概率阈值判断是否开发更新状态后进入下一年。以模拟未来 20 年为例总共跑 20 个时间步。每一年城市像元都在增多邻域密度会变化所以下一年的概率也会随之调整整个过程是动态的。开发概率判断可以采用随机阈值法对每个非城市像元生成一个 0 到 1 的随机数当 P 大于该阈值且 P 超过一个全局最低开发概率常数时才转为城市。也可以用“按概率排序取前 N 个像元”的方式即每轮从候选集合中选概率最高的 N 个像元转成城市这个 N 由历年建设用地增长率决定。随机阈值法更贴近实际因为城市开发不是完全按区位优劣排序的。我常用随机阈值法。2.4 参数标定不能拍脑袋定权重CA 模型最大的争议点在于权重参数。直接用主观经验定权重模拟结果常常不可靠。常见做法是用历史数据反演拿过去某两个年份的真实土地数据按同样的规则模拟中间年份不断调整各因子权重直到模拟结果与真实期末年份数据的 Kappa 系数达到可接受水平一般 0.75 以上。这个过程叫参数率定在第 4 章详细展开。3. 用 Matlab 实现城市增长 CA数据准备与核心代码Matlab 实现 CA 城市模拟我的代码框架包含数据读取、参数初始化、概率计算、状态迭代、结果输出五个模块。为了便于理解我用一个模拟项目 X 的简化数据作为演示实际应用时把数据路径和文件名替换掉即可。3.1 输入数据准备读栅格、统一投影与分辨率首先要准备以下栅格图层初始年份土地利用图如 2010 年城市像元赋值为 1非城市为 0水域为 -1目标年份土地利用图如 2020 年用于参数率定和精度验证DEM 高程数据用于计算坡度道路、河流矢量数据需要转为距离栅格规划管控边界生态红线、基本农田、建设用地管制区转为二值栅格所有栅格必须保持同样的范围、分辨率如 30 米×30 米和行列数。如果投影或分辨率不一致先在 ArcGIS、QGIS 或其他工具里统一重采样和投影否则 Matlab 里叠加后矩阵尺寸不匹配会出现大量 NaN 或错位。这一步最耗时但不可跳过。读取数据我建议读 TIFF 或 IMG 格式% 读取土地利用数据和DEM landuse_2010 readgeoraster(landuse_2010.tif); dem readgeoraster(dem.tif); % 检查行列数不一致时以 landuse_2010 为基准重采样 [rows, cols] size(landuse_2010); dem imresize(dem, [rows, cols], nearest);逻辑说明readgeoraster 是 Matlab 较新版读取栅格的标准函数读入结果是二维矩阵。imresize 的 nearest 方式可以保持原始类别值不变避免像双线性或三次插值那样产生中间值这对土地利用分类数据很重要。矩阵的行列顺序对应地理坐标第 1 行是研究区最北边第 1 列是最西边后面做空间分析时要记住这个约定。3.2 计算约束与适宜性因子以坡度为例从 DEM 用 gradient 函数计算坡度是一个常见快速做法% 计算坡度的近似值并归一化 [gx, gy] gradient(dem, cellsize, cellsize); slope atand(sqrt(gx.^2 gy.^2)); slope_norm 1 - (slope - min(slope(:))) / (max(slope(:)) - min(slope(:)));参数说明cellsize 是像元尺寸比如 30。坡度的归一化采取了反向映射坡度越大、适宜性得分越低。gradient 方法对 DEM 噪声比较敏感正式研究中我推荐在 GIS 里用专业水文工具生成坡度栅格再通读入。Matlab 里的 gradient 方法只适合数据质量较好、分辨率一致的快速原型验证。道路距离因子同样用距离变换函数% 道路二值图1为道路 road_dist bwdist(road_binary); % 归一化距离越近得分越高 max_dist max(road_dist(:)); road_factor (1 - road_dist / max_dist) .* (road_dist 0);3.3 核心模拟主循环模拟部分的核心循环如下。这里约定 land_state 为当前年份的状态矩阵其值为 1 表示城市、0 表示非城市、-1 表示不可开发区域。每次迭代对非城市且非保护的像元计算概率并判断是否转换。% 参数设置 iter_years 20; % 模拟年数 threshold 0.35; % 基础开发阈值经率定确定 neighborhood_weight 0.7; % 邻域影响权重 random_mag 0.03; % 随机扰动幅度 urban_current (land_state 1); % 当前城市像元二值图 excluded (land_state -1); % 不可开发区域 for t 1:iter_years % 计算邻域发展密度 3x3 kernel ones(3,3); kernel(2,2) 0; % 排除自身 neigh_density conv2(double(urban_current), kernel, same) / 8; % 只对非城市、非排除区域计算 candidate ~urban_current ~excluded; % 自然适宜性 slope_norm .* road_factor; 约束 protection_mask suitability slope_norm .* road_factor .* protection_mask; % 邻域影响系数范围0~1.7中心像元非城市时邻域无自身贡献 neighbor_effect 1 neighborhood_weight * neigh_density; % 随机扰动 random_factor 1 random_mag * rand(size(urban_current)); % 综合开发概率 dev_prob neighbor_effect .* suitability .* random_factor; dev_prob(~candidate) 0; % 候选区外概率置0 % 随机阈值比较 random_threshold threshold * rand(size(urban_current)); newly_built (dev_prob random_threshold) candidate; % 更新状态 urban_current urban_current | newly_built; land_state(newly_built) 1; end逻辑说明循环里的步骤其实就是第 2 章公式的代码化实现。conv2 函数用 3×3 邻域核对城市二值图做卷积得到每个像元周围的城市数量除以 8 就得到邻域发展密度。suitability 把自然适宜性统一成一个 0 到 1 的矩阵protection_mask 是二值矩阵生态保护区和基本农田对应的像元为 0该区域开发概率最终为 0。随机扰动项每次迭代重新生成保证每个像元的运气不同。参数说明threshold 是全局开发阈值该值越低开发越容易发生。它的标定依赖于历史数据的逐年变化规律一般初始设为 0.3 到 0.5再在率定阶段精细调优。neighborhood_weight 控制城市集聚效应强度权重越大、城市邻域对开发的影响越显著。在扩张速度过快、城市过度碎片化的场景中优先调高这个权重在模拟结果过于集中时再适当调低。3.4 结果输出与可视化模拟结束时需要把结果保存为 GeoTIFF便于 GIS 软件打开和制图。Matlab 自带 geotiffwrite 函数需要传入坐标参考对象。% 构造地理参考对象沿用读取初始数据时的ref ref readgeoraster(landuse_2010.tif); % 仅取ref信息 geotiffwrite(urban_2030.tif, double(urban_current), ref, CoordRefSysCode, 32650);更稳妥的写法是在一开始读数据时同时获取空间参考信息使用标准格式的读法将 R 对象保存下来[landuse_2010, R] readgeoraster(landuse_2010.tif); % 模拟结束后 geotiffwrite(urban_2030.tif, double(land_state), R);变量 R 来自 readgeoraster 的第二输出是地理栅格参考对象包含了像元分辨率、研究区范围、投影信息写入文件后不会丢空间位置。至此一套完整的最小 CA 城市模拟就跑通了。接下来的问题在于结果的可靠程度如何。4. 参数率定与结果验证把模拟精度调到可用水平不做率定直接跑未来情景结果基本只能作为课堂练习不能用于规划决策。原因很简单初始参数是拍脑袋定的没有依据。要通过历史数据把参数校准到合理值让模型回算历史年份时误差尽量小。4.1 率定思路找出最优参数组合率定的核心思路是“反演验证”已知研究区 2010 年和 2020 年两期真实土地利用图就用 2010 年做初始状态用一组参数模拟 2020 年再把模拟结果与 2020 年真实数据进行对比统计空间一致程度依据误差修正参数反复迭代。需要率定的参数主要是三个基础开发阈值 threshold、邻域影响权重 neighborhood_weight、随机扰动系数 random_mag。坡度、道路距离因子的权重如果还没有把握也需要加入率定但每增加一个参数组合空间都会急剧膨胀。我习惯先固定静态因子权重用熵权法或两两比较法只率定上述三个动态参数效率更高、结果也更可解释。4.2 Kappa 系数与精度检验精度检验两组指标足够像元级一致性指标 Kappa 系数、形态级指标城市斑块数量和平均斑块面积。Kappa 系数是实打实的量化标准公式比较繁琐但 Matlab 代码很短。% predict_map: 模拟期末城市分布二值图 % actual_map: 真实期末城市分布二值图 % 转为逻辑型 pred predict_map 0; act actual_map 0; % 混淆矩阵 n numel(pred); a sum(pred act); % 正确预测为城市 b sum(pred ~act); % 过度预测 c sum(~pred act); % 遗漏预测 d sum(~pred ~act); % 正确预测非城市 p_obs (a d) / n; % 总体一致率 p_exp ((ab)*(ac) (cd)*(bd)) / n^2; % 随机一致率 kappa (p_obs - p_exp) / (1 - p_exp);p_obs 表示的是预测与真实情况完全一致的像元占全部像元的比例p_exp 则是在类别比例不变的前提下随机摆放会产生的期望一致率。Kappa 大于 0.8 说明模拟结果非常好0.6 到 0.8 属于可用水平。如果只有 0.4 以下就需要检查模型结构或数据是否存在严重问题。4.3 参数调节经验先粗调后精调参数粗调看三个方向的信号。如果模拟期末城市像元数远大于真实值说明 threshold 设得太低或 neighborhood_weight 太小城市扩张过快如果城市主要集中在老城区周边呈紧凑圈层蔓延而真实数据却出现了新区跳跃式开发说明随机扰动太小无法模拟远距离新城开发如果模拟结果出现了大量零星孤立像元真实数据却没有那么邻域权重偏低缺少集聚约束。精调阶段可以用网格搜索或粒子群算法自动寻参。网格搜索太费时间粒子群可以较快逼近最优参数组合。这里给一个粒子群调参的最小框架% 目标函数输入参数向量x [threshold, neighbor_weight, random_mag] % 返回负kappa值因为particleswarm是求解最小化 kappa_neg fun(x); options optimoptions(particleswarm, SwarmSize, 20, MaxIterations, 30); best_param particleswarm((x) my_ca_kappa(x, init_state, actual_2020), 3, ... [0.1, 0.3, 0.0], [0.8, 0.9, 0.1], options);逻辑说明particleswarm 是 Matlab 自带的粒子群优化函数它需要目标函数、变量维度、变量下界和上界。这里把模拟 2020 年的整个 CA 过程封装进 my_ca_kappa返回负 Kappa 值作为适应度。每一次粒子移动都会驱动模型完整跑一遍历史反演计算量比较大所以先把栅格降采样到 100 米左右用于率定找到合适参数后再放回全分辨率做最终模拟能省出大量时间。4.4 一个容易被忽视的率定细节时间步对齐历史反演必须保证模拟的时间步长和真实数据间隔一致。如果 2010 到 2020 是 10 年的间隔模型一个时间步代表一年那么就要循环 10 次而不是只循环 1 次。很多新手一开始只跑一个时间步就拿来对比得出的模型几乎必然低估城市扩张总量因为正常扩张应当是累积了 10 年的结果。这个细节看起来简单但实际项目里翻车的概率极高。5. 常见问题与避坑清单从数据到结果的实用排错以下问题都是我实际跑项目时遇到的每条都按现象、原因、解决写清楚方便排查。5.1 数据投影和分辨率不一致导致的错位现象是模拟结果在空间上明显偏离真实分布比如河流位置突然对不上城市像元出现在山上。原因是不同来源的 DEM、土地利用和道路数据没有统一投影和分辨率导致矩阵叠加时位置错位。解决的唯一办法是回到 GIS 里把所有数据统一到一个投影坐标系如 UTM统一重采样到相同像元大小再用 Matlab 读入检查行列数是否一致。这一步是建模前的前置门禁规范的数据整理流程在这类空间模拟里不是锦上添花而是开局的硬性依赖。5.2 城市增长一上来就爆炸式扩散现象是模拟前几年新增城市像元数量巨大远超真实水平随后增长又迅速停滞。原因是基础开发阈值 threshold 设得偏低加上初始城市附近邻域密度高导致第一轮大量像元同时转为城市。解决方法是改用带“最高年增长上限”的约束方案每年允许新增的城市像元数不超过前一年城市像元总数的一个比例比如 2% 到 5%。代码里加一个判断把超过上限的开发像元改回非城市状态。% 每年新增像元上限按当前城市像元数的3% max_new round(sum(urban_current(:)) * 0.03); new_indices find(newly_built); if length(new_indices) max_new % 只在开发概率高的像元中保留前max_new个 [~, order] sort(dev_prob(new_indices), descend); keep new_indices(order(1:max_new)); newly_built false(size(urban_current)); newly_built(keep) true; end这个上限在真实城市增长建模中有实际依据城市建设用地的扩展受年度用地指标、市场去化速度等多重因素约束对一段时间内的扩张速度形成了有形或无形的限制。加入这个限制后模拟结果更符合真实增长节奏。5.3 模拟结果出现“椒盐噪声”般的孤立城市像元现象是模拟图里出现大量零星分布的单一城市像元形态像撒了盐与真实城市的团块状布局差异明显。原因是邻域影响权重太小或随机扰动系数过大导致远离城市中心的位置也有较高概率独立成城。解决方法是提高 neighborhood_weight 至 0.8 以上同时把 random_mag 压缩到 0.02 以下。此外可以在每次迭代后加形态学处理例如开运算去除孤立的单像元但要谨慎使用过度平滑会丢失真实的新城跳跃式增长特征。5.4 水域和保护区像元被模拟成城市现象是模拟结果中水域内部或生态保护区内出现城市像元。原因有两个可能一是初始状态矩阵中水域被赋值为 0没有被标记为 -1二是在状态更新时只判断了非城市没有排除不可开发区。解决方法是初始化时对水域像元强制赋 -1并在 iter 循环更新时把 excluded 区域的可能开发像元再次清除。% 强制水域和不可开发区不变为城市 newly_built newly_built ~excluded;这种问题通常不会制造大的统计偏差但会在图上非常显眼也最容易让结果受到质疑。建模第一件事就是把永久性非建设用地的标记检查一遍不要在事后找补。5.5 不同参数下模拟结果差异巨大现象是略微调整某个参数模拟结果就会产生非常大的变化。这说明模型对某些参数非常敏感在率定时如果参数空间覆盖不足很可能错过最优区间。解决方法是先做敏感性分析在参数率定前逐一改变参数值比如每次只变动一个参数观察 Kappa 值的变化幅度识别敏感参数。对敏感参数加密采样对不敏感参数固定为经验值这样可以大幅压缩率定计算量。敏感性分析本身就是空间模型研究的常见组成部分顺手还可以写进方法部分。5.6 矩阵过大导致 Matlab 内存溢出现象是模型跑一年就报内存不足或者速度极慢。原因是研究区大、分辨率高矩阵运算在每次 conv2 和 rand 操作中产生临时大数组叠加后内存消耗巨大。解决方法是分块计算邻域密度或者用稀疏矩阵表示城市像元。更实用的策略是先降采样到 100 米至 250 米做参数率定最终模拟再恢复到全分辨率。城市增长模拟的尺度效应非常明显不同分辨率下最优参数不同因此率定和最终模拟必须使用同一分辨率。如果条件允许可以把研究区按流域或行政区划切块分别建模、拼合结果。6. 进阶用法多情景模拟、空间自校验与政策落地价值参数率定完成后未来情景模拟才是这个模型的最终舞台。用历史数据校准的参数叠加未来不同的政策情景约束就可以回答“不同政策下城市会怎样长”这类规划问题。6.1 多情景模拟生态保护与城市开发强度的权衡我常用三套情景做对比。基准情景延续过去 10 年的增长趋势不加额外约束。生态优先情景把未划入保护区的林地、草地、水体也加入保护掩膜提高 threshold压低随机扰动模拟在严格生态约束下的增长格局。规划导向情景根据城市总体规划中的增长边界和重点发展片区在适宜性矩阵中加入“规划优势区”加分项让开发在引导区域内更容易发生。每一组模拟使用同一组率定后的基础参数只改变约束和偏好系数。三张结果图叠加显示可以直接得出不同情景下侵占生态用地的面积和空间位置这会成为规划方案对比里非常有说服力的材料。6.2 模拟结果的空间自校验除了 Kappa 系数还要看空间分异是否合理。我通常会把模拟结果按距离城市中心分环带统计城市像元比例再与真实数据对比。城市增长一般符合距离衰减规律距中心越远新增城市比例越小。如果模拟结果在远郊也出现大量新增可能说明空间自相关性没控制好。另一个做法是划定几个重点片区分别统计面积增长然后与真实历史增速对比找到局部高估或低估的原因。这比只看一个总体 Kappa 系数更能定位模型内部的系统误差。6.3 给从业者的实际建议我用这套方法完成某跨平台系统的演练时最后交付物的核心并不是那几张 2030 年城市预测图而是要素齐备的方法论数据说明、率定参数、模型结构、验证结果和不确定性说明。规划师和决策者关心的是你到底有多大把握最坏情况是什么建议后续把每一年的模拟结果全部保存在 GIS 里做年际变化动态演进同时配合 Markov 转移矩阵做数量层面的交叉验证。我一般在迁移到不同地区后会更新距离因子比如港口高速、高铁站位置把新的交通可达性显式加进适宜性选型里。坦白说这些环节里每个都有爬不完的坑。初学时我也曾在参数组合里不设置任何上限结果第一年就爆掉了也曾在数据投影不一致时反复折腾最后才发现是坐标系的问题。现在我的习惯是拿到数据先做五个检查——行列数、投影、像元大小、年份、类别编码然后才开始配置模型参数。希望这条流程能帮你少走一段弯路也希望你做完第一张未来城市模拟图后能体会到这套方法的潜力。本文还有配套的精品资源点击获取
