Matlab实现元胞自动机模拟城市增长

简介：基于元胞自动机的Matlab城市增长模拟项目，适合城乡规划、地理信息科学等方向的师生与研究人员。项目以艾哈迈达巴德地区为案例，通过CA模型模拟未开发土地向住宅、商业、工业等用地的演变，并依据人口迁移、交通便利度等规则迭代更新，覆盖状态转移规则、邻域滤波、道路距离影响等核心算法，并配有实际地图图像与测试图，便于作为课程设计或课题研究基线。压缩包共12个文件，核心为7个.m脚本，包括主程序Main_code.m及多项邻域处理辅助函数；另有3张jpg测试图像、1份txt运行指南和1张png结果示例图，包体仅105KB。已有264人学习下载。配合运行指南和脚本注释，可快速理解从图像到状态矩阵的转换、7×7邻域影响计算及结果可视化全流程，还可在现有规则基础上调整参数，拓展至其他城市增长预测场景。 直接在社区里看到这个标题时，我第一反应是“终于有人把CA模型拿出来聊了”。城市增长模拟这个方向，规划、地理、土木背景的读者应该都不陌生，而元胞自动机（Cellular Automata，简称CA）正是做这类时空模拟最经典的轻量级工具。用Matlab来做，最合适的点在于：你不需要像写C++或Python那样搭一套完整的数据结构，Matlab的矩阵运算天然匹配元胞网格，几十行核心代码就能跑出一个像样的模拟结果。这篇文章我把完整的建模思路、Matlab实现细节、参数标定方法和踩过的坑都整理出来，供正在做毕业设计、论文预实验或实际规划项目的人参考。

1. 项目思路拆解：元胞自动机为什么适合模拟城市增长

城市增长本质上是一个“局部变化驱动全局演化”的过程——一块地是否从非城市用地转为城市用地，往往取决于它周边的状态、交通可达性、自然条件约束和规划政策引导。这种“局部规则决定全局形态”的特性，正好撞在元胞自动机的枪口上。

1.1 元胞自动机的四个核心要素

CA模型听起来玄乎，其实拆开就四样东西：

• 元胞（Cell）：模拟区域划分成的网格单元。在城市增长模拟中，每个元胞通常代表一个固定大小的地块，比如30米乘30米、100米乘100米，具体尺寸取决于你的数据分辨率。
• 状态（State）：每个元胞的属性值。最简化的情况是二值状态——0表示非城市用地，1表示城市用地。也可以扩展成多状态，比如农村、城市、水域、植被等，但多状态会明显增加标定复杂度。
• 邻域（Neighborhood）：决定某个元胞状态变化时参考哪些“邻居”。最常用的是Moore邻域（周围3x3共8个格子），也有用扩展邻域（5x5、7x7）的。邻域越大，模拟出的城市斑块越连续，但计算量也成倍增加。
• 转换规则（Transition Rule）：这是全模型的核心，它决定了一个元胞在给定邻域状态和环境条件下，从非城市变成城市的概率。

城市增长CA的基本逻辑就这么一句话：每个时间步（通常对应现实中的1年），遍历所有非城市元胞，计算它转变为城市的概率，再和随机阈值比较，超过阈值就“城市化”。跑完N个时间步，就得到了N年后的模拟城市格局。

1.2 为什么选择Matlab而不是其他工具

我见过有人用ArcGIS的Model Builder做类似的事，也有用Python的NumPy写的版本，但Matlab在这个场景下有一个独特优势：代码表达与模型逻辑的对应关系非常直观。CA模型里的“全局网格状态更新”直接对应Matlab的矩阵运算，“邻域统计”可以用卷积函数（conv2）一行实现，做可视化更是点几个函数的事。相比之下，ArcGIS操作繁琐且不灵活，Python虽然也行但需要多装不少库。

另外，Matlab的调试体验对学术场景很友好。你可以在循环里中断查看任意元胞的状态变化过程，这对理解和调优CA参数非常有用。很多做城市模拟的论文本身就是用Matlab跑的实验，复现起来也方便。

2. 模型设计：从数据到转换规则

2.1 需要准备哪些数据

城市增长CA模型的输入数据并不复杂，但每一项都直接影响模型可信度。

