基于导向滤波的流固耦合边界数据插值分析

喻 敏 王志鸿 许文辉 章 轩

(武汉理工大学交通学院1) 武汉 430063) (中国舰船研究设计中心2) 武汉 430063)

0 引 言

流固耦合(fluid-structure interaction,FSI)是描述流体动力学定律和结构力学定律之间的多物理场耦合.流固耦合的关键是实现精确的耦合分析，其中分区耦合较适用于数值计算，在CFD数值计算中应用广泛.分区耦合涉及流体到结构的网格节点插值、荷载和变形传输问题，即边界数据传递.利用径向基函数(radial basis function,RBF)插值方法进行数据传递[1]，可有效解决边界网格不匹配而导致节点数据难以直接传递的问题[2]，特别是对复杂非线性流固耦合问题具有很好的适用性，具有良好的工程应用前景.

针对非线性流固耦合问题，基于黏性理论数值求解流体场时，流体网格节点具有较强非线性，使得流体网格节点很容易出现数值噪声(包括异常值).对含有数值噪声的流场网格节点进行径向基函数插值，可能导致插值结果的失真.近年来，关于边界数据传递的径向基函数插值研究主要集中在径向基核函数和算法效率方面[3-4].

常用的数值噪声处理方法有统计分析法、特征提取法和密度聚类法等[5-6].统计分析法，适用于噪声的判断，噪声处理需配合其他算法；特征提取法，如小波去噪，则过分依赖于阈值的选取；密度聚类算法，如DBSCAN等，计算复杂度高，对输入参数敏感.上述方法均不适用于流固耦合边界数据的去噪计算，而图像滤波方法无需先验知识，原理简单，应用广泛、灵活，在三维数据去噪方面效果明显[7-10].文中尝试将图像处理领域的滤波算法应用于CFD流场数据处理，分析高斯滤波、双边滤波和导向滤波三种滤波方法的去噪效果，确定合适的滤波方法，对CFD计算数据进行去噪处理，并将其运用到流固耦合边界数据的插值计算中，进而形成基于径向基函数插值的边界数据传递方法.

1 边界数据的径向基函数插值

1.1 流固耦合数据传递矩阵

在流固耦合问题的流体动力学计算中，由于结构网格与流场网格的差异，需要在流固耦合边界上进行数据传递.流固两相交界处一般满足虚功相等原理，即：

(1)

式中：W为虚功；uf和us分别为流固边界处的虚位移；ff和fs分别为流体表面压力和固体表面压力。以H表示传递矩阵，经推导，可以得到：

(2)

由于流固耦合数据传递是双向的过程，在此为了简化叙述，仅以流体区域信息传递到固体结构区域为例。假定流固耦合边界处两侧分别存在nf(欧拉点数)和ns(拉格朗日点数)个节点，以Zf和Zs分别表示流体和固体节点坐标，数据传递中引入临时传递变量Uf，则径向基函数插值的基本形式为

(3)

(4)

式中：λ、α均为径向基函数中的系数矩阵，式中径向基函数φ的变量是欧拉点(即流体域上的点)之间的距离.则待求系数：

(5)

对于流体传递给固体的情况，有临时传递变量：

(5)

(6)

式中：径向基函数φ的变量是拉格朗点(即固体结构域上的点)与欧拉点(即流体域上的点)之间的距离。由以上公式可得：

(8)

由此，流体区域数据交换到固体结构区域的传递矩阵为

(9)

同理，固体区域数据交换到流体区域的传递矩阵为

(10)

由式(9)～式(10)，结合本节中讨论的数据交换能量守恒，可实现流固耦合边界上速度、压力等高效精确的传递.从推导过程可知，径向基插值对边界数据传递十分必要.

1.2 径向基函数插值

在流固耦合边界数据传递过程中，基于径向基函数的插值不依赖网络节点的拓扑结构，适用于大量散乱数据问题的求解.理论上，径向基函数表示为一类以欧式距离为变量的函数集合，其数学表达式为

(11)

式中：f(x),x∈R为给定高维空间的未知函数；λi为待求的加权系数；Φ(x)为径向基函数；模‖x-xi‖为插值点x到被插值点xi的欧式距离；n为插值点个数.对条件正定的径向基函数，增加一个多项式，以确保方程解唯一.多项式pi(x)需满足如下约束条件.

(12)

由式(11)～式(12)，可得线性方程组：

(13)

有ΘW=Ω，因此为求系数矩阵W，需要求解系数矩阵Θ-1，若插值点x存在噪声，系数矩阵Θ易出现变态，导致径向基函数插值结果的失真，甚至无法求解.因此，有必要对插值点进行去噪处理，以提高径向基函数插值的稳健性.

