薄板样条插值：最小弯曲能量背后的双调和方程

上一篇文章中，我们讨论了 Thin Plate Spline，简称 TPS，为什么可以从稀疏关键点匹配生成连续形变场。它的核心思想可以概括为一句话：在满足控制点匹配的前提下，使整体弯曲能量最小。 在上一篇文章中，我们直接给出了 TPS 的函数形式：

$$\tag{1} f(x,y)=a_0+a_1x+a_2y+\sum_{i=1}^{n}w_i\phi\left(\left\|(x,y)-p_i\right\|\right),$$

其中，径向基函数为：

$$\tag{2} \phi(r)=r^2\log r.$$

但这里其实留下了两个关键问题：

1. 为什么最小化弯曲能量会导出双调和方程？
2. 为什么二维 TPS 的核函数正好是 $r^2\log r$？

本文就继续从这两个问题出发，说明 TPS 函数形式的具体来历。

---

弯曲能量与双调和方程

先考虑一个标量函数 $f(x,y)$。在物理意义上，它可以理解为薄板在二维平面上方的高度；在图像配准中，它也可以理解为二维变换中的某一个坐标分量。TPS 在区域 $\Omega$ 内的弯曲能量定义为：

$$\tag{3} J(f)=\iint_\Omega \left(f_{xx}^2+2f_{xy}^2+f_{yy}^2\right)\,dxdy.$$

其中：

$$\tag{4} f_{xx}=\frac{\partial^2 f}{\partial x^2},\quad f_{xy}=\frac{\partial^2 f}{\partial x\partial y},\quad f_{yy}=\frac{\partial^2 f}{\partial y^2}.$$

这些二阶偏导数刻画的是曲面的弯曲程度。为了找到弯曲能量最小的函数，可以对 $f$ 引入一个很小的扰动：

$$\tag{5} f_\epsilon=f+\epsilon\eta,$$

其中，$\eta(x,y)$ 是任意扰动函数，$\epsilon$ 是一个很小的标量。如果 $f$ 是弯曲能量的极小值点，那么在 $\epsilon=0$ 处，能量的一阶变化应当为 0：

$$\tag{6} \delta J=\left.\frac{d}{d\epsilon}J(f+\epsilon\eta)\right|_{\epsilon=0}=0.$$

将 $f+\epsilon\eta$ 代入式(1)，可得：

$$\begin{equation} \tag{7} \begin{split} J(f+\epsilon\eta) &= \iint_\Omega \Bigl[ (f_{xx}+\epsilon\eta_{xx})^2 \\ &\qquad + 2(f_{xy}+\epsilon\eta_{xy})^2 + (f_{yy}+\epsilon\eta_{yy})^2 \Bigr]\,dxdy. \end{split} \end{equation}$$

接下来对 $\epsilon$ 求导，并在 $\epsilon=0$ 处取值，得到：

$$\tag{8} \delta J=2\iint_\Omega \left(f_{xx}\eta_{xx}+2f_{xy}\eta_{xy}+f_{yy}\eta_{yy}\right)\,dxdy.$$

这一步仍然含有 $\eta$ 的二阶导数。为了得到关于 $f$ 的微分方程，需要通过分部积分把导数从 $\eta$ 上转移到 $f$ 上。经过两次分部积分后，可得：

$$\tag{9} \delta J=2\iint_\Omega \left(f_{xxxx}+2f_{xxyy}+f_{yyyy}\right)\eta\,dxdy+\left.2\mathcal{B}(f,\eta)\right|_{\partial\Omega}.$$

其中，$\left.\mathcal{B}(f,\eta)\right|_{\partial\Omega}$ 表示分部积分后产生的边界项。通常假设扰动函数 $\eta$ 及其一阶导数在边界上为零，或者采用自然边界条件，则边界项可以忽略。因此有：

$$\delta J=2\iint_{\Omega}\left(f_{xxxx}+2f_{xxyy}+f_{yyyy}\right)\eta\,dxdy=0.$$

由于 $\eta$ 在区域内部是任意的，因此必须有：

$$\tag{3} f_{xxxx}+2f_{xxyy}+f_{yyyy}=\left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right)^2 f=\Delta^2 f=0.$$

