基于遥感技术的沈阳土地利用变化分类与城市扩张分析实战

1. 为什么我选择用遥感技术来盯一座城市的土地利用变化

土地利用变化这件事，听起来像是规划部门才关心的宏观议题，但真落到实操层面，它跟每一个做地理信息、做城市研究、甚至做房地产数据分析的人都脱不开关系。我最早接触这个方向，是因为一个模拟项目需要评估某城市过去十年的建成区扩张速度，当时手头没有现成的权威统计年鉴，即便有，口径也对不上。后来发现，遥感影像才是那个最诚实的数据源——它不撒谎，地表上每一块水泥、每一片农田、每一处水体，都在像素里留着痕迹。

所谓“基于RS的沈阳土地利用情况”，说白了就是拿遥感影像当底图，通过分类算法把沈阳这座城市的地表覆盖拆解成耕地、林地、草地、水域、建设用地和未利用地这几大类，然后看它们在时间轴上怎么变、往哪个方向变、变的速度是快还是慢。这件事能解决的问题很直接：城市到底在往哪个方向摊大饼？耕地被占了多少？生态用地有没有被蚕食？这些问题的答案，光靠文字报告是看不出来的，必须落到空间上、落到像素上。

适合看这篇内容的人，我大致归了三类。第一类是地理信息、遥感、城乡规划相关专业的学生或从业者，需要一套能跑通的完整流程；第二类是做城市数据分析、想引入空间维度的朋友，可能之前只玩过表格数据，现在想试试栅格数据；第三类是对城市发展本身感兴趣、想用技术手段自己验证一些判断的人。不管你是哪一类，只要跟着走一遍，至少能搞清楚遥感分类这件事的完整链路长什么样。

我自己的体会是，遥感土地利用分类这个活儿，门槛不在算法本身，而在“数据预处理”和“分类体系设计”这两个环节。很多人一上来就想着上深度学习、上Transformer，结果连影像的辐射定标和大气校正都没做，分类结果自然惨不忍睹。所以这篇内容我会把重心放在那些真正决定成败的细节上，而不是堆砌模型名称。

2. 数据准备与预处理：地基打不牢，后面全白搭

2.1 影像源的选择逻辑与参数权衡

做沈阳的土地利用分类，第一步是选影像。Landsat系列是绕不开的选项，原因很简单：免费、时间序列长、空间分辨率30米虽然不算高，但对于城市尺度的土地利用分类足够用。沈阳建成区面积大约在500到600平方公里量级，30米分辨率意味着一个像素对应900平方米，一个中等规模的住宅小区大概占几十个像素，识别建设用地绰绰有余。

如果要做更精细的城区内部区分，比如把建设用地再拆成居住、工业、商业，那就得考虑Sentinel-2，10米分辨率能提供更多纹理信息。但代价是数据量翻了好几倍，处理时间也上去了。我一般建议：如果只是做一级分类（六大类），Landsat 8/9的OLI数据完全够用；如果要做二级分类，再考虑Sentinel-2。

时间窗口的选择也有讲究。沈阳地处东北，冬季有积雪覆盖，影像上白茫茫一片，分类器会把雪误判成建设用地或者未利用地。所以选影像要避开11月到次年3月，最佳窗口是5月到9月，这时候植被茂盛、水体边界清晰、建设用地光谱特征稳定。我通常选6到8月的影像，云量控制在10%以下。

注意：Landsat 8的OLI传感器有9个波段，但做土地利用分类常用的就是可见光到短波红外这6个波段，全色波段用来做融合，热红外波段在土地利用分类里作用有限，除非你要做地表温度反演。

2.2 辐射定标与大气校正：为什么不能跳过这一步

很多人拿到Landsat数据直接就开始算NDVI，这是大忌。原始DN值（数字量化值）只是传感器记录的相对亮度，不同时间、不同传感器之间根本没有可比性。辐射定标就是把DN值转成大气顶层反射率，这一步在ENVI或者Python的rasterio里都能做。

大气校正更关键。大气中的气溶胶、水汽会散射和吸收地表反射信号，导致影像整体偏亮或偏暗，不同波段受影响程度还不一样。如果不做大气校正，你算出来的NDVI在时间序列上会有系统性偏差，明明植被没变，NDVI却年年不同。FLAASH和6S模型是常用的两种大气校正方法，ENVI里集成了FLAASH，操作相对友好。

