1. 为什么我选择用遥感技术来盯一座城市的土地利用变化土地利用变化这件事听起来像是规划部门才关心的宏观议题但真落到实操层面它跟每一个做地理信息、做城市研究、甚至做房地产数据分析的人都脱不开关系。我最早接触这个方向是因为一个模拟项目需要评估某城市过去十年的建成区扩张速度当时手头没有现成的权威统计年鉴即便有口径也对不上。后来发现遥感影像才是那个最诚实的数据源——它不撒谎地表上每一块水泥、每一片农田、每一处水体都在像素里留着痕迹。所谓“基于RS的沈阳土地利用情况”说白了就是拿遥感影像当底图通过分类算法把沈阳这座城市的地表覆盖拆解成耕地、林地、草地、水域、建设用地和未利用地这几大类然后看它们在时间轴上怎么变、往哪个方向变、变的速度是快还是慢。这件事能解决的问题很直接城市到底在往哪个方向摊大饼耕地被占了多少生态用地有没有被蚕食这些问题的答案光靠文字报告是看不出来的必须落到空间上、落到像素上。适合看这篇内容的人我大致归了三类。第一类是地理信息、遥感、城乡规划相关专业的学生或从业者需要一套能跑通的完整流程第二类是做城市数据分析、想引入空间维度的朋友可能之前只玩过表格数据现在想试试栅格数据第三类是对城市发展本身感兴趣、想用技术手段自己验证一些判断的人。不管你是哪一类只要跟着走一遍至少能搞清楚遥感分类这件事的完整链路长什么样。我自己的体会是遥感土地利用分类这个活儿门槛不在算法本身而在“数据预处理”和“分类体系设计”这两个环节。很多人一上来就想着上深度学习、上Transformer结果连影像的辐射定标和大气校正都没做分类结果自然惨不忍睹。所以这篇内容我会把重心放在那些真正决定成败的细节上而不是堆砌模型名称。2. 数据准备与预处理地基打不牢后面全白搭2.1 影像源的选择逻辑与参数权衡做沈阳的土地利用分类第一步是选影像。Landsat系列是绕不开的选项原因很简单免费、时间序列长、空间分辨率30米虽然不算高但对于城市尺度的土地利用分类足够用。沈阳建成区面积大约在500到600平方公里量级30米分辨率意味着一个像素对应900平方米一个中等规模的住宅小区大概占几十个像素识别建设用地绰绰有余。如果要做更精细的城区内部区分比如把建设用地再拆成居住、工业、商业那就得考虑Sentinel-210米分辨率能提供更多纹理信息。但代价是数据量翻了好几倍处理时间也上去了。我一般建议如果只是做一级分类六大类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_estimators100到200max_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 常见问题速查表问题现象可能原因排查方法解决方案分类结果大片错乱大气校正未做或做错检查影像反射率值范围是否合理重新做大气校正建设用地被误判为耕地训练样本中建设用地样本不足查看混淆矩阵中建设用地的生产者精度增加建设用地训练样本水体边界模糊NDWI阈值设置不当可视化NDWI影像检查水体边界调整NDWI阈值或引入SWIR波段分类结果有规律条带影像条带未去除目视检查原始影像去条带处理Kappa系数低于0.7训练样本质量差或类别混淆严重检查混淆矩阵找出混淆最严重的类别对增加样本量或引入纹理特征面积统计与统计年鉴差异大分类体系不一致或精度不足对比分类体系定义调整分类体系或提升分类精度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, dtypeint) # 训练随机森林 rf RandomForestClassifier(n_estimators150, max_depth25, random_state42) 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(count1, dtyperasterio.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(fKappa系数{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一定要做交叉定标否则分类结果的差异可能来自传感器本身而不是地表变化。这个坑我踩过后来花了整整一周时间重新处理数据教训深刻。 
长河君 · 2026/10/15 13:41:03 · 阅读 1.2k 基于遥感技术的沈阳土地利用变化分类与城市扩张分析实战
1. 为什么我选择用遥感技术来盯一座城市的土地利用变化土地利用变化这件事听起来像是规划部门才关心的宏观议题但真落到实操层面它跟每一个做地理信息、做城市研究、甚至做房地产数据分析的人都脱不开关系。我最早接触这个方向是因为一个模拟项目需要评估某城市过去十年的建成区扩张速度当时手头没有现成的权威统计年鉴即便有口径也对不上。后来发现遥感影像才是那个最诚实的数据源——它不撒谎地表上每一块水泥、每一片农田、每一处水体都在像素里留着痕迹。所谓“基于RS的沈阳土地利用情况”说白了就是拿遥感影像当底图通过分类算法把沈阳这座城市的地表覆盖拆解成耕地、林地、草地、水域、建设用地和未利用地这几大类然后看它们在时间轴上怎么变、往哪个方向变、变的速度是快还是慢。这件事能解决的问题很直接城市到底在往哪个方向摊大饼耕地被占了多少生态用地有没有被蚕食这些问题的答案光靠文字报告是看不出来的必须落到空间上、落到像素上。适合看这篇内容的人我大致归了三类。第一类是地理信息、遥感、城乡规划相关专业的学生或从业者需要一套能跑通的完整流程第二类是做城市数据分析、想引入空间维度的朋友可能之前只玩过表格数据现在想试试栅格数据第三类是对城市发展本身感兴趣、想用技术手段自己验证一些判断的人。不管你是哪一类只要跟着走一遍至少能搞清楚遥感分类这件事的完整链路长什么样。我自己的体会是遥感土地利用分类这个活儿门槛不在算法本身而在“数据预处理”和“分类体系设计”这两个环节。很多人一上来就想着上深度学习、上Transformer结果连影像的辐射定标和大气校正都没做分类结果自然惨不忍睹。所以这篇内容我会把重心放在那些真正决定成败的细节上而不是堆砌模型名称。2. 数据准备与预处理地基打不牢后面全白搭2.1 影像源的选择逻辑与参数权衡做沈阳的土地利用分类第一步是选影像。Landsat系列是绕不开的选项原因很简单免费、时间序列长、空间分辨率30米虽然不算高但对于城市尺度的土地利用分类足够用。沈阳建成区面积大约在500到600平方公里量级30米分辨率意味着一个像素对应900平方米一个中等规模的住宅小区大概占几十个像素识别建设用地绰绰有余。如果要做更精细的城区内部区分比如把建设用地再拆成居住、工业、商业那就得考虑Sentinel-210米分辨率能提供更多纹理信息。但代价是数据量翻了好几倍处理时间也上去了。我一般建议如果只是做一级分类六大类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_estimators100到200max_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 常见问题速查表问题现象可能原因排查方法解决方案分类结果大片错乱大气校正未做或做错检查影像反射率值范围是否合理重新做大气校正建设用地被误判为耕地训练样本中建设用地样本不足查看混淆矩阵中建设用地的生产者精度增加建设用地训练样本水体边界模糊NDWI阈值设置不当可视化NDWI影像检查水体边界调整NDWI阈值或引入SWIR波段分类结果有规律条带影像条带未去除目视检查原始影像去条带处理Kappa系数低于0.7训练样本质量差或类别混淆严重检查混淆矩阵找出混淆最严重的类别对增加样本量或引入纹理特征面积统计与统计年鉴差异大分类体系不一致或精度不足对比分类体系定义调整分类体系或提升分类精度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, dtypeint) # 训练随机森林 rf RandomForestClassifier(n_estimators150, max_depth25, random_state42) 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(count1, dtyperasterio.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(fKappa系数{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一定要做交叉定标否则分类结果的差异可能来自传感器本身而不是地表变化。这个坑我踩过后来花了整整一周时间重新处理数据教训深刻。