从Arrhenius速率到元胞自动机：数学建模连续离散映射-育师

简介：面向备战美国大学生数学建模竞赛的高校学生与指导教师，这份2021年美赛特等奖论文合辑提供获奖队伍的完整论文原文，可用于研读赛题拆解、模型构建与英文写作规范。合辑为1个PDF文件，压缩包约54.28MB，便于集中查阅与打印学习。其中收录的获奖论文以“Fungal Decomposition Modeling via Stochastic Gradient Descent and Cellular Automata”为例，围绕真菌生长与分解效率受温湿度及物种相互作用影响的问题，依次建立扩展速率模型、非线性分解速率优化模型，并用随机梯度下降求解参数，再以细胞自动机构建多物种生长与分解仿真。读者可从中学习阿伦尼乌斯关系建模、优化算法收敛性验证、竞争关系刻画、敏感性分析及122天分解率等结果讨论，理解统计、优化与计算方法如何衔接。已有4027人学习，适合希望对照特等奖范文提升建模深度与论文表达能力的参赛者。

1. 两个数字背后的两次翻译

拿到一份数学建模竞赛特等奖论文，很多人的第一反应是“公式一定很复杂”。把这篇作品逐节拆开会发现，结论只靠四个东西撑着：一个阿伦尼乌斯项、一条对数线性关系、一个 1/2 阶幂次，再加一台格子自动机。122 天分解率 0.68，引入竞争后掉到 0.39、降幅 42.69%，这两个数字之所以站得住，关键不在于拟合得多准，而在于连续方程和离散格子之间的两次转换没有走样。这篇要讲的就是这两次翻译：速率怎么变成概率，概率又怎么变回一个可比的分解率。

2. 速率模型：Arrhenius 与水分权衡怎么拼成一个非线性方程

整篇建模的地基是三段关系，而不是一个大公式。原文的处理顺序很讲究：先单独立扩展速率模型，再在它上面搭分解速率模型，最后才谈两者之间的定量耦合。这个顺序不是行文习惯，而是因为三个量所处的物理层级不同：温度、湿度是环境输入，扩展速率是菌丝的响应，分解速率是底物被转化的结果。层级错位一步，后面 SGD 拟合出来的参数就没有物理含义，只剩一组能让曲线贴合的数值。

2.1 温度项：为什么是指数的指数

扩展速率被写成对温度的 Arrhenius 关系，形式是k = A·exp(-Ea/(RT))。很多模型也会用温度线性项或者二次项，那是把酶活性当成一个可以从外部随意调控的连续量。Arrhenius 的选择意味着把反应速率交给分子碰撞能量分布：只有能量超过活化能Ea的那部分分子能参与反应，而超过阈值的比例随温度按指数增长。

关键在于T出现在分母上。低温区间里T一点点上升，-Ea/(RT)的变化幅度很大，速率爬升很快；到了高温区间，分母变大，同样幅度的升温带来的相对变化被急剧压缩，曲线趋于饱和。这条曲线天然是“两头不对称”的。原文的灵敏度分析说模型对 8% 的温度和湿度波动敏感，其实一半原因就藏在这个分母里：在低温段，8% 的温升换来的速率变化远大于高温段的 8%。如果用线性项，这份敏感性会被抹平，后面所有关于气候带差异的结论都会失真。

这也是为什么模型必须先建速率、再建分解：Arrhenius 只描述“菌丝能长多快”，不负责“长出来之后分解了多少”。把这两件事塞进同一个式子，温度的双重作用就没法分离。

2.2 水分权衡：取对数之后才是一条直线

原文里有一句很容易被略过的话：水分权衡与对数分解速率之间是线性关系。字面看平平无奇，翻成数学语言却是ln D = α·M + β，等价于D ∝ exp(α·M)。也就是说，湿度对分解效率的作用不是“多一点水多一点分解”的线性叠加，而是每一份新增水分都在已有基础上按比例放大。

这符合真菌分解的实际约束。湿度低的时候，菌丝在干基质里移动受限，胞外酶扩散不出去，分解近乎停滞；湿度跨过某个下限之后，水膜连通了菌丝与底物之间的通道，酶的输运效率陡然提升；再往上，水分开始挤占孔隙中的氧气，好氧分解反而被抑制。这条曲线在中间段近似指数上升，在两端被压住。原文只说在“湿度下限”之上讨论，其实就是在避开低湿端的平台区，只拟合有效区间内的指数段。

