王曉宇,賀靜,柯漢兵,林梅

(1.西安交通大學(xué)能源與動力工程學(xué)院,710049,西安;2.武漢第二船舶設(shè)計研究所,430205,武漢)

粒子圖像測速(PIV)技術(shù)是一種誕生于20世紀(jì)80年代的無擾動流場測試技術(shù),它克服了傳統(tǒng)接觸式流場測量中單點測量的局限性,可以對存在變化流速及各類渦結(jié)構(gòu)的復(fù)雜流動實現(xiàn)全流場瞬時非接觸測量[1-2]。在PIV流場圖像后處理中常采用光學(xué)楊氏條紋法、自相關(guān)法或互相關(guān)法,逐點處理PIV圖像,從而獲得流場速度數(shù)據(jù)[3-4]。由于在實驗中可能存在示蹤粒子分布不均勻、PIV圖像變形或扭曲、PIV圖像背景噪聲、成像系統(tǒng)變形等因素,獲得的流速數(shù)據(jù)中會存在干擾正常流場分析的錯誤數(shù)據(jù)。當(dāng)錯誤數(shù)據(jù)造成的誤差超過2%時,便會影響正常的流場分析研判[5]。

針對PIV流場數(shù)據(jù)中錯誤數(shù)據(jù)的消除方法,國內(nèi)諸多學(xué)者都進(jìn)行了相關(guān)研究。王平讓總結(jié)了錯誤數(shù)據(jù)產(chǎn)生的原因,并從PIV圖像濾波處理的角度提出了消除錯誤數(shù)據(jù)的方法[6]。王燦星等建立了基于互相關(guān)技術(shù)的連續(xù)PIV圖像之間的粒子對應(yīng)方法,提出了識別錯誤數(shù)據(jù)的方法和消除準(zhǔn)則[7]。高琪等基于速度場的本征正交分解后處理技術(shù),通過迭代方法有效地實現(xiàn)了錯誤數(shù)據(jù)的識別與修復(fù)[8]。吳龍華等在對人腦判別PIV錯誤數(shù)據(jù)方式進(jìn)行模擬的基礎(chǔ)上,建立了可以進(jìn)行準(zhǔn)確錯誤數(shù)據(jù)識別的多證據(jù)推理Hopfield神經(jīng)網(wǎng)絡(luò)模型[9]。王晉軍等引入不可壓縮連續(xù)性方程對PIV測得的流場數(shù)據(jù)進(jìn)行后處理,使得修正后的流場嚴(yán)格滿足差分形式的連續(xù)性方程,實現(xiàn)了流場數(shù)據(jù)的準(zhǔn)確修正[10]。阮曉東等采用Delaunay網(wǎng)格技術(shù)建立了基于流體連續(xù)性假設(shè)的自動檢測PIV數(shù)據(jù)中錯誤數(shù)據(jù)的算法,實現(xiàn)了錯誤數(shù)據(jù)的識別與替換[11]。張工利用Visual C++開發(fā)了通用的PIV圖像后處理軟件,實現(xiàn)了PIV圖像的增強(qiáng)等功能,可以從PIV圖像中消除部分錯誤數(shù)據(jù)[12]。王元等廣泛研究了各類提高流場測量及數(shù)據(jù)后處理精度的粒子追蹤算法,減少了錯誤數(shù)據(jù)的產(chǎn)生[13-16]。上述方法對PIV流場數(shù)據(jù)中錯誤數(shù)據(jù)的剔除和修正都有較好的效果,但也存在處理流程復(fù)雜、效率不高、引入非物理偏差等不足。

本研究提出了一種基于經(jīng)驗?zāi)B(tài)分解(EMD)的PIV流場數(shù)據(jù)后處理方法,將數(shù)據(jù)分解得到的波形異常的本征模態(tài)分量(IMF)進(jìn)行濾波處理并與其余波形平滑的正常本征模態(tài)分量進(jìn)行反向疊加,重構(gòu)流場數(shù)據(jù),使錯誤數(shù)據(jù)得到有效剔除。當(dāng)分解得到的本征模態(tài)分量較多且錯誤數(shù)據(jù)集中在某一個本征模態(tài)分量時,直接將與原始流場數(shù)據(jù)相關(guān)的本征模態(tài)分量反向疊加重構(gòu)流場數(shù)據(jù),摒棄不相關(guān)的本征模態(tài)分量,實現(xiàn)錯誤數(shù)據(jù)的直接剔除。本方法已經(jīng)在多個模擬或?qū)嶒灚@取的變化較為緩和的平穩(wěn)流場數(shù)據(jù)處理中得到了良好的使用,本文利用該方法對2個典型的人為添加錯誤矢量的模擬標(biāo)準(zhǔn)流場進(jìn)行數(shù)據(jù)處理,同時對某實驗的原始流場數(shù)據(jù)進(jìn)行修正,以上述3個流場為例,說明并驗證該方法的可行性。

