时空贝叶斯非参数聚类

Bayesian nonparametric clustering for spatio-temporal data, with an application to air pollution

时空数据的贝叶斯非参数聚类及其在空气污染中的应用

https://arxiv.org/pdf/2505.24694v1

• 概要 : 本文开发了一种贝叶斯时空模型，用于分析意大利北部有害臭氧空气污染水平。模型结合了时间序列和空间依赖性，利用贝叶斯方法进行参数估计和不确定性量化，为环境政策制定提供了数据支持。
• 推荐理由 : 对于研究环境统计和时空建模的研究者，这篇论文展示了贝叶斯方法在环境科学中的实际应用，具有较强的政策相关性。

摘要 空气污染是一项重大的全球健康危害，细颗粒物（PM10）与严重的呼吸系统和心血管疾病密切相关。因此，分析和聚类时空空气质量数据对于理解污染动态并指导政策干预至关重要。本文综述了贝叶斯非参数聚类方法，特别关注其在环境科学中普遍存在的时空数据上的应用。我们首先介绍了针对点参照时空数据的关键建模方法，强调其在捕捉复杂空间和时间依赖性方面的灵活性。随后，我们回顾了贝叶斯聚类的最新进展，重点关注将空间结构融入聚类过程的空间乘积分区模型。我们以意大利北部的PM10监测数据为例，展示了所提出方法识别有意义污染模式的能力。本综述突显了贝叶斯非参数方法在环境风险评估中的潜力，并为公共健康与环境科学领域时空聚类的未来研究方向提供了见解。 关键词——空气质量；环境风险评估；空间乘积分区模型；时间序列分析

1 引言

空气污染是当今最紧迫的环境挑战之一，直接和间接地威胁着生态系统与公共健康。事实上，空气污染不仅会直接危害吸入污染空气人群的健康，同时也可被视为更广泛的气候变化问题的成因和后果。鉴于这些日益增长的风险，持续的空气质量监测以及对所得数据进行严谨的统计分析，对于识别最脆弱的区域并制定有效的缓解策略至关重要。这要求统计模型能够处理不确定性、捕捉复杂关系，并整合多样化的数据。

空气质量测量本质上具有时空特性，即在空间和时间两个维度上被记录。从空间维度来看，空间数据通常分为三类：点参照（地统计）数据，指在特定位置采集的测量值（例如，不同监测站的空气质量水平）；面域（格网）数据，指在特定空间区域上聚合的测量值（例如，按行政区划统计的土地利用率）；以及点模式数据，指事件在空间上的发生位置（例如，地震或野生动物出没的地点）。关于空间数据的全面概述，可参见 Banerjee 等（2014）和 Cressie 与 Wikle（2015）。时间维度也存在类似的区分，时间可被视为连续的（在 R+ 或其子区间上）或离散的（例如，每小时、每日等）。在后一种情况下，需明确每个测量值代表的是某一时间间隔内的计数或平均值（类似于空间数据中的面域测量），还是某一时刻的单点观测值。

在过去三十年中，同时具有空间和时间索引的数据集日益增多，推动了能够同时考虑空间和时间结构以及二者可能交互作用的随机模型的发展。例如，在空气污染研究中，我们不仅关注污染物的空间分布，还关注这种分布如何随时间变化。贝叶斯时空模型在这一背景下发挥着关键作用，为量化不确定性以及整合空间和时间依赖性提供了连贯的概率框架。贝叶斯模型在聚类分析方面也特别有效，而聚类分析是时空数据建模（更广泛地说，环境数据建模）中的一个著名难题。其应用范围广泛，从早期识别雷区或地震断层的研究（Dasgupta 和 Raftery, 1998），到近年来关于空气质量监测（Cheam 等, 2017）、沙蚤觅食行为（Ranalli 和 Maruotti, 2020）、沿海降雨模式（Paton 和 McNicholas, 2020）以及鱼类生物多样性（Vanhatalo 等, 2021）的研究，仅举数例。在许多应用中，研究重点在于观测值的聚类分配，以及每个聚类内部的模式。

