移动最小二乘法 (Moving Least Squares, MLS)

移动最小二乘法 (Moving Least Squares, MLS)

移动最小二乘法是一种高级的无网格近似方法，其核心是通过在求解域内的每个感兴趣点（不一定是固定节点）处，动态地构建一个局部加权的最小二乘拟合，从而获得该点处函数值及其导数的连续、高精度近似。它不仅是构建无网格方法（如无单元伽辽金法）形函数的关键技术，本身也是一种强大的散乱数据插值与逼近工具。

下面我将为你循序渐进地解析移动最小二乘法的原理与实现。

第一步：从传统最小二乘法到“移动”概念的引入

传统最小二乘法回顾：
对于一个全局数据集，我们选取一组基函数（如多项式基 \(1, x, x^2, ...\)），通过最小化所有数据点的误差平方和，来确定一组全局最优的系数。这得到一个适用于整个区域的单一拟合函数。然而，当数据分布复杂或函数变化剧烈时，一个全局低阶多项式无法精确描述，而提高全局多项式阶数又容易导致龙格现象（在边界处剧烈振荡）。
引入“局部”与“移动”思想：
移动最小二乘法的核心突破在于放弃全局拟合。对于域内任意一个你想要获得函数近似值的点 \(\mathbf{x}\)（称为估值点或兴趣点），它并不使用所有数据点，而是只关注 \(\mathbf{x}\) 附近区域内的数据点。这个“附近”的概念通过一个权函数 \(w(\mathbf{x} - \mathbf{x}_I)\) 来体现，其中 \(\mathbf{x}_I\) 是某个数据节点（或称场节点）的位置。权函数在 \(\mathbf{x}\) 处取得最大值，并随着距离 \(||\mathbf{x} - \mathbf{x}_I||\) 增大而平滑衰减至零。衰减到零的距离定义了该估值点的影响域或支持域。这样，每个估值点的拟合都是基于其局部邻域数据独立进行的，因此是“局部”的。
“移动”的含义：
当你需要计算另一个点 \(\mathbf{y}\) 的近似值时，你会为 \(\mathbf{y}\) 重新定义一个新的权函数中心（即 \(w(\mathbf{y} - \mathbf{x}_I)\)），并基于 \(\mathbf{y}\) 的局部邻域数据再进行一次独立的最小二乘拟合。这个拟合的“中心”随着你关心的点 \(\mathbf{x}\) 的移动而移动，因此称为“移动”最小二乘法。

第二步：移动最小二乘近似的数学构建

考虑在定义域 \(\Omega\) 中有一组离散节点 \(\{\mathbf{x}_I\}_{I=1}^{NP}\)，其对应的已知函数值为 \(\{u_I\}_{I=1}^{NP}\)。我们的目标是构造一个连续的函数近似 \(u^h(\mathbf{x})\)，使其在任意点 \(\mathbf{x}\) 都能很好地逼近真实函数。

局部近似假设：
在估值点 \(\mathbf{x}\) 的一个小邻域内，我们假设函数可以表示为一组基函数的线性组合：

\[ u^h(\mathbf{x}, \bar{\mathbf{x}}) = \sum_{j=1}^{m} p_j(\bar{\mathbf{x}}) a_j(\mathbf{x}) = \mathbf{p}^T(\bar{\mathbf{x}}) \mathbf{a}(\mathbf{x}) \]

这里 \(\bar{\mathbf{x}}\) 是邻域内的局部坐标，\(\mathbf{p}(\bar{\mathbf{x}}) = [p_1(\bar{\mathbf{x}}), p_2(\bar{\mathbf{x}}), ..., p_m(\bar{\mathbf{x}})]^T\) 是 \(m\) 维的基函数向量（常用线性基或二次多项式基），\(\mathbf{a}(\mathbf{x}) = [a_1(\mathbf{x}), a_2(\mathbf{x}), ..., a_m(\mathbf{x})]^T\) 是依赖于估值点 \(\mathbf{x}\) 的系数向量。注意，系数 \(\mathbf{a}(\mathbf{x})\) 不是全局常数，而是随 \(\mathbf{x}\) 变化而变化的，这是 MLS 与全局最小二乘的本质区别。

建立局部加权离散误差泛函：
对于给定的估值点 \(\mathbf{x}\)，我们定义其局部近似与已知节点值的加权误差平方和为：

