Gappy POD 方法重構(gòu)湍流數(shù)據(jù)的研究1)

李天一 Buzzicotti Michele Biferale Luca 萬(wàn)敏平 , 陳十一

* (南方科技大學(xué)力學(xué)與航空航天工程系,廣東深圳 518055)

? (羅馬第二大學(xué)物理系和意大利核物理研究所,意大利羅馬 00133)

引言

本征正交分解(proper orthogonal decomposition,POD)是一種統(tǒng)計(jì)分析、簡(jiǎn)化數(shù)據(jù)集的方法.對(duì)于一個(gè)函數(shù)系綜,例如在實(shí)驗(yàn)中獲得的數(shù)據(jù),POD 可以給出一組模態(tài)分解的基,具有很多應(yīng)用中所需要的良好的性質(zhì).其中,最引人注目的就是最優(yōu)性(optimality),即POD 給出了采用有限模態(tài)來(lái)刻畫(huà)一個(gè)無(wú)限維過(guò)程主要成分的最有效的方式.

POD 方法由Lumley 首次引入到湍流問(wèn)題的背景中[1].在大部分應(yīng)用中,POD 主要用來(lái)分析實(shí)驗(yàn)或計(jì)算數(shù)據(jù)以獲取其中的主要模態(tài)[2-5],即擬序結(jié)構(gòu)(coherent structures).此外,很多研究采用POD 提供的基函數(shù)構(gòu)成的低維子空間來(lái)對(duì)原始問(wèn)題建立降維模型[6-10].

Everson 和Sirovich[11]于1995 年首次提出了采用POD 來(lái)重構(gòu)存在缺失的(gappy)數(shù)據(jù)的方法,即gappy POD.此后很多學(xué)者針對(duì)不同的問(wèn)題使用并發(fā)展了gappy POD 重構(gòu)方法[12-15].Venturi 和Karniadakis[16]對(duì)于圓柱繞流的直接數(shù)值模擬數(shù)據(jù)進(jìn)行了gappy POD 重構(gòu)研究.通過(guò)引入適當(dāng)修正,他們提升了gappy POD 的魯棒性和精度.同時(shí),研究了流場(chǎng)缺失率大小對(duì)gappy POD 重構(gòu)的影響,并給出了相關(guān)解釋.Gunes 等[17]對(duì)比了gappy POD 和克里格(Kriging)插值對(duì)于非穩(wěn)態(tài)圓柱繞流數(shù)據(jù)的重構(gòu)效果,發(fā)現(xiàn)當(dāng)流場(chǎng)的時(shí)間分辨率不高或者空間缺失率較大時(shí),克里格插值比gappy POD 重構(gòu)更加有效.

在前人的工作中,以下兩個(gè)因素對(duì)gappy POD的影響均未被系統(tǒng)研究.第一,數(shù)據(jù)本身的復(fù)雜程度.當(dāng)雷諾數(shù)Re較低時(shí),流動(dòng)處于層流狀態(tài),此時(shí)描述流動(dòng)的POD 模態(tài)較少[16-17];當(dāng)Re很高時(shí),流動(dòng)為復(fù)雜的湍流狀態(tài),需要用來(lái)描述流動(dòng)的POD 模態(tài)往往很多.給定流場(chǎng)的分辨率和破損區(qū)域,流動(dòng)的POD模態(tài)數(shù)量越多,重構(gòu)時(shí)所需要確定的未知系數(shù)就越多.第二,破損區(qū)域的幾何形狀.很多關(guān)于gappy POD 的研究考慮了不同缺失率(缺失面積)對(duì)重構(gòu)效果的影響[15-16],但是其中的缺失區(qū)域均是隨機(jī)生成的,并沒(méi)有對(duì)其幾何形狀的約束.在破損區(qū)域面積相同時(shí),破損區(qū)域的幾何形狀不同會(huì)導(dǎo)致流場(chǎng)損失的相干信息不同,從而可能會(huì)對(duì)gappy POD 的重構(gòu)效果產(chǎn)生影響.

在本文中,選取旋轉(zhuǎn)湍流流場(chǎng)來(lái)進(jìn)行g(shù)appy POD 重構(gòu)的研究.旋轉(zhuǎn)湍流是一個(gè)具有豐富物理現(xiàn)象的湍流的典型范例.其中,通過(guò)外力輸入的能量不僅導(dǎo)致了大尺度的氣旋和反氣旋結(jié)構(gòu),也產(chǎn)生了小尺度的、間歇的、均勻各向同性的高度非高斯的擾動(dòng),這使得旋轉(zhuǎn)湍流表現(xiàn)為分布在不同空間尺度和時(shí)間尺度的混亂流動(dòng)[18-24].此外,旋轉(zhuǎn)湍流對(duì)于很多地球物理現(xiàn)象也有著重要的意義[25-29].

