簡(jiǎn)介遙感影像變化檢測(cè)經(jīng)典算法資源包聚焦 IR-MAD、MAD、CVA、PCA 四種常用方法面向遙感、測(cè)繪、地理信息領(lǐng)域的科研人員和工程師幫助快速建立從多時(shí)相影像預(yù)處理到變化圖生成的完整流程。壓縮包共 97 個(gè)文件、約 10.4MB以 Matlab 腳本為主附有 bmp/tif 示例影像、ENVI hdr 頭文件、fig 結(jié)果圖其中 IRMAD_Update、MADGet、CVADemo、PCADemo 等算法腳本均可直接運(yùn)行便于對(duì)比各方法在光照、大氣條件差異下的表現(xiàn)。已有 2008 人學(xué)習(xí)下載適合系統(tǒng)了解經(jīng)典算法原理、在真實(shí)數(shù)據(jù)上對(duì)比檢測(cè)效果或?yàn)榫唧w課題篩選合適方法的研究者。借助泰州 TM 兩期影像示例可看到多時(shí)相異常檢測(cè)、變化向量分析、主成分降維和比值均值差分各自生成的中間變量與變化圖并可直接替換自己的遙感影像執(zhí)行實(shí)驗(yàn)節(jié)省從零實(shí)現(xiàn)的時(shí)間。1. 兩期影像擺在你面前變化檢測(cè)為什么繞不開這組老算法兩期同一區(qū)域的影像擺在桌上人工對(duì)比通常要小半天而 IR-MAD 這類經(jīng)典算法只需幾分鐘就能把變化斑塊圈出來。遙感影像變化檢測(cè)回答的問題很直接這兩個(gè)時(shí)相之間地表到底哪里變了、變成了什么。標(biāo)題里的四個(gè)算法——IR-MAD、MAD、CVA、PCA——是這一領(lǐng)域最經(jīng)典的組合它們不需要任何標(biāo)注樣本只靠?jī)善谟跋癖旧砭湍茌敵霾町悎D也正因?yàn)檫@個(gè)特性它們至今仍被當(dāng)作深度學(xué)習(xí)變化檢測(cè)方法的對(duì)照基線。適合做耕地監(jiān)測(cè)、違建排查、災(zāi)后評(píng)估的從業(yè)者。一個(gè)反直覺結(jié)論是PCA 這種看似做降維的算法在變化檢測(cè)里恰恰被用來找最有價(jià)值的差異方向而不是壓縮掉它們。2. 算法譜系先立住CVA、PCA 與 MAD 各自在解決什么這組算法經(jīng)常被放在一起比較但它們的定位完全不同。CVA 是最樸素的逐波段差值思路PCA 提供了一種線性變換視角而 MAD 則是為“不變關(guān)系”專門設(shè)計(jì)的統(tǒng)計(jì)建模。把它們拆開看清各自的前提假設(shè)才不會(huì)在實(shí)跑時(shí)用錯(cuò)方向。2.1 CVA逐波段做差后變化強(qiáng)度和方向怎么讀CVAChange Vector Analysis的操作非常簡(jiǎn)單把兩期影像對(duì)應(yīng)波段的灰度值相減得到一組差值向量。差值向量的模長(zhǎng)代表變化強(qiáng)度向量夾角代表變化類型。比如紅光波段從 120 變到 60、近紅外從 90 變到 150這個(gè)差值向量指向的方向就暗示植被可能發(fā)生了變綠或變枯。實(shí)跑時(shí)有兩個(gè)問題繞不開。第一是閾值怎么定常見做法是算所有像元差值模長(zhǎng)的均值與標(biāo)準(zhǔn)差取“均值加 k 倍標(biāo)準(zhǔn)差”作為變化閾值。k 通常在 1.5 到 3 之間具體取多少得看影像方差。第二是輻射差異的干擾如果兩期影像分別來自不同傳感器或不同季節(jié)大氣條件、太陽(yáng)高度角不一樣差值向量里會(huì)混入大量非地表變化成分。CVA 對(duì)這種全局性偏移幾乎沒有抵抗能力這是它最明顯的邊界。不過 CVA 的價(jià)值在于解釋性。它能把“變了”細(xì)化成“從什么方向變到什么方向”在土地覆蓋轉(zhuǎn)移分析里很受歡迎。我一般會(huì)用它做第一步粗篩后續(xù)再用其他算法復(fù)核。2.2 PCA主成分分析不是用來降維而是用來分離差異信息PCA 在變化檢測(cè)里有兩條常見路線。一條是先把兩期影像逐波段做差再對(duì)差值影像做主成分分析取前幾個(gè)主成分作為主要變化分量另一條是把兩期影像的所有波段堆疊成一個(gè)多波段影像一次性做主成分分析認(rèn)為前幾個(gè)主成分捕獲了兩期共有的穩(wěn)定背景信息排在后面的主成分則更多地體現(xiàn)時(shí)相差異。這里要理解一個(gè)關(guān)鍵點(diǎn)PCA 的各個(gè)主成分是沿著方差最大方向排列的。兩期影像里不變的地物通?；叶确植几叨认嚓P(guān)它們會(huì)集中體現(xiàn)在前幾個(gè)主成分里變化像元在統(tǒng)計(jì)上屬于少數(shù)派方差貢獻(xiàn)小反而落到后面的成分中。所以做變化檢測(cè)時(shí)常見的做法是丟掉前幾個(gè)主成分拿后面的成分來合成差異圖。這個(gè)思路和人臉識(shí)別里“特征臉”的典故一致——PCA 的基向量是被數(shù)據(jù)驅(qū)動(dòng)出來的特征方向只是變化檢測(cè)里我們關(guān)心的是那些方差小卻語(yǔ)義明確的尾巴。但 PCA 有一個(gè)結(jié)構(gòu)性弱點(diǎn)它是全局變換。整幅影像共享同一組特征向量局部區(qū)域的差異會(huì)被全局統(tǒng)計(jì)平均掉。如果變化區(qū)很小PCA 的效果就會(huì)變差。2.3 MAD用典型相關(guān)分析給“不變關(guān)系”建模MADMultivariate Alteration Detection的思路比 CVA 更進(jìn)一層。它不是直接在原始波段空間做差而是先對(duì)兩期影像分別做線性組合得到兩個(gè)“典型變量” U a?X、V b?Y然后計(jì)算它們的差值作為變化分量。為什么要繞這一圈因?yàn)?CVA 假設(shè)兩個(gè)時(shí)相同一波段的數(shù)值可以直接相減。但現(xiàn)實(shí)中兩期影像之間存在傳感器定標(biāo)差異、大氣路徑輻射差異直接相減會(huì)把系統(tǒng)性偏差當(dāng)變化。MAD 用典型相關(guān)分析CCA去尋找兩期影像之間最相關(guān)的線性組合找到之后再做差這樣得到的差值在統(tǒng)計(jì)意義上更接近真實(shí)的異常擾動(dòng)。數(shù)學(xué)上MAD 分量的方差等于 2(1?ρ?)其中 ρ? 是第 k 對(duì)典型變量之間的相關(guān)系數(shù)。相關(guān)系數(shù)越低說明這組線性組合在時(shí)序上越不一致也就是變化信息越突出。因此 MAD 分量通常按方差從大到小排列排在前面的分量包含最有辨識(shí)力的變化信號(hào)。相比 CVAMAD 對(duì)輻射偏移和波段間線性關(guān)系的干擾明顯更耐受這是它在實(shí)操中最受歡迎的原因。2.4 四個(gè)算法的適用邊界與選型對(duì)比算法輸入形式核心原理輸出結(jié)果適用場(chǎng)景CVA兩期多波段影像逐波段差值向量變化強(qiáng)度圖 變化方向角土地覆蓋轉(zhuǎn)移分析、快速粗篩PCA差值影像或堆疊影像協(xié)方差矩陣特征分解若干主成分中的差異分量全局變化模式探索、波段壓縮MAD兩期多波段影像典型相關(guān)分析 線性組合差若干 MAD 分量圖多時(shí)相輻射不一致時(shí)的穩(wěn)健檢測(cè)IR-MAD兩期多波段影像MAD 迭代加權(quán)逼近不變像元加權(quán)的 MAD 分量與變化概率圖高精度無監(jiān)督變化檢測(cè)、輻射歸一化選型建議很直接如果兩期影像來自同一傳感器、同季相、經(jīng)過嚴(yán)格輻射定標(biāo)CVA 就夠用如果影像來源復(fù)雜、輻射差異明顯直接用 MAD如果追求更高精度且有耐心調(diào)參數(shù)上 IR-MAD。PCA 更多是作為預(yù)處理或輔助手段出現(xiàn)很少單獨(dú)承擔(dān)最終判定。3. 跑通 MAD 最小實(shí)現(xiàn)數(shù)據(jù)準(zhǔn)備、核心推導(dǎo)與差異圖生成MAD 的原理并不復(fù)雜但真正跑通需要處理影像讀寫、矩陣展平、奇異值分解、差異圖重排這些環(huán)節(jié)。這一章給出一套可以用 GDAL NumPy SciPy 完整復(fù)現(xiàn)的最小實(shí)現(xiàn)不依賴 ArcGIS 或 ENVI 的現(xiàn)成工具箱。3.1 工具選型為什么用 GDAL NumPy 而不是 ArcGIS 工具箱ArcGIS 和 ENVI 里都有變化檢測(cè)工具但對(duì)生產(chǎn)流程有四個(gè)不友好之處一是批處理需要寫 Model Builder調(diào)試不直觀二是中間矩陣不透明出問題難以定位三是許可證環(huán)境對(duì)自動(dòng)化部署不友好四是很難把算法嵌入到自定義的后處理管線里。用 GDAL 負(fù)責(zé)影像 IONumPy 做數(shù)組運(yùn)算SciPy 做統(tǒng)計(jì)分布計(jì)算整個(gè)邏輯全都攤在代碼里每一行都能被檢查。此外公開數(shù)據(jù)集比如 Onera Satellite Change Detection 數(shù)據(jù)集都以 GeoTIFF 格式提供GDAL 是讀取這些數(shù)據(jù)最通用的方式。下面的示例默認(rèn)已安裝 GDAL、NumPy、SciPy 且兩期影像已配準(zhǔn)到同一網(wǎng)格。3.2 數(shù)據(jù)準(zhǔn)備對(duì)齊波段、掩膜無效值、展平樣本矩陣MAD 的輸入是兩期影像各自的像元矩陣。預(yù)處理的目標(biāo)有二一是把影像數(shù)組展平成“樣本×波段”的矩陣方便后續(xù)協(xié)方差計(jì)算二是把無值區(qū)、云區(qū)、邊界區(qū)排除掉避免污染統(tǒng)計(jì)量。from osgeo import gdal import numpy as np def read_to_matrix(path, valid_min0, valid_max65535): ds gdal.Open(path) arr ds.ReadAsArray() # shape: (band, row, col) ds None if arr.ndim ! 3: raise ValueError(需要多波段影像) arr arr.astype(np.float32) b, h, w arr.shape # 構(gòu)建有效像元掩膜所有波段都在有效灰度范圍且不是 NaN mask np.isfinite(arr[0]) for band in arr: mask (band valid_min) (band valid_max) # 展平成 (N, B) 矩陣N 為有效像元數(shù) samples arr[:, mask].T # shape: (N, B) return samples, mask, (h, w) X, mask1, _ read_to_matrix(t1.tif) Y, mask2, _ read_to_matrix(t2.tif) # 兩期影像有效像元交集 mask mask1 mask2 X X[mask] # 嚴(yán)格對(duì)齊后行列號(hào)一致直接用布爾掩膜篩選 Y Y[mask] print(有效像元數(shù):, X.shape[0], 波段數(shù):, X.shape[1])這段代碼里有三個(gè)關(guān)鍵點(diǎn)需要說明。ReadAsArray()返回的通道順序是 (band, row, col)不要和 OpenCV 的 HWC 順序搞混。astype(np.float32)是必需的原始影像常以 UInt16 存儲(chǔ)直接用整型求協(xié)方差時(shí)精度損失很大尤其遇到灰度值 0 和 1 附近的小數(shù)值時(shí)會(huì)嚴(yán)重失真。最后用有效像元交集統(tǒng)一兩期矩陣是為了保證下一步協(xié)方差計(jì)算基于同一批像元位置。3.3 核心推導(dǎo)標(biāo)準(zhǔn)化、SVD 與典型變量的關(guān)系MAD 的數(shù)學(xué)核心在于典型相關(guān)分析。這里用一個(gè)實(shí)現(xiàn)技巧先把兩期矩陣分別標(biāo)準(zhǔn)化為零均值、單位方差這樣自協(xié)方差矩陣變成單位矩陣交叉協(xié)方差的 SVD 結(jié)果就直接給出兩組投影方向。def mad_components(X, Y): # X, Y: (N, B)列對(duì)應(yīng)波段 n, b X.shape # 對(duì)每一波段做 z-score 標(biāo)準(zhǔn)化 x_mean X.mean(axis0) y_mean Y.mean(axis0) x_std X.std(axis0) y_std Y.std(axis0) Xn (X - x_mean) / (x_std 1e-8) # 避免除零 Yn (Y - y_mean) / (y_std 1e-8) # 交叉協(xié)方差矩陣 C (Xn.T Yn) / (n - 1) # SVD左右奇異向量就是 CCA 的投影方向 U, s, Vt np.linalg.svd(C) # 投影得到典型變量 u Xn A.T, v Yn B.T A U.T # shape (b, b) B Vt # 注意 Vt 已經(jīng)是轉(zhuǎn)置后的 V u Xn A.T # (n, b) 每列是一個(gè)典型變量 v Yn B.T # MAD 分量對(duì)應(yīng)列的差 mad u - v # 歸一化每個(gè) MAD 分量的理論方差是 2(1 - rho_k) for k in range(b): rho s[k] mad[:, k] mad[:, k] / np.sqrt(max(2 * (1 - rho), 1e-8)) return mad, s mad, rho mad_components(X, Y)這段代碼是把 CCA 求解簡(jiǎn)化為一次 SVD 的關(guān)鍵所在標(biāo)準(zhǔn)化讓兩個(gè)自協(xié)方差矩陣變成單位陣CCA 的廣義特征分解退化為普通 SVD。U的每一列是 X 側(cè)投影方向的轉(zhuǎn)置Vt的每一行是 Y 側(cè)投影方向。奇異值s[k]就是第 k 對(duì)典型變量的相關(guān)系數(shù)它越接近 1說明這對(duì)變量在兩期影像之間越一致對(duì)應(yīng) MAD 分量的方差越小、變化信息越弱。最后一步對(duì)每個(gè) MAD 分量除以sqrt(2(1 - rho))很微妙。這一步的作用是把分量方差統(tǒng)一歸一化到 1為下一步用卡方分布做統(tǒng)計(jì)檢驗(yàn)做準(zhǔn)備。如果你不做顯著性檢驗(yàn)而只是想肉眼看圖這步可以省略但后續(xù)要算變化概率時(shí)必須保留。3.4 差異圖生成MAD 分量平方和與閾值選取MAD 生成的是若干個(gè)差異分量?？梢园讯鄠€(gè)分量平方求和得到一個(gè)近似服從卡方分布的總統(tǒng)計(jì)量再用生存函數(shù)映射成每個(gè)像元的“變化概率”。這是 MAD 能直接產(chǎn)出概率圖的關(guān)鍵一步。from scipy import stats def mad_change_probability(mad, k_comp3): # 取前 k_comp 個(gè)方差最大的 MAD 分量 m mad[:, :k_comp] chi2_stat np.sum(m ** 2, axis1) # 卡方統(tǒng)計(jì)量 # 越小表示變化越顯著轉(zhuǎn)成“變化概率”便于可視化 p_nochange stats.chi2.sf(chi2_stat, dfk_comp) p_change 1.0 - p_nochange return p_change p_change mad_change_probability(mad) # 把概率向量擺回二維圖 h, w mask.shape change_map np.full((h, w), np.nan) change_map[mask] p_changek_comp取 2 或 3 是常見經(jīng)驗(yàn)值MAD 分量按方差降序排列前幾個(gè)分量集中了主要變化信息太靠后的分量基本是噪聲。卡方檢驗(yàn)的自由度等于分量數(shù)直接決定概率分布的形態(tài)。拿到了概率圖閾值怎么定先用一個(gè)簡(jiǎn)單試驗(yàn)在圖上疊加直方圖觀察概率分布是否呈雙峰。如果像元的概率集中在 0 附近和 1 附近閾值取兩峰之間的谷底即可如果分布連續(xù)那就得人為接受一個(gè)代價(jià)權(quán)衡比如把前 5% 概率最高的像元判為變化。這塊屬于經(jīng)驗(yàn)區(qū)不同影像的閾值差異很大后面避坑章會(huì)展開講。4. 從 MAD 到 IR-MAD迭代加權(quán)的實(shí)現(xiàn)與參數(shù)經(jīng)驗(yàn)MAD 一次性算完得到的結(jié)果往往不夠干凈那些真正的大面積變化區(qū)會(huì)像杠桿點(diǎn)一樣拉扯協(xié)方差矩陣導(dǎo)致投影方向被“帶偏”。IR-MADIteratively Reweighted MAD的補(bǔ)救思路很優(yōu)雅——每次迭代都重新估計(jì)哪些像元更可能沒有變化給這些像元更高權(quán)重再重新計(jì)算投影方向。4.1 為什么 MAD 會(huì)被強(qiáng)變化污染IR-MAD 怎么補(bǔ)救設(shè)想一個(gè)場(chǎng)景某塊區(qū)域發(fā)生了大范圍森林砍伐兩期影像在同一位置的光譜響應(yīng)差異非常大。MAD 在計(jì)算協(xié)方差時(shí)會(huì)給這些像元同等的統(tǒng)計(jì)地位于是協(xié)方差矩陣會(huì)被這批高強(qiáng)度差異牽著走最終算出來的投影方向不再是“不變背景下的最大差異”而是“背景差異與砍伐差異的混合體”。IR-MAD 的做法是給每個(gè)像元一個(gè)權(quán)重 w?權(quán)重代表這個(gè)像元屬于“未變化”的概率。第一輪所有像元等權(quán)等價(jià)于普通 MAD算出 MAD 分量后用卡方分布把每個(gè)像元的差異平方和映射成一個(gè)概率值差異越小的像元概率越高下一輪帶著這些概率重新計(jì)算加權(quán)協(xié)方差和投影方向反復(fù)迭代直到權(quán)重分布不再明顯變化或達(dá)到最大輪數(shù)。這個(gè)過程本質(zhì)上是把“不變像元”的辨識(shí)從一個(gè)一次性的全局假設(shè)變成一個(gè)逐步逼近的迭代估計(jì)。每一輪協(xié)方差矩陣都更偏向于不變背景投影方向也越來越精準(zhǔn)。這是 IR-MAD 比 MAD 精度高的根本原因。4.2 權(quán)重計(jì)算卡方分布尾概率與 mad 歸一化權(quán)重計(jì)算的輸入是當(dāng)前輪的 MAD 分量。每個(gè)像元有 B 個(gè)分量對(duì)分量平方求和得到一個(gè)標(biāo)量 q?。如果該像元真的未變化q? 應(yīng)服從自由度為 B 的卡方分布如果變化強(qiáng)烈q? 會(huì)落在分布的極端尾部。基于這個(gè)假設(shè)用卡方分布的生存函數(shù)把 q? 映射成概率w_i chi2.sf(q_i, dfB)q? 越小生存函數(shù)值越接近 1該像元被視為不變像元的可信度越高q? 極大時(shí)生存函數(shù)值趨近 0其統(tǒng)計(jì)權(quán)重也趨近 0。這輪概率就是下一輪的權(quán)重?！癿ad 歸一化”這個(gè)詞在 IR-MAD 的實(shí)現(xiàn)里通常指兩層操作一層是對(duì)輸入矩陣做加權(quán) z-score 標(biāo)準(zhǔn)化讓協(xié)方差計(jì)算在統(tǒng)一尺度下進(jìn)行另一層是對(duì)每個(gè) MAD 分量除以理論標(biāo)準(zhǔn)差使分量平方和方差匹配卡方分布的自由度。這兩層少了一層卡方檢驗(yàn)的統(tǒng)計(jì)性質(zhì)都會(huì)崩掉。當(dāng)年我第一次實(shí)現(xiàn)時(shí)漏掉分量的方差歸一化概率圖幾乎是一片純黑或純白就是這里出了問題。4.3 IR-MAD 完整迭代代碼與收斂判斷下面給出一個(gè)帶完整迭代邏輯的 IR-MAD 實(shí)現(xiàn)注意看權(quán)重如何在協(xié)方差計(jì)算和標(biāo)準(zhǔn)化兩個(gè)位置同時(shí)生效。def ir_mad(X, Y, n_iter10, tol1e-3): n, b X.shape w np.ones(n) / n # 初始等權(quán) mad_final None for it in range(n_iter): # 1. 用當(dāng)前權(quán)重做加權(quán)標(biāo)準(zhǔn)化mad歸一化 sum_w w.sum() x_mean (X * w[:, None]).sum(axis0) / sum_w y_mean (Y * w[:, None]).sum(axis0) / sum_w Xc X - x_mean Yc Y - y_mean # 2. 加權(quán)協(xié)方差 Cxx (Xc * w[:, None]).T Xc / sum_w Cyy (Yc * w[:, None]).T Yc / sum_w Cxy (Xc * w[:, None]).T Yc / sum_w # 3. 加權(quán) CCACholesky 白化后做 SVD Lx np.linalg.cholesky(Cxx 1e-6 * np.eye(b)) Ly np.linalg.cholesky(Cyy 1e-6 * np.eye(b)) Lx_inv np.linalg.inv(Lx) Ly_inv np.linalg.inv(Ly) K Lx_inv Cxy Ly_inv.T U, s, Vt np.linalg.svd(K) # 4. 投影方向變換回原坐標(biāo) A Lx_inv.T U # 列向量為 X 側(cè)投影方向 V Vt.T B Ly_inv.T V # 列向量為 Y 側(cè)投影方向 # 5. 計(jì)算 MAD 分量并做方差歸一化 u Xc A v Yc B mad u - v for k in range(b): mad[:, k] mad[:, k] / np.sqrt(max(2 * (1 - s[k]), 1e-8)) mad_final mad # 6. 更新權(quán)重卡方分布生存函數(shù) chi2_stat np.sum(mad[:, :b] ** 2, axis1) w_new stats.chi2.sf(chi2_stat, dfb) # 截?cái)鄻O小權(quán)重防止協(xié)方差計(jì)算時(shí)出現(xiàn)病態(tài) w_new np.clip(w_new, 1e-6, 1.0) # 7. 收斂判斷權(quán)重總體變化小于閾值 diff np.abs(w_new - w).sum() / n w w_new if diff tol: print(f第 {it1} 輪收斂, 權(quán)重變化量 {diff:.6f}) break return mad_final, w mad_ir, weight ir_mad(X, Y)收斂邏輯說明第 1 步的加權(quán)標(biāo)準(zhǔn)化確保均值估計(jì)不被強(qiáng)變化像元拉偏這是許多簡(jiǎn)化實(shí)現(xiàn)遺漏的地方。第 3 步加1e-6的正則項(xiàng)是為了防止當(dāng)某個(gè)波段在掩膜后樣本方差趨近 0 時(shí) Cholesky 分解崩潰。第 6 步用當(dāng)前輪 MAD 分量直接計(jì)算權(quán)重實(shí)現(xiàn)的是典型的 EM 式迭代先固定權(quán)重估計(jì)投影方向再固定投影方向更新權(quán)重。參數(shù)方面tol1e-3表示平均每個(gè)像元的權(quán)重相對(duì)變化小于千分之一即停止n_iter10是上限保護(hù)防止異常影像導(dǎo)致不收斂死循環(huán)。實(shí)際經(jīng)驗(yàn)中多數(shù)影像在第 4 到第 6 輪收斂。4.4 三個(gè)必調(diào)參數(shù)迭代次數(shù)、收斂閾值與權(quán)重截?cái)嗟螖?shù)最常見取 5 到 10。超過 10 輪后權(quán)重分布通常會(huì)進(jìn)入“不動(dòng)點(diǎn)”繼續(xù)迭代對(duì)結(jié)果幾乎沒影響。但如果影像中有大量強(qiáng)變化區(qū)前幾輪權(quán)重會(huì)把變化像元壓得極低導(dǎo)致協(xié)方差矩陣只反映純背景反而可能丟失真實(shí)而又微弱的漸變信號(hào)。因此我一般不建議一次拉滿 20 輪而是先用 5 輪看中間結(jié)果必要時(shí)再續(xù)跑。收斂閾值 tol 默認(rèn)1e-3即可。調(diào)大到這個(gè)值的 10 倍會(huì)加快速度但權(quán)重可能沒坐實(shí)調(diào)小可能出現(xiàn)輪次耗盡也不收斂的情況。權(quán)重截?cái)?e-6的作用是防止某些像元權(quán)重被更新成絕對(duì) 0從而徹底退出后續(xù)統(tǒng)計(jì)留下數(shù)值隱患。更保守的做法是截?cái)嗟?e-4這時(shí)邊緣像元對(duì)協(xié)方差仍有微小貢獻(xiàn)適合變化區(qū)特別破碎的影像。還有個(gè)容易被忽略的參數(shù)是卡方自由度。理論上自由度等于參與計(jì)算的波段數(shù) b但有些實(shí)現(xiàn)會(huì)只取前幾個(gè)獨(dú)立 MAD 分量進(jìn)入卡方統(tǒng)計(jì)此時(shí)自由度要改成實(shí)際取用的分量數(shù)。兩者混用會(huì)直接導(dǎo)致權(quán)重整體偏移概率圖閾值完全失真。5. 遙感變化檢測(cè)算法避坑五個(gè)讓結(jié)果翻車的常見問題這一章的價(jià)值直接來自踩坑現(xiàn)場(chǎng)。以下五個(gè)問題是我在多個(gè)項(xiàng)目里反復(fù)遇到的情況幾乎涵蓋了 MAD/IR-MAD 類算法最容易翻車的位置。5.1 兩期影像輻射不一致大面積假變化怎么排查現(xiàn)象結(jié)果圖上某塊連片區(qū)域被整體判為變化但實(shí)地核查發(fā)現(xiàn)地物根本沒變。原因兩期影像大氣條件、太陽(yáng)高度角或傳感器定標(biāo)參數(shù)不同導(dǎo)致同一地物在兩個(gè)時(shí)相的輻射值系統(tǒng)性偏移。IR-MAD 雖然對(duì)線性輻射偏移有抵抗力但如果偏移是非線性的比如大氣水汽對(duì)不同波段影響不同MAD 分量仍會(huì)殘留大量假信號(hào)。解決先做相對(duì)輻射歸一化。常見做法是選擇研究區(qū)內(nèi)的偽不變特征比如深水湖泊的深水區(qū)、大片裸巖、穩(wěn)定不透水面用它們擬合兩期影像的線性回歸關(guān)系再對(duì)第二期影像做校正。IR-MAD 迭代結(jié)束后的權(quán)重圖里權(quán)重接近 1 的像元恰恰可以作為偽不變特征的自動(dòng)候選這是一條閉環(huán)路線。5.2 協(xié)方差矩陣奇異波段數(shù)多于有效樣本數(shù)時(shí)怎么處理現(xiàn)象程序在np.linalg.cholesky或np.linalg.svd處報(bào)錯(cuò)提示矩陣不是正定。原因掩膜后有效像元數(shù)大于波段數(shù)很多時(shí)一般不會(huì)出錯(cuò)真正出錯(cuò)通常是某個(gè)波段動(dòng)態(tài)范圍幾乎為 0或者使用了超高光譜數(shù)據(jù)時(shí)波段間高度共線導(dǎo)致自協(xié)方差矩陣接近奇異。解決在 Cholesky 分解前給協(xié)方差矩陣加上一個(gè)小的正則項(xiàng)代碼里的Cxx 1e-6 * np.eye(b)就是干這個(gè)的。如果加了正則仍然報(bào)錯(cuò)優(yōu)先檢查波段選擇——把相關(guān)性極高的鄰接波段刪掉或用 PCA 預(yù)壓縮是更治本的做法。5.3 椒鹽噪聲嚴(yán)重變化圖里的孤立點(diǎn)怎么清理現(xiàn)象結(jié)果圖上散布大量孤立像元單看每一處都像變化但周圍背景完全穩(wěn)定。原因傳感器噪聲、配準(zhǔn)亞像元誤差、影像拉伸導(dǎo)致的小幅灰度抖動(dòng)都會(huì)在 MAD 分量的高維空間里表現(xiàn)為統(tǒng)計(jì)顯著的變化。IR-MAD 的權(quán)重迭代對(duì)這種高頻噪聲并不敏感因?yàn)樗鼈兊慕y(tǒng)計(jì)特征是空間不相關(guān)。解決在概率圖上做空間后處理最常見的手段是中值濾波、眾數(shù)濾波配合最小圖斑面積約束。給一個(gè)直觀經(jīng)驗(yàn)變化檢測(cè)成果圖斑面積小于 3×3 像元的圖斑在網(wǎng)格尺度下往往沒有制圖意義應(yīng)當(dāng)合并或剔除。5.4 閾值憑經(jīng)驗(yàn)拍腦袋精度虛高的原因與矯正現(xiàn)象取某個(gè)閾值后目視效果和驗(yàn)證點(diǎn)精確率的數(shù)字都很漂亮但模型換到相鄰區(qū)域后立刻失效。原因閾值是在單一影像上優(yōu)化的天然過擬合。MAD 輸出的概率分布形態(tài)在不同影像之間差異很大有的呈 U 形、有的呈單峰長(zhǎng)尾不存在一個(gè)通用閾值。解決用 Otsu 方法在總概率分布上自動(dòng)尋找分割點(diǎn)或使用兩成分高斯混合模型擬合概率分布以兩個(gè)高斯成分的交點(diǎn)為閾值。更穩(wěn)妥的做法是用已有 Ground Truth 上的變化比例反推閾值讓閾值對(duì)應(yīng)的變化面積與先驗(yàn)比例一致。這一條是精度評(píng)估里最后一道保險(xiǎn)也最常被忽略。5.5 配準(zhǔn)誤差導(dǎo)致邊緣條帶后處理掩膜怎么加現(xiàn)象房屋、道路、山脊線的邊緣出現(xiàn)沿地物輪廓的平行條帶看起來像地物“重影”。原因兩期影像配準(zhǔn)存在一兩個(gè)像元的平移誤差導(dǎo)致地物邊界兩側(cè)的灰度差被算法判定為變化。MAD 對(duì)這種高頻空間位移非常敏感這是統(tǒng)計(jì)方法普遍繞不開的結(jié)構(gòu)性缺陷。解決先查看配準(zhǔn)殘余誤差報(bào)告大于 0.5 像元的必須重做配準(zhǔn)后處理階段可以對(duì)變化概率圖進(jìn)行形態(tài)學(xué)開運(yùn)算去除細(xì)長(zhǎng)條帶再用邊緣掩膜把高梯度區(qū)域剔除掉防止線性地物邊緣被誤報(bào)。代價(jià)是真實(shí)沿邊界發(fā)生的小幅擴(kuò)展變化也會(huì)被誤刪這個(gè)取舍要和業(yè)務(wù)方確認(rèn)清楚。6. 結(jié)果驗(yàn)證與進(jìn)階技巧用混淆矩陣和 Kappa 評(píng)估變化檢測(cè)無監(jiān)督變化檢測(cè)算法最容易出現(xiàn)的幻覺是“看起來不錯(cuò)但精度不可證”。IR-MAD 輸出的概率圖必須經(jīng)過量化評(píng)估才能進(jìn)入生產(chǎn)流程。下面給出一個(gè)最小驗(yàn)證腳本搭配兩個(gè)我在實(shí)踐中沉淀下來的技巧。from sklearn.metrics import confusion_matrix, cohen_kappa_score # change_truth: 0/1 驗(yàn)證標(biāo)簽; p_change: 算法輸出的變化概率 threshold 0.7 change_pred (p_change threshold).astype(np.int16) # 全圖二值化 valid (change_truth 0) # 只評(píng)估有標(biāo)簽的像元 cm confusion_matrix(change_truth[valid], change_pred[valid], labels[0, 1]) tn, fp, fn, tp cm.ravel() oa (tn tp) / cm.sum() kappa cohen_kappa_score(change_truth[valid], change_pred[valid]) f1 2 * tp / (2 * tp fp fn) print(fOA{oa:.3f} Kappa{kappa:.3f} F1(change){f1:.3f})閾值threshold不是拍腦袋定的。我會(huì)先直方圖畫一遍p_change的分布如果雙峰明顯閾值放在兩峰谷底如果單峰就取排序后的 95 分位點(diǎn)再微調(diào)。Kappa 系數(shù)比總體精度 OA 更值得盯因?yàn)樵谧兓瘷z測(cè)里“無變化”通常占絕大多數(shù)一個(gè)永遠(yuǎn)全判無變化的模型也能拿到 90% 的 OAKappa 能把這個(gè)幻覺直接打到 0 附近。F1 則代表對(duì)變化類別的查準(zhǔn)與查全的折中業(yè)務(wù)上更符合實(shí)際訴求。另一個(gè)進(jìn)階技巧是把 IR-MAD 迭代結(jié)束時(shí)權(quán)重圖的倒數(shù)當(dāng)作“變化敏感度”的先驗(yàn)把它與概率圖做加權(quán)融合。權(quán)重越低的像元在迭代中被證明越是偏離不變背景的強(qiáng)變化點(diǎn)直接給這些位置的概率加一個(gè)小增量可以緩解閾值分割時(shí)低對(duì)比度變化被漏判的問題。我最早吃過一次虧某次水域淹沒了大片農(nóng)田由于水體信號(hào)強(qiáng)迭代權(quán)重把水邊漸變區(qū)的像元壓得很低概率閾值一卡整片漸淹區(qū)被漏判。后來我就養(yǎng)成了一個(gè)習(xí)慣——只看最終概率圖之前先單獨(dú)檢查 IR-MAD 每一輪權(quán)重圖的變化趨勢(shì)確認(rèn)是否出現(xiàn)局部區(qū)域權(quán)重陡降。實(shí)操上還有一條值得養(yǎng)成慣性在跑 IR-MAD 之前把兩期影像的直方圖和均值差打印出來差距超過一個(gè)標(biāo)準(zhǔn)差時(shí)先做相對(duì)輻射歸一化不要硬跑。這套做法陪我扛過了多次城市擴(kuò)張、森林?jǐn)_動(dòng)和災(zāi)后評(píng)估項(xiàng)目。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取