1 基于經(jīng)驗?zāi)B(tài)分解的流場數(shù)據(jù)重構(gòu)原理

1.1 經(jīng)驗?zāi)B(tài)分解概述

經(jīng)驗?zāi)B(tài)分解是Huang等于1998年提出的一種無需基函數(shù)的自適應(yīng)性信號數(shù)據(jù)處理新方法,在各工程領(lǐng)域得到了廣泛、高效的應(yīng)用[17-18]。該方法利用數(shù)據(jù)自身的尺度特征進(jìn)行數(shù)據(jù)分解,特別適用于非平穩(wěn)、非線性數(shù)據(jù)的平穩(wěn)化處理。經(jīng)驗?zāi)B(tài)分解處理信號的基本步驟如下。

(1)確定隨機(jī)信號x(t)在整個時間尺度上的所有極大值點和極小值點,并用三次樣條差值曲線分別連結(jié)所有極大值點和極小值點,形成上包絡(luò)線a(t)和下包絡(luò)線b(t),求取上、下包絡(luò)線的平均值

(1)

(2)求取原始隨機(jī)信號x(t)與m(t)之間的差值

h(t)=x(t)-m(t)

(2)

(3)如果h(t)的極值點與過零點的個數(shù)最多相差一個且h(t)的上、下包絡(luò)線的均值為0,則h(t)為第一個本征模態(tài)分量c1。求取x(t)與c1的差值

r(t)=x(t)-c1

(3)

將r(t)視為最新的隨機(jī)信號,重新開始上述分解過程(1)～(3)求取剩余的本征模態(tài)分量,直至r(t)變?yōu)閱握{(diào)信號。如果h(t)不能滿足上述條件,則將h(t)視為最新的隨機(jī)信號,重新開始上述分解過程(1)～(3)。

關(guān)于經(jīng)驗?zāi)B(tài)分解更加具體的數(shù)學(xué)化描述及證明過程請參見文獻(xiàn)[19-20]。經(jīng)過經(jīng)驗?zāi)B(tài)分解后,一組數(shù)據(jù)可以被分解為若干個表征不同數(shù)據(jù)尺度和數(shù)據(jù)特征的本征模態(tài)分量,它們是組成原始信號的不同尺度的一系列分量,是對原始信號特征的一系列表達(dá)。理論上連續(xù)的兩個本征模態(tài)分量的標(biāo)準(zhǔn)差在0.2～0.3之間時,則分解得到的本征模態(tài)分量既能滿足線性條件,穩(wěn)定性良好,又能夠較好地反映原始隨機(jī)信號的性質(zhì)。經(jīng)驗?zāi)B(tài)分解具有正交性和完備性的特征,因此將分解得到的本征模態(tài)分量反向疊加即可得到原始數(shù)據(jù)。在處理現(xiàn)代流動測試技術(shù)所獲得的流場數(shù)據(jù)時,數(shù)據(jù)量的大小并不會直接影響所分解出的本征模態(tài)分量數(shù)量,而原始流場的復(fù)雜程度則對本征模態(tài)分量的數(shù)量影響較大。原始數(shù)據(jù)經(jīng)過經(jīng)驗?zāi)B(tài)分解后一般可以得到5～15個本征模態(tài)分量,其中本征模態(tài)分量在5～10個時可以認(rèn)為較少,超過10個時可以認(rèn)為分解得到的本征模態(tài)分量較多。

1.2 相關(guān)系數(shù)概述

在數(shù)據(jù)處理中常用相關(guān)系數(shù)來表征兩組數(shù)據(jù)之間的相關(guān)程度,在本研究中通過求解不同的本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù),可以判斷該分量與原始流場數(shù)據(jù)的相關(guān)程度。相關(guān)系數(shù)的劃分并沒有固定的標(biāo)準(zhǔn),而是根據(jù)所處理數(shù)據(jù)的具體特征來決定的。表1所示的是一種相關(guān)性的劃分方式,本研究根據(jù)這種相關(guān)系數(shù)大小的劃分方式將本征模態(tài)分量與原始流場數(shù)據(jù)的關(guān)系歸類為不相關(guān)、微相關(guān)、實相關(guān)、高度相關(guān)4類,其中微相關(guān)、實相關(guān)、高度相關(guān)統(tǒng)稱為相關(guān)。

