通过基于物理的混合扭曲高斯过程（physics-informed hybrid warped Gaussian Processes）方法提高工业蒸馏过程中结垢情况的预测精度

《Computers & Chemical Engineering》：Improving fouling forecasting in industrial distillation via physics-informed hybrid warped Gaussian Processes

　　弗朗切斯科·内格里（Francesco Negri）| 安德烈亚·加莱亚齐（Andrea Galeazzi）| 弗朗切斯科·加洛（Francesco Gallo）| 弗拉维奥·马内蒂（Flavio Manenti）
Itelyum Regeneration S.p.A.
地址：Via Tavernelle, 19, Pieve Fissiraga, 26854, 意大利

**摘要**
在工艺工厂中，由污染引起的压降呈现出复杂的、非线性的变化模式，这挑战了标准的预测模型，并成为工业4.0背景下的一个重大障碍。本研究介绍了一种基于物理信息的输出变形方法的高斯过程回归（Gaussian Process Regression, GPR）机器学习方法，该方法专门设计用于捕捉污染的渐近行为。通过从真空蒸馏塔及其燃烧加热器的控制方程中解析推导出单调变形函数，将物理约束直接嵌入到GPR观测模型中。这种混合变形机制通过编码操作限制来作为长周期外推的规则化手段。使用工业精炼厂的数据，对三种配置（未变形、对数变形和混合变形）进行了滚动窗口预测评估。结果表明，混合变形在突变情况下（如快速污染或蒸汽吹扫操作）显著优于标准基线方法。通过结合自动核搜索和物理信息变形，该方法提高了关键工业资产预测维护的可靠性。

**1. 引言**
在工业4.0的背景下，数据驱动方法在工艺行业中的应用加速发展，增加的传感和连接性使得复杂操作的预测和决策支持成为可能（Alcácer 和 Cruz-Machado, 2019; Folgado 等, 2024; Venkatasubramanian 等, 2003a, 2003b, 2003c）。在这种环境下，考虑不确定性的时间序列预测对于预测维护尤为重要，因为可以利用操作数据来预测性能退化并支持干预计划（Daoutidis 等, 2024; Torres 等, 2025）。
高斯过程回归（GPR）适用于工业预测，因为它提供了灵活的函数学习能力和概率预测（Kanagawa 等, 2018; Rasmussen 和 Williams, 2006; Roberts 等, 2013）。当与自动核搜索结合使用时，GPR 可以发现可解释的协方差结构，并适应趋势、长距离相关性以及操作状态的变化（Duvenaud 等, 2013; Galeazzi 等, 2024; Kim 和 Teh, 2018; Lloyd 等, 2014; Steinruecken 等, 2019; Wilson 和 Adams, 2013）。在过程系统工程（Process Systems Engineering, PSE）中的最新应用包括预测关键操作变量，如压降、批次时间、温度增量、电池退化、腐蚀及相关工艺性能指标（Dai 等, 2022, 2019; Galeazzi 等, 2024; Hvala 和 Kocijan, 2020; Kaneko, 2023; Li 等, 2025; Negri 等, 2025b, 2025a, 2024a; Yeardley 等, 2020; Zhang 等, 2019; Zhou 等, 2022）。
尽管有这些优势，但如果观测行为偏离常见的建模假设，预测性能可能会受到影响，因为工业时间序列可能表现出偏态、重尾、界限或由物理和操作约束引起的突然变化。这种偏离可能可以通过在GP框架内结合领域知识（包括物理约束、领域边界或关键变量的概率分布知识）来得到补偿（Cross 等, 2024; Swiler 等, 2020）。

科学文献展示了多种在GP回归过程中结合领域知识的策略，这些策略通常根据问题的具体特征而互补。一种方法是通过使用结构化的协方差设计和操作员感知的核函数来强制物理约束（Dalton 等, 2024; Jidling 等, 2017; Lange-Hegermann, 2021）。此外，还可以通过引入非高斯约束因子来施加与单调性、界限、凸性等相关的不等式信息（Agrell, 2019; Riihim?ki 和 Vehtari, 2010; Veiga 和 Marrel, 2012）。此外，领域知识还可以通过基于语法的组合核搜索和在GP优化中的有信息均值/先验选择来指导核函数的选择和基函数组合（Duvenaud, 2014; Duvenaud 等, 2013; Lloyd 等, 2014; Tubbs 和 Mercang?z, 2023; Wilson 和 Adams, 2013）。其他策略则关注GP似然性，例如采用编码定制噪声特征的非高斯似然模型（Algikar 和 Mili, 2023; BinTayyash 等, 2021; Rasmussen 和 Nickisch, 2010; Vanhatalo 等, 2013; Jyl?nki 等, 2011; Kersting 等, 2007; Lázaro-Gredilla 和 Titsias, 2011; Maro?as 等, 2021; M?ller 等, 1998; Shah 等, 2014; Wang 和 Neal, 2012），或者通过变形函数转换输出空间，使高斯性在转换后的表示中更加合理（Oliveira 等, 1997; Snelson 等, 2003; Rios 和 Tobar, 2018, 2019）。在后一种情况下，转换隐式地在原始空间中引入了修改后的似然性，并可能产生更能反映观察到的随机行为的非高斯预测后验。

在这项工作中，我们提出了一种基于物理信息的混合输出变形策略，并将其与Galeazzi等人（2024）的自动核搜索程序相结合。该方法旨在通过引入由污染相关压降的控制关系派生的似然侧结构来提高预测性能，这种结构在概念上与通过核搜索进行的先验侧调整正交。具体来说，混合变形是通过操作真空包装蒸馏塔（T-401）及其进料燃烧加热器（PH-401B）的已知压降关系构建的，编码的物理约束可以规范长时间范围内的外推。该策略通过在这两个单元的操作数据集上进行滚动窗口概率预测来评估，比较了未变形的基线、对数变形的替代方案以及所提出的混合变形。结果表明，混合变形可以改善预测质量，并可能在模式变化后更快地适应，从而支持受污染影响的设备的预测维护决策制定。

本文的其余部分组织如下：第2节总结了关于将领域知识和输出变形结合到GPR中的背景和相关工作。第3节推导了针对蒸馏塔和炉子案例研究的物理信息混合变形函数。第4节描述了实验方案并报告了预测结果。第5节总结了本文并指出了未来研究的方向。

**2. 理论背景**
**2.1 高斯过程回归预测**
高斯过程回归（GPR）是一种非参数贝叶斯回归和预测方法，它提供概率预测以及连贯的不确定性量化（Kanagawa 等, 2018; Rasmussen 和 Williams, 2006; Roberts 等, 2013）。给定观测值集合{(ti, yi)}i=1N，通过协方差核k(t, t′)为潜在函数f(t)分配一个GP先验，而观测值通常通过加性高斯噪声来建模，如公式（1）所示：
(1) ?? = f(t) + ?，?～N(0, σ2n)。
核函数的选择和超参数Θ=(θk, σ2n2)决定了模型的归纳偏差，并在预测应用中影响外推行为和预测不确定性（Rasmussen 和 Williams, 2006）。
在标准GPR中，超参数通常通过最大化对数边际似然（LML）（Type-II最大似然/经验贝叶斯）来学习（Rasmussen 和 Williams, 2006）。在零均值假设下，高斯似然的LML可以表示为公式（2）：
(2) logp(y|t, Θ) = ?1/2Y?K?1Y?12log|K?1|N2log(2π)，
其中Ky是由核函数和噪声模型确定的协方差矩阵。标准GPR过程及其相关LML所涉及的数学实体在补充材料中有所报告。

