简介基于Matlab的遥感图像变化监测课程设计完整源码包面向测绘遥感、计算机、电子信息等专业的课程设计与毕业设计场景适合需要完成变化检测实验或掌握CVA、PCA、NDVI等算法流程的学生与开发者。压缩包共29个文件主体为24个.m源码脚本覆盖影像配准、变化矢量分析、主成分分析、NDVI计算、阈值分割等关键模块既有基础算法实现也包含改进与对比实验另有Markdown说明文档和部署说明便于快速理解工程结构并附带DOCX与PDF格式实验报告可对照学习实验设计与结果分析还有独立zip素材包辅助完整复现。资源仅9.19MB轻量易用目前已有119人学习浏览。项目源码经过完整测试答辩评分较高具备一定工程性可在其基础上扩展修改用于课程答辩、毕业设计或日常算法练习。1. 从两期遥感影像里找真的变了为什么这么难很多人第一次做遥感图像变化监测是直接把两期影像的像素相减差值大于某个固定值就标成变化区结果植被区、水体边缘、云影边界冒出一大片伪变化。原因不复杂卫星两次过境时的太阳高度角、大气条件、传感器视角都不同地面几何位置也有偏移。真正的变化监测流程是配准、辐射归一化、特征变换、变化向量计算、阈值分割。这个基于Matlab的课程设计源码把CVA、PCA、NDVI_BI_CVA三种检测算法、四套配准函数、自动阈值搜索和ROI框选工具都串起来还附带实验报告和部署说明。适合正在做课程设计、大作业的本科生也适合想快速复现一套变化监测测试流程的研究生和工程师。下文我会拆开讲怎么组合、参数怎么调以及哪些文件是你应该自己动手重写的。2. 变化监测的算法选型CVA、PCA与NDVI_BI_CVA2.1 CVA变化向量不是差值图CVA全称Change Vector Analysis。它的输入是同一地理位置两时相影像的像元光谱向量。假设第一时相像元向量是x1[b1,b2,...,bn]第二时相是x2变化向量就是Δx2-x1。Δ的模长反映变化强度方向反映变化类型。比如蓝光和红波段同时变大可能对应土地裸露近红外显著下降而红光上升多半是植被减少。这就是CVA比单纯波段差值图信息量大的地方。源码里的CVA180322A1.m、CVA0323.m、CVA0325A1.m、CVA0325A2.m是不同实验阶段留下的版本核心逻辑基本一致。我一般在课程设计里保留一个主版本把其他版本放进实验记录目录。核心计算可以这样写function [change_mag, change_dir] cva_core(img1, img2) img1 double(img1); img2 double(img2); diff img2 - img1; % 各波段逐像素相减 change_mag sqrt(sum(diff.^2, 3)); % 变化强度向量模长 % 方向角以红波段和近红外波段为例 red diff(:, :, 3); % 假设第3波段为红波段 nir diff(:, :, 4); % 假设第4波段为近红外 change_dir atan2(nir, red); % 方向角单位弧度 end这里img1和img2必须已经完成几何配准否则逐像素减法没有意义。diff是所有波段组成的差值立方体sum(diff.^2,3)沿第三维求和得到每个像素的变化强度。atan2返回的角度范围是[-pi,pi]可以用它给变化类型上色比如植被退化集中在一个角度区间。实际使用中如果波段数不固定不建议写死波段序号最好在调用前读取影像元数据确定红波段和近红外波段的位置。CVA的输出要配合阈值才能得到二值变化图。阈值选多少直接决定结果图是漏检还是爆炸式误检。后面第4章会专门说find_best_threshold怎么用。2.2 PCA把相关波段变成不相关的几个分量PCA法也被称为主成分变化检测。基本思路是把两时相影像的各波段堆叠成一个高维数据集然后做主成分分析让变化信息集中到前几个主成分上。源码里PCA_180321A2.m、PCA_180321A3.m和PCA_analyse.m就是这套流程的不同实现。Find_K_Max_Eigen.m负责从特征值里挑选前K个贡献率最大的分量这个K通常取2到4。PCA的好处是它不需要先验指定哪个波段代表什么地物完全靠数据协方差驱动。坏处是如果两时相影像的辐射差异很大第一主成分往往被整体亮度差占据真实地物变化反而被压到后面。所以很多课程设计会先用直方图匹配或回归归一化再做PCA。下面是一段可跑的示例框架function [pc_img, eig_vals] pca_change(img1, img2) [rows, cols, bands] size(img1); data reshape(double(cat(3, img1, img2)), [], bands * 2); data data - mean(data, 1); % 中心化 cov_m cov(data); % 协方差矩阵 [V, D] eig(cov_m); [eig_vals, idx] sort(diag(D), descend); V V(:, idx); scores data * V(:, 1:4); % 取前4个主成分 pc_img reshape(scores, rows, cols, 4); endcat(3,img1,img2)把两时相波段拼接成2*bands维特征空间中心化后求协方差eig返回特征向量和特征值。注意eig返回的特征值不一定是降序所以要sort。实际项目中主成分得分图还要做阈值化也可以把CVA和PCA结合起来对主成分得分差做CVA源码里的PCA_CVA0322A2.m就是这种混合思路。我的建议是当你的影像只有3到4个波段时用CVA就够了当波段数多比如高光谱或哨兵多波段时先PCA降维再计算变化速度更快也更稳。另外源码中的KT.m是Kauth-Thomas变换和PCA类似也是正交变换但变换系数是固定的专门针对Landsat波段设计。如果答辩时老师问为什么不用KT你可以说KT系数依赖传感器跨传感器数据用PCA更通用。2.3 NDVI_BI_CVA植被指数与亮度指数的组合NDVI_BI_CVA.m把两种指数同时塞进CVA框架。NDVI能突出植被变化BIBrightness Index能突出裸土、建筑物等亮度变化。两时相分别计算NDVI和BI得到两个指数差值再合成一个二维变化向量。这样做的直观理由是只用NDVI只能发现植被变化城市扩展、土壤翻动这类事件虽然NDVI也在变但特征不够明显加上BI后地物变化类型更容易分离。一个常见做法是ndvi1 (nir1 - red1) ./ (nir1 red1 eps); ndvi2 (nir2 - red2) ./ (nir2 red2 eps); bi1 sqrt(sum(img1(:, :, 1:3).^2, 3)); bi2 sqrt(sum(img2(:, :, 1:3).^2, 3)); delta cat(3, ndvi2 - ndvi1, bi2 - bi1); mag sqrt(sum(delta.^2, 3));eps是防止分母为0。BI取可见光波段的平方和的根代表像元整体亮度。如果你处理的是Landsat影像注意先把各波段缩放到同一个尺度比如都转成反射率或归一化到[0,1]否则BI会被数值大的波段带偏。NDVI同理先做辐射归一化再计算不然两期NDVI的差会包含大量大气噪声。2.4 三种算法怎么选算法输入要求适用场景判断依据CVA波段数一致、已配准土地利用/覆盖变化多时相多波段变化向量模长方向角PCA波段数较多高光谱、多光谱批量筛选主成分差NDVI_BI_CVA含红/近红外波段植被减少、城市扩展监测指数差值向量模长实际项目中我建议先用CVA跑一遍如果结果图中变化区域过于破碎再换NDVI_BI_CVA不要一开始就上PCA因为PCA的主成分含义很难向答辩老师解释清楚而CVA和NDVI的物理含义更直观。用于课程设计演示时三种算法都跑一遍把三张结果图放一起对比本身就是很好的加分项。3. 先配准再检测SURF、Nonrigid、MOVINGREG 三条路的取舍3.1 为什么配准脚本占了源码一大半这个项目的文件列表里registerImagesSURF.m、registerImagesNonrigid.m、registerImagesMOVINGREG.m、registerImagesMULTIMODAL.m四个文件都叫registerImages对应同一输入输出接口但内部配准策略不同。之所以要这么多版本是因为遥感影像的时相差异不只是位移可能有旋转、缩放、甚至局部形变。SURF适合整体刚体变换Nonrigid适合地形起伏和传感器角度造成的局部扭曲MOVINGREG和MULTIMODAL适合不同传感器来源的影像比如光学和SAR或者光学与历史地图。3.2 registerImagesSURF 的特征点路线SURF特征点配准步骤在Matlab里可以写成一串函数调用。主流程保留在registerImagesSURF.m中核心逻辑如下function [movingRegistered, tform] registerImagesSURF(moving, fixed) grayMoving im2gray(moving); grayFixed im2gray(fixed); pointsFixed detectSURFFeatures(grayFixed); pointsMoving detectSURFFeatures(grayMoving); [featFixed, validFixed] extractFeatures(grayFixed, pointsFixed); [featMoving, validMoving] extractFeatures(grayMoving, pointsMoving); indexPairs matchFeatures(featFixed, featMoving); matchedFixed validFixed(indexPairs(:, 1), :); matchedMoving validMoving(indexPairs(:, 2), :); [tform, inlierIdx] estimateGeometricTransform2D(... matchedMoving, matchedFixed, similarity); movingRegistered imwarp(moving, tform, ... OutputView, imref2d(size(fixed))); enddetectSURFFeatures返回的特征点数量受MetricThreshold参数控制默认在1000左右。如果影像纹理太少可以把该参数调小到100让更多弱特征进入匹配。extractFeatures之后用matchFeatures默认的最近邻比率匹配比率阈值默认0.6越大匹配对越多但误匹配也越多。estimateGeometricTransform2D中的similarity模型允许平移、旋转、等比例缩放适合大多数Landsat/Sentinel影像的帧间校正如果两张图尺寸差太多改affine会更稳。最后一个注意点是imwarp后图像尺寸要和fixed一致用imref2d(size(fixed))控制输出视窗这是很多新手容易漏掉的一步——不写OutputView输出图和参考图的大小就对不齐。提示如果你的Matlab版本在R2020b之前detectSURFFeatures仍然可用但im2gray需要换成rgb2gray且输入必须是三通道影像。3.3 Nonrigid 与 MULTIMODAL什么时候不要用SURF如果两期影像之间有山谷、山脊这类地形起伏造成的局部偏移整体仿射变换无法消除错位。此时可以用registerImagesNonrigid.m内部通常走disparity估计或配准网络得到逐像素位移场。代价是参数多、跑得慢而且位移场平滑太少会产生空洞。另一个场景是参考图和待配准图来自不同传感器比如光学与SAR或者RGB航拍与历史地形图。SURF的灰度特征在跨模态情况下往往提取不到足够匹配对。registerImagesMULTIMODAL.m里常见的做法是先用互信息作为相似性度量或者用相位相关做粗配准。我给课程设计的建议是同一传感器的两期光学影像优先用SURF不同传感器的影像直接用registerImagesMULTIMODAL只有做完SURF后看到明显局部错位再上Nonrigid。不要三个函数轮流试浪费时间还容易过拟合。3.4 配准后如何肉眼验证配准质量直接影响后续CVA的假阳性。干巴巴看配准前后两幅图很难判断好坏。可以用imshowpair输出叠加图imshowpair(fixed, movingRegistered, falsecolor); title(Falsecolor overlay after registration);falsecolor会把两幅图分别映射到青色和品红色如果地物边缘完全重合图像呈灰色出现彩色边缘就说明还有亚像元级错位。对于课程设计报告我通常会在报告里放3张图配准前overlay、控制点连线、配准后overlay。有这三张图老师就能明白你是真的做了配准而不是只贴一句使用SURF配准。4. 阈值与ROI把变化图变成能交差的结果图4.1 find_best_threshold 在找什么CVA或者PCA输出的是连续变化强度图。要得到二值变化图需要阈值T。T太小噪声全部变成变化区T太大真实变化被滤掉。find_best_threshold.m做的就是这件事。最常见的实现是大津法Otsu它把变化强度直方图分成前景和背景两类让类间方差最大。课程设计如果要求不高直接对change_mag调用graythresh即可如果想体现工作量就在find_best_threshold里写一个遍历算子function bestT find_best_threshold(mag) magVec mag(:); [counts, edges] histcounts(magVec, 256); total sum(counts); omega0 cumsum(counts) / total; % 前景累积占比 mu cumsum(counts .* edges(1:end-1)) / total; % 累积灰度均值 muT mu(end); sigmaB (muT * omega0 - mu).^2 ./ max(omega0 .* (1 - omega0), eps); [~, idx] max(sigmaB); bestT edges(idx); end这段代码遍历256个灰度级计算每个灰度作为阈值时的类间方差sigmaB取最大值对应的灰度作为bestT。eps防止分母为0。实际遥感影像的变化强度直方图往往不是标准双峰分布Otsu可能把阈值往高灰度偏移。这时可以给bestT乘以一个系数比如bestT bestT * 0.8宁可多保留一点候选区再用形态学开运算去掉小斑点。4.2 只统计感兴趣区roicircle 与 inpolygon很多课程设计的数据集里包含水体、云层等不想参与评价的区域。源码中的roicircle.m应该是用来交互式圈圆形感兴趣区的inpolygon.m和myinpolygon.m则是判断点是否落在多边形内。使用时可以先手工画几个多边形把变化检测限制在有效范围内roiMask false(size(changeMask)); for k 1:numel(roiPolygons) [xq, yq] meshgrid(1:cols, 1:rows); in inpolygon(xq, yq, roiPolygons{k}(:, 1), roiPolygons{k}(:, 2)); roiMask roiMask | in; end validChange changeMask roiMask;inpolygon返回与xq/yq同尺寸的逻辑矩阵遍历多个多边形用或运算合并。注意矩阵尺寸是rows*cols但meshgrid的xq对应的是列坐标yq对应行坐标别把x和y传反。另一个常见问题是roiPolygons里的坐标是地理坐标还是图像像素坐标——图像坐标直接用经纬度会全部落在ROI外。处理遥感数据时先用geotiffread读取地理参照信息把经纬度转换成像素行列号再传给inpolygon。其实Matlab的roipoly函数可以交互式选多边形比手写inpolygon更省事源码里的myinpolygon.m我猜是为了避免调用Image Processing Toolbox而写的纯Matlab实现适合在未安装工具箱的机器上跑。同目录的inpolygon_test.m就是用来验证多边形判断是否正确的测试脚本。4.3 变化检测结果的定量评价课程设计答辩时老师大概率会问你怎么证明你的结果是对的。所以除了画图还要给出定量指标。常见做法是把已有土地利用分类图作为真值和你的变化检测结果做混淆矩阵算准确率、误检率、漏检率、Kappa系数。如果没有真值可以人工采样验证点在原始影像上随机生成几百个点人工判断是否真的变化。下面是计算混淆矩阵的快速方法pred changeMask(:); true groundTruth(:); TP sum(pred true); FP sum(pred ~true); FN sum(~pred true); precision TP / (TP FP eps); recall TP / (TP FN eps); F1 2 * precision * recall / (precision recall eps);这里的groundTruth加载后要和changeMask尺寸一致且都用逻辑类型。报告里写清楚采样点数、空间分辨率、评价尺度老师就不会觉得你的结果是拍脑袋。项目自带的任春哲-201511190114-变化检测实验报告.docx/pdf里估计就有类似流程你可以直接照着改但别原样抄——具体影像和处理参数不一样数据是能对出来的。5. 把课程设计改造成批处理脚本从单次运行到多景影像自动监测课程设计归档里的文件大多是单张运行脚本比如script0320A1.m、script0320A2.m、Script180322A3.m一看就是实验过程中每天一个版本。这种命名适合研究记录但不适合产品化。如果你想把基于Matlab遥感图像的变化监测写成简历项目至少要把它改成循环跑多组影像、自动保存结果的工作流。5.1 数据组织固定目录和命名我习惯建立如下目录结构data/ scene1/ before.tif after.tif scene2/ before.tif after.tif result/然后用dir(data/*)遍历场景目录。假设每个场景的before和after都来自同一传感器且已经做了几何精校正批处理主循环可以这样写scenes dir(data); for i 1:length(scenes) if ~scenes(i).isdir || startsWith(scenes(i).name, .) continue; end sceneDir fullfile(data, scenes(i).name); before geotiffread(fullfile(sceneDir, before.tif)); after geotiffread(fullfile(sceneDir, after.tif)); [mag] cva_core(before, after); T find_best_threshold(mag); changeMask mag T; imwrite(changeMask, fullfile(result, [scenes(i).name _mask.tif])); endgeotiffread读入的如果是多波段需要确认波段顺序是BGRN还是RGBfind_best_threshold返回的是灰度阈值如果mag是单精度浮点最好先归一化到[0,1]再比较。imwrite写成TIF时加一个WriteMode,overwrite避免文件已存在报错。每跑完一个场景用fprintf输出场景名、阈值、变化像元数量方便对照日志排查问题。如果数据是ENVI格式把geotiffread替换成源码里的freadenvinew.m逻辑是一样的。5.2 记录参数与边界条件批处理最容易遇到的问题是某个场景尺寸不一致或者某张图片没有地理参考。我在循环里加了一个try-catch捕获异常后把出错场景写入skip.log不中断整体流程。此外每个场景的阈值T会被自动计算但不同时间的影像辐射尺度如果差异很大直接共用全局阈值不合理。建议在批处理前先对每对影像做一次直方图匹配用Matlab的histeq或imhistmatch把待监测影像的直方图映射到参考影像。这样find_best_threshold得到的阈值才具有可比性。最后提醒一个细节批量跑完后用geotiffwrite保存带投影信息的结果栅格而不是用imwrite保存无坐标信息的普通tif。因为你后续可能要导入GIS做叠加分析没有坐标参考的数据还得重新配准一次。geotiffwrite需要传入Rspatial reference object由geotiffread返回这个细节是很多课程设计扣分的地方。把这段流程接到你已有的main脚本末尾整套系统才算闭环。本文还有配套的精品资源点击获取