基于INLA–SPDE模型与校正气象–水文预测因子的湖沼地区钉螺密度空间预测以指导精准防控

《Science in One Health》：Spatial prediction of snail density in lake–marsh regions to inform precision control: an INLA–SPDE model with corrected meteorological–hydrological predictors

　　

背景：湖北钉螺（Oncomelania hupensis）监测依赖劳动密集型调查，空间覆盖有限。存在–缺失模型无法刻画钉螺密度变异。研究人员利用校正后的气象–水文预测因子对活螺密度进行空间建模，以指导监测。方法：研究人员分析了洞庭湖、鄱阳湖及长江安徽段2023年钉螺监测数据。集成嵌套拉普拉斯近似–随机偏微分方程（INLA–SPDE）模型纳入环境协变量与空间效应，并以偏差信息准则指导网格与协变量选择。研究人员采用空间分块交叉验证及9处洲滩900个框的独立调查，将校正模型与未校正模型及普通克里金、随机森林、随机森林–克里金进行比较。结果：校正INLA–SPDE模型的外部均方根误差（RMSE）为0.064，平均绝对误差（MAE）为0.052只/0.11 m²，较未校正模型分别降低23.81%和23.53%。误差以洞庭湖最低，鄱阳湖次之，安徽最高，但9处洲滩预测值均高于观测值。预测潜在高密度区（>0.300只/0.11 m²）面积分别为洞庭湖445.58 km²、鄱阳湖170.02 km²、安徽82.55 km²。温度、淹水时长与植被效应相似，土壤磷与荒野距离效应存在区域差异。结论：校正输入在所述比较中提高了模型性能。空间估计可指导靶向监测与现场核实，但系统性高估及未经验证的高密度分级限制了直接业务化解读。

研究背景方面，血吸虫病曾流行于中国12个省级行政区，主要分布于长江流域及其以南地区，长期构成公共卫生负担。湖北钉螺（Oncomelania hupensis）作为日本血吸虫（Schistosoma japonicum）在中国唯一中间宿主，维持寄生虫传播循环，始终是监测与防控重点。七十余年综合防控已将流行程度降至极低水平，但人–畜–螺–环境构成的复杂传播系统仍留有传播残余风险。消除阶段残余传播仍是挑战：日本血吸虫属人兽共患病原，野生动物保虫宿主感染可能持续并导致再传播；分子方法较常规监测更敏感，可发现低强度感染；气候变异、水文调控、水利工程及旱涝可能影响钉螺栖息地分布、持续与复现，长江流域证据表明钉螺栖息地对水文与气候变化敏感，因此在栖息地、动物保虫宿主与适宜环境条件重叠处可能出现局部复现，持续监测与靶向控制仍不可或缺。中国钉螺主要栖息于长江中下游湖沼地区，密度呈时空异质性，而常规监测依赖劳动密集型调查、覆盖有限。既往模型多以存在–缺失或栖息地适宜性数据及归一化植被指数（NDVI）、地表温度、高程、坡度、湿度、距水距离等遥感指标绘制潜在栖息地，难以区分低密度栖息地与高密度疫点，凸显密度导向预测的必要性。钉螺种群呈空间簇状、带状及不规则斑块分布，需高分辨率环境数据刻画；常用产品如CHIRPS（约5 km）与ERA5（约10 km）可能遗漏洲滩细尺度变异；卫星反演与空间插值虽覆盖连续，但易受大气散射吸收、大气校正不确定性、地形复杂性、尺度不匹配及变异函数或趋势设定影响；粗糙或易错预测因子可能掩盖栖息地异质性与优先区域，校正气象水文输入并以更高空间分辨率预测密度可改进监测规划。贝叶斯地统计模型将环境固定效应与空间随机效应结合，INLA–SPDE框架可实现高效贝叶斯推断与不确定性量化：INLA无需大量马尔可夫链蒙特卡洛抽样即可近似后验边缘分布，SPDE将连续Matérn高斯随机场表示为三角网格上的稀疏高斯马尔可夫随机场（GMRF），二者结合可整合点观测与栅格协变量，估计空间密度、后验可信区间与超限概率，相关应用包括疾病风险制图、空气污染建模与土壤属性预测。为此，研究人员构建并评价洞庭湖、鄱阳湖及长江安徽段活螺密度INLA–SPDE模型，将钉螺监测记录与校正气象–水文及其他环境协变量结合，评估校正输入的预测贡献并比较替代空间模型，所得空间估计在经现场核实前提下可帮助确定现场调查优先级与监测资源分配。

