基于支持向量回归的动网格技术研究

高 翔，廖海翔，徐传福

(国防科技大学 计算机学院，量子信息研究所兼高性能计算国家重点实验室，长沙 410073)

0 引 言

动网格生成技术是指网格做相应调整以适应模拟过程中计算区域的改变，是计算流体力学模拟运动界面流动问题的主要瓶颈之一。动网格方法在气动外形优化设计、流固耦合模拟和气动弹性计算等领域有着广泛的应用。为提升计算精度和结果可靠性，亟需一种高效健壮的动网格方法重新生成或移动网格以满足计算域的变化。第一种方法是直接在各时间步内重构整体或局部的计算网格，但复杂区域的网格生成十分耗时，且新旧网格间的插值传递会引入额外的误差[1]。第二种方法只变化网格节点坐标，但不改变连接关系，从而避免流场的插值传递，并能处理各种类型的网格单元，因此得到了广泛使用，但这种方法很难直接用于多体分离等大尺度运动问题。

在上述第二种方法中，超限插值法（transfinite interpolation，TFI）按网格方向顺序根据边界点与内部点的相对位置更新网格，该方法快速有效但只适用于结构网格[2]。对于没有规则方向的非结构网格运动，一般采用以下两类方法：第一类将网格边类比为弹簧[3]或弹性实体[4]，根据网格拓扑构建力平衡系统，因此需反复求解规模为网格节点数目的线性方程组，计算成本很大。第二类只基于与边界点的距离计算各网格内部点的新坐标，由于计算量少且易于并行，近年来得到了广泛应用。基于Delaunay 背景网格的插值方法先建立计算网格节点与背景单元的映射关系，边界运动后网格节点即可根据相应背景单元的边界点快速计算得到新坐标[5-7]。该方法效率较高且能推广至三维空间，但目前还只能用于凸区域的网格运动。逆距加权插值法（inverse distance weighting，IDW）通过显式函数移动网格内点，不需要求解隐式方程，但对于不同的运动形式和边界类型需采用不同的指数函数，健壮性较差[8-11]。径向基函数插值法（radial basis function，RBF）根据边界点坐标和运动位移拟合一个径向基函数来表示整个计算域内网格点位移的分布。由此需求解一个规模为边界点数目的线性方程组得到函数的待定系数[12]。为进一步提高其计算效率，可只选择部分边界点作为控制点来计算分布函数。主要思路是基于贪心思想，给定初始的控制点集合，由此得到拟合函数后计算剩余边界点的位移，然后再与其已知位移比较并将差异最大的若干个边界点加入该控制点集合，继而继续迭代直至满足设置的误差条件[13-14]。这种方法能够大幅度降低动网格的计算和存储成本，因此被广泛应用于各类问题。

IDW 方法和RBF 方法最初都是加权处理所有边界点对各内部点的影响，因此效率相对较低。而Delaunay 背景网格插值法只基于背景单元的几个边界点影响相应内部点，变形能力又相对较弱。近年来，也出现了一些结合背景网格方法与IDW 方法或RBF 方法的研究[11,15]，主要思路是先将IDW 或RBF方法用于稀疏的背景网格运动变形，然后也维持计算网格节点与背景网格单元的相对位置不变，从而插值计算出网格节点的新坐标。该方法综合了两类算法的优点，但需要另外再生成一套计算域的稀疏背景网格。基于插值或拟合函数，尝试选择部分边界点与所有内部点、甚至不同边界点与不同内部点建立相应映射关系，从而优化和加速动网格方法，这是一个值得继续研究的方向。

本文基于上述思想，进一步放宽贪心选点的径向基函数插值法的约束条件，提出并完善了一种基于支持向量回归（support vector regression，SVR）的动网格方法，该方法将贪心选点的插值问题转化为自适应选点的优化问题。在前期的工作中，已初步基于高斯核函数与贪心RBF 方法进行了对比，在保持网格运动健壮性的同时SVR 方法具有更高的计算效率。本文进一步通过系统性测试分析，以高斯核函数的SVR 方法为基准，针对SVR 动网格方法给出了几个性质良好的核函数及其参数设置策略，并在若干典型二维/三维动网格案例上进行了成功应用，进一步优化和完善了SVR 动网格算法。