2 滤波算法

2.1 高斯滤波

一维高斯函数为

(14)

高斯滤波使用的核函数为x和y两个方向上一维高斯的乘积，且两个维度上的标准差σ相同，其中心点始终为滤波窗口中心坐标，故取均值μ=0，因此，滤波所用到的核函数为

(15)

高斯滤波器宽度是由参数σ控制的，σ越大，高斯滤波器的频带就越宽，平滑程度就越大，因此可以通过调节参数σ的取值，在过平滑和欠平滑之间取得折衷.然而，由于高斯滤波的高斯核只考虑了数据间的空间距离，没有考虑数据取值上的差异，因而容易丢失数据的细节信息.

2.2 双边滤波

双边滤波是一种非线性的滤波方法，综合了高斯滤波器和α-截尾均值滤波器的特点，同时考虑滤波数据的空间信息及其数值上的相似性，以达到保边去噪的目的，一定程度上弥补了高斯滤波的缺陷.双边滤波的核函数为

w(i,j,k,l)=

(16)

式中：(i,j)为滤波窗口中数据的坐标；(k,l)为滤波窗口中心点坐标；I为对应数据的大小；σd、σr均为高斯函数的标准差，但在数值上一般不等.

则w(i,j,k,l)=d(i,j,k,l)×r(i,j,k,l)，即双边滤波的核函数为两个子函数的乘积，d(i,j,k,l)代表数据的空间域信息，决定滤波器模板系数的空间邻近因子，其大小随着坐标点与中心点之间距离的增加而减小；r(i,j,k,l)代表数据取值的范围域信息，决定滤波器模板系数的数据数值相似因子，其大小随着数据之差的增大而减小.

现考虑两种极限，①滤波窗口中数据与中心数据基本相同，邻域内各数据的差趋近于0，有相似因子r(i,j,k,l)≈1，此时空间域权重起主要作用，相当于对数据进行高斯滤波；②窗口中其他数据与中心数据相差极大的情况下，相似因子r(i,j,k,l)→0，此时双边滤波的核函数w(i,j,k,l)≈0，不对数据进行滤波操作，从而达到保留数据细节信息的目的.也正因此，当出现高强度噪声时，滤波算法的范围域核函数权值的计算将受到影响，导致其去噪能力下降，甚至出现“梯度反转”现象[12].

2.3 导向滤波

根据文献[10]，导向滤波需满足以下条件，

1)输出数据Q与原始数据P尽可能的相似，即(Q|P)取得最小.

2)输出数据Q与引导数据I的梯度相似，有Q=aI，即满足Q=aI+b.

现假设在滤波窗口ωk内输出数据Q与引导数据I满足线性关系，有qi=akIi+bk，∀i∈ωk，其中：k表示窗口ωk的中心位置，ak、bk在窗口中为常数.为满足第一个条件，通过求解该线性函数的系数，使得输出值q与输入值p之间差值最小化，即存在

上式可通过最小二乘求解，可得：

(17)

(18)

(19)

(20)

由此可见，在ε≠0的情况下，若引导数据I和输入数据P相同，导向滤波可在处理数据的同时保留数据的细节信息.此外，由于输出数据与引导数据梯度相似，滤波结果不会出现双边滤波的“梯度反转”问题.

2.4 仿真验算分析

为了比较上述三种滤波算法的滤波效果，采用二元函数z=xexp(-x2-y2)进行验算.对函数进行加噪，结果见图1a)，利用三种滤波算法对其进行去噪，比较滤波前后径向基函数插值误差，由此选取最适合的滤波方法对CFD数据进行去噪.

仿真过程中，各滤波器参数的取值经过多次验算获得：高斯滤波使用5×5窗口、标准差σ=1.5的核函数对数据进行滤波处理，结果见图1b).

双边滤波同样取5×5窗口，核函数中取σd=1.5，σr=0.9.此外，由于双边滤波要求输入数据P∈[0,1]，需要对数据进行预处理，预处理过程如下：

1)记Pmin为输入数据P最小值，将数据P中所有元素加上|Pmin|，得新数据P′=P+|Pmin|，记Pmax为其最大值。

2)对数据P′做归一化处理，得最终输入双边滤波器的输入数据P″.

同理，对滤波后数据进行恢复，最终得到双边滤波的结果见图1c).

导向滤波将原始数据同时作为输入数据和引导数据进行滤波，同样选取5×5的滤波窗口，参数ε取1.5，滤波结果见图1d).

图1 加噪数据与三种滤波方法所得数据对比