对数变换在这里的作用是让非线性关系可以用线性工具检验：把观测的分解率取对数，如果对水分权衡量能画出一条直线，就说明指数假设成立。这也解释了为什么模型对湿度 8% 的波动同样敏感：指数关系的相对误差会被放大成绝对误差。

2.3 三条关系怎么合成一个方程

现在可以拼了。原文说扩散律与分解律之间是 1/2 阶关系，常见的解读是：分解只发生在菌丝网络的边界上，被转化的底物量与网络的周长同阶，而周长随面积按平方根增长，面积又正比于扩展速率的累积。于是分解速率对扩展速率取 1/2 次幂。把这一项和上一节的指数水分项相乘，再补一个基线常数，就得到了一个同时满足两条经验关系的结构：

import numpy as np R_GAS = 8.314 # 气体常数, J/(mol*K) def extension_rate(T_c, H, A, Ea, k_h): """扩展速率 = Arrhenius 温度项 * 湿度饱和型修正。 T_c: 摄氏温度; H: 相对湿度 0~1 A: 指前因子; Ea: 活化能; k_h: 湿度半饱和系数""" T_k = T_c + 273.15 temp_term = A * np.exp(-Ea / (R_GAS * T_k)) # 温度在分母，低温段陡、高温段平 moist_term = 1.0 - np.exp(-k_h * H) # 随 H 单调上升并趋于饱和 return temp_term * moist_term def decomposition_rate(E, M, theta): """分解速率 = 扩展速率的 1/2 阶 * 水分权衡的指数项 + 基线。 M: 湿度下限之上的可用水分量; theta = (a, b, c)""" a, b, c = theta return a * np.sqrt(E) * np.exp(b * M) + c

代码里三个系数的分工很明确：np.sqrt(E)负责 1/2 阶尺度关系，np.exp(b * M)负责“取对数后是直线”的水分权衡，c是分解停滞时的基线残留。这个结构的好处是对数一取就现原形：ln(D - c) = ln a + 0.5·ln E + b·M，对E是 0.5 次、对M是 1 次，两条经验关系同时被满足，不需要再额外加修正项。

extension_rate里用1 - exp(-k_h*H)而不是直接乘H，是因为湿度到一定程度后再加水，扩展速率不再线性增长。这个饱和形式在参数估计时还有一个好处：k_h的取值直接反映了“这种菌对干燥的耐受阈值”，是一个可解释的生物学参数，而不是纯数学拟合出来的数。

提示：decomposition_rate里的c不能随意设成 0。分解过程的末端往往有难降解组分（木质素一类），残留量构成一个非零基线，强行置零会让a和b在高分解率区间被系统性拉偏。

3. SGD：非线性参数估计为什么在这里收敛得住

结构定下来之后，剩下的问题是估(a, b, c)。这一步看起来最像“调包”，实际上是最容易翻车的地方：模型对参数是非线性的，梯度不是常数，初值稍偏就可能卡在平坦区或者震荡发散。原文强调算法“稳定的收敛行为证明了有效性”，这句话的分量比表面听起来重——在非线性拟合里，稳定收敛本身就是一条关于目标函数地形的重要证据。

3.1 目标函数：残差平方和为什么够用

拟合的做法是把收集到的分解率观测值和模型预测值的差平方求和，作为待最小化的目标。选择平方误差而不是绝对误差，一个重要原因是它在参数空间里连续可导，梯度处处存在，SGD 的更新式才有定义。另一个原因是平方误差对大偏差的惩罚更重，而分解率数据的噪声主要来自测量误差，极端值往往就是测量出了问题，被重罚一次反而能把参数拉回物理合理区间。

但平方误差在非线性结构下会带来一个副作用：目标函数可能不是凸的。a和b分别是乘性和指数位置上的参数，两者存在耦合——a调大、b调小，在某些数据区间里可以给出几乎一样的预测。这种“参数不可辨识”是黑匣子拟合的经典问题，也是为什么原文必须给参数加上物理约束（活化能为正、湿度系数有界）的原因。不加约束，SGD 会沿着那条平坦的谷底一路滑到没有物理意义的角落。

3.2 随机梯度的“随机”到底随机在哪

批梯度下降每次更新都要扫一遍全部数据，样本量大时单步成本高；SGD 每次只取一小批，用这批数据的梯度近似整体梯度。这个近似引入了噪声，噪声的方向是随机的，但期望等于真实梯度。下面这段更新逻辑把三个参数的偏导都写开了：

