简介一份面向R语言用户的开源高光谱数据分析资源围绕hsdar包提供从数据导入、预处理到特征提取、分类建模及可视化的完整流程。内容涵盖ENVI、HDF、GeoTIFF等多种格式支持以及平滑、大气校正、主成分分析、支持向量机、随机森林等常用方法适合遥感、地学、农业等领域科研人员与工程师参考。压缩包共224个文件以76个R脚本、85个帮助文档rd为主另有25个Fortran源文件、3个C文件等底层算法实现并附带示例数据rdata/rds、PDF手册和图表包体约3.73MB。已有741人学习下载。通过该包可快速搭建高光谱分析环境直接调用封装好的函数完成分类、回归等典型任务源码中Fortran与C实现有助于理解算法底层逻辑配合帮助文档、示例数据及光谱绘图功能可支持二次开发、教学演示或科研复现。1. 用R处理高光谱数据比起ENVI好在哪用Excel给地物光谱画几条折线图容易但当你手里有二三十个样本、几百上千个波段还要区分健康与受胁迫植被时折线图立刻失效。R中高光谱数据的处理和基本分析是把散装光谱变成可复现结论的实用链路hsdar包提供Speclib数据结构prospectr负责平滑去噪caret能搭分类模型全部开源、脚本化、能出图。这个方向适合生态遥感、农业病害监测、土壤属性建模的从业者也适合不想被商业软件黑匣子卡脖子的人。读这篇文章你会拿到从格式解析到分类成图的完整路径以及几个坑了无数人的细节。2. 把ENVI和HDF5读进RSpeclib数据结构与三种存储格式的坑2.1 BSQ/BIL/BIP三种布局读错一个就全乱拿到一份高光谱数据第一件事不是急着跑分析而是搞清楚它是怎么存的。高光谱影像在磁盘上常见三种波段布局BSQ波段顺序、BIL波段行交错、BIP波段像元交错。BSQ最简单所有波段按顺序排列先存完整个波段1再存波段2BIL是几行一个波段来回切换BIP则是每个像元的所有波段连续存放。这三种布局直接影响读取速度和代码写法。用ENVI处理时不感知因为软件帮你做了但在R里如果你拿到的是纯二进制文件外加一个.hdr头文件读错了布局轻则数据错乱重则噪声占比暴涨。hsdar的read.ENVI会读.hdr里的interleave字段一般情况下不用手动管但自己写底层readBin时务必先看head。布局英文全称存储顺序适用场景BSQBand Sequential波段1全部像元 → 波段2全部像元单波段分析、文件压缩率高BILBand Interleaved by Line第1行所有波段 → 第2行所有波段兼顾空间与光谱读取BIPBand Interleaved by Pixel每个像元连续存完所有波段光谱维度操作频繁2.2 用hsdar构造和读取Speclib一段能直接跑的代码hsdar是R里处理高光谱最顺手的开源包核心数据结构叫Speclib它把波长序列、光谱矩阵和辅助信息绑在一起。先看怎么用模拟数据构造一个Speclib这个写法在你手上没有现成文件时用来练手非常合适library(hsdar) wl - seq(400, 2500, by 1) # 波长范围400~2500nm间隔1nm gen_spec - function(wl, chl_depth 0.25, water_depth 0.12) { base - 0.3 0.5 * exp(-((wl - 550) / 80)^2) chl - 1 - chl_depth * exp(-((wl - 680) / 20)^2) water - 1 - water_depth * exp(-((wl - 970) / 30)^2) base * chl * water } set.seed(42) spec_list - lapply(1:30, function(i) { group - ifelse(i 15, health, stress) depth - ifelse(group health, 0.30, 0.15) wd - ifelse(group health, 0.12, 0.10) gen_spec(wl, depth, wd) * (1 rnorm(length(wl), 0, 0.01)) }) spec_mat - do.call(rbind, spec_list) myspec - new(Speclib, wavelength wl, spectra spec_mat, SI data.frame(group rep(c(health, stress), each 15)))这里的关键是gen_spec里三个高斯函数的含义550nm附近的峰模拟反射高峰680nm是叶绿素吸收中心970nm是水分吸收带。chl_depth控制吸收深度健康叶片吸收深胁迫叶片吸收浅。new(Speclib, ...)的三个参数分别是波长、光谱矩阵和辅助数据SISI可以理解为每个样本的标签簿后面做分组、分类全靠它。如果你有现成的野外实测光谱通常是ENVI格式一个.dat文件配一个.hdr文件。读取就一行spec - read.ENVI(grep_files field_spectra.dat, header_file field_spectra.hdr)read.ENVI的grep_files参数填二进制文件header_file填头文件。读完之后立刻用str(spec)看结构用wavelength(spec)确认波长范围。我见过不少人读进来不检查直接跑分析结果发现单位是微米而不是纳米所有指数全部算偏。2.3 高光谱影像的兜底读法terra转Speclib实测点光谱用read.ENVI就够但如果你拿到的是高光谱影像比如机载或星载数据就不能直接用hsdar读了。常见做法是先用terra把影像读进来再转成Speclib因为hsdar的光谱计算函数只认Speclib。library(terra) img - rast(hyperspectral_mosaic.tif) vals - as.matrix(img, wide TRUE) wl - as.numeric(gsub(Band_, , names(img))) img_spec - new(Speclib, wavelength wl, spectra vals)as.matrix(img, wide TRUE)会把栅格对象转成行是像元、列是波段的矩阵这一步恰好是Speclib需要的结构。波长序列如果写在波段名称里就用names(img)拆出来写不进名称就从.hdr里抄千万别自己猜。这条兜底路径能处理绝大多数TIFF格式的高光谱影像但要注意内存波段多、像元多时as.matrix一次全进内存在小机器上大概率崩后面避坑章会专门讲分块。3. R中高光谱预处理实操平滑、掩膜和连续统去除怎么配参数3.1 先画光谱曲线噪声波段一眼就能看出拿到Speclib后的第一个动作永远是plot不画图就直接预处理是自找麻烦。plot输出所有样本的光谱曲线一眼能看出来几件事首尾波段是不是噪声爆表有没有负值曲线整体形态是否合理。地物光谱仪在350~400nm和2400~2500nm两个区间噪声几乎不可避免因为探测器在这两端的信噪比很低。plot(myspec, col grey80)画完图你会看到400nm附近有明显的锯齿2500nm附近则容易乱跳。这里我一般会记录下噪声波段的索引范围后面用掩膜剔掉。掩膜是hsdar里一个重要的概念它标记哪些波段参与计算哪些被排除。设置掩膜用mask函数波段索引从1开始闭区间包含端点mask(myspec) - c(1, 50, 2051, 2101)上面的代码把第1~50和第2051~2101个波段排除掉对应的是400~450nm和2450~2500nm区域。掩膜作用于对象本身后续vegindex、cr这类函数会自动跳过掩膜波段。这个特性很实用因为它不像手动删列那样破坏原始数据想恢复就直接把掩膜设为NULL。3.2 S-G平滑参数怎么设窗口宽度和多项式阶数平滑是高光谱预处理里最容易“好心办坏事”的一步。我用的是prospectr包里的savitzkyGolay它基于Savitzky-Golay卷积平滑核心参数就三个多项式阶数p、窗口宽度w、导数阶数m。library(prospectr) clean_mat - savitzkyGolay(X spectra(myspec), p 3, w 31, m 0) clean_spec - new(Speclib, wavelength wavelength(myspec), spectra clean_mat, SI SI(myspec))p3是多项式阶数对大多数反射光谱够用w31是窗口宽度必须是奇数数值越大平滑力度越强但也会把窄吸收峰抹平m0表示不做导数保留原始反射率形态。如果后续只想做分类很多人会把m设为1或2求一阶导、二阶导来消除基线漂移但代价是噪声被同步放大。窗口宽度怎么选这是经验活。400~2500nm的连续光谱w11到w31都是安全区间。小于11平滑效果约等于没有大于51680nm的叶绿素吸收峰会被削平你后面算的红边参数就失真了。我自己的习惯是先w21跑一版对比w31的结果如果分类精度没有明显变化用更小的窗口保留更多细节。平滑是对所有样本批量做的注意X参数要传矩阵不能直接传Speclib对象。3.3 坏波段剔除与连续统去除让吸收峰站在同一条基线上高光谱分析里有个绕不开的操作叫连续统去除英文是continuum removalhsdar里对应cr函数。它的作用是把光谱的包络线拉平让不同样本的吸收特征在同一个基线上比较。没有这一步两个样本的反射率基线一个高一个低直接比680nm吸收深度对比的是基线而不是真正的吸收强度。cr_spec - cr(clean_spec, method is) plot(cr_spec)methodis是internal standard的缩写用凸包逐段包络原始光谱然后逐波段做除法得到连续统去除后的光谱。结果是0到1之间的曲线吸收峰变成谷谷越深代表吸收越强。另一种methodrms是均方根法适合波形差异较大的场景实际工作中用的少。连续统去除后的光谱还能直接算植被指数hsdar的vegindex一次性能算几十种指数vi - vegindex(cr_spec, index c(NDVI, PRI, NDWI, ARI))vegindex的index参数可以传多个指数名返回一个矩阵行是样本列是指数。这里有个容易忽略的点这些指数的计算公式内部默认波长单位是纳米如果你的数据是微米它不会报错但搜到的波段索引不对算出来的值就错得离谱。所以读数据之后第一件事永远是确认wavelength()返回值这个细节值得刻在脑门上。预处理到这一步光谱已经可以用于分析和建模了。剩下的工作是把几百上千个波段压缩成真正有用的特征这就要进入降维和特征提取环节。4. 降维与特征提取PCA、MNF哪个更适合高光谱4.1 对Speclib做PCA载荷谱里藏着哪些波段高光谱最麻烦的地方是波段太多一个样本动辄几百上千个变量直接丢进分类器不仅慢还容易过拟合。降维第一板斧是主成分分析。hsdar为Speclib重载了princomp方法可以直接对Speclib对象做主成分分析不用手动拆矩阵pca_res - princomp(cr_spec) summary(pca_res) loadings(pca_res) plot(loadings(pca_res)[, 1:3], type l)princomp的内部逻辑是对光谱矩阵做主成分分解前几个主成分集中了绝大部分方差。高光谱数据里前3~5个主成分通常能解释95%以上的方差剩余分量大多是噪声。summary会输出每个主成分的方差贡献率我一般看累计贡献率到95%需要几个分量。loadings的结果是每个波段在主成分上的权重把它画成曲线就可以看到哪些波段贡献最大。比如第一主成分载荷曲线在680nm附近出现高峰说明叶绿素吸收波段是样本间差异的主要来源。这个信息对后续选波段非常有价值它告诉你做特征筛选时优先保留哪些位置。4.2 为什么机载影像推荐MNF它比PCA强在哪PCA的缺陷在于它对噪声太诚实噪声方差大的波段会被当成重要信息前几个分量可能混入大量噪声。对野外点光谱影响还可控对机载高光谱影像就很致命因为传感器噪声、大气残留、条带噪声都叠在里面。这个场景我一般改用MNFMinimum Noise Fraction最小噪声分数变换hsdar里提供了mnf函数mnf_res - mnf(cr_spec) summary(mnf_res)MNF的思路是先估计噪声协方差矩阵把数据白化之后再做一次PCA。它的实际效果是让信息量最大的分量排在前面噪声被压到后面几个分量。在影像上MNF的前三个分量往往是清晰的地物空间分布而PCA的前三个分量经常被条带噪声污染这是两者最直观的差别。用MNF时有个细节对于Speclib它要求数据已经做过掩膜处理因为掩膜外的噪声波段会影响协方差估计。另一件要注意的事情是做MNF之前最好把光谱重采样到统一间隔否则波长分辨率不一致会影响变换结果。MNF的返回值可以直接拿去分类不需要还原到原始光谱。4.3 波段落地vegindex和吸收特征深度怎么选变量降维并不只有PCA和MNF一条路。很多时候你要的不是“压缩后的抽象分量”而是“具体哪个波长、哪个指标能表征这种胁迫”。这也是我认为高光谱分析最有价值的地方找到物理意义明确的特征波段。前面算过的vegindex结果本身就是一组好特征。NDVI对绿度敏感PRI对光合效率敏感NDWI对水分敏感ARI对花青素敏感。把这几个指数拼成特征矩阵加上少量关键波段反射率往往比直接丢几百个波段进模型效果更好也更容易向非遥感背景的人解释。连续统去除之后还有一个利器是吸收特征参数提取hsdar的feature函数可以提取指定波段范围内的吸收峰深度、宽度和面积feat - feature(cr_spec, wavelength c(550, 800))这个调用的含义是在550~800nm区间提取吸收特征对应叶绿素在红光段的吸收。输出的深度、宽度和面积三个参数里深度用得最多它直接反映吸收强弱。参数名在不同版本里可能有差别跑之前先用args(feature)确认一下当前包版本的签名。把feat的结果和vegindex合并成一个数据框就得到了一份“光谱特征表”每行是一个样本每列是一个有物理含义的特征。后面的分类、回归就基于这张表跑比直接拿原始光谱跑要稳得多模型训练时间和过拟合风险也降了一个量级。5. R中高光谱分析血泪避坑五个高频翻车现场5.1 波长单位是微米还是纳米指数算出来全是NaN现象vegindex跑完之后返回矩阵里有大量NaN或者NDVI值异常有的超过1有的为负。 原因数据来源不同波长单位不统一。ASD地物光谱仪导出的ENVI文件波长常以微米为单位数值范围0.4~2.5而hsdar内部所有指数公式都假定波长单位是纳米。单位不匹配时公式去索引特定波段位置找到的全是空气。 解决读入后立刻判断并转换。wl - wavelength(myspec) if (max(wl) 10) { wavelength(myspec) - wl * 1000 }这个判断逻辑是如果最大波长不超过10则几乎可以断定是微米直接乘以1000转成纳米。之后再用wavelength()确认一下范围在400~2500。5.2 read.ENVI报unknown data type先查hdr的data type字段现象read.ENVI读文件直接报错提示“unknown data type”或者读进来的数据全是乱的数值在百万级别。 原因ENVI的.hdr里有data type字段1是Byte、2是Int、4是Float、5是Double等。如果原文件是用非标准工具生成的hdr里data type写错或者写入了R不认识的值hsdar解析就会失败。 解决用文本编辑器打开.hdr看interleave和data type两个字段。如果data type不对最省事的办法是用ENVI或GDAL把数据重存成float32data type4同时确保interleave和实际存储一致。千万不要自己去改.hdr里data type的数值文件里的二进制布局和声明不符时硬改头文件读出来的还是乱码。5.3 整幅高光谱影像塞进内存R直接崩给你看现象用terra读一幅3000行乘3000列、200个波段的高光谱影像as.matrix之后R直接内存不足退出或者等了几分钟然后崩溃。 原因算一笔账3000乘3000乘200乘4字节单精度浮点就是约6.7GB转成矩阵后内存翻倍还没算光谱矩阵在R里的对象开销。把整幅影像变成Speclib在小机器上就是自杀。 解决分块处理按行读取。img - rast(hyperspectral_mosaic.tif) nr - nrow(img) block_size - 512 for (i in 1:ceiling(nr / block_size)) { rows - ((i - 1) * block_size 1):min(i * block_size, nr) block - readValues(img, row rows[1], nrows length(rows), col 1, ncol ncol(img)) block_spec - new(Speclib, wavelength wl, spectra block, SI data.frame(block_id i)) # 这里做预处理或特征提取然后把结果拼到汇总表里 }readValues读出来的是矩阵每行一个像元每列一个波段。分块的关键在于每块处理完只保留你要的结果比如vegindex、SAM角度、分类标签不保留原始光谱这样内存永远只维持在一个block的量级。如果你真的需要全图光谱那建议用一个专门处理大影像的工具先做MNF降维把200波段压缩到10个分量再进R。5.4 重采样把特征波段弄丢了PRI返空值怎么排查现象执行resample把光谱重采样到自定义波长序列之后vegindex的PRI返回空值其他指数正常。 原因PRI需要530nm和570nm两个具体波段。如果重采样时新波长序列没有覆盖这两个点或者间隔太粗跳过了公式找不到对应位置就返回空。多数情况下是重采样时波长范围和间隔没设计好直接把两端截断了。 解决重采样之后先验证波长覆盖范围。rs - resample(cr_spec, seq(400, 2400, by 5), method linear) range(wavelength(rs))检查range是不是还覆盖530和570。如果范围收窄了把新波长序列的起点和终点往两边各扩一点。另一个经验是重采样间隔不要超过5nm超过5nm窄吸收特征会直接消失算出来的指数全都不靠谱。重采样是有损操作能不做就不做必须做时先想清楚你要保留哪些特征波段。5.5 交叉验证精度虚高根源是空间自相关现象用caret跑随机森林五折交叉验证准确率95%样本一到独立地块上验证直接掉到60%。你以为模型稳定其实是交叉验证骗了你。 原因高光谱影像里相邻像元的光谱高度相似。训练集和验证集如果按像元简单随机划分同一个地块或同一棵树冠的光谱会同时出现在训练和验证里模型等于“背过答案”。这就是空间自相关导致的精度虚高。 解决按地块或空间单元分组做交叉验证保证同一个地块不跨训练和验证。library(caret) plot_id - SI(myspec)$plot_id folds - createFolds(unique(plot_id), k 5) fold_idx - lapply(folds, function(f) which(plot_id %in% f)) ctrl - trainControl(method cv, index fold_idx) model - train(x feature_matrix, y group, method rf, trControl ctrl)createFolds把地块ID分成五份fold_idx再把地块映射回像元索引。这样同等地块的像元永远在同一折里模型见不到“同类邻居”。跑出来的精度才是可信的数字很多文章里95%的精度就是这么注水的。6. 进阶技巧光谱角制图加交叉验证一条命令出分类结果6.1 SAM最小实现一条命令算出全图光谱角光谱角制图Spectral Angle Mapper, SAM是高光谱分类里用得最成熟的算法它把每条光谱看作高维空间的一个向量用两向量之间的夹角衡量相似度。夹角越小越相似原理简单抗增益变化、不敏感到光照差异因为这些因素主要影响光谱幅度不影响光谱形状。hsdar提供了sam函数直接把Speclib和参考光谱丢进去就能算ref - avgsp(cr_spec[SI(cr_spec)$group health, ]) angles - sam(cr_spec, ref) hist(angles, breaks 50, main SAM角度分布)avgsp是hsdar的平均函数用于取一组光谱的平均值作为参考光谱这里取健康样本的平均光谱。sam返回每个样本与参考光谱夹角的余弦值转换后就是角度。看分布直方图如果两个类别的角度分布有明显分界直接用阈值就能分类。经验值上0.05~0.1弧度内的光谱可以视为同类0.2以上基本是不同地物。这个阈值对模拟数据很理想实测数据因为有噪声会模糊一些可以先画分布再定阈值。6.2 给分类结果上“后悔药”caret五折交叉验证SAM阈值分类快但阈值怎么定总带点玄学。要对分类结果有把握正规做法还是得跑一遍交叉验证。用caret配randomForest五折交叉验证就能给出客观的准确率估计。有了这个数字你再去定SAM阈值心里就有底了。df - as.data.frame(spectra(cr_spec)) df$group - SI(cr_spec)$group ctrl - trainControl(method cv, number 5) model - train(group ~ ., data df, method rf, trControl ctrl, tuneLength 5) print(model$results)train的group ~ .表示用所有光谱波段预测分组methodrf指定随机森林tuneLength5会自动试5组mtry参数。model$results里能看到不同mtry下的准确率选最高的一行作为最终参数。把这套结果和SAM分类图放在一起看互相验证比单靠任何一种方法都踏实。6.3 出图规范ggplot导出高分辨率jpg和pdf分析做完了最终要落到一张能放进报告里的图。R语言里导出图片用ggplot加ggsave能同时出jpg和pdf两个版本jpg用于快速预览和微信传输pdf用于论文投稿和打印。library(ggplot2) sam_df - data.frame(x rep(1:50, 50), y rep(1:50, each 50), angle as.vector(angles)) p - ggplot(sam_df, aes(x x, y y, fill angle)) geom_raster() scale_fill_viridis_c() theme_minimal() ggsave(sam_surface.pdf, p, width 8, height 6) ggsave(sam_surface.jpg, p, width 8, height 6, dpi 300)dpi300是印刷级分辨率务必写清楚否则默认的72dpi导出的图放大全是马赛克。pdf不用指定dpi它是矢量格式缩放不丢细节。我现在的习惯是每次出图都同时存pdf和jpg避免投稿时临时找不到矢量图再跑一遍数据。整个流程走下来你会发现R这套开源方案真正值钱的不是某一个算法而是整个链路可复现原始光谱进分类图出每一步的参数和代码都在换一批数据也能立刻重跑。这种踏实感是商业软件给不了的。希望帮到你。本文还有配套的精品资源点击获取 R语言高光谱数据分析全流程：从数据读取到分类可视化
简介一份面向R语言用户的开源高光谱数据分析资源围绕hsdar包提供从数据导入、预处理到特征提取、分类建模及可视化的完整流程。内容涵盖ENVI、HDF、GeoTIFF等多种格式支持以及平滑、大气校正、主成分分析、支持向量机、随机森林等常用方法适合遥感、地学、农业等领域科研人员与工程师参考。压缩包共224个文件以76个R脚本、85个帮助文档rd为主另有25个Fortran源文件、3个C文件等底层算法实现并附带示例数据rdata/rds、PDF手册和图表包体约3.73MB。已有741人学习下载。通过该包可快速搭建高光谱分析环境直接调用封装好的函数完成分类、回归等典型任务源码中Fortran与C实现有助于理解算法底层逻辑配合帮助文档、示例数据及光谱绘图功能可支持二次开发、教学演示或科研复现。1. 用R处理高光谱数据比起ENVI好在哪用Excel给地物光谱画几条折线图容易但当你手里有二三十个样本、几百上千个波段还要区分健康与受胁迫植被时折线图立刻失效。R中高光谱数据的处理和基本分析是把散装光谱变成可复现结论的实用链路hsdar包提供Speclib数据结构prospectr负责平滑去噪caret能搭分类模型全部开源、脚本化、能出图。这个方向适合生态遥感、农业病害监测、土壤属性建模的从业者也适合不想被商业软件黑匣子卡脖子的人。读这篇文章你会拿到从格式解析到分类成图的完整路径以及几个坑了无数人的细节。2. 把ENVI和HDF5读进RSpeclib数据结构与三种存储格式的坑2.1 BSQ/BIL/BIP三种布局读错一个就全乱拿到一份高光谱数据第一件事不是急着跑分析而是搞清楚它是怎么存的。高光谱影像在磁盘上常见三种波段布局BSQ波段顺序、BIL波段行交错、BIP波段像元交错。BSQ最简单所有波段按顺序排列先存完整个波段1再存波段2BIL是几行一个波段来回切换BIP则是每个像元的所有波段连续存放。这三种布局直接影响读取速度和代码写法。用ENVI处理时不感知因为软件帮你做了但在R里如果你拿到的是纯二进制文件外加一个.hdr头文件读错了布局轻则数据错乱重则噪声占比暴涨。hsdar的read.ENVI会读.hdr里的interleave字段一般情况下不用手动管但自己写底层readBin时务必先看head。布局英文全称存储顺序适用场景BSQBand Sequential波段1全部像元 → 波段2全部像元单波段分析、文件压缩率高BILBand Interleaved by Line第1行所有波段 → 第2行所有波段兼顾空间与光谱读取BIPBand Interleaved by Pixel每个像元连续存完所有波段光谱维度操作频繁2.2 用hsdar构造和读取Speclib一段能直接跑的代码hsdar是R里处理高光谱最顺手的开源包核心数据结构叫Speclib它把波长序列、光谱矩阵和辅助信息绑在一起。先看怎么用模拟数据构造一个Speclib这个写法在你手上没有现成文件时用来练手非常合适library(hsdar) wl - seq(400, 2500, by 1) # 波长范围400~2500nm间隔1nm gen_spec - function(wl, chl_depth 0.25, water_depth 0.12) { base - 0.3 0.5 * exp(-((wl - 550) / 80)^2) chl - 1 - chl_depth * exp(-((wl - 680) / 20)^2) water - 1 - water_depth * exp(-((wl - 970) / 30)^2) base * chl * water } set.seed(42) spec_list - lapply(1:30, function(i) { group - ifelse(i 15, health, stress) depth - ifelse(group health, 0.30, 0.15) wd - ifelse(group health, 0.12, 0.10) gen_spec(wl, depth, wd) * (1 rnorm(length(wl), 0, 0.01)) }) spec_mat - do.call(rbind, spec_list) myspec - new(Speclib, wavelength wl, spectra spec_mat, SI data.frame(group rep(c(health, stress), each 15)))这里的关键是gen_spec里三个高斯函数的含义550nm附近的峰模拟反射高峰680nm是叶绿素吸收中心970nm是水分吸收带。chl_depth控制吸收深度健康叶片吸收深胁迫叶片吸收浅。new(Speclib, ...)的三个参数分别是波长、光谱矩阵和辅助数据SISI可以理解为每个样本的标签簿后面做分组、分类全靠它。如果你有现成的野外实测光谱通常是ENVI格式一个.dat文件配一个.hdr文件。读取就一行spec - read.ENVI(grep_files field_spectra.dat, header_file field_spectra.hdr)read.ENVI的grep_files参数填二进制文件header_file填头文件。读完之后立刻用str(spec)看结构用wavelength(spec)确认波长范围。我见过不少人读进来不检查直接跑分析结果发现单位是微米而不是纳米所有指数全部算偏。2.3 高光谱影像的兜底读法terra转Speclib实测点光谱用read.ENVI就够但如果你拿到的是高光谱影像比如机载或星载数据就不能直接用hsdar读了。常见做法是先用terra把影像读进来再转成Speclib因为hsdar的光谱计算函数只认Speclib。library(terra) img - rast(hyperspectral_mosaic.tif) vals - as.matrix(img, wide TRUE) wl - as.numeric(gsub(Band_, , names(img))) img_spec - new(Speclib, wavelength wl, spectra vals)as.matrix(img, wide TRUE)会把栅格对象转成行是像元、列是波段的矩阵这一步恰好是Speclib需要的结构。波长序列如果写在波段名称里就用names(img)拆出来写不进名称就从.hdr里抄千万别自己猜。这条兜底路径能处理绝大多数TIFF格式的高光谱影像但要注意内存波段多、像元多时as.matrix一次全进内存在小机器上大概率崩后面避坑章会专门讲分块。3. R中高光谱预处理实操平滑、掩膜和连续统去除怎么配参数3.1 先画光谱曲线噪声波段一眼就能看出拿到Speclib后的第一个动作永远是plot不画图就直接预处理是自找麻烦。plot输出所有样本的光谱曲线一眼能看出来几件事首尾波段是不是噪声爆表有没有负值曲线整体形态是否合理。地物光谱仪在350~400nm和2400~2500nm两个区间噪声几乎不可避免因为探测器在这两端的信噪比很低。plot(myspec, col grey80)画完图你会看到400nm附近有明显的锯齿2500nm附近则容易乱跳。这里我一般会记录下噪声波段的索引范围后面用掩膜剔掉。掩膜是hsdar里一个重要的概念它标记哪些波段参与计算哪些被排除。设置掩膜用mask函数波段索引从1开始闭区间包含端点mask(myspec) - c(1, 50, 2051, 2101)上面的代码把第1~50和第2051~2101个波段排除掉对应的是400~450nm和2450~2500nm区域。掩膜作用于对象本身后续vegindex、cr这类函数会自动跳过掩膜波段。这个特性很实用因为它不像手动删列那样破坏原始数据想恢复就直接把掩膜设为NULL。3.2 S-G平滑参数怎么设窗口宽度和多项式阶数平滑是高光谱预处理里最容易“好心办坏事”的一步。我用的是prospectr包里的savitzkyGolay它基于Savitzky-Golay卷积平滑核心参数就三个多项式阶数p、窗口宽度w、导数阶数m。library(prospectr) clean_mat - savitzkyGolay(X spectra(myspec), p 3, w 31, m 0) clean_spec - new(Speclib, wavelength wavelength(myspec), spectra clean_mat, SI SI(myspec))p3是多项式阶数对大多数反射光谱够用w31是窗口宽度必须是奇数数值越大平滑力度越强但也会把窄吸收峰抹平m0表示不做导数保留原始反射率形态。如果后续只想做分类很多人会把m设为1或2求一阶导、二阶导来消除基线漂移但代价是噪声被同步放大。窗口宽度怎么选这是经验活。400~2500nm的连续光谱w11到w31都是安全区间。小于11平滑效果约等于没有大于51680nm的叶绿素吸收峰会被削平你后面算的红边参数就失真了。我自己的习惯是先w21跑一版对比w31的结果如果分类精度没有明显变化用更小的窗口保留更多细节。平滑是对所有样本批量做的注意X参数要传矩阵不能直接传Speclib对象。3.3 坏波段剔除与连续统去除让吸收峰站在同一条基线上高光谱分析里有个绕不开的操作叫连续统去除英文是continuum removalhsdar里对应cr函数。它的作用是把光谱的包络线拉平让不同样本的吸收特征在同一个基线上比较。没有这一步两个样本的反射率基线一个高一个低直接比680nm吸收深度对比的是基线而不是真正的吸收强度。cr_spec - cr(clean_spec, method is) plot(cr_spec)methodis是internal standard的缩写用凸包逐段包络原始光谱然后逐波段做除法得到连续统去除后的光谱。结果是0到1之间的曲线吸收峰变成谷谷越深代表吸收越强。另一种methodrms是均方根法适合波形差异较大的场景实际工作中用的少。连续统去除后的光谱还能直接算植被指数hsdar的vegindex一次性能算几十种指数vi - vegindex(cr_spec, index c(NDVI, PRI, NDWI, ARI))vegindex的index参数可以传多个指数名返回一个矩阵行是样本列是指数。这里有个容易忽略的点这些指数的计算公式内部默认波长单位是纳米如果你的数据是微米它不会报错但搜到的波段索引不对算出来的值就错得离谱。所以读数据之后第一件事永远是确认wavelength()返回值这个细节值得刻在脑门上。预处理到这一步光谱已经可以用于分析和建模了。剩下的工作是把几百上千个波段压缩成真正有用的特征这就要进入降维和特征提取环节。4. 降维与特征提取PCA、MNF哪个更适合高光谱4.1 对Speclib做PCA载荷谱里藏着哪些波段高光谱最麻烦的地方是波段太多一个样本动辄几百上千个变量直接丢进分类器不仅慢还容易过拟合。降维第一板斧是主成分分析。hsdar为Speclib重载了princomp方法可以直接对Speclib对象做主成分分析不用手动拆矩阵pca_res - princomp(cr_spec) summary(pca_res) loadings(pca_res) plot(loadings(pca_res)[, 1:3], type l)princomp的内部逻辑是对光谱矩阵做主成分分解前几个主成分集中了绝大部分方差。高光谱数据里前3~5个主成分通常能解释95%以上的方差剩余分量大多是噪声。summary会输出每个主成分的方差贡献率我一般看累计贡献率到95%需要几个分量。loadings的结果是每个波段在主成分上的权重把它画成曲线就可以看到哪些波段贡献最大。比如第一主成分载荷曲线在680nm附近出现高峰说明叶绿素吸收波段是样本间差异的主要来源。这个信息对后续选波段非常有价值它告诉你做特征筛选时优先保留哪些位置。4.2 为什么机载影像推荐MNF它比PCA强在哪PCA的缺陷在于它对噪声太诚实噪声方差大的波段会被当成重要信息前几个分量可能混入大量噪声。对野外点光谱影响还可控对机载高光谱影像就很致命因为传感器噪声、大气残留、条带噪声都叠在里面。这个场景我一般改用MNFMinimum Noise Fraction最小噪声分数变换hsdar里提供了mnf函数mnf_res - mnf(cr_spec) summary(mnf_res)MNF的思路是先估计噪声协方差矩阵把数据白化之后再做一次PCA。它的实际效果是让信息量最大的分量排在前面噪声被压到后面几个分量。在影像上MNF的前三个分量往往是清晰的地物空间分布而PCA的前三个分量经常被条带噪声污染这是两者最直观的差别。用MNF时有个细节对于Speclib它要求数据已经做过掩膜处理因为掩膜外的噪声波段会影响协方差估计。另一件要注意的事情是做MNF之前最好把光谱重采样到统一间隔否则波长分辨率不一致会影响变换结果。MNF的返回值可以直接拿去分类不需要还原到原始光谱。4.3 波段落地vegindex和吸收特征深度怎么选变量降维并不只有PCA和MNF一条路。很多时候你要的不是“压缩后的抽象分量”而是“具体哪个波长、哪个指标能表征这种胁迫”。这也是我认为高光谱分析最有价值的地方找到物理意义明确的特征波段。前面算过的vegindex结果本身就是一组好特征。NDVI对绿度敏感PRI对光合效率敏感NDWI对水分敏感ARI对花青素敏感。把这几个指数拼成特征矩阵加上少量关键波段反射率往往比直接丢几百个波段进模型效果更好也更容易向非遥感背景的人解释。连续统去除之后还有一个利器是吸收特征参数提取hsdar的feature函数可以提取指定波段范围内的吸收峰深度、宽度和面积feat - feature(cr_spec, wavelength c(550, 800))这个调用的含义是在550~800nm区间提取吸收特征对应叶绿素在红光段的吸收。输出的深度、宽度和面积三个参数里深度用得最多它直接反映吸收强弱。参数名在不同版本里可能有差别跑之前先用args(feature)确认一下当前包版本的签名。把feat的结果和vegindex合并成一个数据框就得到了一份“光谱特征表”每行是一个样本每列是一个有物理含义的特征。后面的分类、回归就基于这张表跑比直接拿原始光谱跑要稳得多模型训练时间和过拟合风险也降了一个量级。5. R中高光谱分析血泪避坑五个高频翻车现场5.1 波长单位是微米还是纳米指数算出来全是NaN现象vegindex跑完之后返回矩阵里有大量NaN或者NDVI值异常有的超过1有的为负。 原因数据来源不同波长单位不统一。ASD地物光谱仪导出的ENVI文件波长常以微米为单位数值范围0.4~2.5而hsdar内部所有指数公式都假定波长单位是纳米。单位不匹配时公式去索引特定波段位置找到的全是空气。 解决读入后立刻判断并转换。wl - wavelength(myspec) if (max(wl) 10) { wavelength(myspec) - wl * 1000 }这个判断逻辑是如果最大波长不超过10则几乎可以断定是微米直接乘以1000转成纳米。之后再用wavelength()确认一下范围在400~2500。5.2 read.ENVI报unknown data type先查hdr的data type字段现象read.ENVI读文件直接报错提示“unknown data type”或者读进来的数据全是乱的数值在百万级别。 原因ENVI的.hdr里有data type字段1是Byte、2是Int、4是Float、5是Double等。如果原文件是用非标准工具生成的hdr里data type写错或者写入了R不认识的值hsdar解析就会失败。 解决用文本编辑器打开.hdr看interleave和data type两个字段。如果data type不对最省事的办法是用ENVI或GDAL把数据重存成float32data type4同时确保interleave和实际存储一致。千万不要自己去改.hdr里data type的数值文件里的二进制布局和声明不符时硬改头文件读出来的还是乱码。5.3 整幅高光谱影像塞进内存R直接崩给你看现象用terra读一幅3000行乘3000列、200个波段的高光谱影像as.matrix之后R直接内存不足退出或者等了几分钟然后崩溃。 原因算一笔账3000乘3000乘200乘4字节单精度浮点就是约6.7GB转成矩阵后内存翻倍还没算光谱矩阵在R里的对象开销。把整幅影像变成Speclib在小机器上就是自杀。 解决分块处理按行读取。img - rast(hyperspectral_mosaic.tif) nr - nrow(img) block_size - 512 for (i in 1:ceiling(nr / block_size)) { rows - ((i - 1) * block_size 1):min(i * block_size, nr) block - readValues(img, row rows[1], nrows length(rows), col 1, ncol ncol(img)) block_spec - new(Speclib, wavelength wl, spectra block, SI data.frame(block_id i)) # 这里做预处理或特征提取然后把结果拼到汇总表里 }readValues读出来的是矩阵每行一个像元每列一个波段。分块的关键在于每块处理完只保留你要的结果比如vegindex、SAM角度、分类标签不保留原始光谱这样内存永远只维持在一个block的量级。如果你真的需要全图光谱那建议用一个专门处理大影像的工具先做MNF降维把200波段压缩到10个分量再进R。5.4 重采样把特征波段弄丢了PRI返空值怎么排查现象执行resample把光谱重采样到自定义波长序列之后vegindex的PRI返回空值其他指数正常。 原因PRI需要530nm和570nm两个具体波段。如果重采样时新波长序列没有覆盖这两个点或者间隔太粗跳过了公式找不到对应位置就返回空。多数情况下是重采样时波长范围和间隔没设计好直接把两端截断了。 解决重采样之后先验证波长覆盖范围。rs - resample(cr_spec, seq(400, 2400, by 5), method linear) range(wavelength(rs))检查range是不是还覆盖530和570。如果范围收窄了把新波长序列的起点和终点往两边各扩一点。另一个经验是重采样间隔不要超过5nm超过5nm窄吸收特征会直接消失算出来的指数全都不靠谱。重采样是有损操作能不做就不做必须做时先想清楚你要保留哪些特征波段。5.5 交叉验证精度虚高根源是空间自相关现象用caret跑随机森林五折交叉验证准确率95%样本一到独立地块上验证直接掉到60%。你以为模型稳定其实是交叉验证骗了你。 原因高光谱影像里相邻像元的光谱高度相似。训练集和验证集如果按像元简单随机划分同一个地块或同一棵树冠的光谱会同时出现在训练和验证里模型等于“背过答案”。这就是空间自相关导致的精度虚高。 解决按地块或空间单元分组做交叉验证保证同一个地块不跨训练和验证。library(caret) plot_id - SI(myspec)$plot_id folds - createFolds(unique(plot_id), k 5) fold_idx - lapply(folds, function(f) which(plot_id %in% f)) ctrl - trainControl(method cv, index fold_idx) model - train(x feature_matrix, y group, method rf, trControl ctrl)createFolds把地块ID分成五份fold_idx再把地块映射回像元索引。这样同等地块的像元永远在同一折里模型见不到“同类邻居”。跑出来的精度才是可信的数字很多文章里95%的精度就是这么注水的。6. 进阶技巧光谱角制图加交叉验证一条命令出分类结果6.1 SAM最小实现一条命令算出全图光谱角光谱角制图Spectral Angle Mapper, SAM是高光谱分类里用得最成熟的算法它把每条光谱看作高维空间的一个向量用两向量之间的夹角衡量相似度。夹角越小越相似原理简单抗增益变化、不敏感到光照差异因为这些因素主要影响光谱幅度不影响光谱形状。hsdar提供了sam函数直接把Speclib和参考光谱丢进去就能算ref - avgsp(cr_spec[SI(cr_spec)$group health, ]) angles - sam(cr_spec, ref) hist(angles, breaks 50, main SAM角度分布)avgsp是hsdar的平均函数用于取一组光谱的平均值作为参考光谱这里取健康样本的平均光谱。sam返回每个样本与参考光谱夹角的余弦值转换后就是角度。看分布直方图如果两个类别的角度分布有明显分界直接用阈值就能分类。经验值上0.05~0.1弧度内的光谱可以视为同类0.2以上基本是不同地物。这个阈值对模拟数据很理想实测数据因为有噪声会模糊一些可以先画分布再定阈值。6.2 给分类结果上“后悔药”caret五折交叉验证SAM阈值分类快但阈值怎么定总带点玄学。要对分类结果有把握正规做法还是得跑一遍交叉验证。用caret配randomForest五折交叉验证就能给出客观的准确率估计。有了这个数字你再去定SAM阈值心里就有底了。df - as.data.frame(spectra(cr_spec)) df$group - SI(cr_spec)$group ctrl - trainControl(method cv, number 5) model - train(group ~ ., data df, method rf, trControl ctrl, tuneLength 5) print(model$results)train的group ~ .表示用所有光谱波段预测分组methodrf指定随机森林tuneLength5会自动试5组mtry参数。model$results里能看到不同mtry下的准确率选最高的一行作为最终参数。把这套结果和SAM分类图放在一起看互相验证比单靠任何一种方法都踏实。6.3 出图规范ggplot导出高分辨率jpg和pdf分析做完了最终要落到一张能放进报告里的图。R语言里导出图片用ggplot加ggsave能同时出jpg和pdf两个版本jpg用于快速预览和微信传输pdf用于论文投稿和打印。library(ggplot2) sam_df - data.frame(x rep(1:50, 50), y rep(1:50, each 50), angle as.vector(angles)) p - ggplot(sam_df, aes(x x, y y, fill angle)) geom_raster() scale_fill_viridis_c() theme_minimal() ggsave(sam_surface.pdf, p, width 8, height 6) ggsave(sam_surface.jpg, p, width 8, height 6, dpi 300)dpi300是印刷级分辨率务必写清楚否则默认的72dpi导出的图放大全是马赛克。pdf不用指定dpi它是矢量格式缩放不丢细节。我现在的习惯是每次出图都同时存pdf和jpg避免投稿时临时找不到矢量图再跑一遍数据。整个流程走下来你会发现R这套开源方案真正值钱的不是某一个算法而是整个链路可复现原始光谱进分类图出每一步的参数和代码都在换一批数据也能立刻重跑。这种踏实感是商业软件给不了的。希望帮到你。本文还有配套的精品资源点击获取