污水流行病学数学建模：从病毒RNA数据反演社区感染与防控评估

1. 项目概述当数学模型遇见公共卫生危机那次参加2022年认证杯数学建模C题第二阶段的经历至今想起来依然觉得非常“硬核”。题目直接把一个前沿的交叉学科——污水流行病学摆在了我们面前要求我们用它来量化评估新冠疫情防控措施的效果。这可不是纸上谈兵它要求我们把实验室里检测到的病毒RNA浓度通过一系列严谨的数学“翻译”变成能反映社区真实感染水平的“情报”。核心关键词就那几个SPSSPRO我们当时的主要分析工具、污水流行病学原理、校正因子。简单说题目就是给了你一堆从城市污水处理厂采样测出的新冠病毒基因片段数据然后问你怎么通过这些数据反推出这个城市大概有多少人感染了更进一步怎么用这个反推出来的数据去评价封控、筛查这些措施到底有没有用这背后是一套非常巧妙的思路。我们都知道感染者的排泄物里会携带病毒。这些排泄物最终都会汇入城市污水管网。那么理论上检测污水处理厂进水口的病毒载量就能窥见整个服务区内人群的感染情况。这比挨个做鼻咽拭子检测要省时省力得多成本也低还能近乎实时地反映疫情动态甚至能发现那些没有症状的“隐形”感染者。但问题就在于从“污水里的病毒浓度”到“人群中的感染人数”中间隔着千山万水。病毒在污水里会被稀释、会降解不同人的排泄规律不同污水处理厂的流量每天在变检测方法也有其灵敏度极限和误差。校正因子就是搭建这座桥梁的核心数学工具它是一系列系数的集合目的就是把所有干扰因素“熨平”让污水数据能真实地映射回人群疫情。所以这个项目本质上是一次数据驱动的政策评估模拟。它适合所有对数学建模、数据分析、公共卫生交叉领域感兴趣的朋友无论是正在备战数模竞赛的学生还是希望了解如何用数据科学解决实际问题的从业者。接下来我就结合我们当时的解题全过程拆解一下这里面的核心思路、技术细节和那些容易踩坑的地方。2. 核心思路拆解从污水数据到防控评估的数学桥梁面对这样一个题目首要任务是构建清晰的逻辑链条。整个建模过程可以分解为三个核心阶段环环相扣缺一不可。2.1 第一阶段数据理解与预处理题目给出的数据通常包括每日污水处理厂进水口的新冠病毒RNA浓度单位可能是 copies/L、每日污水处理厂的处理流量吨/天、服务区的人口估算数量。我们的第一项工作不是急着建模而是“读懂”数据。数据清洗检查是否有缺失值、异常值。比如某天的流量数据为0或极大这显然不符合实际需要根据前后数据插值或结合实际情况处理。病毒浓度数据也可能因为检测限Limit of Detection, LOD而出现“未检出”的情况通常不能简单记为0可以将其处理为LOD值的一半或采用更复杂的统计方法处理。数据可视化这是至关重要的一步。用SPSSPRO或Python的Matplotlib/Seaborn库绘制病毒浓度随时间变化的曲线、流量变化曲线并将它们与已知的官方疫情报告时间线如封控开始日、大规模筛查日放在一起对比。你可能会发现病毒浓度峰值会滞后于报告病例峰值几天这正好体现了污水监测的“前瞻性”。同时观察流量数据是否呈现周期性如工作日与周末的差异这会影响后续的稀释校正。注意预处理阶段对异常值的处理方式会显著影响最终模型结果。我们当时就发现一次暴雨天后流量激增导致那几天的病毒浓度被异常稀释。如果不加处理直接使用会错误地判断那几天疫情大幅缓解。我们的做法是将暴雨日的数据标记出来在后续计算日均负荷时给予较低的权重或直接使用相邻日期的平滑值替代。2.2 第二阶段构建反演模型——校正因子的核心作用这是整个问题的技术核心目标建立方程估计的日新增感染人数 f(测得的污水病毒浓度 一系列校正因子)。最经典的模型是基于质量平衡原理每日病毒基因拷贝排放总量拷贝/天 污水病毒浓度拷贝/升 × 日污水流量升/天但“排放总量”并不直接等于“当日新增感染者排出的总量”。它包含了来自当日新感染者的排放也包含了前几日尚未康复的感染者的持续排放还可能包括环境残留等。因此需要一个感染者病毒排放函数来描述一个感染者从感染到康复期间每日通过排泄物排出的病毒拷贝数。这个函数通常被简化为一个时间分布模型比如一个右偏的分布感染后几天排放量达到峰值然后缓慢下降。由此我们可以建立一个卷积模型 污水中的日病毒总量 每日新增感染人数序列 * 感染者病毒排放函数 误差这里的“*”是卷积操作。我们的目标是已知左边的“污水病毒总量”时间序列去反推右边的“每日新增感染人数”序列。这是一个典型的反问题通常是不适定的需要引入正则化等方法来求解。而校正因子正是为了让我们能更准确地估计上述方程中的各个组成部分而引入的粪便稀释校正因子 (F1)不是所有人的粪便都进入污水系统例如使用化粪池且每人每日粪便产量不同。F1通常基于文献取一个每人每日平均进入污水系统的粪便量如200克/人/天和粪便中病毒浓度估算值来推算。病毒降解校正因子 (F2)病毒RNA在污水环境中会降解。这个因子与污水在管网中的停留时间、温度、pH值有关。通常需要通过实验或文献确定一个半衰期计算从排泄到检测过程中的损失比例。检测方法回收率校正因子 (F3)病毒从污水样本中提取、浓缩、检测的过程会有损失。这个因子通过向真实污水样品中添加已知量的病毒标准品过程控制来测定是一个介于0到1之间的数如0.6表示回收率为60%。人群贡献校正因子 (F4)污水处理厂服务人口可能波动如旅游城市或者服务区内有医院等特殊污染源。需要根据实际情况调整有效贡献人口。最终经过校正的“估计人群日排放病毒总量”计算公式为P_t (C_t * Q_t) / (F1 * F2 * F3 * F4)其中P_t是t日估计的人群排放总量C_t是测得的浓度Q_t是日流量。2.3 第三阶段防控措施效果评估得到“估计的日新增感染人数”时间序列后第二阶段的问题就转化为一个时间序列干预分析问题。常见的评估模型是分段模型或加入虚拟变量的回归模型。例如我们可以以某项严格封控措施的实施日为分界点将时间序列分为“封控前”和“封控后”两段。分别计算两段序列的统计特征基本再生数 (Rt) 的变化利用反演得到的新增感染序列通过EpiEstim等包或自行实现模型如泊松似然法估算封控前后的Rt值。如果封控后Rt显著且持续小于1说明措施有效。感染增长率的比较计算封控前后序列的日增长率可以通过拟合指数增长/衰减模型得到。封控后增长率应显著下降甚至变为负值即衰减。加入虚拟变量的回归模型构建一个以“估计日新增感染数”的对数为因变量以“时间趋势项”和“封控虚拟变量”封控后为1封控前为0为自变量的线性回归模型。封控虚拟变量的系数大小和显著性直接反映了措施对感染数的瞬时影响水平。实操心得评估时一定要考虑措施的“生效延迟期”。封控不是今天下令明天病毒就消失。通常会有3-7天的滞后因为病毒有潜伏期且措施完全落实需要时间。在建模时可以将虚拟变量的生效日期设为措施开始日滞后天数或者使用分布滞后模型这样评估结果会更符合实际观察。3. 基于SPSSPRO的建模实操全流程解析当时我们团队主要使用SPSSPRO进行核心的统计分析它的图形化界面和丰富的统计模型对于快速验证想法非常友好。以下是我们在SPSSPRO中实现的关键步骤分解。3.1 数据准备与探索性分析首先将清洗好的数据日期、浓度C、流量Q、官方报告病例数等导入SPSSPRO的数据面板。计算日病毒负荷使用“转换 - 计算变量”功能生成新变量Daily_Load C * Q。这就是未经校正的每日病毒基因拷贝排放总量。描述性统计与可视化对Daily_Load和官方病例数做描述性统计均值、标准差、最大值、最小值了解数据范围。使用“图形 - 序列图”绘制Daily_Load的时间序列图。添加参考线标记出重要防控措施的时间点。计算Daily_Load与官方病例数的交叉相关函数CCF观察污水信号是否领先于报告病例。在SPSSPRO中可以通过“分析 - 时间序列预测 - 交叉相关”来完成。我们发现污水负荷峰值普遍领先报告病例峰值5-7天这坚定了我们使用污水数据做早期预警的信心。3.2 校正因子应用与感染人数反演这一步是核心计算部分工作需要在SPSSPRO外进行如卷积反演求解但校正计算可以在SPSSPRO内完成。设定校正因子根据文献调研我们假设了一套初始值F1粪便贡献 基于200克/人/天和文献中粪便病毒浓度中位数计算F2降解损失 0.7假设30%降解F3回收率 0.6F4人口校正 1假设服务区人口稳定。计算校正后负荷再次使用“计算变量”功能生成新变量Corrected_Load Daily_Load / (F1 * F2 * F3 * F4)。这个Corrected_Load就是我们估计的每日从感染人群排出的病毒总量。反演每日新增感染人数这是难点。我们采用了相对简化的方法因为竞赛时间有限。方法A基于病毒排放总量和人均排放模型。假设每个感染者在整个感染期内排出的病毒总量平均为G一个基于文献的固定值例如 10^9 拷贝。那么t日的感染人数I_t可以粗略估计为I_t ≈ Corrected_Load_t / G。但这种方法忽略了感染者排放的时间分布误差较大。方法B使用SPSSPRO的“时间序列预测-ARIMA建模”进行反卷积的近似。我们将Corrected_Load序列视为一个输出序列将“每日新增感染人数”序列视为我们想要求的输入序列并假设一个简单的排放分布如一个持续7天的平均排放。这可以转化为一个时间序列的滤波问题。我们使用SPSSPRO的ARIMA模型对Corrected_Load序列进行拟合然后通过分析其残差和白噪声成分结合设定的排放窗口反向推断出感染人数的相对变化趋势再通过一个缩放因子与总人口和总负荷相关将其绝对化。这种方法得到的更多是感染活跃度的相对指数而非绝对人数但对于评估防控措施的效果已经足够。3.3 防控措施效果评估建模在SPSSPRO中我们主要采用回归模型进行评估。数据准备将反演得到的“估计日新增感染指数”记为E_t取对数以稳定方差得到新变量ln_E。创建一个虚拟变量Lockdown在封控措施开始及之后日期赋值为1之前为0。构建分段回归模型进入“分析 - 回归 - 线性回归”。因变量ln_E自变量引入时间t一个从1开始的连续变量代表趋势Lockdown虚拟变量以及可能的时间与虚拟变量的交互项t * Lockdown用来检验封控是否改变了趋势的斜率。模型解读如果Lockdown的系数为显著负值说明封控措施实施后感染水平有一个立即的、显著的下降水平效应。如果t * Lockdown交互项的系数为显著负值说明封控后感染水平的增长趋势斜率变缓了甚至从增长变为下降。我们还可以分别对封控前和封控后两个子数据集单独进行ln_E对时间t的线性回归比较两个回归线的斜率增长率。在SPSSPRO中可以使用“分析 - 回归 - 曲线估计”或“广义线性模型”中的分组功能来实现。踩坑记录直接使用E_t未取对数进行回归分析残差常常不满足方差齐性的假设导致模型不可靠。取对数后不仅满足了线性回归的许多前提假设其系数的解释也更直观可以解释为百分比变化。例如Lockdown系数为-0.5可以解释为封控措施实施后估计感染指数立即下降了约(1 - e^{-0.5}) * 100% ≈ 39%。4. 模型优化、敏感性分析与常见问题排查一个只跑通一次的模型是脆弱的。在竞赛和实际应用中我们必须检验模型的稳健性并知道哪里可能出问题。4.1 校正因子的敏感性分析我们之前假设的校正因子F1, F2, F3, F4取值有相当大的不确定性。因此进行敏感性分析是必须的。方法对每个校正因子设定一个合理的取值范围例如回收率F3在0.4到0.8之间。然后使用蒙特卡洛模拟的方法在这个范围内随机抽取成千上万套因子组合分别计算对应的Corrected_Load和反演感染人数。分析结果观察最终估计的感染人数峰值、累计数以及评估出的防控措施效果如Rt下降幅度随着不同因子组合的变化范围。这能给出一个估计的不确定性区间。例如你可能会得出结论“在给定的参数不确定性下封控措施使Rt值下降了30%至60%”。这比一个孤零零的数字更有说服力。工具SPSSPRO本身进行大规模蒙特卡洛模拟不太方便我们当时是结合Python用numpy进行随机抽样和循环计算来完成这部分再将结果汇总分析。在SPSSPRO中你可以手动调整几组极端值如高估值、低估值来观察关键输出变量的变化进行简单的敏感性测试。4.2 模型诊断与验证残差分析对最终的回归模型如评估防控措施的模型一定要在SPSSPRO中检查残差图。残差应该随机分布在0附近没有明显的趋势或规律。如果残差图呈现漏斗形或弧形说明模型可能遗漏了重要变量或者因变量需要变换。交叉验证为了防止过拟合可以将时间序列数据按时间顺序分成训练集和测试集例如用前80%的数据建模后20%的数据测试。在SPSSPRO的“时间序列预测”模块中可以设置预测期来进行类似验证。看模型在未见过的数据上表现如何。与外部数据对比虽然我们的目标是评估污水数据的效用但最终反演出的感染趋势应该与其他独立指标存在合理的相关性。例如与官方报告的住院人数、重症监护室ICU入住人数、甚至搜索引擎上“发烧”关键词的搜索量进行滞后相关性分析。如果完全背离就需要回头检查模型假设。4.3 常见问题与排查技巧实录在实际操作中我们遇到了不少典型问题这里整理成一个速查表问题现象可能原因排查思路与解决方案反演出的感染人数出现负值或剧烈震荡。1. 反演模型如卷积反解不适定噪声被放大。2. 校正因子设置严重不合理导致校正后负荷数据出现非物理意义的波动。1.使用正则化方法在反演求解时加入平滑性约束如吉洪诺夫正则化惩罚解的大幅波动。2.检查原始数据回溯Daily_Load的计算检查C和Q数据是否有极端异常值需要处理。3.简化模型竞赛中时间有限可退而求其次采用更稳健但分辨率稍低的方法如计算7日移动平均的感染指数而不是追求每日精确值。污水病毒负荷曲线与官方报告病例曲线完全同步没有领先性。1. 数据时间粒度不对如官方病例是按诊断日期而污水数据是采样日期存在录入延迟。2. 该地区病例报告极其迅速滞后时间极短可能性较小。3. 污水监测点位于管网末端混合停留时间过长失去了预警意义。1.对齐时间基准将官方病例数据按症状发作日期如果可获得或样本采集日期重新对齐再与污水负荷日期对比。2.计算交叉相关系统计算污水负荷与官方病例在不同滞后天数下的相关系数找出最大相关性的滞后天数。可能领先优势体现在更早期的、未被报告的传播上。评估模型中封控虚拟变量不显著p值0.05。1. 封控措施本身效果有限或执行不到位。2. 模型遗漏了重要混杂因素如季节性变化、其他同时实施的干预措施如戴口罩宣传、人群免疫力变化等。3. 反演出的感染指数噪声太大淹没了信号。1.加入控制变量在回归模型中引入“星期几”的虚拟变量控制周末效应、节假日变量等。2.改变模型形式尝试使用中断时间序列分析ITSA的更复杂模型或使用面板数据模型如果有多城市数据。3.平滑数据对因变量感染指数先进行移动平均处理再放入模型以降低随机波动的影响。不同校正因子组合下结论差异巨大。模型对某个或某几个校正因子过于敏感说明模型在该参数取值附近不稳定。1.聚焦关键因子通过敏感性分析识别出对输出影响最大的1-2个因子通常是回收率F3和人均病毒排放量F1。2.文献深耕花更多时间查阅最新、最贴近本地情况的文献为关键因子确定一个更窄、更可靠的先验范围。3.报告不确定性在论文中必须明确报告这种不确定性将结论表述为一个范围而不是一个确定值。这是科学性的体现。5. 从竞赛到现实污水流行病学应用的延伸思考做完这个题目我们最大的收获不是学会了几个SPSSPRO的操作而是深刻理解了数学模型如何作为一个“翻译器”和“放大器”将环境中微弱的生物信号转化为可供决策者参考的宏观态势图。在现实世界的应用中这套流程会更加复杂和严谨。例如现实中的校正因子需要通过本地化的研究来确定比如在目标城市开展专门的病毒排泄率研究、污水病毒衰减实验和检测方法学验证。反演模型也会用到更复杂的贝叶斯统计框架将各种不确定性测量误差、参数先验分布直接整合到模型中最终输出的是一个感染人数的概率分布而不仅仅是一个点估计。对于参加数学建模竞赛的同学来说这个题目提供了一个完美的范本如何将一个新兴的、跨学科的现实问题分解为清晰的数据处理、模型构建、参数估计和效果评估的步骤。关键在于理解每个步骤的物理/生物意义而不是机械地套用公式。敢于对模型做出合理的简化以适配竞赛时间但同时必须清晰地说明这些简化的假设及其可能带来的局限性。最后关于工具选择SPSSPRO在数据预处理、统计检验、回归建模和可视化方面非常高效能节省大量编码时间。但对于更复杂的反演计算、蒙特卡洛模拟或自定义算法结合PythonPandas, NumPy, SciPy或R语言是必要的。在竞赛中灵活运用多种工具发挥各自长处是取得好成绩的关键。我个人在后续的研究中就更倾向于用Python完成全部流程因为它提供了从数据清洗到复杂建模再到生成报告的全链条控制能力但对于快速原型验证和统计分析SPSSPRO的交互式界面依然有其不可替代的优势。