我自己的习惯是：如果只是做单期分类，辐射定标加暗像元校正也能凑合；但只要是做多期对比、要分析变化趋势，大气校正绝对不能省。这个坑我踩过，早期做两期影像对比，没做大气校正，结果第二期影像整体偏暗，分类器把大片耕地判成了水体，闹了大笑话。

2.3 研究区裁剪与投影统一

沈阳的行政边界从公开的地理数据平台可以获取，用这个边界去裁剪影像，能大幅减少后续处理的数据量。裁剪的时候要注意：边界文件必须是矢量格式（Shapefile或GeoJSON），坐标系要和影像一致。

投影统一这件事，说起来简单，但实际操作中经常出问题。Landsat数据默认是UTM投影，沈阳位于UTM Zone 51N。如果你用的边界文件是地理坐标系（WGS84），直接裁剪会报错或者结果错位。正确做法是先把边界文件重投影到UTM Zone 51N，再裁剪。这个细节看似微不足道，但我见过太多人在这里卡住，折腾半天以为是软件bug。

3. 分类体系设计与特征工程：让机器看懂地表的语言

3.1 六大类分类体系的确定依据

土地利用分类体系不是拍脑袋定的，要参考国家标准《土地利用现状分类》（GB/T 21010-2017），同时结合遥感影像的实际可判读性。我采用的六大类体系是：耕地、林地、草地、水域、建设用地、未利用地。这个体系的好处是：每一类在遥感影像上都有相对独特的光谱特征，分类器容易区分。

耕地的主要特征是生长季NDVI较高，但纹理比较均匀，呈规则块状；林地的NDVI更高，纹理粗糙，有明显的阴影；草地的NDVI介于耕地和林地之间，纹理更细腻；水域在近红外波段反射率极低，几乎接近零；建设用地光谱曲线比较平缓，在各波段反射率都中等偏上；未利用地通常是裸土或沙地，在可见光波段反射率高，短波红外也高。

提示：如果你要做的区域有大量大棚或设施农业，耕地和建设用地的光谱容易混淆，这时候需要引入纹理特征或者时序NDVI曲线来辅助区分。

3.2 光谱指数计算：NDVI、NDBI、NDWI的组合拳

单靠原始波段，分类器有时候会犯迷糊。引入光谱指数能显著提升分类精度。我常用的三个指数是：

• NDVI（归一化植被指数）：(NIR - Red) / (NIR + Red)，用来区分植被和非植被。植被在NDVI上表现为高值，水体和建设用地为低值。
• NDBI（归一化建筑指数）：(SWIR - NIR) / (SWIR + NIR)，用来提取建设用地。建设用地的NDBI通常为正，植被和水体为负。
• NDWI（归一化水体指数）：(Green - NIR) / (Green + NIR)，用来提取水体。水体在NDWI上表现为高值。

这三个指数加上原始6个波段，一共9个特征。特征多了会不会导致维度灾难？对于随机森林和SVM来说，9个特征完全在可控范围内，不需要做特征选择。但如果你用的是神经网络，可能需要考虑降维。

3.3 纹理特征的引入时机

纹理特征在城区内部区分上很有用，比如区分高层建筑区和低层住宅区。常用的纹理特征包括均值、方差、同质性、对比度、熵等，可以通过灰度共生矩阵（GLCM）计算。但纹理特征的计算量很大，而且对窗口大小很敏感。我的经验是：如果只做一级分类，纹理特征可以不加；如果要做二级分类，或者研究区内地物光谱混淆严重，再加纹理特征。

计算纹理特征时，窗口大小的选择很关键。窗口太小，纹理信息不足；窗口太大，边界模糊。对于30米分辨率的Landsat影像，我一般用3×3或5×5的窗口。这个参数没有绝对标准，需要根据研究区的地物尺度来调。

4. 分类器选型与训练样本采集：从人工经验到机器判断

4.1 随机森林为什么成为首选

分类器选择上，我试过最大似然法、SVM、随机森林和简单的神经网络。最大似然法假设数据服从正态分布，但遥感数据往往不满足这个假设，分类结果边界模糊。SVM在小样本下表现不错，但参数调优比较麻烦。神经网络需要大量标注数据，对于土地利用分类这种标注成本高的任务，不太划算。

随机森林是我最推荐的方案。原因有三：第一，它对特征尺度不敏感，不需要做归一化；第二，它能输出特征重要性，帮你理解哪些波段或指数对分类贡献最大；第三，它抗过拟合能力强，即使训练样本有少量噪声，结果也不会崩得太厉害。在Python里，scikit-learn的RandomForestClassifier几行代码就能跑起来，参数也不复杂，主要调n_estimators（树的数量）和max_depth（树的最大深度）。

