一种用于标记非均匀泊松点过程强度估计的非参数惩罚方法

《Spatial Statistics》：A nonparametric penalized approach for intensity estimation of marked inhomogeneous Poisson point processes

　　摘要：带标记的点过程为建模发生在空间中的事件及其相关属性提供了自然框架，这些属性被称为标记。该技术应用于多个领域，如地震学、生态学和野火分析，其中的标记往往是连续的，并可能呈现重尾分布特征。一个重要的研究量是联合强度函数，它描述了事件在空间和标记两个维度上的发生情况。现有的估计方法通常依赖于严格的可分性假设，在处理复杂领域时面临困难，且在标记分布右偏的情况下性能较差。我们提出了一种完全非参数的、带有惩罚项的似然估计器，用于连续标记点过程的联合强度函数。该方法避免了可分性假设，能够自然处理非凸和弯曲的空间域，并引入了一个与标记相关的指数型平滑参数，尤其在标记分布为截断重尾型时能提升估计精度。理论分析证明了该估计器的唯一性及渐近收敛性。为便于实际应用，我们采用了基于张量积离散化的数值方法，结合有限元和B样条技术，以及高效的优化算法和基于数据的平滑参数选择方法。最后，通过模拟研究将该方法的性能与现有方法进行对比，并通过加利福尼亚州的野火案例以及全球地震事件的真实数据应用，展示了其实际效用。

引言：带标记的点过程是一组随机点集合，每个点记录了事件的位置及相关属性，即标记。这类模型在多个领域都有应用。例如在地震学中，空间过程由地震震中构成，每次地震的震级即为标记（Ogata, 1988; Schoenberg等, 2002）。在生态学中，则由树木位置表示，树木直径作为标记（Stoyan和Penttinen, 2000; Comas和Mateu, 2007; M?ller等, 2016）。在野火研究中，空间过程由起火点构成，燃烧面积作为标记（Díaz-Avalos等, 2016; Koh等, 2023）。图1展示了一个带标记点过程的示例，显示了2013年5月至2025年3月间加利福尼亚州的野火发生情况及其对应的燃烧面积。野火是一种无法控制的植被火灾，会对生态系统和环境造成威胁，了解起火点及其燃烧面积的空间分布对于有效预防和减轻灾害至关重要。这些数据可从加利福尼亚州林业与消防局获取（CAL FIRE, 2025）。如图左面板所示，起火点并非均匀分布：野火多集中在该州南部、沿海地区以及内陆中央山谷地带。右面板以对数尺度展示了标记的实证分布，可见其具有明显的右偏特征。大多数事件对应的是较小或中等规模的火灾，而少数极端事件却占据了绝大部分的燃烧面积。图2展示了另一个示例，显示了2024年8月至2025年2月间全球范围内震级大于2.5的地震，不同颜色代表不同的震级。在这些情况下，事件发生在二维弯曲的地球表面。大多数地震集中在与构造断层线相对应的狭长区域内，因为这些地方的地震活动最为频繁。受这些应用启发，我们将研究重点放在具有空间不均匀性及定量标记的点模式上，这些标记可能呈截断重尾分布，可视为连续变量进行处理。

从形式上讲，带标记的非均匀泊松过程的实现是一组事件{(si,mi)}i=1n，其中n为事件总数，si∈S?Rd（d=2,3）表示位于可能具有不规则边界或弯曲几何形状的空间区域S内的位置，mi∈M?R+则为相应的定量标记（Cressie, 2015第8.7章）。该过程的计数测度N(S×M)表示位置-标记对落在S×M内的观测数量，其期望值定义了强度测度?Λ(S×M)?E[N(S×M)]。当强度测度关于S和M上的勒贝格测度乘积是绝对连续的时候，就存在一个强度函数λ(s,m)。该函数表示在位置s∈S附近的小区域内，带有标记m∈M的点数期望值（参见例如Chiu等, 2013第4.2章）。将λ(s,m)在区域S×M上积分即可得到Λ(S×M)，即该区域内事件发生的期望数量。本研究的主要目标是开发一种非参数框架，用以精确估计描述事件在空间和标记维度上共同分布情况的λ(s,m)。

