基于導向濾波的流固耦合邊界數據插值分析

喻 敏 王志鴻 許文輝 章 軒

(武漢理工大學交通學院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)對數據進行導向濾波處理后，邊界數據的徑向基函數插值精度明顯提高，且差值精度受邊界數據密度的影響較小.

綜上，利用導向濾波對數據進行去噪后插值，可提高徑向基函數插值精度，有利于實現高效的流固耦合邊界數據傳遞.