def sgd_fit(records, theta0, lr=5e-3, epochs=800, batch=16): theta = np.array(theta0, dtype=float) hist = [] for ep in range(epochs): order = np.random.permutation(len(records)) # 每轮重排，避免批次相关性 for s in range(0, len(records), batch): idx = order[s:s+batch] E = np.array([records[i]['E'] for i in idx]) M = np.array([records[i]['M'] for i in idx]) y = np.array([records[i]['D'] for i in idx]) pred = decomposition_rate(E, M, theta) resid = pred - y # 残差，注意顺序决定梯度符号 base = np.sqrt(E) * np.exp(theta[1] * M) g_a = np.mean(resid * base) # 对 a 的偏导 g_b = np.mean(resid * theta[0] * base * M)# 对 b 的偏导，多一个 M g_c = np.mean(resid) # 对 c 的偏导 theta -= lr * np.array([g_a, g_b, g_c]) hist.append(float(np.mean(resid ** 2))) return theta, hist

batch取 16 而不是 1，是在噪声和效率之间找的平衡点：纯单样本更新的噪声太大，收敛曲线会剧烈抖动，看不出趋势；批次太大又退化成批梯度下降，丧失跳出局部平坦区的能力。lr取 5e-3 这个量级，是因为三个参数的量纲差了好几个数量级——a可能是几十，b是个位数，c接近 0，统一学习率下步子太大会让c直接越过最优值。

g_b的表达式里多乘了一个M，这是链式法则的直接结果：b在指数位置上，求导后会把M带下来。这一项在湿度高时会被显著放大，意味着高湿样本对b的估计影响更大。如果数据集中高湿样本很少，b的梯度信号就弱，收敛会停在离真值较远的地方——这是数据分布问题，换优化器解决不了。

3.3 收敛曲线能证明什么、不能证明什么

原文用“稳定收敛”作为算法有效性的证据，这个论证的边界需要说清楚。一条单调下降最后趋于平的损失曲线，能说明的是：学习率没有设得太大导致发散，梯度方向没有被量纲问题带偏，参数在迭代中没有跑出可行域。它不能说明找到的是全局最优，也不能说明参数有唯一解。

在非线性最小二乘里，判断结果是否可信的常用做法是换几组初值重跑，看最终参数是否收敛到同一片区域。如果多组初值给出差别很大的(a, b, c)但损失几乎相同，问题就不在优化器，而在模型结构本身存在冗余项。另一种验证是把数据分成两段，用一段拟合、一段检验，只有两段的残差量级相当，收敛才真的有意义。

注意：损失曲线的纵轴如果取了绝对值很小的量级（比如 1e-6），下降看起来会很平滑，实际上相对误差可能还很大。比较收敛质量要看相对残差而不是绝对损失值。

4. 元胞自动机：连续速率怎么翻译成离散规则

到这里手上有一个能算速率的连续方程，但它回答不了“40 天后为什么稳定”“竞争之后总数为什么下降三成”。这类关于空间格局和种群动态的问题，连续方程要给答案就得解偏微分方程组，边界条件难定，而且很难体现“某个格子被某个菌种占了”这种离散事实。原文换成元胞自动机，本质是把连续的速率翻译成每个格子上每一步的转移概率。

4.1 格子、邻域与状态：把菌丝网络压成二维场

元胞自动机的三要素是格子、状态、邻域规则。这里的格子对应底物表面的一块微元，状态分两类：空（未被占据）和被某个菌种占据（有多个菌种时状态取值就扩展成 1 到 k）。邻域取摩尔邻域还是冯·诺依曼邻域，直接影响扩展前沿的形状——菌丝是向外辐射生长的，八邻域在二维上更容易还原这种辐射感，四邻域会人为压出方块状的生长前沿。

底物的量不能塞进状态里，需要单独开一个场记录每个格子的剩余养分量。这一点很容易被忽略：如果把“养分量”直接当成状态值，状态空间会变得连续，元胞自动机那种“局部规则反复迭代”的性质就被破坏了。原文的做法是状态只记占据者，消耗量记在旁路数组里，两者同步更新。

时间步长和格子尺寸是一对绑定参数。格子划得越细，单步能跨越的物理时间越短，否则菌丝的扩展速度会被放大成“瞬移”。判断标准是：一个时间步内，扩展前沿推进的距离不应超过一个格子的边长。原文能跑出“前 5 到 10 天数量显著增加、约 10 天后回落、约 40 天稳定”这条曲线，说明步长和格子尺寸的比值是经过校准的，而不是随手填的。