更广义地讲，聚类分析是指将样本推断划分为若干子集（称为聚类），使得同一聚类内的统计单元尽可能相似，而与其他聚类中的单元尽可能不同。换句话说，聚类旨在实现聚类内部的凝聚性和聚类之间的分离性。在基本的聚类形式中，上述划分在数学意义上是一个划分，即一组互不重叠的子集，其并集构成整个样本。也就是说，每个单元仅被分配到一个聚类中。然而，更复杂的聚类方法允许例如同一单元属于多个聚类（重叠聚类，Palla 等, 2005），和/或为每个单元赋予一个隶属度（模糊聚类，Bezdek 等, 1984）或隶属概率（贝叶斯聚类，见第3.1节）。

聚类方法大致可分为启发式（算法型）和基于模型（概率型）两类。一种著名的启发式方法是k均值算法，它通过迭代地将每个数据点分配给最近均值的聚类，将观测值划分为k个聚类。关于k均值及其他基于距离的聚类技术的全面讨论，参见 Sarang（2023）。启发式策略的优点是计算速度快，能够处理大量数据或高维数据集。但其缺点是仅返回上述划分的一个点估计，无法量化每个单元聚类归属的不确定性。相比之下，基于模型的方法（Fraley 和 Raftery, 2002）通过在统计模型中引入聚类结构，能够量化推断划分的不确定性。这还使得我们能够回答诸如需要多少个聚类、如何检测和处理异常值、协变量如何影响聚类等问题。

概率聚类的一种自然方法是基于混合模型，其中观测值由若干密度函数的凸组合（即权重和为1的线性组合）建模，这些密度函数称为混合成分。这些密度通常属于同一参数族，称为混合核。关于混合模型的更多细节，见第3节。然后，可根据观测数据，通过频率学派或贝叶斯方法，对成分权重及各成分密度参数进行推断。关于基于模型聚类的全面综述，参见 Bouveyron 等（2019）；关于贝叶斯方法的深入讨论，参见 Wade（2023）和 Grazian（2023）。

混合模型中的一个关键问题是确定成分的数量。这与确定聚类数量的问题相对应（下文将详细讨论聚类与成分之间的关系）。需注意，尽管对于大小为n的有限样本，最多只需n个组，但在总体层面上，原则上可以允许无限多个聚类。这使得统计聚类问题本质上是非参数的。类似地，在混合模型中，我们可以允许成分数量为无限，或为有限但随机，从而采用非参数方法。尽管无限混合模型在贝叶斯非参数领域已有长期研究，但直到最近，成分数量为随机的情况才被正式纳入非参数框架（Argiento 和 De Iorio, 2022）。在这两种情况下，混合模型都会在划分上诱导一个先验分布，使得在大样本极限下，聚类数量没有上界。此类先验中最著名的例子包括所谓的“有限混合的混合”（参见 Miller 和 Harrison, 2018），它假设成分数量有限但随机；以及狄利克雷过程（DP）混合模型（Ferguson, 1973；Lo, 1984），它允许无限多个成分。文献中还提出了许多其他灵活的推广形式；关于全面综述，可参考 Müller 等（2015）和 Frühwirth-Schnatter 等（2019）。

在对空间点参照数据进行聚类时，一种突出的贝叶斯非参数方法是空间狄利克雷过程（Spatial Dirichlet Process，Gelfand 等，2005），它表现为一系列以概率加权的随机空间曲面的集合。Duan 等（2007）通过引入空间依赖的混合权重对该模型进行了推广。然而，这两种空间狄利克雷过程都需要在每个空间位置上有多次观测。作为替代方案，空间断棒过程（spatial stick-breaking process，Reich 和 Fuentes，2007；Rodríguez 和 Dunson，2011）允许聚类分配概率随空间变化，而无需重复观测。此外，Nguyen 和 Gelfand（2011）提出了狄利克雷标记过程（Dirichlet labeling process），作为一种灵活的空间聚类模型。尽管这些方法非常有吸引力，但当存在可用于辅助聚类识别的协变量时，其性能可能受到限制。