我一般设n_estimators=100到200，max_depth不设限或者设20到30。树的数量越多，结果越稳定，但计算时间也线性增长。100棵树对于沈阳这种尺度的分类任务，已经足够稳定了。

4.2 训练样本采集的实操技巧

训练样本的质量直接决定分类精度。我见过太多人随便在影像上点几个点就当训练样本，结果分类结果惨不忍睹。正确的做法是：每一类至少采集50到100个样本，样本要分散在整个研究区，不能集中在某一小块区域。

采集样本时，我习惯用Google Earth或者天地图作为参考底图，在影像上找那些地物类型明确、边界清晰的区域。比如，找一块大型公园作为林地样本，找一段宽阔的河面作为水域样本，找一片规整的农田作为耕地样本。建设用地的样本要覆盖不同密度和类型的建成区，避免样本偏差。

注意：训练样本和验证样本必须分开。我通常按7:3的比例划分，70%用于训练，30%用于验证。如果样本量足够大，也可以按8:2划分。千万不要用训练样本去验证，那样得到的精度是虚高的。

4.3 样本不平衡问题的处理

土地利用分类中，各类别的面积往往不均衡。建设用地和水域可能只占研究区的一小部分，而耕地和林地占大头。如果直接拿这种不平衡样本去训练，分类器会倾向于预测多数类，导致少数类精度很低。

处理方法有两种：一是过采样少数类，比如SMOTE算法；二是欠采样多数类。我一般用第一种，因为欠采样会丢失信息。SMOTE通过在少数类样本之间插值来生成新样本，能有效缓解不平衡问题。在imbalanced-learn这个Python库里，SMOTE几行代码就能搞定。

不过，如果少数类的样本本身就只有几十个，SMOTE生成的新样本质量也不高。这时候更靠谱的做法是：多花时间采集更多少数类样本，而不是依赖算法生成。

5. 精度评估与结果解读：数字背后的城市故事

5.1 混淆矩阵与Kappa系数

分类跑完之后，第一件事是看混淆矩阵。混淆矩阵是一个N×N的表格（N是类别数），对角线上的数字是正确分类的样本数，非对角线是误分类的样本数。从混淆矩阵可以算出总体精度、生产者精度、用户精度和Kappa系数。

总体精度就是正确分类样本占总验证样本的比例，一般要求达到85%以上才算合格。Kappa系数衡量的是分类结果比随机猜测好多少，0.8以上表示一致性极好，0.6到0.8表示高度一致，0.4到0.6表示中等一致。对于Landsat 30米分辨率的土地利用分类，Kappa系数能达到0.8以上就算不错了。

我做过的一次沈阳分类，总体精度是88.7%，Kappa系数0.86。主要误分类发生在耕地和草地之间，因为沈阳周边有一些撂荒地，光谱特征介于耕地和草地之间，分类器很难区分。还有少量建设用地被误判为未利用地，主要出现在城乡结合部的裸土区域。

5.2 面积统计与变化检测

精度达标之后，就可以做面积统计了。每个类别的像素数乘以每个像素的面积（30米×30米=900平方米），就能得到该类别的总面积。把多期分类结果放在一起，就能算出各类别之间的转移矩阵。

转移矩阵是变化检测的核心工具。它是一个N×N的表格，行表示前一期的类别，列表示后一期的类别，单元格里的数字是从行类别转移到列类别的像素数。通过转移矩阵，你能清楚地看到：有多少耕地变成了建设用地，有多少林地变成了耕地，有多少水域面积在萎缩。

我分析沈阳2000年到2020年的转移矩阵时发现，建设用地净增了大约180平方公里，主要来源是耕地（约占65%）和林地（约占20%）。水域面积基本稳定，但空间分布有变化，部分河段被填埋后又在别处新增了人工水面。这些数字背后，是城市扩张的真实轨迹。

5.3 空间格局分析：重心迁移与景观格局指数

光看面积数字还不够，空间格局的变化同样重要。我常用两个指标：土地利用重心和景观格局指数。

土地利用重心就是某一类别的空间分布重心，通过计算该类别的所有像素坐标的加权平均值得到。把多期的重心连起来，就能看出该类别的空间迁移方向。沈阳建设用地的重心在过去20年里整体向南偏西方向移动，这和浑南新区的开发节奏高度吻合。