由图1可知，高斯滤波对边缘数据野值的处理效果最差，且与双边滤波一样，在梯度变化大的区域(如滤波数据中心部分)滤波效果差，后者主要由“梯度反转”现象所致，而导向滤波并没有出现此类现象，在处理此类数据时滤波效果较为理想.从算法计算效率看，见表1。高斯滤波计算时间最短，这得益于其简单的计算原理；同时，尽管导向滤波和双边滤波均考虑了数据细节信息的特性，导向滤波计算所用时间仍在可接受范围内，而双边滤波颇为耗时，不适用于处理大数据，如流场数值计算.在平均绝对误差(MAE)上，导向滤波数值最大，其原因在于为比较三种算法的滤波效果，而选取了同样大小的窗口，多次仿真计算可知，减小滤波窗口可进一步提高导向滤波的滤波效果，如导向滤波在3×3窗口中MAE=1.894，而高斯滤波和双边滤波已在5×5窗口取得最优滤波效果.

表1 三种滤波方法计算效率对比

将滤波前后的数据分别进行径向基函数插值计算，比较插值精度.插值核函数采用多元二次核函数，径向基函数插值平均绝对误差结果见表2。噪声对插值精度有明显影响；滤波后的插值精度都得到一定程度的提高，其中导向滤波后边界数据的插值效果得到大幅度改善，为三种滤波算法中的最优，因此本文将采用导向滤波对CFD数据进行滤波处理.

表2 五类数据插值误差对比

3 导向滤波流场数据处理应用

3.1 流场流速数据滤波效果对比

文中以CFD数值计算所得某一截面上的速度场数据为例进行滤波计算，见图2a)，速度场以矢量图的形式给出，图2b)为由矢量图转换得到的可视化图形，由此可见流场数据与图像像素值之间存在着相似性.基于这一方面的考虑，本文将导向滤波方法用于流场数据的去噪处理上，以提高径向基函数插值的精度.

文献[13]采用移动最小二乘算法对耦合前CFD数据进行去噪，取得了较好效果，并已应用于实际的自主CFD计算平台.下面将本文的研究对象——导向滤波方法，与移动最小二乘方法在流场数据去噪方面进行对比，结果见图3.结合图2a)速度场数据的变化，可见移动最小二乘算法虽能有效剔除流场数据存在的噪点，但其并不考虑邻域(滤波窗口)内数据的变化趋势，对连续变化物理场(如速度场)的描述显然不如导向滤波，而导向滤波所得结果更符合实际速度场的变化规律，这使得相较于传统的移动最小二乘算法，导向滤波方法在CFD数据去噪方面更具优势.

图2 实验数据及其模信息

图3 移动最小二乘与导向滤波去噪效果对比

3.2 流固耦合边界数据插值验证

为对比导向滤波前后边界数据的径向基函数插值精度，以2.4中导向滤波前后数据为例，分别在不同流体网格密度条件下，对3种密度的边界数据进行插值，以此分析导向滤波方法对流固耦合边界数据插值效果的影响。仿真中，根据流体网格密度分成2组：G1，G2，每组网格中又有着不同密度的待插值边界数据点，见图4，由此可得六种组合方式。为分析导向滤波方法在流固耦合边界数据插值精度上的影响，表3～4为导向滤波前后对应网格密度下边界数据的插值误差对比.

图4 不同边界点密度网格示意图

表3 滤波前数据在不同边界点密度下插值误差

表4 滤波后数据在不同边界点密度下插值误差

由表3～4可知，无论是对滤波前还是滤波后的边界数据进行插值，G2分组所得插值误差都明显大于G1分组。根据数据插值或拟合理论，插值误差和数据密度存在一定的关系，当填充距离很小，即数据密度很大时，插值误差也很小.

一般地，用填充距离h描述数据{xj}的密度：

(21)

可以证明[14]：

(22)

同时，在每一分组内部，不同的边界数据密度对应的插值误差十分接近，这是因为插值精度主要与插值数据密度有关，而与被插值点密度无关。另外，对比表3～4误差，滤波后数据的插值误差明显小于原始数据的插值误差，可知导向滤波处理有效提高了径向基函数插值的精度，对实现高效的流固耦合边界数据传递具有积极作用.

4 结 论

1)导向滤波既能保留数据的细节信息，又具有较高的计算效率，滤波效果较高斯滤波和双边滤波方法更优。

2)相较于移动最小二乘算法，导向滤波在处理流场数据的过程中，综合考虑了邻域内数据的变化，更适用于连续变化物理场数据的处理。

3)对数据进行导向滤波处理后，边界数据的径向基函数插值精度明显提高，且差值精度受边界数据密度的影响较小.

综上，利用导向滤波对数据进行去噪后插值，可提高径向基函数插值精度，有利于实现高效的流固耦合边界数据传递.