简介光谱仪采集到的SVC PSR数据在预处理时常会遇到文件读取格式不统一、平滑参数不好确定、重采样操作较为繁琐等问题。这套使用MATLAB编写的源码正是针对这些环节进行整理适合刚入门的新手也适合有一定经验的开发人员直接改写使用。压缩包内共有两个文件均为点M文件整个压缩包只有两KB结构非常简单没有多余内容。脚本实现了测量数据平均批处理以及光谱重采样两个主要功能可以帮助用户完成多条光谱曲线的批量载入、平滑去噪、按波段重采样以及结果输出大幅减少手动操作时间。作者已对源码进行过实测校正可以稳定运行。目前这个资源已经吸引四百五十四人浏览学习对于遥感、地物光谱测量、实验室光谱分析等方向是一份小巧实用的参考代码。1. SVC PSR 光谱数据读进来第一步就卡在格式上SVC PSR 系列地物光谱仪在野外测量中很常见仪器输出的 .SIG 文件在配套软件里看很正常但要放进 MATLAB 做平滑、重采样和批处理很多人第一步就被挡住了它不是纯文本load和importdata都读不出有效数据直接fread又不知道每个字节代表什么。我最早处理 PSR-3500 数据时也在这个环节上耗了半天后来才把读入流程稳定下来。这篇文章把整条处理链路拆开讲SIG 二进制数据的读入、光谱平滑、重采样到统一波长轴、多文件批处理平均最后给出几个验证脚本有没有跑对的具体方法。适合做植被遥感、土壤光谱和实验室光谱建模的工程师与研究生也适合想摆脱手工 Excel 处理流程、把光谱预处理脚本化的实验室用户。2. SVC PSR SIG 文件读入二进制解析与波长轴重建2.1 SIG 文件结构与读入前的基本判断SVC PSR 的 SIG 文件是文本与二进制混合结构。文件开头是一段可读的 ASCII 头信息记录仪器型号、测量时间、积分时间、增益设置、GPS 坐标等数据区紧随其后一般按float32或float64存放每个波段的响应值。问题在于不同固件版本头部长度不固定有的版本固定 512 字节有的多了扩展配置行导致直接按固定偏移量fread时数据错位读出来要么数量级不对要么波长序列和数值对不上。处理这类文件我建议用先读字节、再找标志的方式而不是写死偏移量。思路是fopen打开后先读一段固定大小的字节转成char后搜索数据区起始标志确定偏移量后再fseek到正确位置读正式数据。下面这个函数框架就是按这个思路写的可以直接套用function [wvl, spc] read_psr_sig(sigfile) % READ_PSR_SIG 读取SVC PSR系列SIG格式光谱文件 % 输出: % wvl - 波长向量(nm) % spc - 光谱值向量 fid fopen(sigfile, rb); if fid -1 error(无法读取文件: %s, sigfile); end % 读入前2048字节, 从中定位数据区起点 headbuf fread(fid, 2048, uint8char); headArea headbuf(1:2048); % 常见做法: 数据区紧跟在包含 DATA 关键字的行之后 tok regexp(headArea, DATA\s*[:]?\s*(\d), tokens, once); if ~isempty(tok) dataStart str2double(tok{1}); else % 兼容没有数值标志的情况: 按关键字出现位置估算 dataStart strfind(headArea, DATA); if isempty(dataStart) error(无法定位数据区起始位置); end dataStart dataStart(end) 4; end fseek(fid, dataStart, bof); % 波段数优先从头部读取, 读不到再用默认值 bandTok regexp(headArea, [Bb]ands\s*[:]?\s*(\d), tokens, once); if ~isempty(bandTok) nbands str2double(bandTok{1}); else nbands 1024; % 兜底值, 请按仪器型号确认 end % 光谱值通常为float32, 读取后转置为行向量 spc fread(fid, nbands, float32); fclose(fid); % 波长轴单独处理, 见2.3节说明 wvl build_wavelength_axis(sigfile, headArea); end逻辑说明regexp去头部文本里搜索Bands字段比直接写死波段数安全。DATA标志在不同固件下的写法有差异所以用正则兼容DATA:与DATA两种格式。如果文件里没有数值标志就回退到字符串定位这种兜底逻辑在野外采集数据里很管用。参数说明fread的uint8char是把字节按无符号整型读出后直接映射为字符避免扩展 ASCII 字符导致乱码。读取光谱数据用float32如果读出来全是极小值或极大值改成float64再试一次。nbands的默认值只是兜底实际波段数以头部字段为准。2.2 头部字段如何影响后续处理参数头部信息里除了能确定数据区起点还能为后续步骤提供依据。比如积分时间长的数据信噪比高平滑窗口可以取小一点增益设置决定量程读数超出预期范围时能快速判断是不是读错了位置。下表是几个关键字段头部字段含义对处理流程的影响MODEL仪器型号决定波段数和波长范围的基本参考Bands通道数传给fread的读取长度IntegrationTime积分时间评估噪声水平影响平滑窗口大小Gain增益设置判断数值量程异常值排查时使用DATA数据区起始位置用于fseek定位偏移量实际项目中我一般会在读完头部后把整段headArea保存到日志文件里方便以后排查。野外数据量大时这样做的价值会体现出来某个文件读出来波形异常翻日志能直接看到当时的积分时间、增益和文件位置不用重新返工。2.3 波长轴的来源与同型号不通用问题波长轴处理是整个读入过程最容易被低估的环节。PSR 系列采用棱镜色散加阵列探测器波长与像元位置不是线性关系而且每台仪器出厂定标都有微小差异同一通道在不同机器上的中心波长可能偏移 0.3~1.5nm。如果直接拿另一台机器的波长表来用后续的重采样、吸收特征位置对比都会引入系统性偏差。我常用的做法是分两步第一次拿到数据时把 SIG 文件与仪器配套软件导出的 ASCII 文本做一次波长匹配将波长轴保存成chs.mat之后处理同批次数据直接加载不需要重复对波长。如果仪器送修或更换过 CCD再重新做一次定标匹配。% 首次处理: 从仪器软件导出的ASCII文件建立波长轴 wvl_export load(psr3500_wavelength.txt); % 第一列序号, 第二列波长 save(chs.mat, wvl_export); % 后续处理: 直接加载保存的定标波长 S load(chs.mat); wvl S.wvl_export(:, 2);这样既避免了每次处理都重复对波长也不会因为 SIG 文件里没有波长信息而把流程卡住。3. 光谱平滑的工程选择Savitzky-Golay 参数与边界控制3.1 移动平均为什么不适合地物光谱对光谱做平滑最朴素的思路是移动平均即窗口内所有点等权平均。PSR 光谱的分辨率本身有限移动平均的代价是会把突变的吸收谷整体拉低峰谷被抻平窄吸收特征被吞掉。高光谱数据预处理里更常用的是 Savitzky-Golay 卷积平滑在一个滑动窗口内用最小二乘多项式拟合用拟合值替代中心点相当于一个能按信号局部形态调整权重的低通滤波器。一句话区别移动平均是纯低通Savitzky-Golay 是保留局部形态的低通。后面做吸收特征提取、植被指数计算时保形性比看起来平滑重要得多。3.2 sgolayfilt 的窗口、阶数与判断依据MATLAB 中的调用格式是y sgolayfilt(x, order, framelen)。order是拟合多项式阶数framelen是滑动窗口长度必须是奇数且大于order1。这两个参数的配合主要看两点波段间隔和噪声水平。使用场景orderframelen说明一般野外测量, 噪声较低25~9平滑适中, 峰形保留好噪声较高的老型号设备2~311~15平滑强, 注意过平滑风险吸收特征很窄的矿物光谱35尽量保留窄吸收峰批量建模前的统一预处理27兼顾速度与稳定性表格里的范围是我处理 PSR-1100 和 PSR-3500 数据时的经验值并非严格规定。判断窗口是否过大一个有效的定量指标是看一阶差分信号的标准差变化。原始光谱的一阶差分包含噪声和吸收特征引起的高频变化平滑后差分标准差会下降如果降幅过大说明真实特征也被压掉了% 计算一阶差分, 评估平滑对高频分量的影响 d_orig diff(x); % 原始数据一阶差分 x_s sgolayfilt(x, 2, 9); d_s diff(x_s); ratio std(d_s) / std(d_orig); fprintf(平滑后高频分量剩余 %.1f%%\n, ratio * 100); if ratio 0.15 warning(高频信息损失过多, 建议减小窗口); end逻辑说明一阶差分反映相邻波段间的变化率随机噪声在差分后标准差与噪声水平成正比。平滑后比值低于 15% 时窗口往往已经把一部分真实吸收特征平掉了需要调小窗口或降低阶数。3.3 边界效应与低信噪比波段的掩膜处理sgolayfilt对序列两端采用非对称窗口补全边缘几个点会出现明显失真。PSR 数据在 350~400nm 和 2400~2500nm 两端信噪比很低平滑这两段只会产出假纹理没有实际物理意义。常见做法不是对整个波长范围一刀切平滑而是先按可信波段范围做掩膜% 按物理上可信的波段范围做掩膜 valid (wvl 400) (wvl 2400); wvl_v wvl(valid); x_v x(valid); % 平滑只在有效波段内进行 x_smooth_v sgolayfilt(x_v, 2, 7); % 填回整谱, 无效波段置NaN x_smooth NaN(size(x)); x_smooth(valid) x_smooth_v; % 后续计算用 omitnan 忽略无效区域 avg_refl mean(x_smooth, omitnan);掩膜处理后两端噪声大的波段不会拖累全谱统计量后面做平均和多文件对比时也更稳。4. 光谱重采样与插值方法选择重采样_PSR.m 的核心思路4.1 为什么要做重采样什么时候必须做SVC PSR 的波长间隔本身是非均匀的可见光波段可能 1.5nm 一个点到短波红外区间可能拉开到 6~8nm。另外野外连续测量时环境温漂可能使同一台仪器的波长刻度发生轻微偏移。所以跨文件平均、并列对比、放入建模框架之前通常需要把所有光谱重采样到统一的等间隔波长网格比如 1nm 步长。4.2 重采样_PSR.m 的四个核心步骤这种脚本的核心逻辑可以拆解成四步确定目标波长轴清理原始波长轴上的重复点和 NaN用插值函数在目标网格上计算光谱值把超出原始范围的位置置 NaN。function [wvl_new, spc_new] resample_psr(wvl_orig, spc_orig, step) % 重采样_PSR.m 核心流程 if nargin 3, step 1; end % 第一步: 定义统一目标波长轴 lo ceil(min(wvl_orig)); hi floor(max(wvl_orig)); wvl_new (lo : step : hi); % 第二步: 排序并去重, interp1要求x单调递增且不重复 [wvl_orig, idx] sort(wvl_orig); spc_orig spc_orig(idx); dup [false; diff(wvl_orig) 0]; wvl_orig(dup) []; spc_orig(dup) []; % 第三步: 插值, linear最稳, spline平滑但可能引入负值 spc_new interp1(wvl_orig, spc_orig, wvl_new, linear, extrap); % 第四步: 超出原始覆盖范围的位置置NaN out (wvl_new min(wvl_orig)) | (wvl_new max(wvl_orig)); spc_new(out) NaN; end逻辑说明先排序再插值是因为interp1要求 x 单调递增二进制读入顺序如果反了这一步能避免运行期报错。去重步骤处理的是同一波长位置上出现两个测量值的偶发情况不去除会导致interp1报网格不唯一的错。参数说明linear对含噪光谱最稳spline出来的曲线更顺滑但在吸收谷附近可能插出负反射率物理上不合理。pchip保持单调性不会像 spline 那样振荡不过计算量大一些。批量重采样时我默认用linear只有需要制图平滑效果时才换pchip。4.3 插值方法对比与步长选择插值方法平滑度主要风险适用场景linear一般保留原始噪声批量处理默认选择spline高可能产生负值或振荡仅用于可视化pchip高计算量较大要求严格单调且保形时步长方面1nm 是很多地物光谱分析框架的标准网格。如果后续要与多光谱卫星波段响应函数做卷积模拟可能需要更密的 0.5nm 网格但步长设得比原始分辨率更细不会增加真实信息量只是让文件体积变大。5. 批处理与平均光谱把单文件脚本串成自动化流程5.1 目录遍历与文件名排序的细节批处理第一步是确定读哪些文件、按什么顺序读。直接用dir(*.SIG)得到的结果按字节顺序排文件名带数字时可能出现SITE01_002.SIG排在SITE01_010.SIG后面的情况。常见做法是用 File Exchange 上的natsortfiles做自然排序或者用regexp从文件名提取序号再排序% 遍历目录下所有SIG文件 folder D:\psr_data\field0602; files dir(fullfile(folder, *.SIG)); n length(files); % 提取文件名中的数字序号, 按测量顺序排序 numIdx zeros(n, 1); for k 1 : n tok regexp(files(k).name, (\d), tokens, once); if ~isempty(tok) numIdx(k) str2double(tok{1}); end end [~, order] sort(numIdx); files files(order);逻辑说明把序号提取出来再排序比直接对字符串排序更符合野外测量顺序。如果同一个测点包含白板参考测量和目标测量文件名里一般有ref或target标识需要先分组再分别处理。参数说明regexp的tokens返回捕获组内容once只取第一个匹配。如果文件名里有多组数字例如日期和序号要按实际情况调整正则表达式。5.2 多条曲线平均前为什么要先重采样批处理流程里平均这一步有一个容易忽略的前提所有参与平均的光谱必须落在同一波长网格上。如果只在原始波长上直接平均各条光谱的波长偏移会被平均操作掩盖掉后期建模时又会被重新放大。所以正确顺序是读入 → 重采样到统一网格 → 平滑 → 平均。先平均再重采样是常见错误。平均时的异常剔除也很关键。野外测量常因云遮挡、镜头晃动或其他人为因素产生异常光谱逐条目检不现实。常用做法是计算每条光谱与中值光谱的偏差再按统计量判断% 每条光谱与中值光谱的绝对偏差 med_spc median(spectra, 2); % 逐波段中值 diff_mat abs(spectra - med_spc); total_dev sum(diff_mat, 1); % 每条光谱的总偏差 % 阈值: 中位数 3倍MAD threshold median(total_dev) 3 * mad(total_dev, 1); valid_idx total_dev threshold; % 剔除异常后平均 spectra_clean spectra(:, valid_idx); avg_spc mean(spectra_clean, 2, omitnan);逻辑说明异常光谱往往在个别波段出现强烈偏离比如云遮挡导致 1400nm 和 1900nm 水汽吸收带异常加深全波段总偏差会显著大于正常测量。用中位数加 3 倍 MAD 做阈值能降低异常值本身对阈值的干扰。5.3 批处理主脚本的组织方式与输出规范批处理脚本建议拆成独立函数读入、平滑、重采样、平均各一个文件主脚本只负责调度和输出。这样后期更换仪器型号或调整平滑参数时只改对应函数即可。输出时建议按处理日期建目录避免中间产物把原始数据目录搞乱输出文件内容用途processed.mat所有处理后的光谱矩阵供后续 MATLAB 建模average.csv平均光谱表报告配图或外部软件使用stats.mat每条光谱的标准差与异常标记便于复查6. 跑完脚本后必须做的几个验证从数值和波形上确认没跑偏6.1 与仪器软件导出的 ASCII 数据交叉验证读入脚本是否正确最直接的办法是把同一份 SIG 文件用 SVC 官方软件导出为 ASCII 文本再在 MATLAB 里与脚本读出的结果做比较。比较的内容包括三个值域范围是否一致吸收峰位置是否相同长度是否匹配。% 与官方导出的文本文件对比 ref load(exported_from_SVC.txt); % 第一列波长, 第二列反射率 plot(ref(:, 1), ref(:, 2), k-, wvl, spc, r--); legend(官方导出, 脚本读取);两者差别应只来自定标波长表的来源差异。如果趋势一致但整体有细微偏移多半是波长轴没对齐优先检查chs.mat里的波长表来源。6.2 重采样前后积分面积比检查重采样不应该大幅改变光谱的总能量。计算重采样前后在有效波长范围内的积分面积比值应接近 1。偏离超过 5% 时要检查插值方法是否引入了震荡或者目标波长轴范围是不是截多了% 原始与重采样光谱的积分面积比 area_orig trapz(wvl_orig_valid, spc_orig_valid); area_new trapz(wvl_new(~isnan(spc_new)), spc_new(~isnan(spc_new))); ratio area_new / area_orig; fprintf(重采样前后积分面积比: %.3f\n, ratio);trapz在非均匀网格下也能算但重采样前后用同一栅格更公平。如果比值偏差明显优先检查是否把原始范围外的 NaN 算进了积分。6.3 平滑效果的残差检查判断平滑是否把噪声去掉而没有伤及信号可以把平滑后的曲线从原始曲线中减去观察残差光谱。正常残差应该基本围绕零线随机分布且幅值与仪器噪声水平相当。如果残差里能看到明显的吸收峰或谷形说明窗口过大真实信号被当成噪声平滑掉了residual x - x_smooth; figure; plot(wvl, residual, b-); % 计算残差均方根, 与仪器噪声水平对比 rms_res rms(residual, omitnan); fprintf(残差均方根: %.4f\n, rms_res);残差里出现系统性峰谷时把窗口调小或降一阶重新平滑再重复这一步即可。本文还有配套的精品资源点击获取 SVC PSR光谱数据MATLAB处理：从SIG读取到平滑重采样与批处理
简介光谱仪采集到的SVC PSR数据在预处理时常会遇到文件读取格式不统一、平滑参数不好确定、重采样操作较为繁琐等问题。这套使用MATLAB编写的源码正是针对这些环节进行整理适合刚入门的新手也适合有一定经验的开发人员直接改写使用。压缩包内共有两个文件均为点M文件整个压缩包只有两KB结构非常简单没有多余内容。脚本实现了测量数据平均批处理以及光谱重采样两个主要功能可以帮助用户完成多条光谱曲线的批量载入、平滑去噪、按波段重采样以及结果输出大幅减少手动操作时间。作者已对源码进行过实测校正可以稳定运行。目前这个资源已经吸引四百五十四人浏览学习对于遥感、地物光谱测量、实验室光谱分析等方向是一份小巧实用的参考代码。1. SVC PSR 光谱数据读进来第一步就卡在格式上SVC PSR 系列地物光谱仪在野外测量中很常见仪器输出的 .SIG 文件在配套软件里看很正常但要放进 MATLAB 做平滑、重采样和批处理很多人第一步就被挡住了它不是纯文本load和importdata都读不出有效数据直接fread又不知道每个字节代表什么。我最早处理 PSR-3500 数据时也在这个环节上耗了半天后来才把读入流程稳定下来。这篇文章把整条处理链路拆开讲SIG 二进制数据的读入、光谱平滑、重采样到统一波长轴、多文件批处理平均最后给出几个验证脚本有没有跑对的具体方法。适合做植被遥感、土壤光谱和实验室光谱建模的工程师与研究生也适合想摆脱手工 Excel 处理流程、把光谱预处理脚本化的实验室用户。2. SVC PSR SIG 文件读入二进制解析与波长轴重建2.1 SIG 文件结构与读入前的基本判断SVC PSR 的 SIG 文件是文本与二进制混合结构。文件开头是一段可读的 ASCII 头信息记录仪器型号、测量时间、积分时间、增益设置、GPS 坐标等数据区紧随其后一般按float32或float64存放每个波段的响应值。问题在于不同固件版本头部长度不固定有的版本固定 512 字节有的多了扩展配置行导致直接按固定偏移量fread时数据错位读出来要么数量级不对要么波长序列和数值对不上。处理这类文件我建议用先读字节、再找标志的方式而不是写死偏移量。思路是fopen打开后先读一段固定大小的字节转成char后搜索数据区起始标志确定偏移量后再fseek到正确位置读正式数据。下面这个函数框架就是按这个思路写的可以直接套用function [wvl, spc] read_psr_sig(sigfile) % READ_PSR_SIG 读取SVC PSR系列SIG格式光谱文件 % 输出: % wvl - 波长向量(nm) % spc - 光谱值向量 fid fopen(sigfile, rb); if fid -1 error(无法读取文件: %s, sigfile); end % 读入前2048字节, 从中定位数据区起点 headbuf fread(fid, 2048, uint8char); headArea headbuf(1:2048); % 常见做法: 数据区紧跟在包含 DATA 关键字的行之后 tok regexp(headArea, DATA\s*[:]?\s*(\d), tokens, once); if ~isempty(tok) dataStart str2double(tok{1}); else % 兼容没有数值标志的情况: 按关键字出现位置估算 dataStart strfind(headArea, DATA); if isempty(dataStart) error(无法定位数据区起始位置); end dataStart dataStart(end) 4; end fseek(fid, dataStart, bof); % 波段数优先从头部读取, 读不到再用默认值 bandTok regexp(headArea, [Bb]ands\s*[:]?\s*(\d), tokens, once); if ~isempty(bandTok) nbands str2double(bandTok{1}); else nbands 1024; % 兜底值, 请按仪器型号确认 end % 光谱值通常为float32, 读取后转置为行向量 spc fread(fid, nbands, float32); fclose(fid); % 波长轴单独处理, 见2.3节说明 wvl build_wavelength_axis(sigfile, headArea); end逻辑说明regexp去头部文本里搜索Bands字段比直接写死波段数安全。DATA标志在不同固件下的写法有差异所以用正则兼容DATA:与DATA两种格式。如果文件里没有数值标志就回退到字符串定位这种兜底逻辑在野外采集数据里很管用。参数说明fread的uint8char是把字节按无符号整型读出后直接映射为字符避免扩展 ASCII 字符导致乱码。读取光谱数据用float32如果读出来全是极小值或极大值改成float64再试一次。nbands的默认值只是兜底实际波段数以头部字段为准。2.2 头部字段如何影响后续处理参数头部信息里除了能确定数据区起点还能为后续步骤提供依据。比如积分时间长的数据信噪比高平滑窗口可以取小一点增益设置决定量程读数超出预期范围时能快速判断是不是读错了位置。下表是几个关键字段头部字段含义对处理流程的影响MODEL仪器型号决定波段数和波长范围的基本参考Bands通道数传给fread的读取长度IntegrationTime积分时间评估噪声水平影响平滑窗口大小Gain增益设置判断数值量程异常值排查时使用DATA数据区起始位置用于fseek定位偏移量实际项目中我一般会在读完头部后把整段headArea保存到日志文件里方便以后排查。野外数据量大时这样做的价值会体现出来某个文件读出来波形异常翻日志能直接看到当时的积分时间、增益和文件位置不用重新返工。2.3 波长轴的来源与同型号不通用问题波长轴处理是整个读入过程最容易被低估的环节。PSR 系列采用棱镜色散加阵列探测器波长与像元位置不是线性关系而且每台仪器出厂定标都有微小差异同一通道在不同机器上的中心波长可能偏移 0.3~1.5nm。如果直接拿另一台机器的波长表来用后续的重采样、吸收特征位置对比都会引入系统性偏差。我常用的做法是分两步第一次拿到数据时把 SIG 文件与仪器配套软件导出的 ASCII 文本做一次波长匹配将波长轴保存成chs.mat之后处理同批次数据直接加载不需要重复对波长。如果仪器送修或更换过 CCD再重新做一次定标匹配。% 首次处理: 从仪器软件导出的ASCII文件建立波长轴 wvl_export load(psr3500_wavelength.txt); % 第一列序号, 第二列波长 save(chs.mat, wvl_export); % 后续处理: 直接加载保存的定标波长 S load(chs.mat); wvl S.wvl_export(:, 2);这样既避免了每次处理都重复对波长也不会因为 SIG 文件里没有波长信息而把流程卡住。3. 光谱平滑的工程选择Savitzky-Golay 参数与边界控制3.1 移动平均为什么不适合地物光谱对光谱做平滑最朴素的思路是移动平均即窗口内所有点等权平均。PSR 光谱的分辨率本身有限移动平均的代价是会把突变的吸收谷整体拉低峰谷被抻平窄吸收特征被吞掉。高光谱数据预处理里更常用的是 Savitzky-Golay 卷积平滑在一个滑动窗口内用最小二乘多项式拟合用拟合值替代中心点相当于一个能按信号局部形态调整权重的低通滤波器。一句话区别移动平均是纯低通Savitzky-Golay 是保留局部形态的低通。后面做吸收特征提取、植被指数计算时保形性比看起来平滑重要得多。3.2 sgolayfilt 的窗口、阶数与判断依据MATLAB 中的调用格式是y sgolayfilt(x, order, framelen)。order是拟合多项式阶数framelen是滑动窗口长度必须是奇数且大于order1。这两个参数的配合主要看两点波段间隔和噪声水平。使用场景orderframelen说明一般野外测量, 噪声较低25~9平滑适中, 峰形保留好噪声较高的老型号设备2~311~15平滑强, 注意过平滑风险吸收特征很窄的矿物光谱35尽量保留窄吸收峰批量建模前的统一预处理27兼顾速度与稳定性表格里的范围是我处理 PSR-1100 和 PSR-3500 数据时的经验值并非严格规定。判断窗口是否过大一个有效的定量指标是看一阶差分信号的标准差变化。原始光谱的一阶差分包含噪声和吸收特征引起的高频变化平滑后差分标准差会下降如果降幅过大说明真实特征也被压掉了% 计算一阶差分, 评估平滑对高频分量的影响 d_orig diff(x); % 原始数据一阶差分 x_s sgolayfilt(x, 2, 9); d_s diff(x_s); ratio std(d_s) / std(d_orig); fprintf(平滑后高频分量剩余 %.1f%%\n, ratio * 100); if ratio 0.15 warning(高频信息损失过多, 建议减小窗口); end逻辑说明一阶差分反映相邻波段间的变化率随机噪声在差分后标准差与噪声水平成正比。平滑后比值低于 15% 时窗口往往已经把一部分真实吸收特征平掉了需要调小窗口或降低阶数。3.3 边界效应与低信噪比波段的掩膜处理sgolayfilt对序列两端采用非对称窗口补全边缘几个点会出现明显失真。PSR 数据在 350~400nm 和 2400~2500nm 两端信噪比很低平滑这两段只会产出假纹理没有实际物理意义。常见做法不是对整个波长范围一刀切平滑而是先按可信波段范围做掩膜% 按物理上可信的波段范围做掩膜 valid (wvl 400) (wvl 2400); wvl_v wvl(valid); x_v x(valid); % 平滑只在有效波段内进行 x_smooth_v sgolayfilt(x_v, 2, 7); % 填回整谱, 无效波段置NaN x_smooth NaN(size(x)); x_smooth(valid) x_smooth_v; % 后续计算用 omitnan 忽略无效区域 avg_refl mean(x_smooth, omitnan);掩膜处理后两端噪声大的波段不会拖累全谱统计量后面做平均和多文件对比时也更稳。4. 光谱重采样与插值方法选择重采样_PSR.m 的核心思路4.1 为什么要做重采样什么时候必须做SVC PSR 的波长间隔本身是非均匀的可见光波段可能 1.5nm 一个点到短波红外区间可能拉开到 6~8nm。另外野外连续测量时环境温漂可能使同一台仪器的波长刻度发生轻微偏移。所以跨文件平均、并列对比、放入建模框架之前通常需要把所有光谱重采样到统一的等间隔波长网格比如 1nm 步长。4.2 重采样_PSR.m 的四个核心步骤这种脚本的核心逻辑可以拆解成四步确定目标波长轴清理原始波长轴上的重复点和 NaN用插值函数在目标网格上计算光谱值把超出原始范围的位置置 NaN。function [wvl_new, spc_new] resample_psr(wvl_orig, spc_orig, step) % 重采样_PSR.m 核心流程 if nargin 3, step 1; end % 第一步: 定义统一目标波长轴 lo ceil(min(wvl_orig)); hi floor(max(wvl_orig)); wvl_new (lo : step : hi); % 第二步: 排序并去重, interp1要求x单调递增且不重复 [wvl_orig, idx] sort(wvl_orig); spc_orig spc_orig(idx); dup [false; diff(wvl_orig) 0]; wvl_orig(dup) []; spc_orig(dup) []; % 第三步: 插值, linear最稳, spline平滑但可能引入负值 spc_new interp1(wvl_orig, spc_orig, wvl_new, linear, extrap); % 第四步: 超出原始覆盖范围的位置置NaN out (wvl_new min(wvl_orig)) | (wvl_new max(wvl_orig)); spc_new(out) NaN; end逻辑说明先排序再插值是因为interp1要求 x 单调递增二进制读入顺序如果反了这一步能避免运行期报错。去重步骤处理的是同一波长位置上出现两个测量值的偶发情况不去除会导致interp1报网格不唯一的错。参数说明linear对含噪光谱最稳spline出来的曲线更顺滑但在吸收谷附近可能插出负反射率物理上不合理。pchip保持单调性不会像 spline 那样振荡不过计算量大一些。批量重采样时我默认用linear只有需要制图平滑效果时才换pchip。4.3 插值方法对比与步长选择插值方法平滑度主要风险适用场景linear一般保留原始噪声批量处理默认选择spline高可能产生负值或振荡仅用于可视化pchip高计算量较大要求严格单调且保形时步长方面1nm 是很多地物光谱分析框架的标准网格。如果后续要与多光谱卫星波段响应函数做卷积模拟可能需要更密的 0.5nm 网格但步长设得比原始分辨率更细不会增加真实信息量只是让文件体积变大。5. 批处理与平均光谱把单文件脚本串成自动化流程5.1 目录遍历与文件名排序的细节批处理第一步是确定读哪些文件、按什么顺序读。直接用dir(*.SIG)得到的结果按字节顺序排文件名带数字时可能出现SITE01_002.SIG排在SITE01_010.SIG后面的情况。常见做法是用 File Exchange 上的natsortfiles做自然排序或者用regexp从文件名提取序号再排序% 遍历目录下所有SIG文件 folder D:\psr_data\field0602; files dir(fullfile(folder, *.SIG)); n length(files); % 提取文件名中的数字序号, 按测量顺序排序 numIdx zeros(n, 1); for k 1 : n tok regexp(files(k).name, (\d), tokens, once); if ~isempty(tok) numIdx(k) str2double(tok{1}); end end [~, order] sort(numIdx); files files(order);逻辑说明把序号提取出来再排序比直接对字符串排序更符合野外测量顺序。如果同一个测点包含白板参考测量和目标测量文件名里一般有ref或target标识需要先分组再分别处理。参数说明regexp的tokens返回捕获组内容once只取第一个匹配。如果文件名里有多组数字例如日期和序号要按实际情况调整正则表达式。5.2 多条曲线平均前为什么要先重采样批处理流程里平均这一步有一个容易忽略的前提所有参与平均的光谱必须落在同一波长网格上。如果只在原始波长上直接平均各条光谱的波长偏移会被平均操作掩盖掉后期建模时又会被重新放大。所以正确顺序是读入 → 重采样到统一网格 → 平滑 → 平均。先平均再重采样是常见错误。平均时的异常剔除也很关键。野外测量常因云遮挡、镜头晃动或其他人为因素产生异常光谱逐条目检不现实。常用做法是计算每条光谱与中值光谱的偏差再按统计量判断% 每条光谱与中值光谱的绝对偏差 med_spc median(spectra, 2); % 逐波段中值 diff_mat abs(spectra - med_spc); total_dev sum(diff_mat, 1); % 每条光谱的总偏差 % 阈值: 中位数 3倍MAD threshold median(total_dev) 3 * mad(total_dev, 1); valid_idx total_dev threshold; % 剔除异常后平均 spectra_clean spectra(:, valid_idx); avg_spc mean(spectra_clean, 2, omitnan);逻辑说明异常光谱往往在个别波段出现强烈偏离比如云遮挡导致 1400nm 和 1900nm 水汽吸收带异常加深全波段总偏差会显著大于正常测量。用中位数加 3 倍 MAD 做阈值能降低异常值本身对阈值的干扰。5.3 批处理主脚本的组织方式与输出规范批处理脚本建议拆成独立函数读入、平滑、重采样、平均各一个文件主脚本只负责调度和输出。这样后期更换仪器型号或调整平滑参数时只改对应函数即可。输出时建议按处理日期建目录避免中间产物把原始数据目录搞乱输出文件内容用途processed.mat所有处理后的光谱矩阵供后续 MATLAB 建模average.csv平均光谱表报告配图或外部软件使用stats.mat每条光谱的标准差与异常标记便于复查6. 跑完脚本后必须做的几个验证从数值和波形上确认没跑偏6.1 与仪器软件导出的 ASCII 数据交叉验证读入脚本是否正确最直接的办法是把同一份 SIG 文件用 SVC 官方软件导出为 ASCII 文本再在 MATLAB 里与脚本读出的结果做比较。比较的内容包括三个值域范围是否一致吸收峰位置是否相同长度是否匹配。% 与官方导出的文本文件对比 ref load(exported_from_SVC.txt); % 第一列波长, 第二列反射率 plot(ref(:, 1), ref(:, 2), k-, wvl, spc, r--); legend(官方导出, 脚本读取);两者差别应只来自定标波长表的来源差异。如果趋势一致但整体有细微偏移多半是波长轴没对齐优先检查chs.mat里的波长表来源。6.2 重采样前后积分面积比检查重采样不应该大幅改变光谱的总能量。计算重采样前后在有效波长范围内的积分面积比值应接近 1。偏离超过 5% 时要检查插值方法是否引入了震荡或者目标波长轴范围是不是截多了% 原始与重采样光谱的积分面积比 area_orig trapz(wvl_orig_valid, spc_orig_valid); area_new trapz(wvl_new(~isnan(spc_new)), spc_new(~isnan(spc_new))); ratio area_new / area_orig; fprintf(重采样前后积分面积比: %.3f\n, ratio);trapz在非均匀网格下也能算但重采样前后用同一栅格更公平。如果比值偏差明显优先检查是否把原始范围外的 NaN 算进了积分。6.3 平滑效果的残差检查判断平滑是否把噪声去掉而没有伤及信号可以把平滑后的曲线从原始曲线中减去观察残差光谱。正常残差应该基本围绕零线随机分布且幅值与仪器噪声水平相当。如果残差里能看到明显的吸收峰或谷形说明窗口过大真实信号被当成噪声平滑掉了residual x - x_smooth; figure; plot(wvl, residual, b-); % 计算残差均方根, 与仪器噪声水平对比 rms_res rms(residual, omitnan); fprintf(残差均方根: %.4f\n, rms_res);残差里出现系统性峰谷时把窗口调小或降一阶重新平滑再重复这一步即可。本文还有配套的精品资源点击获取