简介针对SVC PSR光谱数据的处理需求这套MATLAB源码实现了数据读入、光谱平滑、重采样与测量数据平均批处理等核心功能面向遥感、地物光谱分析领域的新手及有一定经验的开发人员。压缩包内共2个.m脚本整体大小仅2KB体量虽小但代码精炼覆盖了从单文件读入到批量处理的完整链条便于学习、复用与二次开发。脚本命名直观与处理步骤对应便于定位功能模块。已有454人学习下载源码经过校正测试可确保直接运行。读者从中既可以看到PSR数据读取与预处理的实现思路也能借助平均批处理脚本快速整理多组光谱测量结果对于正在搭建光谱预处理流程的科研场景这份资源既能辅助理解算法细节也可作为实用工具直接嵌入现有项目减少重复编码提升效率。1. SVC PSR 光谱数据读入别让第一行代码卡在文件格式上野外光谱仪 SVC PSR 系列在植被遥感、土壤属性反演、矿物填图里很常见拿到的数据多半是.sig文本文件。很多人在 MATLAB 里直接用load或者textread读结果报错或者读进来全是乱的原因在于.sig不是纯两列数据而是按「暗电流段、白参考段、目标光谱段」分段存储的 ASCII 文件不同段的数据长度和含义完全不一样。这篇文章就围绕「读入、平滑、重采样、批处理」这四步把 SVC PSR 数据在 MATLAB 里的完整处理路径讲清楚。适合刚接触光谱数据处理的研究生也适合要把上千条光谱批量整理成训练集的工程师。读完你能自己写一个不依赖 SVC 自带后处理软件的 MATLAB 流水线把.sig文件变成干净、平滑、波长统一的光谱矩阵。2. 用 MATLAB 解析 SVC PSR 的 .sig 文件段标记识别与三种读入方案2.1 先看懂 .sig 文件的段结构HR 光谱段、DR 暗电流段、WR 白参考段SVC PSR 系列输出的.sig文件多数情况下是制表符或空格分隔的文本文件里除了文件头信息还会出现若干段标记。常见的段标记包括HR目标光谱段、DR暗电流段、WR白参考段每个标记行后面跟着波长和 DN 值两列数据。不同版本仪器写出的文件略有差异但「标记行 数据块」的结构基本一致。在写读入代码之前先用任意文本编辑器打开一个.sig文件观察前 20 行。重点确认三件事段标记是独立一行还是和数据在同一行波长列和 DN 列之间的分隔符是制表符、逗号还是空格每段数据的起始波长和结束波长是否一致。这些信息直接决定用textscan还是readtable。我处理过的 SVC PSR 文件里大多是独立标记行加两列数值但也碰到过把反射率和 DN 值混在一行输出的导出文件所以读入逻辑要写成「按标记行切换状态」而不是「固定行号读取」。% 查看文件前 20 行确认段标记格式 fid fopen(sample.sig, r); for i 1:20 tline fgetl(fid); if ~ischar(tline), break; end fprintf(%03d: %s\n, i, tline); end fclose(fid);这段代码先把文件头结构摸清楚避免后续解析时踩「标记行没识别到」的坑。fgetl按行读取返回字符向量如果返回 -1 说明已经到文件末尾。fprintf里的%s直接打印当前行内容方便肉眼核对。2.2 方案一按标记行分段读入的通用解析函数确认段标记后写一个按状态机方式读入的解析函数。核心思路是读到HR、DR、WR等标记时切换当前段读到数值行时把数据追加到当前段的数组里。function T read_sig_file(filename) % 读取 SVC PSR .sig 文件返回包含波长和 DN 值的 table fid fopen(filename, r); if fid -1 error(无法打开文件: %s, filename); end wl []; dn []; % 当前段的数据缓存 sections struct(); % 存放各段数据 while ~feof(fid) tline fgetl(fid); if ~ischar(tline) break; end tline strtrim(tline); % 去掉首尾空格 if isempty(tline) continue; end % 判断是否为段标记行 if startsWith(tline, HR) sections.HR begin_section(sections, HR); wl []; dn []; elseif startsWith(tline, DR) sections.DR begin_section(sections, DR); wl []; dn []; elseif startsWith(tline, WR) sections.WR begin_section(sections, WR); wl []; dn []; else nums sscanf(tline, %f); if numel(nums) 2 wl(end1) nums(1); dn(end1) nums(2); end end end % 把最后一段写入结构体 if ~isempty(wl) sections.(current_section) [wl(:), dn(:)]; end fclose(fid); T struct2table(sections); end function s begin_section(s, secname) % 上一段数据尚未保存时先写入 if isfield(s, secname) s.(secname) s.(secname); end end这段代码里需要注意几个细节。startsWith是 R2016b 之后才有的函数如果你的 MATLAB 版本较老改成strncmp(tline, HR, 2) 1。sscanf(tline, %f)会一次性解析一行里的所有数值SVC 的.sig数据行通常只有波长和 DN 两个数但有些文件会带第三列标记numel(nums) 2的判断能兼容这种情况。begin_section函数写得比较保守实际应用中你可以在切换段标记时直接把上一段的[wl(:), dn(:)]存入sections这段代码只是示范「数据结构如何切换」。真实场景里一定要处理「遇到标记行时先把上一段缓存存起来」否则最后一段数据会丢失。2.3 方案二用 readtable 省事的适用场景如果你的 SVC 软件已经把.sig导出成了标准 CSV 或 TXT第一行是列名后面是波长和反射率两列那么不需要自己解析readtable一步到位。% 读取已导出的标准格式光谱文件 opts detectImportOptions(reflectance.csv); opts.VariableNamingRule preserve; T readtable(reflectance.csv, opts); % 检查列名和数据类型 disp(T.Properties.VariableNames); disp(class(T.Wavelength));detectImportOptions会自动推断分隔符、变量名和数据格式适合列结构固定的导出文件。注意VariableNamingRule要设为preserve否则 MATLAB 会把列名里的空格、横杠替换成下划线容易对不上原文件字段。这种方案不适用于原始.sig因为原始文件里三种段的数据行混在一起readtable会把暗电流和白参考数据也当光谱数据读进来。2.4 读入后的数据组织统一放进 table 还是分开存读进来的 SVC PSR 数据我一般推荐用table组织而不是分开的wl、dn两个向量。原因有两个第一批处理时需要把每条光谱的波长范围、点数、缺失值数量等元数据一起管理table的列式存储天然适合做这种「每条记录一行」的结构第二后面重采样和拼接矩阵时table可以直接用outerjoin或rows2vars完成宽表转换比手动维护索引省事很多。% 把 HR 段数据存成标准光谱表 specTable table(sections.HR(:,1), sections.HR(:,2), ... VariableNames, {Wavelength, Reflectance}); % 去掉 NaN 和负值 valid isfinite(specTable.Reflectance) specTable.Reflectance 0; specTable specTable(valid, :);这段代码先构造标准两列表再用isfinite配合逻辑索引过滤异常值。SVC PSR 在波段边缘经常出现 NaN 或者负的 DN 值这属于正常现象但不能带着它们做平滑和重采样否则interp1会传播这些无效值。3. 光谱平滑与光谱重采样的参数陷阱窗口宽度、多项式阶数、插值核3.1 为什么原始 SVC PSR 光谱先平滑还是先重采样顺序不能乱SVC PSR 的原始波长间隔并不是均匀的VNIR 波段约 1.4 nmSWIR 波段约 2 nm 甚至更粗。直接对这样的非等间隔数据做 Savitzky-Golay 平滑MATLAB 的sgolayfilt默认按等间隔假设滤波结果会引入偏差。常见做法是先把原始光谱重采样到统一的 1 nm 网格再在这个网格上做平滑。反过来先平滑再重采样也不是不行但要自己写非等间隔的滑动窗口逻辑得不偿失。另外一个容易被忽略的问题重采样本身会改变噪声的统计特性。线性插值得出的光谱噪声方差比原始数据要小样条插值可能会放大局部的抖动。所以「平滑」和「重采样」不是两个独立的步骤而是一条流水线的两段。我一般按这个顺序走读入.sig→ 去 NaN/负值 → 重采样到 1 nm 网格 → Savitzky-Golay 平滑 → 导出。这个顺序写出来的批处理脚本对同一条光谱换不同参数重跑结果差异最小。3.2 Savitzky-Golay 平滑sgolayfilt 的窗口宽度与多项式阶数怎么配SG 平滑是光谱处理里最常用的方法MATLAB 里直接调sgolayfilt(x, k, f)k是多项式阶数f是窗口宽度必须为奇数且f k。参数选择直接影响平滑效果阶数越高保形越好去噪能力越差窗口越宽去噪越强但窄的吸收峰容易被抹平。处理目标多项式阶数 k窗口宽度 f适用场景强去噪不关心细节221 或 31土壤光谱建模关注整体趋势平衡去噪与保峰315 或 19植被光谱保留红边和叶绿素吸收精细结构分析49 或 11矿物光谱保留 2.2 μm 附近的小吸收峰% 在 1nm 等间隔网格上做 SG 平滑 k 3; % 多项式阶数 f 15; % 窗口宽度必须为奇数 smoothRefl sgolayfilt(specTable.Reflectance, k, f); % 对比平滑前后的差异 figure; plot(specTable.Wavelength, specTable.Reflectance, -, DisplayName, 原始); hold on; plot(specTable.Wavelength, smoothRefl, -, DisplayName, SG平滑); legend; xlabel(波长/nm); ylabel(反射率);参数里的关键约束是f必须为奇数且f k1。如果你设f 3、k 3MATLAB 会直接报错。记住一个经验公式窗口宽度约等于光谱中最小吸收峰半高宽的 23 倍太窄了起不到平滑作用太宽了把吸收峰抹成平台。对于 SVC PSR 的 1 nm 重采样光谱我一般从f15, k3起步如果平滑后 RMSE 变化小于仪器噪声水平说明窗口可以再大一点。3.3 用 interp1 把光谱重采样到统一波长网格重采样的核心是把 SVC PSR 的不规则波长映射到规则网格。最常用的函数是interp1但它有好几个插值核选错了会出问题。% 生成统一的 1nm 网格从 350nm 到 2500nm target_wl (350:1:2500); % 三种插值方法对比 refl_linear interp1(specTable.Wavelength, specTable.Reflectance, target_wl, linear); refl_pchip interp1(specTable.Wavelength, specTable.Reflectance, target_wl, pchip); refl_spline interp1(specTable.Wavelength, specTable.Reflectance, target_wl, spline); % 检查 spline 是否引入了负值 fprintf(spline 负值数量: %d\n, sum(refl_spline 0));linear是最稳妥的选择不会引入超过原始数据范围的过冲但会在波长间隔变化大的区域产生轻微的斜率不连续。pchip分段三次 Hermite 插值能保持单调性适合有陡峭吸收峰的光谱同时不会像spline那样在峰谷附近上下振铃。spline数学上最平滑但对于测量噪声大的光谱会在吸收峰两侧产生伪峰甚至负值。我给出的代码里专门检查负值数量因为反射率是物理量出现负值意味着插值方法引入了不可信的伪振荡。对于 SVC PSR 数据默认用pchip最稳。3.4 重采样网格怎么定起点、终点和步长网格设计决定了批处理后期所有样本是否对齐。很多人直接350:1:2500但这未必合理。SVC PSR 的实测范围通常到 2500 nm但有些仪器型号到 2400 nm 就截止了直接插值到 2500 会在尾部产生一段外推数据interp1默认返回 NaN如果没处理后面的矩阵拼接会报错。% 动态生成网格而不是硬编码 wl_min floor(min(specTable.Wavelength)); wl_max ceil(max(specTable.Wavelength)); target_wl (wl_min:1:wl_max);这段代码先计算当前样本的波长覆盖范围再生成网格。批处理时所有样本的wl_min和wl_max可能不一样统一的处理策略是取全体样本的交集。比如 100 条光谱里最短的只到 2450 nm那所有样本都只保留到 2450 nm避免矩阵尾部堆满 NaN。如果建模时用的是固定波段范围比如 4002400 nm直接截取这个范围内数据更干净。4. 文件批处理把读入、平滑、重采样串成一条 MATLAB 流水线4.1 用 dir 批量定位 .sig 文件而不是让用户一个个选取批处理的第一步是让 MATLAB 自己遍历文件夹里的所有.sig文件而不是用uigetfile手动点选。dir返回的文件列表里包含name、folder、bytes等字段配合fullfile拼接完整路径就可以在for循环里逐个处理。% 获取指定目录下所有 .sig 文件 data_dir D:\spectra\field_site_A; file_list dir(fullfile(data_dir, *.sig)); % 按文件名排序确保输出顺序稳定 [~, idx] sort({file_list.name}); file_list file_list(idx); fprintf(共找到 %d 个 SVC 文件\n, numel(file_list));dir(fullfile(data_dir, *.sig))返回的是 struct 数组每个元素对应一个文件。sort({file_list.name})是把文件名按字典序排列因为dir返回的顺序不保证后续如果按位置写日志会对不上号。4.2 从文件名提取站点和时间信息regexp 的常用匹配模式野外光谱的文件名往往自带信息比如SVC_PSR_20230615_0935_siteA.sig站点和时间都在文件名里。处理时把这些元数据提取出来保存到汇总表里后面按站点分组建模会很方便。% 匹配文件名中的日期、时间和站点标识 pattern (?year\d{4})(?month\d{2})(?day\d{2})_(?hh\d{2})(?mm\d{2})_(?site\w)\.sig; info regexp(file_list(1).name, pattern, names); % 把提取结果拼成时间向量 t datetime(str2double(info.year), str2double(info.month), str2double(info.day), ... str2double(info.hh), str2double(info.mm), 0);regexp配合命名分组(?name...)可以直接把年、月、日、时、分、站点名一次性抓出来。日期时间用datetime构造之后批量按时间排序或者按站点分组就非常方便。如果文件名格式不统一建议先用disp打印几个不匹配的样本再调整正则表达式。4.3 批处理里的异常隔离try-catch 保住整批数据批处理最怕的是第 50 个文件格式不对整个脚本崩溃前面 49 个结果白算。用try-catch把每个文件的处理包起来出错的样本记录到日志里其他样本继续跑。results table(); error_log {}; for i 1:numel(file_list) filename fullfile(data_dir, file_list(i).name); try % 读入 .sig 文件并检查 HR 段是否存在 T read_sig_file(filename); if ~isfield(T, HR) || isempty(T.HR) error(缺少 HR 光谱段); end % 重采样到 1nm 网格 target_wl (400:1:2400); refl interp1(T.HR(:,1), T.HR(:,2), target_wl, pchip); % SG 平滑 smooth_refl sgolayfilt(refl, 3, 15); % 当前样本的结果存入 table results [results; table( ... string(file_list(i).name), ... target_wl(:), smooth_refl(:), ... VariableNames, {FileName, Wavelength, Reflectance})]; catch ME % 记录错误信息继续下一个文件 error_log{end1, 1} file_list(i).name; error_log{end, 2} ME.message; fprintf(处理失败: %s, 原因: %s\n, file_list(i).name, ME.message); end end核心是try块里每一步都可能抛出异常catch ME把ME.message存下来这样整批任务跑完错误日志里一目了然。注意循环里不要用outerjoin或频繁拼接table来累积数据样本量上千时性能很差可以先存成 cell 数组最后一次性cell2table。4.4 批量导出保存 .mat 文件并生成汇总 CSV处理完所有.sig文件后通常需要两种产出一种是每个样本单独一个.mat文件方便后续按需加载另一种是一张宽表 CSV行是样本列是波长值是平滑反射率这是机器学习建模最常用的输入格式。% 把长表转成宽表每行一个样本每列一个波长 wideTable unstack(results, Reflectance, Wavelength); writetable(wideTable, fullfile(data_dir, all_smooth_spectra.csv)); % 另外保存原始长表作为详细记录 save(fullfile(data_dir, spectra_processing_results.mat), results, error_log);unstack是 MATLAB 里制作宽表的高效函数第一个参数是长表数据第二个参数是要展开的列第三个参数是分类依据的列。生成的宽表里列名就是波长值如果列数太多比如 2001 个波长CSV 文件会很大但 MATLAB 处理起来没问题。保存.mat文件时把error_log也存进去方便复盘哪些样本没处理成功。5. 批量结果质量检查用诊断表快速定位坏样本与过平滑样本处理完批量数据后别急着拿去训练模型先跑一遍质量检查。光看一眼平滑曲线是看不出问题的 光谱数量一大必须用数值诊断量把可疑样本筛出来。我一般会为每个样本计算下面这几个指标整理成诊断表然后用isoutlier批量标记异常。诊断字段计算方式异常判定标准WavelengthCount重采样后的有效波段数明显少于预期如少于 1800说明网格截断异常NegativeCount光滑后反射率小于 0 的个数大于 0 说明插值或平滑引入了非物理值NaNCount重采样后 NaN 的数量大于 0说明原始数据末端覆盖不足SmoothRMSEsqrt(mean((raw - smooth).^2))过大说明窗口参数不匹配或原始噪声异常SignalNoiseRatiomedian(smooth) / median(abs(raw - smooth))低于 10 说明光谱信噪比太差慎用% 为每个样本计算诊断量写入诊断表 diagTable table(); for i 1:height(wideTable) wl str2double(wideTable.Properties.VariableNames(2:end)); refl table2array(wideTable(i, 2:end)); NaNCount sum(isnan(refl)); NegativeCount sum(refl 0); SmoothRMSE sqrt(mean((refl - smooth(refl, 21, moving)).^2)); SignalNoiseRatio median(refl) / median(abs(refl - smooth(refl, 21, moving))); diagTable [diagTable; table( ... wideTable.FileName(i), NaNCount, NegativeCount, SmoothRMSE, SignalNoiseRatio, ... VariableNames, {FileName, NaNCount, NegativeCount, SmoothRMSE, SignalNoiseRatio})]; end % 标记异常样本 diagTable.Abnormal isoutlier(diagTable.SmoothRMSE) | ... diagTable.NegativeCount 0 | ... diagTable.SignalNoiseRatio 10;isoutlier是 R2017a 之后引入的函数默认用四分位距法判断离群值会自动算出上下边界。对于SmoothRMSE这种分布不规则的指标四分位距法比均值加减三倍标准差更稳健。计算SmoothRMSE时用smooth(refl, 21, moving)做对照是为了避免用和批处理流水线相同的 SG 参数从而独立评估平滑效果。批量质量检查最忌讳的是「看起来差不多就跳过」把诊断表存成 CSV 后用 Excel 透视表按站点或日期分组往往能发现某个仪器档位或者某一天采集的数据系统性偏差这些信息比调平滑参数值钱得多。本文还有配套的精品资源点击获取