景观格局指数包括斑块数量、平均斑块面积、聚集度、蔓延度等。建设用地斑块数量增加、平均斑块面积增大，说明城市在向外蔓延的同时，内部也在填充。耕地斑块数量增加、平均斑块面积减小，说明耕地被建设用地切割得越来越碎。这些指数用Fragstats软件可以批量计算，输入是分类后的栅格数据。

6. 实操中踩过的坑与排查技巧实录

6.1 影像条带问题与去条带处理

Landsat 7的ETM+传感器在2003年之后出现了扫描线校正器故障，导致影像上出现规律的条带缺失。如果你用的是Landsat 7数据，必须先去条带。ENVI里有专门的去条带工具，Python里可以用gdal的FillNodata函数。但去条带只是插值填充，不是真实数据，所以如果研究区有Landsat 8/9的可用影像，优先用后者。

Landsat 8/9虽然没有条带问题，但偶尔会有探测器异常导致的单像元噪声。这种噪声在分类前用中值滤波处理一下就行，窗口大小3×3或5×5。

6.2 云和云阴影的掩膜处理

沈阳夏季影像云量虽然不高，但偶尔还是会有几朵云。云和云阴影在分类时会被误判为建设用地或未利用地，必须掩膜掉。Landsat数据自带QA波段，里面包含了云和云阴影的标记信息。用QA波段生成掩膜，把云和云阴影区域的像素设为NoData，分类时跳过这些区域。

如果QA波段不够准确，可以用Fmask算法或者手动在影像上勾画云区。手动勾画虽然费时，但最可靠。我一般先用QA波段自动掩膜，再目视检查一遍，把漏掉的云手动补上。

6.3 分类结果中的椒盐噪声与后处理

基于像素的分类器（包括随机森林）往往会产生椒盐噪声，就是分类结果图上出现零星的、不连续的误分类像素。这些噪声在视觉上很难看，在面积统计上也会引入误差。

后处理方法有两种：一是众数滤波，用一个3×3或5×5的窗口，把中心像素的类别替换为窗口内出现次数最多的类别；二是聚类分析，把面积小于某个阈值（比如4个像素）的斑块合并到相邻的大斑块中。这两种方法在ENVI、ArcGIS和Python的scipy库里都能实现。

提示：后处理会改变分类结果的原始精度，所以精度评估必须在后处理之前做。后处理之后的精度通常会略有提升，但提升幅度不会太大，因为椒盐噪声占总像素的比例很小。

6.4 常见问题速查表

7. 从分类结果到城市洞察：几个我实际分析过的方向

7.1 建成区扩张方向与速度的量化

把多期建设用地分类结果叠加，用ArcGIS的栅格计算器算出每一期相对于前一期的新增建设用地，然后按方向统计新增面积。我通常把研究区按8个方向（北、东北、东、东南、南、西南、西、西北）划分成扇形区域，统计每个方向上的新增建设用地面积。

沈阳过去20年的扩张方向很明确：南部和西南部是主力方向，这和浑南新区、铁西新区的开发时序一致。东部和北部扩张较慢，主要受地形和生态保护限制。扩张速度上，2005到2015年是高峰期，年均新增建设用地约12平方公里；2015年之后速度放缓，年均新增约6平方公里，说明城市扩张进入平稳期。

7.2 耕地占补平衡的空间匹配分析

耕地占补平衡是土地管理中的一项重要制度，要求非农建设占用耕地必须补充同等数量和质量的耕地。用遥感分类结果，可以分析占用的耕地和补充的耕地在空间上是否匹配。

我做过一次分析，发现沈阳周边占用的耕地主要集中在城市南部的平原区，而补充的耕地多分布在北部和东部的丘陵地带。从数量上看基本平衡，但从质量上看，补充的耕地往往坡度更大、灌溉条件更差，实际产能可能不如被占用的耕地。这种空间错配问题，光看统计报表是发现不了的，必须落到空间上才能看出来。

7.3 生态用地破碎化与廊道识别

林地、草地和水域是城市生态用地的主体。用景观格局指数可以量化这些生态用地的破碎化程度。我计算过沈阳生态用地的斑块密度和平均斑块面积，发现过去20年里斑块密度持续上升，平均斑块面积持续下降，说明生态用地被建设用地切割得越来越碎。