其中 $\Delta = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ 是二维 Laplace 算子。上式就是双调和方程。也就是说，TPS 的最小弯曲能量问题，经过变分法会导出四阶双调和方程。

---

引入控制点约束

如果没有控制点约束，弯曲能量最小化只会导出：

$$\Delta^2f=0.$$

但 TPS 并不是任意找一张最平滑的薄板，而是要求薄板经过给定控制点。假设有 $n$ 个控制点：

$$\tag{16} p_i=(x_i,y_i),\quad i=1,\dots,n,$$

并且每个控制点对应的目标标量值为：

$$\tag{17} f(p_i)=z_i,\quad i=1,\dots,n.$$

于是，TPS 的优化问题可以写成带等式约束的形式：

$$\tag{18} \min_f J(f),\quad \text{s.t.}\quad f(p_i)=z_i,\quad i=1,\dots,n.$$

为了处理这些控制点约束，可以引入拉格朗日乘子 $\lambda_i$，构造增广泛函：

$$\tag{19} \mathcal{L}(f,\lambda)=J(f)-2\sum_{i=1}^{n}\lambda_i\left(f(p_i)-z_i\right).$$

其中，$\lambda_i$ 表示第 $i$ 个控制点约束对应的拉格朗日乘子。对 $f$ 做扰动 $f+\epsilon\eta$，并在 $\epsilon=0$ 处取一阶变分。这里简化下记号，用$\mathbf{x}=(x,y)$ 表示二维空间位置。忽略了分部积分产生的边界项，可得：

$$\tag{20} \delta\mathcal{L}=2\iint_\Omega \Delta^2 f(\mathbf{x})\eta(\mathbf{x})\,d\mathbf{x}-2\sum_{i=1}^{n}\lambda_i\eta(p_i).$$

这里的关键是，控制点约束是离散的。为了把离散点上的取值写成积分形式，可以利用 Dirac delta 函数构造：

$$\tag{21} \eta(p_i)=\iint_\Omega \eta(\mathbf{x})\delta(\mathbf{x}-p_i)\,d\mathbf{x}.$$

因此：

$$\tag{22} \begin{aligned} \sum_{i=1}^{n}\lambda_i\eta(p_i) &=\sum_{i=1}^{n}\lambda_i\iint_\Omega \eta(\mathbf{x})\delta(\mathbf{x}-p_i)\,d\mathbf{x}\\ &=\iint_\Omega \eta(\mathbf{x})\sum_{i=1}^{n}\lambda_i\delta(\mathbf{x}-p_i)\,d\mathbf{x}. \end{aligned}$$

代回一阶变分式，可得：

$$\tag{23} \delta\mathcal{L}=2\iint_\Omega \left[\Delta^2 f(\mathbf{x})-\sum_{i=1}^{n}\lambda_i\delta(\mathbf{x}-p_i)\right]\eta(\mathbf{x})\,d\mathbf{x}.$$

由于扰动函数 $\eta(\mathbf{x})$ 在区域内部是任意的，要使 $\delta\mathcal{L}=0$，必须满足：

$$\tag{24} \Delta^2 f(\mathbf{x})=\sum_{i=1}^{n}\lambda_i\delta(\mathbf{x}-p_i).$$

这说明，加入控制点约束以后，薄板不再只是满足自由弯曲的齐次双调和方程，而是受到多个离散点源的作用。换句话说，每一个控制点都可以理解为作用在薄板上的一个点源，而 TPS 曲面就是这些点源响应叠加后的结果。

---

齐次解、非齐次解与 Green 函数

到这里，我们已经从最小弯曲能量出发，得到了 TPS 对应的非齐次双调和方程：

$$\tag{25} \Delta^2 f(\mathbf{x})=\sum_{i=1}^{n}w_i\delta(\mathbf{x}-p_i).$$

这是一个带点源项的非齐次双调和方程。 根据线性微分方程的基本思想，它的解通常可以分为两部分：

$$\tag{26} f(\mathbf{x})=f_h(\mathbf{x})+f_p(\mathbf{x}).$$