研究方法方面，研究人员采用2023年洞庭湖、鄱阳湖及长江安徽段国家血吸虫病监测数据，在登记的历史或现症钉螺栖息地开展系统抽样，记录栖息地坐标并将钉螺带回实验室作存活评估与压碎镜检，分析前对行政位置、栖息地名称与类型、植被类型、流行分类、活螺数及调查框数进行标准化；原始数据集分别含441、350、537条记录，剔除位置不完整、调查框数少于50或位置重复记录后最终纳入1,316处栖息地（435、346、535处）。以活螺总数除以系统调查的0.11 m²框数计算平均活螺密度，符合纳入标准但2023年无活螺者密度记为零；为适应零值与右偏分布，密度加1后作λ=0的Box–Cox变换，即z=ln(density+1)，变换密度为响应变量，连续协变量以训练折均值与标准差标准化并应用于验证数据，预测反变换为density=exp(z)?1，MAE与RMSE在原尺度计算；洞庭湖、鄱阳湖、安徽残差Moran's I分别为0.0078、0.0092、0.0023（均P>0.05），指定空间权重下无显著残差空间自相关。环境协变量涵盖直接或间接影响钉螺存活、生长、发育、繁殖与扩散的因素：淹水时长、温度与降水经误差校正并降尺度至100 m，纳入2020（滞后三年）、2021（滞后两年）、2022（滞后一年）年淹水数据以考虑滞后效应；土壤属性包括湿度、pH及全氮、磷、钾；地形与人为输入包括SRTM-DEM（2000）、距河流、道路与荒野距离及夜间灯光强度（WorldPop，100 m）；土壤养分（全磷、钾、氮及pH，0–5 cm）来自时空三极环境大数据平台并由90 m重采样至100 m；土壤湿度来自地理遥感生态网（100 m），以3–8月季节均值代表钉螺活动期；植被协变量包括SPOT/VEGETATION与PROBA-V 100 m产品的NDVI及现场调查植被类型。以ArcGIS 10.2提取各采样点环境协变量值，连续协变量以z分数标准化，偏态或极端离群变量标准化前作对数等变换。INLA–SPDE模型包含截距、环境协变量效应、空间结构随机效应及独立高斯观测误差，空间效应为零均值Matérn协方差高斯随机场，由边际方差σ²与有效范围刻画，经SPDE近似以网格节点权重与稀疏精度矩阵表示，INLA估计固定效应、空间效应与超参数的后验边缘分布；SPDE参数在对数尺度重参数化并设高斯先验，截距与固定效应系数亦设高斯先验。模型选择与比较采用DIC、WAIC与边际对数似然（MLL），初含全部环境协变量，标准化后绝对后验均值系数<0.001者视为可忽略并在剔除降低DIC时移除，直至DIC最小且保留系数绝对后验均值≥0.001；校正与未校正气象–水文协变量模型按同一标准比较，选定后以投影矩阵将预测网格链接至数据栈内网格节点，由INLA生成后验空间效应与钉螺密度预测分布。五种模型在同一数据集、验证记录与100 m输出分辨率下评价：普通克里金（OK）、随机森林（RF）、RF–克里金（RF+OK）及校正或未校正输入的INLA–SPDE；未校正模型将原始气象水文栅格直接重采样至100 m生成预测协变量，校正模型则经环境误差校正与空间降尺度生成同分辨率预测因子，二者分辨率相同而环境数据处理流程不同；响应变换为ln(density+1)，标准化、变量选择与调参仅用训练数据，MAE与RMSE反变换后计算；OK采用指数半变异函数与最近12个观测，RF用500棵树、每次分裂约三分之一预测因子、最小叶节点5、随机种子42，RF+OK先拟合RF再对训练残差克里金；两种INLA–SPDE模型均用高斯似然、Matérn协方差、Mesh 3、相同协变量集与相同先验；固定效应以后验均值与95%可信区间（CrIs，取边缘后验分布2.5与97.5百分位）汇总，基于后验均值与标准差的近似双侧后验尾概率以正态近似计算，95% CrIs排除零视为具后验支持，尾概率不解读为频率学派P值；校正与未校正模型比较评价环境误差校正与空间降尺度相对直接重采样的联合贡献，相对误差降低按（比较模型误差?提出模型误差）/比较模型误差×100%计算。验证方面，采用空间分块交叉验证避免训练与验证记录空间依赖导致性能高估：各区域划分20 km×20 km区块并分配至折，十折交叉验证重复十次（每区域100次训练–验证评价），约90%观测训练、10%验证并尽量平衡折大小，区块尺寸超过估计区域空间范围（洞庭湖5.39 km、鄱阳湖4.82 km、安徽2.39 km）以降低训练与验证数据依赖。外部预测精度在9处独立湖沼洲滩评价（每区域3处，按区域内相对高、中、低高程目的性选取且不共享高程阈值），2023年春季每处按20 m间距布设10×10网格100框，观测平均密度为活螺总数除以全部100框（含零计数框），模型预测在原密度尺度于框坐标提取并在各洲滩内平均，由9对洲滩水平观测–预测均值计算MAE与RMSE。不确定性评估采用研究特定流程：对INLA–SPDE模型作空间分块自助重拟合并以逆RMSE权重聚合预测，100个子模型共享结构与先验，约90%发育数据区块训练、10%留出验证，训练区块有放回重采样，各子模型在共同网格预测密度，集合权重为验证RMSE归一化倒数，加权集合预测为各子模型点预测加权和，子模型点预测间变异以加权自助预测标准差汇总，度量区块重采样与模型重拟合所致点预测加权离散度。密度分级方面，空间预测掩膜至湖南张家界、常德、益阳、岳阳、长沙、株洲，江西九江、南昌、上饶，安徽安庆、池州、铜陵、芜湖、马鞍山；预测密度四舍五入至三位小数后分为零（0）、极低（0.001–0.020）、低（0.021–0.100）、中（0.101–0.300）、高（0.301–0.500）、极高（>0.500只/0.11 m²），预测栖息地面积为掩膜内预测密度>0的网格单元面积之和；分析建模用R 4.2.3，制图用ArcGIS 10.2。

