简介：面向环境科学与水文学研究者，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(f=1..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-360 | 400-460 |
| 类富里酸 | 230-260, 305-350 | 370-450 |
| 类色氨酸 | 270-280 | 330-380 |
| 类酪氨酸 | 220-230, 270-280 | 300-330 |
4.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 结果一起打包存档，这是同行评审阶段最省事的做法。
本文还有配套的精品资源，点击获取
