基于本征正交分解(POD)的PIV數(shù)據(jù)壞點剔除方法

1、基于本征正交分解（POD）的PIV數(shù)據(jù)壞點辨識方法摘要：1.引言：粒子圖像測速技術是一種基于圖像互相關的激光測速技術。在進行PIV實驗的過程的由于示蹤粒子濃度不均勻、激光強度分布不均、粒子成像質量差等原因容易造成互相關峰值不確定1，出現(xiàn)所謂的“壞點”。最近幾年發(fā)展起來的3D3CPIV2, 3，由于在重構空間粒子場的時候會不可避免的引入虛假粒子，進一步增加了流場中了壞矢量。因此，壞矢量的剔除和重建是PIV數(shù)據(jù)后處理的一個重要內(nèi)容。一般通過對比該點誤差與周圍3*3或5*5鄰域內(nèi)平均誤差來確定該點是不是壞點。目前比較通用并且效果很好的是Westerweel和Scarano 提出的歸一化中值檢測方法4

2、（normalized median test）?？紤]一個位移矢量和其周圍的相鄰矢量。用表示的中值，用表示殘差，定義為（i=1,2,8）。的殘差用的中值歸一化后得到下面的公式：其中是與流場的平均噪音有關的量，一般設為0.1-0.2 pixel。對于平面PIV，的閾值設為24，大于該閾值的點被認為是壞點。其他壞點識別的方法與該方法原理相同，基本上都是基于相鄰點的統(tǒng)計特性判斷該點的性質。這種方法有缺陷，一是只能用到少量相鄰點的值，忽略了流場的全場特性；二是局部相鄰點的速度矢量很有可能也是壞點，導致判斷的結果不可靠；三是的閾值需要根據(jù)流中壞點的多少確定，在實際運用中并不能精確的知道壞點所占的比例。因

3、此，本文作者提出了基于POD的壞點檢測方法。該方法對周期性或類周期性流場PIV數(shù)據(jù)的后處理有很好的效果。本文的結構安排如下：首先在第二節(jié)簡單的介紹了本征正交分解的原理，推導了POD與流場湍動能的關系；在第三節(jié)對流場的誤差進行了數(shù)值模擬，討論了誤差對POD分解的影響以及詳細的介紹了基于POD的壞點剔除方法；一個真實PIV計算出的流場運用該方法進行了壞點剔除并和歸一化中值檢測方法進行了比較，這在第四節(jié)給出；第五節(jié)是結論。2.POD原理POD是一種從統(tǒng)計意義上提取流場中主要流動結構的方法，可以實現(xiàn)對復雜非線性系統(tǒng)的線性降維處理。這種方法最早由Lumley（1967年）引入湍流的研究，用于辨識大尺度擬

4、序結構，后來Berkooz5等人對POD方法給出了系統(tǒng)的介紹。POD方法就是要找出和原流場最相似的流動模態(tài)5，對于復雜流場這樣的流動模態(tài)不止一個，于是可以假設流場是一系列POD基模態(tài)的線性組合，即：其中是不同時刻的速度場，代表POD的第k階模態(tài)，代表與第k階模態(tài)相關的時間系數(shù)。POD的模態(tài)滿足正交關系：從數(shù)學上可以從下面的公式解出： （1）其中C是和的兩點互相關系數(shù)，即。這里我們著重討論一下特征值的取值。對于兩點互相關矩陣C中的元素Cij有計算公式：一般都是用速度的脈動場進行POD分解的，式中代表的是速度脈動。當擁有足夠多的不同時刻流場的數(shù)據(jù)時，不論是對于classical POD還是snap

5、shot POD，都代表流場的湍動能，而且不會有太大的變化，即在Cij的公式中分母基本都相同，這也正是可以直接用Kij替換Cij帶入公式8計算的原因。因此，當對歸一化后，特征值就是一個與該階模態(tài)的湍動能有關的量。越大，該階模態(tài)所含的湍動能越多，也就是說POD可以按能量提取流場的主要特征。詳細的介紹可以參考文獻5。3.誤差模擬分析3.1 模擬方法評價PIV計算結果好壞的一個重要指標就是壞點所占的百分比。一般來說，在PIV中存在兩類誤差。一類是以單個點形式出現(xiàn)的錯誤矢量，這類誤差大小和分布都帶有明顯的隨機性，主要是由于查詢窗內(nèi)互相關信噪比太低峰值不明確造成的；另一類誤差成片形式出現(xiàn)，通常有好幾個誤