虽然核函数的设计和组合主要影响潜在函数的先验，但预测性能也可能受到假设似然性的影响，即观测噪声和约束的表示方式。在工业时间序列中，可能会出现偏离高斯加性噪声假设的情况（例如偏态、重尾、界限），而这些效应并不总能仅通过先验侧调整来解决。在这项工作中，先验结构是通过自动核搜索算法选择的（Duvenaud 等, 2013; Kim 和 Teh, 2018; Lloyd 等, 2014; Steinruecken 等, 2019; Wilson 和 Adams, 2013），该算法基于Galeazzi等人（2024）提出的程序，并已应用于工业环境中的污染相关压降预测，旨在进行预测维护（Negri 等, 2025a, 2025b, 2024a）。有关采用的核搜索的更多细节将在第3.2节提供。

为了补充先验侧调整，通过输出变形来修改似然性，即通过对观测值应用确定性的单调双射T(?)，以便在变换后的空间中进行标准GP回归（Snelson 等, 2003; Maro?as 等, 2021）。表示z=T(y)时，观测模型变为公式（3）：
(3) ?? ～ GP(m(t), k(t, t′), z=T(y) = f? + ?，?～N(0, σ2n2)，
这样在变形空间中的推理仍然具有解析可行性，而预测可以通过对y=T?1(z)进行反变换映射回原始尺度。这种反变换通常在原始空间中产生非高斯预测分布，这可能更好地反映观察到的过程的约束和随机特性（Snelson 等, 2003; Maro?as 等, 2021）。对于没有额外自由参数的确定性变形，变量变换的雅可比矩阵只对θ的目标函数贡献了一个加法常数，因此不影响公式（2）中的标准超参数优化。然而，变形函数通常包括变量超参数，这些超参数通过雅可比变量变换项对似然性进行修改（Snelson 等, 2003）。因此，寻找最优超参数的GP训练问题变为最大化公式（4）中的雅可比修改后的LML：
(4) logp(y|t, Θ, ψ) = ?1/2Z?Kz?1Z?12log|Kz| ? N2log(2π) + ∑i=1N log|?T(yi; ψ)?yi|

**2.2 结合领域知识的高斯过程相关研究**
在复杂的工程应用中，从业者可能对系统具有先验信息，包括物理约束、操作限制或关键变量的预期分布特征。大量文献表明，在高斯过程模型中结合此类领域知识可以提高可靠性和预测性能，这些方法往往是互补的而非相互排斥的（Cross 等, 2024; Swiler 等, 2020）。特别是，可以通过（i）通过核函数和均值结构设计引入领域信息，（ii）通过额外因素或虚拟观测进行受限推理，以及（iii）似然侧建模（包括非高斯似然和输出空间变换）来引入领域信息（Cross 等, 2024; Swiler 等, 2020）。

在先验方面，可以通过结构化协方差设计和操作员感知的核函数来强制物理约束，包括满足已知线性关系的线性受限GP（Jidling 等, 2017），用于线性PDE/ODE约束的算法构建（Lange-Hegermann, 2021），以及精确实施线性边界条件的边界受限GP先验（Dalton 等, 2024）。在推理方面，可以通过引入非高斯约束因子来施加不等式信息（例如单调性、界限或凸性），通常通过虚拟观测来实现；一个典型的例子是使用虚拟导数观测和期望传播进行单调GP建模（Riihim?ki 和 Vehtari, 2010），并扩展到更一般的不等式约束（Agrell, 2019; Veiga 和 Marrel, 2012）。此外，领域知识还可以指导核函数的选择和基函数组合，例如通过组合核搜索和基于GP的探索和优化中的有信息均值/先验选择（Duvenaud, 2014; Duvenaud 等, 2013; Lloyd 等, 2014; Wilson 和 Adams, 2013）。

在似然方面，一种直接策略是用更好地反映观察到的随机特征的分布替换高斯观测模型。对于重尾行为下的稳健回归，常见的选择包括学生t分布和Huber型似然（Algikar 和 Mili, 2023; Jyl?nki 等, 2011; Tang 等, 2017）。对于计数值数据，常采用离散似然，如泊松分布和负二项分布（BinTayyash 等, 2021; M?ller 等, 1998）。当输出空间严格为正时，适合使用右偏的连续似然，如对数正态分布、伽马分布、指数分布、逆高斯分布和韦伯分布（Limpert 等, 2001; Rasmussen 和 Nickisch, 2010; Vanhatalo 等, 2013），而在[0,1]范围内的有界输出则适合使用贝塔分布（Jensen 等, 2013）。这些非高斯似然选择可以提高真实性，但可能会破坏共轭性，因此需要近似推理方案（Boukouvalas 等, 2012; Nickisch 等, 2018）。

另一种似然侧选项是输出变形，其中观测变量通过单调双射T(?)进行变换，以便在假设高斯性有效的变形空间中进行标准GP回归（Snelson 等, 2003; Maro?as 等, 2021; Chen 和 Gopinath, 2000; Meng 等, 2020; Saon 等, 2004）。这种策略在概念上与经典的方差稳定变换（如Box-Cox变换，Box 和 Cox, 1964; Oliveira 等, 1997）是一致的。在通用概率（GP）文献中，早期的贡献包括扭曲高斯过程（WGP），其中T是参数化的（例如，通过tanh分量的总和）并与GP超参数一起学习（Snelson等人，2003年）。后续的发展包括受Box-Cox启发的GP扭曲（BCGP）（Rios和Tobar，2018年）和组合扭曲GP模型（CWGP），其中T是通过简单的双射（例如，log、arsinh、Box-Cox、sinh-arcsinh）构建的（Rios和Tobar，2019年）。在严格单调的扭曲下，预测分位数可以通过逆变换映射回原始空间，而原始尺度上的矩通常需要数值近似，如蒙特卡洛采样或矩展开（Harrison，2010年，Hosmer等人，2008年，Mooney，1997年，Oehlert，1992年）。在这项工作中，我们专注于通过基于物理信息的混合输出扭曲进行似然侧建模，建立在扭曲GP公式（Snelson等人，2003年）的基础上。与Maro?as等人（2021年）对输出扭曲领域的观点一致，目标是通过提出反映原始输出尺度上高斯性合理偏离的变换来减轻潜在的模型错误规范，并将这些偏差编码在扭曲函数本身中。在这里，候选扭曲是通过操纵表示与污垢相关的压降物理约束的控制关系获得的，从而将可行性和其他领域信息特征直接嵌入到观测模型中。这种似然侧的修改旨在补充Galeazzi等人（2024年）提出的通过自动核搜索进行的先验侧调整，因为先验和似然在GP回归中扮演概念上正交的角色（Maro?as等人，2021年，Rasmussen和Williams，2006年，Snelson等人，2003年）。第3.1.3节详细介绍了所提出的扭曲策略及其在预测流程中的集成。

3.1. 设计基于物理信息的扭曲策略
本节为受污垢影响的真空蒸馏装置的压力降预测开发了基于物理信息的输出扭曲函数，重点关注填料蒸馏塔及其进料加热器。该策略从控制压降关系开始，并引入取决于污垢的项（在适当的情况下进行参数化），以获得编码物理约束和原始ΔP尺度上高斯性合理偏离的扭曲。然后使用这些扭曲来补充通过自动核搜索进行的GP先验调整，同时保持所引入的似然修改的物理解释。

3.1.1. 受污垢影响的填料蒸馏中的压降
案例研究的第一部分（第4.1节）涉及真空填料蒸馏塔中的压降，这里的压降按单位质量使用过的油料进行了归一化。作为参考，从Bravo等人（1986年）那里改编的（干）结构填料压降关系在方程（5）中报告，其中时间依赖性被明确表示并始终以随机方式处理：
(5) ?ΔP(t) = ΔP?(t)Foil(t) = F(t)ρg(t)UgS2(t)ΔZFoil(t)?2(t)sin2(ω)S，
其中ΔP(t)和ΔP?(t)分别是填料高度ΔZ上的质量比压降和实际压降，Foil(t)是使用过的油的进料流速。参数S和ω是固定的几何量，分别对应于填料侧边长度和波纹角度；ρg(t)和UgS(t)表示气体密度和表观速度；F(t)是摩擦系数；??(t)是填料空隙率。尽管方程（5）是严格为干填料推导出来的，但在这里用作一个简化的代理，因为它保留了当前扭曲构建所需变量的功能依赖性（Bravo等人，1986年，Rocha等人，1993年）。对于基线扭曲来说，方程（5）表现出乘法结构，表明ΔP可能受到多个随机因素的乘积/比率的影响，因此提出了一个将在本节后面讨论的第一个简化扭曲策略。
在高度污染的操作中，由于沉积物的积累，填料几何形状可能会随时间变化，因此将??(t)视为随时间变化而不是常数是合理的（Moore等人，1984年）。为了构建一个明确的考虑污垢的代理模型，我们假设ΔP(t)的时间变化主要由填料通道内的污垢层厚度控制，而其余贡献被合并为一个正参数。采用的几何抽象在图1中示出，其中填料几何形状以无量纲形式表示（为了保密，通过参考边长S0进行了缩放，不影响推导）。
将污垢与空隙率联系起来，本工作中使用的关系在方程（6）中报告（Negri等人，2025a）：
(6) ??(t) = Vvoid(t)Vtotal = ΔZ?Avoid(t)ΔZ?Atotal = λ?2(θ+δf(t))λ2，
其中面积/体积比率形式化了随着沉积物增长而减少的可用流动面积。将方程（6）代入方程（5）并将剩余项合并为一个正参数α，得到的明确考虑污垢的代理模型在方程（7）中给出（Negri等人，2025a）：
(7) ΔP(t) = αλ1?2λ(θ+δf(t))?4，
其中λ和θ是固定的与填料相关的几何参数，而α∈R+汇总了其余的操作和流体动力学贡献。在当前的构建中，δf(t)被建模为一个高斯过程；因此，方程（7）可以解释为一个混合建模层，它将机械代理关系与GPR联系起来（Sansana等人，2021年，Schweidtmann等人，2024年）。

下载：下载高分辨率图像（382KB）
下载：下载全尺寸图像
图1. 逐渐被污垢沉积物阻塞的填料通道的示意图。参数包括无量纲填料侧边长度λ和厚度θ，以及时间依赖的无量纲污垢厚度δf(t)。照片来自本工作中分析的工业装置中安装的填料。
方程（7）中的映射对于δf(t)∈?∞,λ2?θ是单调的，并在δf(t)=λ2?θ处表现出垂直渐近线。在这项工作中，渐近线是一个具有物理意义的可行性特征：它代表了与严重通道阻塞相关的操作极限，并且可以在预测结果偏离物理上不太可能的范围时作为正则化约束。

两个额外的考虑支持使用方程（7）作为扭曲设计的代理。首先，该构建不需要沿着塔均匀的污染增长：因为填料层是串联排列的，所以总体ΔP预计将由最具限制性和污染最快的部分主导，因此异质增长仍然可以与总体行为保持一致。其次，参数λ和θ使得可以明确区分不同的填料类型；这一点在这里很重要，因为分析的数据集涵盖了使用不同填料模型的多次运行，如第4节后面所讨论的。
关注提出的扭曲策略，第一个基线扭曲来自方程（5）的乘法结构。在简化假设下，认为贡献的随机因素可以视为独立同分布（i.i.d.），对数空间的中心极限定理（CLT）论证为ΔP提供了一个对数正态表示（Limpert和Stahel，2011年，Limpert等人，2001年，Lyon，2014年）；补充材料中提供了进一步的讨论。因此，第一个扭曲策略是对数扭曲T(ΔP)=log(ΔP)，它没有额外的超参数，因此除了变量变换之外不改变标准的LML优化。
为塔提出的主要扭曲是从方程（7）中的明确考虑污垢的代理推导出来的。通过反转ΔP和δf之间的双射，得到了在方程（8）中报告的混合、基于物理信息的扭曲（Negri等人，2025a）：
(8) ?T(ΔP)?δf(ΔP) = λ21?αλΔP1/4?θ，
其中λ和θ是固定的填料参数，α>0是一个扭曲超参数。方程（8）中的变换的雅可比矩阵在方程（9）中给出：
(9) Tα′(ΔP) = λ2?14αλ1/4?ΔP?5/4=α1/4λ3/48?ΔP?5/4。
由于变量变换，修改后的对数边际似然（LML）获得了对数雅可比贡献；本工作中使用的目标在方程（10）中显示：
(10) logp(ΔP∣t,Θ,α) = ?12δf?Kδf?1δf?12log|Kδf|?N2log(2π)+N4logα+3N4logλ?Nlog8?54∑i=1NlogΔPi，
其中α>0，λ>0，ΔPi>0。最大化方程（10）涉及标准GP超参数Θ和扭曲超参数α的联合优化。一个实际的初始化是设置α(0)，使得ΔP|t=0=ΔP0在δf|t=0=0时成立，其中ΔP0表示清洁条件下的设计参考压降。
为了方便，表1总结了塔的扭曲选项及其建模假设。

总的来说，这些扭曲在似然侧变换中明确了ΔP中不同来源的非高斯性。对数扭曲对应于与乘法随机结构一致的基线假设，而混合扭曲编码了ΔP与近似高斯潜在污垢厚度之间的机械联系，以及与可行性相关的特征，如由通道阻塞暗示的操作极限。

表1. 受污垢影响的真空填料蒸馏的扭曲策略。
扭曲类型 扭曲方程 假设
未扭曲 T(ΔP)=I(ΔP) ΔP可以被视为标准GP
对数扭曲 T(ΔP)=log(ΔP) ΔP源自几个近似独立同分布（i.i.d.）随机因素的乘法组合
混合扭曲 T(ΔP)=λ21?αλΔP1/4?θ ΔP主要由方程（7）中的代理控制的GP建模的污垢厚度δf决定

3.1.2. 受污垢影响的加热器中的压降
案例研究的第二部分将相同的构建逻辑应用于服务于真空塔的进料加热器，重点关注作为维护操作指标的质量比压降ΔPPH(t)。起点是Darcy-Weisbach关系（改编自Brown（2012年）），在方程（11）中报告，其中时间变化量被视为随机过程：
(11) ΔPPH(t) = ΔPPH?(t)Foil(t) = fD(t)ρ(t)Q(t)2L2DAw(t)2Foil(t)，
其中L和D表示管道长度和内径，fD(t)是Darcy摩擦系数，ρ(t)是流体密度，Q(t)是体积流量，Aw(t)是湿润截面积。在当前分析中，假设分布摩擦占主导，忽略了次要损失（入口、弯曲、接头）；如果需要，可以通过标准损失系数纳入这些贡献。
为了建模污垢，关键的随时间变化的几何因素被认为是Aw(t)，反映了减少可用流动面积的沉积物增长。采用的几何抽象在图2中显示，其中几何形状以无量纲形式表示（为了保密，通过参考直径D0进行了缩放，不影响推导），并且炉膛污垢厚度δfPH(t)被视为高斯过程。
将湿润面积表示为污垢厚度的函数得到方程（12）：
(12) Aw(t)D02 = π4[??2δfPH(t)]2，
其中?是无量纲内径。将方程（12）代入方程（11）并将剩余因素合并为一个正参数β，得到的明确考虑污垢的代理模型在方程（13）中报告：
(13) ΔPPH(t) = β??2δfPH(t)?5
与方程（7）中的填料塔代理类似，方程（13）中的模型意味着对污垢厚度的逆多项式依赖，并表现出垂直渐近线（这里在δfPH(t)=?/2处），可以解释为在严重阻塞下的可行性相关操作极限。

下载：下载高分辨率图像（530KB）
下载：下载全尺寸图像
图2. 逐渐被污垢沉积物阻塞的炉管示意图。参数包括无量纲内径?和时间依赖的无量纲炉膛污垢厚度δfPH(t)。照片来自本工作中分析的工业装置。
在定义扭曲策略时，基线对数扭曲再次受到方程（11）的乘法结构的启发：在独立同分布（i.i.d.）假设下，对数空间的CLT论证表明ΔPPH可以对数正态表示（Limpert和Stahel，2011年，Limpert等人，2001年，Lyon，2014年）。因此，第一个炉膛扭曲策略是T(ΔPPH)=log(ΔPPH)。
混合 furnace 拥扭扭曲是根据方程（13）中的代理关系反转得到的。得到的变换在方程（14）中报告：
(14) ?T(ΔPPH)?δfPH(ΔPPH) = 12??βΔPPH1/5，
其中?是固定的，β>0是一个进入修改后的LML的扭曲超参数。方程（14）中的雅可比矩阵在方程（15）中给出：
(15) Tβ′(ΔPPH) = 110β1/5ΔPPH?6/5，
这产生了方程（16）中报告的修改后的LML目标：
(16) logp(ΔPPH∣t,Θ,β) = ?12δfPH?KδfPH?1δfPH?12log|KδfPH|?N2log(2π)+N5logβ?Nlog10?65∑i=1NlogΔPPH,i，
其中β>0且ΔPPH,i>0。与塔的情况（方程（10）一样，通过联合优化标准GP超参数Θ和扭曲超参数β来进行推理。一个实际的初始化是选择β(0)，使得ΔPPH|t=0=ΔPPH,0在δfPH|t=0=0时成立，其中ΔPPH,0表示清洁条件下的设计参考压降。
为了方便，表2总结了加热器的扭曲选项及其假设。

表2. 受污垢影响的加热器的扭曲策略。
扭曲类型 扭曲方程 假设
未扭曲 T(ΔPPH)=I(ΔPPH) ΔPPH可以被视为标准GP
对数扭曲 T(ΔPPH)=log(ΔPPH) ΔPPH源自几个近似独立同分布（i.i.d.）随机因素的乘法组合
混合扭曲 T(ΔPPH)=12??βΔPPH1/5 ΔPPH主要由方程（13）中的代理控制的GP建模的炉膛污垢厚度δf决定

3.2. 将扭曲纳入自动核搜索
本工作中提出的基于物理信息的输出扭曲被集成到Galeazzi等人（2024年）的自动核搜索流程中，用于通过协方差结构的组合来选择GP先验。目标是将似然侧调整纳入相同的工作流程，以便可以共同调整先验结构和似然规范，以适应此处考虑的预测任务。此算法的未扭曲版本之前已在同一进料加热器的操作数据集上进行了测试（Galeazzi等人，2024年），后来作为预测维护策略的一部分在真空蒸馏塔上进行了评估（Negri等人，2025b）。相关工作进一步探索了在同一填料塔上的简化混合扭曲（Negri等人，2025a）。与这个已建立的预测和预测性维护流程一致，我们将问题表述为单变量时间序列预测，并使用运行时间作为高斯过程（GP）的输入，以便预测结果可以直接用维护调度的决策相关领域来表示。在这些研究的基础上，当前工作将基于物理的变形技术整合到核搜索中，以检验在污染条件下预测质量和鲁棒性是否可以得到提高。自动核搜索建立在基于语法的高斯过程协方差发现方法之上（Duvenaud等人，2013年；Lloyd等人，2014年）。在这类方法中，候选核是通过简单的运算符（例如求和与乘积）组合基本核生成的，然后根据证据或预测性能标准进行排名。典型的例子包括组合核搜索（CKS）（Duvenaud等人，2013年）和自动贝叶斯协方差发现（ABCD）（Lloyd等人，2014年；Steinruecken等人，2019年），这些方法在保留已知核组件显式求和与乘积的可解释性的同时，探索了丰富的核空间。这里采用的核搜索算法遵循Galeazzi等人（2024年）提出的贪心搜索结构，这是Duvenaud等人（2013年）CKS的一种变体。简而言之，该方法逐层构建候选协方差结构，在时间顺序的训练/验证划分上评估候选方案，并在每次迭代中保留得分最高的模型。贪心搜索树如图3所示。在第一层，候选模型包括所有基本核及其两两求和与乘积，这样可以降低过早收敛到不良局部最优解的可能性。在后续层次中，通过使用求和（+）和乘积（?）运算符将当前最佳核与每个基本核组合来形成新的候选模型。对于每个候选协方差模型，在训练子集上通过最大化（可能经过修改的）对数边际似然来估计核超参数。在未变形和对数变形的情况下，训练使用标准的高斯过程似然（LML）（第2.1节，公式（2））；在混合变形的情况下，训练使用从变形雅可比矩阵导出的修改后的目标函数（第3.1节，公式（10）用于列和公式（16）用于 furnace）。训练完成后，每个候选模型在变形空间中使用基于分位数的目标函数进行评估。在Galeazzi等人（2024年）的原始实现中，验证目标是平均绝对误差（MAE），对应于预测中位数的评估。在这里，我们改用加权区间评分（WIS）作为验证目标，以便同时考虑中位数准确性和区间行为（见下面的公式（18））。在每个搜索层次中，保留获得最低验证分数的候选模型作为该层次的最佳模型，并从此核开始进一步搜索。一旦达到最大深度，就会在所有层次中选择全局最佳模型，因为最深的模型不一定是最优的。最后，在进行预测之前，将选定的核重新拟合到完整的数据集（训练+验证）上，以减少因忽略最新观测数据而造成的信息损失（Galeazzi等人，2024年；Negri等人，2025b）。值得注意的是，当前工作对Galeazzi等人（2024年）的原始实现进行了少量调整：（i）扩大了核库；（ii）相应增加了搜索深度；（iii）用变形空间中计算的WIS替换了验证目标，如下所述。本工作中使用的核库总结在表3中。相对于Galeazzi等人（2024年），它增加了三个具有不同平滑参数ν的Matérn核，这些核可能更适合粗糙的潜在函数（Duvenaud等人，2013年；Porcu等人，2023年）。所的结果集包括线性（LIN）、径向基函数（RBF）、有理二次（RQ）和Matérn核，其ν∈{5/2,3/2,1/2}，分别表示为M52、M32和EXP（ν=1/2的Matérn通常称为指数核）。搜索深度增加了两级，以适应扩大的基础集。

本工作对验证空间中的目标函数进行了轻微调整，用WIS替换了MAE，认为WIS是其自然的泛化形式。区间依赖的组成部分是区间评分ISα，在公式（17）中为名义覆盖率1?α定义（Bracher等人，2021年）：(17) ISα(P,y) = (u?l) + 2α(l?y)1(yu)，其中l和u是（1?α）预测区间的下限和上限，y是观测值，1(?)是指示函数。在这里，WIS仅使用95%区间（即α=0.05）以及预测中位数来计算。相应的WIS表达式如公式（18）所示（Bracher等人，2021年）：(18) WISα(P,y) = 11/2 + 12/2|y?m| + α2ISα(P,y)，α=0.05，其中m表示预测中位数。公式（17）–（18）中的所有量都在变形尺度上评估，即对于z的预测分布和变形后的观测值zval=T(yval;ψ)。选择WIS作为调整选项的动机是认为预测平均可以从公式（17）中引入的ISα成分中受益。正如在相同算法框架内使用MAE进行时间序列预测和预测性维护的先前结果所示（Galeazzi等人，2024年；Negri等人，2025a；Negri等人，2025b），使用MAE作为核搜索的目标仍然是合适的。应该注意的是，确定WIS是否在外推设置中始终产生更精确的区间需要跨多个数据集进行专门的比较。

表3. 核搜索算法中的高斯过程核库。
| 核 | 方程 | 描述 |
|------------|---------------------------------|--------------------------------------|
| LIN | k(t,t′) = σb2 + σv2(t?t′) | 线性趋势，非平稳行为 |
| RBF | k(t,t′) = σ2exp?(t?t′)2 | 非常平滑的函数，无限可微 |
| RQ | k(t,t′) = σ21+(t?t′)22α?2?α | RBF的尺度混合，多尺度平滑 |
| M52 | k(t,t′) = σ21+5|t?t′|?+5(t?t′)23?2exp?5|t?t′|? | 中等平滑的路径，可二次微分 |
| M32 | k(t,t′) = σ21+3|t?t′|?exp?3|t?t′|? | 比M52粗糙，可一次微分 |
| EXP | k(t,t′) = σ2exp?|t?t′|? | 不可微分路径，变化急剧 |

最后，整个工作流程在算法1中进行了总结。对于每个候选核M和选定的变形T(y;ψ)，通过在变形空间中最大化相应的（修改后的）LML来进行训练，当ψ被优化时包括对数雅可比项∑ilog|?T/?y|（混合变形），而当T没有自由参数时则简化为标准LML最大化（未变形或对数变形）。算法1通过变量变换应用输出变形，而GP训练和推断在变形空间z=T(y;ψ)下仍按照标准的高斯过程精确回归（GPR）进行。因此，与未变形情况相比，渐近成本保持不变，每次拟合的规模为O(n3)，其中n=|Ttrain|（Rasmussen和Williams，2006年；Snelson等人，2003年）。变形仅增加了对T（及其在优化ψ时的雅可比项）的逐点评估，以及最终的逆变形映射；因此，它最多只影响与优化器重启和迭代相关的常数因子。原始输出尺度上的预测使用适当的评分规则进行评估，以总结校准和精度。主要的指标是WIS，现在直接在输出变量y上进行评估，并考虑了一个缩放变体sWIS来减轻尺度效应。本工作中使用的sWIS由作者在公式（19）中提出：(19) ? sWISα{0:K}(P,y) ? WISα{0:K}(P,y) × (|m| + |y|)/2，其中m是预测中位数，y是观测值。这种缩放受到了SMAPE中使用的归一化的启发（Armstrong，1985年；Flores，1986年；Makridakis，1993年）。由于WIS（和sWIS）仅依赖于选定的分位数，还报告了两种依赖于完整预测分布的额外评分规则：连续排名概率评分（CRPS）和缩放变体（sCRPS），在公式（20）、（21）中定义。（20）CRPS(P,y) = ∫RFP(x)?H(x?y)2dx，（21）sCRPS(P,y) = EX～P[|X?y|] + EX,X′～P[|X?X′|] + 1/2logEX,X′～P[|X?X′|]，其中FP(x)是预测P的累积分布函数，H(x?y)是以y为中心的Heaviside阶跃函数，X,X′是从P中独立抽取的样本（Bolin和Walin，2022年；Epstein，1969年；Gneiting和Raftery，2007年；Matheson和Winkler，1976年；Pic等人，2025年）。

4. 案例研究：TDA过程
本节介绍了本工作的核心工业案例研究：一个用于废油精炼厂的热脱沥青（TDA）部分，包括一个真空包装的蒸馏塔（T-401）及其相关的进料燃烧加热器（PH-401B）。这两个单元都受到污染的影响，对压降和定期清洁需求有可测量的影响。将未变形的自动核搜索程序与为这两个单元提出的对数和混合变形策略进行比较，以评估在预测性维护方面概率预测精度的潜在提升。

4.1. 工业过程和数据收集
所提出的基于物理的自动化高斯过程（GPR）预测工作流程使用Itelyum Regeneration S.p.A.（意大利洛迪的Pieve Fissiraga）拥有的工业废油再生设施的运营数据进行了验证。该工厂处理矿物润滑油，并生产再生的润滑油基础油以及真空煤油和沥青，通过顺序预闪蒸（PF）、热脱沥青（TDA）和加氢精制（HDF）阶段进行操作（Gallo，2016年）。由于本研究关注TDA部分，图4展示了该子系统的简化表示。TDA阶段在高温和深度真空下运行，将脱水进料分离成半成品馏分，随后将其送往加氢精制（Negri等人，2024b）。该图突出了本研究中考虑的两个单元（T-401和PH-401B）以及用于预测的信号。运营时间序列是从Yokogawa Exaquantum工厂信息管理系统（PIMS）中获取的，该系统与工厂的分布式控制系统（DCS）（Yokogawa，2020年）相连。对于T-401，获取了废油流量（FIC-4058）和塔底部/顶部压力（PI-4075/PI-4076）；对于PH-401B，获取了相同的废油流量和炉子入口/出口压力（PI-4301/PI-4304）。本工作中分析的关键目标是质量-specific压降，定义为按废油流量缩放的压力差：ΔP=(Pbot?Ptop)/Foil（对于塔）和ΔPPH=(Pin?Pout)/Foil（对于炉子）。

该工业场所之前曾作为未变形环境中相同自动核搜索GPR流程的测试平台。具体来说，Galeazzi等人（2024年）将该方法应用于PH-401B的运营数据集，而Negri等人（2024a）、Negri等人（2025b）则在其上进行了预测和预测性维护研究。鉴于与异质废料进料相关的污染严重性（Negri等人，2025b），图5提供了从T-401提取的清洁和污染结构填料的视觉比较，展示了该系统面临的挑战。

4.2. 数据集和滚动窗口预测
本研究中分析的数据集在图6中总结，T-401塔的运行标记为A–E（图6(a)），PH-401B的运行标记为PH_F–PH_G（图6(b)）。引入了无量纲时间变量，以确保每个系统中最长的运行时间为10个单位（塔的符号为τ，炉子的符号为τPH）。压降以缩放形式报告，以保护工业机密性；相应地，时间相关单位被省略，不影响比较分析。这些运行的持续时间和轨迹各不相同，反映了不同运行中塔填充类型的差异以及废油进料本身的成分变化，这影响了随时间的污染倾向（Negri等人，2025b）。

在所有运行中，自动核搜索用于预测在三种输出变形选择下的ΔP和ΔPPH的演变：（i）未变形，（ii）对数变形，以及（iii）混合变形（第3.1节）。未变形模型是GPR中的标准默认值。对数变形作为物理驱动的基线，与在独立同分布（i.i.d.）近似下的乘法控制结构一致。混合变形是通过反转一个明确污染的代理关系得出的，在该关系中，污染厚度被视为压降演变的主要驱动因素。用于混合变形的参数包括填充依赖的塔参数和变形超参数的初始值，见表4。在每次迭代中，都会计算预测范围内的评分规则（CRPS、sCRPS、WIS、sWIS），并取其平均值以获得每次迭代的单一值。具体来说，设第k个训练终点为τk=0.5k，则在[0,τk]区间内训练模型，并在预测范围（τk,τk+2）内对其进行评估。选择的范围长度是Negri等人（2025b）在同一研究流程中使用的两倍，以探究更长时间范围的外推行为，并突出不同变形策略之间的潜在差异。表4列出了不同数据集的混合变形参数。λ和θ表示无量纲的包装边缘长度和厚度；?表示无量纲的内管直径。α(0)和β(0)是用于初始化修改后的LML优化的混合变形超参数的初始值，分别用于柱子和炉子。数据集λθα(0)?β(0)分别为A2.0590.02941.477、B2.0590.02941.477、C1.0000.00290.799、D1.0000.00290.799、E1.4710.02940.963、PH_F–1.0000.1933和PH_G–1.0000.1667。

4.3. 不同变形策略下的预测评估
本节展示了针对每个调查过程的单元进行的预测评估。每个小节都以一个基准开始，该基准将提出的GP变体（未变形、对数变形、混合变形）与几个基线进行比较：经典时间序列模型（线性趋势OLS、ARIMA、ETS）、神经序列模型（GRU）以及带有滞后输入的自回归GP（AR-GP）。基准测试的完整结果以表格和图表的形式提供在补充材料中，随时间变化。然后展示了每个单元最长运行时间的提出的GP方法变体（未变形、对数变形、混合变形）的评分规则概况（T-401的案例E和PH-401B的案例PH_F）。这个选择主要是为了简洁起见，因为下面讨论的定性行为在每个单元内的多次运行中大体一致。未变形、对数变形和混合变形GP的所有运行的完整图形输出和扩展数值结果提供在附录A中，其中附录A.2.4.3.1中以表格形式报告了按预测来源划分的详细评分规则值。

4.3.1 T-401真空包装柱的预测评估
提出的预测框架针对一组常用的单变量时间序列外推基线进行了基准测试，以提供预测误差大小的背景（Hyndman和Athanasopoulos，2021; Lim和Zohren，2021; Makridakis等人，2020）。考虑的参考包括：（i）基于最近时间窗口拟合的滚动普通最小二乘（OLS）线性趋势外推；（ii）适用于ΔP序列的经典ARIMA模型；（iii）无季节性的Holt-Winters指数平滑（ETS）；（iv）自回归高斯过程（AR-GP），作为对滞后ΔP输入的GP回归实现；以及（v）紧凑的门控循环单元（GRU）网络。所有基线都在本工作中使用相同的滚动窗口协议进行评估，并使用WIS作为测试评分规则。

这项基准测试的目的并不是声称提出的混合变形GP在所有数据集上普遍优于所有其他方法，这需要更深入的调查，包括对每个基线进行详尽的调整。相反，目的是验证提出的GP变体获得的误差水平与同一任务上的典型预测方法大致相当。为了简洁起见，比较仅使用WIS来总结。这里报告了值得注意的配置选择：ARIMA(1,1,1)包含漂移项，而ARIMA(0,2,2)则不包含；GRU-ΔP直接在ΔP上训练，而GRU-ΔΔP首先在对ΔP的一阶差分上进行训练，然后进行积分以恢复ΔP预测；AR-GP在静态RBF核和非静态-静态组合（LIN+RBF）下进行测试，两种情况下GP超参数都是通过LML最大化估计的。WIS的平均值和标准差在表5中进行了总结。

总体而言，提出的GP变体产生的WIS值与所考虑基线的范围大致一致，表明该方法适合预测任务。混合变形GP在大多数运行中取得了最低的平均WIS值，并且其预测来源的分散度较低，但在本数据集中有两个例外：案例B中，对数变形稍微有利；案例C中，OLS线性趋势产生了最低的平均WIS值。在这两种情况下，混合变形仍然接近表现最好的模型。案例C中的ARIMA(1,1,1)条目用星号(*)标出，因为该模型在某个点上出现了分歧。因此，这里显示的值是通过移除相应的异常值获得的。有关所有方法和案例的WIS概况的额外基准测试细节在补充材料中提供。由于本工作的主要目标是量化自动核搜索GP流程中似然侧变形的边际效益，因此分析的其余部分集中在三种提出的GP配置上：未变形、对数变形和混合变形。

表5. T-401滚动窗口性能总结，通过每次运行的WIS平均值和标准差（μWIS±σWIS）来表示。提出的GP策略（未变形、对数变形、混合变形）与典型的基准进行了评估。
模型 | 案例A | 案例B | 案例C | 案例D | 案例E |
| ---- | ---- | ---- | ---- | ---- |
| 未变形GP | 0.0545±0.0414 | 0.0562±0.0305 | 0.1586±0.1038 | 0.3070±0.4202 |
| 对数变形GP | 0.0442±0.0355 | 0.0371±0.0232 | 0.1084±0.0784 | 0.2879±0.4172 |
| 混合变形GP | 0.0388±0.0248 | 0.0399±0.0229 | 0.0959±0.0411 | 0.2303±0.3777 |
| 线性趋势（OLS） | 0.0417±0.0206 | 0.0375±0.0099 | 0.0878±0.0445 | 0.2838±0.4043 |
| ETS | 0.0411±0.0245 | 0.0381±0.0129 | 0.2405±0.4907 | 0.2445±0.3420 |
| ARIMA (1,1,1) | 0.0570±0.0551 | 0.0412±0.0295 | (0.1794±0.1128)* | 0.3848±0.4928 |
| ARIMA (0,2,2) | 0.0532±0.0333 | 0.0715±0.0886 | 0.2109±0.3675 | 0.2595±0.3406 |
| AR-GP (RBF) | 0.0789±0.0724 | 0.0640±0.0402 | 0.2670±0.1647 | 0.4825±0.6080 |
| AR-GP (LIN+RBF) | 0.0712±0.0670 | 0.0583±0.0348 | 0.1683±0.0969 | 0.2635±0.3230 |

分析T-401柱的情况，图7展示了案例E在三种变形策略（未变形、对数变形和混合变形）下的滚动窗口预测性能。每个点对应一个预测来源τ?，并报告了在预测范围（τ?,τ?+2]内的评分规则平均值（例如，τ?=5是在τ∈(5,7)内评估的）。实际上，未缩放的（CRPS, WIS）和缩放的（sCRPS, sWIS）版本显示出类似的趋势；缩放主要减少了运行过程中由于幅度变化引起的漂移。此外，CRPS和WIS通常表现出类似的趋势，这表明WIS可以是一个方便的总结指标，因为它只需要计算选定的分位数，而CRPS需要完整的预测分布。图A.1中报告了其余案例（A–D）中评分规则随时间变化的图表，显示出类似的趋势。