在现有文献中，估计强度函数λ(s,m)的方法多种多样，但往往依赖于一些简化假设。其中最严格的是一阶可分性假设，根据该假设，强度可分解为独立的空间成分和标记成分，即λ(s,m)=λg(s)λM(m)（Cressie, 2015第8.6.3章）。此处，λg(s)和λM(m)分别表示边缘空间强度和标记强度。在这种情况下，仅通过空间位置和标记即可独立估计这两个分量。这一假设已被应用于一些领域，比如在林业统计中，点代表树木位置，标记为胸高处的树干直径（Stoyan和Penttinen, 2000）；在地震学中，则假设地震震级与震中位置相互独立（Ogata, 1998）。

相比之下，许多估计方法依赖于较弱的条件可分性假设，根据该假设，联合强度可表示为λ(s,m)=λg(s)fM(m∣s)，其中fM(m∣s)表示给定位置s时标记的条件分布。在这种框架下，地面点模式{si}i=1n被建模为具有强度λg(s)的非均匀泊松过程，而标记则在给定位置的情况下被视为具有密度fM(m∣s)的独立随机变量（Daley等, 2003定义6.4.I）。虽然地面强度λg(s)可以使用标准的空间强度估计技术来估算，比如核平滑法（Diggle, 1985）、对数高斯Cox过程建模法（M?ller等, 1998），或是基于惩罚样条和偏微分方程的正则化方法（Gu, 1993; Ferraccioli等, 2021），但对条件标记分布fM(m∣s)的直接估计方法则发展得较少。现有方法往往依赖于较强的分布假设，这些假设通常是在地质统计框架下提出的。例如，阈值克里金法已被视为处理带标记点过程的一种通用方法（Mase, 1996），而标准克里金法则被用于林业领域，以插值计算诸如树干尺寸之类的连续标记值（Biondi等, 1994）。这些方法具有很强的适用场景限制，仅适用于地质统计标记情况，即当m(s)=U(s)，且U(s)为与空间点过程{si}i=1n无关的随机场时（W?lder和Stoyan, 1996; Stoyan和W?lder, 2000）。

即便在条件可分性假设下，这两个分量也并非总能被独立估计。例如，虽然fM(m∣s)中存在空间依赖性是可行的（Ogata和Katsura, 1988），但Myllym?ki和Penttinen（2009）提出的那种基于强度依赖性的标记形式fM(m∣λg(s))则是不被允许的，因为它会阻碍分离空间效应和标记效应。尽管已有方法可用于检验一阶可分性（Schoenberg, 2004; Chang和Schoenberg, 2011）以及标记与位置的依赖关系（Schlather等, 2004; Guan和Afshartous, 2007; Zhang, 2017），但据我们所知，目前还不存在正式的标准来判断λg(s)和fM(m∣s)是否可以独立估计。此外，条件可分性假设还隐含着要求空间强度结构在标记值变化时保持不变。正如Goulard等人（1996）所指出的，实际情况可能并非如此，因为点的空间分布本身就可能受到标记结构的影响。这种现象有时被称为“纠缠效应”（Renshaw等, 2007），它会破坏条件可分性，使得仅依据位置条件来分析标记变得不可行。正因如此，通过在乘积空间S×M上进行分析来估计联合强度函数λ(s,m)，可以避免那些可能简化潜在生成机制并降低估计精度的可分性假设。

目前用于估计联合强度函数的方法主要以非参数方法为主。核估计器已被直接应用于具有连续标记的带标记点过程，尤其是在林业研究中（Penttinen等, 1992）。然而，核密度估计存在诸多局限性。已知它在处理具有大量尾部质量的严重偏态标记分布时表现不佳，而这类分布在许多涉及连续标记的应用中十分常见。例如，按体积或蔓延距离测量的山体滑坡、按水位峰值或流量测量的洪水、按燃烧面积测量的野火，以及按震级测量的地震等（Katz等, 2002; Malamud和Turcotte, 1999）。此外，核密度估计还存在难以修正的边缘效应，尤其是在非凸空间域中（Baddeley等, 2016第6.5章）。最后，随着维度增加，其性能会迅速下降，因为维度灾难会导致需要使用更宽的带宽，进而造成过度平滑现象（Scott, 2015第7章）。在贝叶斯框架下也提出了核密度估计的替代方法。特别是，基于狄利克雷过程混合模型的非参数方法被用于建模联合标记-位置分布，从而能够同时对强度和条件标记分布进行推断（Taddy和Kottas, 2012）。不过，这一框架并未配备公开可用的代码，限制了其可重复性和实际应用前景。此外，还在分层贝叶斯框架下开发了用于具有协变量的带标记非均匀泊松过程的有限高斯混合模型（Christou Micheas, 2014），但目前的相关实现并不支持连续标记（Micheas和Chen, 2018）。

