迭代加权多元变化检测(IRMAD) python实现

原论文标题：The Regularized Iteratively Reweighted MAD Method for Change Detection in Multi and Hyperspectral Data 

import numpy as np
from numpy.linalg import inv, eig # 求逆，求特征值
from scipy.stats import chi2 # 卡方分布

def irmad(X, Y, max_iter=150, epsilon=1e-6):
    assert(X.shape==Y.shape)
    c, n = X.shape

    # 正则化项，避免不同波段线性相关的情况，一般是对于高光谱影像
    # L = np.eye(c-2, c) - 2*np.eye(c-2, c, k=1) +np.eye(c-2, c,k=2)
    # Ohm = L.T @ L

    w = np.ones((1, n)) # 权重
    lambd_last = 100*np.ones((c, 1))
    for _iter in range(max_iter):
        # 加权平均
        X_mean = np.sum(w*X, axis=1, keepdims=True) / w.sum()
        Y_mean = np.sum(w*Y, axis=1, keepdims=True) / w.sum()

        # 加权协方差
        X_centered = X - X_mean
        Y_centered = Y - Y_mean
        V = np.concatenate([X_centered, Y_centered], axis=0)
        S = ( w*V @ V.T )/( (n-1) *w.sum()