鉴于不同的评分规则表现出大致相似的时间趋势，下面的讨论重点使用WIS作为柱子运行的代表性指标。提出的GP变体（未变形、对数变形、混合变形）的详细来源WIS（以及其他评分规则）表格在附录A.2中提供。

下载：下载高分辨率图片（745KB）
下载：下载全尺寸图片
图7. 案例E的预测评估，滚动窗口长度为2τ。评分规则比较了三种输出变形策略：未变形（蓝色圆圈）、对数变形（橙色方形）和混合变形（紫色菱形）。
在多次运行中，WIS的概况显示出一种重复的定性模式：（i）早期阶段的误差相对较高，（ii）中间阶段的误差较低且更稳定，在某些运行中（案例C–E），还有两个阶段（iii）误差可能会急剧增加（异常阶段），然后在增长中的训练窗口内发生变化后（后期阶段）误差会减少。初始阶段可能受到相对于固定预测范围的训练数据有限的影响，而中间阶段表明选定的协方差结构通常能够充分捕捉到主导动态。案例C–E中观察到的急剧增加与ΔP轨迹的突然变化一致，这在滚动窗口设置中是预期的，因为变化被吸收进训练集后性能会改善。

对于这种制度转变的过程方面解释可能是结垢老化，可能与外部干扰相结合。老化表示从最初较软的、粘稠的层到更硬的、类似焦炭的结构的形态转变，这种结构具有更高的抗侵蚀性，可能成为进一步沉积的基底并加速净结垢的生长（Coletti等人，2010；Ishiyama等人，2010；Moore等人，1984）。虽然可用数据并未唯一确定这一机制，但这种解释可能与在某些窗口期间ΔP的突然加速一致。