表1 相關(guān)性劃分

1.3 流場數(shù)據(jù)處理流程

圖1是本研究所提出的基于經(jīng)驗?zāi)B(tài)分解的PIV流場數(shù)據(jù)后處理方法流程圖,其中實線表示的流程是通用方法,虛線表示的流程是在原始流場較為復(fù)雜、原始流場數(shù)據(jù)經(jīng)過經(jīng)驗?zāi)B(tài)分解后得到的本征模態(tài)分量較多且錯誤數(shù)據(jù)集中在某一個本征模態(tài)分量時對通用方法的一種簡化。

圖1 基于EMD的PIV流場數(shù)據(jù)后處理流程圖

錯誤數(shù)據(jù)在數(shù)值大小和方向上與正確數(shù)據(jù)有明顯的差異,在原始流場數(shù)據(jù)沒有經(jīng)過經(jīng)驗?zāi)B(tài)分解時,錯誤數(shù)據(jù)摻雜在正確數(shù)據(jù)中難以識別。當(dāng)用經(jīng)驗?zāi)B(tài)分解將原始流場數(shù)據(jù)分解為若干個本征模態(tài)分量時,錯誤數(shù)據(jù)就會造成本征模態(tài)分量曲線出現(xiàn)明顯的局部不平滑。由于錯誤數(shù)據(jù)造成本征模態(tài)分量曲線局部的不平滑十分明顯,因此常用的濾波方法就可以將錯誤數(shù)據(jù)點在本征模態(tài)分量中濾除而不誤刪正確數(shù)據(jù)。根據(jù)經(jīng)驗?zāi)B(tài)分解的正交性和完備性,波形異常的本征模態(tài)分量經(jīng)濾波處理后與波形平滑正常的本征模態(tài)分量反向疊加就可以重構(gòu)剔除了錯誤數(shù)據(jù)的流場數(shù)據(jù)。上述的原始流場數(shù)據(jù)處理思路可表示為圖1中實線所示的流程。

原始流場數(shù)據(jù)經(jīng)過經(jīng)驗?zāi)B(tài)分解后得到的本征模態(tài)分量越多,每個本征模態(tài)分量攜帶的原始流場數(shù)據(jù)就越少,所占原始流場數(shù)據(jù)的比例就越小。在有些實驗中,粒子圖像測速技術(shù)所獲得的流場數(shù)據(jù)相對比較復(fù)雜且錯誤數(shù)據(jù)比較明顯,這時分解得到的本征模態(tài)分量較多且錯誤數(shù)據(jù)會集中在某一個本征模態(tài)分量。摻雜了錯誤數(shù)據(jù)的本征模態(tài)分量波形異常且與原始流場數(shù)據(jù)的相關(guān)性較低(呈不相關(guān)),因此在該流程中可以通過求解各本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)性大小來判斷錯誤數(shù)據(jù)所在的本征模態(tài)分量。為了提高處理速度和效率,可以按照圖1中虛線所示的處理流程,省略相關(guān)性低(呈不相關(guān))的本征模態(tài)分量的濾波過程,直接將該分量刪除并將其余相關(guān)性較高的本征模態(tài)分量反向疊加直接重構(gòu)流場數(shù)據(jù)。隨著錯誤數(shù)據(jù)的刪除,也會損失該本征模態(tài)分量中的正確數(shù)據(jù),但在分解得到的本征模態(tài)分量較多的情況下,損失的小部分正確數(shù)據(jù)不會對整體流場數(shù)據(jù)造成明顯的影響,流場的特性仍然可以清楚地表達(dá)。

2 算例驗證與應(yīng)用分析

2.1 模擬流場中的驗證與分析

