尧图网站建设 
简介面向环境科学与水文学研究者DOMFluor.zip是一个MATLAB环境下的平行因子分析PARAFAC工具箱专用于溶解有机物DOM荧光光谱数据的成分解析与来源识别。压缩包共116个文件大小6.42MB内含100个m函数、11个mat数据文件以及csv、txt、xml等配置与文档覆盖数据预处理、模型拟合、结果评估和三维可视化等完整分析流程。已有459人浏览学习。借助内置示例数据和脚本用户可在MATLAB中直接开展背景扣除、噪声去除、因子数确定、模型质量验证与负荷图绘制等操作省去从零搭建算法的时间。对希望揭示DOM组成结构、追踪污染源或评估气候变化影响的研究者这份工具箱提供了完整分析链路和可扩展的代码框架只需具备一定MATLAB和统计学基础即可上手。1. 拿到 DOMFluor.zip 之后平行因子分析才算落地水环境样品的三维荧光光谱EEM堆起来是一个“样本 × 激发波长 × 发射波长”的三维数组常规二维分析只能对每个荧光峰做手动的区域积分十几个样品还可以上百个样品就会发现峰重叠、基线漂移、人工判读不一致的问题。平行因子分析PARAFAC这类三线性分解方法就是为了解决这个问题被引入的DOMFluor 则是 MATLAB 里做 EEM 平行因子分析最常用的工具箱。它把主成分的思想推广到三维把混合光谱拆成若干个稳定的荧光组分以及每个样本在每组组分上的相对浓度。适合人群是环境、水质、化学计量学方向的学生和工程师手头有批量 EEM 数据想在不用自己重写迭代算法的情况下得到可发表、可复现的组分模型。下面从工具箱装配、数据整理开始一直到组分数筛选与结果验证把整个流程按最省事的路径走一遍。2. 为什么平行因子分析能拆开 EEM三线性模型与 DOMFluor 装配2.1 三线性模型从“两张谱”到“三个载荷矩阵”EEM 数据的物理基础是荧光叠加原理在吸光度足够低的条件下每个纯组分的荧光强度近似等于该组分浓度与激发谱、发射谱三者的乘积。对单个样本来说这是一个二维矩阵对一批样本来说所有样本堆叠在一起就构成三维数组。主成分分析PCA只能在一个二维矩阵上做双线性分解所以面对混合样品时PCA 的载荷没有唯一方向必须做旋转才能给出物理解释而旋转往往又会破坏非负性。平行因子分析把三维数组直接写成三线性模型X(i,j,k) sum(f1..F) A(i,f) * B(j,f) * C(k,f) E(i,j,k)其中X(i,j,k)表示第 i 个样本、第 j 个激发波长、第 k 个发射波长处的荧光强度A是样本方向上的得分矩阵对应各组分的相对浓度B和C分别对应激发模式和发射模式上的载荷谱E是残差。只要组分间光谱形状有足够差异样本组成有足够变异这个三线性模型的解在大多数情况下是唯一的不需要旋转这也是 PARAFAC 被化学计量学界广泛接受的根本原因。求解过程通常用交替最小二乘ALS先固定 B、C用线性最小二乘更新 A再固定 A、C 更新 B再固定 A、B 更新 C如此循环直到残差平方和不再明显下降。DOMFluor 真正替你做的正是这部分迭代但它不是从零重写算法而是把数据规范化、模型初始化、约束设置、诊断输出包装成更适合 EEM 分析的工作流底层求解仍然依赖 N-way Toolbox 的parafac函数。这一点直接影响你后面怎么配环境只解压 DOMFluor 而不装 N-way模型根本跑不起来。2.2 在 MATLAB 里配齐 DOMFluor 与 N-way 依赖网上能下载到的 DOMFluor 通常是一个 zip 压缩包解压后里面有若干.m文件也可能带示例数据N-way Toolbox 是另一个独立压缩包。安装顺序上没有硬性要求但路径一定要确认能访问到。常见做法是把两个目录放在同一个工具目录下然后用addpath挂载而不是把文件直接拷进 MATLAB 安装目录后者在 MATLAB 升级时容易被覆盖。% 解压后把两个工具箱目录都加进 MATLAB 搜索路径 addpath(genpath(D:\tools\DOMFluor)); addpath(genpath(D:\tools\nway30)); % 保存路径避免下次启动 MATLAB 重新配置 status savepath; if status 0 disp(路径保存成功); else disp(路径保存失败请检查 MATLAB 工作目录写入权限); endaddpath(genpath(...))里的genpath会把指定目录下所有子目录一并加入搜索路径这对工具目录结构不透明的压缩包最稳妥少写一层子目录也不会漏。savepath是把你当前会话里的路径设置写到pathdef.m这样下次启动就不用再执行一遍addpath。如果提示保存失败多半是 MATLAB 安装目录只读可以用userpath下的自定义目录或者每次都先运行一个startup.m。挂载完成后的第一件事是验证依赖完整直接查两个关键函数是否在路径上which parafac which corcond which domfluorwhich命令会返回函数的完整路径。如果parafac返回“未找到”说明 N-way 没挂载成功如果corcond找不到说明 N-way 版本里没有核一致性诊断函数后面做模型阶数筛选时会卡住。domfluor能找到则说明 DOMFluor 本体可执行。再用help parafac看一眼函数签名能正常弹出帮助文本就代表版本基本兼容。提示路径中尽量不要出现中文和$、这类特殊字符。旧版 MATLAB 对 Unicode 路径支持不佳工具箱文件一旦放到中文目录下容易在启动时出现“找不到文件或函数”的报错排查起来非常费时间。3. 把 CSV 和原始 EEM 整理成平行因子分析的输入数组3.1 三维数组排布样本 × 激发 × 发射的顺序规则平行因子分析的输入必须是一个三维数组三个维度的顺序没有数学上的硬性要求但建议统一固定为“样本 × 激发 × 发射”。这样在所有后续操作里你只要一看到size(X)是[nSample, nEx, nEm]就能立刻知道每个维度代表什么避免在掩膜、画图、解释载荷时把 Ex 和 Em 轴搞反。仪器导出的 EEM 通常是一个“激发波长 × 发射波长”的二维矩阵每个样本一个文件。堆叠前先要把所有样本插值到统一的波长网格上因为不同批次样品的扫描范围或步长可能不一致直接堆叠会形成大量NaN缺口PARAFAC 迭代时这些缺口会拖慢收敛甚至导致模型不稳定。EmCommon 250:2:600; % 统一的发射波长轴 X nan(nSample, length(Ex), length(EmCommon)); for i 1:nSample % F_raw 是当前样本的原始 EEM行激发列原始发射波长 F_interp interp1(Em_raw, F_raw, EmCommon, pchip); X(i, :, :) F_interp; end size(X) % 输出: nSample nEx nEminterp1默认沿着数组第一个非单一维度插值所以原始矩阵要先转置让发射波长变成行方向插值完成后转置回来这样行列语义不会乱。pchip是分段三次 Hermite 插值荧光光谱通常是平滑峰形用pchip比线性插值更自然也不会像样条插值那样出现过度振荡。堆叠之前还要检查波长坐标是否单调递增。个别仪器导出的 Excel 表格会把激发波长按降序排列这种情况下直接interp1会报错。先执行issorted(Ex)和issorted(Em_raw)返回 0 就手动flipud或fliplr调整。3.2 去除瑞利散射与拉曼散射再喂给模型三维荧光谱上最刺眼的“假峰”是瑞利散射一阶瑞利满足发射波长约等于激发波长二阶瑞利满足发射波长约等于两倍激发波长。散射带强度往往比荧光信号高一个数量级如果不处理PARAFAC 会把散射区域当成一个“组分”提取出来不仅掩盖真实组分峰还会让核一致性诊断崩溃。常见做法是先对散射区域做掩膜把异常值置为NaN再决定是否插值填补。建模前是否要填补取决于你用哪个求解函数部分 PARAFAC 实现能直接跳过NaN但更多实现要求完整数组。保险起见我们同时准备两套版本。[EmG, ExG] meshgrid(Em, Ex); M ones(nEx, nEm); M(EmG - ExG 30) NaN; % 一阶瑞利附近 M(EmG - 2*ExG 25) NaN; % 二阶瑞利附近 % 乘上掩膜散射区域全部变成 NaN X_masked X; for i 1:nSample X_masked(i, :, :) squeeze(X(i, :, :)) .* M; end这里的30和25是经验带宽对应仪器狭缝宽度和波长偏移量。实际使用中可以先拿一个纯水空白样看看散射带实际宽度再调整阈值不要照抄文献数值。掩膜裁出来的凹陷区域如果要补洞就在发射方向做线性插值X_fill X_masked; for i 1:nSample for j 1:nEx row squeeze(X_fill(i, j, :)); row fillmissing(row, linear, EndValues, nearest); X_fill(i, j, :) row; end endfillmissing在 MATLAB R2016b 及以后版本可用linear表示沿发射波长方向线性插值EndValues, nearest表示散射带两端用最近邻值补齐避免数组两端出现NaN。这个补洞策略对散射带这种窄带缺失够用不建议用高阶多项式否则会把散射残余带进真实荧光区域。3.3 从 CSV 批量导入 EEM 数据的脚本模板很多荧光仪器支持把数据导出成 CSV常见有两种格式长表格式每行是一个测量点列包括激发波长、发射波长、荧光强度、样本编号宽表格式则是每行一个激发波长每列一个发射波长。这里给一套兼容长表格式的批量导入模板。files dir(eem/*.csv); for k 1:length(files) T readmatrix(fullfile(files(k).folder, files(k).name)); ExV T(:, 1); % 表头: Excitation EmV T(:, 2); % 表头: Emission val T(:, 3); % 表头: Value if k 1 % 用第一个文件构造三维数组 Ex unique(ExV); Em unique(EmV); X nan(length(files), length(Ex), length(Em)); end % 通过索引映射到三维数组 [~, ix] ismember(ExV, Ex); [~, iy] ismember(EmV, Em); sub sub2ind([length(files), length(Ex), length(Em)], ... k*ones(size(ix)), ix, iy); X(sub) val; endreadmatrix是 R2018b 以后推荐的读取函数能自动识别数值列旧版可以用csvread代替但对列顺序敏感。unique取出波长坐标后ismember把每个测量点的波长映射到网格下标sub2ind一次性把长表数据灌进三维数组避免了双层循环。整套脚本的代价是要求所有文件使用完全一致的波长表所以前面的统一插值步骤不能省。CSV 列名含义处理去向Excitation激发波长nm三维数组第 2 维坐标Emission发射波长nm三维数组第 3 维坐标Value荧光强度三维数组元素SampleID样本编号第 1 维索引4. 用 DOMFluor 跑平行因子分析组分数选择与结果解读4.1 图形界面与脚本的取舍在 MATLAB 命令窗口输入domfluor如果路径配置正确会弹出 DOMFluor 主界面。界面里通常提供数据导入、组分数量设置、约束选择、运行与诊断等入口适合第一次上手时直观感受“跑一个模型需要设置哪些东西”。但图形界面不适合批量处理比如你要对比 2 到 6 个组分模型的核一致性在界面上挨个点会非常低效。我的习惯是用 GUI 读一遍示例数据确认工具箱确实能跑通正式分析全部走脚本。这样组分数扫描、残差计算、图表导出都能一键复现也方便以后换一批数据时直接改路径重跑。脚本的核心参数是组分数量一般先扫描 2 到 6 个组分。组分数量太少不同荧光团会被强行合并数量太多模型会把噪声和散射残余拟合进来出现负载荷或者峰位漂移。4.2 用核一致性诊断确定组分数量核一致性core consistency是目前最常用的阶数判别指标原理是把 PARAFAC 解出的载荷反推回一个超对角核再比较这个核与理想单位核的接近程度。接近 100 表示模型和数据的三线性结构吻合明显偏低或为负说明数据被过度分解。在 N-way Toolbox 里对应的函数是corcond。rng(2024); % 固定随机种子保证试验可复现 X X_fill; % 上一步处理好的三维数组 for n 2:6 [A, B, C] parafac(X, n); % 跑 n 组分模型 corco corcond(X, A, B, C); % 核一致性诊断 fprintf(组分数量 %d - 核一致性 %.1f%%\n, n, corco); endparafac的前两个参数是数据数组和组分数量 F返回值 A、B、C 分别是三个模式上的载荷矩阵。迭代结果对初始值敏感而初始值往往是随机生成的所以循环开始前先用rng(2024)固定随机数生成器否则同一次扫描每次跑出来的核一致性都会略有波动无法判断是数据问题还是随机性导致。parafac后面还可以接约束参数和选项结构体但默认无约束模型能跑通之前不建议一上来就加复杂约束。核一致性只适合横向比较不同 F 值的相对好坏不能单独作为“最优模型”的最终裁决。我一般遵循下面这套判断逻辑核一致性值模型状态后续动作大于 90%结构非常符合三线性可以继续做 split-half 验证50% ~ 90%可接受但需谨慎对比相邻组分数的载荷图接近 0 或负数阶数过高或数据本身非三线性降低组分数量或检查散射处理4.3 载荷图与荧光峰位鉴定组分数量初步确定后把 B、C 载荷分别画到发射和激发波长轴上就能看到每个组分的光谱形状。荧光峰的峰位是鉴定组分的核心依据。比如类腐殖质 C 峰通常出现在激发 320-360 nm、发射 400-460 nm 的位置类色氨酸则更靠蓝移。下面这段代码生成标准的两段式载荷图figure; subplot(2, 1, 1); plot(Ex, A(:, 1), LineWidth, 1.5); hold on; plot(Ex, A(:, 2), LineWidth, 1.5); plot(Ex, A(:, 3), LineWidth, 1.5); hold off; xlabel(激发波长 (nm)); ylabel(激发载荷); legend(组分1, 组分2, 组分3); subplot(2, 1, 2); plot(Em, B(:, 1), LineWidth, 1.5); hold on; plot(Em, B(:, 2), LineWidth, 1.5); plot(Em, B(:, 3), LineWidth, 1.5); hold off; xlabel(发射波长 (nm)); ylabel(发射载荷);画完图先看每个载荷是否光滑、是否只有一个明显峰。如果某个组分的发射载荷出现多个尖锐毛刺或者像镜像一样成对出现基本可以判定模型阶数不对。下图是文献中常见荧光峰位的经验区间不同仪器会有 5-10 nm 的偏差做比对时允许一定容差。组分类型激发峰nm发射峰nm类腐殖质 C 峰320-360400-460类富里酸230-260, 305-350370-450类色氨酸270-280330-380类酪氨酸220-230, 270-280300-3304.4 用 split-half 验证模型稳定性载荷图看着合理还不够还要验证模型的稳定性。最通用的验证方法是 split-half把样本随机分成两半分别用相同组分数量建模比较两个模型在相同模式上的载荷谱是否高度一致。理想情况下两半数据应该还原出几乎相同的 B、C 载荷。rng(7); perm randperm(size(X, 1)); half floor(length(perm) / 2); Xa X(perm(1:half), :, :); Xb X(perm(half 1:end), :, :); [Aa, Ba, Ca] parafac(Xa, 3); [Ab, Bb, Cb] parafac(Xb, 3); for f 1:3 r corr(Ba(:, f), Bb(:, f)); fprintf(组分 %d 发射载荷相关系数: %.3f\n, f, r); end这里用相关系数做快速对照学术上更严谨的做法是计算 Tucker congruence阈值通常在 0.95 以上。需要注意 split-half 对样本量有要求每个子集至少要保留几十个样本否则子模型本身就不稳定验证意义不大。样本量不足时可以用残差平方和的肘部图辅助判断模型阶数。5. DOMFluor 平行因子分析的验证技巧与常见报错排查5.1 用随机种子锁住可复现结果PARAFAC 的 ALS 迭代用到随机初始化每次运行结果会有细微差异。正式报告里的最终模型一定要固定随机种子后再运行并在方法部分写明随机种子和迭代终止条件。推荐在建模脚本最顶部写一行rng(42)不要把它埋在循环内部。如果复现时发现两次运行结果差异很大优先怀疑数据里存在大量NaN或散射残余而不是随机性问题。最终结果用save命令存成.mat文件save(parafac_result.mat, A, B, C, X_fill, rng_state, -v7.3);.mat文件体积超过 2 GB 时需要-v7.3格式普通小数据不加也没问题。rng_state可以通过rng函数无参调用取得存下来是为了以后精确复现建模环境。如果要在 MATLAB 之外继续分析可以用 Python 的scipy.io读取这个.mat文件A、B、C 三个矩阵会原样导出成 NumPy 数组这一步在换工具链做作图或统计时很常用。5.2 轴顺序与波长区间不一致的排查跑完模型发现激发载荷峰位整体偏移先别急着怀疑模型。检查插值时用的Ex、Em顺序以及绘制载荷图时用的是哪一个矩阵。常见错误是把发射模式载荷当成激发模式画了或者两个波长轴的顺序一个是递增一个是递减导致图像左右翻转。确认方法很简单取一个已知组分的载荷峰位对照 4.3 的表格峰位偏差超过 15 nm 就应该回到数据整理阶段重新检查。散射线掩膜也可能把真实荧光切掉。掩膜带宽设得过大时蓝移组分如类酪氨酸的短波长侧会被吃掉载荷图出现截断峰。遇到这种情况调小带宽后重新运行模型对比峰位是否移动。5.3 生成可直接插入论文的矢量图画载荷谱和 EEM 图时先用set(gcf, Color, w)把图形背景设为白色再用exportgraphics导出高分辨率图片这样线条不会发虚set(gcf, Color, w); exportgraphics(gcf, loading_3comp.eps, Resolution, 300);exportgraphics支持-v7.3无关的eps、png、pdf等格式Resolution参数控制导出分辨率。如果画 EEM 等值线图时波长点太多横轴刻度会挤成一条黑带用xticks手动指定显示位置例如xticks(250:50:600)同时配合xlim把有效波长区间卡紧。最后在投稿前把最终模型脚本、随机种子、核一致性数值和 split-half 结果一起打包存档这是同行评审阶段最省事的做法。本文还有配套的精品资源点击获取 MATLAB中DOMFluor工具箱实现EEM平行因子分析全流程指南
简介面向环境科学与水文学研究者DOMFluor.zip是一个MATLAB环境下的平行因子分析PARAFAC工具箱专用于溶解有机物DOM荧光光谱数据的成分解析与来源识别。压缩包共116个文件大小6.42MB内含100个m函数、11个mat数据文件以及csv、txt、xml等配置与文档覆盖数据预处理、模型拟合、结果评估和三维可视化等完整分析流程。已有459人浏览学习。借助内置示例数据和脚本用户可在MATLAB中直接开展背景扣除、噪声去除、因子数确定、模型质量验证与负荷图绘制等操作省去从零搭建算法的时间。对希望揭示DOM组成结构、追踪污染源或评估气候变化影响的研究者这份工具箱提供了完整分析链路和可扩展的代码框架只需具备一定MATLAB和统计学基础即可上手。1. 拿到 DOMFluor.zip 之后平行因子分析才算落地水环境样品的三维荧光光谱EEM堆起来是一个“样本 × 激发波长 × 发射波长”的三维数组常规二维分析只能对每个荧光峰做手动的区域积分十几个样品还可以上百个样品就会发现峰重叠、基线漂移、人工判读不一致的问题。平行因子分析PARAFAC这类三线性分解方法就是为了解决这个问题被引入的DOMFluor 则是 MATLAB 里做 EEM 平行因子分析最常用的工具箱。它把主成分的思想推广到三维把混合光谱拆成若干个稳定的荧光组分以及每个样本在每组组分上的相对浓度。适合人群是环境、水质、化学计量学方向的学生和工程师手头有批量 EEM 数据想在不用自己重写迭代算法的情况下得到可发表、可复现的组分模型。下面从工具箱装配、数据整理开始一直到组分数筛选与结果验证把整个流程按最省事的路径走一遍。2. 为什么平行因子分析能拆开 EEM三线性模型与 DOMFluor 装配2.1 三线性模型从“两张谱”到“三个载荷矩阵”EEM 数据的物理基础是荧光叠加原理在吸光度足够低的条件下每个纯组分的荧光强度近似等于该组分浓度与激发谱、发射谱三者的乘积。对单个样本来说这是一个二维矩阵对一批样本来说所有样本堆叠在一起就构成三维数组。主成分分析PCA只能在一个二维矩阵上做双线性分解所以面对混合样品时PCA 的载荷没有唯一方向必须做旋转才能给出物理解释而旋转往往又会破坏非负性。平行因子分析把三维数组直接写成三线性模型X(i,j,k) sum(f1..F) A(i,f) * B(j,f) * C(k,f) E(i,j,k)其中X(i,j,k)表示第 i 个样本、第 j 个激发波长、第 k 个发射波长处的荧光强度A是样本方向上的得分矩阵对应各组分的相对浓度B和C分别对应激发模式和发射模式上的载荷谱E是残差。只要组分间光谱形状有足够差异样本组成有足够变异这个三线性模型的解在大多数情况下是唯一的不需要旋转这也是 PARAFAC 被化学计量学界广泛接受的根本原因。求解过程通常用交替最小二乘ALS先固定 B、C用线性最小二乘更新 A再固定 A、C 更新 B再固定 A、B 更新 C如此循环直到残差平方和不再明显下降。DOMFluor 真正替你做的正是这部分迭代但它不是从零重写算法而是把数据规范化、模型初始化、约束设置、诊断输出包装成更适合 EEM 分析的工作流底层求解仍然依赖 N-way Toolbox 的parafac函数。这一点直接影响你后面怎么配环境只解压 DOMFluor 而不装 N-way模型根本跑不起来。2.2 在 MATLAB 里配齐 DOMFluor 与 N-way 依赖网上能下载到的 DOMFluor 通常是一个 zip 压缩包解压后里面有若干.m文件也可能带示例数据N-way Toolbox 是另一个独立压缩包。安装顺序上没有硬性要求但路径一定要确认能访问到。常见做法是把两个目录放在同一个工具目录下然后用addpath挂载而不是把文件直接拷进 MATLAB 安装目录后者在 MATLAB 升级时容易被覆盖。% 解压后把两个工具箱目录都加进 MATLAB 搜索路径 addpath(genpath(D:\tools\DOMFluor)); addpath(genpath(D:\tools\nway30)); % 保存路径避免下次启动 MATLAB 重新配置 status savepath; if status 0 disp(路径保存成功); else disp(路径保存失败请检查 MATLAB 工作目录写入权限); endaddpath(genpath(...))里的genpath会把指定目录下所有子目录一并加入搜索路径这对工具目录结构不透明的压缩包最稳妥少写一层子目录也不会漏。savepath是把你当前会话里的路径设置写到pathdef.m这样下次启动就不用再执行一遍addpath。如果提示保存失败多半是 MATLAB 安装目录只读可以用userpath下的自定义目录或者每次都先运行一个startup.m。挂载完成后的第一件事是验证依赖完整直接查两个关键函数是否在路径上which parafac which corcond which domfluorwhich命令会返回函数的完整路径。如果parafac返回“未找到”说明 N-way 没挂载成功如果corcond找不到说明 N-way 版本里没有核一致性诊断函数后面做模型阶数筛选时会卡住。domfluor能找到则说明 DOMFluor 本体可执行。再用help parafac看一眼函数签名能正常弹出帮助文本就代表版本基本兼容。提示路径中尽量不要出现中文和$、这类特殊字符。旧版 MATLAB 对 Unicode 路径支持不佳工具箱文件一旦放到中文目录下容易在启动时出现“找不到文件或函数”的报错排查起来非常费时间。3. 把 CSV 和原始 EEM 整理成平行因子分析的输入数组3.1 三维数组排布样本 × 激发 × 发射的顺序规则平行因子分析的输入必须是一个三维数组三个维度的顺序没有数学上的硬性要求但建议统一固定为“样本 × 激发 × 发射”。这样在所有后续操作里你只要一看到size(X)是[nSample, nEx, nEm]就能立刻知道每个维度代表什么避免在掩膜、画图、解释载荷时把 Ex 和 Em 轴搞反。仪器导出的 EEM 通常是一个“激发波长 × 发射波长”的二维矩阵每个样本一个文件。堆叠前先要把所有样本插值到统一的波长网格上因为不同批次样品的扫描范围或步长可能不一致直接堆叠会形成大量NaN缺口PARAFAC 迭代时这些缺口会拖慢收敛甚至导致模型不稳定。EmCommon 250:2:600; % 统一的发射波长轴 X nan(nSample, length(Ex), length(EmCommon)); for i 1:nSample % F_raw 是当前样本的原始 EEM行激发列原始发射波长 F_interp interp1(Em_raw, F_raw, EmCommon, pchip); X(i, :, :) F_interp; end size(X) % 输出: nSample nEx nEminterp1默认沿着数组第一个非单一维度插值所以原始矩阵要先转置让发射波长变成行方向插值完成后转置回来这样行列语义不会乱。pchip是分段三次 Hermite 插值荧光光谱通常是平滑峰形用pchip比线性插值更自然也不会像样条插值那样出现过度振荡。堆叠之前还要检查波长坐标是否单调递增。个别仪器导出的 Excel 表格会把激发波长按降序排列这种情况下直接interp1会报错。先执行issorted(Ex)和issorted(Em_raw)返回 0 就手动flipud或fliplr调整。3.2 去除瑞利散射与拉曼散射再喂给模型三维荧光谱上最刺眼的“假峰”是瑞利散射一阶瑞利满足发射波长约等于激发波长二阶瑞利满足发射波长约等于两倍激发波长。散射带强度往往比荧光信号高一个数量级如果不处理PARAFAC 会把散射区域当成一个“组分”提取出来不仅掩盖真实组分峰还会让核一致性诊断崩溃。常见做法是先对散射区域做掩膜把异常值置为NaN再决定是否插值填补。建模前是否要填补取决于你用哪个求解函数部分 PARAFAC 实现能直接跳过NaN但更多实现要求完整数组。保险起见我们同时准备两套版本。[EmG, ExG] meshgrid(Em, Ex); M ones(nEx, nEm); M(EmG - ExG 30) NaN; % 一阶瑞利附近 M(EmG - 2*ExG 25) NaN; % 二阶瑞利附近 % 乘上掩膜散射区域全部变成 NaN X_masked X; for i 1:nSample X_masked(i, :, :) squeeze(X(i, :, :)) .* M; end这里的30和25是经验带宽对应仪器狭缝宽度和波长偏移量。实际使用中可以先拿一个纯水空白样看看散射带实际宽度再调整阈值不要照抄文献数值。掩膜裁出来的凹陷区域如果要补洞就在发射方向做线性插值X_fill X_masked; for i 1:nSample for j 1:nEx row squeeze(X_fill(i, j, :)); row fillmissing(row, linear, EndValues, nearest); X_fill(i, j, :) row; end endfillmissing在 MATLAB R2016b 及以后版本可用linear表示沿发射波长方向线性插值EndValues, nearest表示散射带两端用最近邻值补齐避免数组两端出现NaN。这个补洞策略对散射带这种窄带缺失够用不建议用高阶多项式否则会把散射残余带进真实荧光区域。3.3 从 CSV 批量导入 EEM 数据的脚本模板很多荧光仪器支持把数据导出成 CSV常见有两种格式长表格式每行是一个测量点列包括激发波长、发射波长、荧光强度、样本编号宽表格式则是每行一个激发波长每列一个发射波长。这里给一套兼容长表格式的批量导入模板。files dir(eem/*.csv); for k 1:length(files) T readmatrix(fullfile(files(k).folder, files(k).name)); ExV T(:, 1); % 表头: Excitation EmV T(:, 2); % 表头: Emission val T(:, 3); % 表头: Value if k 1 % 用第一个文件构造三维数组 Ex unique(ExV); Em unique(EmV); X nan(length(files), length(Ex), length(Em)); end % 通过索引映射到三维数组 [~, ix] ismember(ExV, Ex); [~, iy] ismember(EmV, Em); sub sub2ind([length(files), length(Ex), length(Em)], ... k*ones(size(ix)), ix, iy); X(sub) val; endreadmatrix是 R2018b 以后推荐的读取函数能自动识别数值列旧版可以用csvread代替但对列顺序敏感。unique取出波长坐标后ismember把每个测量点的波长映射到网格下标sub2ind一次性把长表数据灌进三维数组避免了双层循环。整套脚本的代价是要求所有文件使用完全一致的波长表所以前面的统一插值步骤不能省。CSV 列名含义处理去向Excitation激发波长nm三维数组第 2 维坐标Emission发射波长nm三维数组第 3 维坐标Value荧光强度三维数组元素SampleID样本编号第 1 维索引4. 用 DOMFluor 跑平行因子分析组分数选择与结果解读4.1 图形界面与脚本的取舍在 MATLAB 命令窗口输入domfluor如果路径配置正确会弹出 DOMFluor 主界面。界面里通常提供数据导入、组分数量设置、约束选择、运行与诊断等入口适合第一次上手时直观感受“跑一个模型需要设置哪些东西”。但图形界面不适合批量处理比如你要对比 2 到 6 个组分模型的核一致性在界面上挨个点会非常低效。我的习惯是用 GUI 读一遍示例数据确认工具箱确实能跑通正式分析全部走脚本。这样组分数扫描、残差计算、图表导出都能一键复现也方便以后换一批数据时直接改路径重跑。脚本的核心参数是组分数量一般先扫描 2 到 6 个组分。组分数量太少不同荧光团会被强行合并数量太多模型会把噪声和散射残余拟合进来出现负载荷或者峰位漂移。4.2 用核一致性诊断确定组分数量核一致性core consistency是目前最常用的阶数判别指标原理是把 PARAFAC 解出的载荷反推回一个超对角核再比较这个核与理想单位核的接近程度。接近 100 表示模型和数据的三线性结构吻合明显偏低或为负说明数据被过度分解。在 N-way Toolbox 里对应的函数是corcond。rng(2024); % 固定随机种子保证试验可复现 X X_fill; % 上一步处理好的三维数组 for n 2:6 [A, B, C] parafac(X, n); % 跑 n 组分模型 corco corcond(X, A, B, C); % 核一致性诊断 fprintf(组分数量 %d - 核一致性 %.1f%%\n, n, corco); endparafac的前两个参数是数据数组和组分数量 F返回值 A、B、C 分别是三个模式上的载荷矩阵。迭代结果对初始值敏感而初始值往往是随机生成的所以循环开始前先用rng(2024)固定随机数生成器否则同一次扫描每次跑出来的核一致性都会略有波动无法判断是数据问题还是随机性导致。parafac后面还可以接约束参数和选项结构体但默认无约束模型能跑通之前不建议一上来就加复杂约束。核一致性只适合横向比较不同 F 值的相对好坏不能单独作为“最优模型”的最终裁决。我一般遵循下面这套判断逻辑核一致性值模型状态后续动作大于 90%结构非常符合三线性可以继续做 split-half 验证50% ~ 90%可接受但需谨慎对比相邻组分数的载荷图接近 0 或负数阶数过高或数据本身非三线性降低组分数量或检查散射处理4.3 载荷图与荧光峰位鉴定组分数量初步确定后把 B、C 载荷分别画到发射和激发波长轴上就能看到每个组分的光谱形状。荧光峰的峰位是鉴定组分的核心依据。比如类腐殖质 C 峰通常出现在激发 320-360 nm、发射 400-460 nm 的位置类色氨酸则更靠蓝移。下面这段代码生成标准的两段式载荷图figure; subplot(2, 1, 1); plot(Ex, A(:, 1), LineWidth, 1.5); hold on; plot(Ex, A(:, 2), LineWidth, 1.5); plot(Ex, A(:, 3), LineWidth, 1.5); hold off; xlabel(激发波长 (nm)); ylabel(激发载荷); legend(组分1, 组分2, 组分3); subplot(2, 1, 2); plot(Em, B(:, 1), LineWidth, 1.5); hold on; plot(Em, B(:, 2), LineWidth, 1.5); plot(Em, B(:, 3), LineWidth, 1.5); hold off; xlabel(发射波长 (nm)); ylabel(发射载荷);画完图先看每个载荷是否光滑、是否只有一个明显峰。如果某个组分的发射载荷出现多个尖锐毛刺或者像镜像一样成对出现基本可以判定模型阶数不对。下图是文献中常见荧光峰位的经验区间不同仪器会有 5-10 nm 的偏差做比对时允许一定容差。组分类型激发峰nm发射峰nm类腐殖质 C 峰320-360400-460类富里酸230-260, 305-350370-450类色氨酸270-280330-380类酪氨酸220-230, 270-280300-3304.4 用 split-half 验证模型稳定性载荷图看着合理还不够还要验证模型的稳定性。最通用的验证方法是 split-half把样本随机分成两半分别用相同组分数量建模比较两个模型在相同模式上的载荷谱是否高度一致。理想情况下两半数据应该还原出几乎相同的 B、C 载荷。rng(7); perm randperm(size(X, 1)); half floor(length(perm) / 2); Xa X(perm(1:half), :, :); Xb X(perm(half 1:end), :, :); [Aa, Ba, Ca] parafac(Xa, 3); [Ab, Bb, Cb] parafac(Xb, 3); for f 1:3 r corr(Ba(:, f), Bb(:, f)); fprintf(组分 %d 发射载荷相关系数: %.3f\n, f, r); end这里用相关系数做快速对照学术上更严谨的做法是计算 Tucker congruence阈值通常在 0.95 以上。需要注意 split-half 对样本量有要求每个子集至少要保留几十个样本否则子模型本身就不稳定验证意义不大。样本量不足时可以用残差平方和的肘部图辅助判断模型阶数。5. DOMFluor 平行因子分析的验证技巧与常见报错排查5.1 用随机种子锁住可复现结果PARAFAC 的 ALS 迭代用到随机初始化每次运行结果会有细微差异。正式报告里的最终模型一定要固定随机种子后再运行并在方法部分写明随机种子和迭代终止条件。推荐在建模脚本最顶部写一行rng(42)不要把它埋在循环内部。如果复现时发现两次运行结果差异很大优先怀疑数据里存在大量NaN或散射残余而不是随机性问题。最终结果用save命令存成.mat文件save(parafac_result.mat, A, B, C, X_fill, rng_state, -v7.3);.mat文件体积超过 2 GB 时需要-v7.3格式普通小数据不加也没问题。rng_state可以通过rng函数无参调用取得存下来是为了以后精确复现建模环境。如果要在 MATLAB 之外继续分析可以用 Python 的scipy.io读取这个.mat文件A、B、C 三个矩阵会原样导出成 NumPy 数组这一步在换工具链做作图或统计时很常用。5.2 轴顺序与波长区间不一致的排查跑完模型发现激发载荷峰位整体偏移先别急着怀疑模型。检查插值时用的Ex、Em顺序以及绘制载荷图时用的是哪一个矩阵。常见错误是把发射模式载荷当成激发模式画了或者两个波长轴的顺序一个是递增一个是递减导致图像左右翻转。确认方法很简单取一个已知组分的载荷峰位对照 4.3 的表格峰位偏差超过 15 nm 就应该回到数据整理阶段重新检查。散射线掩膜也可能把真实荧光切掉。掩膜带宽设得过大时蓝移组分如类酪氨酸的短波长侧会被吃掉载荷图出现截断峰。遇到这种情况调小带宽后重新运行模型对比峰位是否移动。5.3 生成可直接插入论文的矢量图画载荷谱和 EEM 图时先用set(gcf, Color, w)把图形背景设为白色再用exportgraphics导出高分辨率图片这样线条不会发虚set(gcf, Color, w); exportgraphics(gcf, loading_3comp.eps, Resolution, 300);exportgraphics支持-v7.3无关的eps、png、pdf等格式Resolution参数控制导出分辨率。如果画 EEM 等值线图时波长点太多横轴刻度会挤成一条黑带用xticks手动指定显示位置例如xticks(250:50:600)同时配合xlim把有效波长区间卡紧。最后在投稿前把最终模型脚本、随机种子、核一致性数值和 split-half 结果一起打包存档这是同行评审阶段最省事的做法。本文还有配套的精品资源点击获取