多元变化检测算法(Multivariate Alteration Detection，MAD)由A. A. Nielsen等人在论文Multivariate Alteration Detection (MAD) and MAF Postprocessing in Multispectral, Bitemporal Image Data: New Approaches to Change Detection Studies中首次提出，而后在论文The Regularized Iteratively Reweighted MAD Method for Change Detection in Multi- and Hyperspectral Data中，A. A. Nielsen提出了MAD的迭代版本IRMAD。
MAD基于典型相关分析CCA，是以投影特征差值方差最大化为准则提出的变化检测算法。
MAD算法的基础是典型相关分析(Canonical Correlation Analysis, CCA)，因此下面先介绍CCA。CCA是简单相关、多重相关的推广，是研究两组变量之间相关性的一种统计分析方法，CCA也是一种降维技术。CCA是将高维的两组数据分别降维到1维，然后用相关系数分析相关性。降维的原则是降到1维后，两组数据的相关系数最大。
{ C = [ X , Y ] T E ( C ) = [ μ X , μ Y ] T Σ = V a r ( C ) = [ Σ X X , Σ X Y Σ Y X , Σ Y Y ] \left\{\begin{matrix}C=\left[X, Y\right]^T \\ E(C) =\left[\mu_X, \mu_Y\right]^T \\ \Sigma= Var(C)=\left[\begin{matrix}\Sigma_{XX},\Sigma_{XY}\\\Sigma_{YX},\Sigma_{YY}\end{matrix}\right] \end{matrix}\right. ⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ C = [ X , Y ] T E ( C ) = [ μ X , μ Y ] T Σ = V a r ( C ) = [ Σ X X , Σ X Y Σ Y X , Σ Y Y ]
将X ,Y投影后得到的一维向量f X f_X f X f Y f_Y f Y
f X = a T X f Y = a T Y f_X=a^TX \\ f_Y=a^TY f X = a T X f Y = a T Y
f X f_X f X f Y f_Y f Y
{ v a r ( f X ) = a T Σ X X a v a r ( f Y ) = b T Σ Y Y b c o v ( f X , f Y ) = a T Σ X Y b \left\{\begin{matrix}var(f_X)=a^T\Sigma_{XX}a \\ var(f_Y)=b^T\Sigma_{YY}b \\ cov(f_X,f_Y)=a^T\Sigma_{XY} b \\ \end{matrix}\right. ⎩ ⎨ ⎧ v a r ( f X ) = a T Σ X X a v a r ( f Y ) = b T Σ Y Y b c o v ( f X , f Y ) = a T Σ X Y b
CCA的优化目标是最大化ρ ( f X , f Y ) \rho(f_X,f_Y) ρ ( f X , f Y ) a , b a,b a , b
a r g m a x a , b ρ = c o v ( f x , f y ) D ( f X ) D ( f Y ) \mathop {argmax} \limits_{a,b} \rho=\frac{cov(f_x,f_y)}{\sqrt {D(f_X)D(f_Y)}} a , b a r g m a x ρ = D ( f X ) D ( f Y ) c o v ( f x , f y )
由于分子分母增大相同的倍数，优化目标结果不变，可以采用和SVM类似的优化方法，固定分母，优化分子，优化问题转化为
a r g m a x a , b a T Σ X Y b s . t . { a T Σ X X a = 1 b T Σ Y Y b = 1 \mathop {argmax} \limits_{a,b} \ a^T\Sigma_{XY} b \\
s.t.\left\{\begin{matrix}a^T\Sigma_{XX} a=1
\\ b^T\Sigma_{YY} b=1
\end{matrix}\right. a , b a r g m a x a T Σ X Y b s . t . { a T Σ X X a = 1 b T Σ Y Y b = 1
求解方法是构造Lagrange等式
L = a T Σ X Y b − λ 2 ( a T Σ X X a ) − θ 2 ( b T Σ Y Y b ) L=\ a^T\Sigma_{XY} b-\frac{\lambda}{2}\left(a^T\Sigma_{XX} a \right)-\frac{\theta}{2}\left(b^T\Sigma_{YY} b \right) L = a T Σ X Y b − 2 λ ( a T Σ X X a ) − 2 θ ( b T Σ Y Y b )
L L L a a a b b b
Σ X X − 1 Σ X Y b = λ a Σ Y Y − 1 Σ Y X a = λ b \Sigma^{-1}_{XX}\Sigma^{}_{XY}b=\lambda a \ \ \ \Sigma^{-1}_{YY}\Sigma^{}_{YX}a=\lambda b Σ X X − 1 Σ X Y b = λ a Σ Y Y − 1 Σ Y X a = λ b
将第二个式子带入第一个式子可得
Σ X X − 1 Σ X Y Σ Y Y − 1 Σ Y X = λ 2 a \Sigma^{-1}_{XX}\Sigma^{}_{XY}\Sigma^{-1}_{YY}\Sigma^{}_{YX}=\lambda^2a Σ X X − 1 Σ X Y Σ Y Y − 1 Σ Y X = λ 2 a
对该式求特征值和特征向量即可得到λ \lambda λ a a a λ \lambda λ a a a b b b a a a b b b λ \lambda λ ρ ( f X , f Y ) \rho(f_X,f_Y) ρ ( f X , f Y )
假设有两景不同时相的多波段影像X X X Y Y Y p p p
a r g m a x a , b v a r ( a T X − b T Y ) s . t . { a T Σ X X a = 1 b T Σ Y Y b = 1 \mathop {argmax} \limits_{a,b} \ var \left(a^TX-b^TY \right) \\
s.t.\left\{\begin{matrix}a^T\Sigma_{XX} a=1
\\ b^T\Sigma_{YY} b=1
\end{matrix}\right. a , b a r g m a x v a r ( a T X − b T Y ) s . t . { a T Σ X X a = 1 b T Σ Y Y b = 1
根据约束条件我们可以得到v a r ( a T X − b T Y ) = 2 ( 1 − ρ ( a T X , b T Y ) ) var \left(a^TX-b^TY \right)=2 \left(1- \rho \left(a^TX,b^TY\right) \right) v a r ( a T X − b T Y ) = 2 ( 1 − ρ ( a T X , b T Y ) ) a a a b b b ρ \rho ρ ρ \rho ρ ρ \rho ρ a a a b b b M A D = a T x ^ − b T y ^ MAD=a^{T}\widehat{x}-b^{T}\widehat{y} M A D = a T x − b T y
由于MAD变量是X X X Y Y Y p p p
T j = ∑ i = 1 P ( M A D i j 2 σ M A D i ) ∈ χ 2 ( p ) T_{j}=\sum_{i=1}^{P}\left(\frac{MAD_{ij}^{2} }{\sigma_{MAD_i}}\right) \in \chi^{2}\left(p\right) T j = i = 1 ∑ P ( σ M A D i M A D i j 2 ) ∈ χ 2 ( p )
根据T j T_{j} T j
w j = P { T j > t } = P { χ 2 ( p ) > T j } w_{j}=P\left\{T_{j}>t \right\}=P\left\{\chi^{2}\left(p\right)>T_{j} \right\} w j = P { T j > t } = P { χ 2 ( p ) > T j }
然后让w j w_j w j
X i ‾ = ∑ j = 1 N w j X j i ∑ j = 1 N w j S k l = ∑ j = 1 N w j ( X j k − X k ‾ ) ( X j l − X l ‾ ) ( N − 1 ) ∑ j = 1 N w j / N \overline{X_{i}}=\frac{\sum_{j=1}^{N}w_{j}X_{ji}}{\sum_{j=1}^{N}w_{j}} \\
S_{kl}=\frac{\sum_{j=1}^{N}w_{j}\left( X_{jk}-\overline{X_{k}}\right) \left( X_{jl}-\overline{X_{l}}\right)}{(N-1)\sum_{j=1}^{N}w_{j}/N} X i = ∑ j = 1 N w j ∑ j = 1 N w j X j i S k l = ( N − 1 ) ∑ j = 1 N w j / N ∑ j = 1 N w j ( X j k − X k ) ( X j l − X l )
反复迭代直到典型相关系数收敛。
首先是MAD算法在泰州数据上的六个波段，如图3.1所示。
图3.1 MAD 6个波段
IRMAD在泰州数据上的六个波段。
图3.2 IRMAD 6个波段
图中非常亮的或者非常暗的区域代表变化、灰色区域代表未变化。可以看出在后两个波段，变化区域从未变化区域中被分离出来，前两个波段只包含少量变化信息。
由于T j T_{j} T j T j T_{j} T j
图3.3 binary change map(左边为MAD，右边为IRMAD)
相较于MAD，IRMAD通过迭代更新权值，得到的binary change map更精确，噪声更少，道路等变化区域被更好地从非变化区域中提取出来。
图3.4 四种方法的ROC曲线和AUC
从ROC曲线和计算得到的AUC可以看出MAD算法和SFA算法都是非常好的非监督变化检测算法，它们都可以达到很高的精度。
提出IRMAD的论文 The Regularized Iteratively Reweighted MAD Method for Change Detection in Multi and Hyperspectral Data https://ieeexplore.ieee.org/document/4060945 https://blog.csdn.net/weixin_41278720/article/details/80150098 https://blog.csdn.net/mbx8x9u/article/details/78824216
