耦合粒子圖像測速誤差的連續伴隨數據同化技術

鄧志文，何創新，劉應征,*

1. 上海交通大學 機械與動力工程學院 動力機械與工程教育部重點實驗室，上海 200240

2. 上海交通大學 燃氣輪機研究院，上海 200240

近年來，隨著實驗流體力學與計算流體力學的迅猛發展，人們對湍流現象的認知和了解亦不斷深入。常規認識湍流的途徑主要是通過實驗測量和數值模擬兩種手段。然而，單純依賴實驗測量或數值模擬的手段已無法滿足學界與工業界日益深入的研究需求。這和兩者各自的局限問題是分不開的。實驗測量方面，無論是傳統的熱線風速儀(Hot-wire Anemometer)，還是先進的粒子圖像測速(Particle Image Velocimetry，PIV)技術，依然受限于實驗測量成本高昂，測量區域有限，時空分辨率較低，實驗數據類型單一，多物理同步測量難度大以及實驗測量誤差難避免。數值模擬方面，高精度的直接數值模擬(Direct Numerical Simulation, DNS)和大渦模擬(Large-Eddy Simulation, LES)需要消耗大量的計算資源，而計算效率高的雷諾時均(Reynolds-Averaged Navier-Stokes, RANS)模型大多基于半經驗半理論公式，受限于模型結構形成缺陷和經驗常數難精確給出，對復雜流場的模擬精度不高；此外，數值模擬的初始條件或邊界條件在實際工程計算中常常難以準確獲得。

為了彌補實驗測量與數值模擬各自的不足，實現對高精度湍流流場的精準描述，數據同化(Data Assimilation, DA)作為一種深度融合實驗數據和數值模擬的技術，近年來愈發受到國內外流體力學領域研究學者的關注。數據同化，最早應用于氣象領域，隨后拓展到水文、海洋及地質領域。數據同化的三要素有實驗觀測、系統模型和數據同化算法，其核心思想在于利用實驗數據對系統模型進行結構上或參數上的調整修正，進而優化系統模型，提升預測性能。

數據同化算法主要可以分為兩類：① 基于數據統計算法的數據同化，如卡爾曼濾波、集合卡爾曼濾波(Ensemble Kalman Filter, EnKF)等；②基于變分算法的數據同化，如3DVar (Three Dimensional Variation)、4DVar等。Kato等通過采用EnKF及其變種的同化算法，優化了翼型繞流的來流攻角、馬赫數以及湍流渦黏分布，成功預測了翼型表面壓力系數分布。本文作者團隊用EnKF同化算法優化了不同RANS模型的模型經驗參數，成功預測了自由射流速度場。然而，大部分濾波算法主要針對湍流模型的邊界條件、初始條件、模型經驗參數進行修正，因而同化效果有限。相對而言，變分算法則通過添加源項的方式優化模型結構，設置與實驗觀測相關的目標損失函數，推導伴隨變分方程，通過求解目標損失函數極小值，便可修正優化湍流模型預測路徑。He等在Spalart-Allmaras(SA)模型中添加了模型修正因子,通過連續伴隨算法同化得到的分布，進而提升了SA模型對流動分離的預測能力。Foures和Symon等使用連續伴隨方法直接對雷諾應力項進行同化，準確還原了翼型繞流的平均流場。He等繼而提出了基于連續伴隨的順序數據同化算法，結合了順序數據同化與變分數據同化的優勢：既保留了變分算法的精度，又采用順序同化的方式實時動態耦合實驗觀測值修正模型預測結果；這種方式避免了變分算法在時間上正向逆向的積分，算法計算效率更高。值得提出的是，連續伴隨算法在實驗觀測準確可靠時，具有很好的同化效果，而實驗測量誤差這一重要因素如何在同化算法中準確考慮，無疑是很重要的。

本文在He等的研究基礎上，從PIV互相關原理的不確定度出發，提出了耦合PIV實驗誤差的連續伴隨同化算法，優化了原變分算法的目標損失函數以增強該同化算法的魯棒性。為了驗證算法的可行性，本文先對已知PIV流場植入合成誤差，進行同化對比測試，繼而對PIV互相關算法的不同參數設置所獲得的流場進行同化研究。結果表明，相比原變分算法，耦合PIV實驗誤差的變分算法能夠對實驗觀測數據去偽存真，抗誤差干擾能力明顯提升，魯棒性更強，實現對湍流流場的準確預測。