\[ J(\mathbf{a}(\mathbf{x})) = \sum_{I=1}^{NP} w_I(\mathbf{x}) \left[ \mathbf{p}^T(\mathbf{x}_I - \mathbf{x}) \mathbf{a}(\mathbf{x}) - u_I \right]^2 \]

其中 \(w_I(\mathbf{x}) = w(\mathbf{x} - \mathbf{x}_I)\) 是权重函数，它赋予了距离 \(\mathbf{x}\) 更近的节点更高的权重。只有当 \(\mathbf{x}_I\) 在 \(\mathbf{x}\) 的影响域内时，\(w_I(\mathbf{x}) > 0\)，否则为 0。因此求和实际上只在有限个节点上进行。

最小化误差求解系数：
通过令泛函 \(J\) 对系数向量 \(\mathbf{a}(\mathbf{x})\) 的导数为零，即 \(\partial J / \partial \mathbf{a} = 0\)，我们得到：

\[ \mathbf{A}(\mathbf{x}) \mathbf{a}(\mathbf{x}) = \mathbf{B}(\mathbf{x}) \mathbf{u} \]

其中：
- \(\mathbf{A}(\mathbf{x}) = \sum_{I} w_I(\mathbf{x}) \mathbf{p}(\mathbf{x}_I - \mathbf{x}) \mathbf{p}^T(\mathbf{x}_I - \mathbf{x})\) 是一个 \(m \times m\) 的矩陣，称为矩陣。
- \(\mathbf{B}(\mathbf{x}) = [w_1(\mathbf{x})\mathbf{p}(\mathbf{x}_1 - \mathbf{x}), w_2(\mathbf{x})\mathbf{p}(\mathbf{x}_2 - \mathbf{x}), ..., w_{NP}(\mathbf{x})\mathbf{p}(\mathbf{x}_{NP} - \mathbf{x})]\) 是一个 \(m \times NP\) 的矩陣。
- \(\mathbf{u} = [u_1, u_2, ..., u_{NP}]^T\) 是节点值向量。
因此，估值点 \(\mathbf{x}\) 的局部近似系数为：\(\mathbf{a}(\mathbf{x}) = \mathbf{A}^{-1}(\mathbf{x}) \mathbf{B}(\mathbf{x}) \mathbf{u}\)。

第三步：获得形函数与近似表达式

推导形函数：
将求解得到的 \(\mathbf{a}(\mathbf{x})\) 代回最初的局部近似表达式。注意，此时我们关心的是在 \(\mathbf{x}\) 点本身的近似值，所以令局部坐标 \(\bar{\mathbf{x}} = 0\)（即以估值点自身为局部坐标系原点）。于是：

\[ u^h(\mathbf{x}) = \mathbf{p}^T(\mathbf{0}) \mathbf{a}(\mathbf{x}) = \mathbf{p}^T(\mathbf{0}) \mathbf{A}^{-1}(\mathbf{x}) \mathbf{B}(\mathbf{x}) \mathbf{u} \]

通常，基向量 \(\mathbf{p}(0) = [1, 0, 0, ...]^T\)（对于常数项为1的基）。上式可以清晰地写成节点值的线性组合：

\[ u^h(\mathbf{x}) = \sum_{I=1}^{NP} \phi_I(\mathbf{x}) u_I \]

其中，\(\phi_I(\mathbf{x})\) 就是对应于节点 \(I\) 的 MLS 形函数：

\[ \phi_I(\mathbf{x}) = \mathbf{p}^T(\mathbf{0}) \mathbf{A}^{-1}(\mathbf{x}) \mathbf{b}_I(\mathbf{x})， \quad \mathbf{b}_I(\mathbf{x}) = w_I(\mathbf{x}) \mathbf{p}(\mathbf{x}_I - \mathbf{x}) \]