关于变形的效果，通常在每次运行的更规则阶段差异不大，而在制度转变前后可能会出现更大的分离。在表5中，混合变形在五次运行中的四次（A、C、D、E）中产生了最低的平均WIS值，而对数变形在案例B中略微有利，这是由于运行长度和包装类型的组合，导致在A–B中的两个运行中没有遇到异常和后期阶段。这种行为在图7的案例E中也可见，在运行的后期部分，混合模型通常显示出比其他方案更低的WIS值（例如，在τ?≈5及以后）。这些结果表明，当观察到的动态偏离训练窗口捕获的早期制度时，由变形嵌入的似然侧约束可能是有益的。

为了提供变化后制度下更直接的定性比较，图8展示了案例E的代表性预测，分别对应于异常窗口和后期窗口，在每种变形策略下。附录A.3中报告了案例A–D的类似图表。在图表中，训练/验证部分和保留的测试部分之间的分界由实线与虚线观察轨迹表示，预测均值由虚线橙色线表示，预测不确定性通过浅绿色阴影区域表示68%和95%的置信区间。虽然所有三种方法在包含足够的变更后观察结果后都可以适应，但变形模型（特别是混合规格）通常产生的预测分布具有与原始尺度上的观察数据更一致的偏度，能够更好地响应变化后的制度，这与案例E后期观察到的较低WIS水平是一致的。