6、差向量集中在同一個區(qū)域，這種誤差很有可能是由于圖像在這一區(qū)域粒子太少或質量太差造成的。設誤差向量的連通區(qū)域的大小為，當時表示單個誤差矢量，表示兩個誤差矢量連在一起，以此類推。為了得到的分布，對一個PIV計算出來的真實的流場用歸一化中值檢測方法檢識別壞點，然后統(tǒng)計壞點連通區(qū)域的分布。得到圖 1中的統(tǒng)計結果。從直方圖中可以看出單個誤差矢量出現(xiàn)的概率最高，連通區(qū)域越大出現(xiàn)的概率越小。對數(shù)據(jù)進行高斯擬合得到光滑的擬合曲線：代表該類型誤差占所有向量的比例。誤差的大小在查詢窗內(nèi)隨機分布。當流場誤差總的個數(shù)確定后，按該概率密度函計算不同類型的誤差并添加到基本流場之上，完成誤差場的模擬。圖 1誤差向量連通區(qū)域

7、的概率密度分布圖 2Q=3%時模擬結果基本流場取自一套平滑過后的圓柱擾流的平面PIV數(shù)據(jù)。該數(shù)據(jù)在北京航空航天大學低速水槽完成，自由來流速度為35mm/s。圓柱水平放置在水槽中間直徑為10mm。雷諾數(shù)為250。在流場中撒播直徑為5m ，密度為1.05gmm-3的空心玻璃微珠作為示蹤粒子，用一臺2w的激光器照亮測量區(qū)域。相機的分辨率為640*480，采樣頻率為80Hz。在進行互相關計算時采用了窗口變形算法，查詢窗最終大小為16*16,50%重疊區(qū)。本文模擬了8中不同的工況，分別是Q=0.5%,1%,2%,3%,4%,5%,7%,10%。對其中Q=3%的情況進行了詳細的分析。圖 2給出了Q=3%時

8、模擬的結果。3.2 誤差向量對POD分解的影響POD分解按流場的湍動能提取擬序結構，而誤差向量直接影響流場湍動能。直觀上判斷流場中壞點越多，流場的湍動能就會越大，而且這種湍動能的增加并不是流動結構的變化引起的。設Q代表誤差向量的個數(shù)占流場全部向量個數(shù)的百分比，代表引入誤差后流場湍動能的值與真實流場湍動能值的比值，當Q=0時=1。圖 3給出了隨的變化規(guī)律。從圖中可以看出，明顯隨Q呈線性變化，增加1%的誤差向量流場的湍動能增加23.27%。但是值得注意的是這種湍動能是由于隨機誤差造成的，并不包含流動結構。圖 3流場湍動能隨Q的變化規(guī)律接下來我們考慮誤差向量對POD分解的影響。在前面一節(jié)我們已經(jīng)知道

9、基于兩點互相關的POD按湍動能提取流場的流動特征。那么，由于引入了誤差導致湍動能增加會影響POD分解的模態(tài)嗎？怎樣影響？為此，我們選擇包含0%，0.5%，1%，3%，5%，10%的誤差的流場進行本征正交分解，根據(jù)擬合的公式可知流場的湍動能分別為：100%，110.75%，121.5%，164.5%，207.5%，315%。圖 4給出了分解后的能量譜。由于流場整體能量的增加，各階模態(tài)所占的相對能量百分比下降。設各階模態(tài)的絕對能量為，為原始沒有誤差向量的流場所含總湍動能的絕對值。那么，原始場各階模態(tài)的相對能量為：當引入誤差矢量后，流場的湍動能會按照線性規(guī)律增加。增加的比例用表示。用表示該流場模態(tài)的

10、相對能量，于是得到：為了和原始模態(tài)進行比較，在上面的公式兩端都乘以，進行歸一化處理。得到的結果如圖所示。從圖中可以看出，歸一化后的各階模態(tài)的相對能量曲線在低階時很好的重合在一起。雖然流場的湍動能在不斷增加，但是主要流動結構（低階模態(tài)）的能量并沒有發(fā)生變化。這說明誤差向量并沒有影響低階模態(tài)的提取，用較少的低階模態(tài)同樣能重構出和原流場非常近似的流場。但是，從圖中我們發(fā)現(xiàn)，隨著誤差向量的越來越多，對流場的影響逐漸深入到低階高能模態(tài)，而且模態(tài)的階數(shù)越來越多。這與流場的相關性降低有關。高階模態(tài)反映了流場的隨機誤差，雖然隨機誤差場占有很大一部分湍動能，但這些湍動能并不包含流動結構并且每階模態(tài)雜亂無章。低階

