搞過遙感的人對ENVI應(yīng)該都不陌生但能把這個軟件里的主成分分析PCA真正用明白的人其實不算多。我最早接觸PCA是在做多光譜影像分類的時候——9個波段一股腦扔進去分類精度反而比只用3個波段還差后來才明白問題就出在波段信息嚴(yán)重冗余上。從那時起PCA就成了我處理多光譜數(shù)據(jù)的默認(rèn)第一步。這篇教程按我自己的實操習(xí)慣來寫從原理講到ENVI里的具體操作再到紋理特征提取和后續(xù)踩坑記錄全文沒有廢話適合想真正把PCA用起來的同學(xué)無論是做分類、變化檢測還是影像融合都能從中找到可以直接照做的套路。1. 主成分分析到底在解決什么問題1.1 多光譜數(shù)據(jù)的“信息冗余”困境遙感影像和普通照片最大的區(qū)別是它不只有RGB三個波段。Landsat 8 OLI有9個波段Sentinel-2有13個波段高光譜甚至動輒上百個波段。波段多當(dāng)然信息量更大但很多波段之間高度相關(guān)比如紅波段和近紅外波段雖然數(shù)值差異很大但在植被覆蓋區(qū)域它們的變化趨勢幾乎同步這就是“信息冗余”。冗余帶來的直接后果是數(shù)據(jù)維度膨脹但有效信息沒有同比例增加分類算法會把重復(fù)的計算量和噪聲都吃進去輕則訓(xùn)練時間變長重則出現(xiàn)維度災(zāi)難模型越復(fù)雜精度反而越差。我見過很多新手一拿到影像就急急忙忙去做監(jiān)督分類結(jié)果分出十幾個類驗證精度卻不到60%問題很多時候就出在特征沒有提前做降維和去相關(guān)。PCA要解決的正是這個問題把多個相關(guān)波段通過線性變換壓縮成一組互不相關(guān)的新變量這些新變量叫主成分。第一個主成分承載原始數(shù)據(jù)中最大的方差信息第二個主成分承載剩余信息中最大的方差以此類推。實際操作下來通常前三個主成分就能扛起原始影像絕大部分信息剩下的則是壓縮后的噪聲和冗余。1.2 PCA的數(shù)學(xué)核心與直觀理解PCA的數(shù)學(xué)原理其實不復(fù)雜核心就是求協(xié)方差矩陣的特征值和特征向量。給定一個n維數(shù)據(jù)矩陣X先計算各波段之間的協(xié)方差矩陣C然后對C做特征分解求出特征值λ和對應(yīng)的特征向量v滿足關(guān)系式C v λ v特征值λ越大說明對應(yīng)的特征向量方向上數(shù)據(jù)方差越大也就是信息越多。每個特征向量其實就是一組權(quán)重系數(shù)主成分PC_i可以理解為原始所有波段的加權(quán)線性組合PC_i v1 * Band1 v2 * Band2 ... vn * Band_n用生活里的例子類比一個班50個學(xué)生每個人有語文、數(shù)學(xué)、英語三門成績這三門課高度相關(guān)成績好的通常三門都好。如果只允許用一個分?jǐn)?shù)給學(xué)生排名最好的辦法不是隨便挑一門課而是按“綜合得分”排名。這個綜合得分相當(dāng)于把三門課按各自權(quán)重加起來得到一個最能區(qū)分學(xué)生水平的新分?jǐn)?shù)——這其實就是PCA干的事情。在ENVI里PC1就是那個“綜合得分”它盡最大可能把數(shù)據(jù)間的差異集中到一個維度上。理解了這一層你就知道PCA為什么能做數(shù)據(jù)壓縮和噪聲抑制了既然是按方差從大到小排列主成分保留前幾個主成分就相當(dāng)于保住了“大信號”丟棄后面的主成分就相當(dāng)于丟掉了“小波動”這些波動往往就是噪聲或者波段間的隨機干擾。1.3 協(xié)方差矩陣和相關(guān)矩陣ENVI里怎么選ENVI的Forward PC Rotation運行時會讓你做一個選擇用協(xié)方差矩陣Covariance Matrix還是相關(guān)矩陣Correlation Matrix這個選擇題我見過很多人隨手就點了協(xié)方差矩陣其實里面有講究。協(xié)方差矩陣對波段本身的量綱和動態(tài)范圍非常敏感如果某個波段的數(shù)值范圍比其他波段大很多它就會在協(xié)方差矩陣中占主導(dǎo)地位算出來的主成分會過度偏向這個波段。相關(guān)矩陣本質(zhì)上是先對每個波段做了標(biāo)準(zhǔn)化處理讓所有波段處在同一個尺度上再算相關(guān)關(guān)系這樣每個波段對主成分的貢獻相對均衡。對多光譜影像來說如果各個波段之間的輻射定標(biāo)比較統(tǒng)一、數(shù)值范圍接近用協(xié)方差矩陣沒問題。但如果影像里混了熱紅外、短波紅外這類數(shù)值范圍差異很大的波段或者數(shù)據(jù)來自不同傳感器拼接我會建議用相關(guān)矩陣更穩(wěn)妥。高光譜數(shù)據(jù)做PCA時也建議優(yōu)先考慮相關(guān)矩陣因為波段間的量綱差異通常非常大。提示實際判斷方法很簡單——先看一眼各波段的統(tǒng)計值最大值、最小值、標(biāo)準(zhǔn)差的量級差別在三倍以內(nèi)用協(xié)方差矩陣沒有大問題明顯差出一兩個數(shù)量級的還是乖乖選相關(guān)矩陣吧。2. ENVI主成分分析完整操作流程2.1 數(shù)據(jù)準(zhǔn)備與軟件版本說明我用的是ENVI 5.6和5.7這兩個版本在Toolbox的菜單路徑上基本一致如果你還在用ENVI Classic經(jīng)典界面操作入口是Transform菜單下的Principal Components本質(zhì)上是一樣的算法只是入口和界面風(fēng)格不同。在操作前建議先確認(rèn)兩件事一是影像是否已經(jīng)做了輻射定標(biāo)和大氣校正PCA算的是波段之間的統(tǒng)計關(guān)系如果輸入數(shù)據(jù)本身有問題主成分結(jié)果也會帶著同樣的毛病二是如果影像存在明顯的無效值區(qū)域比如邊緣的黑邊、云和陰影最好先做一次掩膜處理不然這些像素會把協(xié)方差矩陣帶偏。打開文件的方式我也不啰嗦了File → Open As → Optical Sensor → Landsat Geometric或直接Open External File把影像先加載進來就可以開始操作。整個流程不需要任何第三方擴展ENVI自帶模塊就能完成。2.2 正向主成分旋轉(zhuǎn) Forward PC Rotation 的詳細(xì)操作正向主成分旋轉(zhuǎn)就是把原始波段變換成主成分序列操作路徑是Toolbox → Transform → Principal Components → Forward PC Rotation → Forward PC Rotation New Statistics and Rotate。在彈出的文件選擇框里選中你要處理的影像點擊OK后進入?yún)?shù)設(shè)置。這里有幾個關(guān)鍵參數(shù)要注意每個都有實際意義Stats Filename統(tǒng)計文件輸出路徑這個文件會記錄特征值和特征向量后邊Inverse反向旋轉(zhuǎn)和查看貢獻率時還要用到千萬別刪也別用中文路徑和中文文件名ENVI對中文路徑支持仍然不友好容易報錯。Spatial Subset只在影像某個子區(qū)域做統(tǒng)計并旋轉(zhuǎn)。如果你只想對研究區(qū)中心區(qū)域做處理或者需要剔除大量噪聲邊緣在這里框選范圍。Spectral Subset選擇參與計算的波段子集。比如Landsat 8有9個波段但不想把沿海氣溶膠波段和卷云波段放進來就可以在這里只選可見光到短波紅外的7個反射率波段。Covariance Matrix / Correlation Matrix前面提到的矩陣類型選擇。設(shè)置完成后點擊OKENVI會先計算統(tǒng)計信息再輸出一個多波段結(jié)果文件這個文件從PC1到PCn排列n就是輸入波段的個數(shù)。運行過程中如果數(shù)據(jù)量很大軟件界面可能會有幾秒到幾十秒的無響應(yīng)這是正常的不是卡死了。注意Output Result選項里默認(rèn)會生成一個臨時文件建議改成“Memory”或指定到本地磁盤路徑。如果數(shù)據(jù)量很大且內(nèi)存吃緊一定要存在磁盤上否則處理到一半內(nèi)存占滿整個ENVI都會崩潰我為此丟過好幾次沒保存的結(jié)果。2.3 如何讀懂PCA輸出的特征值表運行結(jié)束后很多人盯著生成的PC圖像不知道下一步該干嘛關(guān)鍵是要看懂那個.sta統(tǒng)計文件。你可以在文件管理器里用記事本打開也可以用ENVI的Layer Manager右鍵點擊結(jié)果文件查看Statistics但最直接的方式還是打開.sta文件內(nèi)容類似這樣Eigenvalues PC1 0.452317 82.343 82.343 PC2 0.061082 11.118 93.461 PC3 0.021553 3.922 97.383 PC4 0.008377 1.525 98.908 PC5 0.003648 0.664 99.572 PC6 0.001371 0.249 99.821 PC7 0.000982 0.179 100.000三列數(shù)字分別是特征值、單波段貢獻率百分比、累計貢獻率百分比。這是我手頭一個Landsat 8影像7個反射率波段的典型結(jié)果PC1貢獻率82.34%PC2貢獻率11.12%兩者累計已經(jīng)達到93.46%也就是說前兩個主成分就保留了原始7個波段超過93%的信息量。判斷保留多少個主成分我不建議死記“前三個”這種口訣而是看累計貢獻率。一般做分類累計貢獻率超過90%就可以了做數(shù)據(jù)壓縮存儲想盡量保留細(xì)節(jié)就取到95%以上做去噪反而可以適當(dāng)少留幾個主成分把后面的高頻噪聲直接扔在重建過程之外。還有一個需要留意的點特征值越大對應(yīng)的主成分圖像細(xì)節(jié)越豐富但并不是說后面那些貢獻率小的PC就毫無用處。在個別應(yīng)用中比如提取線性構(gòu)造、檢測地表異常信息這些低方差的PC往往會給出意想不到的線索因為它們?yōu)V掉了共性背景留下了特殊差異。2.4 反向主成分旋轉(zhuǎn) Inverse PC Rotation 的妙用反向旋轉(zhuǎn)的作用是從選定主成分中重建原始波段路徑是Toolbox → Transform → Principal Components → Inverse PC Rotation。這個操作看似冷門實際上非常實用。最典型的場景是基于PCA的影像去噪。處理流程是對原始影像做正向旋轉(zhuǎn)得到從PC1到PCn的序列把貢獻率很低的那些PC直接丟棄只選擇前面幾個高貢獻率PC作為輸入再執(zhí)行Inverse PC Rotation選擇正向旋轉(zhuǎn)時生成的.sta特征值文件ENVI就會用這幾個主成分的線性組合反算出一組新的波段圖像。這組重建出來的波段在視覺上和原始影像幾乎一樣但細(xì)節(jié)上的隨機噪聲明顯減少因為噪聲主要集中在那幾個被丟棄的低貢獻率PC里。我用這個方法處理過Sentinel-2影像再做后續(xù)分類整體精度比直接拿原始影像分類高出差不多3到5個百分點。另一個常見用途是數(shù)據(jù)壓縮存儲。如果原始影像有40個波段需要長期保存或者傳輸可以把正向旋轉(zhuǎn)后的結(jié)果只保留前8個PC輸出這樣存儲空間直接少了80%等到需要分析時再用反向旋轉(zhuǎn)重建。當(dāng)然這是有損壓縮對精度要求高的正式成果不建議長時間只保留壓縮版本至少要給自己留一份完整原始數(shù)據(jù)。3. 把PCA用在紋理特征提取上更香3.1 紋理特征與PCA有什么關(guān)系很多人提到PCA第一反應(yīng)是光譜降維但其實PCA在紋理特征提取上也是個神兵利器。紋理特征描述的是像素在空間上的灰度變化規(guī)律比如相干矩陣、反差、熵、同質(zhì)性、相異性等這些特征需要通過灰度共生矩陣GLCM來計算。問題在于GLCM紋理特征往往不止一個高分辨率影像或雷達影像提取出來動輒十幾個、幾十個紋理特征波段波段之間同樣存在嚴(yán)重的相關(guān)性。比如“均值”和“同質(zhì)性”在很多區(qū)域高度相關(guān)“對比度”和“相異性”也經(jīng)常聯(lián)動。這種情況下對紋理特征影像再做一次PCA效果立竿見影。PCA可以把幾十個紋理波段壓縮成少數(shù)幾個能夠衡量“紋理強度”“紋理復(fù)雜度”“紋理方向性”的綜合特征特征數(shù)量大幅減少但分類器拿到的紋理信息反而更純。我在做城市高分辨率影像分類時最常用的就是光譜波段PCA與紋理特征PCA的組合輸入。3.2 基于PCA的紋理特征提取實操步驟操作流程分四步每一步都有需要注意的參數(shù)細(xì)節(jié)。第一步計算GLCM紋理特征。打開影像后進入Toolbox → Texture → Co-occurrence Measures選擇需要計算紋理的波段。這里不是所有波段都要算選一個最具有代表性的波段往往效果最好比如近紅外波段對植被和建筑區(qū)分度就比紅波段更好。窗口大小建議選5x5或7x7窗口太小紋理噪聲大窗口太大又會平滑掉細(xì)節(jié)。第二步設(shè)置灰度量化級別Quantization Levels。這個參數(shù)控制灰度級數(shù)16級計算速度快但紋理細(xì)節(jié)損失明顯32級是均衡選擇大多數(shù)場景我都用它64級最精細(xì)但計算量和文件大小都直線上升小范圍研究可以用。步長Distances一般取1方向選All Directions。第三步生成紋理特征影像并做PCA。把計算出來的所有紋理特征波段合并成一個多波段文件然后按照第二部分的Forward PC Rotation流程對這個紋理特征文件做PCA。這一步的參數(shù)選擇和光譜PCA完全一致仍然要關(guān)注特征值表中的累計貢獻率。第四步選取紋理主成分參與后續(xù)建模。這一步我一般會做一個波段組合實驗把光譜主成分和紋理主成分放到一起再計算最佳指數(shù)因子OIF來挑選參與分類的最佳波段組合。通常紋理PCA的前兩個主成分就夠用加多了反而引入紋理噪聲。實操心得紋理PCA的PC1更多反映的是整體紋理強度比如建筑密集區(qū)和整齊農(nóng)田在PC1上往往差異巨大PC2則更多反映紋理的方向性和空間異質(zhì)性。當(dāng)你發(fā)現(xiàn)PC1和PC2區(qū)分度不夠時可以試試PC3甚至PC4不要急著否定紋理特征的有效性。3.3 PCA特征與原始波段的組合思路還有一種常見思路是把PCA壓縮后的特征和原始波段混合使用這在高分辨率影像分類里很流行。比如用WorldView-3做土地利用分類你可以保留原始4個多光譜波段的PC1和PC2再疊加紋理PCA的PC1構(gòu)成一個三維輸入特征空間。組合的關(guān)鍵問題是怎么判斷該保留哪些特征我自己的經(jīng)驗是分三步走第一步先觀察每個候選特征與已知地物類別之間的相關(guān)性計算各特征的類間距離第二步通過逐步判別或者隨機森林特征重要性排序篩選貢獻度高的特征第三步用OIF或J-M距離定量評估候選組合的可分性選得分最高的組合。這么說可能有點抽象簡單講不要一次性把所有特征全部塞進分類器先縮減到最多15到20個特征再靠特征重要性排序一步步淘汰。PCA在這里的價值就是保證你最終留下的特征相關(guān)性低、信息量高分類器跑得快還不容易過擬合。4. 常見問題排查與避坑指南4.1 前兩個主成分累計貢獻率偏低怎么辦我見過有人做完P(guān)CA后PC1貢獻率只有40%PC2也只有20%前兩個主成分加起來還不到70%原以為是軟件出了問題實際上多數(shù)情況下是這幾個原因。一是原始波段之間的相關(guān)性本身就很低。比如你輸入的數(shù)據(jù)里既有光學(xué)波段又有DEM、坡度等非遙感數(shù)據(jù)它們之間本來就沒有強相關(guān)性PCA自然擠不出一個主導(dǎo)性主成分。這種情況建議把數(shù)據(jù)按來源分組分別做PCA后再把主成分合起來不要強行混在一起。二是影像中存在大量無效值或異常像素。比如大范圍云覆蓋、水體表面太陽耀斑這些異常像元會干擾協(xié)方差統(tǒng)計。解決辦法是在運行Forward PC Rotation前先做一次像元篩選或?qū)τ跋褡鲅谀ぷ寘⑴c統(tǒng)計的像素更干凈。三是你在選擇輸入文件時混入了一個噪聲特別大的波段。比如某些熱紅外波段或受傳感器影響嚴(yán)重的波段它們本身方差很大但信息價值低擠占了主成分的權(quán)重。處理辦法是查看每個波段的直方圖和標(biāo)準(zhǔn)差把標(biāo)準(zhǔn)差異常大且分布發(fā)散的波段剔除后再試。4.2 ENVI里的SARscape工具包沒有GACOS怎么處理這個問題的出現(xiàn)頻率很高尤其是做InSAR時序分析的同學(xué)經(jīng)常會搜到“SARscape做大氣延遲校正需要GACOS數(shù)據(jù)”然后在ENVI的SARscape菜單里怎么翻都找不到GACOS相關(guān)的模塊心里就開始懷疑是不是自己的SARscape安裝不完整。先說結(jié)論SARscape菜單里沒有GACOS入口并不是軟件安裝問題而是GACOS數(shù)據(jù)需要通過在線服務(wù)單獨獲取SARscape本身只是一個數(shù)據(jù)處理框架它不會替你把這種外部氣象數(shù)據(jù)下載下來。GACOS的完整名稱是Generic Atmospheric Correction Online Service用于InSAR大氣延遲相位校正數(shù)據(jù)以網(wǎng)格文件形式提供給用戶。標(biāo)準(zhǔn)的處理流程是先到GACOS在線服務(wù)平臺注冊并申請覆蓋研究區(qū)域和對應(yīng)成像日期的數(shù)據(jù)文件下載后得到的是經(jīng)緯度網(wǎng)格格式的大氣延遲數(shù)據(jù)然后在SARscape的InSAR處理流程里找到與大氣校正相關(guān)的模塊通過讀取外部數(shù)據(jù)的方式把GACOS文件導(dǎo)入再進行相位校正。具體到Envisat或Sentinel-1數(shù)據(jù)的處理中需要先把GACOS文件轉(zhuǎn)換成SARscape能識別的格式這一步通常在SARscape的數(shù)據(jù)導(dǎo)入工具中完成。提示如果你在SARscape里確實找不到大氣校正或GACOS的相關(guān)子模塊先確認(rèn)自己安裝的是不是完整版SARscape模塊包括InSAR擴展而且不是所有版本和授權(quán)級別都開放了全部工具可以先查看Help里的模塊列表。數(shù)據(jù)下載請走官方申請渠道不要輕信網(wǎng)上打包好的第三方數(shù)據(jù)來源不明的數(shù)據(jù)質(zhì)量和時效性都沒保障。4.3 ENVI下載和安裝時容易踩的坑很多人搜“ENVI下載”是想找個免費包這個我只能給一個非常明確的建議ENVI作為商業(yè)軟件最好從官方渠道下載試用版或者通過所在單位、學(xué)校購買的正版授權(quán)來使用。網(wǎng)上那些來路不明的安裝包不僅可能帶病毒而且破解過程中經(jīng)常出現(xiàn)許可過期、模塊缺失反而更浪費時間。如果你已經(jīng)裝了正版但許可出現(xiàn)問題常見原因是許可服務(wù)器地址沒配對或者License過期。在ENVI啟動時會讀取許可配置建議檢查環(huán)境變量和許可文件路徑確認(rèn)服務(wù)器地址寫的是你單位許可服務(wù)器的IP而不是默認(rèn)的localhost。還有一個很常見的問題是安裝后Toolbox里某些工具是灰色的這說明當(dāng)前許可類型沒有包含對應(yīng)模塊比如SARscape模塊就是獨立的擴展授權(quán)ENVI基礎(chǔ)版裝好了也不能直接用。4.4 用Python復(fù)核ENVI的PCA結(jié)果最后分享一個我自己常用的交叉驗證方法拿Python的sklearn跑一遍同樣數(shù)據(jù)的PCA和ENVI的結(jié)果對比驗證操作有沒有出錯也方便做批量處理。下面這段代碼可以讀入ENVI導(dǎo)出的影像數(shù)據(jù)完成與ENVI幾乎相同的PCA計算并輸出各主成分的貢獻率。import numpy as np from osgeo import gdal from sklearn.decomposition import PCA # 讀取ENVI格式影像 ds gdal.Open(landsat8_subset.dat) arr ds.ReadAsArray() # shape: [波段數(shù), 行數(shù), 列數(shù)] rows, cols arr.shape[1], arr.shape[2] # 轉(zhuǎn)成二維每個像素一行每列是一個波段 data arr.reshape(arr.shape[0], -1).T # shape: [像素數(shù), 波段數(shù)] # 剔除無效值像素 data data[np.all(np.isfinite(data), axis1)] # 標(biāo)準(zhǔn)化到零均值等價于使用協(xié)方差矩陣 data_mean data - data.mean(axis0) # sklearn PCA pca PCA(n_componentsdata.shape[1]) scores pca.fit_transform(data_mean) # 輸出特征值和貢獻率 print(特征值方差:, pca.explained_variance_) print(貢獻率:, pca.explained_variance_ratio_) print(累計貢獻率:, np.cumsum(pca.explained_variance_ratio_)) # 查看第一主成分圖像 pc1 np.full((rows * cols, 1), np.nan) pc1[np.all(np.isfinite(arr.reshape(arr.shape[0], -1).T), axis1)] scores[:, 0] pc1_img pc1.reshape(rows, cols)對比ENVI輸出的.sta文件特征值兩者的差異應(yīng)該非常小一般在小數(shù)點后三位以內(nèi)。如果差得多優(yōu)先檢查預(yù)處理步驟是否一致比如是否做了標(biāo)準(zhǔn)化、是否排除了相同的無效像元。有一點要特別留意ENVI默認(rèn)用的是協(xié)方差矩陣sklearn的PCA也是基于協(xié)方差矩陣因為會先去中心化但如果數(shù)據(jù)量綱差異大ENVI里選了相關(guān)矩陣那Python這邊就要先用StandardScaler標(biāo)準(zhǔn)化數(shù)據(jù)再跑PCA不然兩邊的結(jié)果對不上。最后再分享一個使用技巧說回PCA本身我目前的固定習(xí)慣是拿到任何多光譜影像第一步先跑一次PCA看一眼特征值表花不了兩分鐘但能讓你對數(shù)據(jù)信息分布有個整體把握。如果PC1占比超過80%說明數(shù)據(jù)冗余度很高后續(xù)分類不用那么多波段如果PC1比較低說明各波段獨立性較強需要更謹(jǐn)慎地篩選特征。另外一個小技巧是PCA在影像融合中的應(yīng)用。很多人在做高分辨率全色影像和多光譜影像融合時只想到Brovey、GS變換其實把多光譜波段做PCA后用高分辨率全色波段替換PC1再反向旋轉(zhuǎn)回原始波段空間這種融合方式的色彩保真度在很多情況下優(yōu)于傳統(tǒng)方法值得一試。PCA是一個被講濫了但實際應(yīng)用仍然非常廣的工具希望這篇教程能幫你避開我踩過的那些坑少走點彎路。