1 方法描述

1.1 贪心选点RBF 方法约束条件的放宽

RBF 动网格方法在各坐标方向分别采用一个RBF 函数，以描述整个计算域内任意点在该方向将要移动的位移量，一般的表达形式如下：

求解上述线性方程组确定RBF 函数的系数后，从而得到了网格的运动函数。

为进一步提高RBF 方法在求解大规模问题时的计算效率，贪心选点的思想是逐次新增若干边界点作为控制点，由此重新计算得到RBF 插值函数，最终将其他非控制点的边界点位移量误差控制在一定范围内，从而减少待求方程组的规模。其数学内涵可以描述如下：

若进一步将公式（6）的约束条件放宽，使得对所有边界点都允许计算得到的位移与已知位移存在一定误差，即：

并在此基础上求解某个最值问题，则将贪心选点的RBF 插值方法转换成了一个约束优化问题。而支持向量回归模型正是解决这类问题的高效方法，且其回归函数的表达形式也是由若干径向基函数线性组合而成。因此，完全可以尝试采用支持向量回归作为网格运动的位移分布函数。

1.2 支持向量回归

支持向量回归的机器学习方法发展于统计学习的理论，提出近三十年来已广泛应用于诸多领域的模式识别、分类以及回归分析。如人脸识别、文本分类以及复杂工程分析的近似等[16-18]。SVR 的核心原理是基于一个非线性核函数将样本空间映射到高维特征空间，并将原问题转化为在此高维空间求解一个线性拟合问题[19]。基于支持向量的概念，SVR 所得到的拟合函数即要求所有样本点都落在给定的误差范围内，并使得函数在特征空间的斜率达到最小。

将SVR 应用于求解网格运动时，本文提出的方法把所有边界点坐标及其某一坐标方向的已知位移量作为训练样本集，令其为：

其中假设线性情况下的拟合函数为：

如图1 所示的样本集拟合函数，所有训练样本点都落在函数向两侧平移 ε而形成的管道内。从中可知该函数实际上仅取决于部分样本点，这些实质决定回归函数的训练样本被称为支持向量，也就相当于贪心选点RBF 方法中选择出来的控制点。

图1 SVR 在样本空间的示意图Fig. 1 Schematic of SVR in the sample space

公式（9）的二次规划问题一般采用拉格朗日乘数法（Lagrange Multiplier）转化为对偶问题再更有效地求解。对偶后的函数可表示为如下形式：

此外，采用非线性核函数k(x,xi)代 替内积 〈x,xi〉则可求解更高维特征空间的线性回归。这等价于将样本数据高维扩展后再做内积，其结果等同于将样本代入核函数得到的值[20]。通过这种方式从而极大提高了SVR 的拟合能力，扩展后的回归函数为：

关于支持向量回归对偶转化和非线性扩展更详细的推导，可参阅文献[21]。

1.3 SVR 针对网格运动的适配

针对网格运动，本方法与贪心选点的RBF 插值法类似，也需要寻找合适的核函数k(x,xi)和误差阈值ε，使之具有较强的拟合能力和变形能力。此外，适应各类运动的核函数和通用的参数设置方法对于动网格算法的推广和应用也非常重要。前期相关研究初步采用了应用广泛的高斯核函数处理网格运动，相比于贪心选点的RBF 方法，在保持相当网格质量的基础上计算效率提高了数倍[21]。本文的一个重要研究内容是通过系统全面地测试，分析各类核函数在SVR 动网格方法框架下的性能表现，具体内容见第2 节。

对于SVR 方法的误差阈值ε，为了减少累积误差并避免网格发生交错，阈值 ε一般应不大于网格点间最短距离Δmin与 变形迭代次数s的比值，即有：

在此基础上，误差阈值 ε越小，SVR 方法的拟合能力越强，但相应的建模效率可能会降低。一般可根据网格运动的复杂程度，在0 到1 之间粗略设置阈值系数 λ即可。