观察到的混合变形在后期区域平均而言往往优于对数变形，可能反映了控制ΔP演变的主导因素的变化。由于每种变形策略编码了不同的物理假设，当相应的假设与当前的操作制度更加吻合时，可以预期性能会有所提高。在初始和中间区域，由于污染相对较轻，可以合理假设多个随机因素共同影响压降。在这种设置中，与对数变形（表1）的乘法组合一致。相比之下，异常和后期区域可能对应于更严重的污染，其中沉积物的增长成为压降的主要决定因素。这种解释与混合公式的假设相符，其中污染厚度δf被视为ΔP的主要驱动因素，而残余效应被汇总到超参数α中。

值得注意的是，即使在δf尚未成为主要贡献者的情况下，混合变形通常也能在初始和中间区域实现与对数变形相当的性能。在异常和后期区域，对数变形通常优于未变形基线，但通常不如混合模型。一种可能的解释是对数变形最适合于多个因素具有相似影响的制度；随着污染厚度的增加，这种近似可能不再适用。在这种情况下，混合变形可能需要δf不是始终主要的驱动因素，而混合变形则能够在单一似然侧规格中适应混合驱动条件和污染主导的行为。总体而言，T-401的结果表明，自动核搜索可以在相对规则的演变期间提供准确的预测，这在很大程度上独立于变形策略。当过程显示出ΔP动态的突然变化时，物理信息驱动的变形似乎在几次运行中提供了额外的鲁棒性，特别是混合规格在后期窗口中通常与较低的WIS相关联。这些观察结果与变形作为似然侧调味的预期作用一致，特别是在观察到的随机行为偏离早期良好代表的操作制度时。PH-401B喂料燃烧加热器的预测评估

