尧图网页设计 
简介面向遥感影像变化检测研究的经典算法代码包完整实现了IR-MAD、MAD、CVA、PCA四种方法均为Matlab可直接运行的脚本可解决多时相影像中地物变化识别、异常区域探测和光谱特征降维等问题适用于土地利用、环境监测、城市规划等研究方向需要读者具备基础遥感与矩阵运算知识。压缩包共97个文件、10.4MB以15个m源码为核心搭配28个bmp、24个tif测试影像以及hdr头文件和fig结果图源码中既包含IRMAD_Update、MADGet、CVADemo、PCADemo等主程序也有DAcom、covw、eigen2等辅助函数并附有泰州地区真实遥感样例数据可直接执行得到强度图和二值变化图。四种算法各有侧重IR-MAD和MAD对光照、大气影响更鲁棒CVA输出直观PCA擅于去噪和降维实际应用中可根据影像质量、变化类型和结果解释复杂度灵活选择包内另有Change Result Comparing.m等对比脚本方便输出定量比较结果。截止目前已有2008人学习适合遥感方向学生、科研人员及工程师快速上手经典算法并开展实验验证。1. 为什么“老三样”变化检测算法到今天还没过时遥感影像变化检测这几年被深度学习刷了很多屏但你去看实际生产的生态督查、违建监测、灾害评估项目跑在服务器上最稳的仍是MAD、IR-MAD、CVA、PCA这批经典算法。它们不吃海量训练样本不需要GPU你把两期影像丢进来半小时内就能出一张可交代的变化图。这篇文章就把这四个算法串起来讲清楚每个方法在做什么数学假设、代码怎么落地、参数怎么调、以及我实际踩过的那些坑。2. MAD与IR-MAD变化检测里的“差分增强器”2.1 MAD为什么比直接做差靠谱典型相关分析的直觉直接做差值有个根治不了的毛病——两期影像的辐射条件很难完全一致。传感器响应、太阳高度角、大气水汽任何一个变量变了差值图上就会出现大面积伪变化。MADMultivariate Alteration Detection的思路是不直接比较灰度而是先对两个多波段影像各做一次线性组合让组合后的分量之间相关性最强再去比较这两个组合分量的差。用数学语言描述就是找两组典型变量U a^T XV b^T Y最大化corr(U, V)。然后取M_i U_i - V_i作为MAD分量。为什么这样能压制辐射差异因为典型相关只保留了两个影像中“共同变化”的模式而共同变化以外的部分——比如传感器噪声、光照不均、大气影响——被踢到了残差里。理论上MAD分量服从卡方分布所以后续用卡方分布做显著性检验有天然的统计依据。这是它和直接差值的本质差别直接做差是像素级操作MAD是多元统计变换。我见过不少新手拿到两期Landsat影像上来就做波段差值结果变化图上全是耕地翻耕的条带和大气差异造成的噪声块。换成MAD之后这些噪声被自动过滤掉剩下的变化斑块才是真正需要关注的地表状态改变。2.2 最小可复现的MAD用numpy实现最简实现不需要遥感软件把两期影像的像元光谱当成两个随机向量就行。先把影像重排成像素数×波段数矩阵然后做典型相关分析import numpy as np def mad_transform(x, y, max_components3): # x, y: 像素数 x 波段数 的二维矩阵波段数通常3~8 # 返回mad分量矩阵每一列是一个mad分量 n_pixels x.shape[0] # 去均值典型相关分析要求零均值数据 xc x - x.mean(axis0) yc y - y.mean(axis0) # 计算协方差矩阵加正则项防止奇异 sxx np.cov(xc, rowvarFalse) np.eye(x.shape[1]) * 1e-6 syy np.cov(yc, rowvarFalse) np.eye(y.shape[1]) * 1e-6 sxy np.cov(xc, yc, rowvarFalse)[:x.shape[1], x.shape[1]:] syx sxy.T # 解广义特征值问题Sxy * Syy^-1 * Syx mat np.linalg.pinv(sxx) sxy np.linalg.pinv(syy) syx eigval, eigvec np.linalg.eig(mat) # 按特征值降序排序特征值越大相关性越强 idx np.argsort(eigval.real)[::-1] eigvec eigvec[:, idx] # 取前max_components个典型向量 a_vec eigvec[:, :max_components].real # 计算对应的b向量 b_vec np.linalg.pinv(syy) syx a_vec b_vec b_vec / np.sqrt((b_vec ** 2).sum(axis0)) # 计算典型变量并做差 u xc a_vec v yc b_vec mad_components u - v return mad_components这段代码的核心在广义特征值求解。mat矩阵的特征向量的意义是把x投影到某个方向后和y的某个方向相关性最大。注意我在sxx和syy对角线上加了1e-6这是防止协方差奇异。对实际遥感影像来说波段间相关性极强协方差矩阵很容易接近奇异不加这个正则项特征值求解会直接报错或者给出发散的结果。参数说明max_components一般取波段数或波段数减1。对4波段影像取3个MAD分量就够第四个分量对应的典型相关特征值趋近于1数学上的意义是“两期影像该方向上几乎无变化”直接舍弃。eigvec前几个分量的特征值应当接近1但又小于1如果你解出来的特征值等于1.0000000多半是正则化没加或者影像有大量无效值比如黑色边界没剔除。2.3 从MAD到IR-MAD迭代权重是怎么提精度的MAD有个明显的缺陷——它假设两期影像所有像素都是“没变化”的并在这个假设下估计典型相关结构。但真实场景里有变化的像素虽然不多可总有那么几个如果变化区域过大比如一场洪水淹了几万平方公里协方差矩阵会被真实变化污染典型相关计算出来的方向就不再是“共同稳定模式”了。IR-MADIteratively Reweighted MAD就是为这个设计的每次迭代根据当前MAD分量计算每个像素属于“no-change”的权重然后带权重新计算协方差矩阵再做一次MAD。变化像素的权重会越来越小no-change像素的主导地位就被加强了。按文献惯例权重按下式计算w_i P(卡方统计量 z_i)其中z_i sum((mad_k / sigma_k)^2) 是各MAD分量归一化后的平方和卡方分布自由度取分量数。from scipy.stats import chi2 def ir_mad(x, y, max_iter10, tol1e-4): # x, y: 像素数 x 波段数, 必须用有效值掩膜剔除nodata n_bands x.shape[1] # 初始权重全为1即第一轮等同于普通MAD w np.ones(x.shape[0]) for i in range(max_iter): # 带权计算均值和协方差 xm (x * w[:, None]).sum(axis0) / w.sum() ym (y * w[:, None]).sum(axis0) / w.sum() xc (x - xm) * np.sqrt(w[:, None]) yc (y - ym) * np.sqrt(w[:, None]) sxx xc.T xc / w.sum() np.eye(n_bands) * 1e-8 syy yc.T yc / w.sum() np.eye(n_bands) * 1e-8 sxy ((x - xm).T (w[:, None] * (y - ym))) / w.sum() syx sxy.T # 解广义特征值问题与普通MAD相同 mat np.linalg.pinv(sxx) sxy np.linalg.pinv(syy) syx eigval, eigvec np.linalg.eig(mat) order np.argsort(eigval.real)[::-1] a_vec eigvec[:, order].real # 用全部特征向量计算所有MAD分量 u (x - xm) a_vec v (y - ym) (np.linalg.pinv(syy) syx a_vec) mad u - v # 对方差归一化每个分量尺度统一 sigma mad.std(axis0) 1e-12 z ((mad / sigma) ** 2).sum(axis1) # 卡方分布计算no-change权重 new_w 1 - chi2.cdf(z, dfn_bands) # 收敛判断 if np.abs(new_w - w).max() tol: w new_w break w new_w return mad, w这里有个细节MAD分量计算出来后要对每个分量做方差归一化因为典型相关只保证相关性最大不保证分量尺度一致。z实际上是一个卡方统计量在no-change假设下应服从自由度等于分量数的卡方分布所以用1 - chi2.cdf(z)作为变化概率的补数也就是no-change权重。权重小说明该像素大概率真的变了这个权重向量本身就是一张很好的变化概率图。迭代停止条件用max_iter和tol双保险。我一般设max_iter10实际三五次迭代权重就收敛了。如果你发现迭代到第十次权重还在明显变化先别急着加迭代次数回去查一下影像预处理大概率是两期影像配准上有系统偏移或者有一条波段没做辐射定标只做了快速大气校正。2.4 IR-MAD的关键参数与停止条件IR-MAD名义上参数少实际上有两个地方非常影响结果。第一个是协方差矩阵的奇异值处理。遥感波段数通常48但像素数动辄几千万协方差矩阵肯定能求逆。坑不在这坑在无效值。影像四周的黑边、云掩膜产生的0值如果不做掩膜直接参与统计协方差会被拉向零值典型相关方向就是错的。我处理高分一号、高分二号的经验是先做一次有效值掩膜波段值0且不是nodata掩膜后的像素才进IR-MAD。第二个是迭代里的正则项系数。上面代码里sxx加的是1e-8这个值看数据情况。如果你发现前后两次迭代的权重分布出现振荡——第一次权重集中在A区域第二次却跑到B区域那多半是1e-8太小导致广义特征值求解数值不稳定。这时把正则项调到1e-6甚至1e-5振荡基本能消掉。代价是权重图会变钝一点但对最终二值化影响有限。3. CVA与PCA两个“先做变换再定阈值”的经典流派3.1 CVA的变化向量方向比大小更值钱CVAChange Vector Analysis的思路比MAD直白得多把每个像素看成多维光谱空间里的一个点两期影像的同一像素构成两条光谱向量直接做差得到一个“变化向量”。变化向量的模长就是变化强度变化向量的方向就是变化类型。这个算法最值得称道的地方是它把“变没变”和“变成了什么”分开了。模长超过阈值的像素被判定为变化区域而方向向量可以用来区分变化类别植被变裸地方向角一般在某个范围水体变建设用地又在另一个范围。实际项目里我常拿CVA当第一道粗筛先出整景变化强度图然后按方向角把变化像素聚类成几类再做分类后处理。但CVA的缺陷同样明显——它对辐射归一化的要求极高。直接做差意味着任何大气条件差异都会以“变化”的形式出现在模长里。所以用CVA前两期影像必须做至少一次直方图匹配或相对辐射归一化。这跟MAD不同MAD在数学上自带对公共模式的过滤CVA完全裸奔。3.2 CVA的Python实现从波段差到变化强度import numpy as np def cva_change_detection(x1, x2): # x1, x2: 已配准的两期多光谱影像, 形状为 h x w x n_bands # 先做相对辐射归一化: 用线性回归把x2匹配到x1的辐射范围 from sklearn.linear_model import LinearRegression h, w, nb x1.shape x1_flat x1.reshape(-1, nb).astype(np.float64) x2_flat x2.reshape(-1, nb).astype(np.float64) # 逐波段做线性归一化 x2_norm np.zeros_like(x2_flat) for b in range(nb): lr LinearRegression().fit(x2_flat[:, b:b1], x1_flat[:, b:b1]) x2_norm[:, b] lr.predict(x2_flat[:, b:b1]) diff x1_flat - x2_norm # 变化强度 欧氏距离, 变化方向 归一化的差值向量 magnitude np.sqrt((diff ** 2).sum(axis1)) direction diff / (magnitude[:, None] 1e-10) # 重排回影像形状 mag_img magnitude.reshape(h, w) dir_img direction.reshape(h, w, nb) return mag_img, dir_img这段代码最容易被忽略的是前面的线性归一化。很多初学CVA的帖子不写这一步直接用原始DN值做差结果在山区影像上几乎全花——因为地形阴影的辐射差远大于真实的地表变化。线性回归把x2整体调整到x1的辐射水平虽然不是严格意义上的大气校正但对付两期影像间的系统辐射偏移是管用的。参数说明LinearRegression逐波段做意思是一个波段一个增益和一个偏置。更高档的做法是考虑邻近像元光谱关系做相对辐射归一化但线性回归已经能消掉80%的系统辐射差异。如果你发现归一化后变化强度图还是噪点密布把影像先做一次3×3的均值滤波信噪比会明显改善——这是CVA最便宜的去噪手段。另外要注意direction的计算模长接近0的像素方向不稳定实际生产里我会把模长低于阈值的像素方向直接置为0避免后面的方向聚类被噪声主导。3.3 PCA做变化检测对差值影像降维取主要变化PCA进入变化检测的方式有点特殊不是对两期影像直接做主成分而是先做差值影像通常是CVA的差值然后对差值影像的各波段做主成分分析。PCA在这里的角色是“信息浓缩器”把n个波段的变化信息压缩到前几个主成分里第一主成分就是变化的主导方向方差贡献率越高说明变化越集中在一两种模式里。这个思路在很多项目里被用来区分真变化和噪声噪声在主成分上的分布通常均匀而真实的地表变化往往集中在第一、第二主成分上。所以操作上就是取前两个主成分的模长作为变化强度。打个比方这就跟PCA做特征脸类似——先用主成分提纯信号把分散在多维空间里的主要模式抽出来剩下的尾项当成噪声丢掉。对变化检测来说丢掉的那些尾项里确实藏着大量传感器噪声。def pca_change_detection(diff_flat, n_components2): # diff_flat: 像素数 x 波段数 的差值矩阵 # 去均值 d diff_flat - diff_flat.mean(axis0) # 计算协方差并做特征分解 cov np.cov(d, rowvarFalse) eigval, eigvec np.linalg.eigh(cov) # 特征值升序排列, 取最大的n_components个 idx np.argsort(eigval)[::-1][:n_components] proj d eigvec[:, idx] # 变化强度 前n_components个主成分的模长 change_mag np.sqrt((proj ** 2).sum(axis1)) return change_mag, eigval[idx]这里用np.linalg.eigh而不是np.linalg.eig因为协方差矩阵是对称阵eigh更稳定也更快。返回值里的eigval[idx]是前两个主成分的特征值它们占总特征值和的比率就是方差解释率。我看到很多人在报告里只写“PCA提取了主要变化”连解释率都不给——这个数字其实特别重要如果前两个主成分只解释了不到60%的方差说明变化模式很散PCA路线不适合这个场景别硬用。3.4 PCA的代码与“主成分数”怎么定主成分数在变化检测里一般取2到3个就够。不是越多越好多取一个主成分多带进一部分噪声。判断依据是特征值贡献率曲线也就是“陡坡图”。理想情况下曲线在前两三个点上急速下降长尾平缓拐点处就是该截断的位置。如果曲线一直平缓没有明显拐点说明差值影像本身没有主导的变化方向这时PCA和直接取差值模长没有本质区别。实际使用中还有一个变体——把PCA和MAD结合先做MAD得到几个分量再对这些分量做PCA压缩成一维或两维的变化指标。这个做法适合波段多的数据比如Sentinel-2的10米波段MAD已经把辐射差异压制了PCA再降一次维可视化效果比直接看卡方统计量更直觉。我最近处理一个城市扩张项目就是这么干的MAD出4个分量PCA压缩到2维变化强度图比单用MAD的卡方值干净而且能隐约看出不同扩张方向在颜色上的区分。4. 四个算法参数怎么设一张表和四类场景建议4.1 参数对比表MAD/IR-MAD/CVA/PCA算法核心步骤主要参数对辐射归一化的要求输出计算开销CVA逐像素光谱差阈值方向聚类角度极高必须先归一化变化强度图方向图最低纯逐像素运算PCA差值影像去均值特征分解主成分数解释率阈值高归一化与否影响提取方向主成分变化图低一次特征分解MAD典型相关分析差分max_components正则项低共同模式被自动过滤MAD分量卡方统计量中等需构造协方差矩阵IR-MAD带权迭代典型相关迭代次数正则项收敛容差低迭代中动态调整权重MAD分量权重图高多次协方差与特征分解这张表的第一列和第二列值得细读。CVA看着最“朴素”但对预处理的依赖最重MAD系列数学上最优雅但代码和调参的细节最多。PCA居中它能不能用取决于你差值影像的质量而差值影像的质量又取决于CVA那步归一化做得好不好。所以实际项目里这四个算法不是互相替代的关系而是流水线里的不同环节先归一化再做差值或MAD再用PCA或CVA出强度图最后二值化。4.2 不同数据场景的选型建议场景一只有两期中低分辨率多光谱Landsat级30m。优先IR-MAD。这类数据辐射一致性较好但云和阴影的多时相差异明显IR-MAD的权重迭代能在迭代过程中逐渐压掉云影伪变化。我做过一个2000到2010年Landsat的城市扩张监测IR-MAD比MAD直接出的变化斑块干净很多尤其是山区阴影带。场景二两期高分辨率影像0.52m有云阴影和建筑物阴影。CVA不乐观PCA也不乐观——高分辨率影像的地物阴影变化太剧烈。正确顺序是先做阴影掩膜掩膜外区域再做IR-MAD。如果实在想用CVA必须配合逐块本地阈值而不是全局阈值不然阴影边缘全是伪变化。场景三多光谱波段特别多8个以上比如Sentinel-2部分波段或WorldView-2。MAD系列的问题在于波段间相关性会迅速变强广义特征值求解的数值稳定性会出问题。务必把正则项加大到1e-5量级且只取前46个波段参与计算宁可用波段子集也不要一股脑把12个波段全塞进去。场景四只有单波段雷达数据或高程数据。MAD和CVA都需要多波段单波段可以用CVA的退化版——直接做差取绝对值或者用时空主成分分析把多期影像的时间维当成波段维做PCA。这个方法在地表沉降、水体面积监测里很常见。5. 变化检测常见坑从配准误差到阈值玄学5.1 影像没配准变化检测全白做现象变化检测结果图上一半检测出来的“变化”沿着道路、河流、田埂的边界成双线分布而且集中在线状地物上。原因两期影像虽然大致对齐但存在0.51个像素的系统偏移。线状地物只要偏移半个像素边缘灰度差就能轻松超过阈值。这在MAD里尤其隐蔽——MAD分量拖尾上的像素大部分来自这种边缘错位。解决进变化检测之前用控制点做一个自动配准精度检查。至少保证均方根误差低于0.5个像素。我对Landsat数据会额外做一次影像互相关配准用前期影像的某块区域做模板在后期影像对应区域搜索把残余偏移降到0.2像素以内再跑IR-MAD。5.2 辐射归一化没做透伪变化比真变化还多现象变化强度图上是整片整片的低频分布块跟地貌的阴坡阳坡高度相关而不是地物边界。原因两期影像的大气条件差异没有消除太阳高度角差异导致的整体辐射差被CVA和PCA当成变化。MAD系列对此免疫较好但如果你先用CVA算出了差值再交给PCA处理那误差就进PCA了。解决严格的生产流程应该是先做辐射定标和大气校正对MAD和IR-MAD可以跳过对CVA和PCA不能跳再做直方图匹配到参考影像。直方图匹配的参考影像我一般选云量少、时相接近、传感器一致性好的那一期。注意直方图匹配会改变光谱特征如果你后续要做变化方向分类方向角数据要在匹配前保存不然光谱信息就被破坏了。5.3 阈值不是玄学但Otsu不是万能的现象变化强度图出来后用Otsu自动阈值二值化结果把一大片水域全部标成变化而真正新增的建筑斑块反而没进去。原因Otsu基于双峰假设但真实的变化强度直方图往往是偏态单峰——绝大多数像素聚集在低强度区长尾延伸到高强度区。Otsu双峰假设不成立时它的分位点会落在偏高处小目标变化就丢了。解决这种分布用“均值若干倍标准差”或卡方分布的97.5百分位做阈值更稳。IR-MAD直接给卡方统计量用卡方分布找阈值是原则透明且可复现的CVA和PCA的强度值没有统计分布我一般取98分位或99分位数先看看哪些是高置信变化再用区域生长向外扩出完整地块。如果项目要求自动跑就把阈值设为97.5分位数人工再复核一遍避免一刀切。5.4 千万像素影像的内存爆炸现象处理一景GF-2影像约1.2万×1.2万像素IR-MAD的权重迭代里内存占用飙升最终程序被杀掉。原因IR-MAD每轮迭代要构造n×nn为像素数的中间矩阵例如(x - xm).T (w[:, None] * (y - ym))这一行如果用np.diag(w)来写直接生成n×n矩阵n是几千万时当然爆炸。即使在python里用广播写法如上面代码所示中间结果仍然需要容纳整个影像的加权差值。解决换成逐块处理。按瓦片划分影像比如每块256×256或512×512每块独立跑IR-MAD块与块之间保留20像素overlap最后用overlap区域的平均变化强度拼缝。overlap就是后悔药——如果拼缝出现明显跳变说明两块的协方差估计差异大回看该块是否有云或无效像素剔除后重跑。这个方案对生产级数据基本是必要的不要试图让整个协方差矩阵扛住全部数据。5.5 迭代不收敛时的“翻车”自查清单现象IR-MAD跑到第10次迭代权重图跟第一次比还在波动甚至出现了交替翻转。原因常规原因有两个。一是影像中存在大量异常值坏线、条纹、无效值协方差估计被污染权重在每个迭代重新计算后产生系统性偏移二是波段间相关性太高广义特征值求解的排序不稳定顺序一换权重就跟着换。解决先做严格的质量控制掩膜把所有已知的坏线、条纹、云、阴影、nodata全部掩掉。然后检查特征值排序是否在迭代中稳定打印每轮迭代前两个特征值如果某轮排序发生变化把掩膜外的像素比例降到5%以下通常就稳了。这个比例是我在实践中摸出来的——掩膜外的像素一旦超过20%IR-MAD基本就不收敛。6. 用模拟变化数据验证算法替换贴出的第一手经验最后这个技巧是我自己用的验证方法。真实影像做变化检测没有ground truth你永远不知道算法报出来的“变化”到底准不准。所以我的做法是随手找一景干净影像人工制造已知变化再跑算法看它能不能找回我埋的变化。构造方法很简单取一景影像X用本地编辑器把几个区域的值改掉比如把一个地块的波段值整体加10%再把另一个地块的波段全部换成另一景影像对应区域的值。这样我就有了“变化前X”和“变化后Y”而且知道每个变化区域的确切位置。跑MAD、IR-MAD、CVA、PCA四种算法对生成的强度图做阈值二值化和真值对比计算三个指标召回率检测出的变化像素占真值变化像素的比例。这个指标看漏检漏检多说明阈值太高或算法敏感性不足。误检率检出但真值没变化的像素占检出像素的比例。误检多说明伪变化压制不住。F1综合得分推荐在调参时用它当目标函数而不是只看召回率。我做的模拟实验里IR-MAD的F1通常比CVA高1015个百分点差距主要来自辐射差异干扰下的误检压制能力。但这只是我自己的数据你的数据不同结论可能翻转所以强烈建议把这个验证流程固化到你的流水线里。每次接新数据源先花半天做这个模拟验证后面处理全景数据心里就有底。这比直接拿真实数据磨半天却不知道哪里错要靠谱得多。希望帮到你。本文还有配套的精品资源点击获取 遥感变化检测必知：MAD、IR-MAD、CVA与PCA实战详解
简介面向遥感影像变化检测研究的经典算法代码包完整实现了IR-MAD、MAD、CVA、PCA四种方法均为Matlab可直接运行的脚本可解决多时相影像中地物变化识别、异常区域探测和光谱特征降维等问题适用于土地利用、环境监测、城市规划等研究方向需要读者具备基础遥感与矩阵运算知识。压缩包共97个文件、10.4MB以15个m源码为核心搭配28个bmp、24个tif测试影像以及hdr头文件和fig结果图源码中既包含IRMAD_Update、MADGet、CVADemo、PCADemo等主程序也有DAcom、covw、eigen2等辅助函数并附有泰州地区真实遥感样例数据可直接执行得到强度图和二值变化图。四种算法各有侧重IR-MAD和MAD对光照、大气影响更鲁棒CVA输出直观PCA擅于去噪和降维实际应用中可根据影像质量、变化类型和结果解释复杂度灵活选择包内另有Change Result Comparing.m等对比脚本方便输出定量比较结果。截止目前已有2008人学习适合遥感方向学生、科研人员及工程师快速上手经典算法并开展实验验证。1. 为什么“老三样”变化检测算法到今天还没过时遥感影像变化检测这几年被深度学习刷了很多屏但你去看实际生产的生态督查、违建监测、灾害评估项目跑在服务器上最稳的仍是MAD、IR-MAD、CVA、PCA这批经典算法。它们不吃海量训练样本不需要GPU你把两期影像丢进来半小时内就能出一张可交代的变化图。这篇文章就把这四个算法串起来讲清楚每个方法在做什么数学假设、代码怎么落地、参数怎么调、以及我实际踩过的那些坑。2. MAD与IR-MAD变化检测里的“差分增强器”2.1 MAD为什么比直接做差靠谱典型相关分析的直觉直接做差值有个根治不了的毛病——两期影像的辐射条件很难完全一致。传感器响应、太阳高度角、大气水汽任何一个变量变了差值图上就会出现大面积伪变化。MADMultivariate Alteration Detection的思路是不直接比较灰度而是先对两个多波段影像各做一次线性组合让组合后的分量之间相关性最强再去比较这两个组合分量的差。用数学语言描述就是找两组典型变量U a^T XV b^T Y最大化corr(U, V)。然后取M_i U_i - V_i作为MAD分量。为什么这样能压制辐射差异因为典型相关只保留了两个影像中“共同变化”的模式而共同变化以外的部分——比如传感器噪声、光照不均、大气影响——被踢到了残差里。理论上MAD分量服从卡方分布所以后续用卡方分布做显著性检验有天然的统计依据。这是它和直接差值的本质差别直接做差是像素级操作MAD是多元统计变换。我见过不少新手拿到两期Landsat影像上来就做波段差值结果变化图上全是耕地翻耕的条带和大气差异造成的噪声块。换成MAD之后这些噪声被自动过滤掉剩下的变化斑块才是真正需要关注的地表状态改变。2.2 最小可复现的MAD用numpy实现最简实现不需要遥感软件把两期影像的像元光谱当成两个随机向量就行。先把影像重排成像素数×波段数矩阵然后做典型相关分析import numpy as np def mad_transform(x, y, max_components3): # x, y: 像素数 x 波段数 的二维矩阵波段数通常3~8 # 返回mad分量矩阵每一列是一个mad分量 n_pixels x.shape[0] # 去均值典型相关分析要求零均值数据 xc x - x.mean(axis0) yc y - y.mean(axis0) # 计算协方差矩阵加正则项防止奇异 sxx np.cov(xc, rowvarFalse) np.eye(x.shape[1]) * 1e-6 syy np.cov(yc, rowvarFalse) np.eye(y.shape[1]) * 1e-6 sxy np.cov(xc, yc, rowvarFalse)[:x.shape[1], x.shape[1]:] syx sxy.T # 解广义特征值问题Sxy * Syy^-1 * Syx mat np.linalg.pinv(sxx) sxy np.linalg.pinv(syy) syx eigval, eigvec np.linalg.eig(mat) # 按特征值降序排序特征值越大相关性越强 idx np.argsort(eigval.real)[::-1] eigvec eigvec[:, idx] # 取前max_components个典型向量 a_vec eigvec[:, :max_components].real # 计算对应的b向量 b_vec np.linalg.pinv(syy) syx a_vec b_vec b_vec / np.sqrt((b_vec ** 2).sum(axis0)) # 计算典型变量并做差 u xc a_vec v yc b_vec mad_components u - v return mad_components这段代码的核心在广义特征值求解。mat矩阵的特征向量的意义是把x投影到某个方向后和y的某个方向相关性最大。注意我在sxx和syy对角线上加了1e-6这是防止协方差奇异。对实际遥感影像来说波段间相关性极强协方差矩阵很容易接近奇异不加这个正则项特征值求解会直接报错或者给出发散的结果。参数说明max_components一般取波段数或波段数减1。对4波段影像取3个MAD分量就够第四个分量对应的典型相关特征值趋近于1数学上的意义是“两期影像该方向上几乎无变化”直接舍弃。eigvec前几个分量的特征值应当接近1但又小于1如果你解出来的特征值等于1.0000000多半是正则化没加或者影像有大量无效值比如黑色边界没剔除。2.3 从MAD到IR-MAD迭代权重是怎么提精度的MAD有个明显的缺陷——它假设两期影像所有像素都是“没变化”的并在这个假设下估计典型相关结构。但真实场景里有变化的像素虽然不多可总有那么几个如果变化区域过大比如一场洪水淹了几万平方公里协方差矩阵会被真实变化污染典型相关计算出来的方向就不再是“共同稳定模式”了。IR-MADIteratively Reweighted MAD就是为这个设计的每次迭代根据当前MAD分量计算每个像素属于“no-change”的权重然后带权重新计算协方差矩阵再做一次MAD。变化像素的权重会越来越小no-change像素的主导地位就被加强了。按文献惯例权重按下式计算w_i P(卡方统计量 z_i)其中z_i sum((mad_k / sigma_k)^2) 是各MAD分量归一化后的平方和卡方分布自由度取分量数。from scipy.stats import chi2 def ir_mad(x, y, max_iter10, tol1e-4): # x, y: 像素数 x 波段数, 必须用有效值掩膜剔除nodata n_bands x.shape[1] # 初始权重全为1即第一轮等同于普通MAD w np.ones(x.shape[0]) for i in range(max_iter): # 带权计算均值和协方差 xm (x * w[:, None]).sum(axis0) / w.sum() ym (y * w[:, None]).sum(axis0) / w.sum() xc (x - xm) * np.sqrt(w[:, None]) yc (y - ym) * np.sqrt(w[:, None]) sxx xc.T xc / w.sum() np.eye(n_bands) * 1e-8 syy yc.T yc / w.sum() np.eye(n_bands) * 1e-8 sxy ((x - xm).T (w[:, None] * (y - ym))) / w.sum() syx sxy.T # 解广义特征值问题与普通MAD相同 mat np.linalg.pinv(sxx) sxy np.linalg.pinv(syy) syx eigval, eigvec np.linalg.eig(mat) order np.argsort(eigval.real)[::-1] a_vec eigvec[:, order].real # 用全部特征向量计算所有MAD分量 u (x - xm) a_vec v (y - ym) (np.linalg.pinv(syy) syx a_vec) mad u - v # 对方差归一化每个分量尺度统一 sigma mad.std(axis0) 1e-12 z ((mad / sigma) ** 2).sum(axis1) # 卡方分布计算no-change权重 new_w 1 - chi2.cdf(z, dfn_bands) # 收敛判断 if np.abs(new_w - w).max() tol: w new_w break w new_w return mad, w这里有个细节MAD分量计算出来后要对每个分量做方差归一化因为典型相关只保证相关性最大不保证分量尺度一致。z实际上是一个卡方统计量在no-change假设下应服从自由度等于分量数的卡方分布所以用1 - chi2.cdf(z)作为变化概率的补数也就是no-change权重。权重小说明该像素大概率真的变了这个权重向量本身就是一张很好的变化概率图。迭代停止条件用max_iter和tol双保险。我一般设max_iter10实际三五次迭代权重就收敛了。如果你发现迭代到第十次权重还在明显变化先别急着加迭代次数回去查一下影像预处理大概率是两期影像配准上有系统偏移或者有一条波段没做辐射定标只做了快速大气校正。2.4 IR-MAD的关键参数与停止条件IR-MAD名义上参数少实际上有两个地方非常影响结果。第一个是协方差矩阵的奇异值处理。遥感波段数通常48但像素数动辄几千万协方差矩阵肯定能求逆。坑不在这坑在无效值。影像四周的黑边、云掩膜产生的0值如果不做掩膜直接参与统计协方差会被拉向零值典型相关方向就是错的。我处理高分一号、高分二号的经验是先做一次有效值掩膜波段值0且不是nodata掩膜后的像素才进IR-MAD。第二个是迭代里的正则项系数。上面代码里sxx加的是1e-8这个值看数据情况。如果你发现前后两次迭代的权重分布出现振荡——第一次权重集中在A区域第二次却跑到B区域那多半是1e-8太小导致广义特征值求解数值不稳定。这时把正则项调到1e-6甚至1e-5振荡基本能消掉。代价是权重图会变钝一点但对最终二值化影响有限。3. CVA与PCA两个“先做变换再定阈值”的经典流派3.1 CVA的变化向量方向比大小更值钱CVAChange Vector Analysis的思路比MAD直白得多把每个像素看成多维光谱空间里的一个点两期影像的同一像素构成两条光谱向量直接做差得到一个“变化向量”。变化向量的模长就是变化强度变化向量的方向就是变化类型。这个算法最值得称道的地方是它把“变没变”和“变成了什么”分开了。模长超过阈值的像素被判定为变化区域而方向向量可以用来区分变化类别植被变裸地方向角一般在某个范围水体变建设用地又在另一个范围。实际项目里我常拿CVA当第一道粗筛先出整景变化强度图然后按方向角把变化像素聚类成几类再做分类后处理。但CVA的缺陷同样明显——它对辐射归一化的要求极高。直接做差意味着任何大气条件差异都会以“变化”的形式出现在模长里。所以用CVA前两期影像必须做至少一次直方图匹配或相对辐射归一化。这跟MAD不同MAD在数学上自带对公共模式的过滤CVA完全裸奔。3.2 CVA的Python实现从波段差到变化强度import numpy as np def cva_change_detection(x1, x2): # x1, x2: 已配准的两期多光谱影像, 形状为 h x w x n_bands # 先做相对辐射归一化: 用线性回归把x2匹配到x1的辐射范围 from sklearn.linear_model import LinearRegression h, w, nb x1.shape x1_flat x1.reshape(-1, nb).astype(np.float64) x2_flat x2.reshape(-1, nb).astype(np.float64) # 逐波段做线性归一化 x2_norm np.zeros_like(x2_flat) for b in range(nb): lr LinearRegression().fit(x2_flat[:, b:b1], x1_flat[:, b:b1]) x2_norm[:, b] lr.predict(x2_flat[:, b:b1]) diff x1_flat - x2_norm # 变化强度 欧氏距离, 变化方向 归一化的差值向量 magnitude np.sqrt((diff ** 2).sum(axis1)) direction diff / (magnitude[:, None] 1e-10) # 重排回影像形状 mag_img magnitude.reshape(h, w) dir_img direction.reshape(h, w, nb) return mag_img, dir_img这段代码最容易被忽略的是前面的线性归一化。很多初学CVA的帖子不写这一步直接用原始DN值做差结果在山区影像上几乎全花——因为地形阴影的辐射差远大于真实的地表变化。线性回归把x2整体调整到x1的辐射水平虽然不是严格意义上的大气校正但对付两期影像间的系统辐射偏移是管用的。参数说明LinearRegression逐波段做意思是一个波段一个增益和一个偏置。更高档的做法是考虑邻近像元光谱关系做相对辐射归一化但线性回归已经能消掉80%的系统辐射差异。如果你发现归一化后变化强度图还是噪点密布把影像先做一次3×3的均值滤波信噪比会明显改善——这是CVA最便宜的去噪手段。另外要注意direction的计算模长接近0的像素方向不稳定实际生产里我会把模长低于阈值的像素方向直接置为0避免后面的方向聚类被噪声主导。3.3 PCA做变化检测对差值影像降维取主要变化PCA进入变化检测的方式有点特殊不是对两期影像直接做主成分而是先做差值影像通常是CVA的差值然后对差值影像的各波段做主成分分析。PCA在这里的角色是“信息浓缩器”把n个波段的变化信息压缩到前几个主成分里第一主成分就是变化的主导方向方差贡献率越高说明变化越集中在一两种模式里。这个思路在很多项目里被用来区分真变化和噪声噪声在主成分上的分布通常均匀而真实的地表变化往往集中在第一、第二主成分上。所以操作上就是取前两个主成分的模长作为变化强度。打个比方这就跟PCA做特征脸类似——先用主成分提纯信号把分散在多维空间里的主要模式抽出来剩下的尾项当成噪声丢掉。对变化检测来说丢掉的那些尾项里确实藏着大量传感器噪声。def pca_change_detection(diff_flat, n_components2): # diff_flat: 像素数 x 波段数 的差值矩阵 # 去均值 d diff_flat - diff_flat.mean(axis0) # 计算协方差并做特征分解 cov np.cov(d, rowvarFalse) eigval, eigvec np.linalg.eigh(cov) # 特征值升序排列, 取最大的n_components个 idx np.argsort(eigval)[::-1][:n_components] proj d eigvec[:, idx] # 变化强度 前n_components个主成分的模长 change_mag np.sqrt((proj ** 2).sum(axis1)) return change_mag, eigval[idx]这里用np.linalg.eigh而不是np.linalg.eig因为协方差矩阵是对称阵eigh更稳定也更快。返回值里的eigval[idx]是前两个主成分的特征值它们占总特征值和的比率就是方差解释率。我看到很多人在报告里只写“PCA提取了主要变化”连解释率都不给——这个数字其实特别重要如果前两个主成分只解释了不到60%的方差说明变化模式很散PCA路线不适合这个场景别硬用。3.4 PCA的代码与“主成分数”怎么定主成分数在变化检测里一般取2到3个就够。不是越多越好多取一个主成分多带进一部分噪声。判断依据是特征值贡献率曲线也就是“陡坡图”。理想情况下曲线在前两三个点上急速下降长尾平缓拐点处就是该截断的位置。如果曲线一直平缓没有明显拐点说明差值影像本身没有主导的变化方向这时PCA和直接取差值模长没有本质区别。实际使用中还有一个变体——把PCA和MAD结合先做MAD得到几个分量再对这些分量做PCA压缩成一维或两维的变化指标。这个做法适合波段多的数据比如Sentinel-2的10米波段MAD已经把辐射差异压制了PCA再降一次维可视化效果比直接看卡方统计量更直觉。我最近处理一个城市扩张项目就是这么干的MAD出4个分量PCA压缩到2维变化强度图比单用MAD的卡方值干净而且能隐约看出不同扩张方向在颜色上的区分。4. 四个算法参数怎么设一张表和四类场景建议4.1 参数对比表MAD/IR-MAD/CVA/PCA算法核心步骤主要参数对辐射归一化的要求输出计算开销CVA逐像素光谱差阈值方向聚类角度极高必须先归一化变化强度图方向图最低纯逐像素运算PCA差值影像去均值特征分解主成分数解释率阈值高归一化与否影响提取方向主成分变化图低一次特征分解MAD典型相关分析差分max_components正则项低共同模式被自动过滤MAD分量卡方统计量中等需构造协方差矩阵IR-MAD带权迭代典型相关迭代次数正则项收敛容差低迭代中动态调整权重MAD分量权重图高多次协方差与特征分解这张表的第一列和第二列值得细读。CVA看着最“朴素”但对预处理的依赖最重MAD系列数学上最优雅但代码和调参的细节最多。PCA居中它能不能用取决于你差值影像的质量而差值影像的质量又取决于CVA那步归一化做得好不好。所以实际项目里这四个算法不是互相替代的关系而是流水线里的不同环节先归一化再做差值或MAD再用PCA或CVA出强度图最后二值化。4.2 不同数据场景的选型建议场景一只有两期中低分辨率多光谱Landsat级30m。优先IR-MAD。这类数据辐射一致性较好但云和阴影的多时相差异明显IR-MAD的权重迭代能在迭代过程中逐渐压掉云影伪变化。我做过一个2000到2010年Landsat的城市扩张监测IR-MAD比MAD直接出的变化斑块干净很多尤其是山区阴影带。场景二两期高分辨率影像0.52m有云阴影和建筑物阴影。CVA不乐观PCA也不乐观——高分辨率影像的地物阴影变化太剧烈。正确顺序是先做阴影掩膜掩膜外区域再做IR-MAD。如果实在想用CVA必须配合逐块本地阈值而不是全局阈值不然阴影边缘全是伪变化。场景三多光谱波段特别多8个以上比如Sentinel-2部分波段或WorldView-2。MAD系列的问题在于波段间相关性会迅速变强广义特征值求解的数值稳定性会出问题。务必把正则项加大到1e-5量级且只取前46个波段参与计算宁可用波段子集也不要一股脑把12个波段全塞进去。场景四只有单波段雷达数据或高程数据。MAD和CVA都需要多波段单波段可以用CVA的退化版——直接做差取绝对值或者用时空主成分分析把多期影像的时间维当成波段维做PCA。这个方法在地表沉降、水体面积监测里很常见。5. 变化检测常见坑从配准误差到阈值玄学5.1 影像没配准变化检测全白做现象变化检测结果图上一半检测出来的“变化”沿着道路、河流、田埂的边界成双线分布而且集中在线状地物上。原因两期影像虽然大致对齐但存在0.51个像素的系统偏移。线状地物只要偏移半个像素边缘灰度差就能轻松超过阈值。这在MAD里尤其隐蔽——MAD分量拖尾上的像素大部分来自这种边缘错位。解决进变化检测之前用控制点做一个自动配准精度检查。至少保证均方根误差低于0.5个像素。我对Landsat数据会额外做一次影像互相关配准用前期影像的某块区域做模板在后期影像对应区域搜索把残余偏移降到0.2像素以内再跑IR-MAD。5.2 辐射归一化没做透伪变化比真变化还多现象变化强度图上是整片整片的低频分布块跟地貌的阴坡阳坡高度相关而不是地物边界。原因两期影像的大气条件差异没有消除太阳高度角差异导致的整体辐射差被CVA和PCA当成变化。MAD系列对此免疫较好但如果你先用CVA算出了差值再交给PCA处理那误差就进PCA了。解决严格的生产流程应该是先做辐射定标和大气校正对MAD和IR-MAD可以跳过对CVA和PCA不能跳再做直方图匹配到参考影像。直方图匹配的参考影像我一般选云量少、时相接近、传感器一致性好的那一期。注意直方图匹配会改变光谱特征如果你后续要做变化方向分类方向角数据要在匹配前保存不然光谱信息就被破坏了。5.3 阈值不是玄学但Otsu不是万能的现象变化强度图出来后用Otsu自动阈值二值化结果把一大片水域全部标成变化而真正新增的建筑斑块反而没进去。原因Otsu基于双峰假设但真实的变化强度直方图往往是偏态单峰——绝大多数像素聚集在低强度区长尾延伸到高强度区。Otsu双峰假设不成立时它的分位点会落在偏高处小目标变化就丢了。解决这种分布用“均值若干倍标准差”或卡方分布的97.5百分位做阈值更稳。IR-MAD直接给卡方统计量用卡方分布找阈值是原则透明且可复现的CVA和PCA的强度值没有统计分布我一般取98分位或99分位数先看看哪些是高置信变化再用区域生长向外扩出完整地块。如果项目要求自动跑就把阈值设为97.5分位数人工再复核一遍避免一刀切。5.4 千万像素影像的内存爆炸现象处理一景GF-2影像约1.2万×1.2万像素IR-MAD的权重迭代里内存占用飙升最终程序被杀掉。原因IR-MAD每轮迭代要构造n×nn为像素数的中间矩阵例如(x - xm).T (w[:, None] * (y - ym))这一行如果用np.diag(w)来写直接生成n×n矩阵n是几千万时当然爆炸。即使在python里用广播写法如上面代码所示中间结果仍然需要容纳整个影像的加权差值。解决换成逐块处理。按瓦片划分影像比如每块256×256或512×512每块独立跑IR-MAD块与块之间保留20像素overlap最后用overlap区域的平均变化强度拼缝。overlap就是后悔药——如果拼缝出现明显跳变说明两块的协方差估计差异大回看该块是否有云或无效像素剔除后重跑。这个方案对生产级数据基本是必要的不要试图让整个协方差矩阵扛住全部数据。5.5 迭代不收敛时的“翻车”自查清单现象IR-MAD跑到第10次迭代权重图跟第一次比还在波动甚至出现了交替翻转。原因常规原因有两个。一是影像中存在大量异常值坏线、条纹、无效值协方差估计被污染权重在每个迭代重新计算后产生系统性偏移二是波段间相关性太高广义特征值求解的排序不稳定顺序一换权重就跟着换。解决先做严格的质量控制掩膜把所有已知的坏线、条纹、云、阴影、nodata全部掩掉。然后检查特征值排序是否在迭代中稳定打印每轮迭代前两个特征值如果某轮排序发生变化把掩膜外的像素比例降到5%以下通常就稳了。这个比例是我在实践中摸出来的——掩膜外的像素一旦超过20%IR-MAD基本就不收敛。6. 用模拟变化数据验证算法替换贴出的第一手经验最后这个技巧是我自己用的验证方法。真实影像做变化检测没有ground truth你永远不知道算法报出来的“变化”到底准不准。所以我的做法是随手找一景干净影像人工制造已知变化再跑算法看它能不能找回我埋的变化。构造方法很简单取一景影像X用本地编辑器把几个区域的值改掉比如把一个地块的波段值整体加10%再把另一个地块的波段全部换成另一景影像对应区域的值。这样我就有了“变化前X”和“变化后Y”而且知道每个变化区域的确切位置。跑MAD、IR-MAD、CVA、PCA四种算法对生成的强度图做阈值二值化和真值对比计算三个指标召回率检测出的变化像素占真值变化像素的比例。这个指标看漏检漏检多说明阈值太高或算法敏感性不足。误检率检出但真值没变化的像素占检出像素的比例。误检多说明伪变化压制不住。F1综合得分推荐在调参时用它当目标函数而不是只看召回率。我做的模拟实验里IR-MAD的F1通常比CVA高1015个百分点差距主要来自辐射差异干扰下的误检压制能力。但这只是我自己的数据你的数据不同结论可能翻转所以强烈建议把这个验证流程固化到你的流水线里。每次接新数据源先花半天做这个模拟验证后面处理全景数据心里就有底。这比直接拿真实数据磨半天却不知道哪里错要靠谱得多。希望帮到你。本文还有配套的精品资源点击获取