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