对于面域数据的聚类，Fernández 和 Green（2002）的开创性工作提出了一种具有随机成分数量的有限混合模型，并通过马尔可夫随机场对具有空间依赖性的混合权重进行建模。在时间序列聚类方面，Nieto-Barajas 和 Contreras-Cristán（2014）在状态空间框架内提出了一种贝叶斯非参数混合模型，而 Frühwirth-Schnatter 和 Kaufmann（2008）则提出了一种动态回归模型的有限混合模型。

基于上述文献，本文展示了贝叶斯非参数模型在聚类时空数据方面的潜力，这类数据在环境科学中极为普遍，包括但不限于空气质量监测。具体而言，此类方法使我们能够：量化对估计划分及聚类特异性参数的不确定性；避免预先指定聚类数量；并整合地理坐标信息——这在环境数据集中通常也是可用的。特别是，引入协变量有助于生成更简洁、空间上更连贯（因而更具可解释性）的聚类。

本文其余部分结构如下：第2节和第3节概述了用于聚类分析的混合模型，重点关注时空数据和贝叶斯非参数方法。特别地，在第4节中，我们详细讨论了乘积分区模型（Product Partition Models, PPM，Hartigan, 1990）及其在时空场景下的扩展。随后，我们在真实世界空气质量数据上展示了该方法的应用（第5节），旨在突出此类技术在所有涉及时空数据的应用中的广泛适用性。最后，第6节提出总结性评述及未来研究方向，强调贝叶斯非参数方法在环境数据科学中的巨大潜力。

2 时空数据模型

连续时空数据的一种常用模型是：

3 用于聚类分析的混合模型

在混合模型中，每个观测值被认为来自 M ≤ ∞ 个混合成分中的某一个。即，该模型表示为：

混合核的选择同样至关重要——尽管常常被忽视——它取决于数据的性质以及分析的目标。对于连续型数据，虽然高斯核是最流行的选择，但可通过使用多元偏正态（skew-Normal）或偏t分布（Frühwirth-Schnatter and Pyne, 2010; Lee and McLachlan, 2014）、平移后的非对称拉普拉斯分布（Franczak et al., 2013）或正态逆高斯分布（O'Hagan et al., 2016）来实现偏态或对异常值的鲁棒性。对于分类数据，可以采用潜在类别模型，这类模型依赖于伯努利或多项分布的混合（Goodman, 1974; Argiento et al., 2024）。最后，对于计数数据，泊松混合（Karlis and Xekalaki, 2005; Krnjajić et al., 2008）、负二项分布（Liu et al., 2024）或零膨胀泊松、零膨胀负二项分布对于处理稀疏计数数据非常有用（Wu and Luo, 2022）。

3.1 贝叶斯聚类

将方程 (3) 中的抽样模型置于贝叶斯框架下，需要对 P设置一个先验分布，即对一个几乎必然离散的概率测度设置先验。这一任务通常需要高级的概率工具，但由于引入了分配向量 c，模型设定可以被方便地简化如下：

请注意，EPPF 仅依赖于各聚类的大小

。此外，注意：对于一个合适的 EPPF，存在一个几乎必然离散的随机概率测度 P，称为物种采样模型（Species Sampling Model，Pitman, 1996），使得该模型可以表示为公式 (4) 的形式。

3.2 后验推断

当成分数量固定且有限时，使用马尔可夫链蒙特卡洛（MCMC）方法从聚类分配和聚类特定参数的后验分布中采样相对简单。在这些情况下，算法通常可以简化为具有共轭完全条件分布的吉布斯抽样器。例如，当基测度 P0与混合核共轭，且权重来自狄利克雷分布时，就会出现这种情况。

相比之下，当成分数量为无限或有限但随机（如有限混合的混合模型）时，后验采样变得更加具有挑战性。特别是后者情况需要一种跨维度的算法，以适应不同维度的参数空间之间的跳跃（根据混合成分的数量）。早期的重要突破由 Green (1995) 和 Richardson 与 Green (1997) 实现，他们提出了用于单变量高斯混合模型的可逆跳跃 MCMC 算法。该算法非常流行，但它要求设计良好的可逆跳跃移动，这在实际应用中——尤其是在高维参数空间中——带来了挑战。