本文通過(guò)旋轉(zhuǎn)湍流場(chǎng)數(shù)據(jù),研究了不同流場(chǎng)復(fù)雜程度以及破損區(qū)域的面積大小和幾何形狀對(duì)gappy POD 重構(gòu)結(jié)果的影響.通過(guò)更加嚴(yán)格地表述gappy POD 重構(gòu)過(guò)程,其推導(dǎo)了gappy POD 重構(gòu)誤差的公式,并解釋了它與流場(chǎng)復(fù)雜程度以及破損區(qū)域的關(guān)系.本工作期待為gappy POD 的應(yīng)用和發(fā)展提供參考.

1 POD 理論與gappy POD 重構(gòu)

1.1 POD 簡(jiǎn)介

為了在1.2 節(jié)中系統(tǒng)、清楚地描述gappy POD重構(gòu)過(guò)程,本節(jié)簡(jiǎn)要介紹POD 分解的過(guò)程.在進(jìn)行g(shù)appy POD 重構(gòu)之前,需要由Nc個(gè)已知的完整訓(xùn)練數(shù)據(jù)來(lái)求得POD 模態(tài).

對(duì)于第c個(gè)訓(xùn)練數(shù)據(jù)uc(x) ,可計(jì)算其在位置x處相對(duì)所有訓(xùn)練數(shù)據(jù)系綜平均的偏移

1.2 Gappy POD 重構(gòu)

本節(jié)介紹gappy POD 的重構(gòu)過(guò)程,相比Everson和Sirovich[11]中的描述,這里進(jìn)一步在數(shù)學(xué)上嚴(yán)格分析了不同因素對(duì)重構(gòu)效果的影響,包括數(shù)據(jù)復(fù)雜程度、采用的POD 模態(tài)個(gè)數(shù)以及缺失區(qū)域.

假定考慮的流動(dòng)數(shù)據(jù)最多可以由Nflow個(gè)POD模態(tài)描述,那么Nflow的大小反映了流動(dòng)的自由度或復(fù)雜程度.當(dāng)僅用前Ncomp個(gè)POD 分量描述系統(tǒng)時(shí),對(duì)任意一個(gè)數(shù)據(jù)樣本 φ (x),有

式中的截?cái)鄽埐顁取決于Ncomp,為

一個(gè)存在部分缺失的樣本可以被表示為

其中在缺失的位置x處有m=0,而在其他位置處有m=1.對(duì)于該樣本,有

Everson 和Sirovich[11]中通過(guò)最小二乘擬合,即最小化如下誤差

2 旋轉(zhuǎn)湍流流場(chǎng)重構(gòu)

2.1 湍流數(shù)據(jù)簡(jiǎn)介和重構(gòu)問(wèn)題描述

本工作選用旋轉(zhuǎn)湍流流場(chǎng)來(lái)進(jìn)行g(shù)appy POD 重構(gòu)的研究.TURB-Rot 是一個(gè)公開(kāi)的旋轉(zhuǎn)湍流數(shù)據(jù)集[30],包含105 600 張垂直于旋轉(zhuǎn)方向的二維流場(chǎng)切片,切片的分辨率為 6 4×64,所有切片均來(lái)自同一個(gè)三維旋轉(zhuǎn)湍流的直接數(shù)值模擬.由于數(shù)據(jù)集中的旋轉(zhuǎn)湍流場(chǎng)高度二維化,所有切片均可看作來(lái)自同一系綜.數(shù)值模擬中采用了周期邊界條件,因此可對(duì)二維流場(chǎng)切片進(jìn)行如下傅里葉譜展開(kāi)

由于旋轉(zhuǎn)流場(chǎng)的各向異性,gappy POD 對(duì)u的不同分量的重構(gòu)效果并不相同.不失一般性,本文考慮速度模u=‖u‖2的gappy POD 重構(gòu).在由速度u計(jì)算速度模u時(shí),需要注意消除混淆誤差[31],使得u可以寫(xiě)為