除上述方法之外，更广泛的空间点过程研究领域中也发展出了若干方法论框架，但其中只有少数被直接用于估计具有连续标记的点过程的强度。这些方法包括基于广义加性模型和张量积平滑技术的半参数方法（Baddeley等, 2012; Wood, 2017）、用于泊松过程的参数化对数线性模型（M?ller和Waagepetersen, 2003），以及如对数高斯Cox过程之类的潜场建模方法（M?ller等, 1998; Diggle等, 2013），这类方法通常通过集成嵌套拉普拉斯近似法来实现（Rue等, 2009; Simpson等, 2016）。这些方法大多是为无标记过程或具有分类标记的过程开发的（Illian等, 2012; Illian等, 2013; Diggle等, 2013; Gomez-Rubio等, 2015）。在少数涉及连续标记的情况下，人们通常通过假设各变量相互独立来避免对完整联合强度λ(s,m)的估计（参见例如Díaz-Avalos等, 2016）。将这些框架扩展为直接针对乘积空间S×M上的完整联合分布，仍是一个有待研究的课题。

此外，准确估计完整的强度λ(s,m)以及由此衍生的地面强度λg(s)，为量化该过程的联合空间-标记结构奠定了基础。由于二阶汇总统计量是通过强度重新加权来定义的，因此强度估计的误差会直接影响到这些统计量；反之，可靠的强度估计则有助于更准确地描述聚类和相互作用模式。尤其是像非均匀K函数和J函数这样的量，它们依赖于强度重新加权，被广泛用于检测空间聚类或抑制现象（Cronie和Van Lieshout, 2016）。在带标记的情境下，Ghorbani等人（2021）提出的标记加权非均匀K函数，以及D’Angelo等人（2024）提出的其局部扩展形式，进一步拓展了这些工具的功能，它们能够考虑标记与空间位置之间的相互作用。此外，针对带标记的时空过程（Iftimi等, 2019）以及网络上的带标记点过程（Eckardt和Moradi, 2024），也有了进一步的推广。相关的频域分析框架，比如用于多变量带标记点过程的图形建模方法（Eckardt和Mateu, 2019），同样依赖于准确的强度估计，以便量化子模式之间的条件依赖关系。总体而言，这些发展凸显了基于强度的方法的广泛应用价值，同时也强调了准确强度估计方法的重要性。另外，对于在如地球表面这样的弯曲空间域上估计带标记点过程的强度，目前仍几乎未被探索，据我们所知，尚不存在适用于此类数据场景的方法。

正是基于这些考虑，我们有必要开发新的方法，能够灵活地估计完整的联合强度函数，同时兼顾重尾标记以及形状不规则或弯曲的空间域。我们提出的方法正好解决了这些问题。它提供了一种完全非参数的、基于似然的λ(s,m)估计器，无需遵循严格的可分性假设，而且能够自然处理在非凸和弯曲空间域中观察到的点模式。该方法在那些不可分离性不成立且空间区域的几何结构可能影响研究现象的实际应用场景中具有广泛的适用性。此外，还引入了指数型的标记依赖平滑参数，以便更好地处理具有强右偏特征的、定义在有限支撑集上的标记分布。这项工作为起源于单变量场景（Good和Gaskins，1971年；Silverman，1982年）的研究方向做出了贡献，后来该研究被扩展到多变量点模式（Scott等人，1978年；Scott等人，1980年）以及发生在空间区域内的点模式（Gu和Qiu，1993年）。近期的发展还将这些方法应用于非凸空间域（Ferraccioli等人，2021年）、曲面（Arnone等人，2022年）以及受网络约束的区域（Clemente等人，2026年；Panzeri等人，2025年），同时还用于时空点过程分析（Begu等人，2024年）。在这一研究领域中，当前的工作提出了首个专门用于估计连续标记点过程联合强度的惩罚似然框架。我们确立了相应估计问题的适定性，并证明了所得估计量的一致性。通过模拟研究以及两个真实数据应用——加利福尼亚州的野火事件和全球的地震事件——对该方法的性能进行了评估。本文的结构如下：第2节介绍所提出的方法；第3节阐述函数嵌入方法，并探讨估计量的理论性质，包括唯一性和一致性；第4节描述数值实现方法，涵盖张量积离散化、优化算法以及平滑参数的选择；第5节通过在平面空间域上进行的两次模拟研究来评估该方法的准确性，并将其与现有的多元核密度估计以及在条件可分离性假设下的核平滑方法进行比较。此外，还在曲面上进行了模拟，以展示所提框架的灵活性，目前尚无其他竞争性方法可用于此类场景；第6节介绍了针对野火事件的真实数据应用，第7节则将分析扩展到全球的地震事件；第8节最后总结了研究成果并提出了未来研究的方向。