通常而言，一般的SVR 应用会采用归一化方法将样本各个维度的数据缩放至[0, 1]范围内。由于网格坐标各维度采用同一度量单位，且归一化还会降低数据的精度，因此针对网格运动不需要进行这一操作。除此之外，因为所有网格边界点都是有效数据，并且边界在运动形变中通常是连续的，故而不存在所谓的噪声数据，也就无需在算法模型中采取松弛变量处理例外[18]。

最后，针对SVR 动网格方法对应的特殊二次规划问题，可采用序列最小最优化（sequential minimal optimization，SMO）算法[22]进行高效求解。SMO 算法是在每次迭代中只选择两个变量而令其他为常数，从而将原问题分解为若干尽可能小的二次规划问题，再采用解析方法直接求解。因此在每个时间步内进行大规模网格运动时，相对于迭代求解规模递增方程组的贪心选点RBF 方法，SVR 方法有着很大的效率优势。此外，SMO 算法从训练样本集全局考虑，自适应选择最合适的边界点作为支持向量，且无需选择若干边界点作为初始的控制点集合，进一步减少了人工干预。

2 基准测试

SVR 动网格方法总体分为两步，一是利用已知位移的边界点建立SVR 模型，二是利用该模型预测计算内部点的位移。其中模型建立的时间开销通常更大，且直接决定第二步的网格质量，因此本节主要关注第一步的计算效率和对边界点位移值的拟合效果。对于内部点运动后的最终网格质量评估将在第3 节通过案例应用测试分析。

2.1 测试案例设计

为测试不同核函数针对网格边界运动规律的拟合能力、计算效率以及不同疏密网格下的健壮性，设计如图2 所示的三套不同密度的二维笛卡尔网格。这些网格x、y方向的坐标范围为[-2, 2]，z方向坐标为0，网格单元数分别为 2 0×20、 4 0×40和 8 0×80。

图2 三套不同密度的二维笛卡尔网格Fig. 2 Three 2D Cartesian grids with different mesh densities

在实际动网格案例中SVR 方法是分别在三个方向独立建模，因此研究其拟合能力时只需针对某一方向测试即可。根据不同类型的网格运动，本节采用四个相对复杂的解析函数指导测试网格在z方向的变形，其z坐标的位移值等于用该点坐标代入相应函数所得。正值表示网格点沿z轴正向运动，负值表示沿z轴反向运动。4 个运动函数如下：

变形后的网格如图3 所示，经F1 函数变形后的网格特点是平缓延展、曲面光滑；F2 函数变形后如对称的涡流；F3 函数变形后网格呈现周期性上下起伏；F4 函数变形则同时具备周期性、激变性、对称性等多种运动特点。因此，这些网格运动具有很好的代表性。

图3 四种运动函数对应变形后的网格Fig. 3 Deformed grids for four motion functions

将初始网格的坐标点和运动变形后的位移量作为不同核函数下SVR 方法的训练样本，采用SMO 算法训练得到相应的SVR 模型。再将初始网格点代入这些SVR 模型计算得到预测的位移量，从而可通过比较预测值和真实值的误差、选择的支持向量（控制点）考察不同核函数的拟合能力及训练开销等性能。

本文采用均方差（mean square error，MSE）来衡量回归拟合的精确度，其计算公式如下：

其中，zi表 示网格点i的 真实位移值，～zi为预测位移值，Nb为边界样本点数目。MSE 值越大，说明对应SVR模型对网格运动的表示能力越差。

2.2 核函数类型

对于某个确定的回归问题，Vapnik 等证明必然存在可行的核函数，且满足Mercer 定理的函数均能作为核函数[16]。由于满足Mercer 定理的函数较多，对于特定类型的问题，通常需要根据经验或大量实验对比才能给出更适合的核函数及其参数，从而能显著提高SVR 的建模效率和质量。径向基核函数是处理无先验知识数据集首选的核函数类型，其以空间中两点的欧氏距离为自变量，影响周围区域的样本取值。

按照作用范围的不同，径向基函数可分为全支撑（Global Support）基函数和紧支撑（Compact Support）基函数两类。全支撑基函数在整个定义域内均为非零值，其常见的函数类型如表1 所示，其中MQB 和IMQB 函数使用参数a控制基函数的形状，本文设置其常用值a=10-3。参考标准高斯核函数（表1 中的Gauss 函数），本文进一步将全支撑基函数的自变量通过系数r进行缩放，从而控制核函数对样本差异的敏感程度，即：

