简介这是一份基于R语言的高光谱数据处理工具包hsdar的完整源码与文档面向遥感、农业、环境科学等领域的科研人员及数据分析者解决高光谱数据从读取到建模的全流程问题。压缩包内含224个文件以R脚本76个、Rd帮助文档85个、Fortran源码25个为主并附有PDF说明、示例数据和可视化图片整体大小3.73MB目录结构规整便于离线安装、查阅源码与二次开发。资源覆盖hsdar包的全部核心功能支持通过rgdal接口读取ENVI、HDF、GeoTIFF等格式数据提供smooth平滑、atmcorr大气校正等预处理方法实现PCA、ICA、PLSR等特征提取算法集成SVM、随机森林分类及plsr回归建模工具并包含plotSpectra光谱曲线绘制、falseColor假彩色合成等可视化函数。通过研读这些源码与文档读者可深入理解高光谱处理算法细节快速搭建自己的分析流程。已有741人学习下载适合希望掌握R高光谱分析或扩展科研工具的中高级用户。1. R 高光谱数据分析怎么开始从一个能复现的流程说起真正做过一次高光谱数据的人都会承认这个方向的前半小时最容易让人放弃一架无人机飞完拿回来的是几百波段的二进制文件不是可以直接画图的图片想要做白板校正、暗电流扣除却发现手边的遥感软件要么授权过期要么菜单里的算法是个黑匣子。R 的开源生态其实早就把这条链路补齐了hsdar、hyperSpec、hsi这些包覆盖了从 ENVI 数据读取、反射率校正、平滑去噪到 PCA 降维、光谱角制图和混合像元分解的完整流程。这套做法能解决一个很具体的诉求在不依赖商业授权、每个步骤都可见可改的前提下把高光谱数据变成可解释的分类结果和定量参数。适合遥感、农业、地质、材料检测方向的工程师也适合正在做课程设计和论文复现的学生。2. 把高光谱读进 R 之前开源包选型与数据导入的第一次成功2.1 先想清楚你要处理的是“光谱”还是“立方体”R 里的高光谱开源包并不少但它们的定位差异很大。hyperSpec面向的是“一束一束的光谱”数据对象围绕光谱矩阵和波长展开适合处理地物光谱仪测出来的点光谱做聚类、判别和化学成分反演hsdar的数据对象叫speclib更像一个“光谱库”每个样本是一行、每个波段是一列同时对影像像元和点光谱都友好内置了SAM、光谱重采样和植被指数计算hsi则是专门面向“影像立方体”的包可以加载公开数据集、按像元提取光谱、做图像级可视化再往下还有terra和raster它们不是高光谱专用但负责最后一步的大栅格分块写入。选型不必纠结“哪个包最好”而是先回答两个问题数据是点光谱还是影像立方体最终要出分类图还是出参数表包数据对象擅长做的事不适合做的事hyperSpechyperSpec 对象点光谱处理、可视化和化学计量分析空间影像的整块读写hsdarspeclibENVI 读取、光谱重采样、SAM、植被指数超大型栅格的空间平滑hsiHSI 影像对象高光谱影像可视化、公开数据集加载细致的光谱预处理terra / rasterSpatRaster / RasterLayer大栅格分块、并行写入、空间滤波光谱角度计算我在实际项目里的一般取法是如果拿到的是机载或星载的 ENVI 标准格式影像立方体主轴用hsdar数据可视化用hsi辅助如果是对地物光谱仪导出的 CSV 点光谱直接走hyperSpec最后要做大规模制图把光谱矩阵算好后导给terra分块落盘。这套组合在 R 里只需要四个开源包不依赖任何商业授权。新手最容易踩的坑是“一个包想干所有事”结果装了半天发现它连自己常用的数据格式都不支持。2.2 安装与第一次成功读取ENVI 文件的正确打开方式安装这组包本身没什么难的真正让新手卡住的是hsdar在 Windows 上编译依赖。先把基础依赖装好代码顺序如下options(repos c(CRAN https://cloud.r-project.org)) install.packages(c(hsdar, hsi, prospectr, nnls))prospectr提供 Savitzky-Golay 平滑和标准正态变换给预处理备用nnls后面做混合像元分解要用hsdar负责读取 ENVI 格式hsi负责公开数据集加载。安装hsdar时如果提示缺少编译工具Windows 用户需要先安装 Rtools并在弹窗里把“Add to PATH”勾上这是最常见的安装失败原因。数据读取是最需要耐心的环节。一个标准的 ENVI 高光谱影像通常由.dat或.img和.hdr两个文件组成.hdr里写明了波段数、行数、列数、数据类型、存储顺序interleave和波长列表。只要两个文件同名且在同一个目录readSpeclib就能自动关联library(hsdar) # 假设拿到的文件是 scene.dat 和 scene.hdr img - readSpeclib(path/to/scene.dat, type ENVI) img读取成功后务必马上做三件事看波段数对不对、看像元数对不对、看波长有没有异常。nBands(img) # 影像波段数 nSpectra(img) # 像元总数 wavelength(img) # 每个波段的中心波长 summary(spectra(img))speclib对象的核心其实就是一个光谱矩阵加一个波长向量。矩阵的每一行是一个像元每一列是一个波段这个结构会让后面的 PCA、SAM 和回归都特别顺手。summary(spectra(img))能快速暴露数据质量问题比如某个波段的数值全为 NA或者最大最小值相差几个量级这些都是采集成像阶段留下的伤疤。如果读取失败或者波段数不匹配第一件事不是换函数而是打开.hdr检查三个字段samples、lines、bands和interleave。常见错误是 ENVI 文件里的波段数是 400读进来只有 100多半是interleaveBIL/BIP/BSQ被误判了需要在读取函数的参数里显式指定或者先按错位的波长次序读进来再重排。R 的索引从 1 开始ENVI 头文件里的波段名从 1 开始这部分通常对齐问题一般出在存储顺序而不是下标。如果你用的是公开数据集比如 ICVL 高光谱数据集省去手工拼矩阵的麻烦library(hsi) icvl - loadICVL() # 首次运行会从网络下载并缓存到本地 plot_image(icvl)hsi的下载类函数把图像立方体、波长和元数据一并封装好了读进来以后可以直接spectra()提取矩阵也可以配合自己的预处理流程走。离线环境下最稳妥的方式仍然是拿到原始 ENVI 文件按前面readSpeclib的路径走一遍公开数据集本身的格式不一定和 ENVI 完全一致别为了省事硬套函数。3. 高光谱预处理是一条流水线暗电流、白板校正和去噪的必调参数3.1 一套可以直接改的预处理脚本高光谱数据的预处理最忌讳零敲碎打。今天扣一个暗电流明天做一次白板校正参数口口相传等换了一台设备就全部失效。我一般把预处理固定成五个步骤写在同一个脚本里裁剪噪声波段、扣除暗电流、白板校正、平滑去噪、回写对象。下面这段代码可以直接套在自己的数据上# 1. 裁掉两端信噪比崩掉的波段 wl - wavelength(img) keep - wl 420 wl 2350 img - img[, keep] # 2. 暗电流扣除每个波段取暗像元均值广播到所有像元 dark_current - readSpeclib(path/to/dark.dat, type ENVI) dc - colMeans(spectra(dark_current), na.rm TRUE) dn - spectra(img) dn_corr - sweep(dn, 2, dc, -) # 3. 白板反射率校正得到相对反射率 white - readSpeclib(path/to/white.dat, type ENVI) white_ref - colMeans(spectra(white), na.rm TRUE) reflectance - dn_corr / (white_ref - dc) reflectance[reflectance 0] - 0 reflectance[reflectance 1.5] - 1.5 # 4. Savitzky-Golay 平滑去噪 library(prospectr) reflectance_smooth - savitzkyGolay(reflectance, m 0, p 1, w 11) # 5. 回写到原来的 speclib 对象 spectra(img) - reflectance_smooth这段脚本的逻辑很直白暗电流是传感器在没有光进入时仍然读到的响应必须在所有像元上先减掉白板校正则是用标准反射白板的像元均值做分母得到相对反射率。sweep(dn, 2, dc, -)的意思是按列广播也就是每个波段用同一个暗电流均值去减这个操作比apply快得多。reflectance里的负值基本来自传感器噪声直接置 0 是为了避免取对数或算夹角时出现 NaN上限截到 1.5 而不是 1是因为白板反射率管不到百分百轻微超过 1 是正常的过校正现象。第 4 步的savitzkyGolay参数最容易被人误解。m 0表示不平导数只做平滑p 1是拟合多项式阶数对高光谱反射率数据取 1 就够w 11是窗口宽度必须传奇数。窗口越大曲线越平滑但吸收特征也会被抹掉。一个比较稳的参考是传感器光谱分辨率在 10 nm 左右时窗口取 5分辨率在 2 nm 左右时窗口取 11 到 15。取 11 是从“肉眼看起来平滑同时保留叶绿素吸收谷”这个经验来的拿到自己的数据后可以先对 3 到 5 条光谱反复试看到 680 nm 附近的红边没有被拖平这组参数基本可以用。3.2 预处理完先别急着分析三步确认法预处理做没做对不能只看脚本有没有报错要用“三步确认法”验证第一画光谱曲线看形态第二检查特征吸收谷位置第三对比校正前后的统计量。画图用以下代码set.seed(2024) idx - sample(nSpectra(img), size min(15, nSpectra(img))) wl_plot - wavelength(img) spc_plot - spectra(img)[idx, ] matplot(wl_plot, t(spc_plot), type l, lty 1, col rgb(0, 0.4, 0.8, 0.6), xlab Wavelength / nm, ylab Reflectance, main Preprocessed spectra)正常反射率曲线应该是缓变的波段之间没有剧烈的锯齿。如果你看到 940 nm、1140 nm、1380 nm 附近有明显的下凹那不是错误那是水汽吸收带反而是判断波长位置有没有对齐的好标志。如果 1380 nm 处原本应该有个吸收谷画出来却是个凸起就要怀疑波长标定反了或者参考光谱没有重采样。白板校正解决的是“相对反射率”不是绝对辐亮度。在 R 的开源生态里我们不硬凑大气校正。如果你的目标只是做分类、丰度反演、同一区域不同日期的比对白板校正已经够用但如果你想把自己测的光谱和某个公开光谱库逐波段比对两个数据源的波长网格必须一致。hsdar提供了resampleSpectra可以把一个 0.1 nm 分辨率的点光谱重采样到影像的 10 nm 网格上library(hsdar) image_wl - wavelength(img) field_spec - readSpeclib(field_spectrum.txt, type Synthetic) field_resampled - resampleSpectra(field_spec, image_wl)这里的关键参数是目标波长向量必须来自影像本身而不是随意生成一个等差数列否则光谱角计算会偏差特别大。还有人喜欢在预处理阶段顺手做一阶导数我要泼一盆冷水对反射率直接做导数会放大高频噪声R 里没有所谓“自动帮你做好平滑再做导数”的函数必须先平滑再求导。这种顺序问题比任何参数都致命。4. 基本分析做到能出图PCA 降维、SAM 分类与像元分解的 R 实现4.1 降维先标准化PCA 和它的内存坑高光谱影像动辄几十万个像元、几百个波段直接丢进分类器之前必须先降维。很多人第一反应是用 MNF但 MNF 在 R 里没有稳定的一站式实现我常用的做法是“先做 PCA再做可视化验证”效果差不多关键是参数别错。PCA 的代码非常短但有两个隐藏细节标准化和内存控制。X - spectra(img) Xc - scale(X, center TRUE, scale TRUE) # 每个波段标准化 pca - prcomp(Xc, center FALSE, scale. FALSE, rank. 10) scores - pca$x[, 1:5]scale(X, center TRUE, scale TRUE)这一步不是可选项。高光谱波段之间的能量差异极大某些波段方差可能是其他波段的几百倍如果不标准化前几个主成分会被这些高方差波段牵着走最后画出来的得分图几乎就是一张噪声图。prcomp里的center FALSE, scale. FALSE不是说不需要中心化而是告诉prcomp输入已经标准化过了别再重复处理一遍。rank. 10是为了限制只计算前 10 个主成分能省下非常可观的内存。但prcomp的内部实现仍然需要把整个影像矩阵读进内存参与 SVD 分解。如果你碰上一块 400 波段、有几百万像元的影像prcomp很可能跑到一半就把 R 会话卡死。更稳妥的做法是用irlba包做截断 SVD它专门为主成分个数远小于矩阵维度的情况设计# install.packages(irlba) pca_irlba - irlba::prcomp_irlba(Xc, n 10, center FALSE, scale. FALSE) scores_irlba - pca_irlba$x[, 1:5]两个版本的得分结果在趋势上没有本质区别但irlba的内存占用低很多。降维后要看的不是主成分图本身而是方差解释率。一个合格的高光谱 PCA 结果前三个主成分至少解释 90% 以上的方差否则说明数据里噪声太大或者波段裁剪得太保守。用barplot(pca_irlba$sdev^2 / sum(pca_irlba$sdev^2))扫一眼如果第一个主成分就占到 80% 以上通常不是因为信息集中而是因为你忘了标准化。4.2 光谱角制图和线性混合像元分解光谱角制图SAM把每个像元的光谱看作高维空间里的一个向量用夹角衡量相似度。这个方法的优点是对光照强度不敏感适合直接用在白板校正后的相对反射率上。hsdar里提供了SAM函数# 假设已经有三个端元光谱植被、土壤、阴影行数代表端元数 endmember_matrix - rbind(soil_spec, veg_spec, shade_spec) colnames(endmember_matrix) - wavelength(img) ref - speclib(spectra endmember_matrix, wavelength wavelength(img)) sam_result - SAM(img, ref)SAM返回的矩阵行数是影像像元数列数是端元数每个值是弧度制的夹角越接近 0 说明越像。拿到结果后用apply(sam_result, 1, which.min)给每个像元分配类别标签就能出一张初级的分类图。这里最容易被忽略的细节是参考光谱的波长网格必须和影像一致否则参考光谱在 750 nm 处有个反射峰影像光谱在 755 nm 处才是峰值夹角会被人为拉大。用前面提到的resampleSpectra对齐一次再丢进SAM。混合像元分解解决的是“一个像元里混了多少种地物”的问题。最简单也最好用的线性混合模型假设每个像元是端元光谱的线性组合系数就是丰度。R 里可以用nnls包做非负最小二乘约束每个系数不小于零library(nnls) pixel - spectra(img)[500, ] # 任取一个像元 A - rbind(soil_spec, veg_spec) # 波段数 x 端元数 fit - nnls(A, pixel) abundance - fit$x / sum(fit$x) # 归一化到总和为 1nnls(A, pixel)的返回对象里x就是各端元的非负丰度值。之所以要除以总和是因为nnls本身只保证非负不保证和为 1。使用这个函数前要注意A必须是“波段数 × 端元数”的方向如果你的端元光谱矩阵是反的nnls会直接报维度错误这也是新手最容易翻车的地方。端元的选取不要拍脑袋一般做法是先在 PCA 得分图上找聚类的顶点把每个类里最“纯”的几个像元光谱平均一下作为端元光谱。这一步做完你和“遥感专业的高光谱分析”之间已经没有多少距离了。5. R 高光谱分析避坑这五个现象最像玄学其实都有确切原因5.1 数据读取阶段的三个坑现象一ENVI 文件读进来了但波段顺序是反的画出来的光谱曲线整体左右颠倒。原因通常是.hdr文件里的波长列表wavelength顺序和二进制体里的存储顺序不一致或者你在拷贝文件时只拷了.dat忘带了.hdrR 自动生成了一份错误的头。解决办法是先用readLines(scene.hdr)人工看一眼波长列表检查第一行波长的数值是否大于最后一行如果反了就执行img - img[, nBands(img):1]把波段顺序重排回来。不要相信任何自动修正必须亲眼确认 680 nm 和 750 nm 的位置再往下走。现象二白板校正之后图像出现整行整列的条纹还有些像元反射率变成几十倍。原因是白板不是在同一批次测量条件下采集的或者白板波段响应在紫外和短波红外端接近于零除法把噪声放大了。解决办法是先把白板光谱画出来找到白板响应低于暗电流响应 10 倍的波段区间这些波段在反射率计算中直接剔除而不是硬算完再筛。同时白板校正必须在暗电流扣除之后进行顺序反了分母出现负值结果就是一片乱码。现象三读取几百 MB 的.dat时 R 直接卡死连 RStudio 都无响应。原因是readSpeclib会把整个影像一次性读成矩阵几十万像元乘以几百波段的内存占用轻松超过几个 GB。解决办法不是在 R 里硬扛而是在读取前先用 GDAL 或 ENVI 把影像重采样比如空间上每隔一个像元采样一次把数据量先降到可以放进内存的规模。R 里也可以临时用terra::rast读入小范围子集再转成 speclib但主流程上我更推荐“先降采样再进分析”。高光谱数据做算法验证时期不需要原始全分辨率等流程定型后再上大机器跑全图。5.2 分析与内存阶段的两个坑现象四PCA 做完前两个主成分画出来不是地物分布而是一条一条的水汽噪声带。原因是波段标准化被跳过或者裁剪波段时只裁了两端没有剔除 1380 nm 和 1900 nm 附近的强水汽吸收带。解决办法分两步先按wl 420 wl 2350做一次粗裁再把 1350 到 1450 nm、1800 到 1950 nm 这两段明确剔除然后再scale(center TRUE, scale TRUE)。如果你发现标准化之后 PCA 还是很脏检查一下Xc里是否有个别波段的方差为 0那通常是因为原始数据里那几个波段全是同一个数值这种波段要整体删掉否则scale会得到一堆 NaN。现象五SAM 分类结果像椒盐图每个像元的标签在相邻像元之间跳来跳去。原因是单个像元的光谱噪声太大或者参考光谱来自不同波段网格。解决办法是先对影像做一次空间平滑再把平滑后的影像送进SAM。空间平滑可以用terra::focal对每个波段单独操作也可以用最笨的 3×3 均值滤波R 里面最容易理解的做法是把影像先转成terra的 SpatRaster然后library(terra) r - rast(scene_smooth.tif) r_smooth - focal(r, w matrix(1, 3, 3), fun mean, na.rm TRUE)注意focal是对每个图层分别做空间滤波所以波段之间不会互相污染。做完平滑后光谱角分类结果的斑块会明显连续。另一个隐藏问题是参考光谱没有重采样到影像波长网格这个前面强调过实际操作里比噪声更容易造成“整片错分”必须检查参考光谱和影像光谱的波长最大值是否一致。6. 让它落地再快一点分块计算、并行与结果验证分析流程跑通之后接下来面对的就是效率和可信度问题。一整块高光谱影像全量跑一次 PCA 和 SAM 可能很慢但我通常不会一上来就开并行而是先做一件事把像元按批次处理每个批次只保留特征波段上的数据。具体做法是把spectra(img)按 2 万像元一批切出来每批算出 PCA 得分和 SAM 夹角后把标签写进一个预分配的整数向量里批次之间互不影响自然就可以并行。Windows 下我习惯用future.applyLinux 下用mclapply但注意不要在 RStudio 里嵌套并行否则内存会被迅速吃光。验证这一步比优化速度更重要。我的习惯是永远留出 5% 的像元不参与任何分析和训练专门用来做“结果复核”。SAM 分类跑完后把这 5% 像元的类别归属和原始光谱曲线逐条对照看标签是不是落在合理的端元附近。PCA 降维也有个值得做的后验指标重建误差。把前五个主成分的得分乘回载荷矩阵再和标准化后的原始光谱比 RMSE如果误差在 0.1 个反射率单位以内说明降维没有伤到主要地物信息如果误差很大说明前五个主成分撑不住这个数据集需要增加主成分数量。记得重建时把scale的均值和标准差还原回来不然算出来的 RMSE 没有实际含义。这套流程走到这里已经覆盖了开源 R 做高光谱数据分析和处理的主干数据读取、预处理、降维、光谱角分类、混合像元分解以及关键位置的避坑。真到项目上线我习惯把每一批数据预处理后的统计量存成一个 RDS 文件保留所有参数和波长信息这样将来回看结果时能知道“这组分类图是在什么参数下生成的”。数据、参数、结果三件套齐全比任何花哨算法都让人安心。希望这些经验能帮你在 R 里少走一段弯路把时间留给真正要解释的地物问题。本文还有配套的精品资源点击获取