11、模態(tài)包含流場的大尺度流動結構，這些模態(tài)所含的能量越低，越容易受隨機誤差的影響。但是值得注意的是，隨機誤差在這些流動結構中所含的能量是十分有限的，大部分誤差的能量仍然集中在高階模態(tài)，并且各階模態(tài)之間能量差別很小。正是這樣的原因導致了POD模態(tài)的急劇增加，達到上千階。對于含誤差3%的PIV數(shù)據(jù)，如果選擇前15階和100階模態(tài)進行重構得到的流場如圖 5所示。從圖中可以看出，前15階得到的流場受隨機誤差的影響要小很多，重構階數(shù)越多流場受到的擾動越大。而且隨機誤差會傳播到流場中的每一個點，這種誤差均勻分布于全流場，通過對比很容易判斷出誤差的分布。圖 4不同Q下流場POD分解的能量譜圖 5Q=3%流場與P

12、OD重構對比通過上面的討論得知，隨機誤差會快速的增加流場的湍動能，但是在一般PIV誤差量級下，其能量在POD分解時還不足以影響流場的大尺度結構，即低階高能模態(tài)。雖然各個點的誤差會傳播到整個流場，導致POD的模態(tài)不光滑，但并不影響對流場的近似重構。3.3 基于POD的壞點剔除通過上面的分析我們知道，運用較少的低階高能模態(tài)重構可以引入較少的隨機誤差。對比重構的流場和原始流場，在存在壞點的地方速度矢量的差別很大，而不存在壞點的地方絕對誤差與周圍的絕對誤差相近。由于重構采用了固定的階數(shù)，不可避免的產(chǎn)生了階段誤差，但是從局部3*3 或 5*5 的區(qū)域來看，這種重構的階段誤差應該是非常相近的，這為我們判斷

13、壞點提供了方便。與Westerweel和Scarano4相同，考慮速度矢量和其周圍3*3區(qū)域的相鄰矢量，與之相對應的前m階重構的速度矢量用和。在平面PIV中每個速度矢量有兩個分量u,v。重構的絕對誤差可以看出為兩矢量差的長度。該誤差與重構階數(shù)、流場特性有關。一般而言在強剪切區(qū)或湍流度比較大的區(qū)域該誤差較大，這與POD分解的特性有關，這些區(qū)域需要更多的重構階數(shù)才能達到一定的精度。另外，該誤差還與該點速度的大小有關。由此可見用絕對誤差作為誤差判斷的標準并不能適應各種流動工況，一種簡單的解決方法就是用周圍點的絕對誤差對其歸一化。具體公式為：式中代表點誤差歸一化后的殘差，分母上為周圍8個點的絕點誤差的

14、最小值。和Westerweel和Scarano文獻不同，之所以選擇最小值，是加速收斂和對成片壞點有更好的剔除效果。為一常數(shù)，反映流場中整個誤差的平均水平?？梢灶A見在光滑的區(qū)域應該是趨近于1的，而在存在壞點的地方殘差應該很大。需要說明的是，之所以分母上選擇中值而非平均，是因為中值濾波可以減少壞點的干擾，算法上更穩(wěn)定。這樣可以通過迭代剔除壞點不斷的逼近真實流場。算法的實現(xiàn)過程如下：a) 準備第k次迭代b) 對前一次迭代得到的流場進行POD分解，若k等于1，則對原始PIV數(shù)據(jù)進行POD分解。這里要求原流場擁有足夠多的幀數(shù)和周期，在時間上能夠很好的解析整個流動過程。本文選擇snapshot POD方法

15、。c) 選擇前4階模態(tài)的相對能量作為收斂條件的判斷。設代表第k次迭代后前4階模態(tài)的能量。隨著壞點的不斷剔除，低階模態(tài)的能量會越來越大，最終趨于穩(wěn)定。因此可以認為當時，程序已收斂，停止迭代；否則繼續(xù)下一步。的取值范圍為12%。d) 計算重構階數(shù)并用前階模態(tài)對流場進行重構。當k=1時，=4。當k>=2時，取原始數(shù)據(jù)能量譜和第次迭代完成并插值后數(shù)據(jù)能量譜的交點對應的階數(shù)。如圖XX所示，該點位置在迭代過程基本不變。e) 按照公式計算殘差。式中，即該次迭代所有誤差的平均值。并計算每一個時刻的殘差的均值和脈動。f) 設置殘差的閾值，若殘差大于該閾值則認為是壞點被剔除。這里3是較為合理的，若太小則會剔