圖2是用計算機(jī)模擬的流速為1 m/s的單向均勻流場,在4 000個原始流場數(shù)據(jù)中加入了明顯的錯誤數(shù)據(jù),人為添加的誤差約為1.7%,因此速度矢量場中出現(xiàn)了數(shù)值大小和方向明顯錯誤的速度矢量。圖3是原始流場數(shù)據(jù)按圖1中實線流程處理后得到的速度矢量場。圖4和圖5是對該原始流場x、y方向全部的4 000個流場數(shù)據(jù)進(jìn)行經(jīng)驗?zāi)B(tài)分解后得到的本征模態(tài)分量,分別為6個和8個,不符合使用簡化流程的要求,因此采用的是通用方法。由圖可知x、y方向的錯誤數(shù)據(jù)都集中在了第1個本征模態(tài)分量中,因此該本征模態(tài)分量的曲線中出現(xiàn)了不平滑的跳變點,也正是錯誤數(shù)據(jù)存在的位置。對x、y方向上錯誤數(shù)據(jù)集中的第1個本征模態(tài)分量進(jìn)行濾波處理,即可把錯誤數(shù)據(jù)存在的跳變點濾除而不誤刪正確數(shù)據(jù)。將x、y方向上濾波處理后的第1個本征模態(tài)分量與其他本征模態(tài)分量反向疊加,即可重構(gòu)流場數(shù)據(jù)。圖2與圖3對比表明,實線表示的處理流程能有效實現(xiàn)該單向均勻流場錯誤數(shù)據(jù)的剔除和速度矢量場的修正,修正并重構(gòu)后的流場數(shù)據(jù)與正確流場數(shù)據(jù)的誤差約為0.002%。

圖2 原始單向流場數(shù)據(jù)得到的速度矢量場

圖3 處理后的單向流場數(shù)據(jù)得到的速度矢量場

圖4 原始單向流場x方向數(shù)據(jù)的本征模態(tài)分量

圖5 原始單向流場y方向數(shù)據(jù)的本征模態(tài)分量

為了驗證所提出的方法在處理相對復(fù)雜的流場數(shù)據(jù)時的有效性,本研究模擬了速度大小不定的單向旋轉(zhuǎn)流場,并在原始流場數(shù)據(jù)中加入了一些明顯的錯誤數(shù)據(jù),人為添加的誤差約為3.3%。圖6和圖7是對該原始流場x、y方向的流場數(shù)據(jù)進(jìn)行經(jīng)驗?zāi)B(tài)分解后得到的本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù),本征模態(tài)分量的數(shù)量和特點符合上文提出的虛線表示的簡化處理流程的要求,因此處理該流場數(shù)據(jù)時采用的是簡化流程。圖8和圖9分別是由原始流場數(shù)據(jù)和按圖1中虛線表示的簡化流程處理后的流場數(shù)據(jù)得到的速度矢量場。對比表明,圖1中虛線表示的簡化處理流程能夠有效地實現(xiàn)該速度大小不定的單向旋轉(zhuǎn)流場中錯誤數(shù)據(jù)的剔除和速度矢量場的修正,修正并重構(gòu)后的流場數(shù)據(jù)與正確流場數(shù)據(jù)的誤差約為0.18%。

圖6 旋轉(zhuǎn)流場x方向流場數(shù)據(jù)的各本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù)

圖7 旋轉(zhuǎn)流場y方向流場數(shù)據(jù)的各本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù)

圖8 原始旋轉(zhuǎn)流場數(shù)據(jù)得到的速度矢量場

圖9 處理后的旋轉(zhuǎn)流場數(shù)據(jù)得到的速度矢量場

2.2 實驗流場中的應(yīng)用與分析

為了驗證本文所提出的基于經(jīng)驗?zāi)B(tài)分解的PIV流場數(shù)據(jù)后處理方法在實驗流場數(shù)據(jù)處理中的準(zhǔn)確性與實用性,本研究對低速風(fēng)洞中某一T型通道的PIV流場原始數(shù)據(jù)進(jìn)行了處理與分析。關(guān)于該T型通道PIV流動測試實驗的具體細(xì)節(jié)請參見文獻(xiàn)[21-22]。