标记非均匀泊松过程的强度估计
我们考虑一个标记非均匀泊松过程的实现形式{(si,mi):si∈S,mi∈M}i=1n，其中n表示观测到的事件数量。这些空间位置si位于紧凑的空间域S?Rd中，d可取2或3，该空间域可能代表一个非凸的平面区域，也可能是嵌入在R3中的光滑黎曼流形。与之对应的标记mi则取值于有界区间M=[M0,M1]?R+内。我们的目标是在完全非参数的框架下估计联合强度函数λ(s,m)，而不需要施加任何额外的假设。

理论性质
在本节中，我们为所提出的方法奠定了理论基础。在3.1节中，我们研究了在适当功能空间内该估计问题的适定性。在3.2节中，我们分析了估计量的渐近行为，包括其一致性，以及偏差和方差的大小。

数值估计
在本节中，我们提出了一种数值程序，用于近似解决(1)式中目标函数的无限维最小化问题。在4.1节中，我们基于有限元和B样条引入了对空间域的离散化处理，从而得到了(1)式的可行离散形式。随后，利用Begu等人（2024年）描述的优化和初始化策略来求解这个凸优化问题。最后，在4.2节中，我们讨论了平滑参数的选择问题。

模拟研究
为了评估所提方法的性能，我们进行了三项模拟研究。5.1节介绍了标准的竞争性方法，即基于条件可分离性假设的多元核密度估计和核平滑方法。这些方法仅适用于平面空间域；关于标记点过程中基于核的强度估计的详细讨论，可参考Baddeley等人（2016年，第15章）。5.2节展示了一项基于……的模拟研究。

加利福尼亚州的野火发生情况及其燃烧面积
野火指的是自然植被发生的失控火灾，它们会带来严重的环境和生态风险。要有效预防和减轻野火危害，就需要了解点火点的空间分布及其造成的影响。燃烧面积是衡量影响程度的关键指标，它可以作为生物量损失、温室气体排放以及生态系统受影响程度的替代指标（Gajendiran等人，2024年）。该数据集包含了2013年5月期间记录的野火点火点的经度-纬度坐标信息。

全球不同震级的地震发生情况
地震是具有空间分布特征的事件，其震级用于量化其所释放的能量。准确建模地震的空间-震级分布对于理解全球地震活动模式以及事件位置与震级之间的关联至关重要。在此背景下，强度函数λ(s,m)是一种重要的分析工具，可用于描述地震发生率在地理位置和震级方面的共同变化规律。估计这一函数对于……

结论
在这项工作中，我们提出了一种用于连续标记点过程联合强度的非参数惩罚似然估计器，该估计器基于偏微分方程驱动的正则化方案，并采用了新颖的标记依赖型正则化方法。该方法无需遵循严格的可分离性假设，能够处理强右偏特征的标记分布，还能自然地适配几何结构复杂的区域。第5节中的模拟实验表明，与基于核的方法相比，该方法具有明显优势。

致谢
我们感谢责任编辑和审稿人提出的极具建设性且富有价值的意见。B. Begu的研究得到了爱尔兰研究委员会（SFI）以及NexSys战略合作伙伴计划下的资助（项目编号21/SPP/3756）。L.M. Sangalli得益于意大利大学与研究部（MUR）以及欧盟NextGenerationEU计划资助的PRIN研究项目CoEnv - 复杂环境数据与建模（项目编号CUP 2022E3RY23）。L.M. Sangalli还……

Blerta Begu|Eleonora Arnone|Laura M. Sangalli|Janine Illian|Michelle Carey