16、除正確的矢量，若太大則要迭代更多的次數(shù)。g) 對已經(jīng)剔除的點進行插值。插值點附近依然可能存在壞點，故插值的方法盡量對壞點不敏感。線性插值、樣條插值等一般插值方法都是基于周圍數(shù)據(jù)的統(tǒng)計信息進行插值計算的，很容易受到噪音的影響。Gappy POD是一種基于全流場信息的插值方法，在這里非常適用。但由于插值并不是本文的重點，這里并不引入該方法。本文選擇了一種基于POD重構流場的插值方法?？紤]第k次迭代并剔除壞點后的速度矢量和其周圍3*3區(qū)域的速度，與之相對應的前階重構的速度矢量用和表示。用表示誤差，即。若被當做壞點剔除，則可用下面的公式計算其插值：h) 插值后的流場作為新的原始流場，進入第a步。3.4

17、 數(shù)值評估按照前面介紹的誤差模擬方法，對基本流場添加5%的誤差向量并運用POD方法識別壞點。按照上面所描述的方法，設置收斂條件=2%，6次迭代后滿足收斂條件。由于基本流場是一個周期性很好的流動，前兩階模態(tài)的能量占流場湍動能的90%左右。因此計算出的重構階數(shù)較小，大概為16，并且隨迭代次數(shù)的變化很小，如圖XX中國黑色圓圈所示。紅色線代表真實流場的能譜，從圖XX可以看出，隨著迭代次數(shù)的增加，流場中的壞點被剔除，流場的能譜逐漸逼近真實解。同時由隨機誤差造成的高階低能模態(tài)越來越少，最終POD模態(tài)的階數(shù)也會逼近真實解。從中也可以看出，前幾階模態(tài)在真實流場中的能量應該是最大的。圖XX給出了前四階能量的收斂

18、曲線和對應壞點比例曲線。在迭代剔除壞點的過程中，壞點會逐漸被判斷出來，對應的前四階能量會逐漸變大，直到達到收斂條件。最終整個流場識別出5.2%的壞點，非常接近預設的5%。通過上面的分析，我們從整體上判斷出流場的收斂性，但是并不清楚流場壞點識別和填充的細節(jié)。從流場取出某一幀的數(shù)據(jù)，用紅色標示出識別出的壞點，得到圖XX?；旧纤袎狞c都被識別。按照POD插值后得到圖XX的結果。圖中顏色標示速度矢量的大小。通過計算插值結果與原數(shù)據(jù)的相對誤差，給出了相對誤差的概率密度分布（PDF），如圖XX。相對誤差落在0,2%范圍內(nèi)的概率為85%，落在0,3%的概率為91%。和一般的基于周圍點的數(shù)學插值方法不同，本

19、文所用的基于POD的插值方法利用了流場所有的時空信息，能夠有效的還原流動結構，對成片的壞點和強剪切區(qū)域的插值有很好的效果。為了說明這一點，在真實流場中隨機剔除速度矢量然后用兩種方法插值，比較相對誤差的大小。得到的結果如圖XX所示。橫坐標表示插值區(qū)域的大小，例如5，則代表5*5的插值區(qū)域?？v坐標是統(tǒng)計了上萬個點后相對誤差的平均值?？梢郧宄目吹剑S著區(qū)域的增大，基于POD的插值方法的相對誤差明顯低于線性插值。圖中仍然取前16階進行插值。需要說明的是，之所以進行了統(tǒng)計分析，是因為插值的精度與插值的位置有關，速度梯度變化大的地方精度會變低。同時本文所用的插值方法與重構階數(shù)的選擇有關，盡量選擇沒有被誤