1 算法模型原理

1.1 算法框架設計

圖1 耦合PIV實驗誤差的連續伴隨算法流程Fig.1 Flow chart of continuous-adjoint algorithm coupled with PIV error

1.2 PIV實驗誤差分析

為了定量表征PIV實驗誤差的分布，本文提出了全場各點實驗權重分布函數。本節將從粒子圖像測速的互相關算法原理出發，詳細闡述實驗權重分布函數的提取實現過程。

粒子圖像測速技術是通過在流場中布撒示蹤粒子跟隨流體運動，利用激光加強粒子反射，經由相機在極短的時間內雙次或多次曝光，拍攝記錄粒子的運動軌跡，分析相鄰兩幀的粒子圖像從而獲得流場的速度信息。其中，流動速度矢量的獲取通常采用圖像互相關算法。該算法是假設在相鄰兩幀圖像中，被測流體的速度保持不變，并且粒子之間的相對位置在空間上保持一致，然后將粒子圖像劃分成若干個子區域(即問詢區域)，通過逐個區域粒子群的互相關計算，獲得其與臨近區域的相關性。傳統相關性的計算表達式為

=?(,)·(+Δ,+Δ)dd

(1)

式中：(,)和(,)分別對應兩幀圖像某一問詢區域的灰度分布函數；Δ、Δ分別為粒子在、方向的移動距離。如圖2(a)所示，當Δ、Δ等于粒子在、方向的真實位移時，互相關計算平面就會出現一個峰值，表明此處的相關性最大，可作為該區域粒子的平均運動位移，將其除以運動時間即可獲得速度矢量，即該子區域內粒子群的平均速度。為了節省計算資源，提高運算效率，一般采用基于傅里葉變換的互相關函數計算，如式(2)～式(4)所示:

圖2 PIV實驗誤差分析Fig.2 Procedure of PIV error analysis

(,)=((,))

(2)

(,)=((,))

(3)

(4)

如上所述，通過式(2)～式(4)計算便可獲得互相關函數的分布，如圖2(b)所示。經由圖像分析可知，該相關性分布總體呈現出二維高斯分布特征。因而，本文采用二維高斯分布函數對各個問詢區域的互相關函數進行獨立的擬合表達(～(,Δ,Δ,,))。設其概率密度函數(,)為

(5)

式中:為高斯分布的峰值大小;、為、方向所對應的位移方差。和的數值大小與PIV實驗的誤差密切相關，、數值越大表明此處的實驗值不確定度越高，誤差越大，實驗數值越不可信。

為了更好地擬合高斯分布規律，排除虛假相關干擾的影響，本文設置截斷濾波，濾波閥值為0.1倍的峰值，即相關性小于濾波閥值的點將置于0，如圖2(c)所示。在經過濾波處理過后，采用最小二乘法對高斯分布中的未知參數進行擬合，擬合結果如圖2(d)所示。最終分析對比高斯分布重構的互相關函數分布(圖2(e))與原始粒子圖像測算的互相關分布(圖2(a))情況，可以發現排除虛假相關干擾后，二者分布規律基本一致，采用此算法可以很好地重構的分布規律。

(6)

(7)

1.3 連續伴隨同化算法原理

本文的研究目標主要是針對平均場的數據同化，因此采用的實驗觀測為PIV平均場數據，同化的模擬結果為穩態結果。本文所采用的連續伴隨同化算法是由何創新等在Foures等的基礎上進行改進的，如式(8)所示，算法的核心思想是在穩態Navier-Stokes方程中加入同化控制源項，通過優化源項的分布，進而不斷減小實驗觀測與模型預測之間的誤差。

(8)

(9)

式中：、、分別為速度矢量、壓力和流體密度；為等效黏度，即流體分子黏度與湍流各項同性渦黏部分之和，可以通過RANS模型直接計算獲得;為各項異性渦黏，可通過伴隨變分求解。原連續伴隨的目標損失函數為

(10)

式中：為實驗觀測速度；為無窮遠處的來流速度；代表數值模擬的計算區域；為單位無量綱轉換系數；為Navier-Stoke方程和連續性方程約束；作為掩模函數指定了實驗觀測所在的位置：當=1時，表明此處有實驗觀測約束，反之為0，則表明無實驗觀測。