表1 常见的全支撑基函数Table 1 Common global support basis functions

表2 常见的紧支撑基函数Table 2 Common compact support basis functions

新的SVR 动网格方法将基于以上核函数进行测试，更多有关径向基函数的信息及其误差和收敛特性可参阅文献[23]。

2.3 结果分析

除了前期研究采用的高斯核函数外，本节通过大范围调整核参数，对其余13 个核函数进行了大量拟合测试，得到了7 个能较好拟合公式（14）描述的四类网格运动的核函数及其核参数设置范围，结果如表3 所示。其余核函数均不适用于SVR 动网格方法，后续不再讨论。

表3 成功拟合网格运动的核函数Table 3 Kernel functions for successfully fitting mesh motions

当SVR 方法采用原来的高斯核函数时，基于测试案例 40×40的中等网格拟合四类运动，求解回归方程占总的时间开销均在85%以上。这应该是由于其为指数函数，运算复杂度较高。为评估不同核函数的拟合效率，本文进一步测试了它们在 40×40网格上的计算开销，结果如图4 所示。可以看出，高斯核函数的计算速度最慢，尤其是在拟合运动特性复杂的F4 函数时，表现明显落后于除 CTPSC0以外的其他核函数。此外，IMQB、IQB 和 C PC2的效率最高。

图4 不同核函数在网格上的时间开销Fig. 4 Time cost on the grid for different kernel functions

为考察边界点规模对SVR 方法拟合性能的影响，采用运动特性较全面的F4 函数测试了核函数在不同密度网格下（分别用20、40 和80 表示）完成变形的时间开销、控制点选择和均方差。图5 给出了各类核函数在不同密度网格上的时间开销，从中可知边界点规模增大会明显降低建模的速度，但不同核函数的影响差异显著。高斯核函数和 CPC6函数在边界点增加时效率下降严重，而其他核函数的效率下降相对缓慢。

图5 不同密度网格上各类核函数时间开销Fig. 5 Time cost on different density grids for various kernel functions

图6 展示了不同核函数选择控制点占边界点数的比例，可以看出在处理较密两套网格时高斯核函数和 CPC6函数依赖的控制点数更多，由于SVR 建模速度主要取决于选择的支持向量数目，这也解释了为何这两个函数的时间开销更大。其他核函数则都能在不改变核参数的情况下将控制点数量维持在一定范围内，这一特性在处理高密度的大规模网格变形时将具有极大优势。

图6 不同密度网格上各类核函数选择的控制点数量Fig. 6 Number of control points selected by various kernel functions on different density grids

按照径向基函数分类，全支撑组和紧支撑组选择的控制点分布分别如图7（a）和图7（b）所示。从图中可明显看出，对于所有核函数，控制点的总体分布不会随网格密度增加而有明显改变。此外，不同核函数所选择的控制点虽然分布半径和密度有所区别，但基本都吻合同心圆样式的排列，这正是图3（d）中F4 函数运动变化最剧烈的区域。由此反映出SVR 方法能够自适应地选择位移最突出的边界点来表示边界运动。

图7 不同密度网格上各类核函数的控制点分布Fig. 7 Distribution of control points for various kernel functions on different density grids

图8 进一步给出了不同核函数拟合的均方差，说明所有被测核函数均能有效拟合边界的运动规律，但高斯核函数的拟合能力在密度最大的网格上有所下降，说明其健壮性较差。同时注意到 C PC4和 C PC6相比紧支撑组的其他核函数，其特性更接近全支撑函数，且自身阶数较高，较难体现紧支撑核函数高计算效率的优势。

图8 不同密度网格上各类核函数拟合的均方差Fig. 8 MSE of various kernel functions fitting on different density grids

综合以上分析，针对SVR 动网格方法，全支撑核函 数IQB、IMQB 和 紧 支 撑 核 函 数 CPC0、 CPC2、CTPSC0在建模效率和健壮性方面都要明显优于高斯核函数，且 CPC2核函数的效率最高。这些核函数对网格质量的影响将在第3 节做进一步评估。

3 案例应用