4.2 两个映射：速率变概率，概率变面积

从连续方程到离散规则有两次转换，两次都必须单调、可逆，否则累积出来的分解率就没有意义。

第一次是把扩展速率映射成定殖概率。做法通常是拿扩展速率除以一个参考速率做归一化，得到 0 到 1 之间的数，作为空格子被邻胞定殖的概率。这里要注意，概率不能等于速率本身，速率再大也不该保证一定定殖，否则扩展前沿会变成刚性推进，失去随机性带来的粗糙边缘。

第二次是把分解速率映射成每步的底物消耗量。这一步需要单位换算：连续模型给出的分解速率是“每单位时间转化多少比例的底物”，元胞自动机需要的是“每一步扣掉多少绝对量”。两者之间差一个时间步长因子，漏掉这个因子，122 天的分解率就会差出一个数量级。

def ca_step(grid, nutrient, ext, dec, rng, species=3): """grid: 状态矩阵, 0 为空, 1..species 为占据者 nutrient: 底物场，与 grid 同形 ext: 各菌种归一化后的定殖概率; dec: 各菌种每步的底物消耗量""" h, w = grid.shape new = grid.copy() for i in range(1, h - 1): for j in range(1, w - 1): if grid[i, j] == 0: # 空格子：看邻胞能否定殖 for s in range(1, species + 1): if has_neighbor(grid, i, j, s) and rng.random() < ext[s]: new[i, j] = s break else: # 已占据：消耗底物 s = grid[i, j] nutrient[i, j] -= dec[s] if nutrient[i, j] <= 0: new[i, j] = 0 # 底物耗尽，退出该格 return new

ext[s]是概率，dec[s]是绝对消耗量，两个数组的量纲完全不同，不能混用。nutrient用独立数组记录而不是拼进grid，是为了让状态空间保持离散有限，迭代不会因为浮点累积而漂移。整体分解率则由初始底物总量和当前剩余量做归一化得到——122 天 0.68 这个数，就是累积消耗量除以初始总量。

提示：has_neighbor遍历八邻域时要注意边界，越界索引在 Python 里会静默回绕成负索引，导致对侧边缘被错误地当成邻胞，生长前沿会莫名其妙地在对面出现。

4.3 竞争关系：33.65% 这个数字是被定义出来的

引入多物种之后，规则里要加入竞争。原文的做法很克制：没有去建一套资源争夺的动力学，而是先给“竞争的强度”下一个可计算的度量——用共存时某物种的扩展速率相对于它单独生长时的扩展速率的下降幅度，来刻画它被压制了多少。

这个定义的好处是把一个模糊的生态学概念压缩成一个标量，可以直接塞进元胞规则：被压制的物种，其ext[s]按这个比例打折。打折之后，同格子数下的定殖概率下降，前沿推进变慢，长期累积下来占据的面积就少。原文给出的结果是多物种共存时数量平均下降 33.65%，122 天分解率从 0.68 降到 0.39，降幅 42.69%。

两个降幅不一样，这本身是个信息量很大的细节。物种数量降了三成多，分解率却降了四成多，说明竞争的影响不是均匀摊在每个格子上的，而是集中在底物富集、生长最快的区域。这也解释了原文的另一个结论：短期内相互作用不明显，长期才显著——早期底物充足，各物种都在自己的邻域里扩张，重叠区小；后期养分紧张，重叠区扩大，打折规则才开始生效。

从机制上看，这套竞争建模牺牲了“物种为什么竞争”的细节，换来的是“竞争会导致什么”的可计算性。对竞赛类建模来说这是合理取舍，但如果要追问竞争优势的来源（是酶分泌能力、还是菌丝穿透力），这个定义给不出答案。

5. 常见问题排查：这套模型最容易翻车的四处

模型结构不复杂，真正让人卡住的是实现层面的细节。下面四条是按“现象 → 原因 → 解决”的顺序整理的高频问题，每一条都对应一个具体可查的环节。

5.1 现象：SGD 收敛曲线看着很平，但预测值系统性偏高

原因通常出在decomposition_rate里那个基线常数c。如果c被初值固定在一个偏大的位置，同时又没有对a加约束，优化器会倾向于用c去吸收残差，因为改c的梯度是常数、最容易下降。结果是c一路涨到接近观测均值，a被压到接近 0，模型退化成“常数预测”，损失看着降了，物理意义全丢了。