所提出的预测框架也对PH-401B炉进行了基准测试，使用的是第4.3.1节中为柱子介绍的相同标准基线集合，目的是为了提供预测误差幅度的背景信息。对于柱子的测试而言，基准测试的目的并不是为了确立所提出的混合扭曲高斯过程（GP）相对于所有替代方案的系统性优势，而是为了验证所提出的GP变体在相同任务上的误差水平是否与典型的预测方法大致相当。比较结果同样仅使用WIS（Weighted Integrated Squared Error）进行总结，并且这里保留了第4.3.1节中为ARIMA、GRU和AR–GP选择的相同配置。表6报告了WIS在各个预测起点上的平均值和标准差。

总体而言，所提出的GP变体产生的WIS值与考虑的基线范围大致一致，表明该方法也适用于炉子预测任务。混合扭曲GP在PH_F案例中获得了最低的平均WIS，并且与表现最佳的GP变体（对数扭曲）基本无法区分，同时在两次运行中都展示了最低的WIS分散度。PH_G案例中标记有星号的ARIMA(0,2,2)条目反映了该拟合在某个起点产生了一个异常预测；报告的值是在移除了相应的异常值后计算得出的。包括所有方法和案例的WIS概况在内的更多基准测试细节在补充材料中提供。由于这项工作的主要目标是量化自动核搜索GP pipeline 中似然侧扭曲的边际收益，因此分析的其余部分集中在三种GP配置上：未扭曲、对数扭曲和混合扭曲。