為了求解目標損失函數的極小值，引入拉格朗日函數進行變分方程的推導求解：

(11)

式中：、為伴隨速度變量和伴隨壓力變量，經過變分推導可得伴隨變分方程組為

(12)

(13)

連續伴隨的數據同化過程，可以看作是求解目標損失函數的極小值過程。通過不斷優化迭代的分布，最終使得速度場既滿足實驗觀測規律又滿足Navier-Stokes方程以及變分方程的約束。平均場的同化計算迭代流程圖如圖3所示。圖中、為初始流場的速度和壓力；、分別為計算相對誤差和人為設定的收斂誤差值；為松弛因子，控制的迭代速度。本文通過引入實驗觀測權重分布函數修正模型的目標損失函數，如式(14)所示，從而將PIV各個速度矢量的實驗誤差與連續伴隨同化算法迭代計算相耦合，進而提升該算法的抗誤差干擾能力。

圖3 平均場的連續伴隨數據同化計算迭代流程Fig.3 Iteration procedure of continuous-adjoint based DA for mean flow field

(14)

2 實驗設計驗證

為了驗證所提同化算法的有效性和魯棒性，本文采用有限長平板繞流的PIV實驗數據作為同化過程中的實驗觀測輸入。對于此類分離再附流動，大多RANS模型難以精準預測，與實驗結果相差甚遠，因而需要數據同化技術進行優化。該實驗的原始數據來源于文獻[28]。

該實驗是在一個300 mm (高)×300 mm (寬)×2 000 mm (長)的低速開式風洞中進行的。為了消除風洞壁面對流動二維性的影響，平板厚度設計為=24 mm，弦厚比為=9，并垂直安裝于風洞測試段的正中間位置。本實驗中自由來流的流速保持在=10 m/s，依據平板厚度計算對應的雷諾數為=15 800。如圖4所示，為了捕捉平板繞流的分離與再附現象，實驗數據的觀測區域(Region of Observation，ROB)選定于平板正上方的流場(0≤≤9,0≤≤25)。本實驗采用Di-Ethyl-Hexyl-Sebacate (DEHS)油滴作為示蹤粒子，粒徑為1 μm左右，對實驗流體具有良好的跟隨性。激光平面由Nd:YAG激光器(135 mJ/pulse, 532 nm, 8 ns, Litron, UK)產生，厚度為1 mm。采用CCD相機(IPX 16M, IMPERX, USA)對粒子的運動圖像進行采集記錄。激光脈沖頻率與相機采樣速率由同步器控制，均保持在1 Hz。最終共采集6 000幅圖像，圖像分辨率為4 872 pixel×3 248 pixel。通過PIV互相關算法可計算獲得3 000 個速度矢量場，求平均可得PIV實驗的統計平均流場。關于此實驗更詳細的描述和進一步的結果可在文獻[28]中獲得。

圖4 實驗裝置示意圖Fig.4 Schematic diagram of experimental setup

3 結果分析討論

為了全面充分地驗證本文算法的優化效果，本文對比了原連續伴隨算法以及改進后的算法在不同誤差環境下的同化效果。首先對不同程度的合成誤差干擾進行了測試驗證，隨后進一步對不同PIV互相關算法設置引入的算法誤差進行同化效果測試。

圖5為PIV實驗粒子圖所獲得的全場速度矢量的誤差分布云圖，顏色越藍代表此處的互相關性越高，不確定度越低，即實驗誤差越小，實驗觀測結果可信度高。如圖5(a)所示，平板上方的主流流動區域顏色較藍，表明實驗誤差小，較為可信，因而實驗權重分布值較大。偏離藍色的區域主要集中于平板繞流的剪切層附近以及近壁面區域，這是因為剪切層附近的流動比較紊亂，而近壁區域存在壁面反射，因而實驗結果誤差較大。圖5(b) 是剪切層附近的局部細節放大圖，每一塊子區域都表示一個速度矢量的互相關函數分布情況，能夠通過1.2節所述方法進行高斯分布擬合。圖5(c)為經過誤差分析后轉換得到的實驗權重分布函數。對于實驗誤差相對較高的壁面和剪切層區域權重分配較小，而誤差較小的自由來流區域的權重分配較大。

圖5 鈍體平板繞流PIV實驗誤差分析Fig.5 Analysis of PIV error of flow around blunt plate

3.1 合成誤差對同化效果的影響