形函数的关键特性：
- 非插值性：通常 \(\phi_I(\mathbf{x}_J) \neq \delta_{IJ}\)（克罗内克δ）。即，形函数在自身节点处的值不一定为1，在其他节点处也不一定为0。因此 \(u^h(\mathbf{x}_I) \neq u_I\)。MLS 近似是一种拟合而非插值。
- 紧支性：由于权函数 \(w_I(\mathbf{x})\) 的紧支性，形函数 \(\phi_I(\mathbf{x})\) 也具有紧支性，其支撑域就是该节点权函数的支持域。这导致了带状稀疏矩阵，有利于数值计算。
- 高连续性：形函数 \(\phi_I(\mathbf{x})\) 的连续性由权函数 \(w(\mathbf{r})\) 和基函数 \(\mathbf{p}(\mathbf{r})\) 的连续性共同决定。即使使用低阶多项式基（如线性基），只要权函数 \(C^k\) 连续，形函数也是 \(C^k\) 连续的。这天然地能构造出高次连续性的近似空间，优于许多基于单元的插值。

第四步：权函数的选择与影响

权函数是 MLS 的灵魂，它直接决定了近似的局部性和连续性。

常见形式：最常用的是基于距离 \(r_I = ||\mathbf{x} - \mathbf{x}_I||\) 的紧支正函数。
- 样条权函数：如三次样条、四次样条。例如四次样条：\(w(r) = \begin{cases} 1 - 6s^2 + 8s^3 - 3s^4, & 0 \le s \le 1 \\ 0, & s > 1 \end{cases}\)，其中 \(s = r / d_m\)，\(d_m\) 是节点 \(I\) 的影响域半径。
- 高斯权函数： \(w(r) = \begin{cases} \exp(-(r/c)^2k) - \exp(-(d_m/c)^2k), & r \le d_m \\ 0, & r > d_m \end{cases}\)，其中 \(c\) 和 \(k\) 是形状参数。
影响域半径 \(d_m\) 的确定：
- 常取为 \(d_m = \alpha \cdot h_I\)，其中 \(h_I\) 是节点 \(I\) 的邻近节点间距，\(\alpha\) 是一个缩放因子（通常为 2.0~4.0）。
- 影响域必须足够大，以包含足够多的节点，使得矩陣 \(\mathbf{A}(\mathbf{x})\) 可逆（即局部有至少 \(m\) 个节点）。但同时，过大的影响域会削弱局部性，增加计算量，并使近似过于平滑。

第五步：导数的计算与在无网格法中的应用

导数近似：
由于 MLS 形函数 \(\phi_I(\mathbf{x})\) 是显式已知的，函数近似 \(u^h(\mathbf{x})\) 的导数可以直接通过对形函数求导得到：

\[ u^h_{,i}(\mathbf{x}) = \sum_{I=1}^{NP} \phi_{I,i}(\mathbf{x}) u_I \]

其中 \(\phi_{I,i}(\mathbf{x})\) 表示 \(\phi_I(\mathbf{x})\) 对 \(x_i\) 的偏导数。这可以通过对 \(\phi_I(\mathbf{x}) = \mathbf{p}^T(0) \mathbf{A}^{-1}(\mathbf{x})\mathbf{b}_I(\mathbf{x})\) 应用链式法则和矩阵求导公式得到，涉及 \(\mathbf{A}^{-1}\) 的导数。导数的连续性阶数比形函数本身低一阶。

在无网格伽辽金法中的应用：
MLS 最主要的应用之一是作为无网格伽辽金法（如 EFGM）的构形函数。其流程为：
- 用 MLS 构造出近似函数 \(u^h(\mathbf{x}) = \sum \phi_I(\mathbf{x}) u_I\)。
- 将 \(u^h(\mathbf{x})\) 代入控制方程的弱形式（伽辽金积分方程）。
- 由于 MLS 形函数不具备克罗内克δ性质，本质边界条件的施加不能像有限元那样直接赋值，需要采用拉格朗日乘子法、罚函数法或配点法等方法进行处理。
- 最终形成以节点参数 \(\mathbf{u}\) 为未知量的线性或非线性方程组进行求解。

总结：
移动最小二乘法通过为每个估值点独立进行局部加权最小二乘拟合，巧妙地结合了局部性（通过紧支权函数）、高连续性（通过权函数光滑性）和灵活性（基函数可选），生成了精度高、光滑性好的近似函数。它是连接离散节点数据与连续场变量的强大桥梁，是无网格方法得以摆脱网格束缚、处理大变形和动态裂纹等问题的核心数学工具。理解 MLS 的关键在于把握其“移动”的局部拟合思想，以及权函数在平衡局部精度与全局稳定性中所起的决定性作用。