正如上文所述,實驗流場的復(fù)雜性一般相對較高,流場數(shù)據(jù)經(jīng)過經(jīng)驗?zāi)B(tài)分解之后得到的本征模態(tài)分量較多。該實驗流場x、y方向的原始流場數(shù)據(jù)經(jīng)過經(jīng)驗?zāi)B(tài)分解之后分別得到了10個本征模態(tài)分量和12個本征模態(tài)分量,可以用圖1中的簡化流程進(jìn)行數(shù)據(jù)處理。圖10和圖11分別表示了x方向流場數(shù)據(jù)、y方向流場數(shù)據(jù)的本征模態(tài)分量與原始流場數(shù)據(jù)之間的相關(guān)程度。由圖可知,x方向流場數(shù)據(jù)的第1個本征模態(tài)分量和y方向流場數(shù)據(jù)的第6個本征模態(tài)分量與原始流場數(shù)據(jù)的整體相關(guān)性很小(呈不相關(guān)),是錯誤數(shù)據(jù)集中的具體表現(xiàn)。該流場數(shù)據(jù)可以按照圖1中虛線表示的流程進(jìn)行處理,直接刪除這兩個本征模態(tài)分量,將x、y方向上剩余的9個、11個本征模態(tài)分量反向疊加重構(gòu)流場數(shù)據(jù)。圖12、圖13分別表示了該流場數(shù)據(jù)處理前后的速度矢量場,可知原始流場數(shù)據(jù)中位于速度矢量場上部的錯誤數(shù)據(jù)(圖12中區(qū)域1)得到了準(zhǔn)確的剔除,一些明顯的錯誤矢量得到了有效的消除,速度矢量場更加平滑,品質(zhì)得到了顯著改善,速度矢量場上部的流動特征(圖13中區(qū)域1＇)也更加趨近于文獻(xiàn)[22]報道的各工況下的真實流動情況,但在個別不存在錯誤矢量的區(qū)域出現(xiàn)了流場特征的輕微變化。這是隨著簡化流程中呈不相關(guān)的本征模態(tài)分量的直接刪除,造成少部分正確數(shù)據(jù)損失所導(dǎo)致的。如果對局部流場特性的準(zhǔn)確性要求較高,可以不采用簡化流程,而采用相對較為煩瑣的通用方法來進(jìn)行處理。

圖10 實驗流場x方向流場數(shù)據(jù)的各本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù)

圖11 實驗流場y方向流場數(shù)據(jù)的各本征模態(tài)分量與原始流場數(shù)據(jù)的相關(guān)系數(shù)

圖12 原始實驗流場數(shù)據(jù)得到的速度矢量場

圖13 處理后的實驗流場數(shù)據(jù)得到的速度矢量場

3 結(jié) 論

本研究提出了一種基于經(jīng)驗?zāi)B(tài)分解的PIV流場數(shù)據(jù)后處理方法,對存在錯誤數(shù)據(jù)的本征模態(tài)分量進(jìn)行濾波處理或直接刪除,將其余本征模態(tài)分量反向疊加,實現(xiàn)了錯誤數(shù)據(jù)的剔除與流場數(shù)據(jù)的重構(gòu)。利用該方法處理了計算機(jī)模擬流場和某實驗流場的原始數(shù)據(jù),結(jié)果表明錯誤數(shù)據(jù)得到了有效剔除,重構(gòu)的流場數(shù)據(jù)可更加真實地反映流場特征。

目前任何一種PIV流場數(shù)據(jù)后處理方法都不能完全將錯誤數(shù)據(jù)剔除而不影響正確數(shù)據(jù),本研究所提出的方法也是如此,但本方法具備以下兩個優(yōu)點:首先,該方法最大限度地從原始流場數(shù)據(jù)中剔除了錯誤數(shù)據(jù),并非通過各類算法在計算時“避開”錯誤數(shù)據(jù),具備處理的徹底性;其次,該方法簡單易行,沒有復(fù)雜的錯誤數(shù)據(jù)識別流程與標(biāo)準(zhǔn),可以快捷簡便地提高流場分析的準(zhǔn)確性。在PIV技術(shù)中,錯誤數(shù)據(jù)超過2%就已經(jīng)會影響流場分析,證明實驗中存在的誤差太大,此時解決誤差大的最好辦法是重新進(jìn)行實驗,而不是消除已有的錯誤數(shù)據(jù)。就這一點而言,本方法對誤差在2%以內(nèi)的流場數(shù)據(jù)是有效的,但本方法也有自身的缺點,當(dāng)所處理的流場變化非常劇烈、雜亂(例如射流撞擊壁面)時,本方法容易把部分真實反映流場的數(shù)據(jù)視為錯誤數(shù)據(jù),因此本方法只適用于處理變化較為緩和的平穩(wěn)流場。

受該方法的啟發(fā),也可以將原始的PIV圖像視為一個二維數(shù)組,利用經(jīng)驗?zāi)B(tài)分解的方法通過上述流程來對其進(jìn)行處理。這時處理的對象是PIV原始粒子圖像,各個本征模態(tài)分量也是以圖像或數(shù)組的形式來呈現(xiàn)的,最終通過本征模態(tài)分量的反向疊加得到處理后的PIV圖像,并從中提取流場數(shù)據(jù)。此外,本方法在三維流場數(shù)據(jù)中的應(yīng)用及可行性還有待驗證,在接下來的研究工作中,上述兩部分內(nèi)容將是研究的重點。