圖6為添加不同級別的合成誤差干擾后所得的流向速度平均場()云圖對比。圖6(a)為采用最優的PIV算法參數設置所得的統計平均場結果，亦作為本次同化實驗中的真實速度。合成誤差的干擾程度主要由參數控制，如式(15)所示。

(15)

圖6(b)～圖6(c)中為添加不同程度的實驗干擾后的平均速度云圖，對應=01,02,05的情況。合成誤差的添加規律與PIV實驗誤差分布相關，即在壁面及流動剪切區域的合成誤差較高，同時由圖6(d)可知當=05時，合成誤差已對實驗結果產生了較大的影響，因而本研究中的最大值為0.5。將帶有合成誤差的PIV結果(圖6(b)～圖6(d))作為實驗觀測輸入進行同化對比測試，最終同化效果如圖7～圖9所示。通過對比研究發現，當采用帶有合成誤差的實驗結果進行同化時，原連續伴隨同化模型(如圖7(a)～圖7(c) Adjoint結果，圖8中藍色實線)為了盡可能地同化帶有誤差的實驗值(對應圖8 中PIV+結果)，出現了極不合理的預測結果，違背了真實的流動分布規律(對應圖8 中Ref結果)。尤其在近壁面的速度容易被該區域內誤差較大的速度分布所干擾(圖8中綠色圓點)，進而導致整體同化效果很差，流動計算發散。而當采用耦合PIV實驗誤差的同化模型時(如圖7(a)～圖7(c) Adjoint+結果，圖8紅色實線)，則可降低合成誤差的干擾，提升計算收斂的穩定性，同時相比傳統的SA模型預測結果有很大的改善，對流動剪切的預測更為精準，與真實的流動速度分布曲線吻合良好。

圖6 添加不同合成誤差所對應的流向速度云圖Fig.6 Contours of streamwise velocity adding different levels of synthetic error

圖7 不同合成誤差下采用不同同化方法的結果對比云圖Fig.7 Contour plots of assimilated results with different synthetic errors using different data assimilation methods

圖8 不同位置處流向速度對比曲線(α=0.2)Fig.8 Transverse profiles of streamwise velocity at different locations (α=0.2)

圖9為平板流動再附點位置以及流線圖對比結果。從圖中可以看出，耦合PIV實驗誤差的數據同化模型(圖9(b)～圖9(d))與PIV結果吻合良好，對分離再附點的位置預測結果準確，且即使在高合成誤差場景下(圖9(d)，=05)，依然有很高的預測精度。

圖9 不同合成誤差下流場流線圖以及再附點位置分布Fig.9 Streamlines and location distribution of reattachment point obtained using data assimilation method with different synthetic errors

3.2 PIV算法參數設置對同化效果的影響

為了進一步驗證所提同化算法的有效性和魯棒性，本文采用不同的PIV算法設置所獲得的流場作為同化的觀測輸入進行研究。采用不同的PIV互相關算法參數設置對同一流場的粒子圖進行分析所得的計算結果也有所不同。如當問詢窗口的尺寸發生變化時，PIV計算的速度矢量數量及其對應的誤差大小也隨之發生改變。為了便于理解和討論，本文中將這部分僅由不同PIV算法參數設置所引起的誤差，簡稱為“PIV算法誤差”。

如圖10所示，當不對原粒子圖進行任何預處理(諸如圖像對比度增強以及噪聲濾波算法等)，并采用不同的問詢窗口設置時(24×24,32×32,48×48)，PIV算法的計算結果會產生較大的計算誤差，尤其在近壁面區域。圖10白色圓圈內的流動與真實的流動現象并不相符，分離泡對壁面的依附出現失真，無法準確捕捉再附點位置。

圖10 采用不同問詢窗口計算所得的流向速度云圖Fig.10 Contours of streamwise velocity calculated using different sizes of interrogation windows

圖11為當采用不同問詢窗口時，實驗權重的分布函數云圖。如圖所示，針對本文的驗證案例，當采用較大尺寸的問詢窗口(48×48)，PIV實驗結果的誤差會有所減小，相應地，其所對應的實驗權重分布函數的權重值會有所增大(圖11(c))。而當問詢窗口尺寸較小時(24×24)，所對應的實驗權重值較小。觀察實驗權重的分布可得，近壁面及剪切流動區域的實驗誤差較大，實驗權重取值小。當問詢窗口較小時(圖11(a)和圖11(b))，近壁及剪切區域的權值分布接近于0，證明此處的實驗結果誤差很大，實驗觀測不可信。