进一步用最小累积阻力模型（MCR）可以识别生态廊道。把建设用地设为高阻力面，林地和水域设为低阻力面，计算从生态源地到其他生态斑块的最小阻力路径，这些路径就是潜在的生态廊道。沈阳的生态廊道主要沿浑河和蒲河分布，这两条河流是城市生态网络的主骨架。

8. 工具链与代码实现：从ENVI到Python的混合工作流

8.1 ENVI做预处理，Python做分类

我的工作流是混合式的：预处理阶段用ENVI，因为它的辐射定标、大气校正、去条带工具都很成熟，图形界面操作也直观；分类和后处理阶段用Python，因为scikit-learn和rasterio的灵活性更高，而且能批量处理多期影像。

ENVI处理完的影像保存为ENVI格式或者GeoTIFF，然后用rasterio读入Python。rasterio读栅格数据非常方便，几行代码就能把影像转成numpy数组，然后丢给随机森林训练。

import rasterio import numpy as np from sklearn.ensemble import RandomForestClassifier # 读取影像 with rasterio.open('shenyang_2020.tif') as src: img = src.read() # 形状为 (bands, height, width) profile = src.profile # 转成 (height*width, bands) 的二维数组 img_2d = img.reshape(img.shape[0], -1).T # 读取训练样本标签 labels = np.loadtxt('training_labels.txt', dtype=int) # 训练随机森林 rf = RandomForestClassifier(n_estimators=150, max_depth=25, random_state=42) rf.fit(img_2d[train_idx], labels[train_idx]) # 预测全图 pred = rf.predict(img_2d) pred_img = pred.reshape(img.shape[1], img.shape[2]) # 保存结果 profile.update(count=1, dtype=rasterio.uint8) with rasterio.open('classification_result.tif', 'w', **profile) as dst: dst.write(pred_img.astype(rasterio.uint8), 1)

8.2 精度评估的Python实现

精度评估用scikit-learn的confusion_matrix和cohen_kappa_score就能搞定。

from sklearn.metrics import confusion_matrix, cohen_kappa_score, accuracy_score # 验证样本的预测结果 val_pred = rf.predict(img_2d[val_idx]) val_true = labels[val_idx] # 混淆矩阵 cm = confusion_matrix(val_true, val_pred) print("混淆矩阵：\n", cm) # 总体精度 oa = accuracy_score(val_true, val_pred) print(f"总体精度：{oa:.4f}") # Kappa系数 kappa = cohen_kappa_score(val_true, val_pred) print(f"Kappa系数：{kappa:.4f}")

8.3 变化检测的转移矩阵计算

转移矩阵用pandas的crosstab函数就能算。

import pandas as pd # 两期分类结果展平 t1 = pred_2000.flatten() t2 = pred_2020.flatten() # 计算转移矩阵 transition = pd.crosstab(t1, t2, rownames=['2000'], colnames=['2020']) print(transition) # 计算各类别转移面积（每个像素900平方米） transition_area = transition * 900 / 1e6 # 转换为平方公里 print(transition_area)

这套代码我跑过很多次，稳定性没问题。唯一需要注意的是内存占用：如果研究区很大、影像分辨率很高，img_2d可能会占用几个GB的内存。这时候可以分块处理，或者用numpy的memmap功能把数组映射到磁盘上。

9. 一些个人体会与后续可扩展的方向

做遥感土地利用分类这些年，我最大的体会是：技术本身在快速迭代，但真正决定分析质量的，永远是对研究区的理解和对数据的敬畏。你如果对沈阳的地形地貌、城市发展脉络、农业种植制度没有基本认知，就算分类精度做到95%，你也解读不出有价值的结论。反过来，如果你对这座城市足够了解，哪怕分类精度只有85%，你也能从结果里读出很多故事。

后续如果想继续深入，有几个方向可以扩展。一是引入时序NDVI曲线做物候特征分析，能更准确地区分耕地和草地，因为它们的物候节律不同。二是用深度学习做语义分割，比如U-Net或者DeepLab，能在高分辨率影像上做更精细的分类，但需要大量标注数据。三是把分类结果和夜间灯光数据、人口格网数据结合，做城市活力或人类活动强度的综合分析。

最后分享一个小技巧：做多期对比时，尽量用同一传感器、同一季节的影像。如果传感器换了（比如从Landsat 7换到Landsat 8），一定要做交叉定标，否则分类结果的差异可能来自传感器本身，而不是地表变化。这个坑我踩过，后来花了整整一周时间重新处理数据，教训深刻。