本节通过三个典型案例测试分析核函数及其参数设置对SVR 动网格方法的全面影响。对于变形运动后的网格质量，本文采用文献[12]提出的方法对三角形/四面体单元进行评估，并统计网格平均质量和最差单元质量。该方法通过函数fss(0 ≤fss≤1)综合相对尺寸度量和形状斜交度量评估变形后的网格单元，当函数值为0 时说明单元已退化、为1 时说明单元为正多边形/正多面体且其面积/体积保持不变，且对应函数值越大说明单元质量越高。

3.1 二维矩形块旋转平移

本案例主要考察在较大形变运动中采用SVR 方法进行多步变形后的网格质量以及核函数的影响。初始网格如图9 所示，其中外边界尺寸为 25×25，正中的矩形边界尺寸为 5×1。实验设置矩形块经由多步沿x和y轴负方向平移5 个单位长度，并绕中心逆时针旋转 60°。考虑网格变形程度较为严重，设置较小的回归阈值系数 λ =0.4。

图9 二维矩形块旋转平移案例的初始网格Fig. 9 Initial mesh for the 2D case of rectangular block rotation and translation

在进行实际动边界问题数值模拟时，从初始到最终状态往往需要计算多个时间步，因此也需要经过多步将初始网格变形到最终状态。首先设总的变形步数为20，研究不同核参数对网格变形后质量的影响。图10 给出了最终网格的最差单元质量度量值，从中可知 CPC0和 CTPSC0函数的变形质量始终较差，予以排除。三组核函数各自的核参数可行范围并不类似，其中高斯核可行取值范围最小，全支撑组的范围最大，此外 CPC2函数的最差质量在可行核参数范围内最稳定。

图10 分组对比不同核参数对最差单元质量的影响Fig. 10 Group comparison of the worst cell quality by different kernel parameters

为考察网格整体质量，测试了不同步数下变形到最终状态的网格单元最差质量和平均质量，结果如图11 所示。从中可知 CPC2函数和所有全支撑核函数的最差质量均随着变形步数增加而增大，且超过一定步数后质量趋于稳定。此外， CPC2函数的最差质量值是所有函数最高的，且在步数超过4 和8 后，IMQB 和IQB 函数的最差质量也分别超过了对照的高斯核函数。所有核函数变形后的网格平均质量度量值都较高，但需注意由于 CPC0和 CTPSC0等函数未能很好地将边界运动扩散至周边网格，导致矩形块附近的网格单元质量差而其余单元受影响较小保持了原有高质量，所以反而其网格平均质量更高。

图11 不同核函数下多步变形后最终的网格质量对比Fig. 11 Comparison of the final mesh quality after multi-step deformation by different kernel functions

图12 直观展示了采用IMQB 和高斯核函数通过4 步变形到最终状态的网格及其质量度量分布。在矩形块经过大幅度旋转平移后，其周围的网格单元仍与之保持着较好的连接关系，且采用IMQB 核函数变形后的网格质量明显更优。

图12 采用4 步变形的最终网格Fig. 12 Final mesh deformed after four steps by IMQB and Gaussian kernel function

3.2 二维多段翼型旋转

本案例通过多段翼的多体相对运动，进一步研究不同密度的网格变形质量受核函数的影响。初始网格的远场边界为圆形，半径长为25L，其中L为尾翼弦长，令内部边界中三段翼型的尾翼沿顺时针方向旋转1 5°，旋转中心为尾翼根部，其余部件保持静止。两套不同密度的初始网格如图13 所示，其中尾翼与最近的边界距离仅约0.1L。低密度网格的单元数为13 120，高密度网格的单元数为131 772。

图13 两套三段翼初始网格（局部）Fig. 13 Two initial grids of the three-segment wing (partial in detail)

本案例网格变形较缓和，取回归阈值系数 λ=1。此外，考虑到动静边界距离很近，应当取较小的核参数。结合前面实验的分析，设置高斯核宽度 σ为0.125L，全支撑组的核参数r为150L，紧支撑组的支撑半径R为5L。在两套网格上应用相同的核参数，不同变形步数下最终网格的最差单元质量如图14 所示。首先可以发现全支撑组函数和 CPC2的变形能力明显优于高斯核函数，且紧支撑组的其余两个函数的表现明显好于之前案例，可以推断是由于动静边界距离较小的缘故。