圖11 采用不同問詢窗口所得的實驗權重分布函數云圖Fig.11 Contours of experimental weighting function distribution obtained using different sizes of interrogation windows

將不同PIV算法所得的實驗結果(圖10(b)～圖10(d))作為實驗觀測進行同化測試，最終同化結果如圖12～圖14所示。圖12為流向速度的平均場云圖，由圖可得，相比于實驗結果(圖10 (b)～圖10 (d)，圖13綠色圓點)，兩種數據同化算法計算收斂較好，未觀察到分離泡依附失真現象。然而，仔細對比無誤差引入的連續伴隨模型(圖12(a)～圖12(c)中Adjoint結果，圖13 藍色實線)與耦合PIV誤差的同化模型(圖12(a)～圖12(c) 中Adjoint+結果，圖13紅色實線)，可以發現，原伴隨變分模型的預測結果在剪切層附近存在速度變化不光滑不連續現象，且對近壁面的流動結果預測并不準確。而耦合實驗誤差修正的同化模型所得的流場速度曲線連續光滑，不存在較大的偏離值，與真實速度曲線分布(黑色圓點)具有良好的一致性。圖14為平板流動再附點位置以及流線圖對比結果，從圖中可以看出，由于采用不同的PIV互相關算法參數設置，引入了除實驗誤差影響外的算法誤差，導致流場的再附點位置存在較大的偏差(圖14(a)～圖14(c)中PIV計算結果)，與真值4.43(圖9(a))相去甚遠。采用原伴隨變分同化算法所得再附點的位置(圖14(a)～圖14 (c)中Adjoint結果)會貼合含有實驗誤差和PIV算法誤差的觀測數據(圖14(a)～圖14(c)PIV計算結果)，并無明顯的改善優化。而采用耦合PIV誤差的同化模型時(圖14(a)～圖14(c)中Adjoint+結果)，同化結果受PIV算法誤差干擾較小，其對再附點的位置預測更為準確，并且對流動模式預測良好，同化效果更好。

圖12 采用不同同化模型在不同PIV算法誤差下的同化結果云圖Fig.12 Contours of assimilated results with different errors induced by different PIV algorithms using different data assimilation methods

圖13 不同位置處流向速度對比曲線 (問詢窗口：32×32)Fig.13 Transverse profiles of streamwise velocity at different locations (interrogation window: 32×32)

圖14 不同問詢窗口設置下PIV和同化結果的流場流線圖及再附點位置分布Fig.14 Streamlines and location distribution of reattachment points determined by PIV and assimilated results obtained using different settings of interrogation windows

4 結 論

本文從粒子圖像互相關的基本原理出發，推導獲得反映PIV實驗誤差分布規律的實驗權重分布函數，通過將其與連續伴隨算法耦合，優化同化算法的目標損失函數以提升算法的抗誤差干擾能力。為了驗證所提算法的有效性和魯棒性，本文先對已知PIV流場植入合成誤差進行同化對比測試，繼而對PIV互相關算法的不同參數設置所獲得的流場進行對比研究，可以得到以下結論：

1) 從粒子圖像互相關算法推導得到的PIV實驗權重分布函數能定量反映PIV實驗誤差分布規律，能與連續伴隨同化算法緊密耦合。

2) 在合成誤差同化測試中，原伴隨變分模型易受合成誤差干擾，無法進行有效地同化，預測結果有悖物理規律。而耦合PIV實驗誤差的同化模型能夠顯著地降低合成誤差的干擾，提升求解收斂的穩定性，準確預測流場的真實信息。

3) 在由PIV互相關算法的不同參數設置所獲得的流場同化測試中，原伴隨變分模型無法精準預測流場的細節，而耦合PIV實驗誤差的同化模型可以降低PIV算法誤差的影響，對流場的預測結果較為準確，貼近流場真實分布規律，細節還原度高。

在后續的研究工作中，本文將進一步對所提算法的魯棒性和普適性進行探究，如針對不同問詢窗口、粒子濃度和粒子直徑等多參數耦合作用下的PIV流場進行同化對比測試。同時開展對瞬態流場的同化對比研究，分析速度脈動、雷諾應力等統計量的同化精度。