基于软化技术的任意平滑度保形插值

《Computer Physics Communications》：Shape preserving interpolation with arbitrary smoothness based on mollification

　　米凯莱·盖里尼·罗科 | 费德里科·斯乔斯科利 | 洛伦佐·菲吉尼
意大利米兰20126，米兰-比科卡大学物理系

摘要
在物理学和工程学中，一个常见的问题是在保留数据整体形状的同时对离散样本进行插值：例如，保持温度曲线的正性，或累积概率的单调性。现有的常用方法并不具备这一特性，虽然文献中有一些解决方案，但它们可能对于某些应用来说平滑度不足，或者计算成本过高或难以实现。本文我们提出了一种能够保留数据形状的插值方法，所得插值函数的可微性等级可任意设定。该方法适用于微分方程中的数据建模，且实现快速简单。我们的方法基于一种称为“软化”的卷积操作，通过将该操作应用于更简单的形状保持插值函数，从而获得所需的平滑程度。借助适当的多项式软化函数，该卷积具有封闭形式，因此可以快速计算，从而得到Cn阶插值函数。我们为任意软化函数定义了通用方案，给出了误差分析，并提供了C2阶的具体实现方式。最后，我们展示了该方法如何成功应用于电子回旋波射线追踪程序GRAY中等离子体密度的建模。

引言
在基于测量数据开展计算机模拟以建模物理量时，通常需要进行插值。除了构建平滑函数外，我们可能还需要保留数据的形状，这里所说的形状指的是某些导数的符号。例如，为了避免出现不符合物理规律的结果，某些温度读数的插值结果应始终为正；对于累积概率而言，其值应始终非递减，这样其导数才能表示分布函数。最常见的插值方法——样条插值和厄米插值——通常不具备形状保持功能：例如，用于插值正数值的立方样条，如果样本值接近零，就很容易出现负值。多年来，人们提出了多种解决此问题的插值方案[1]、[2]、[3]、[4]。然而，我们发现这些方案都不适合我们的应用场景：即为需要用高阶积分器求解的常微分方程组构建数据。这要求数据至少具有C2阶连续性，包括在边界外进行外推时也是如此。（更多细节及各种方案的比较见第4节）

为了解决这些缺陷，我们开发了一种新的插值方案，采用了函数分析领域中的一种称为“软化”的技术。我们首先对数据进行分段线性插值，这种方法本身就能保留数据的形状，但显然平滑度不够。不过，通过将数据与软化函数——即具有紧支集的平滑函数——进行卷积，初始的粗糙插值函数就会变得平滑，同时仍保持其形状保持特性。特别是，软化后的插值函数会继承软化函数的可微性等级。通常，软化函数是具有紧支集的C∞阶函数，仅作为理论工具使用，但若将其放宽为在边界处具有Cn阶连续性的多项式，则该卷积可以通过解析方法计算，从而便于实际数值实现。

在信号处理等领域，利用卷积来平滑数据是一种常见做法，其中一些滤波器，尤其是有限脉冲响应（FIR）滤波器，可视为离散卷积。不过，我们的方法基于连续卷积，在这方面更类似于非参数统计中的核方法。通过将软化函数的标准定义从C∞阶函数放宽为在边界处具有Cn阶可微性的多项式，使得我们可以构建具有简单封闭形式的插值函数，而且关键的是，这类插值函数具备所需的连续性、形状保持性和紧支集特性。

与厄米插值不同，但与平滑样条类似，这种方法并非严格的插值方法：插值函数与样本点之间的距离，即误差，并非为零，但通过调整软化函数的宽度，可以将误差缩小到任意小的程度。作为一个附带效果，软化操作实际上起到了低通滤波器的作用，能够有效消除噪声。当数据来自不确定的测量结果且在插值前未经过处理时，这一特性十分有用。

在下一节中，我们将简要介绍软化函数并证明其主要性质。本节的所有结论都是众所周知的，不过通常只在C∞阶情况下被阐述。我们的插值方法的平滑性和形状保持性可直接从这些结论中得出。在第3节中，我们将为具有任意可微性等级的通用软化函数定义该插值方案，给出误差分析，并提供使用C2阶软化函数的详细实现方式。最后，在第4节中，我们将探讨其他形状保持方法，并对比它们与我们方法的优劣。

章节摘录
**软化函数理论**
我们将软化函数定义为具有以下性质的连续函数g:Rd→C [5]：
1. supp(g)={x∈Rd∣g(x)≠0}是紧集；
2. ∫g(x)dx=1；
3. 对于几乎所有的x0以及所有属于Lloc1(Rd)的函数f，都有limw→0+w?d∫g((x0?x)/w))f(x)dx=f(x0)。

考虑具有紧支集的Cn阶函数构成的空间D≡Ccn(Rd)。我们取一个属于D的函数g（满足性质1），如有必要，对其进行归一化处理，使其也满足性质2，然后……

**定义**
我们考虑一组样本点(xi,yi)∈R2，其中i=1...N，且满足x1< />• χ为特征函数。