另一方面，无限混合模型的后验采样器可分为两类：(i) 条件算法，能够对混合参数和聚类结构进行完整的贝叶斯推断（Papaspiliopoulos 和 Roberts, 2008；Kalli 等, 2011）；以及 (ii) 边缘算法，通过将混合参数积分掉来简化计算，从而专注于聚类本身（MacEachern 和 Müller, 1998）。早期进展包括引入 Pólya urn 吉布斯抽样器（Escobar, 1994；MacEachern, 1994；Escobar 和 West, 1995；MacEachern, 1998），随后 Ishwaran 和 James (2001) 引入了狄利克雷过程的断棒表示，使得可以通过吉布斯抽样实现后验推断。

最近，利用有限混合模型与无限混合模型之间的联系，为贝叶斯非参数模型开发的算法已被推广到有限混合模型。例如，Miller 和 Harrison (2018) 提出的中国餐馆过程抽样器、Frühwirth-Schnatter 等 (2021) 开发的望远镜式抽样，以及 Argiento 和 De Iorio (2022) 提出的两种增强型吉布斯抽样器。这些进展提高了有限混合模型的计算效率和实际应用能力。

4 空间乘积分区模型

4.1 sPPM 下的后验采样

5 空气质量数据的分析

根据世界卫生组织（World Health Organization, 2022）的数据，每年约有700万人的过早死亡可归因于室内和室外空气污染的暴露。颗粒物（PM）是衡量空气污染的标准代理指标。与其它污染物不同，PM并非单一化学物质，而是一种由固体和液体微粒组成的复杂混合物，这些微粒能在大气中长期悬浮，从而经历扩散和传输过程。这些颗粒物来源于自然因素，如土壤侵蚀、火山活动和花粉传播，也来源于人类活动，包括工业生产、取暖和机动车交通。暴露于PM会对生物体造成严重的健康风险，并对环境产生显著影响，加剧气候变化，同时污染土壤和水源。因此，持续且广泛地监测PM水平对于防止达到被认为有害健康的浓度至关重要。PM水平通常通过分布在目标区域内的监测站，在特定时间间隔内进行测量。因此，PM数据具有固有的时空特性。在本节中，我们将说明前几节所介绍的贝叶斯非参数聚类方法如何应用于这类数据。

在贝叶斯文献中，已有多种方法被提出用于分析PM数据。例如，Sahu等（2006）提出了一种包含两个随机效应的贝叶斯模型，分别用于农村和城市地区的PM水平；Cocchi等（2007）则提出了一个贝叶斯时空模型，用于预测每日PM浓度。Cameletti等（2011）基于拟合优度和可解释性比较了不同的贝叶斯模型用于PM数据分析。Hamm等（2015）提出了一种具有空间变化系数的模型，用于在大范围空间区域内绘制PM水平分布图。最后，尽管这绝非详尽无遗的列表，Datta等（2016）引入了一种基于时空高斯过程近似的贝叶斯模型，用于分析和预测大尺度PM数据。

5.1 数据描述

在本研究中，我们关注的数据集包含2019年1月1日至12月31日期间，意大利北部162个监测站每日测量的PM10浓度水平（即直径小于10μm的颗粒物）。该数据集由欧洲环境署（EEA）收集，并可免费获取于 https://eeadmz1-downloads-webapp.azurewebsites.net/。图1展示了这些监测站中有多少个超过了50μg/m³的浓度限值，以及超过该限值的天数。需要注意的是，根据EEA的规定（European Commission, 2008），每个监测站每年不应超过上述浓度限值超过35天。然而，在2019年，意大利北部162个监测站中有77个站点在超过35天的时间内超出了该限值。在一些受影响特别严重的区域，某些站点甚至记录到超过80天的超标情况。这凸显了意大利北部空气质量的严峻性，亟需进行细致的监测。该地区的空气污染——尤其是波河谷地区——是一个有充分文献记载的现象，其成因包括高密度的工业活动、城市区域以及位于被山脉环绕的平坦盆地中的密集畜牧业，这些因素阻碍了空气的有效交换（例如，Larsen等，2012）。