图14 两套网格变形后的最差单元质量Fig. 14 Worst cell quality after deformation for two grids

低密度和高密度网格的初始最差单元质量度量值分别为0.71 和0.69，对比两套网格变形后的质量，可知在同一核参数下，网格密度的增加令高斯核函数的变形能力明显下降，但对其他核函数影响不大。图15 分别展示了两套网格最终状态下动边界附近的网格单元质量分布，结果表明采用 CPC2和IMQB 函数进行变形后的网格质量明显优于高斯核，且质量较差单元也更少。综上所述，在多体相对运动的案例中，全支撑组和紧支撑组的核函数比对照的高斯核函数具有更好的变形能力和参数鲁棒性。

图15 三段翼动边界附近的最终网格质量分布Fig. 15 Final mesh quality distribution near the moving boundary of the three-segment wing

3.3 三维返回舱俯仰振荡

为进一步测试SVR 方法在三维动网格上的表现，本案例对三维周期性运动后的网格质量稳定性进行研究。初始的四面体网格如图16 所示，计算域的外边界成圆柱形，底部半径为5L，高约16L，L为返回舱中轴线长，返回舱最大半径为0.25L。共有61 514个边界点，303 927 个内部点。

图16 返回舱初始网格Fig. 16 Initial mesh of the return capsule

图17 给出了采用本文方法变形返回舱一个周期的运动示意图，其俯仰中心位于中轴线尾部五分之一处，并在y轴方向以 20°的振幅做往复运动。为测试网格质量的稳定性，令返回舱运动3 个周期，分120 步完成，每10 步完成1/4 个周期并运动至最大俯仰角。

本案例动边界结构复杂，为减小潜在的插值累积误差，取较小的回归阈值系数 λ=0.1，且由于是做周期运动，计算回归阈值时应取运动到最大形变时的步数计算，即s=10。此外，由于动静边界距离大于运动部件长度L的10 倍，应取较大的核参数，这里将高斯核宽度 σ设为1.5L，全支撑组的核参数r为1 200L，紧支撑组的支撑半径R为5L。

图18 给出了每一变形步后四面体网格单元的最差质量和平均质量。当达到最大俯仰角时，高斯核函数相应的网格质量略好于其他核函数，但相差很小且网格质量均较好。另一方面从图17 也可直观看出运动过程中保持了良好的网格质量。此外，之前实验表现最好的全支撑函数和 CPC2函数，其对应网格质量能很好地跟随返回舱周期运动而表现出周期性变化。不仅各周期对应变形步质量相当，且在三个运动周期后最差网格质量均能回到与初始几乎相等的水平，说明在整个变形过程中网格质量并未衰退。

图17 返回舱的周期俯仰运动Fig. 17 Periodic pitching motion of the return capsule

图18 返回舱各变形步完成后的网格质量Fig. 18 Mesh quality after each deformation step of the return capsule

表4 进一步给出了采用不同核函数进行变形的总时间开销，并以高斯核函数为基准列出了其他函数所节省的时间比例。其中采用 CPC2核函数的计算效率比高斯核提高了68%，其余全支撑函数也能节省超过40%的时间。因此，在大规模网格变形时，应用 CPC2、IMQB 和IQB 核函数具有更好的性能优势。

表4 不同核函数针对返回舱案例的时间开销Table 4 Time cost of different kernel functions for the return capsule case

4 结 论

本文通过进一步放宽贪心选点RBF 插值法的约束条件，提出了一种更高效的基于支持向量回归的机器学习动网格技术，并以前期研究采用的高斯核函数为基准，总结了所提方法的核参数设置规律。对于常见运动形式，一般可将 CPC2函数支撑半径设为3L～5L，IMQB 函数缩放系数设为300L～500L，其中L表示运动部件的特征长度。

在包含网格运动的大规模数值模拟中，动网格方法的计算效率始终是影响模拟周期的重要因素。在下一步的工作中，可以尝试将本文算法进行不同模式的并行化改造，使之适用于大规模复杂问题的工程实践。