其中，$f_h(\mathbf{x})$ 是齐次方程的解，满足：

$$\tag{27} \Delta^2 f_h(\mathbf{x})=0.$$

而 $f_p(\mathbf{x})$ 是由右侧点源项引起的一个特解，满足：

$$\tag{28} \Delta^2 f_p(\mathbf{x})=\sum_{i=1}^{n}w_i\delta(\mathbf{x}-p_i).$$

在 TPS 的标准形式中，齐次部分通常写成全局仿射项：

$$\tag{29} f_h(x,y)=a_0+a_1x+a_2y.$$

需要注意的是，严格来说，双调和方程 $\Delta^2 f=0$ 的齐次解空间并不只包含仿射函数。例如，$x^2$、$y^2$、$xy$ 等多项式也满足双调和方程。但在 TPS 的标准表达中，通常只保留仿射项来描述整体趋势，例如平移、旋转、缩放或剪切；而局部非刚性弯曲则交给径向基函数项来表示。

接下来的问题就是：非齐次部分 $f_p(\mathbf{x})$ 应该如何求？ 关键就在于 Green 函数。

Green 函数可以理解为单位点源作用下系统的基本响应。只要知道双调和算子对一个单位点源的响应，就可以通过线性叠加得到多个控制点共同作用下的非齐次解。这一步也解释了 TPS 核函数的来源：二维双调和算子的 Green 函数。

---

双调和算子与 Green 函数

既然控制点对应的是点源，那么我们先考虑最基本的单位点源问题：

$$\tag{30} \Delta^2 G(\mathbf{x})=\delta(\mathbf{x}).$$

这里的 $G(\mathbf{x})$ 称为双调和算子的 Green 函数，表示单位点源作用下系统产生的响应。由于双调和方程是线性的，如果知道单位点源的响应 $G(\mathbf{x})$，那么多个点源的响应可以直接叠加：

$$\tag{31} f_p(\mathbf{x})=\sum_{i=1}^{n}w_iG(\mathbf{x}-p_i).$$

也就是说，TPS 的非刚性形变项本质上就是若干个 Green 函数响应的线性叠加。因此，只要得到 $G(\mathbf{x})$ 的具体形式，就可以进一步写出 TPS 的函数形式。

---

二维双调和 Green 函数的径向形式

在二维空间中，双调和算子的 Green 函数具有径向对称性，因此它只与当前位置到点源的距离有关。令 $r=\sqrt{x^2+y^2}$, 于是有：$G(x,y)=G(r)$。需要注意的是，下面的常规函数推导默认在 $r>0$ 的区域进行。在 $r=0$ 处，点源项需要从分布意义理解，不能简单地把 $r=0$ 直接代入 $\log r$。对于单位点源问题：

$$\tag{32} \Delta^2 G(\mathbf{x})=\delta(\mathbf{x}),$$

在远离点源的位置，也就是 $r>0$ 的区域，右侧没有点源，因此有：

$$\tag{33} \Delta^2G(r)=0.$$

对于二维径向函数 $G(r)$，Laplace 算子可以写为：

$$\tag{34} \Delta G(r)=G''(r)+\frac{1}{r}G'(r).$$

为了求解 $\Delta^2G=0$，令：$H(r)=\Delta G(r)$， 那么在 $r>0$ 时，$\Delta H(r)=0$，从而有：