20、差影響的模態(tài)，如步驟d所示。只要階數(shù)選擇合理，基于POD的插值方法在精度和穩(wěn)定性上更具有優(yōu)勢，能夠加速迭代過程的收斂。 相對誤差的概率密度分布 兩種插值方法對比從上面的結論中可以看出，在流場中壞點較少的情況下，該方法能夠很好的識別壞矢量并且能同時插值出很好的結果。為了進一步發(fā)掘該方法的優(yōu)勢以及測試該方法的有效性，需要建立更加苛刻的模擬條件。對基本流場添加15%的隨機誤差和一個5*5大小的朝一個方向的壞點區(qū)。依然按照給定的參數(shù)進行迭代。迭代7次后流場收斂，收斂過程如圖XX和XX所示。能譜依然能夠收斂到真實解。最終識別的壞點有15.2%，稍微高于設定的15%。這一次，我們與歸一化中值檢測方法做了對

21、比。圖XX是POD識別壞點的結果，圖XX是歸一化中值檢測的結果。總體上看歸一化檢測方法只識別出13%的壞點，仍有2%壞點被漏掉。注意圖中綠色圓圈內(nèi)的區(qū)域，POD識別出全部的壞點，但中值檢測方法在這里失效，這是因為基于相鄰區(qū)域統(tǒng)計的方法無法計算出正確的參考向量和誤差，導致判斷失效?？傊瑢τ趩蝹€壞點兩種方法識別的結果相同，但對于成片的壞點，基于POD的壞點剔除明顯優(yōu)于中值檢測。對于插值的結果，圖給出了相對誤差分布的云圖。圖XX是原始流場，圖XX是POD插值后的結果。在成片誤差區(qū)域相對誤差較大，單點誤差較小，與前面的分析一致。原始速度場（矢量長度） 插值后速度場 相對誤差通過上面的分析，基于POD

22、的壞點識別和其他方法4, 6相比有如下優(yōu)勢：a) 該方法基于全流場的POD分解，主要流動結構的重構能夠很好的揭示流動的物理含義，不再單純的依靠數(shù)學統(tǒng)計上的意義判斷壞點；b) 將兩個分量耦合求解，改進了歸一化中值檢測方法，更加容易識別出壞點，對成片的壞點也能很好的識別；c) POD插值方法使得在剔除壞點的同時能夠高精度的填充壞點，得到的流場既可以直接用作數(shù)據(jù)分析，也可以作為流場平滑的參考。不過該方法不能用在多層互相關的程序中，比較適合對數(shù)據(jù)進行后處理。4.實例本節(jié)將給出用該方法識別真實PIV原始數(shù)據(jù)中壞點的例子。本實驗在北京航空航天大學流體力學研究所的回流式水槽中完成，研究低速條帶在層流邊界層中

23、的演化過程。條帶由固定位置的圓柱干擾絲產(chǎn)生，在下游垂直于干擾絲放置光滑的有機玻璃板。條帶進入邊界層后，由于行波擾動失穩(wěn)，將產(chǎn)生Sinuous模式和Varicose模式兩種模態(tài)7。粒子圖像的大小為640*480，如圖XX所示。該粒子圖像的粒子濃度偏低，而且出現(xiàn)了較多的單個像素，容易導致互相關計算時峰值鎖定。從粒子圖像上就可以判斷互相關計算的結果應該含有較多的壞點。在實際計算的過程中，不做圖像的前處理，采用了圖像變形算法計算速度場，最終查詢窗口的大小為16*16,50%的重疊區(qū)。每幀圖像有60*80個速度矢量，總共有8497幀。為了便于顯示，只取其中40*40個向量。圖XX給出本實例最重要的收斂曲

24、線。由于層流邊界層中的條帶的演化是一個類周期的過程，故前四階模態(tài)所占的能量比較低，迭代10次后收斂到38%左右。而壞點的比例達到了13%。這么高的壞點比例，正如前面所講，是因為粒子圖像的質量不高同時又沒進行圖像前處理造成的。壞點識別以及插值填充的結果在圖XX中給出。針對以上問題，可以發(fā)展出一種迭代POD壞點識別程序。即在前一步識別壞點的基礎上對流場進行高精度數(shù)學插值，也可以運用 gappy POD方法對缺失的點進行高精度重構8, 9，得到新的含壞點較少的流場。在該流場的基礎上繼續(xù)前一步剔除壞點的工作，循環(huán)直到流場中無法檢測到壞點。1.Hart DP. PIV error correction. Exp Fluids. 2000; 29(1): 13-22.2.Elsinga GE, van Oudheusden BW, Scarano F. Experimental assessment of Tomographic-PIV accuracy. 2006.3.Gao Q, Wang H, Wang J. A single camera volumetric particle image velocimetry and its application. Scie