式中kf為32.可以看到,kf的大小確定了u所包含的傅里葉模態(tài)個(gè)數(shù),從而確定了流場(chǎng)的復(fù)雜程度.因此可以通過(guò)僅保留 ‖k‖＜kf=8 的模態(tài)來(lái)獲得一個(gè)低復(fù)雜度的流場(chǎng)數(shù)據(jù).值得指出的是,對(duì)于現(xiàn)在所考慮的均勻湍流系統(tǒng),可以證明其傅里葉模態(tài)就是POD模態(tài)[1].

根據(jù)式(1)～ 式(5),選取Nc=84 480 個(gè)速度模的二維切片(訓(xùn)練集),可以得到kf=8和kf=32 時(shí)POD 模態(tài)的特征值曲線如圖1 所示.圖1 表明對(duì)于kf=8和kf=32 的數(shù)據(jù)分別有193 和3205 個(gè)特征值非零的POD 模態(tài),這與相應(yīng)的傅里葉模態(tài)數(shù)目相等.

圖1 POD 模態(tài)的特征值曲線Fig.1 Eigenvalues of the POD modes

獲得系統(tǒng)的POD 模態(tài)后,選取20 480 個(gè)非訓(xùn)練集的數(shù)據(jù)作為測(cè)試集,來(lái)研究不同破損區(qū)域?qū)appy POD 重構(gòu)效果的影響.圖2(a)～ 圖2(e) 給出了擁有相同破損面積Agap=32×32 和不同幾何的破損區(qū)域示意圖.其中,圖2(a)～ 圖2(d) 中的破損區(qū)域?yàn)榫鶆蚺挪嫉姆叫螀^(qū)域,但是不同圖片中的方形邊長(zhǎng)和數(shù)量不同;圖2(e) 中的破損區(qū)域?yàn)殡S機(jī)分布的像素點(diǎn).真實(shí)的完整流場(chǎng)如圖2(f)所示.

圖2 (a)-(e) 破損面積相同但破損區(qū)域幾何不同的部分缺失流場(chǎng),(f) 真實(shí)的完整流場(chǎng)Fig.2 (a)-(e) Damaged flow fields with gaps of the same area but different geometries.(f) The complete flow field

圖2 (a)-(e) 破損面積相同但破損區(qū)域幾何不同的部分缺失流場(chǎng),(f) 真實(shí)的完整流場(chǎng) (續(xù))Fig.2 (a)-(e) Damaged flow fields with gaps of the same area but different geometries.(f) The complete flow field (continued)

2.2 湍流場(chǎng)復(fù)雜程度的影響

本節(jié)研究當(dāng)分辨率相同時(shí),湍流場(chǎng)數(shù)據(jù)的復(fù)雜程度對(duì)gappy POD 重構(gòu)的影響.下面以Nflow代表湍流場(chǎng)數(shù)據(jù)的復(fù)雜程度,即特征值非零的POD 模態(tài)數(shù)目;同時(shí)以Ncomp代表gappy POD 重構(gòu)采用的POD模態(tài)數(shù)目.實(shí)際計(jì)算中往往會(huì)選取Ncomp≤Nflow,以在獲得一定重構(gòu)精度的同時(shí)降低計(jì)算量.下面對(duì)于復(fù)雜程度不同的數(shù)據(jù),考慮采用不同的Ncomp對(duì)gappy POD 重構(gòu)誤差的影響.

對(duì)于2.1 節(jié)中kf=8和kf=32 的數(shù)據(jù),分別有Nflow=193和Nflow=3205 .為了定量描述重構(gòu)效果,選取Nc=20 480 個(gè)測(cè)試數(shù)據(jù),它們與訓(xùn)練數(shù)據(jù)來(lái)自數(shù)值模擬中時(shí)間相關(guān)性很弱的不同時(shí)間段.定義破損區(qū)域中的均方重構(gòu)誤差為

其中積分域?yàn)闊o(wú)破損的完整區(qū)域,A代表完整區(qū)域面積,u0為平均速度

對(duì)于一個(gè)位于流場(chǎng)中心的 8 ×8 的正方形破損(圖2(a)),圖3對(duì)于kf=8 和kf=32的數(shù)據(jù)給出了破損區(qū)域中的均方重構(gòu)誤差MSEgap關(guān)于NcompNflow的變化曲線.注意圖中Ncomp的取值范圍為 [ 1,Nflow-1] .

