企业官网建设流程全解析

简介：面向遥感影像变化检测研究的经典算法代码包，完整实现了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为什么比直接做差靠谱：典型相关分析的直觉

直接做差值有个根治不了的毛病——两期影像的辐射条件很难完全一致。传感器响应、太阳高度角、大气水汽，任何一个变量变了，差值图上就会出现大面积伪变化。MAD（Multivariate Alteration Detection）的思路是不直接比较灰度，而是先对两个多波段影像各做一次线性组合，让组合后的分量之间相关性最强，再去比较这两个组合分量的差。

用数学语言描述就是找两组典型变量：

U = a^T X，V = 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_components=3): # x, y: 像素数 x 波段数 的二维矩阵，波段数通常3~8 # 返回mad分量矩阵，每一列是一个mad分量 n_pixels = x.shape[0] # 去均值，典型相关分析要求零均值数据 xc = x - x.mean(axis=0) yc = y - y.mean(axis=0) # 计算协方差矩阵，加正则项防止奇异 sxx = np.cov(xc, rowvar=False) + np.eye(x.shape[1]) * 1e-6 syy = np.cov(yc, rowvar=False) + np.eye(y.shape[1]) * 1e-6 sxy = np.cov(xc, yc, rowvar=False)[: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(axis=0)) # 计算典型变量并做差 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-MAD（Iteratively 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_iter=10, tol=1e-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(axis=0) / w.sum() ym = (y * w[:, None]).sum(axis=0) / 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(axis=0) + 1e-12 z = ((mad / sigma) ** 2).sum(axis=1) # 卡方分布计算no-change权重 new_w = 1 - chi2.cdf(z, df=n_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_iter=10，实际三五次迭代权重就收敛了。如果你发现迭代到第十次权重还在明显变化，先别急着加迭代次数，回去查一下影像预处理：大概率是两期影像配准上有系统偏移，或者有一条波段没做辐射定标只做了快速大气校正。

2.4 IR-MAD的关键参数与停止条件

IR-MAD名义上参数少，实际上有两个地方非常影响结果。

第一个是协方差矩阵的奇异值处理。遥感波段数通常4～8，但像素数动辄几千万，协方差矩阵肯定能求逆。坑不在这，坑在无效值。影像四周的黑边、云掩膜产生的0值，如果不做掩膜直接参与统计，协方差会被拉向零值，典型相关方向就是错的。我处理高分一号、高分二号的经验是：先做一次有效值掩膜（波段值>0且不是nodata），掩膜后的像素才进IR-MAD。

第二个是迭代里的正则项系数。上面代码里sxx加的是1e-8，这个值看数据情况。如果你发现前后两次迭代的权重分布出现振荡——第一次权重集中在A区域，第二次却跑到B区域，那多半是1e-8太小导致广义特征值求解数值不稳定。这时把正则项调到1e-6甚至1e-5，振荡基本能消掉。代价是权重图会变钝一点，但对最终二值化影响有限。

3. CVA与PCA：两个“先做变换再定阈值”的经典流派

3.1 CVA的变化向量：方向比大小更值钱

CVA（Change 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:b+1], x1_flat[:, b:b+1]) x2_norm[:, b] = lr.predict(x2_flat[:, b:b+1]) diff = x1_flat - x2_norm # 变化强度 = 欧氏距离, 变化方向 = 归一化的差值向量 magnitude = np.sqrt((diff ** 2).sum(axis=1)) 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_components=2): # diff_flat: 像素数 x 波段数 的差值矩阵 # 去均值 d = diff_flat - diff_flat.mean(axis=0) # 计算协方差并做特征分解 cov = np.cov(d, rowvar=False) 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(axis=1)) 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看着最“朴素”，但对预处理的依赖最重；MAD系列数学上最优雅，但代码和调参的细节最多。PCA居中，它能不能用取决于你差值影像的质量，而差值影像的质量又取决于CVA那步归一化做得好不好。所以实际项目里这四个算法不是互相替代的关系，而是流水线里的不同环节：先归一化，再做差值或MAD，再用PCA或CVA出强度图，最后二值化。

4.2 不同数据场景的选型建议

场景一：只有两期中低分辨率多光谱（Landsat级，30m）。优先IR-MAD。这类数据辐射一致性较好，但云和阴影的多时相差异明显，IR-MAD的权重迭代能在迭代过程中逐渐压掉云影伪变化。我做过一个2000到2010年Landsat的城市扩张监测，IR-MAD比MAD直接出的变化斑块干净很多，尤其是山区阴影带。

场景二：两期高分辨率影像（0.5～2m），有云阴影和建筑物阴影。CVA不乐观，PCA也不乐观——高分辨率影像的地物阴影变化太剧烈。正确顺序是先做阴影掩膜，掩膜外区域再做IR-MAD。如果实在想用CVA，必须配合逐块本地阈值而不是全局阈值，不然阴影边缘全是伪变化。

场景三：多光谱波段特别多（8个以上，比如Sentinel-2部分波段或WorldView-2）。MAD系列的问题在于波段间相关性会迅速变强，广义特征值求解的数值稳定性会出问题。务必把正则项加大到1e-5量级，且只取前4～6个波段参与计算，宁可用波段子集也不要一股脑把12个波段全塞进去。

场景四：只有单波段雷达数据或高程数据。MAD和CVA都需要多波段，单波段可以用CVA的退化版——直接做差取绝对值，或者用时空主成分分析，把多期影像的时间维当成波段维做PCA。这个方法在地表沉降、水体面积监测里很常见。

5. 变化检测常见坑：从配准误差到阈值玄学

5.1 影像没配准，变化检测全白做

现象：变化检测结果图上一半检测出来的“变化”沿着道路、河流、田埂的边界成双线分布，而且集中在线状地物上。

原因：两期影像虽然大致对齐，但存在0.5～1个像素的系统偏移。线状地物只要偏移半个像素，边缘灰度差就能轻松超过阈值。这在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×n（n为像素数）的中间矩阵，例如(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高10～15个百分点，差距主要来自辐射差异干扰下的误检压制能力。但这只是我自己的数据，你的数据不同，结论可能翻转，所以强烈建议把这个验证流程固化到你的流水线里。每次接新数据源，先花半天做这个模拟验证，后面处理全景数据心里就有底。这比直接拿真实数据磨半天却不知道哪里错，要靠谱得多。希望帮到你。

本文还有配套的精品资源，点击获取