需要注意的是，f仅属于cc0([x1,xn])，但它能保持数据的形状，至少在一阶导数层面如此。实际上，假设x∈uj且所有yj的符号相同，则：
1. signf(x)=sign(δj/hj(x?xi)+yj)=signyj；
2. signdf(x)=signδj。

**与现有方案的比较**
我们简要介绍了一些文献中提出的形状保持插值方案，以及它们与我们所提方法的差异。标准厄米立方插值的一个著名扩展是fritsch-carlson算法[1]，该算法通过调整数据估计导数的值来保持数据的单调性。由于厄米插值是局部操作，因此该方法既可用于分段单调数据，也可用于全局数据。不过，该方案仅能保证c1阶连续性。

**结论**
我们提出了一种基于软化技术的新型插值方案，该方案具有形状保持功能——即能够保持数据的符号、单调性和凸性——并且可以实现任意程度的平滑。这一方法非常适合用于物理模拟，尤其是为微分方程构建数据，因为在这些应用中往往需要具有cn阶连续性和紧支集的函数。我们还提供了使用六阶多项式软化函数的c2阶实现方案，并证明了实现更高阶连续性是可行的。

**作者贡献声明**
米凯莱·盖里尼·罗科：写作——审稿与编辑、写作——初稿撰写、可视化、验证、软件开发、方法设计、研究分析、形式化分析、概念构思。
费德里科·斯乔斯科利：写作——审稿与编辑、验证。
洛伦佐·菲吉尼：写作——审稿与编辑、写作——初稿撰写、验证。

**利益冲突声明**
作者声明不存在任何已知的财务利益或个人关系，这些利益或关系可能会影响本文的研究工作。 • χ为特征函数。 需要注意的是，f仅属于cc0([x1,xn])，但它能保持数据的形状，至少在一阶导数层面如此。实际上，假设x∈uj且所有yj的符号相同，则： 1. signf(x)=sign(Δj/hj(x?xi)+yj)=signyj； 2. signdf(x)=signΔj。 **与现有方案的比较** 我们简要介绍了一些文献中提出的形状保持插值方案，以及它们与我们所提方法的差异。标准厄米立方插值的一个著名扩展是fritsch-carlson算法[1]，该算法通过调整数据估计导数的值来保持数据的单调性。由于厄米插值是局部操作，因此该方法既可用于分段单调数据，也可用于全局数据。不过，该方案仅能保证c1阶连续性。 **结论** 我们提出了一种基于软化技术的新型插值方案，该方案具有形状保持功能——即能够保持数据的符号、单调性和凸性——并且可以实现任意程度的平滑。这一方法非常适合用于物理模拟，尤其是为微分方程构建数据，因为在这些应用中往往需要具有cn阶连续性和紧支集的函数。我们还提供了使用六阶多项式软化函数的c2阶实现方案，并证明了实现更高阶连续性是可行的。 **作者贡献声明** 米凯莱·盖里尼·罗科：写作——审稿与编辑、写作——初稿撰写、可视化、验证、软件开发、方法设计、研究分析、形式化分析、概念构思。 费德里科·斯乔斯科利：写作——审稿与编辑、验证。 洛伦佐·菲吉尼：写作——审稿与编辑、写作——初稿撰写、验证。 **利益冲突声明**>• χ为特征函数。

需要注意的是，f仅属于cc0([x1,xn])，但它能保持数据的形状，至少在一阶导数层面如此。实际上，假设x∈uj且所有yj的符号相同，则：
1. signf(x)=sign(δj/hj(x?xi)+yj)=signyj；
2. signdf(x)=signδj。

**与现有方案的比较**
我们简要介绍了一些文献中提出的形状保持插值方案，以及它们与我们所提方法的差异。标准厄米立方插值的一个著名扩展是fritsch-carlson算法[1]，该算法通过调整数据估计导数的值来保持数据的单调性。由于厄米插值是局部操作，因此该方法既可用于分段单调数据，也可用于全局数据。不过，该方案仅能保证c1阶连续性。

**结论**
我们提出了一种基于软化技术的新型插值方案，该方案具有形状保持功能——即能够保持数据的符号、单调性和凸性——并且可以实现任意程度的平滑。这一方法非常适合用于物理模拟，尤其是为微分方程构建数据，因为在这些应用中往往需要具有cn阶连续性和紧支集的函数。我们还提供了使用六阶多项式软化函数的c2阶实现方案，并证明了实现更高阶连续性是可行的。

**作者贡献声明**
米凯莱·盖里尼·罗科：写作——审稿与编辑、写作——初稿撰写、可视化、验证、软件开发、方法设计、研究分析、形式化分析、概念构思。
费德里科·斯乔斯科利：写作——审稿与编辑、验证。
洛伦佐·菲吉尼：写作——审稿与编辑、写作——初稿撰写、验证。

**利益冲突声明**
作者声明不存在任何已知的财务利益或个人关系，这些利益或关系可能会影响本文的研究工作。>