研究结果方面，三角网格选择：各区域构建四个三角网格（Mesh 1–4）并外扩至调查边界外以限制边缘效应，网格密度增加使顶点数增多（如洞庭湖由Mesh 1的7,370增至Mesh 3的70,170）；DIC由Mesh 1至Mesh 3下降，Mesh 4仅边际改善而计算成本显著增加；Mesh 3跨区域约58,273–70,170个顶点，兼顾拟合与效率，在16–32 GB RAM的Intel i7处理器上拟合一个INLA–SPDE模型约5–10 min，Mesh 4超过20 min，故选用Mesh 3；超参数后验中块金方差、边际方差与范围近似正态，Mesh 3与Mesh 4块金后验相似，安徽段空间范围较短约2–3 km，洞庭湖与鄱阳湖分别为5.39 km与4.82 km。环境协变量选择：纳入环境协变量在三个区域均改善拟合（DIC与WAIC更低、MLL更高），基于DIC的选择分别保留洞庭湖13个、鄱阳湖12个、安徽13个协变量，选择后三区域DIC与WAIC分别下降约19%、12%、6%，表明较少协变量下拟合改善；相同协变量集下，误差校正气象–水文数据模型较未校正输入DIC与WAIC更低，三区域分别下降约39%、19%、15%，MLL相应升高，信息准则提示拟合更佳，但不足以单独确立样本外预测改善或将校正效应与输入处理变化分离。空间预测模型比较：外部验证以RMSE与MAE比较，普通克里金误差最高（RMSE=0.103；MAE=0.093只/0.11 m²），其次RF（0.090；0.071）与RF+OK（0.081；0.061），校正输入INLA–SPDE误差最低（RMSE=0.064；MAE=0.052只/0.11 m²），相对OK降低37.86%与44.09%，相对RF+OK降低20.99%与14.75%；相对未校正INLA–SPDE，RMSE由0.084降至0.064、MAE由0.068降至0.052只/0.11 m²，分别降低23.81%与23.53%。钉螺密度空间预测：洞庭湖区域预测钉螺栖息地约2,029.72 km²，环绕湖泊并沿澧、沅、资、湘水系分布，潜在较高密度区集中于西部与西北部尤其澧水与沅水流域，高与极高密度栖息地分别覆盖226.21 km²（11.14%）与219.37 km²（10.81%），呈带状与斑块状，中至极高密度区主要位于石门、临澧、澧县、安乡、桃源、岳阳与华容，中密度383.29 km²（18.88%），极低与低密度主要环绕较高密度斑块，分别736.27 km²（36.27%）与464.58 km²（22.89%）；鄱阳湖区域预测栖息地约992.07 km²，高与极高密度分别94.78 km²（9.55%）与75.24 km²（7.58%），主要围绕鄱阳湖北部、通江河道与修水流域呈带状与斑块状，中至极高密度集中于九江的星子、德安、彭泽、瑞昌、庐山、湖口、永修与都昌及南昌的安义与进贤，中、低、极低密度在髙密度区与五河水系附近更广泛，分别212.16 km²（21.39%）、227.52 km²（22.93%）、382.37 km²（38.54%）；安徽预测湖沼钉螺栖息地约281.62 km²，高与极高密度分别46.82 km²（16.63%）与35.73 km²（12.69%），主要位于安庆的大观、太湖、怀宁、望江，池州的东至、贵池，铜陵的枞阳、铜陵，芜湖的无为及马鞍山的当涂、和县、含山，中密度43.67 km²（15.51%），低与极低密度分别73.23 km²（26.00%）与82.17 km²（29.18%），较低密度区一般环绕较高密度斑块；三区域预测栖息地面积合计超2023年报告钉螺栖息地面积16.29%。验证与不确定性评估：空间分块交叉验证总体RMSE为0.059、MAE为0.048只/0.11 m²，误差以洞庭湖最低、鄱阳湖次之、安徽最高；外部验证用9处独立洲滩900框，总体RMSE为0.064、MAE为0.052只/0.11 m²，洞庭湖最低（RMSE=0.046；MAE=0.037），鄱阳湖次之（0.048；0.043），安徽最高（0.089；0.076），9处洲滩均呈高估；子模型预测加权自助标准差范围0.001–0.060只/0.11 m²，洞庭湖、鄱阳湖及安徽长江主廊道采样密集区变异较低，远离主要湖河水系处较高。模型参数解读：空间随机效应参数上，三区域块金方差、边际方差与空间范围差异明显，提示未解释变异与空间依赖的区域异质性；块金方差洞庭湖最低（均值0.314；95% CrI：0.241–0.385），鄱阳湖居中（0.601；0.522–0.679），安徽最高（0.788；0.714–0.862），提示安徽局部异质性或观测噪声更大；边际方差区域间相近（0.441–0.514）；空间范围洞庭湖最长（5.388 km；95% CrI：3.288–7.484），鄱阳湖次之（4.816 km；3.559–6.071），安徽最短（2.392 km；1.921–2.865），表明安徽钉螺密度空间连续性弱于两湖区域。固定效应参数上，后验估计提示淹水时长、极端温度、NDVI、植被类型与距水道距离与钉螺密度的关联在三区域大体相似，土壤属性与部分距离变量存在区域差异；土壤全磷在洞庭湖（β=?0.047；95% CrI：?0.073，?0.021）与鄱阳湖（β=?0.059；?0.085，?0.033）与密度负相关，安徽则不显著（β=?0.027；?0.073，0.019）；距荒野距离仅在安徽与密度负相关（β=?0.033；?0.051，?0.015），与荒野越近预测密度越高一致；上述格局表明共享水文与植被预测因子之外存在区域变异。

