做高分辨率全色圖與多光譜影像融合我最早嘗試的方法不是IHS也不是Brovey而是PCA。原因很現(xiàn)實當時手里只有Matlab基礎(chǔ)工具箱沒有專業(yè)遙感軟件PCA用自帶的cov、eig、reshape就能完整搭起來從讀取數(shù)據(jù)到輸出融合圖、算完評價指標一個腳本搞定。后來用這套流程處理了好幾組不同傳感器的影像反復(fù)調(diào)過參數(shù)、踩過不少坑今天把完整的實現(xiàn)思路、代碼邏輯和評價指標的計算方法一并整理出來。這篇內(nèi)容適合剛接觸圖像融合的研究生也適合工作需要快速出圖的工程師。如果你已經(jīng)會基本的Matlab操作但對“PCA融合到底在做什么”“為什么替換第一主成分”“評價指標應(yīng)該看哪幾個”這些問題還沒完全理清楚這篇文章可以直接幫你落地。1. 為什么選PCA做全色與多光譜融合1.1 全色圖和多光譜圖的“互補關(guān)系”遙感影像里全色圖PanchromaticPAN通常是單波段覆蓋的波長范圍很寬空間分辨率高拍出來的是一張灰度圖細節(jié)紋理清楚但沒有任何顏色信息。多光譜圖MultispectralMS正好相反通常包含紅、綠、藍、近紅外等多個波段每個波段的空間分辨率低但光譜信息豐富能還原地物的色彩差異。高分辨率全色圖與多光譜圖像融合目的就是把兩者的優(yōu)勢合并讓輸出影像既有多光譜的顏色判定能力又有全色圖的細節(jié)清晰度。這個需求在自然資源監(jiān)測、城市規(guī)劃、農(nóng)林調(diào)查里都很常見。比如一塊區(qū)域里多光譜圖能區(qū)分植被和水體但邊界模糊全色圖邊界非常清楚卻沒有類別顏色。融合之后邊界和顏色都有了后續(xù)再做分類、目標提取精度會明顯提升。1.2 主成分分析在這里到底做了什么PCA主成分分析的核心邏輯不復(fù)雜多光譜影像的幾個波段之間往往是高度相關(guān)的比如紅光和綠光在某些地物上變化趨勢接近。PCA把這幾個相關(guān)波段重新組合成一組互不相關(guān)的分量按方差從大到小排列第一主成分集中了最大的信息量本質(zhì)上就是多光譜影像共有的“亮度/結(jié)構(gòu)”成分。高分辨率全色圖恰恰也是以空間結(jié)構(gòu)信息為主的灰度圖。于是思路就順了把多光譜影像做PCA用全色圖去替換第一主成分再把替換后的主成分變換回原始波段空間。這樣全色圖的空間結(jié)構(gòu)信息被帶進了每個波段而其余主成分保留了多光譜特有的光譜差異融合結(jié)果就同時具備兩者的優(yōu)點。我第一次看到這個思路的時候覺得挺巧妙的但實現(xiàn)中有一個非常關(guān)鍵的細節(jié)全色圖不能直接替換PC1必須先做直方圖匹配。否則替換進去的灰度分布和原PC1差異太大逆變換后會出現(xiàn)整體偏色和亮度失真這個問題后面會專門講。1.3 與IHS、Brovey、小波方法對比PCA融合不是唯一方案甚至在某些場景下不是最優(yōu)方案但它的適用性和Matlab實現(xiàn)的便捷性很好。IHS融合把多光譜影像從RGB空間變換到亮度、色相、飽和度空間用全色圖替換亮度分量再反變換。優(yōu)點是速度快、實現(xiàn)簡單缺點是只適合三個波段而且光譜畸變比較大尤其在植被和陰影區(qū)域會出現(xiàn)明顯的顏色偏離。Brovey融合本質(zhì)是多光譜波段歸一化后乘以全色圖屬于比值融合。它能很好地保持亮度信息但色彩失真和噪聲放大問題也很明顯不適合做定量分析。小波融合把全色圖和多光譜圖都做小波分解在近似分量和細節(jié)分量上分別融合。光譜保真度最高但參數(shù)多、計算量大需要選擇小波基和分解層數(shù)對新手不友好。PCA融合在光譜保真度上介于IHS和Wavelet之間實現(xiàn)難度遠低于小波也不需要額外工具箱只要數(shù)據(jù)完成幾何配準Matlab幾行核心代碼就能跑通。這就是我推薦把它作為入門方案的原因。2. PCA圖像融合原理與Matlab實現(xiàn)細節(jié)2.1 PCA變換的數(shù)學(xué)基礎(chǔ)PCA變換本質(zhì)是一種線性變換把原始波段數(shù)據(jù)投影到新的正交空間中。假設(shè)多光譜影像有B個波段每個像素有B個灰度值構(gòu)成一個B維向量。所有像素的均值向量記為μ協(xié)方差矩陣記為C。PCA就是求C的特征值和特征向量特征值λ表示對應(yīng)特征向量方向上的方差反映該方向包含的信息量特征向量矩陣V的每一列代表一個主成分的投影方向。把原始像元向量減去均值后左乘V的轉(zhuǎn)置就得到主成分影像。第一主成分就是特征值最大對應(yīng)的那個方向上的投影它的方差最大信息含量最高。整個變換可以用一個矩陣乘法完成Matlab里不需要手動寫協(xié)方差分解直接用cov和eig就行。有一點容易被忽略eig函數(shù)返回的特征值默認按升序排列也就是最后一個才是最大特征值。如果不加處理直接取第一行當?shù)谝恢鞒煞謺玫降姆讲钭钚?、信息最少的分量整個融合結(jié)果就廢了。排序這步一定不能省。2.2 第一主成分替換的關(guān)鍵步驟把多光譜影像重采樣到與全色圖相同的行列數(shù)然后按像素矩陣形式排列(rows*cols, bands)的尺寸。對中心化后的數(shù)據(jù)求協(xié)方差矩陣計算特征向量和特征值按特征值降序?qū)μ卣飨蛄颗判蛴成涞街鞒煞挚臻g得到每個主成分的影像。第一主成分是一個二維矩陣尺寸與全色圖一致。替換不是簡單賦值而是先做直方圖匹配讓全色圖的灰度統(tǒng)計特性接近PC1再把匹配后的全色圖放回第一主成分的位置。最后逆變換這里用的逆變換矩陣是特征向量矩陣的轉(zhuǎn)置因為PCA變換矩陣是正交矩陣逆矩陣等于轉(zhuǎn)置。加上均值向量后reshape回多波段影像融合完成。整個過程里從“主成分空間”回到“原始波段空間”這一步的方向不要寫反。如果變換是pc (ms_centered * V)那么逆變換就是fused (pc_fused * V) ms_mean。我第一次實現(xiàn)的時候把這個轉(zhuǎn)置關(guān)系搞反了結(jié)果輸出圖像都是花屏排查了很久才發(fā)現(xiàn)問題。2.3 直方圖匹配為什么要做直方圖匹配這一步的官方名稱叫histogram matching目的是將全色圖灰度直方圖的形狀調(diào)整為與PC1一致。直方圖匹配的本質(zhì)是做一個非線性灰度映射讓兩張圖的灰度分布盡可能接近。為什么必須做因為PC1是多光譜波段共同亮度結(jié)構(gòu)的加權(quán)結(jié)果它的均值和方差反映了多光譜影像的輻射水平。全色圖雖然也是灰度圖但其灰度范圍、均值和方差跟PC1大概率不同直接把全色圖塞進PC1的位置逆變換后各波段的灰度值會被整體抬高或者壓低融合圖出現(xiàn)嚴重偏色光譜信息被破壞。用直方圖匹配之后全色圖的灰度分布被限制在與PC1相近的范圍里融合圖在提升空間細節(jié)的同時整體色調(diào)保持穩(wěn)定。Matlab里做直方圖匹配有兩種方式老版本的histeq(pan, imhist(pc1))和新版本的imhistmatch(pan, pc1)。后者在R2018a之后比較常用接口更直觀直接指定參考圖像。實測下來對于尺寸相同的兩張圖imhistmatch效果更穩(wěn)推薦優(yōu)先使用。2.4 核心代碼實現(xiàn)下面這段代碼是完整融合流程的Matlab實現(xiàn)我在多個數(shù)據(jù)集上驗證過按順序執(zhí)行即可。% 讀取數(shù)據(jù)ms為多光譜影像pan為全色影像 ms_img imread(multispectral.tif); pan_img imread(panchromatic.tif); % 多光譜重采樣到全色圖尺寸 [rows, cols] size(pan_img); ms_resized imresize(ms_img, [rows, cols]); % 轉(zhuǎn)double避免uint8運算精度丟失 ms_double im2double(ms_resized); pan_double im2double(pan_img); % 按像素矩陣排列每一行是一個像元每一列是一個波段 [rows, cols, bands] size(ms_double); ms_flat reshape(ms_double, rows*cols, bands); % 中心化 ms_mean mean(ms_flat); ms_centered ms_flat - ms_mean; % 協(xié)方差矩陣與特征分解 cov_mat cov(ms_centered); [V, D] eig(cov_mat); % 特征值降序排序防止取錯主成分 [lambda, idx] sort(diag(D), descend); V V(:, idx); % 投影到主成分空間 pc (ms_centered * V); % 每一行是一個主成分影像 % 提取第一主成分并還原成二維圖像 pc1 reshape(pc(1, :), rows, cols); % 直方圖匹配全色圖匹配到PC1的灰度分布 pan_matched imhistmatch(pan_double, pc1); % 替換第一主成分 pc_fused pc; pc_fused(1, :) reshape(pan_matched, 1, rows*cols); % PCA逆變換回到原始波段空間 fused_flat (pc_fused * V) ms_mean; % 還原為圖像尺寸 fused reshape(fused_flat, rows, cols, bands); fused im2uint8(fused);代碼里有一個細節(jié)值得注意im2double把圖像歸一化到0~1范圍能避免uint8計算溢出問題但最后要用im2uint8轉(zhuǎn)回原始位深否則保存成圖片時會丟失信息或變成全黑。如果原始數(shù)據(jù)是uint16最后要轉(zhuǎn)回uint16im2uint16即可。3. 融合評價指標怎么算、怎么解讀3.1 空間信息類指標融合后的影像必須回答兩個問題空間細節(jié)是否增強了光譜信息是否保住了針對前者常用指標有標準差、信息熵、平均梯度。標準差反映圖像的對比度值越大說明灰度分布越分散圖像層次越豐富。信息熵衡量圖像包含的信息量熵值越高代表細節(jié)越多、不確定性越大。平均梯度反映圖像的清晰度和紋理變化強度平均梯度大說明邊緣更銳利、細節(jié)更明顯。這三個指標都是單波段計算的多波段影像要分別算各波段再取平均。它們的共同特點是“數(shù)值越大越好”但注意它們只是單方面描述空間信息不能反映顏色是否正確所以必須結(jié)合光譜保真指標一起看。3.2 光譜保真類指標光譜保真的核心是“融合后的顏色與原始多光譜盡量一致”。最常用的指標是相關(guān)系數(shù)CC、光譜扭曲度、均方根誤差RMSE和相對全局維度綜合誤差ERGAS。相關(guān)系數(shù)計算融合影像與原始多光譜影像對應(yīng)波段的相關(guān)系數(shù)越接近1說明融合過程對光譜信息的破壞越小。光譜扭曲度先算融合影像與原始多光譜逐波段差的絕對值平均值再除以原始多光譜均值得到一個相對畸變比例。值越小越好。RMSE融合影像和原始多光譜影像逐像素差的均方根反映了整體輻射偏差。ERGAS遙感領(lǐng)域用得比較多的綜合指標它在RMSE的基礎(chǔ)上考慮了各波段的均值差異和分辨率比例數(shù)值越低代表整體融合質(zhì)量越好。這里有一個容易踩的坑光譜指標的計算基準不是原始未經(jīng)重采樣的低分辨率多光譜影像而是重采樣到全色分辨率后的多光譜影像。因為融合后的影像分辨率與全色圖一致你要保證對比的是同一空間分辨率下的圖像否則算出來的誤差包含了一部分重采樣造成的信息差異不能反映融合算法本身的光譜保真度。3.3 綜合評價指標怎么選實際寫論文或者做方案報告時不需要把所有指標都堆上去。我的建議是空間信息選信息熵、平均梯度、標準差里挑兩個光譜保真選相關(guān)系數(shù)、光譜扭曲度、RMSE里挑兩個這樣一組實驗下來就有四個指標足夠支撐結(jié)論。但如果只讓我用一個綜合指標我會選ERGAS因為它同時兼顧了多個波段和分辨率比例。對于Matlab實現(xiàn)需要先確定一個融合比例全色圖分辨率除以多光譜分辨率。這個比例只用于ERGAS計算數(shù)值上等于多光譜影像尺寸放大到全色尺寸的倍數(shù)。3.4 指標計算代碼% 假設(shè) fused 是融合結(jié)果uint8ms_resized 是重采樣后的原多光譜double fused_double im2double(fused); % 1. 標準差各波段平均 std_vals zeros(1, bands); for k 1:bands std_vals(k) std2(fused_double(:, :, k)); end mean_std mean(std_vals); % 2. 信息熵各波段平均 ent_vals zeros(1, bands); for k 1:bands p imhist(fused(:, :, k)) / numel(fused(:, :, k)); p(p 0) []; ent_vals(k) -sum(p .* log2(p)); end mean_ent mean(ent_vals); % 3. 平均梯度各波段平均 grad_vals zeros(1, bands); for k 1:bands [gx, gy] gradient(fused_double(:, :, k)); grad_vals(k) mean(sqrt(gx.^2 gy.^2), all); end mean_grad mean(grad_vals); % 4. 相關(guān)系數(shù)各波段與原始多光譜對應(yīng)波段 cc_vals zeros(1, bands); for k 1:bands cc_vals(k) corr2(fused_double(:, :, k), ms_double(:, :, k)); end mean_cc mean(cc_vals); % 5. 光譜扭曲度 spec_dist mean(abs(fused_double - ms_double), all) / mean(ms_double, all); % 6. 均方根誤差RMSE rmse_val sqrt(mean((fused_double - ms_double).^2, all));這些指標計算方式都比較常規(guī)重點是要統(tǒng)一數(shù)據(jù)類型fused和ms_double都必須是double類型取值范圍一致否則算出的誤差會異常偏大或偏小。另外各波段指標算完后取平均是通常做法報告中也可以把各波段指標單獨列出來這樣能看出哪些波段光譜失真更明顯。4. 完整實操從數(shù)據(jù)準備到結(jié)果分析4.1 數(shù)據(jù)準備實操之前先明確使用什么數(shù)據(jù)。如果只是驗證流程可以用Matlab自帶的pout.tif或cameraman.tif簡單模擬把彩色圖降采樣模擬低分辨率多光譜把灰度圖當全色圖。但更好的辦法是用公開的遙感數(shù)據(jù)集比如GeoEye、QuickBird、WorldView系列的樣例圖這些數(shù)據(jù)原尺寸就是多光譜與全色同時獲取融合結(jié)果更有說服力。數(shù)據(jù)準備好之后最耗時間的反而是預(yù)處理。首先要做幾何配準確保多光譜圖和全色圖在空間上完全對齊。如果兩張圖有偏移融合結(jié)果會出現(xiàn)重影和邊緣虛化指標再好也是假的。其次要檢查投影信息和行列分辨率記錄全色圖與多光譜圖的尺寸比例這對后續(xù)ERGAS計算很重要。4.2 融合過程在Matlab里我習(xí)慣把融合和指標計算分兩個腳本寫pca_fusion.m負責載入數(shù)據(jù)、重采樣、PCA融合、保存結(jié)果eval_metrics.m負責讀入原多光譜、全色圖和融合結(jié)果計算所有指標。pca_fusion.m里除了前面展示的核心代碼外還應(yīng)該把關(guān)鍵中間變量存出來比如PC1影像、直方圖匹配后的全色圖。這個習(xí)慣幫我避免過很多次調(diào)試困難如果融合結(jié)果異常可以先看PC1和匹配后的全色圖長什么樣判斷問題出在PCA環(huán)節(jié)還是替換環(huán)節(jié)。4.3 結(jié)果分析與指標解讀拿一組實驗數(shù)據(jù)舉例原始多光譜圖分辨率2m全色圖分辨率0.5m尺寸比例是4倍波段數(shù)為4R、G、B、NIR。重采樣多光譜到0.5m后PCA融合得到的影像在目視效果上邊緣銳利度明顯提升道路、建筑邊界不再模糊。此時看指標標準差和平均梯度比原多光譜提升明顯信息熵基本持平或略有提升說明空間細節(jié)增強有效相關(guān)系數(shù)在0.85到0.95之間光譜扭曲度在0.05到0.15之間說明光譜信息有一定損失但在可接受范圍。如果相關(guān)系數(shù)低于0.8或者光譜扭曲度超過0.2就要考慮是不是直方圖匹配環(huán)節(jié)出了問題或者原始影像配準精度不夠。目視檢查也不能少。指標再好如果融合圖出現(xiàn)明顯的偏色、重影或者邊緣震蕩說明算法流程里有bug。我一般會把原多光譜、PCA融合圖、全色圖三張圖放在一起對比看重點觀察植被區(qū)域和陰影區(qū)域這些區(qū)域最容易暴露光譜畸變。5. 排坑實錄與心得5.1 常見問題速查表問題現(xiàn)象可能原因解決辦法融合圖整體偏色嚴重直方圖匹配未做或匹配后分布差異大檢查imhistmatch是否生效對比PC1與匹配后全色圖直方圖結(jié)果圖出現(xiàn)花屏特征向量排序錯誤或逆變換方向不對檢查sort排序確認逆變換用V而非V空間細節(jié)提升不明顯全色圖被過度平滑或重采樣核設(shè)置不當嘗試用雙三次插值重采樣檢查全色圖本身是否清晰指標算出來數(shù)值異常大數(shù)據(jù)類型混用uint8與double直接做差統(tǒng)一轉(zhuǎn)換為double最后再轉(zhuǎn)回uint8輸出相關(guān)系數(shù)低到0.5以下圖像未配準或多光譜波段數(shù)太少先做精確配準確保參考影像合理融合圖有重影或邊界虛化多光譜與全色圖存在空間位移用控制點法或互相關(guān)法重新配準5.2 幾個容易踩的細節(jié)第一個細節(jié)多光譜影像和全色圖的行列數(shù)在很多數(shù)據(jù)里不是整數(shù)倍關(guān)系。比如多光譜是4500×4500全色圖是9000×9000比例正好是2倍這種情況比較順利但有的傳感器數(shù)據(jù)經(jīng)過處理全色圖尺寸不是嚴格整數(shù)倍此時重采樣時不能直接按imresize默認比例縮放要先明確目標尺寸是[size(pan,1), size(pan,2)]否則會出現(xiàn)輕微的對不齊。第二個細節(jié)PCA對波段間的相關(guān)性敏感。如果輸入的多光譜波段沒有做輻射歸一化比如各波段DN值范圍差異極大協(xié)方差矩陣會被高方差波段主導(dǎo)第一主成分主要是這個波段的信息融合效果也會受影響。穩(wěn)妥的做法是先把各波段歸一化到0~1或z-score標準化再進PCA。第三個細節(jié)直方圖匹配后的全色圖雖然灰度分布接近PC1但空間紋理信息是保留原樣的。如果原全色圖本身有傳感器噪聲匹配過程會把噪聲也放大并代入融合結(jié)果。遇到這種情況可以預(yù)先對全色圖做一個適度的高斯濾波或去噪處理但不要過度平滑否則細節(jié)增強效果就沒了。第四個細節(jié)關(guān)于特征向量符號。PCA分解出的特征向量方向并不是唯一的正負方向在不同版本、不同環(huán)境下可能不同。這不是bug因為最終結(jié)果經(jīng)過逆變換后會一致但如果你在中途直接查看主成分影像可能會發(fā)現(xiàn)某一次運行PC1的灰度方向是反的此時需要對比PC1和全色圖的明暗趨勢必要時乘以-1統(tǒng)一方向再使用。我在實際跑過幾組數(shù)據(jù)之后最深的體會是PCA融合的效果上限取決于輸入數(shù)據(jù)的預(yù)處理質(zhì)量。配準不準后面所有環(huán)節(jié)都白搭重采樣方式不同會直接影響指標數(shù)值。做實驗時盡量保證所有對比方案在相同預(yù)處理條件下運行這樣評價指標才公平。另外PCA融合雖然參數(shù)少但它畢竟是一種全局變換對于局部地物差異很大的場景融合結(jié)果可能會出現(xiàn)局部過亮或過暗這時候要考慮分塊處理或者換用小波融合。這套Matlab流程的優(yōu)勢在于改動成本低想對比方法只需要把PCA替換部分換成其他變換代碼就行。如果你手頭有現(xiàn)成的多光譜和全色數(shù)據(jù)建議直接拿一份出來跟著跑一遍。融合的結(jié)果可以先用目視檢查再算指標指標的計算代碼本身也可以當成一個小工具箱后續(xù)做其他融合方法時直接復(fù)用。這樣一套流程走下來既理解了PCA融合的原理也建立了一個可復(fù)用的融合評估框架。