1. 项目背景与问题引入当数学模型遇上公共卫生危机2022年初我参与了一场特殊的数学建模竞赛——认证杯SPSSPRO杯。当时全球仍笼罩在新冠疫情的阴影下我们团队抽到的C题第二阶段直接将我们拉入了一个极具现实意义的战场如何利用“污水流行病学”的原理为疫情防控提供量化的决策支持这个题目让我印象深刻因为它完美地诠释了数学建模如何从象牙塔走向街头巷尾解决真实世界的棘手问题。你可能听说过通过检测污水来追踪毒品滥用情况这就是污水流行病学的经典应用。其核心逻辑很简单人体代谢后的许多物质会随排泄物进入城市污水系统。通过检测污水中这些特定生物标志物比如病毒RNA片段、药物残留的浓度结合污水流量、服务人口等数据就能反推出一个社区内某种物质的消费量或某种病原体的流行情况。这种方法具有匿名、客观、覆盖面广、可早期预警的巨大优势。在新冠疫情中它被用来监测社区感染水平甚至早于临床报告发现疫情苗头。然而题目给我们的不是一篇综述而是一系列冰冷的数据和几个尖锐的问题。它模拟了一个城市的污水监测数据要求我们建立数学模型估算真实感染人数分析数据中的异常波动并评估不同防控措施的效果。这听起来像是一个完美的“黑箱”预测问题但实际操作中从污水浓度到感染人数中间隔着人口流动、病毒载量差异、检测灵敏度、管道沉降、水力停留时间等无数个“灰箱”。题目提供的“校正因子”关键词正是打开这些灰箱的一把关键钥匙。整个解题过程实际上是一场关于“如何用数学语言描述并量化不确定性”的思维体操。2. 核心问题拆解从污水数据到管理决策的三重挑战面对题目给出的数据集和问题描述我们首先需要摒弃“直接套公式”的幻想。污水流行病学应用于新冠监测其建模链条比传统应用更为复杂。我们将赛题的核心诉求拆解为三个环环相扣的层次这构成了我们整个论文的骨架。2.1 第一层反演问题——如何从浓度估算感染人数这是最根本的数学模型问题。输入是污水中新冠病毒RNA的浓度单位拷贝数/升输出是监测区域内可能的活跃感染人数。这个过程绝非简单的线性比例计算。我们需要构建一个“反演模型”。基本公式框架如下[ N_{infected} \frac{C_{wastewater} \times F \times V}{f \times R \times S \times \eta} ]其中(N_{infected})估算的感染人数。(C_{wastewater})污水中测得的病毒RNA浓度拷贝数/L。(F)污水稀释因子考虑地下水渗入、雨水等。(V)污水厂日均处理流量L/天。(f)每人每天平均粪便排放量L/天/人。这是一个生理常数但个体差异大。(R)感染者日均通过粪便排出的病毒RNA拷贝数。这是最不确定、最核心的参数它取决于病程阶段、个体免疫状况、病毒毒株等多种因素。(S)病毒RNA在污水输送和处理过程中的存活率/回收率。病毒RNA会降解也会被吸附在管道污泥上。(\eta)实验室检测方法的回收率/效率。题目中提到的“校正因子”本质上就是针对 (R)、(S)、(\eta) 等这些高度不确定的参数进行整体性或分项式的标定与修正。我们的首要任务就是设计一种算法或模型来合理地确定或校准这个综合校正因子 (K 1/(f \times R \times S \times \eta))使得模型估算结果与已知的、有限的其他数据如同期临床报告病例数在统计意义上最为吻合。注意这里切忌直接使用文献中的单一数值。不同城市、不同季节、不同毒株、不同污水管网材质这些参数都可能不同。建模的价值就在于根据本地数据动态地确定或校准这些参数。2.2 第二层时空分析——如何识别和解释数据异常题目提供的污水监测数据通常包含时间序列每日/每周浓度可能还有空间序列来自城市不同采样点的浓度。第二层挑战在于分析这些数据的波动模式。趋势分析感染人数是持续上升、平台期还是下降这对应疫情的发展阶段。异常点检测某天的病毒浓度突然飙升或骤降是真实的疫情暴发/管控起效还是采样误差、检测异常、或大型活动如音乐节、体育赛事导致的人口临时聚集造成的我们需要利用统计方法如3σ原则、移动平均法、时间序列分解识别异常点并结合收集到的“上下文信息”如是否举办大型活动、是否发布新的防控政策进行归因分析。这一步是将数据转化为情报的关键。空间异质性分析如果有多点数据可以分析不同区域如商业区、居民区、大学城的病毒浓度差异。这能帮助定位高风险区域实现精准防控。2.3 第三层措施评估——如何量化防控策略的效果这是模型的最终应用层。题目可能会假设几种情景例如情景A保持现有措施。情景B实施全员核酸筛查。情景C对特定区域进行封控。情景D提倡居家办公、关闭娱乐场所。我们需要利用建立好的模型来预测不同措施下污水病毒浓度及反演的感染人数随时间的变化。这通常需要引入机制模型如经典的传染病模型SIR, SEIR与我们的反演模型进行耦合。例如封控措施降低了人群接触率这会在SEIR模型中体现为传播系数β的下降进而改变每日新增感染人数I(t)最终映射到模型估算的污水病毒浓度上。通过比较不同情景下的预测曲线如感染峰值、达峰时间、总感染人数可以定量评估各种措施的成本效益。3. 模型构建与求解一个融合了统计与机制的混合框架基于以上拆解我们团队没有采用单一的模型而是设计了一个“分阶段、混合式”的建模框架。这个框架的核心思想是先利用数据驱动方法校准基础模型再引入机制模型进行推演。3.1 阶段一基于贝叶斯推断的校正因子动态校准模型面对反演公式中巨大的参数不确定性我们选择了贝叶斯统计模型作为核心工具。它的优势在于能将我们对参数的不确定性先验知识和观测数据的不确定性似然函数结合起来得到参数完整的概率分布后验分布而不仅仅是一个点估计。模型设定观测数据时间序列的污水病毒浓度 (C_t)以及可能不可靠但可作为参考的同期临床新增报告病例数 (I_t^{reported})。待估参数综合校正因子 (K_t)。我们允许它随时间缓慢变化以捕捉病毒排毒规律、检测效率等随时间的变化。我们假设 (\log(K_t)) 服从一个随机游走过程。状态空间模型状态方程描述参数如何演化(\log(K_t) \log(K_{t-1}) \epsilon_t, \quad \epsilon_t \sim N(0, \sigma^2))。观测方程描述数据如何生成对于污水数据(C_t (N_t^{infected} / K_t) \nu_t)。这里 (N_t^{infected}) 需要从一个简单的传播模型如SIR初步估算或者初期可假设与 (I_t^{reported}) 成比例。(\nu_t) 是观测误差。对于临床数据如果可用(I_t^{reported} \sim \text{Binomial}(N_t^{infected}, \rho))。其中 (\rho) 是报告率也是一个待估参数。求解与实现 我们使用马尔可夫链蒙特卡洛MCMC方法通过Python的PyMC3或Stan库来实现贝叶斯推断。代码的核心是定义上述概率模型然后让采样器去探索参数的后验分布。# 伪代码示例展示PyMC3模型结构思路 import pymc3 as pm with pm.Model() as wastewater_model: # 先验分布对log(K)的初始值和波动性进行假设 sigma_K pm.HalfNormal(sigma_K, sigma0.5) # K随时间变化的波动大小 K_log pm.GaussianRandomWalk(K_log, sigmasigma_K, shapelen(days)) K pm.Deterministic(K, pm.math.exp(K_log)) # 转换为正数 # 报告率的先验 reporting_rate pm.Beta(reporting_rate, alpha2, beta8) # 假设报告率较低先验集中在0.2左右 # 连接感染人数与污水浓度 (简化版假设感染人数已知为I_estimated) # I_estimated 可以从一个独立的SEIR模型生成作为输入 expected_concentration I_estimated / K # 似然函数假设观测浓度服从正态分布方差未知 conc_sigma pm.HalfNormal(conc_sigma, sigma10) observed_conc pm.Normal(observed_conc, muexpected_concentration, sigmaconc_sigma, observedreal_wastewater_data) # 连接感染人数与临床报告病例 (如果数据可用) reported_cases pm.Binomial(reported_cases, npm.math.ceil(I_estimated), # 感染人数取整 preporting_rate, observedclinical_reported_data) # 运行MCMC采样 trace pm.sample(2000, tune1000, cores2, return_inferencedataFalse)运行这个模型后我们得到的不是单一的K值而是K随时间变化的一系列可能轨迹后验分布样本。我们可以用其中位数轨迹作为最佳估计同时其分位数如95%置信区间清晰地展示了估算的不确定性。这比给出一个孤零零的数字要有力得多。3.2 阶段二耦合SEIR模型进行预测与措施模拟在阶段一我们可能用一个简单的趋势模型来生成 (I_estimated)。为了进行措施评估我们需要一个能反映传染病动力学的机制模型。我们引入一个经典的SEIR模型[ \begin{aligned} \frac{dS}{dt} -\frac{\beta(t) I S}{N} \ \frac{dE}{dt} \frac{\beta(t) I S}{N} - \sigma E \ \frac{dI}{dt} \sigma E - \gamma I \ \frac{dR}{dt} \gamma I \end{aligned} ]其中(S, E, I, R) 分别代表易感者、潜伏者、感染者、康复者(N)为总人口。(\sigma)是潜伏期倒数(\gamma)是感染期倒数。关键参数是传播率 (\beta(t))它是一个可以受防控措施影响的变量。措施如何量化我们将不同的防控措施映射为对 (\beta(t)) 的调制封控/静默大幅降低接触使 (\beta) 降至一个很低的基础水平如原来的20%。社交距离/口罩令中等程度降低传播效率使 (\beta) 减少一定比例如原来的40%-60%。疫苗接种通过降低易感者比例S和可能降低感染者传染力来间接影响有效传播。可以在模型中加入疫苗覆盖率变量或直接调整有效接触率。模型耦合与预测流程参数校准利用历史数据包括阶段一校准后的感染人数估算值使用最小二乘法或MCMC来拟合SEIR模型的参数主要是初始的 (\beta_0)。情景模拟设定一个时间点 (t_{intervention})从此处开始将 (\beta(t)) 的值根据预设措施进行调整。运行模型数值求解SEIR微分方程组得到未来一段时间内每日的感染人数 (I(t))。反演回污水浓度利用阶段一得到的最优校正因子 (K)或其后验均值通过公式 (C_{predicted}(t) I(t) / K)预测未来污水病毒浓度的变化趋势。效果对比绘制不同情景下的 (I(t)) 和 (C_{predicted}(t)) 曲线比较峰值高度、达峰时间、曲线下面积近似总感染人数等指标。4. 求解全过程的关键细节与避坑指南纸上谈兵总是容易但将上述框架实现并得出可靠结果过程中布满了“坑”。以下是我们在解题和后续复盘时总结的几个关键细节和避坑点。4.1 数据预处理异常值处理与平滑原始的污水监测数据噪声极大。某天浓度奇高可能是因为采样时正好有一段高浓度污水团流过某天浓度奇低可能是检测失败。直接使用原始数据拟合模型会导致灾难性后果。我们的做法滑动窗口中位数滤波对于时间序列数据我们首先采用滑动窗口中位数例如窗口大小为7天进行平滑。中位数比均值对异常值更不敏感。这能有效滤除短暂的尖峰或低谷保留趋势。基于模型的异常检测在建立了初步的模型如简单的趋势模型贝叶斯反演后计算每个数据点的标准化残差。将残差绝对值大于3倍标准差的数据点标记为“待审查异常点”。上下文归因查阅这些异常点对应日期的“上下文信息表”题目若提供或根据常识推断。如果发现某天有大型集会则该高点可能是真实的如果无特殊事件则很可能是技术误差可以考虑用前后数据的插值进行替换或在模型中赋予该点一个更大的观测误差方差降低其对参数估计的影响。4.2 校正因子K的不确定性传递这是整个建模中最精妙也最易出错的部分。K不是一个固定常数它的不确定性会直接“传递”到感染人数估算 (N C \times K) 中。在贝叶斯框架下我们得到了K的后验分布那么N的后验分布自然就是 (C \times K) 的分布。在报告结果时必须报告置信区间例如“估算今日感染人数为 1250人95% CI: 800 - 1800”。只报告点估计1250人是严重不完整的会误导决策者。在将K用于SEIR模型预测时更稳妥的做法是进行不确定性传播分析。例如从K的后验分布中抽取1000个样本对每个样本K_i运行一次SEIR模型得到一条预测曲线。最终你会得到1000条预测曲线形成一个“预测带”。这个带子的宽度直观地展示了由于基础参数不确定性导致的预测不确定性。这比只用K的均值跑一次模型要科学得多。4.3 SEIR模型参数化与过拟合陷阱SEIR模型看似只有几个参数(\beta, \sigma, \gamma)但在拟合短期、嘈杂的数据时极易过拟合。特别是初始条件E(0), I(0), R(0)的设定影响巨大。我们的经验固定部分参数潜伏期 (\sigma) 和感染期 (\gamma) 可以根据医学研究设定为固定值如潜伏期5天则 (\sigma 1/5)感染期7天则 (\gamma 1/7)。这减少了待估参数。谨慎设定初始感染者I(0)不要简单设为临床报告数。可以利用反演模型估算的初始感染人数作为I(0)的强先验信息。在贝叶斯框架下可以将I(0)也作为一个待估参数并为其设定一个以反演估算值为中心、有一定宽度的先验分布。使用正则化或先验信息对传播率 (\beta) 施加合理的先验分布如对数正态分布防止其拟合出一些物理意义上不合理的极端值。验证如果数据量允许采用时间序列交叉验证。用前80%的数据拟合模型预测后20%看预测效果。避免使用全部数据拟合后直接宣称模型完美。4.4 措施效果评估的“反事实”框架评估措施效果时最大的挑战是我们无法同时观测到“实施措施”和“不实施措施”两个平行世界的结果。我们只能观测到实施措施后的实际数据。因此我们的评估本质上是构建一个“反事实”场景如果没有这项措施情况会怎样我们的模型校准后的SEIR就是用来生成这个反事实场景的工具。具体步骤用措施实施前的所有数据将模型参数校准到最佳。假设措施未实施即 (\beta) 参数保持措施前的水平不变运行模型预测措施实施后的疫情发展。这是“反事实”曲线。将反事实预测曲线与实际观测到的数据或经过处理后的污水反演感染人数进行对比。效果量化计算两条曲线之间的差异。例如可以计算在措施实施期T内反事实曲线预测的总感染人数与实际估算的总感染人数之差这个差值就是模型估计的“措施避免的感染人数”。也可以比较峰值降低的百分比。必须强调这个评估的可靠性完全依赖于模型本身的质量和我们对“若无措施β保持不变”这一假设的信心。任何模型评估都必须附带对假设和不确定性的充分讨论。5. 论文写作与结果呈现如何将复杂计算转化为清晰故事数学建模竞赛的论文本质上是向评委讲述一个用数学工具解决实际问题的完整故事。写作和图表呈现与模型本身同等重要。5.1 论文结构逻辑我们的论文大致遵循了以下结构但这并非模板而是逻辑的自然流动问题重述与分析用自己的语言精炼地复述问题并立即展示我们拆解出的“反演-分析-评估”三层框架图。让评委一眼看懂你的解题思路。模型准备说明数据预处理步骤平滑、异常值处理、列出所有符号假设。介绍SEIR模型和贝叶斯反演模型的基本原理并重点阐述为何选择它们以及它们在本问题中如何连接。模型求解首先展示贝叶斯反演的结果给出校正因子K的时间序列图带置信区间以及由此反演的感染人数估算图与临床报告数对比。这张图是第一个亮点直观展示了模型校准的合理性和不确定性。然后展示SEIR模型的拟合效果用校准后的感染人数数据拟合SEIR给出拟合曲线图。最后进行措施模拟用不同的β值模拟不同措施绘制未来感染人数和污水浓度的预测对比图。用阴影表示不确定性带。结果分析敏感性分析展示关键参数如先验分布的设定、报告率ρ的变化如何影响最终感染人数估算。这体现了模型的稳健性思考。异常点分析列出检测到的异常点并结合假设的“上下文信息”给出合理解释。措施评估表用表格清晰对比不同措施下的关键指标预测值峰值人数、达峰时间、避免感染人数等。模型评价与推广客观讨论模型的优点如能量化不确定性、整合多源数据、缺点如对先验信息依赖、假设简化以及未来改进方向如引入空间网络模型、考虑病毒变异等。5.2 图表可视化技巧一图胜千言在建模论文中尤其如此。图1数据概览将污水浓度时间序列和临床报告病例数画在同一张图上双Y轴第一时间揭示数据的特征和关联。图2模型原理示意图绘制一个流程图展示“原始数据 - 预处理 - 贝叶斯反演 - SEIR模型 - 措施模拟 - 输出决策支持”的完整逻辑链。图3贝叶斯反演结果这是核心。用一条深色线表示K或感染人数的中位数估计用浅色带状区域表示其95%置信区间。务必在图上标注出识别出的异常点。图4SEIR模型拟合与预测用散点表示历史数据或反演估算值用实线表示模型拟合曲线。预测部分用虚线延伸并用不同的颜色和线型来区分不同措施情景下的预测曲线预测不确定性用半透明色带表示。表1措施效果对比设计一个简洁的表格列是不同措施情景行是关键指标如传播率β降低比例、预测感染峰值、达峰时间、总感染人数估算、避免感染人数等。5.3 代码与可复现性虽然论文正文不展示大量代码但在附录或提交的附件中清晰、注释良好的代码是加分项。我们当时使用了Python主要库包括pandas数据处理、numpy数值计算、matplotlib/seaborn绘图、pymc3贝叶斯建模、scipy微分方程求解。代码文件按功能模块组织并有一个主脚本清晰地按顺序调用各个模块确保评委或任何人能够一键复现我们的主要结果。回顾整个解题过程从最初面对“污水数据”和“校正因子”这些陌生概念的茫然到最终构建出一个能自圆其说、量化不确定性的混合模型最大的收获不是学会了某个特定算法而是锻炼了一种“结构化定义问题”和“用概率思维拥抱不确定性”的能力。在实际的公共卫生决策中数据从来都是不完美的模型也永远是现实的简化。一个优秀的模型其价值不在于做出精准无比的预测而在于能清晰揭示关键驱动因素、量化不同选择的风险与收益并将结论的不确定性坦诚地呈现给决策者。这次数学建模的经历正是对这一理念的一次深刻演练。