讨论部分，研究人员构建了结合钉螺监测记录、滞后水文变量及包括降尺度误差校正温度与淹水数据在内的多种环境预测因子的INLA–SPDE钉螺密度模型，外部RMSE较未校正输入降低23.81%；尽管9处验证洲滩预测均高于观测，但预测与观测正相关，支持进一步评价基于密度的调查优先级排序；在“同一健康”背景下，这些估计描述血吸虫病监测的环境与中间宿主组分，结合钉螺感染、传染源及人畜接触疫水数据可助力感染风险评估与协同防控，目前主要用途为识别候选监测点供现场调查。模型旨在优先安排调查而非替代现场观察：洞庭湖报告密度梯度为澧水>沅水>资水>湘水，潜在高密度栖息地集中于西北部澧水与西部沅江流域，与既往湖岸、堤垸与洲滩带研究一致；预测亦提示沅江上游，既往调查钉螺密度与分布有限，可能因复杂支流妨碍发现，此差异需现场核查；2,029.72 km²预测栖息地中445.58 km²超过0.300只/0.11 m²，为三区域最大；河网复杂性、堤垸、三峡大坝后径流减少与洲滩扩张可能利于钉螺持续或回升，而生态修复与退田还湖在某些条件下可能增加栖息地适宜性。鄱阳湖预测栖息地992.07 km²，其中170.02 km²超过0.300只/0.11 m²，高密度预测聚集于鄱阳湖北部与通江河道，湖水位与洲滩高程影响栖息地；三峡工程后较早枯水与较长枯季可能使适宜区北移，现场研究亦发现北部感染面积与密度高于南部，故北部及通江河道值得靶向核实。安徽预测钉螺栖息地面积最小但超过0.300只/0.11 m²比例最高（29.31%），主要沿长江及其支流呈带状与斑块状；虽2004–2020年人、牛、螺感染率下降，但钉螺密度与感染面积下降有限，2004–2020年安徽年度灭螺仅覆盖21.40%感染栖息地，水位波动进一步限制江滩效果；操作者依赖的调查检测可能遗漏后回升的低密度种群，密度预测可帮助靶向现场检查并在核实后指导防控。环境关联部分跨区域共享：淹水时长、极端温度与植被效应大体相似，土壤与景观距离变量不同；土壤全磷在两湖与密度负相关而安徽无关，此差异可能反映区域养分累积与微生境结构，缓流湿地磷累积可致富营养化、氧耗竭与水质变差而影响钉螺，安徽冲积洪泛平原流速较强与沉积不稳定可能削弱该关联；仅安徽距荒野越远密度越低，可能因荒野与连续河岸植被及潮湿岸坡重叠而非直接影响密度，两湖荒野边界较不清晰降低该变量解释力；提示可将共享气候、水文与植被预测因子作为跨区域模型起点，而土壤与景观距离变量需本地校准。空间与残差后验参数区域不同：安徽较短估计空间范围提示控制协变量后残差空间相关在较短距离衰减，可能反映带状栖息地、河流整治、圩区农业与景观破碎化；两湖较长范围提示空间依赖覆盖更大距离；空间方差描述空间随机场变异，块金方差描述该场与纳入协变量未解释的残差变异，鄱阳湖与安徽较高块金估计可能提示ln(density+1)尺度上更多非结构变异，包括测量误差或未解析微尺度变异；淹水估计不确定性可能贡献但无法与其他来源区分，估计应结合后验区间解读，更大校准样本与长期监测尤其淹水时长有助于评估改善协变量精度、时间一致性、残差变异与预测。预测钉螺栖息地面积超现场报告感染面积16.29%，此差异可能反映现场检测不完善、抽样代表性或时间尺度不匹配而非单纯模型高估；地形、可达性与劳力限制可妨碍发现隐蔽或低密度栖息地，而空间栅格模型覆盖未调查区域；年度或聚合环境预测因子与季节现场调查代表不同时期，钉螺活动具季节性，覆盖、可检测性与时间尺度差异可能促成面积差异，但其相对贡献未确立。局限性包括：外部验证使用9处洲滩各100框的洲滩水平平均密度而非单个100 m网格结局，最大观测密度低于高密度阈值（>0.300只/0.11 m²），该类别内性能未经验证，预测高密度区需现场确认，需更多覆盖更宽密度范围的洲滩；一年钉螺监测与时间聚合环境数据不能完全刻画季节或年际密度变化；钉螺防控干预记录不可得，未记录干预可能解释部分高估但未证实；淹水时长校正主要依据2023–2025年洞庭湖视频观测校准后转移重建三区域2020–2022年状况，安徽缺乏独立视频观测，转移校正及其对安徽高估的可能贡献无法验证，需更多区域与年份独立观测评估时间稳定性、空间可迁移性及重建淹水不确定性对密度预测的影响。研究结论部分翻译如下：将钉螺监测与校正并降尺度的环境数据整合，在所述模型比较中提高了预测精度。预测的潜在高密度区可帮助确定现场核实与监测点分配的优先级，但未经进一步验证不应视为确诊的高密度栖息地。气候、水文与植被可提供共享的跨区域预测因子，而土壤与景观因素需本地校准。该方法可在数据驱动的钉螺监测与防控中补充而非替代现场调查。