在Matlab中，统一用geotiffread或imread读取这些栅格，然后利用经纬度或投影坐标系将它们配准到同一网格。配准这一步是搞遥感、GIS的人的老熟人——如果分辨率不一致，用imresize重采样到统一尺寸；如果投影不同，最好在ArcGIS或QGIS里先统一投影再导出。

2.2 转换概率如何计算

这是模型的“灵魂”。我采用了主流文献里常见的逻辑回归（Logistic Regression）形式：

P_growth = 1 / (1 + exp(-z))

z = a0 + a1 * dist_road + a2 * dist_cbd + a3 * density_neighbor + a4 * slope + ...

核心思路是：用历史数据（比如2010年和2020年两期土地利用）来拟合这个回归方程。把那些“从非城市变成城市”的元胞标记为1，“始终未变”的标记为0，把它们对应的距离、密度、坡度等变量提取出来做逻辑回归，得到的系数就是转换规则里的权重。

实际操作中，我建议至少包含三个变量：道路距离（反映交通引力）、城市核密度K（反映邻域城市化强度）、坡度（反映建设适宜性）。这三个变量已经能模拟出比较像样的城市扩张形态。想更精细就再加距CBD距离、到水体距离等。

城市核密度K的计算是这样：

K = sum(neighborhood_cells_urban) / total_neighbor_count

在Matlab里，可以直接用conv2对当前时刻的城市状态矩阵做卷积：

K = conv2(double(cityState), ones(3)/9, 'same');

这一行就完成了每个元胞3x3邻域内城市占比的计算，向量化效率极高，比写双层for循环快两个数量级，这也是Matlab实现CA的核心技巧。

2.3 随机性和约束条件

纯粹的逻辑回归概率会让城市增长显得“太确定”——历史上明明是随机扩散，模型却产出规则图案。所以我在每个时间步引入一个随机扰动项，用蒙特卡洛思想过滤转换概率：

cell_to_urban = (P_growth > rand(size(cityState))) & (developable == 1);

这里rand生成0到1之间的均匀随机数，与概率比较就等于以P_growth的概率发生转换。你每次运行的结果会有细微差别，这是正常的，CA模拟本身就是一种随机模拟方法。实验中要复现结果就固定随机种子。

developable是约束矩阵，坡度大于25度或属于水域的元胞就赋值为0，永远不允许开发。这一步放在概率计算之后，能够高效排除不可建设区域，也比在回归里强行加哑变量更可控。

3. Matlab完整实现过程

3.1 代码框架总览

整个模拟程序可以分为四个模块：数据预处理、参数标定、循环模拟、可视化输出。核心的模拟循环大概五六十行，一个上午基本能写完调试通过。

3.2 模拟主循环实现

下面是一段精简但完整的模拟循环代码，可以直接在Matlab里跑通：

% 参数设置 width = size(cityInit, 1); height = size(cityInit, 2); years = 20; % 模拟20年 developable = ones(width, height); developable(slope > 25 | water == 1) = 0; % 逻辑回归系数（示意，实际需用历史数据标定） coeff.road = -0.008; coeff.cbd = -0.003; coeff.density = 2.5; coeff.slope = -0.05; coeff.const = -1.2; % 随机种子（保证结果可复现） rng(42); % 初始化状态 cityState = cityInit; growthRecord = zeros(width, height, years); for t = 1:years % 1. 计算距离因子（每个时间步可更新，这里假设静态） % 距离栅格distRoad和distCBD在循环外已计算好 % 2. 计算邻域城市密度 K = conv2(double(cityState), ones(3)/9, 'same'); % 3. 计算综合发展概率 z = coeff.const ... + coeff.road * distRoad ... + coeff.cbd * distCBD ... + coeff.density * K ... + coeff.slope * slope; P = 1 ./ (1 + exp(-z)); % 4. 加约束和随机性 P = P .* developable; newUrban = (P > rand(width, height)) & ~cityState; % 5. 更新城市状态 cityState = cityState | newUrban; growthRecord(:, :, t) = double(newUrban); fprintf('第%d年模拟完成，新增城市元胞：%d\n', t, sum(newUrban(:))); end

3.3 关键点解释

代码里有几个细节值得展开说。

第一，conv2计算邻域密度时，ones(3)/9的窗口对应Moore邻域。如果你想让道路更近的城区扩张更快，可以把窗口改成带权重的核，比如给中心元胞更大的权重。但做模拟实验时不要随便改窗口，因为窗口形态本身就是一个需要论证的模型假设。