$$\tag{35} \begin{array}{rcll} \text{径向 Laplace 方程：} & H''(r)+\dfrac{1}{r}H'(r) &=& 0,\\[6pt] \text{两边乘以 } r：& rH''(r)+H'(r) &=& 0,\\[6pt] \text{写成全导数形式}： & \left(rH'(r)\right)' &=& 0,\\[6pt] \text{积分一次}： & rH'(r) &=& C_1,\\[6pt] \text{整理}： & H'(r) &=& \dfrac{C_1}{r},\\[6pt] \text{再积分一次}： & H(r) &= & C_1\log r+C_2.\\[6pt] \text{代入$\Delta G(r)$}： & G''(r)+\frac{1}{r}G'(r)&=&C_1\log r+C_2, \\[6pt] \text{两边乘以}r： & rG''(r)+G'(r)&=&C_1r\log r+C_2r,\\[6pt] \text{写成全导数形式}： & \left(rG'(r)\right)'&=&C_1r\log r+C_2r,\\[6pt] \end{array}$$

对 $r$ 积分一次，得到：

$$\tag{36} rG'(r)=\frac{C_1}{2}r^2\log r-\frac{C_1}{4}r^2+\frac{C_2}{2}r^2+C_3.$$

因此：

$$\tag{37} G'(r)=\frac{C_1}{2}r\log r-\frac{C_1}{4}r+\frac{C_2}{2}r+\frac{C_3}{r}.$$

继续积分，并将常数重新整理为 $A,B,C,D$，可以得到：

$$\tag{38} G(r)=Ar^2\log r+Br^2+C\log r+D.$$

这就是二维双调和方程在 $r>0$ 区域的径向解形式。 但 TPS 中最终只保留了核心项 $r^2\log r$。由于 $G(r)$ 需要在 $r=0$ 处有限，因此必须有 $C=0$。此外，$Br^2+D$ 是一个二次多项式，它的四阶导数为零，因此不会对 $\Delta^2G$ 产生贡献，可以纳入到齐次解中。于是可以令 $B=D=0$，从而得到：

$$\tag{39} G(r)=Ar^2\log r.$$

接下来只需要确定常数 $A$。直接计算 $G(r)=Ar^2\log r$ 的 Laplace 算子：

$$\tag{40} \begin{aligned} G'(r)&=A\left(2r\log r+r\right),\\ G''(r)&=A\left(2\log r+3\right),\\ \Delta G(r)&=4A\log r+4A,\\ \Delta^2G(r)&=4A\Delta\log r. \end{aligned}$$

二维 Laplace 算子的基本解满足：

$$\tag{41} \Delta\log r=2\pi\delta(\mathbf{x}).$$

所以：

$$\tag{42} \Delta^2G(r)=4A\Delta\log r=8\pi A\delta(\mathbf{x}).$$

为了满足单位点源方程 $\Delta^2G=\delta(\mathbf{x})$，必须有 $8\pi A=1$，即 $A=\frac{1}{8\pi}$。于是，二维双调和算子的 Green 函数为：

$$\tag{43} G(r)=\frac{1}{8\pi}r^2\log r.$$

由于常数系数 $\frac{1}{8\pi}$ 可以吸收到权重参数 $w_i$ 中，因此 TPS 中通常直接取径向基函数为：

$$\tag{44} \phi(r)=r^2\log r.$$

到这里，$r^2\log r$ 的来源就比较清楚了：它正是二维双调和算子在单位点源作用下的基本响应。

小结

将齐次解和非齐次解合在一起，二维空间中 TPS 的函数形式可以写成：

$$\tag{45} \begin{aligned} f(\mathbf{x}) &= f(x,y) \\ &=\underbrace{a_0+a_1x+a_2y}_{\text{全局仿射项}}+\sum_{i=1}^{n}\underbrace{w_i\phi\left(\left\|\mathbf{x}-p_i\right\|\right)}_{\text{局部非刚性弯曲项}}\\ &=a_0+a_1x+a_2y+\sum_{i=1}^{n}w_i\left\|\mathbf{x}-p_i\right\|^2\log\left\|\mathbf{x}-p_i\right\|. \end{aligned}$$

其中，$a_0+a_1x+a_2y$ 是全局仿射项，用来描述整体平移、旋转、缩放或剪切趋势。

而：

$$\tag{46} \sum_{i=1}^{n}w_i\phi\left(\left\|\mathbf{x}-p_i\right\|\right)$$

则来自多个 Dirac 点源的 Green 函数响应叠加，用来描述由控制点诱导的局部非刚性弯曲。所以，TPS 的核心并不是简单地“选择了一个看起来合适的径向基函数”，而是有一条清晰的数学来源。这也解释了为什么 TPS 能够从稀疏控制点生成连续、平滑的形变场：

控制点决定薄板必须经过哪里，而最小弯曲能量决定薄板在这些点之间如何自然地弯过去。