遥感影像滑坡场景分类:特征工程到SVM全流程Python实现
简介这是一套基于Python实现的遥感影像滑坡场景分类项目代码与文档面向毕业设计、课程设计及项目开发人员适合具备一定机器学习基础的读者参考。项目采用SVM分类器完成滑坡场景识别完整覆盖特征提取光谱特征、GLCM纹理统计、K-Means视觉词聚类、LDA主题抽象以及Libsvm分类四个阶段。压缩包内含21个文件包括4个Python脚本、4个pkl数据文件、4个txt样本文件、模型文件、LDA状态文件及项目说明文档等整体仅1.86MB结构紧凑便于快速部署。目前已有57人学习下载。代码经过严格测试附带详细md说明可从数据生成、样本标签制作到训练与预测全流程进行复用适合在此基础上做算法对比或扩展应用。1. 一个能直接跑通流程的滑坡场景分类骨架这套代码到底帮你省了什么山体滑坡的遥感影像识别最难的不是最后那个分类器而是中间那一大段特征工程怎么从一张影像里把光谱、纹理信息量化出来怎么把高维特征抽象成分类器能吃的输入再怎么组织训练样本让结果不翻车。这套基于 Python 的滑坡场景分类任务代码把整条管线完整串了一遍——光谱特征提取、GLCM 纹理统计、K-Means 视觉词袋、LDA 主题抽象、Libsvm 训练与预测每个环节都有对应的脚本和中间产物还有一份项目文档把流程和参数交代清楚。对正在做毕业设计或课程设计、需要快速搭起一个可运行的场景分类任务的人来说这套代码相当于把最耗时的特征工程和样本组织部分替你趟平了。这次任务做到的是场景级分类判断影像块是不是滑坡场景不是逐像素的语义分割先把这个边界搞清楚再往下看。2. 数据组织从哪下手训练样本怎么切、标签怎么对齐、文件列表怎么读拿到项目先别急着跑 classify.py先把数据管线理顺。遥感影像分类和普通图像分类最大的区别在于原始影像通常尺寸大、地物杂不能整张喂给分类器必须先切成带标签的样本块。这套项目里 data_tool.py 和 generate_samples_and_labels.py 就是干这件事的理解了这两个脚本后面所有环节都不会懵。2.1 滑窗切割样本块的尺寸与步长是第一个需要决策的参数遥感影像的处理一般是先人工标注出滑坡区域再把影像切成固定大小的样本块。常见做法是用一个固定尺寸的滑窗在影像上移动窗口落进标注区域的样本标记为正类落在非滑坡区域的标记为负类。项目里 data_tool.py 的核心职责就是这类批量切割和整理工作。import numpy as np import cv2 import os def crop_image_with_labels(image_path, label_path, window_size, stride, save_dir): 对遥感影像做滑窗切割并根据标注文件生成对应的样本标签 image_path: 原始影像路径大幅遥感影像 label_path: 标注文件路径二值掩膜1表示滑坡区域 window_size: 窗口大小比如(64, 64) stride: 滑动步长控制相邻窗口的重叠程度 save_dir: 保存切割结果的目录 img cv2.imread(image_path) label_map cv2.imread(label_path, cv2.IMREAD_GRAYSCALE) h, w img.shape[:2] samples, labels [], [] for y in range(0, h - window_size[1] 1, stride): for x in range(0, w - window_size[0] 1, stride): img_patch img[y:y window_size[1], x:x window_size[0]] label_patch label_map[y:y window_size[1], x:x window_size[0]] # 当标注区域内滑坡像素占比超过阈值时视为正样本 pos_ratio np.sum(label_patch 0) / (window_size[0] * window_size[1]) if pos_ratio 0.5: samples.append(img_patch) labels.append(1) elif pos_ratio 0.1: samples.append(img_patch) labels.append(0) # 保存样本和对应标签 for idx, (sample, label) in enumerate(zip(samples, labels)): class_dir os.path.join(save_dir, fclass_{label}) os.makedirs(class_dir, exist_okTrue) cv2.imwrite(os.path.join(class_dir, f{idx}.jpg), sample) return len(samples) # 用64x64窗口、32像素步长切割一张遥感影像 sample_count crop_image_with_labels( ./data/slope_img.tif, ./data/slope_label.png, window_size(64, 64), stride32, save_dir./train_samples ) print(f生成 {sample_count} 个样本块)这段逻辑的精髓在于两个参数。一是 window_size窗口太小滑坡地物的纹理特征展不开窗口太大正样本里混入大量非滑坡地物标签就不纯了。项目里的样本块尺寸我不确定是多少但如果你要自己调64×64 或 128×128 是比较稳妥的起步区间。二是 stride它直接决定样本量和重叠程度。stride 比 window 小相邻窗口有重叠样本量增大但容易产生相似的重复样本stride 等于 window 则完全不重叠样本独立性好但可能漏掉滑坡边界。2.2 正负样本失衡这个坑在数据组织阶段就要堵住滑坡场景分类的典型困境是正样本稀缺。一张遥感影像里滑坡区域可能只占百分之几你切出来的正样本块数量远少于负样本块。如果不处理后面训练出的 SVM 会无脑把所有样本判为负类因为这样准确率也很高。看项目里的文件列表training.txt、training2.txt 这些应该就是样本清单文件每行记录样本路径和类别标签。我的习惯是在 generate_samples_and_labels.py 里加一个负样本下采样逻辑import random def balance_samples(sample_list_file, output_file, max_neg_ratio1.5): 对样本列表做负样本下采样控制正负样本比例 sample_list_file: 原始样本清单每行格式为 路径 标签 max_neg_ratio: 允许的最大负样本/正样本比例 pos_samples, neg_samples [], [] with open(sample_list_file, r) as f: for line in f: path, label line.strip().split() if int(label) 1: pos_samples.append(f{path} {label}\n) else: neg_samples.append(f{path} {label}\n) # 随机抽样负样本使正负比不超过 max_neg_ratio max_neg int(len(pos_samples) * max_neg_ratio) if len(neg_samples) max_neg: neg_samples random.sample(neg_samples, max_neg) with open(output_file, w) as f: f.writelines(pos_samples neg_samples) print(f正样本 {len(pos_samples)} 张负样本 {len(neg_samples)} 张)为什么这里强调下采样而不是上采样影像样本维度高用 SMOTE 这类插值方法容易在特征空间里造出虚假样本对 SVM 这种依赖样本分布的分类器反而有害。直接随机丢弃部分负样本让正负比控制在 1:1.5 到 1:2 之间在绝大多数遥感场景分类任务里效果都够用。看着 training.txt 里那些路径和标签数字如果你发现正负样本数量悬殊优先考虑这个方案。3. 特征提取不是玄学光谱统计与 GLCM 纹理的代码级拆解样本切好了下一步是 feature.py 干的事。这个阶段的目标是把每个样本块从像素矩阵变成一行有意义的特征向量。项目明确写了提取两类特征光谱特征和 GLCM 纹理特征这两类特征在整个流程里是拼接在一起用的。3.1 光谱特征均值与标准差为什么够用滑坡区域在遥感影像上的典型表现是土壤裸露、植被覆盖减少在光谱上会呈现特定的反射特性。光谱特征计算的是影像块在每个波段上的像素值统计量最常见的就是均值和标准差。均值反映整体亮度水平标准差反映区域内像素值的离散程度。import numpy as np def extract_spectral_features(image): 从影像块中提取光谱特征 image: HWC格式的影像块channel在最后一维 返回每个波段的均值、标准差、归一化差值 features [] h, w, c image.shape for channel in range(c): band image[:, :, channel] # 均值描述整体灰度水平滑坡裸露地表的灰度普遍偏高 mean_val np.mean(band) # 标准差描述灰度波动滑坡区域内部灰度变化通常比较剧烈 std_val np.std(band) features.extend([mean_val, std_val]) # 如果是多光谱影像可以额外计算波段间关系 if c 3: # 归一化植被指数近似值红波段和近红外波段计算 red image[:, :, 0].astype(np.float32) nir image[:, :, 2].astype(np.float32) ndvi (nir - red) / (nir red 1e-6) features.append(np.mean(ndvi)) return np.array(features) # 假设一个64x64x3的样本块 sample_patch np.random.randint(0, 256, (64, 64, 3), dtypenp.uint8) spec_feat extract_spectral_features(sample_patch) print(f光谱特征维度: {spec_feat.shape[0]}) # 三波段影像输出7维特征3个均值3个标准差1个NDVI光谱特征提取最大误区是贪婪堆波段。RGB 三波段影像你提 6 维统计量就够了不需要把每个波段的直方图也拉进来。加 NDVI 是个不错的策略但前提是数据源里有近红外波段如果只是普通 RGB 航空影像近红外没有就别硬算。3.2 GLCM 纹理特征灰度共生矩阵的窗口、方向和统计量滑坡区域的光谱特征容易与裸土、道路混淆真正能区分开的是纹理。滑坡区域的纹理通常是破碎的、方向杂乱的这与人工建筑的规则纹理有本质区别。GLCM灰度共生矩阵通过统计像素对在特定方向和距离上的灰度联合分布来描述纹理。import numpy as np from skimage.feature import graycomatrix, graycoprops def extract_glcm_features(gray_image, distances[1], angles[0, np.pi/4, np.pi/2, 3*np.pi/4]): 提取GLCM纹理特征 gray_image: 单通道灰度图值域0-255 distances: 像素对距离滑坡纹理的特征通常在短距离内体现 angles: 四个方向的夹角弧度覆盖所有纹理朝向 返回每个方向上的对比度、相关性、能量、同质性 # 灰度级压缩到8级0-7压缩能抑制噪声干扰并加速计算 gray_norm (gray_image / 32).astype(np.uint8) glcm graycomatrix( gray_norm, distancesdistances, anglesangles, levels8, symmetricTrue, normedTrue ) props [contrast, correlation, energy, homogeneity] features [] for i in range(len(angles)): for prop in props: val graycoprops(glcm, prop)[0, i] features.append(val) return np.array(features) # 灰度压缩到8级后计算 gray cv2.cvtColor(sample_patch, cv2.COLOR_BGR2GRAY) glcm_feat extract_glcm_features(gray) print(fGLCM特征维度: {glcm_feat.shape[0]}) # 4个方向×4个统计量16维特征GLCM 里有几个参数值得细说。灰度级压缩到 8 或 16 是常规做法原始 256 级灰度计算出的共生矩阵极其稀疏统计量几乎没有判别力。方向的选择上四个方向全上比只算单个方向更稳健滑坡纹理没有固定朝向如果只取 0 度方向会漏掉垂直方向的纹理信息。距离参数设 1 即可因为灰度共生矩阵捕捉的是像素间的局部关系距离太远反而引入不相关的宏观纹理。最终每个样本块的完整特征向量是光谱特征和 GLCM 特征拼接而成这也是后面聚类和主题分析的输入。4. 从特征到视觉词袋再到 LDA 主题高级特征是如何炼成的把每个样本块变成一行特征向量后整个流程才走了不到一半。接下来的三个步骤——K-Means 聚类、视觉词袋构建、LDA 主题分析——层层递进把底层特征抽象成分类器真正能用的高级特征。这部分是整套代码里最能体现设计思路的环节。4.1 K-Means 聚类与视觉单词为什么要对特征做聚类直接拿 23 维特征向量去训练 SVM 不是不行但效果通常不够好。原因在于底层特征维度低、表达粗糙样本间区分度不足。K-Means 聚类在这里的作用是把全部训练样本的特征向量聚成若干簇每个簇中心就是一个视觉单词。当处理新样本时计算它与各聚类中心的距离得到的是一个离散化的词频分布相当于把连续特征映射到了更高层级的语义空间。from sklearn.cluster import KMeans import numpy as np def build_visual_vocabulary(feature_matrix, n_clusters100, random_state42): 用K-Means对全部训练样本的特征向量聚类生成视觉词典 feature_matrix: 形状为(n_samples, n_features)的特征矩阵 n_clusters: 聚类簇数即视觉单词数量 返回: 训练好的KMeans模型和每个样本的聚类标签 # 标准化后聚类效果更稳定避免量纲差异大的特征主导聚类 from sklearn.preprocessing import StandardScaler scaler StandardScaler() feat_scaled scaler.fit_transform(feature_matrix) kmeans KMeans( n_clustersn_clusters, initk-means, n_init10, max_iter300, random_staterandom_state ) cluster_labels kmeans.fit_predict(feat_scaled) # 统计每个样本的视觉单词频率构成词袋向量 vocab_size n_clusters bow_matrix np.zeros((len(feature_matrix), vocab_size)) for i, label in enumerate(cluster_labels): bow_matrix[i, label] 1 return kmeans, scaler, bow_matrix, cluster_labels # 假设有2000个训练样本每个样本是23维特征向量 train_features np.random.rand(2000, 23) kmeans_model, scaler_model, bow_feat, labels build_visual_vocabulary( train_features, n_clusters100 )聚类簇数的选择是比较关键的参数。簇数太少视觉单词过于笼统区分度不足簇数太多每个簇覆盖的样本太少统计不稳定。100 到 300 之间是这类任务的常见区间如果你样本量大可以往 300 以上靠。项目文件里的 visual_dict 和 zhang_spec_clt.pkl、zhang_text_clt.pkl 应该就是聚类器的序列化产物分别对应不同特征类型的视觉词典。4.2 LDA 主题分析从词袋到主题分布词袋向量的维度等于聚类簇数通常上百维且大量维度是稀疏的。LDA隐含狄利克雷分配在这里的作用是对词袋向量做第二次抽象把稀疏的高维词袋向量映射成低维的主题分布。注意LDA 传统上用于文本主题建模这里是把视觉单词当作词影像当作文档来用是跨领域迁移的典型做法。from sklearn.decomposition import LatentDirichletAllocation def extract_topic_features(bow_matrix, n_topics20, max_iter50): 用LDA对视觉词袋做主题分析得到每个样本的主题分布 bow_matrix: 形状为(n_samples, n_vocab)的词袋矩阵整数值 n_topics: 主题数量即高级特征维度 返回: 主题分布矩阵(n_samples, n_topics)和训练好的LDA模型 lda LatentDirichletAllocation( n_componentsn_topics, learning_methodbatch, max_itermax_iter, random_state42 ) # 主题分布矩阵每行代表一个样本的主题概率分布 topic_dist lda.fit_transform(bow_matrix) return lda, topic_dist # 用聚类得到的词袋特征做LDA bow_int bow_feat.astype(int) lda_model, topic_features extract_topic_features(bow_int, n_topics20) print(f主题特征维度: {topic_features.shape[1]}) print(f第一个样本的主题分布: {topic_features[0]}) # 输出是一个概率分布所有主题概率之和为1每一行主题分布是一个概率向量所有维度相加等于 1这意味着特征被归一化到了同一量纲对 SVM 训练是友好的。主题数量 n_topics 我一般设置在 10 到 30 之间太小语义区分不够太大容易出现噪声主题。项目的 generate_samples_and_labels.py 里应该有对应逻辑把样本路径、类别标签和主题分布三者组合成 Libsvm 标准的训练样本格式# Libsvm格式label index1:value1 index2:value2 1 1:0.023 2:0.081 3:0.012 ... 20:0.005 0 1:0.001 2:0.003 3:0.132 ... 20:0.0424.3 可视化地理解中间产物visual_lda.model 里装了什么项目根目录下的 visual_lda.model 文件看起来就是用 Gensim 训练的 LDA 模型。Gensim 的 LDA 模型保存后会生成多个文件.expElogbeta.npy 保存主题-词分布矩阵.id2word 保存词典映射.state 保存训练状态。如果你不想重新跑一遍 K-Means 和 LDA可以像下面这样直接加载from gensim.models import LdaModel import numpy as np # 加载训练好的LDA模型 lda LdaModel.load(./visual_lda.model) # expElogbeta是主题-视觉单词矩阵形状为(n_topics, n_vocab) topic_word_matrix np.exp(lda.expElogbeta) print(f主题数: {topic_word_matrix.shape[0]}, 视觉单词数: {topic_word_matrix.shape[1]}) # 对新的词袋向量预测主题分布 new_bow [(0, 1), (5, 2), (12, 1)] # 格式为(id, count) topic_dist lda.get_document_topics(new_bow, minimum_probability0)在调试阶段把主题数、视觉单词数打印出来比对能快速确认管线是否跑通了。如果不一致多半是先后步骤的 K-Means 配置对不上了。5. 避坑总览五条高频翻车记录与排查路径这套代码在毕设场景下算是跑得通的但真正自己动手时你大概率会在下面这几个地方栽跟头。每条都是实际运行过的路径按现象 → 原因 → 解决记录如下。5.1 训练集准确率接近 100%测试集崩到 60%这是最容易踩的坑。现象是模型在训练集上表现极好换成测试数据就大幅缩水。核心原因是样本切分时的数据泄露如果滑窗的 stride 小于 window_size相邻样本块存在重叠同时这些重叠样本被分到了训练集和测试集SVM 严格记忆了相似样本泛化能力自然差。解决方法是切样本时要保证同一位置的影像块不跨训练测试集或者对整张影像按空间位置划分数据集而不是对全部样本随机打乱。具体做法是先用一张完整影像的所有滑窗组成一个样本组整组进入训练集或测试集不混分。5.2 GLCM 计算报错all values in a row are equal现象是运行特征提取时抛出类似错误集中在全黑的影像块上。根本原因是灰度值完全相同的影像块无法计算灰度共生矩阵或者压缩到 8 级后所有像素落在一个灰度桶里。滑坡影像通常包含大面积暗色阴影区域这种全黑影块很常见。解决方法是两步走在特征提取之前先过滤方差过低的样本块比如灰度方差小于 5 的直接丢弃如果不想丢样本检查灰度压缩逻辑适当增大 levels 参数从 8 调到 16可以让灰度分布更稀疏但不至于退化到单桶。5.3 主题分布全是均匀值SVM 训练出来像猜拳现象是 LDA 输出的主题分布每一行都差不多比如全都是 0.05、0.05 这样均匀排列导致 SVM 根本找不到分类边界。原因通常是词袋矩阵过于稀疏K-Means 聚类后每个样本只有一两个视觉单词有词频大量维度为零LDA 在这种数据上学不出有区分力的主题。解决方法是回到 K-Means 阶段检查每个样本的聚类标签是否集中在少数几个簇里。如果你发现 80% 的样本落在 10 个簇内说明 n_clusters 设置过大需要减小同时考虑增加 LDA 的 max_iter从默认 50 加到 100 或 200让主题分布收敛得更充分。5.4 Libsvm 格式行列对应错位样本全部判成一类现象是训练正常但预测结果全为同一个标签检查中间文件时发现每一行的 index 和 value 对不上。原因是生成 Libsvm 格式文件时主题分布的特征顺序和标签没有对齐比如某行写了1 2:0.5 1:0.3index 没有按升序排列或者 index 从 0 开始而不是从 1 开始。Libsvm 对输入格式要求严格index 必须从 1 开始且递增。解决方法是生成文件后写个几行的校验脚本读取文件检查每行的 index 是否严格递增标签是否等于 0 或 1样本数量是否和源文件一致。项目里的 training2.txt 和 test2.txt 应该就是经过这一步处理的产物。5.5 预测阶段加载模型后维度对不上直接报错现象是用训练好的模型预测新样本时提示特征维度不一致或者聚类器输出和模型预期维度对不上。原因是预测时只加载了 SVM 模型却没有加载训练阶段的 StandardScaler 和 K-Means 模型。整个预测链路包含三个序列化对象scaler标准化、kmeans聚类、lda主题分析、svm分类器缺一个都会断裂。解决方法是把四个对象打包保存到一个 pickle 文件里或者用字典统一存储import pickle def save_pipeline(scaler, kmeans, lda, svm, save_path): 把完整推理链路统一序列化避免预测时缺组件 scaler: 特征标准化器 kmeans: 聚类模型 lda: LDA主题模型 svm: 训练好的Libsvm分类器 with open(save_path, wb) as f: pickle.dump({ scaler: scaler, kmeans: kmeans, lda: lda, svm: svm }, f) print(f完整pipeline已保存到{save_path})6. 复现后的下一步模型验证、整图预测与调参习惯跑通项目流程只是起点真正要把这套东西用在自己数据上还需要做三件事验证模型可靠度、对大幅影像做预测、理解参数调整的边界。交叉验证怎么落地。项目自带 test2.txt 测试集但这只能反映一次划分的结果说服力有限。我一般会用 5 折交叉验证重新评估主题分布加 SVM 的组合。具体做法是把训练样本按空间位置划分成 5 份每折用 4 份训练、1 份验证最终取 5 折平均准确率、精确率和召回率三个指标。需要特别注意交叉验证的划分单位是影像而不是样本块否则重叠样本会虚高评估结果。整图预测的滑窗策略。拿到一张新的遥感影像没有标注信息无法直接切样本这时需要滑动窗口逐块预测。窗口尺寸、步长与训练时完全一致每个窗口生成一个预测概率而不是硬标签。用概率而不是 0/1 硬标签的好处是可以设定阈值调整灵敏度比如滑坡是重点关注对象把阈值从 0.5 降到 0.3宁可误报也不漏报。可视化时把每个窗口的预测概率映射到色带就能生成热力图一眼看出模型认为哪些区域存在滑坡风险。调参的几个优先级。整套管线里参数优先级从高到低分别是样本块大小和步长、主题数量、聚类簇数、SVM 的 C 和 gamma。样本块大小影响特征质量是首要优化对象主题数量决定高级特征维度影响分类边界复杂度SVM 的 RBF 核参数最后才用网格搜索调优。完整跑一遍管线的耗时主要花在 K-Means 和 LDA 训练上如果样本量大先用小数据集验证 pipeline 通顺再全量训练。LDA 的训练还有个小习惯先用少量训练数据粗跑一次观察主题分布是否有区分度再全量训。从那以后我每次做这类场景分类任务都强制走一遍中间产物可视化聚类前输出特征分布直方图LDA 后输出主题分布的热力图SVM 训练后打印交叉验证矩阵。所有中间产物看一眼再进下一步看起来慢实际节省的时间比直接跑完三版多得多。希望帮到你。本文还有配套的精品资源点击获取