表6. PH-401B滚动窗口性能总结，通过每次运行的WIS平均值和标准差（τPH?，μWIS±σWIS）来表示。提出的GP策略（未扭曲、对数扭曲、混合扭曲）与典型的基准进行了评估。

| 模型 | PH_F案例 | PH_G案例 |
| --- | --- | --- |
| 未扭曲GP | 0.0341±0.0252 | 0.0355±0.0149 |
| 对数扭曲GP | 0.0252±0.0146 | 0.0186±0.0088 |
| 混合扭曲GP | 0.0202±0.0105 | 0.0187±0.0079 |
| 线性趋势（OLS） | 0.0339±0.0249 | 0.0432±0.0319 |
| ETS | 0.0375±0.0358 | 0.0395±0.0218 |
| ARIMA (1,1,1) | 0.0517±0.0611 | 0.0362±0.0120 |
| ARIMA (0,2,2) | 0.0292±0.0261 | (0.0255±0.0088)* |
| GRU–ΔP | 0.0740±0.0204 | 0.1424±0.0690 |
| GRU–ΔΔP | 0.0325±0.0148 | 0.0721±0.0483 |
| AR-GP (RBF) | 0.0647±0.0412 | 0.0713±0.0415 |
| AR-GP (LIN+RBF) | 0.0749±0.0898 | 0.0462±0.0458 |