图2展示了2019年所有监测站点平均PM10浓度的时间序列，同时附有四分位距和90%置信带，以说明整体污染水平及其季节性模式。具体而言，PM10浓度在较冷月份达到峰值，数值频繁超过50 μg/m³，偶尔接近甚至超过100 μg/m³。这种季节性模式与冬季常见的已知因素相符，如取暖排放增加、大气扩散能力下降以及静稳空气条件（Pietrogrande等，2022）。相比之下，夏季则表现出较低的PM10浓度，通常低于50 μg/m³。图2中的置信带也显示出各监测站之间存在显著的变异性。

相反，图3展示了每个监测站点PM10浓度的均值和标准差。具体而言，图(a)显示最高均值（通常超过30 μg/m³）出现在都灵、米兰、布雷西亚、维罗纳和威尼斯等城市和工业区，而较低的均值浓度（低于20 μg/m³）则出现在北部和沿海地区。此外值得注意的是，南部波河谷地区的浓度处于中等水平，通常在20至30 μg/m³之间。通过对比图(a)和图(b)，后者展示了每个站点的标准差，我们可以观察到：平均浓度较高的地区，特别是城市和工业区，往往也经历更大的时间波动。

5.2 贝叶斯非参数聚类

在聚类过程中引入相似性函数（即在sPPM框架下）所获得的结果如图6所示。图(a)中的共聚类矩阵反映出，与未包含空间协变量的PPM相比（即与图5的(a)部分相比），划分结构的不确定性进一步降低：聚类数量更少，且站点分配到聚类的不确定性也更低。图(b)展示了根据估计聚类结果着色的监测站分布，说明了引入相似性函数后，使得原本因靠近而被分配到不同聚类的监测站更倾向于被归入同一聚类。由此生成的地图揭示出三个明显不同的区域：波河谷中心区域（红色聚类）、阿尔卑斯山和亚平宁山脉周边地区（橙色聚类），以及沿海和山区（黄色至绿色聚类）——这些区域的持续性和变异性依次递减（仅有一个例外）。

最后，图7中的河流图进一步展示了基于PPM（左侧）和sPPM（右侧）所得到的VI损失点估计划分之间的差异。颜色代表在PPM下形成的聚类，而框内标注的数字是通过以估计划分条件运行算法1计算得到的聚类特定参数。再次注意到，利用空间信息可得到更少的聚类（5个而非9个），从而形成更简洁、更具可解释性的聚类结构。

为了进一步说明我们聚类方法的结果，在附录（第A4节）中，我们展示了从PPM和sPPM分别得到的每个聚类的时间序列均值，以及四分位距和90%置信带。这些图表明，具有较高持续性和变异性特征的聚类对应于PM10污染水平更严重的区域，而持续性和变异性较低的聚类则通常记录到较为适中的PM10浓度。

6 总结与未来研究方向

本文综述了用于聚类时空环境数据的贝叶斯非参数方法。我们首先回顾了时空数据的建模方法，然后讨论了贝叶斯非参数聚类方法，并展示了如何纳入个体协变量（例如地理坐标）。随后，我们将该方法专门应用于环境科学中普遍存在的时空数据。最后，我们以意大利北部波河谷地区的PM10监测数据为例，展示了所提出方法在识别有意义空间模式方面的有效性。

具体而言，在本文的应用研究中，我们采用了一个多层次的分层线性混合模型，引入一阶自回归过程以刻画时间效应，并将具有相同自回归系数和时间方差的时间序列聚类在一起。随后，我们通过sPPM框架引入地理坐标，对非参数模型进行了扩展，以在聚类过程中诱导空间同质性。

未来关于时空环境数据的贝叶斯非参数聚类研究，可聚焦于探索sPPM框架内不同的凝聚函数和相似性函数。此外，将sPPM方法扩展至面域数据（areal data）将是一项具有重要实际意义的原创性贡献。另一个重要的研究方向可能是优化超参数的选择，以提升模型性能。最后，现有文献在这些模型的预测能力方面尚存空白。尽管这些模型已展现出强大的推断性能，但关于其预测表现的研究仍十分有限，而预测能力对于理解环境变化至关重要。

原文链接：https://arxiv.org/pdf/2505.24694v1