简介Ecopath 食物网质量平衡算法的 MATLAB 实现面向生态学、渔业资源管理与海洋或淡水生态系统研究者用于构建食物网模型定量分析物种生物量、摄食关系与能量流动在生态建模与资源评估中具有实用价值。压缩包共一百零四个文件大小约一点零六兆字节核心为四十四个算法源文件辅助材料包含二十四个数据表、若干脚本与数据库文件以及可视化图像和说明文档覆盖数据预处理、能量平衡迭代计算、食物网构建、营养级统计与结果输出的完整流程。另有测试脚本和示例数据可供复现与二次开发并有助于开展后续研究。目前已有三百二十三人学习浏览代码组织清晰、目录划分明确适合希望快速上手生态模型、开展生态系统能量流动与稳定性分析的 MATLAB 使用者。1. 当食物网模型遇上矩阵求逆Ecopath 质量平衡算法的 MATLAB 落地生态学里有个挺反直觉的事实一个看上去物种丰富的生态系统往往只被少数几条能量通路支配。Ecopath 模型要做的就是把这个“谁吃谁、吃多少、剩多少”的复杂关系压缩成一组线性方程用质量平衡去反推那些无法直接观测的摄食率。早年做这类分析基本靠 Ecopath with Ecosim 的桌面软件但当你需要批量跑场景、把模型嵌入蒙特卡洛模拟、或者动态调整功能组划分时GUI 就成瓶颈了。这个 MATLAB 实现的价值在于它把质量平衡的核心算法拆成了可读、可改的函数包配合 REco_Flatfish1.csv、tb_diet.csv 这类数据文件你可以在脚本里完成从数据加载、平衡求解到网络指标计算的全流程。适合有三类需求的读者想搞懂质量平衡算法内部迭代逻辑的需要把 Ecopath 接入自定义工作流的以及被桌面版 License 限制但手头有 MATLAB 环境的研究者。2. 从 CSV 到系数矩阵Ecopath 方程组的数据建模与 MATLAB 数据结构Ecopath 的数学模型核心是 Master Equation每个功能组 i 满足Bi * (P/B)i * EEi - Σj [Bj * (Q/B)j * DCji] - EXi 0其中 Bi 是生物量P/B 是生产量/生物量比EE 是生态效率ecotrophic efficiencyQ/B 是消费量/生物量比DCji 表示 j 组对 i 组的摄食比例EXi 是输出捕捞加迁移。这里关键在于当模型假设系统处于平衡态时方程组中每个功能组都有一个未知数通常是某几个组的生物量或 EE并且方程组本身是线性的。这也解释了为什么 Ecopath 被称作“质量平衡”模型——它不模拟时间动态它求解的是稳态下的流量分配。2.1 数据文件的读取策略与字段映射拿到资源包后第一件事是检查 CSV 的结构。以 REco_Flatfish1.csv 和 REco_Roundfish1.csv 为例这类文件通常以功能组为行单位列包含 biomass、pb即 P/B 值、qbQ/B 值、diet 矩阵等。MATLAB 中处理 CSV 数据最稳妥的方式是用 readtable因为它保留列名且能识别混合类型。下面是我在加载食物网数据时常用的模式% 加载模型参数表每个功能组一行 T readtable(tb_model.csv, PreserveVariableNames, true); % 加载摄食矩阵行捕食者列猎物值为DCji D readtable(tb_diet.csv); dietMat D{:, 2:end}; % 第一列是功能组名称其余列是摄食比例 bio T.biomass; pb T.pb; % P/B 生产量与生物量之比 qb T.qb; % Q/B 消费量与生物量之比 % 捕获量渔业输出若没有则置0 land T.landing;readtable 返回的表对象在处理生态学数据时有一个隐性坑它会自动将列名中的括号或数字前缀转为合法字段名。我在处理 REco_Roundfish2.csv 时就遇到过列名被改写成x_PB的情况。所以 PreserveVariableNames 参数建议始终打开字段访问用T.(P/B)而不是T.P_B。加载完数据后需要验证摄食矩阵的行和与 qb 列的关系sum(dietMat, 2)应该接近 1否则说明食性数据是原始观测值还没有归一化。2.2 生物量缺失组的处理连续方程补充在真实的渔业数据里很难做到所有功能组都有独立的生物量调查。底栖鱼类的生物量可以从底拖网估算但浮游动物呢这时候需要用连续方程来反解。Ecopath 的处理方式是为缺失生物量的组预设 EE 值通常取 0.95然后把 Bi 作为未知数解出来。在 MATLAB 里实现时需要把未知数分成两组已知 Bi 求 EE和已知 EE 求 Bi。这也是质量平衡矩阵组装时最容易出错的地方——你必须把方程组重新组织成矩阵形式 A * x b。其中 x 是未知量向量A 的每一行对应一个功能组的平衡方程。对于未知 Bi 的功能组将其从系数矩阵中转移到 x 侧对于未知 EE则把 Bi 乘进系数。这个转换我一般写成一个独立的子函数function [A, rhs] build_ecopath_matrix(bio, pb, qb, dietMat, knownBio, targetEE) n length(bio); A zeros(n); rhs zeros(n, 1); for i 1:n % 图中已知生物量的功能组方程为 B_i*(P/B)_i sum(B_j*(Q/B)_j*DC_ji) EX_i if knownBio(i) A(i, i) pb(i); % 未知数为 EE_i 时B_i*(P/B)_i*EE_i 中 B_i*(P/B)_i 是系数 A(i, 1:n) A(i, 1:n) - (qb .* bio) .* dietMat(:, i); rhs(i) -bio(i) * pb(i) * targetEE(i) land(i); else % 未知生物量的功能组EE 给定把生物量移到未知量一侧 A(i, i) pb(i) * targetEE(i); for j 1:n if j ~ i A(i, j) A(i, j) - qb(j) * dietMat(j, i); end end rhs(i) land(i); end end end这段代码的思路是每个功能组的平衡方程都写成待求变量的一次式。已知生物量时EE 是待求量落在对角线未知生物量时对角线上是 pb 乘预设 EE其余列是被该捕食者摄食带来的生物量流出。需要注意的是dietMat(:, i)代表所有捕食者对 i 的摄食比例按列索引。实际运行前建议打印 A 和 rhs 的规模再检查cond(A)——如果状态数过大通常意味着存在几乎不与其他组相连的孤立功能组或某两组的食性矩阵高度共线。这个预检能省掉后面很多调平衡的时间。3. 主迭代与 EE 收敛质量平衡求解器的 MATLAB 实现线性方程组 A * x b 的求解在 MATLAB 里一行代码就能完成x A \ b但 Ecopath 的真实求解过程没有那么直接。原因在于模型中有一部分功能组的 EE 值在生物学上并不自由——EE 不能大于 1。如果一个功能组的 EE 算出来是 1.3说明该系统无法支撑该组的当前生物量这组数据在平衡意义上是不成立的。所以求解器要做的不仅是解方程还要检查 EE 是否落在 (0, 1] 区间内。超过上限的组需要调整系统输入参数通常是降低其 P/B 或生物量然后重新求解。这是一个迭代逼近的过程不是一次性解方程。3.1 先固定预平衡再求全系统解资源包中REco_Flatfish1.csv与REco_Roundfish1.csv的数据差异恰好能说明迭代的必要性不同渔获物代表的资源调查往往给出相互矛盾的摄食压力估计。我的求解器采用了两阶段策略。第一阶段把所有捕食者的 Q/B 值和生物量视为已知对每个被捕食功能组独立解方程求初始 EE第二阶段将初值代入全局矩阵针对残差最大的功能组做参数松弛比如微调其 qb 或 diet 比例然后重新组装 A 矩阵求解。两个阶段交替直到所有 EE 均小于等于 1 且系统总残差低于阈值maxIter 20; tol 1e-6; for iter 1:maxIter [A, rhs] build_ecopath_matrix(bio, pb, qb, dietMat, knownBio, fixedEE); x A \ rhs; % 将解映射回生物量和EE ee ones(n,1); bioNew bio; for i 1:n if knownBio(i) ee(i) x(i) / bio(i) / pb(i); % 标准化 else bioNew(i) x(i); end end if all(ee 1.0 1e-3) all(ee 0) res max(abs(A * x - rhs)); if res tol break; end else % 定位越界组 badIdx find(ee 1.0); % 对这些组的qb做向下微调幅度2% qb(badIdx) qb(badIdx) * 0.98; end end这里 EE 标准化那一步是较容易出错的细节。因为 A 矩阵对角线在 knownBio 分支上直接放的是pb(i)而非bio(i) * pb(i)所以解出的 x(i) 实际是bio(i) * EE(i)需要除以 bio 还原。这种方式比直接构造大矩阵编码 EE 更不容易写错。另外很多 MATLAB 入门教程会用inv(A)求逆来解方程组这里必须强调A \ b使用的是 LU 分解加列主元数值稳定性远高于显式求逆尤其在 A 不是正定矩阵时。矩阵求逆是 O(n³) 且对舍入误差敏感而左除运算符会根据矩阵结构自动选择最优分解算法。3.2 收敛判据与实际迭代效果的对照迭代过程中我把残差记录到数组里便于观察收敛行为residuals(iter) max(abs(A * x - rhs)); fprintf(iter%d maxEE%.3f resid%.2e\n, iter, max(ee), residual(iter));实际跑的时候我观察到残差一般会在前四次迭代内跌到 1e-7 以下说明线性求解本身没有数值问题。真正的瓶颈是 EE 越界带来的 qb 调参。如果某个功能组的 EE 反复大于 1单靠 qb 整体缩放效率很低。更好的做法是修改dietMat中该组作为捕食者那一行的摄食比例把摄食压力分散到其他猎物上。所以程序里应该在 badIdx 分支中同时调整两个向量而不是只降 qb。如果调整了 20 轮还不收敛基本可以断定是输入数据存在系统性问题比如某猎物的生物量低于该生态系统实际水平或某捕食者的 Q/B 值超出了该类群的生理上限。这时候我会把该组标记为 “unbalanced”导出 JSON 诊断信息而不是强行用数值手段把结果凑出来。4. 输出重算与网络图可视化营养级、混合营养影响与食物网连边平衡表收敛之后Ecopath 的输出不只是 X 向量本身。更常用的指标包括有效营养级Effective Trophic Level、混合营养影响矩阵Mixed Trophic ImpactMTI、以及各功能组之间的能量流矩阵。其中营养级的计算依赖摄食矩阵和流量TL_i 1 Σ_j DC_ij * TL_j。这里需要注意公式里用的 DC 是捕食者 i 对猎物 j 的摄食比例方向与原始 dietMat 互为转置初学者容易弄反。MATLAB 里计算过程如下TL ones(n, 1); for iter 1:20 TLOld TL; for i 1:n prey find(dietMat(i, :) 0); if isempty(prey) TL(i) 1; % 生产者或碎屑组 else TL(i) 1 sum(dietMat(i, prey) .* TLOld(prey)); end end if max(abs(TL - TLOld)) 1e-4, break; end end这段代码采用迭代法解线性方程组而不是直接求逆。原因在于 TL 的定义本质上是一个马尔可夫链的稳态分布迭代法天然满足收敛性且能容忍 DC 矩阵中少量非归一化行。收敛后通常浮游植物营养级为 1滤食性鱼类在 2.2~2.6顶级捕食者能到 4 以上。我用 REco_Roundfish1.csv 跑出来的结果显示圆尾鲀类功能组的有效营养级在 3.8 左右说明这个区食物网中它对中级鱼类的依赖很强这个数字可以直接写进论文方法部分。4.1 从邻接矩阵到网络图有向加权图的 MATLAB 绘制食物网的可视化重点是展示能量流方向需要按流量加权设置线宽和透明度。MATLAB 的digraph函数能直接接受邻接矩阵并生成有向图对象然后通过plot控制布局。我习惯用layout参数设为force这样捕食者会被自动推向外圈生产者聚在中心视觉上类似 Ecopath with Ecosim 的网络图% 流量矩阵 FF(i,j) 表示从 j 流向 i捕食者 i 消费猎物 j F (qb .* bio) .* dietMat; % 流量 消费量 * 生物量 * 摄食比例 G digraph(F, speciesName); p plot(G, Layout, force, LineWidth, 1 5 * normalize(edges(G).Weight, range)); p.MarkerSize 7; labeledge(p, 1:numedges(G), round(edges(G).Weight, 2));这段代码里 normalize 函数将权重映射到 [0,1] 再乘以缩放系数避免大流量边把图撑爆。运行时我建议先看summary(G)确认节点数量和孤立节点如果有孤立节点说明 CSV 里存在没有被任何摄食关系连接的功能组通常是碎屑组或浮游生物组没做好关联。这个检查很重要因为 Ecopath 模型允许把碎屑组设为独立功能组但碎屑组几乎不被任何捕食者引用时网络会出现断环很多后续指标会失真。4.2 生态网络指标系统吞吐量与平均传输效率MTI 矩阵不仅能看出直接影响还能揭示间接效应。MTI 的定义基于 Leontief 逆矩阵而 Leontief 逆矩阵的计算本质是凯恩斯乘数在生态学中的翻版这也是为什么经济学家转行做生态模型往往上手很快。MATLAB 中实现方式足够直接I eye(n); S (I - F ./ sum(F, 2)) \ I; % 结构矩阵 % 按捕食者-猎物对计算混合营养影响 mti -S .* (bio * (1 ./ sum(F, 2)));计算后输出为热图时我常把对角线强制设为零因为 MTI 对角线表示的是密度依赖效应在大多数分析中会被单独解读。热图色标建议用parula它对色盲读者更友好且黑色背景下的对比度好于 jet。除了这些指标也要合理地引入系统总吞吐量 TSTTotal System Throughput。其计算是sum(F(:)) sum(bio .* qb)它能反映整个食物网的规模是不同生态系统间对比的基准量。用表格整理核心输出参数时一般会列出营养级、EE、P/Q 比值和流量占比。下表是 REco_Flatfish1.csv 预处理后某一功能组平衡成功后的输出功能组生物量 (t/km²)P/BQ/B营养级EE鳎科4.281.857.623.420.86牙鲆2.172.108.153.880.92虾类9.364.2016.82.510.79这张表直接回答了“平衡后的系统长什么样”这个问题。注意 EE 值是质量平衡后的核心输出它代表了该功能组被上层捕食和被渔业利用的总比例。EE 接近 1 说明该组几乎被“吃干榨净”后续做管理策略模拟时需要重点关注。5. 平衡失败时的定位技巧从 EE 越界追踪到食性矩阵的污染行调试 Ecopath 模型最花时间的不是算法本身而是弄清楚为什么某个功能组的 EE 死活大于 1。这部分技巧可以直接用于任何质量平衡类生态模型包括营养盐收支模型。dietMat的污染行是最隐蔽的坑。食物网数据来自文献汇总不同文献对同一物种的食性描述可能差异很大。资源包里的REco_Flatfish2.csv是从另一种渔获物调查中整合的其中某个捕食者的猎物比例总和达到 1.13——这显然是有机碎屑和浮游植物的比例加超了。但注意问题在于这个捕食者本身可能不是系统中最关键的那个但它会通过摄食压力传导影响多个猎物的 EE。所以排查顺序要反过来先列所有 EE 1 的功能组再列出这些组的所有捕食者检查这些捕食者的摄食比例是否归一化。dietRowSum sum(dietMat, 2); badPred find(abs(dietRowSum - 1) 1e-3); for i badPred fprintf(捕食者 %d: 摄食比例和%.3f, i, dietRowSum(i)); % 找出是哪个列猎物贡献最多 [val, idx] max(dietMat(i, :)); fprintf( 最大项: 猎物%d 占比%.2f\n, idx, val); end如果发现确实有多个组摄食比例超标修复时不要盲目归一化。归一化会改变功能组间的相对摄食压力导致原本平衡的组反而失衡。正确做法是检查原始数据来源看超标是否来自碎屑组。碎屑组在 Ecopath 中往往被同时当作猎物和汇很多模型为了闭合能量平衡会在碎屑组上叠加额外的呼吸损失。这类数据不能简单归一化需要把多余比例转移到碎屑组的“不可利用”系数里相当于增加系统呼吸。另一个高价值调试指标是功能组的呼吸量估算R Q - P - 未同化部分。MATLAB 里可以按 P/Q 比值直接推算。如果算出的呼吸量为负通常说明该组的 Q/B 值低于维持代谢所需说明系统的摄食数据偏低。这个负呼吸的检查比 EE 的检查更早暴露问题因为呼吸量是从热力学约束入手的不需要等矩阵求解结束。最后提一个个人习惯我会在平衡表求解完成后把总生物量、总消费量、总呼吸量分别打印出来然后手动算一下系统总净生产力是否大于总呼吸消耗。如果不大于说明这个系统处于衰败状态——这在数学上模型照样能收敛但生态学上没有意义。这种交叉验证能帮你筛掉很多“数字上平衡、理论上荒谬”的模型结果避免在后续基于此模型的渔业管理建议中翻车。本文还有配套的精品资源点击获取