分析PH-401B炉的情况时，图9报告了使用第4.3.1节中为柱子采用的相同符号，对PH_F案例在三种扭曲策略下的滚动窗口预测评估。每个点对应一个预测起点τPH?，并报告了预测范围内的评分规则平均值（τPH?,τPH?+2]。与柱子的情况一样，不同的评分规则倾向于传达一致的信息，缩放主要影响轮廓的视觉稳定性（在这里，sCRPS在运行的中部似乎更敏感）。附录A.1提供了两个炉子案例在整个时间范围内的完整评分规则图表。

鉴于评分规则之间的相似性，讨论的重点是WIS作为PH-401B的代表性指标。所有评分规则的详细起点表格在附录A.2中提供。与T-401柱相比，PH-401B表现出不同的模式：WIS值倾向于保持在一个基线水平附近，偶尔会出现突发峰值，随后迅速回归到基线。在两次炉子运行中，扭曲策略通常相对于未扭曲基线降低了WIS，其中混合扭曲在PH_F中产生了最低的平均误差，在PH_G中的表现与对数扭曲相当。

下载：下载高分辨率图像（766KB）
下载：下载全尺寸图像

图9. PH_F案例的预测评估，使用长度为2τPH的滚动窗口。评分规则比较了三种输出扭曲策略：未扭曲（蓝色圆圈）、对数扭曲（橙色方块）和混合扭曲（紫色菱形）。

PH-401B趋势中错误激增的一个合理解释是偶尔采用了炉子蒸汽吹扫操作。在这种除焦过程中，油流暂时停止，并注入过热蒸汽以部分去除污垢沉积物；当操作恢复时，通常会观察到ΔPPH的下降，这可能会延长运行时间而无需完全停机。由于蒸汽吹扫是一个外生操作决策，因此没有任何版本的预测算法仅凭历史轨迹就能预测到它的发生。相反，相关问题是不同策略在干预数据可用后适应得有多快。

为了说明这种响应，图10比较了PH_F案例在蒸汽吹扫事件前后的代表性预测。PH_G案例的类似图表在附录A.3中报告。在干预之前，所有策略提供的平均轨迹都与观察到的演变接近；然而，未扭曲的模型可能显得相对过于自信，其区间较窄，可能会覆盖不到实际测试点，而对数和混合扭曲则倾向于产生更宽、正偏态的预测分布，这些分布与实际结果更吻合。蒸汽吹扫后立即，适应的差异变得更加明显：混合扭曲通常比其他策略更紧密地跟踪后续轨迹，而未扭曲的模型可能会将最近的下降外推为持续下降。这种行为与之前的观察结果一致，即未扭曲版本的算法可能需要额外的干预后观测才能恢复（Galeazzi等人，2024年）。在当前结果中，似然侧扭曲的一个实际优势似乎是达到误差接近基线阶段所需的干预后观测次数较少，其中混合扭曲在误差激增后最常与最低的WIS相关联。这可能是由于核搜索在更规则的扭曲空间中进行，其中ΔPPH的下降可能不如更广泛的趋势有影响力。

考虑到一个假设的预测维护场景，值得注意的是，与对数和混合扭曲相关的更宽、正偏态预测区间似乎与吹扫后的响应相当一致。因此，在权衡蒸汽吹扫与完全停机和清洁时，可以合理预期吹扫后的ΔPPH会保持在较低的95%预测范围内，这可以为这一决策提供一些指导。

下载：下载高分辨率图像（1013KB）
下载：下载全尺寸图像

图10. PH_F案例在蒸汽吹扫事件前后的预测，展示了三种扭曲策略。干预后的窗口对应于由ΔPPH的外生不连续性引起的误差峰值。特别是混合扭曲形式，在下降后往往适应得更快，并提供与原始尺度更对齐的偏态预测分布。

5. 结论

本研究调查了基于物理信息的输出扭曲作为高斯过程回归的似然侧扩展，结合了自动核搜索，以在一家使用过的石油精炼厂中的污垢相关压力降预测作为工业案例研究。所提出的方法从真空填充柱（T-401）及其喂料燃烧加热器（PH-401B）的压力降控制关系中推导出确定性的、单调的变换，旨在在观察模型中引入物理上有意义的约束，同时将先验规范留给核搜索过程。在混合公式中，得到的映射表现出操作极限行为，包括在全阻塞时的渐近线，这可以作为长时段外推的物理驱动约束。

使用有限的基准与标准基线进行比较，以了解误差幅度的背景情况，而主要分析集中在扭曲相对于未扭曲核搜索GP的边际收益上。在两个单元的多次运行中，扭曲模型通常相对于未扭曲基线降低了概率预测误差，其中混合扭曲最常与最佳的整体性能相关联。在制度中断期间，改进最为明显：对于T-401，ΔP的突然增加与可能由老化引起的更严重污垢的转变一致；对于PH-401B，在外生蒸汽吹扫程序后ΔPPH的突然下降。在这些制度变化的阶段，混合扭曲通常需要更少的干预后观测次数即可恢复到较低的误差水平，并且其在原始尺度上的预测分布倾向于正偏态，这种形状与未扭曲模型暗示的对称后验更为一致。

从应用的角度来看，将基于物理信息的混合扭曲纳入核搜索管道可能会在预测维护设置中提高预测的稳健性，尤其是在发生意外转变的情况下。对于炉子来说，得到的不确定性界限也有助于在将其与停机进行比较时，理解蒸汽吹扫的预期影响。

未来的工作可以扩展管道的两个组成部分。在先验方面，核搜索可以超越贪婪搜索（例如，束搜索、有限差异搜索或遗传策略），以减少错过竞争结构的风险。在似然方面，扭曲可以推广到输入依赖的形式（例如，ψ=ψ(t)），使用参数模型或关于扭曲超参数的贝叶斯先验。

作者贡献声明：
Francesco Negri：写作——审查与编辑，写作——原始草稿，可视化，验证，软件，方法论，调查，形式分析，数据管理，概念化。
Andrea Galeazzi：写作——审查与编辑，写作——原始草稿，验证，监督，软件，项目管理，方法论，调查，形式分析，数据管理，概念化。
Francesco Gallo：写作——原始草稿，验证，监督，资源管理，项目管理，资金获取，概念化。
Flavio Manenti：写作——原始草稿，可视化，验证，监督，资源管理，项目管理，资金获取，概念化。