简介：面向城市规划与地理信息建模场景，基于元胞自动机的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 = ((a+b)*(a+c) + (c+d)*(b+d)) / 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 转移矩阵做数量层面的交叉验证。我一般在迁移到不同地区后，会更新距离因子（比如港口高速、高铁站位置），把新的交通可达性显式加进适宜性选型里。
坦白说，这些环节里每个都有爬不完的坑。初学时我也曾在参数组合里不设置任何上限，结果第一年就爆掉了；也曾在数据投影不一致时反复折腾，最后才发现是坐标系的问题。现在我的习惯是：拿到数据先做五个检查——行列数、投影、像元大小、年份、类别编码，然后才开始配置模型参数。希望这条流程能帮你少走一段弯路，也希望你做完第一张未来城市模拟图后，能体会到这套方法的潜力。
本文还有配套的精品资源，点击获取