第二，~cityState这个条件保证了已经城市化的元胞不会重复经历“转变”过程。城市用地一旦变成城市，在模型里就是永久性的——这在大多数城市增长模拟中是合理的假设，毕竟实际中“去城市化”的案例非常罕见。

第三，约束条件P .* developable用一个矩阵相乘就实现了不可开发区域的屏蔽。如果某个区域的规划政策是限制开发但未完全禁止，可以把developable的值设置成0到1之间的系数，比如0.2，表示该区域开发概率打两折。这种设置在政策模拟中很好用——你可以在规划方案A和方案B之间来回切换比较扩张结果。

3.4 可视化与结果输出

模拟完成之后，我习惯做两张图：一张是最终城市格局，另一张是逐年新增城市密度变化图。

% 最终城市状态地图 figure; imagesc(cityState); colormap([0.85 0.85 0.85; 0.2 0.3 0.5]); axis equal; axis tight; set(gca, 'YDir', 'normal'); title('20年后模拟城市格局'); % 逐年新增面积变化 figure; annualGrowth = squeeze(sum(sum(growthRecord, 1), 2)); bar(1:years, annualGrowth); xlabel('年份'); ylabel('新增城市元胞数'); grid on;

如果需要出GeoTIFF文件放到ArcGIS里叠加分析，用geotiffwrite输出即可：

geotiffwrite('simulation_result_2040.tif', uint8(cityState), R_ref);

其中R_ref是参考栅格的空间参考对象，从原始数据里用geotiffread读出来时可以直接获取。

4. 参数标定与结果验证

4.1 历史数据标定法

直接拍脑袋定系数是做CA模拟的大忌。合理的做法是用两期历史数据来标定——比如你有2010年和2020年两期土地利用图，那么：

1. 从2010年城市状态出发模拟一年；
2. 对比2020年真实城市状态的差异，调整参数让模型预测结果最接近真实；
3. 用另外一期数据（比如2015年的）做验证。

逻辑回归的具体做法是：在2010年非城市元胞中随机抽样（因为栅格数量太大，全量跑逻辑回归很慢），记录每个样本元胞的2010年特征变量（距离道路、距离CBD、邻域密度、坡度）和2020年标签（是否变为城市），然后用fitglm拟合：

% sampleData每一行是特征变量，label是0/1响应变量 model = fitglm(sampleData, label, 'Distribution', 'binomial');

拟合完成后，coeff就来自model.Coefficients.Estimate。这一步既解决了“参数从哪来”的问题，也能通过显著性检验判断哪些变量应该保留、哪些应该舍弃。

我实测下来的经验是：道路距离和邻域密度几乎总是显著变量；坡度在平原城市作用弱，在山区城市作用极强；距CBD距离在单中心城市的模型中很关键，在多中心城市中通常不显著。这些结论不是算法问题，而是城市发展规律的体现，拿出来写论文就是很好的分析点。

4.2 精度评估

模型跑完后，用混淆矩阵来评估模拟结果和真实格局的吻合程度。常用指标包括：

• 总体精度（Overall Accuracy）：预测正确的元胞占比；
• Kappa系数：排除随机一致性的精度指标；
• 城市元胞命中率：真实新增城市元胞中有多少被模型正确预测。

Kappa计算比较繁琐，可以直接用Matlab自带的混淆矩阵工具，或者自己写一个简单的函数：

% simMap和realMap是相同尺寸的二值矩阵 confMat = confusionmat(realMap(:), simMap(:)); OA = trace(confMat) / sum(confMat(:));

如果OA在0.85以上，基本可以说模型表现良好。低于0.8就需要检查变量设置和数据配准问题了。

4.3 模型局限与理性看待

CA模型虽然是城市模拟的经典工具，但它本质上是“混合了随机性的经验模型”，缺乏对城市发展内在经济驱动力的解释——比如就业增长、人口迁移、产业集聚这些深层原因并不直接出现在模型里。CA模拟的结果更适合作为“趋势外推”的参考，不能说它是“预测未来”的绝对答案。在论文或项目报告中，这个定性表述要格外注意，否则容易引起争议。

5. 常见问题与排查技巧实录

5.1 模拟结果全是噪点，没有连片城市