解决的办法是给参数加边界，并且在训练开始前先做一个粗略的量级估计：把观测数据里最大分解率和最小分解率之差作为a的参考量级，c的上界设成最小观测值的一半左右。另外可以监控a的轨迹，如果它在训练初期就单调跌向 0，说明结构或者数据有问题，继续跑没有意义。

5.2 现象：分解率算出来是 0.9 以上，明显不合常理

先查时间步长因子有没有漏。连续模型给出的分解速率是单位时间比例，元胞自动机每一步要乘上步长对应的物理时间。如果步长设成 1 天但实际按 0.1 天记录，累积消耗就会被放大十倍。这个错误不会报错，只会让结果看起来“效率极高”。

其次查底物场的初始化。如果nutrient每个格子的初值不统一，归一化时分母用的是初始总量，分子却按实际消耗算，两者口径不一致，结果就会失真。稳妥的做法是初始化时把总量记下来，之后所有比例都除以这个固定值，而不是动态求和。

还有一种情况是边界处理导致重复计算：遍历时用了range(1, h-1)排除边界，但初始化却把整个矩阵都填满了底物，边界格子永远不会被消耗，分母偏大，比例反而偏小——如果同时叠加了步长错误，两个偏差会部分抵消，结果看起来“差不多对”，实际上两个因子都错了。

5.3 现象：加了竞争之后，物种数量反而变多了

这个反直觉的结果多半是竞争规则的施加方向搞反了。竞争度量的定义是“扩展速率下降的幅度”，作用到ext[s]上应该是乘(1 - 降幅)。如果写成了乘(1 + 降幅)，被压制的物种反而获得了更高的定殖概率，长期下去它会挤占其他物种的空间。

另一个可能是竞争只在某些条件下触发，但条件判断写在了更新之后，导致同一轮里有的格子用了打折前的概率、有的用了打折后的，统计口径不统一。解决办法是把竞争系数的计算提到每轮迭代的最前面，做成一个只读的数组，整轮共用，避免迭代中出现状态不一致。

5.4 现象：换一组初值重跑，参数差很多但损失差不多

这不是 bug，是结构性的参数不可辨识。a和b在a·sqrt(E)·exp(b·M)里是乘性和指数位置的组合，当数据集中在某个窄区间时，不同的(a, b)组合可以给出几乎相同的预测。

处理办法有两条：一是检查数据覆盖范围，如果M的取值只在一个很窄的带里，b的梯度信号本来就弱，任何优化器都估不准，需要补数据而不是换算法；二是对b施加物理上合理的约束区间，把搜索空间压缩到有意义的范围，让最优解在约束内唯一。判断是否属于这种情况，最直接的方法是画(a, b)的损失等高线，如果能看到一条狭长的谷，就说明是结构问题。

6. 设计取舍：为什么用元胞自动机而不是偏微分方程

把连续速率换成格子规则，代价和收益都很明确。收益是离散性：每个格子记着“被谁占了”“还剩多少底物”，物种竞争、局部耗尽、前沿粗糙这些现象是规则的自然产物，不需要额外建模。代价是尺度依赖：格子尺寸和时间步长成了新的自由参数，结论对这两个数敏感，而连续方程里不存在这个问题。原文给出的 40 天稳定、122 天 0.68 这些具体数值，严格来说是“在这套离散参数下的结果”，换成更细的格子会得到略有差异的曲线。

真正值得琢磨的是竞争建模那一层的取舍。用“扩展速率下降比例”来刻画竞争，等于把复杂的资源争夺压缩成一个标量系数。这让模型可解、可比较、可跑出 33.65% 这样干净的结论，但也让模型失去了回答“谁为什么赢”的能力。原文最后能给出不同环境下的优势种，靠的是不同温湿度组合改变了各物种的扩展速率，竞争系数只是在这个基础上做缩放——优势种是“在打折之后仍然更快”的那个，不是“更擅长争夺资源”的那个。

验证这类模型，我一般看三件事：换初值后参数是否稳定收敛到同一区域；灵敏度和原文报告的量级是否一致（温度湿度 8% 的扰动应该带来可观测但不爆炸的变化）；以及去掉竞争项之后，结果是否能退化回单物种情形的数量级。第三点最容易被跳过，但它是最硬的一条——一个加了交互的模型如果退化不回无交互的基准，说明交互项里混进了别的东西。