ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

ENVI主成分分析实战:从原理到多光谱影像降维应用

ENVI主成分分析实战:从原理到多光谱影像降维应用 搞过遥感的人对ENVI应该都不陌生但能把这个软件里的主成分分析PCA真正用明白的人其实不算多。我最早接触PCA是在做多光谱影像分类的时候——9个波段一股脑扔进去分类精度反而比只用3个波段还差后来才明白问题就出在波段信息严重冗余上。从那时起PCA就成了我处理多光谱数据的默认第一步。这篇教程按我自己的实操习惯来写从原理讲到ENVI里的具体操作再到纹理特征提取和后续踩坑记录全文没有废话适合想真正把PCA用起来的同学无论是做分类、变化检测还是影像融合都能从中找到可以直接照做的套路。1. 主成分分析到底在解决什么问题1.1 多光谱数据的“信息冗余”困境遥感影像和普通照片最大的区别是它不只有RGB三个波段。Landsat 8 OLI有9个波段Sentinel-2有13个波段高光谱甚至动辄上百个波段。波段多当然信息量更大但很多波段之间高度相关比如红波段和近红外波段虽然数值差异很大但在植被覆盖区域它们的变化趋势几乎同步这就是“信息冗余”。冗余带来的直接后果是数据维度膨胀但有效信息没有同比例增加分类算法会把重复的计算量和噪声都吃进去轻则训练时间变长重则出现维度灾难模型越复杂精度反而越差。我见过很多新手一拿到影像就急急忙忙去做监督分类结果分出十几个类验证精度却不到60%问题很多时候就出在特征没有提前做降维和去相关。PCA要解决的正是这个问题把多个相关波段通过线性变换压缩成一组互不相关的新变量这些新变量叫主成分。第一个主成分承载原始数据中最大的方差信息第二个主成分承载剩余信息中最大的方差以此类推。实际操作下来通常前三个主成分就能扛起原始影像绝大部分信息剩下的则是压缩后的噪声和冗余。1.2 PCA的数学核心与直观理解PCA的数学原理其实不复杂核心就是求协方差矩阵的特征值和特征向量。给定一个n维数据矩阵X先计算各波段之间的协方差矩阵C然后对C做特征分解求出特征值λ和对应的特征向量v满足关系式C v λ v特征值λ越大说明对应的特征向量方向上数据方差越大也就是信息越多。每个特征向量其实就是一组权重系数主成分PC_i可以理解为原始所有波段的加权线性组合PC_i v1 * Band1 v2 * Band2 ... vn * Band_n用生活里的例子类比一个班50个学生每个人有语文、数学、英语三门成绩这三门课高度相关成绩好的通常三门都好。如果只允许用一个分数给学生排名最好的办法不是随便挑一门课而是按“综合得分”排名。这个综合得分相当于把三门课按各自权重加起来得到一个最能区分学生水平的新分数——这其实就是PCA干的事情。在ENVI里PC1就是那个“综合得分”它尽最大可能把数据间的差异集中到一个维度上。理解了这一层你就知道PCA为什么能做数据压缩和噪声抑制了既然是按方差从大到小排列主成分保留前几个主成分就相当于保住了“大信号”丢弃后面的主成分就相当于丢掉了“小波动”这些波动往往就是噪声或者波段间的随机干扰。1.3 协方差矩阵和相关矩阵ENVI里怎么选ENVI的Forward PC Rotation运行时会让你做一个选择用协方差矩阵Covariance Matrix还是相关矩阵Correlation Matrix这个选择题我见过很多人随手就点了协方差矩阵其实里面有讲究。协方差矩阵对波段本身的量纲和动态范围非常敏感如果某个波段的数值范围比其他波段大很多它就会在协方差矩阵中占主导地位算出来的主成分会过度偏向这个波段。相关矩阵本质上是先对每个波段做了标准化处理让所有波段处在同一个尺度上再算相关关系这样每个波段对主成分的贡献相对均衡。对多光谱影像来说如果各个波段之间的辐射定标比较统一、数值范围接近用协方差矩阵没问题。但如果影像里混了热红外、短波红外这类数值范围差异很大的波段或者数据来自不同传感器拼接我会建议用相关矩阵更稳妥。高光谱数据做PCA时也建议优先考虑相关矩阵因为波段间的量纲差异通常非常大。提示实际判断方法很简单——先看一眼各波段的统计值最大值、最小值、标准差的量级差别在三倍以内用协方差矩阵没有大问题明显差出一两个数量级的还是乖乖选相关矩阵吧。2. ENVI主成分分析完整操作流程2.1 数据准备与软件版本说明我用的是ENVI 5.6和5.7这两个版本在Toolbox的菜单路径上基本一致如果你还在用ENVI Classic经典界面操作入口是Transform菜单下的Principal Components本质上是一样的算法只是入口和界面风格不同。在操作前建议先确认两件事一是影像是否已经做了辐射定标和大气校正PCA算的是波段之间的统计关系如果输入数据本身有问题主成分结果也会带着同样的毛病二是如果影像存在明显的无效值区域比如边缘的黑边、云和阴影最好先做一次掩膜处理不然这些像素会把协方差矩阵带偏。打开文件的方式我也不啰嗦了File → Open As → Optical Sensor → Landsat Geometric或直接Open External File把影像先加载进来就可以开始操作。整个流程不需要任何第三方扩展ENVI自带模块就能完成。2.2 正向主成分旋转 Forward PC Rotation 的详细操作正向主成分旋转就是把原始波段变换成主成分序列操作路径是Toolbox → Transform → Principal Components → Forward PC Rotation → Forward PC Rotation New Statistics and Rotate。在弹出的文件选择框里选中你要处理的影像点击OK后进入参数设置。这里有几个关键参数要注意每个都有实际意义Stats Filename统计文件输出路径这个文件会记录特征值和特征向量后边Inverse反向旋转和查看贡献率时还要用到千万别删也别用中文路径和中文文件名ENVI对中文路径支持仍然不友好容易报错。Spatial Subset只在影像某个子区域做统计并旋转。如果你只想对研究区中心区域做处理或者需要剔除大量噪声边缘在这里框选范围。Spectral Subset选择参与计算的波段子集。比如Landsat 8有9个波段但不想把沿海气溶胶波段和卷云波段放进来就可以在这里只选可见光到短波红外的7个反射率波段。Covariance Matrix / Correlation Matrix前面提到的矩阵类型选择。设置完成后点击OKENVI会先计算统计信息再输出一个多波段结果文件这个文件从PC1到PCn排列n就是输入波段的个数。运行过程中如果数据量很大软件界面可能会有几秒到几十秒的无响应这是正常的不是卡死了。注意Output Result选项里默认会生成一个临时文件建议改成“Memory”或指定到本地磁盘路径。如果数据量很大且内存吃紧一定要存在磁盘上否则处理到一半内存占满整个ENVI都会崩溃我为此丢过好几次没保存的结果。2.3 如何读懂PCA输出的特征值表运行结束后很多人盯着生成的PC图像不知道下一步该干嘛关键是要看懂那个.sta统计文件。你可以在文件管理器里用记事本打开也可以用ENVI的Layer Manager右键点击结果文件查看Statistics但最直接的方式还是打开.sta文件内容类似这样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三列数字分别是特征值、单波段贡献率百分比、累计贡献率百分比。这是我手头一个Landsat 8影像7个反射率波段的典型结果PC1贡献率82.34%PC2贡献率11.12%两者累计已经达到93.46%也就是说前两个主成分就保留了原始7个波段超过93%的信息量。判断保留多少个主成分我不建议死记“前三个”这种口诀而是看累计贡献率。一般做分类累计贡献率超过90%就可以了做数据压缩存储想尽量保留细节就取到95%以上做去噪反而可以适当少留几个主成分把后面的高频噪声直接扔在重建过程之外。还有一个需要留意的点特征值越大对应的主成分图像细节越丰富但并不是说后面那些贡献率小的PC就毫无用处。在个别应用中比如提取线性构造、检测地表异常信息这些低方差的PC往往会给出意想不到的线索因为它们滤掉了共性背景留下了特殊差异。2.4 反向主成分旋转 Inverse PC Rotation 的妙用反向旋转的作用是从选定主成分中重建原始波段路径是Toolbox → Transform → Principal Components → Inverse PC Rotation。这个操作看似冷门实际上非常实用。最典型的场景是基于PCA的影像去噪。处理流程是对原始影像做正向旋转得到从PC1到PCn的序列把贡献率很低的那些PC直接丢弃只选择前面几个高贡献率PC作为输入再执行Inverse PC Rotation选择正向旋转时生成的.sta特征值文件ENVI就会用这几个主成分的线性组合反算出一组新的波段图像。这组重建出来的波段在视觉上和原始影像几乎一样但细节上的随机噪声明显减少因为噪声主要集中在那几个被丢弃的低贡献率PC里。我用这个方法处理过Sentinel-2影像再做后续分类整体精度比直接拿原始影像分类高出差不多3到5个百分点。另一个常见用途是数据压缩存储。如果原始影像有40个波段需要长期保存或者传输可以把正向旋转后的结果只保留前8个PC输出这样存储空间直接少了80%等到需要分析时再用反向旋转重建。当然这是有损压缩对精度要求高的正式成果不建议长时间只保留压缩版本至少要给自己留一份完整原始数据。3. 把PCA用在纹理特征提取上更香3.1 纹理特征与PCA有什么关系很多人提到PCA第一反应是光谱降维但其实PCA在纹理特征提取上也是个神兵利器。纹理特征描述的是像素在空间上的灰度变化规律比如相干矩阵、反差、熵、同质性、相异性等这些特征需要通过灰度共生矩阵GLCM来计算。问题在于GLCM纹理特征往往不止一个高分辨率影像或雷达影像提取出来动辄十几个、几十个纹理特征波段波段之间同样存在严重的相关性。比如“均值”和“同质性”在很多区域高度相关“对比度”和“相异性”也经常联动。这种情况下对纹理特征影像再做一次PCA效果立竿见影。PCA可以把几十个纹理波段压缩成少数几个能够衡量“纹理强度”“纹理复杂度”“纹理方向性”的综合特征特征数量大幅减少但分类器拿到的纹理信息反而更纯。我在做城市高分辨率影像分类时最常用的就是光谱波段PCA与纹理特征PCA的组合输入。3.2 基于PCA的纹理特征提取实操步骤操作流程分四步每一步都有需要注意的参数细节。第一步计算GLCM纹理特征。打开影像后进入Toolbox → Texture → Co-occurrence Measures选择需要计算纹理的波段。这里不是所有波段都要算选一个最具有代表性的波段往往效果最好比如近红外波段对植被和建筑区分度就比红波段更好。窗口大小建议选5x5或7x7窗口太小纹理噪声大窗口太大又会平滑掉细节。第二步设置灰度量化级别Quantization Levels。这个参数控制灰度级数16级计算速度快但纹理细节损失明显32级是均衡选择大多数场景我都用它64级最精细但计算量和文件大小都直线上升小范围研究可以用。步长Distances一般取1方向选All Directions。第三步生成纹理特征影像并做PCA。把计算出来的所有纹理特征波段合并成一个多波段文件然后按照第二部分的Forward PC Rotation流程对这个纹理特征文件做PCA。这一步的参数选择和光谱PCA完全一致仍然要关注特征值表中的累计贡献率。第四步选取纹理主成分参与后续建模。这一步我一般会做一个波段组合实验把光谱主成分和纹理主成分放到一起再计算最佳指数因子OIF来挑选参与分类的最佳波段组合。通常纹理PCA的前两个主成分就够用加多了反而引入纹理噪声。实操心得纹理PCA的PC1更多反映的是整体纹理强度比如建筑密集区和整齐农田在PC1上往往差异巨大PC2则更多反映纹理的方向性和空间异质性。当你发现PC1和PC2区分度不够时可以试试PC3甚至PC4不要急着否定纹理特征的有效性。3.3 PCA特征与原始波段的组合思路还有一种常见思路是把PCA压缩后的特征和原始波段混合使用这在高分辨率影像分类里很流行。比如用WorldView-3做土地利用分类你可以保留原始4个多光谱波段的PC1和PC2再叠加纹理PCA的PC1构成一个三维输入特征空间。组合的关键问题是怎么判断该保留哪些特征我自己的经验是分三步走第一步先观察每个候选特征与已知地物类别之间的相关性计算各特征的类间距离第二步通过逐步判别或者随机森林特征重要性排序筛选贡献度高的特征第三步用OIF或J-M距离定量评估候选组合的可分性选得分最高的组合。这么说可能有点抽象简单讲不要一次性把所有特征全部塞进分类器先缩减到最多15到20个特征再靠特征重要性排序一步步淘汰。PCA在这里的价值就是保证你最终留下的特征相关性低、信息量高分类器跑得快还不容易过拟合。4. 常见问题排查与避坑指南4.1 前两个主成分累计贡献率偏低怎么办我见过有人做完PCA后PC1贡献率只有40%PC2也只有20%前两个主成分加起来还不到70%原以为是软件出了问题实际上多数情况下是这几个原因。一是原始波段之间的相关性本身就很低。比如你输入的数据里既有光学波段又有DEM、坡度等非遥感数据它们之间本来就没有强相关性PCA自然挤不出一个主导性主成分。这种情况建议把数据按来源分组分别做PCA后再把主成分合起来不要强行混在一起。二是影像中存在大量无效值或异常像素。比如大范围云覆盖、水体表面太阳耀斑这些异常像元会干扰协方差统计。解决办法是在运行Forward PC Rotation前先做一次像元筛选或对影像做掩膜让参与统计的像素更干净。三是你在选择输入文件时混入了一个噪声特别大的波段。比如某些热红外波段或受传感器影响严重的波段它们本身方差很大但信息价值低挤占了主成分的权重。处理办法是查看每个波段的直方图和标准差把标准差异常大且分布发散的波段剔除后再试。4.2 ENVI里的SARscape工具包没有GACOS怎么处理这个问题的出现频率很高尤其是做InSAR时序分析的同学经常会搜到“SARscape做大气延迟校正需要GACOS数据”然后在ENVI的SARscape菜单里怎么翻都找不到GACOS相关的模块心里就开始怀疑是不是自己的SARscape安装不完整。先说结论SARscape菜单里没有GACOS入口并不是软件安装问题而是GACOS数据需要通过在线服务单独获取SARscape本身只是一个数据处理框架它不会替你把这种外部气象数据下载下来。GACOS的完整名称是Generic Atmospheric Correction Online Service用于InSAR大气延迟相位校正数据以网格文件形式提供给用户。标准的处理流程是先到GACOS在线服务平台注册并申请覆盖研究区域和对应成像日期的数据文件下载后得到的是经纬度网格格式的大气延迟数据然后在SARscape的InSAR处理流程里找到与大气校正相关的模块通过读取外部数据的方式把GACOS文件导入再进行相位校正。具体到Envisat或Sentinel-1数据的处理中需要先把GACOS文件转换成SARscape能识别的格式这一步通常在SARscape的数据导入工具中完成。提示如果你在SARscape里确实找不到大气校正或GACOS的相关子模块先确认自己安装的是不是完整版SARscape模块包括InSAR扩展而且不是所有版本和授权级别都开放了全部工具可以先查看Help里的模块列表。数据下载请走官方申请渠道不要轻信网上打包好的第三方数据来源不明的数据质量和时效性都没保障。4.3 ENVI下载和安装时容易踩的坑很多人搜“ENVI下载”是想找个免费包这个我只能给一个非常明确的建议ENVI作为商业软件最好从官方渠道下载试用版或者通过所在单位、学校购买的正版授权来使用。网上那些来路不明的安装包不仅可能带病毒而且破解过程中经常出现许可过期、模块缺失反而更浪费时间。如果你已经装了正版但许可出现问题常见原因是许可服务器地址没配对或者License过期。在ENVI启动时会读取许可配置建议检查环境变量和许可文件路径确认服务器地址写的是你单位许可服务器的IP而不是默认的localhost。还有一个很常见的问题是安装后Toolbox里某些工具是灰色的这说明当前许可类型没有包含对应模块比如SARscape模块就是独立的扩展授权ENVI基础版装好了也不能直接用。4.4 用Python复核ENVI的PCA结果最后分享一个我自己常用的交叉验证方法拿Python的sklearn跑一遍同样数据的PCA和ENVI的结果对比验证操作有没有出错也方便做批量处理。下面这段代码可以读入ENVI导出的影像数据完成与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: [波段数, 行数, 列数] rows, cols arr.shape[1], arr.shape[2] # 转成二维每个像素一行每列是一个波段 data arr.reshape(arr.shape[0], -1).T # shape: [像素数, 波段数] # 剔除无效值像素 data data[np.all(np.isfinite(data), axis1)] # 标准化到零均值等价于使用协方差矩阵 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文件特征值两者的差异应该非常小一般在小数点后三位以内。如果差得多优先检查预处理步骤是否一致比如是否做了标准化、是否排除了相同的无效像元。有一点要特别留意ENVI默认用的是协方差矩阵sklearn的PCA也是基于协方差矩阵因为会先去中心化但如果数据量纲差异大ENVI里选了相关矩阵那Python这边就要先用StandardScaler标准化数据再跑PCA不然两边的结果对不上。最后再分享一个使用技巧说回PCA本身我目前的固定习惯是拿到任何多光谱影像第一步先跑一次PCA看一眼特征值表花不了两分钟但能让你对数据信息分布有个整体把握。如果PC1占比超过80%说明数据冗余度很高后续分类不用那么多波段如果PC1比较低说明各波段独立性较强需要更谨慎地筛选特征。另外一个小技巧是PCA在影像融合中的应用。很多人在做高分辨率全色影像和多光谱影像融合时只想到Brovey、GS变换其实把多光谱波段做PCA后用高分辨率全色波段替换PC1再反向旋转回原始波段空间这种融合方式的色彩保真度在很多情况下优于传统方法值得一试。PCA是一个被讲滥了但实际应用仍然非常广的工具希望这篇教程能帮你避开我踩过的那些坑少走点弯路。
返回列表