这是最容易遇到的问题。如果模拟出的城市格局像撒胡椒面一样零散分布，最可能的原因是邻域密度项的系数太高或太低，或者邻域窗口太小。

解决办法：把邻域窗口从3x3扩展到5x5，或者增大density系数。在Matlab里把conv2的核改成ones(5)/25即可，实际效果立竿见影。另外，检查一下distRoad是不是没有做归一化——距离值如果动辄几千上万，而逻辑回归的常数项只有个位数，z值就会被距离项主导，概率完全由距离决定，形态自然不对。建议把所有连续变量都做min-max归一化到0到1。

5.2 城市增长慢得像蜗牛，20年才长了几个元胞

这种情况通常是因为逻辑回归的常数项过于负向，导致全区域的基准概率都非常低。检查一下coeff.const的量级，如果小于-5，那么即使所有正向变量都取最大值，概率也上不去。

另外，随机数种子的影响比你想象的大。如果你rng固定在一个不好的种子上，可能在某个时间步恰好所有概率比较都没通过，整个模拟就停滞了。建议跑3到5个随机种子，取平均结果或者选中间值代表情景。

5.3 边界元胞的邻域密度计算偏小

所有基于卷积的计算都有边界效应——图像边缘的元胞做3x3卷积时，卷积核有一部分超出图像范围，Matlab的conv2默认在边界处补零，导致边界区域的邻域密度被低估。如果你的研究区域是完整行政区划，而城市边缘又恰好在行政边界附近，这种偏差就会影响模拟结果。

处理方案有几种：一是用conv2(..., 'same')之后，直接忽略边界向外扩一圈得到的结果；二是给研究区域周围扩展一圈缓冲带，缓冲带内的状态从真实数据补全。我一般倾向第二种，因为城市规划模拟的研究区往往不到整个行政区的边界，留一圈缓冲区不仅修正卷积误差，也更符合城市发展连续性的实际。

5.4 模拟结果在相邻时间步之间反复震荡

这种“闪烁”现象通常是因为每一次状态更新用了异步更新的方式——也就是说，循环内部边更新城市状态、边用已经更新过的新值去计算后续元胞的概率。CA模型严格来说要求同步更新：当前时间步所有元胞的状态变化都基于上一时间步的全局状态，更新全部计算完之后，再一次性替换。

我在代码里用newUrban = (P > rand) & ~cityState然后cityState = cityState | newUrban，就是先算完整张概率矩阵，再统一更新状态矩阵，这就是标准同步更新。如果你发现模拟结果有震荡，先检查代码里是否无意中在同一个时间步内修改了cityState却又拿来参与后面的计算。

6. 一点拓展建议

CA模型的框架搭好之后，可以扩展的方向很多。

一个值得尝试的方向是“多情景模拟”。设置高位增长、基准增长、紧凑发展三种情景，分别调整约束矩阵和系数，就可以对比不同开发策略下的城市格局差异，这是规划决策中非常常见的应用方式。实际操作时，我一般把约束矩阵做成一个可交互变量，用Matlab的live script或者直接写个inputdlg对话窗口，切换起来很方便。

另一个方向是引入更细粒度的分区控制，比如工业区、居住区、商业区各自用不同的转换规则。这样模型就从单状态提升到了多状态，虽然标定工作量翻倍，但模拟结果的解释力会明显增强。

还有一个比较有意思的点：把CA模型和未来土地利用需求总量约束结合。比如你已知规划目标到2040年城市总建设用地需要增加多少平方公里，这在CA框架内可以换算成需要新增多少个元胞，然后在时间步循环中加一个全局判断条件——达到这个总量就停止新增。这种方法能保证模拟的城市总面积符合宏观规划，比纯粹的规则驱动更有说服力。我已经把相关代码放在自己的项目仓库里，标题提到的这个Matlab实现也增加了这个功能模块，感兴趣可以直接测试一下。

最后再分享一个小经验：做这类模型，数据不追求多，但配准必须严格。曾经有一次我用了两份分辨率不一致的栅格直接跑模型，结果出来的城市形态全挤在投影错位的边界上，排查了很久才发现是数据问题，白白浪费了几天时间。建议所有输入栅格在进入模拟前，统一用imresize重采样到同一尺寸，并用sum(water + slope > 0)这类矩阵运算快速检查配准情况。数据对了，模型就跑顺了。

本文还有配套的精品资源，点击获取