圖3 破損區(qū)域中的均方重構(gòu)誤差 M SEgap 關(guān)于Ncomp/Nflow的變化曲線Fig.3 Normalized mean square error in thegap,M SEgap,as a function of Ncomp/Nflow

圖4 的最小奇異值 σ min 關(guān)于 N comp / N flow 的變化曲線Fig.4 The minimum singular value of ,σ min,as a function of N comp /Nflow

本節(jié)討論了湍流場(chǎng)的復(fù)雜程度對(duì)gappy POD 重構(gòu)的影響.湍流場(chǎng)的復(fù)雜程度由對(duì)應(yīng)的POD 模態(tài)決定,它們?cè)谝阎c(diǎn)上的值構(gòu)成了矩陣(式(13)).當(dāng)破損區(qū)域使得為列滿(mǎn)秩時(shí),重構(gòu)誤差大小約等于流場(chǎng)POD 展開(kāi)的截?cái)嗾`差除以的最小奇異值σmin(式(27)).當(dāng)流場(chǎng)復(fù)雜度較低時(shí),σmin的值并不十分小,重構(gòu)誤差主要由截?cái)嗾`差所影響,可以根據(jù)對(duì)重構(gòu)精度和計(jì)算量的要求采用較少的POD 模態(tài)來(lái)進(jìn)行重構(gòu).當(dāng)流場(chǎng)復(fù)雜度較高時(shí),由于截?cái)嗾`差會(huì)被極小的 σmin放大,此時(shí)必須采用全部POD 模態(tài)來(lái)進(jìn)行有效的重構(gòu).

2.3 破損區(qū)域面積大小和幾何形狀的影響

本節(jié)研究破損區(qū)域大小和幾何形狀對(duì)gappy POD 重構(gòu)的影響.考慮kf=32 的高復(fù)雜度的數(shù)據(jù),根據(jù)2.2 節(jié)中的討論,始終選取Ncomp=Nflow來(lái)進(jìn)行重構(gòu).

給定破損區(qū)域面積Agap,定義破損尺寸(gap size)為圖5 給出了破損區(qū)域中的均方重構(gòu)誤差MS Egap關(guān)于破損尺寸的變化曲線,其中的誤差棒代表重構(gòu)誤差在Nc=20 480 個(gè)測(cè)試數(shù)據(jù)上的標(biāo)準(zhǔn)差.可以看到,對(duì)任何一個(gè)破損形狀,總是存在一個(gè)臨界的破損尺寸,使得破損小于等于該尺寸時(shí)重構(gòu)誤差很小,而大于該尺寸時(shí)重構(gòu)誤差很大..對(duì)于圖2(a)～ 圖2(e) 中所示的破損形狀,

圖5 對(duì)于不同幾何形狀的破損,破損區(qū)域中的均方重構(gòu)誤差MS Egap關(guān)于破損尺寸的變化曲線Fig.5 Normalized mean square error in the gap,M S Egap,as a function of the gap size for different gap geometries

由圖5 和圖6 還可以看出,對(duì)于不同的幾何形狀,相應(yīng)的臨界破損尺寸也不同.例如,對(duì)于圖2(a)所示的一個(gè)方形的破損,臨界破損尺寸為12;而對(duì)于圖2(e) 所示的隨機(jī)點(diǎn)破損,臨界破損尺寸為30.這表明對(duì)于同樣的破損尺寸,當(dāng)破損形狀不同時(shí),gappy POD 的重構(gòu)效果不同.圖7 給出了破損尺寸為16 時(shí),不同破損幾何的重構(gòu)結(jié)果圖.其中,前3 列分別為破損區(qū)域、重構(gòu)結(jié)果和原始區(qū)域的示意圖;第4 列為紅色參考線上重構(gòu)結(jié)果與真實(shí)結(jié)果的分布;第5 列顯示了在固定的x1處的平均重構(gòu)誤差

圖6 對(duì)于不同幾何形狀的破損,N flow-r 關(guān)于破損尺寸的變化曲線Fig.6 N flow-r as a function of the gap size for different gap geometries

圖7 相同破損面積,不同破損幾何的gappy POD 重構(gòu)結(jié)果(對(duì)應(yīng)不同行).第1 列:破損的流場(chǎng);第2 列:重構(gòu)的流場(chǎng);第3 列:原始流場(chǎng);第4 列:第一列中紅色參考線上重構(gòu)結(jié)果(虛線)與真實(shí)結(jié)果(實(shí)線)的分布;第5 列:平均重構(gòu)誤差(紅色),Δ u(x1),及其在測(cè)試集上的平均(黑色),〈 Δu(x1)〉 . x1 和 x2 分別代表水平和豎直方向Fig.7 Gappy POD reconstruction results for the same gap area and different gap geometries (one for each row).1st column:damaged image in input.2nd column:image generated in output.3rd column:ground truth.4th column:generated (dashed) and ground truth (solid) profiles along the vertical line shown in the 1st column.5th column:mean reconstruction error,Δ u(x1) ,for each image (red line) and the average error,〈 Δu(x1)〉,over the test images (black curve).Note that x1 and x2 denotes the horizontal and vertical directions,respectively

式中的積分區(qū)域?yàn)榇嬖谄茡p的x2位置,Lgap為相應(yīng)的破損長(zhǎng)度. 〈 Δu(x1)〉 代表 Δu(x1) 在測(cè)試數(shù)據(jù)上的平均.可以看到,除了一個(gè)大的正方形破損,gappy POD 對(duì)于其他幾何形狀的破損均可精確重構(gòu)流場(chǎng).雖然破損區(qū)域的總面積相同,相比于其他幾何形狀,一個(gè)大的正方形破損會(huì)導(dǎo)致?lián)p失的信息具有更大的相關(guān)性,因此導(dǎo)致gappy POD 無(wú)法重構(gòu)流場(chǎng).

本節(jié)通過(guò)考慮kf=32 的高復(fù)雜度的湍流數(shù)據(jù),并選取Ncomp=Nflow,研究了破損區(qū)域大小和幾何形狀對(duì)gappy POD 重構(gòu)的影響.對(duì)于任何一種幾何形狀的破損,都對(duì)應(yīng)著一個(gè)保證重構(gòu)成功的最大尺寸,稱(chēng)為臨界破損尺寸.破損的幾何形狀影響著損失的信息具有的相關(guān)性,從而影響著上述臨界破損尺寸的大小.

3 結(jié)論

本文研究了gappy POD 在湍流數(shù)據(jù)重構(gòu)上的應(yīng)用,主要考慮了湍流場(chǎng)復(fù)雜程度以及破損區(qū)域的大小和幾何形狀的影響.

通過(guò)更嚴(yán)格地表述了gappy POD 重構(gòu)的過(guò)程,本文指出重構(gòu)誤差由兩部分構(gòu)成.第一部分來(lái)自流場(chǎng)POD 展開(kāi)的截?cái)嗾`差,經(jīng)由POD 基函數(shù)在已知點(diǎn)上組成的矩陣的最小特征值放大.當(dāng)湍流場(chǎng)的復(fù)雜度較低時(shí),這一項(xiàng)隨著采用的POD 模態(tài)數(shù)目增大而減小,因此實(shí)際應(yīng)用中可根據(jù)精度和計(jì)算量要求合適選取POD 模態(tài)數(shù)目.然而,當(dāng)湍流場(chǎng)復(fù)雜度較高時(shí),POD 基函數(shù)在已知點(diǎn)上組成的矩陣的最小特征值非常小,此時(shí)即使有很小的POD 截?cái)嗾`差也會(huì)被極大地放大.因此,為了對(duì)高復(fù)雜度湍流場(chǎng)進(jìn)行g(shù)appy POD 重構(gòu),必須采用其所有POD 模態(tài)使得截?cái)嗾`差為零.

Gappy POD 重構(gòu)誤差的第二部分來(lái)自POD 基函數(shù)在已知點(diǎn)上組成的矩陣的非列滿(mǎn)秩性,它主要取決于破損區(qū)域的面積大小和幾何形狀.對(duì)于固定幾何形狀的破損,當(dāng)破損區(qū)域增大時(shí),流場(chǎng)損失的信息逐漸增多,使得上述矩陣逐漸變?yōu)榉橇袧M(mǎn)秩的,從而使得重構(gòu)失效.而對(duì)于相同的破損面積,當(dāng)破損導(dǎo)致的損失信息包含的相關(guān)性越大(如一個(gè)大的方形破損),上述矩陣越容易變?yōu)榉橇袧M(mǎn)秩的,從而使得重構(gòu)失效;而當(dāng)破損導(dǎo)致的損失信息包含的相關(guān)性不大時(shí)(如隨機(jī)點(diǎn)組成的破損),上述矩陣易于維持列滿(mǎn)秩性質(zhì),此時(shí)流場(chǎng)可以被精確重構(gòu).