專(zhuān)利名稱(chēng)：粒子圖像測(cè)速處理方法
技術(shù)領(lǐng)域：
本發(fā)明涉及圖像處理，特別是涉及一種粒子圖像測(cè)速處理方法。
背景技術(shù)：
自然界圖像的運(yùn)動(dòng)形式千差萬(wàn)別，流體運(yùn)動(dòng)是一種典型的非剛體運(yùn)
動(dòng)，流體運(yùn)動(dòng)圖像計(jì)算與分析即粒子圖像測(cè)速(Particle Image Velocimitry,簡(jiǎn)稱(chēng)PIV)是一種新型的非接觸式的測(cè)量技術(shù)。通常PIV是 在流體中投入示蹤粒子，在激光片光的照射下，在與片光垂直方向用攝像 機(jī)拍攝隨流體運(yùn)動(dòng)的粒子圖像，再對(duì)圖像進(jìn)行處理分析計(jì)算，最終得到流 場(chǎng)橫切面上二維速度場(chǎng)的一種測(cè)量方法。作為一種全流場(chǎng)、無(wú)接觸、無(wú)擾 動(dòng)、高精度的流動(dòng)可視化方法，PIV適用于湍流、非定常流等時(shí)變復(fù)雜流 場(chǎng)的測(cè)量。如今，PIV已是一門(mén)跨學(xué)科的交叉綜合技術(shù)，其結(jié)合了激光技 術(shù)、視頻圖像處理技術(shù)、計(jì)算機(jī)技術(shù)、流體力學(xué)和近代光學(xué)技術(shù)的最新成 果，廣泛應(yīng)用于實(shí)驗(yàn)流體力學(xué)、生物醫(yī)學(xué)、工業(yè)制造等許多領(lǐng)域。PIV將整 個(gè)運(yùn)動(dòng)矢量場(chǎng)作為一個(gè)整體連續(xù)函數(shù)進(jìn)行估計(jì)，其分辨率可以達(dá)到圖像中 的每一個(gè)像素，為微觀、精確地進(jìn)行局部圖像運(yùn)動(dòng)計(jì)算和分析提供了可能， 同時(shí)克服了傳統(tǒng)相關(guān)方法中的一些固有缺陷。PIV系統(tǒng)在實(shí)驗(yàn)流體力學(xué)、 空氣動(dòng)力學(xué)、制造業(yè)、醫(yī)學(xué)、工業(yè)、飛機(jī)制造、水利水電、氣象學(xué)等領(lǐng)域 都有著重要的科學(xué)和經(jīng)濟(jì)價(jià)值。
PIV計(jì)算的方法主要分為基于楊氏干涉條紋法、基于相關(guān) (cross-correlation)的方法和粒子跟蹤測(cè)速(Particle Image Tracing, 簡(jiǎn)稱(chēng)PTV)，后兩者則是目前PIV計(jì)算的主流方法。PIV系統(tǒng)的精度與運(yùn)動(dòng) 估計(jì)算法的性能、實(shí)際應(yīng)用場(chǎng)合和物理器件本身特性有關(guān)，拋開(kāi)后兩者就 算法而言，傳統(tǒng)的算法雖逐漸成熟、實(shí)用，但這些方法本身存在一些先天 的缺陷?；谙嚓P(guān)法的不足歸納如下
1. 由于相關(guān)計(jì)算時(shí)所使用的局部鄰域有限，在前后時(shí)間內(nèi)鄰域里的 粒子出現(xiàn)丟失或出現(xiàn)未知新粒子時(shí)，相關(guān)計(jì)算結(jié)果的峰值不一定對(duì)應(yīng)真實(shí) 的運(yùn)動(dòng)位移。
2. 根據(jù)統(tǒng)計(jì)原理，基于相關(guān)法得到的運(yùn)動(dòng)矢量只是概率上的最大可能位移。
3.當(dāng)粒子運(yùn)動(dòng)的速度很慢時(shí)(位移在亞像素級(jí))或有些區(qū)域粒子密度
過(guò)于集中，當(dāng)采用較大的窗口進(jìn)行相關(guān)搜索時(shí)，會(huì)出現(xiàn)多峰值的現(xiàn)象。為 了克服這個(gè)問(wèn)題常采用尺寸和形狀都不固定的局部窗。
與相關(guān)方法類(lèi)似的PTV是通過(guò)跟蹤單個(gè)粒子來(lái)完成對(duì)運(yùn)動(dòng)速度的估 計(jì)。PTV也存在一些先天的缺陷如當(dāng)跟蹤過(guò)程中粒子消失或別的粒子出 現(xiàn)時(shí)會(huì)帶來(lái)計(jì)算上的不穩(wěn)定；雖然PTV能得到粒子在一段時(shí)間內(nèi)的運(yùn)動(dòng)軌 跡(在可靠跟蹤的前提下)，但PTV對(duì)粒子圖像相對(duì)相關(guān)法要求更苛刻，而 且PTV不能得到100%密度的運(yùn)動(dòng)場(chǎng)，通常需要插值運(yùn)算。
近年來(lái)，基于變分(如光流:物體亮度模式變化的反映，通常指運(yùn)動(dòng)矢 量場(chǎng))的PIV計(jì)算受到人們?cè)絹?lái)越多的關(guān)注，由于該方法是將整個(gè)運(yùn)動(dòng)矢量 場(chǎng)作為一個(gè)整體連續(xù)函數(shù)進(jìn)行估計(jì)，所以從理論上講可以克服上面所提到 的問(wèn)題。絕大多數(shù)光流計(jì)算的分辨率可以達(dá)到圖像中的每一個(gè)像素，這為 微觀、精確地進(jìn)行局部圖像運(yùn)動(dòng)計(jì)算和分析提供了可能。同時(shí)粒子運(yùn)動(dòng)圖 像大多在時(shí)間和空間上都是連續(xù)變化的，這與變分方法中圖像序列的局部 時(shí)空可微假設(shè)在物理意義上是一致的。光流場(chǎng)計(jì)算是通過(guò)圖像序列獲得其 運(yùn)動(dòng)信息的一種方法，是計(jì)算機(jī)視覺(jué)領(lǐng)域一項(xiàng)重要技術(shù)問(wèn)題。光流不僅攜 帶了被觀察物體的運(yùn)動(dòng)信息，還攜帶有被觀察物的三維結(jié)構(gòu)、深度、傳感 器參數(shù)、非剛性物體的局部形變，甚至流體運(yùn)動(dòng)的矢量結(jié)構(gòu)特征等豐富信 息。光流計(jì)算在商業(yè)領(lǐng)域、科學(xué)研究、工業(yè)生產(chǎn)中都有廣泛的應(yīng)用，如 視頻分割、運(yùn)動(dòng)估計(jì)、運(yùn)動(dòng)目標(biāo)跟蹤、場(chǎng)景的三維結(jié)構(gòu)重建、微觀血流運(yùn) 動(dòng)分析。光流計(jì)算的研究真正起始于80年代初Horn和chunck以及Lucas 和Kanade等人的奠基性的工作，在其后的30多年時(shí)間里出現(xiàn)了許多新方 法、新思想。近幾年隨著偏微分方程(PDE)，張量分析，微分幾何等數(shù)學(xué)方 法和計(jì)算手段在圖像分析中的不斷滲透，光流計(jì)算又有了快速發(fā)展，無(wú)論 在計(jì)算精度、可靠性、算法實(shí)時(shí)性等方面都取得了較大的飛躍。
雖然基于光流技術(shù)的PIV計(jì)算有著重要的學(xué)術(shù)和應(yīng)用價(jià)值，但其作為 一個(gè)新興的交叉研究領(lǐng)域還存在許多有待解決及改進(jìn)的地方，歸納如下
1) 對(duì)于PIV系統(tǒng)來(lái)講運(yùn)動(dòng)矢量計(jì)算需要結(jié)合實(shí)際的流體力學(xué)物理意 義，如散度和旋度場(chǎng)的數(shù)據(jù)約束等。另外除矢量場(chǎng)計(jì)算外，如何精確計(jì)算 無(wú)散場(chǎng)、無(wú)旋場(chǎng)、粘性區(qū)域、提取拓?fù)潼c(diǎn)(吸收點(diǎn)、發(fā)散點(diǎn))等還存在許多 沒(méi)有解決好的問(wèn)題。
2) 許多變分光流計(jì)算需要解決帶有邊值問(wèn)題的橢圓偏微分方程，并最終通過(guò)求解線性或非線性方程組實(shí)現(xiàn)，其需要大量的迭代計(jì)算。如何提 高光流的計(jì)算速度使之達(dá)到實(shí)時(shí)性的要求，一直是一個(gè)具有挑戰(zhàn)性的課題。
3)對(duì)于流體運(yùn)動(dòng)圖像，除運(yùn)動(dòng)矢量本身外矢量的散度和旋度場(chǎng)也是 重要的計(jì)算與分析對(duì)象。運(yùn)動(dòng)矢量場(chǎng)的散度和旋度常常被用來(lái)計(jì)算與分析
流體矢量場(chǎng)的結(jié)構(gòu)特征和進(jìn)行流場(chǎng)可視化(如分析擴(kuò)散源、匯聚點(diǎn)、渦旋， 以及進(jìn)行流場(chǎng)的拓?fù)浣Y(jié)構(gòu)分割等)。通常，會(huì)將這些特征單獨(dú)分析處理再融 合判斷，而如何基于高維圖像的觀點(diǎn)對(duì)流體運(yùn)動(dòng)場(chǎng)進(jìn)行綜合描述和分析仍 是一個(gè)空白。

發(fā)明內(nèi)容
本發(fā)明的主要目的就是針對(duì)現(xiàn)有技術(shù)的不足，提供一種精度和可靠性 高的粒子圖像測(cè)速處理方法。
本發(fā)明還提供一種能進(jìn)一步有效加快處理過(guò)程的粒子圖像測(cè)速處理方法。
為實(shí)現(xiàn)上述目的，本發(fā)明采用以下技術(shù)方案 一種粒子圖像測(cè)速處理方法，包括以下步驟
A、 對(duì)于NXM像素范圍區(qū)域Q構(gòu)建水平集函數(shù)d)，水平集函數(shù)(fr為定義 在NXM像素范圍內(nèi)的連續(xù)函數(shù)，用水平集函數(shù)())將區(qū)域Q劃分為第一子區(qū) 域^和第二子區(qū)域02，滿(mǎn)足運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)在第一子區(qū)域Q,和 第二子區(qū)域02內(nèi)光滑連續(xù)，但在交界處不連續(xù)，其中，(j^O的區(qū)域?qū)?yīng) 第一子區(qū)域Q,， (|)<0的區(qū)域?qū)?yīng)第二子區(qū)域02，
其中，第一子區(qū)域q對(duì)應(yīng)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)[V,A，"」，第二 子區(qū)域02區(qū)域?qū)?yīng)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)[V2，^，r!2]，運(yùn)動(dòng)矢量
V=[vx,vy]，散度卜div(V^^ + ^，旋度^curl(V)-^k—、,、為 l 」 x y y x y
沿x軸和y軸方向的兩個(gè)速度分量；
構(gòu)建指示函數(shù)H，定義為H(x)=l,x20and H(x)-0,x〈0;
B、 構(gòu)造區(qū)域Q的能量函數(shù)E-E,+XE2，其中E為全局能量泛函，Ei為 數(shù)據(jù)約束能量，E2光滑約束能量；
C、 最小化全局能量泛函，實(shí)現(xiàn)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)的最佳逼近。 優(yōu)選地，數(shù)據(jù)約束能量和光滑約束能量分別采用如下形式
6^ = Jll(X+V,，t) —I(X，t)|2H(*) + |l(X + V2，t)-I(X，t)|2HH))dX + J]div(VJ-^12 HW + |div(V2)-g2 HH))dx + J|Curl(V1)-Ti1|2H(<|)) + |Curl(V2)-"2fH(-<|))dX
E2 二aflVVj2H(<|)) + |VV2|2H(-()))dx +ajKH((|)) + |vg2HH))dX
+ccJ|V"1|2H(())) + |Vt12|2H(—<|))dX
n
+(5jlVH(<l))l2dX
優(yōu)選地，全局能量泛函的最小化通過(guò)求解以下線性方程組實(shí)現(xiàn)
AU = F
U-[^V2,^，^,r^， ,(1)]， A為系數(shù)矩陣，F(xiàn)為向量，表示線性方程組的 奇次項(xiàng)。
可以通過(guò)經(jīng)典的變分計(jì)算方法獲得和處理線性方程組AU = F 。 優(yōu)選通過(guò)多網(wǎng)格方法處理線性方程組AU-F 。 所述處理優(yōu)選地包括
設(shè)置第一網(wǎng)格，獲得方程組AU-F對(duì)應(yīng)在第一網(wǎng)格h下的方程組 AhUh = Fh ;
迭代求解方程組AhUh-Fh ，根據(jù)所得實(shí)際解tJh得到方程誤差 Rh=Fh_Ah^jh>
設(shè)置尺寸較第一網(wǎng)格更粗的第二網(wǎng)格H，將在第一網(wǎng)格h下對(duì)方程組 AV = Rh的求解尺度映射為在第二網(wǎng)格H下對(duì)方程組AHeH = RH的求解； 迭代求解在第二網(wǎng)格H下的高頻計(jì)算誤差eH ;
對(duì)第二網(wǎng)格H下的高頻計(jì)算誤差eH進(jìn)行尺度映射，得到在第一網(wǎng)格h 下的低頻計(jì)算誤差eN和
用第一網(wǎng)格h下的低頻計(jì)算誤差eh更新實(shí)際解ti;;ew =eh + t!h 。 所述處理還可以包括
在更新實(shí)際解后，再在第一網(wǎng)格h下迭代求解AhUh-Fh，來(lái)消除G^ 的高頻計(jì)算誤差。
優(yōu)選地，循環(huán)多次進(jìn)行所述處理，其中第二網(wǎng)格的尺度逐次遞增。本發(fā)明有益的技術(shù)效果是-
與灰度圖像類(lèi)似，矢量場(chǎng)的結(jié)構(gòu)信息表現(xiàn)為局部矢量場(chǎng)的邊緣、多矢 量交匯、分叉等結(jié)構(gòu)特點(diǎn)，即可以看作運(yùn)動(dòng)估計(jì)中的多運(yùn)動(dòng)和運(yùn)動(dòng)邊界。
本發(fā)明粒子圖像測(cè)速處理方法，對(duì)于2維PIV運(yùn)動(dòng)矢量場(chǎng)、散度和旋度場(chǎng) 獲得，采用了基于光流運(yùn)動(dòng)方程的計(jì)算思想，以整體函數(shù)的形式計(jì)算運(yùn)動(dòng) 矢量，并通過(guò)構(gòu)建水平集(levelset)函數(shù)，將水平集的策略和思想集成于 現(xiàn)有的PIV計(jì)算框架，解決了運(yùn)動(dòng)邊界和多運(yùn)動(dòng)問(wèn)題，同時(shí)提高了流體運(yùn) 動(dòng)特征描述的精度，在抑止噪聲同時(shí)又能保持矢量場(chǎng)、散度和旋度場(chǎng)的局 部結(jié)構(gòu)信息，精度高，可靠性好。
進(jìn)一步地，在圖像計(jì)算中采用多網(wǎng)格方法進(jìn)行處理，可以有效加快線 性方程組迭代的收斂速度，大大加快圖像處理速度，同時(shí)不會(huì)降低圖像處 理的可靠性和精度。


圖1為粒子圖像測(cè)速處理系統(tǒng)的結(jié)構(gòu)框圖2為本發(fā)明粒子圖像測(cè)速處理方法一種實(shí)施例的流程圖3為用水平集函數(shù)^來(lái)劃分內(nèi)部光滑連續(xù)而邊界不連續(xù)的場(chǎng)的示意
圖4為一種實(shí)施例中利用多網(wǎng)格進(jìn)行加速處理的流程圖； 圖5a和圖5b分別為V型和W型多網(wǎng)格結(jié)構(gòu)的示意圖； 圖6為一種優(yōu)選實(shí)施例的多網(wǎng)格結(jié)構(gòu)的示意圖。 本發(fā)明的特征及優(yōu)點(diǎn)將通過(guò)實(shí)施例結(jié)合附圖進(jìn)行詳細(xì)說(shuō)明。
具體實(shí)施例方式
一個(gè)完整的粒子圖像測(cè)速處理系統(tǒng)其關(guān)鍵部分如圖l所示，主要包括 數(shù)據(jù)采集模塊、圖像預(yù)處理模塊、運(yùn)動(dòng)矢量場(chǎng)(光流場(chǎng))計(jì)算模塊和后續(xù)數(shù) 據(jù)分析模塊，后續(xù)數(shù)據(jù)分析模塊包括旋度計(jì)算模塊和散度計(jì)算模塊。其中 的運(yùn)動(dòng)矢量場(chǎng)計(jì)算模塊是整個(gè)PIV系統(tǒng)的核心模塊。
一種粒子圖像測(cè)速處理方法請(qǐng)參考圖2。在計(jì)算2維運(yùn)動(dòng)矢量場(chǎng)、散 度和旋度場(chǎng)時(shí)，采用了基于光流運(yùn)動(dòng)方程的計(jì)算思想，以整體函數(shù)的形式 計(jì)算運(yùn)動(dòng)矢量。為了能夠抑止噪聲同時(shí)又能保持矢量場(chǎng)、散度和旋度場(chǎng)的 局部結(jié)構(gòu)信息，還結(jié)合了水平集的思想。
2維運(yùn)動(dòng)矢量場(chǎng)是指在大小為NXM像素范圍內(nèi)(如通常熟知的PAL制-768 X 576或640 X 480個(gè)像素)能夠描述每一像素位置(NX M內(nèi))下該像素當(dāng) 前時(shí)刻的位移矢量，該矢量為一浮點(diǎn)數(shù)二維矢量，表示為V = [vx,Vy]，
8其中v.,vy為x軸和y軸方向的兩個(gè)速度分量。為了統(tǒng)一描述，用大寫(xiě)粗體
表示矢量，小寫(xiě)細(xì)體表示標(biāo)量。速度場(chǎng)V實(shí)際是坐標(biāo)位置變量X-[x,y]的函數(shù)。 流體運(yùn)動(dòng)圖像計(jì)算出了獲得瞬時(shí)運(yùn)動(dòng)矢量場(chǎng)V外，還需計(jì)算V的散度場(chǎng)g和
旋度場(chǎng)n，如(1)(2)式所示。<formula>formula see original document page 9</formula>由于實(shí)際的流體運(yùn)動(dòng)可能會(huì)有多矢量流交界的現(xiàn)象，即運(yùn)動(dòng)邊界和多 運(yùn)動(dòng)問(wèn)題，為了解決運(yùn)動(dòng)邊界和多運(yùn)動(dòng)問(wèn)題，同時(shí)為了提高流體運(yùn)動(dòng)特征 描述的精度，將水平集的策略和思想集成于原有的光流計(jì)算框架。
如圖3所示，在某一區(qū)域D內(nèi)運(yùn)動(dòng)矢量、散度和旋度有可能不連續(xù)， 為此可以將Q劃分為兩個(gè)子區(qū)域Q,A， ^區(qū)域內(nèi)對(duì)應(yīng)的運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和 旋度場(chǎng)為[Hn,；Q:區(qū)域內(nèi)對(duì)應(yīng)的運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)為[V2,^,"2]。
假設(shè)在區(qū)域n,內(nèi)[v",，n,光滑連續(xù)，在區(qū)域^內(nèi)[V2，^，nd光滑連續(xù)，但在n,A 的交界處運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)并不連續(xù)。為了保證計(jì)算時(shí)[v,，^n、l和和^2，《2,112]時(shí)同時(shí) 要不斷計(jì)算更新^。
為了估計(jì)出[Hr^和[v"^]，構(gòu)造以下能量函數(shù)。<formula>formula see original document page 9</formula>
其中E為全局能量泛函，E,，E,分別為以下式(4)、 (5)形式的數(shù)據(jù)約束 能量和光滑約束能量。<formula>formula see original document page 9</formula>E2 = a J|VV! |2 H④+1W212 H(-小)dx o
+a<[|V^|2H(<|)) + |V^|2H(—小)dx
n (5) +ocJlVri1|2H((t)) + |VTl2|2HH))dX
ft
+|5j|VH((|))|2dX
這里i(x,t)代表粒子圖像序列。最終估計(jì)的光流場(chǎng)可以描述為-
]_V2 else
通過(guò)全局能量泛函的最小化來(lái)實(shí)現(xiàn)光流場(chǎng)的最佳逼近。
能量泛函(3)的最小化可以通過(guò)求解線性方程組(6)實(shí)現(xiàn)，。
AU = F (6) 其中11 = ^1，^，^,^，111,112,小]，A為9X9的系數(shù)矩陣，F(xiàn)為9X1的向量
表示線性方程的奇次項(xiàng)。
線性方程組(6)的獲得和計(jì)算可以通過(guò)經(jīng)典的變分計(jì)算方法實(shí)現(xiàn)，這 里不詳細(xì)闡述。
在優(yōu)選的實(shí)施例中，可以將基于多網(wǎng)格加速策略與實(shí)際的PIV計(jì)算相 結(jié)合，通過(guò)多網(wǎng)格方法來(lái)處理線性方程組(6)，加速PIV計(jì)算。
雖然光流計(jì)算無(wú)論在計(jì)算精度、可靠性等方面都取得了較大的飛躍， 但其數(shù)值方法和手段并沒(méi)有重大突破。已在工程數(shù)學(xué)和有限元數(shù)值計(jì)算中 廣泛使用的多網(wǎng)格方法，憑借其優(yōu)越的加速收斂特性開(kāi)始得到圖像處理領(lǐng) 域人員的關(guān)注。
多網(wǎng)格計(jì)算的一個(gè)中心任務(wù)是求解各種各樣來(lái)自實(shí)際問(wèn)題的偏微分 方程及線性方程組。由于許多圖像分析問(wèn)題需要解決帶有邊值問(wèn)題的橢圓 偏微分方程，并最終通過(guò)求解線性或非線性方程組實(shí)現(xiàn)，如基于變分方 法的圖像非線性擴(kuò)散、主動(dòng)輪廓計(jì)算、圖像恢復(fù)等，在大尺寸網(wǎng)格上線性 方程的高頻誤差很容易消除，相反即使經(jīng)過(guò)大量迭代低頻誤差卻很難消除， 但在小尺度網(wǎng)格上這種低頻誤差只需經(jīng)少量迭代就能明顯消除。因此多網(wǎng) 格策略可通過(guò)不同網(wǎng)格的迭代計(jì)算來(lái)加速收斂過(guò)程。對(duì)于圖像分析而言， 網(wǎng)格的定義比較簡(jiǎn)單，原始圖像的每一個(gè)像素表示一個(gè)網(wǎng)格，即最細(xì)的網(wǎng) 格，計(jì)算都是基于這樣的離散網(wǎng)格點(diǎn)上的值運(yùn)算完成的。
請(qǐng)參考圖4， 一種優(yōu)選的實(shí)施例的粒子圖像測(cè)速處理還包括如下步驟:可以將式(6)寫(xiě)為式(7):
AhUh=Fh (7) 其中指標(biāo)h代表了網(wǎng)格尺寸為、xhy的第一網(wǎng)格，Ah，TJh分別表示在網(wǎng)
格尺寸為h^hy (不同網(wǎng)格尺寸對(duì)應(yīng)不同的圖像分辨率)下的方程系數(shù)矩
陣和解。
假設(shè)方程(7)經(jīng)過(guò)A次迭代后得到的解為tJh，并令lJh為理想的真實(shí)解，
則運(yùn)動(dòng)矢量的計(jì)算誤差為
eh=Uh-t!h (8) 由于計(jì)算結(jié)果的高頻誤差經(jīng)過(guò)少量迭代就可以基本消除，因此eh在這 里主要包含計(jì)算誤差的低頻能量。如果能知道^則可以對(duì)tJh進(jìn)行補(bǔ)償達(dá)到 減少計(jì)算誤差的目的。雖然直接得到eh并不容易，但可以通過(guò)式(9)和式(10) 進(jìn)行間接求解。
通過(guò)已知的!Jh我們可以計(jì)算出方程(7)的誤差Rh ，有式(9)。
Rh=Fh-Ah(Jh (9)
由于Ah是線性的因此有
AV=Rh (10) 這里的eh主要是在網(wǎng)格h下的低頻計(jì)算誤差。根據(jù)前面的分析，知道 僅從方程(10)很難得到低頻誤差eh。多網(wǎng)格計(jì)算的巧妙之處就在于^不是 在原先網(wǎng)格h下計(jì)算得到的，而是在更粗的網(wǎng)格上計(jì)算得到的。
對(duì)于較粗網(wǎng)格H的尺寸為H^Hy，并有H^^h,，Hy〉hy，在粗網(wǎng)格上可
以將(10)寫(xiě)為:
AHeH=RH (11) 其中AH，RH為Ah，Rh經(jīng)過(guò)尺度映射后得到的新方程系數(shù)，映射關(guān)系如
下
AH=Ph—H(Ah)，RH=Ph—H(Rh) P"h為映射函數(shù)，eH反映的是網(wǎng)格H下的高頻計(jì)算誤差，同時(shí)又是 細(xì)網(wǎng)格h下的低頻計(jì)算誤差。由于eH和eb不在一個(gè)尺度上，因此還需要將eh 映射為eH，如下
eh=Ph—H(eH)
在網(wǎng)格H下的高頻計(jì)算誤差eH只需通過(guò)很少的迭代就能得到，通過(guò)映 射PH—h (□)最終可以得到eh ，有了 eh就可以通過(guò)(12)來(lái)更新tlh 。
D;;ew=eh+tJh (12)
pew由于尺度之間的映射也會(huì)帶來(lái)一定的高頻誤差，因此最后還可以在網(wǎng) 格h下再進(jìn)行n,次迭代來(lái)消除C匕的高頻誤差。
可將多網(wǎng)格結(jié)構(gòu)與圖像計(jì)算中的多尺度思想相結(jié)合應(yīng)用于PIV計(jì)算。 采用多尺度(有可能尺度因子小于2)逐級(jí)逼近的策略，可以降低方程模型 誤差帶來(lái)的影響，改善計(jì)算精度。每次尺度變化都會(huì)進(jìn)行一次基于當(dāng)前估 計(jì)結(jié)果的補(bǔ)償更新。將圖像尺度的變化看作是不同尺寸網(wǎng)格的變化，本質(zhì) 上多尺度逼近計(jì)算可以看作是多網(wǎng)格算法的一種特例。
因此，根據(jù)實(shí)際計(jì)算需要也可以采用更多尺寸的網(wǎng)格。多網(wǎng)格算法在 實(shí)際中通常分為三類(lèi)V型多網(wǎng)格計(jì)算、W型多網(wǎng)格計(jì)算和完全型多網(wǎng)格計(jì) 算，圖5a和圖5b給出了 V型和W型多網(wǎng)格計(jì)算的示意圖。
優(yōu)選地，采用如圖6所示的多網(wǎng)格結(jié)構(gòu)。圖6中的圓點(diǎn)代表不同尺度 的網(wǎng)格(從下至上網(wǎng)格尺寸依次增大)，上升箭頭表示從小尺度網(wǎng)格向大尺 度網(wǎng)格過(guò)渡，向下箭頭則相反。每一完成一次V型多網(wǎng)格計(jì)算后，便將式 (7)中的方程系數(shù)A和F進(jìn)行更新，進(jìn)行多次循環(huán)處理。通過(guò)多網(wǎng)格計(jì)算可 以有效加快線性方程組迭代的收斂速度，同時(shí)不會(huì)降低計(jì)算的可靠性和精 度。
以上內(nèi)容是結(jié)合具體的優(yōu)選實(shí)施方式對(duì)本發(fā)明所作的進(jìn)一步詳細(xì)說(shuō) 明，不能認(rèn)定本發(fā)明的具體實(shí)施只局限于這些說(shuō)明。對(duì)于本發(fā)明所屬技術(shù) 領(lǐng)域的普通技術(shù)人員來(lái)說(shuō)，在不脫離本發(fā)明構(gòu)思的前提下，還可以做出若 干簡(jiǎn)單推演或替換，都應(yīng)當(dāng)視為屬于本發(fā)明的保護(hù)范圍。
1權(quán)利要求
1.一種粒子圖像測(cè)速處理方法，其特征在于包括以下步驟A、對(duì)于N×M像素范圍區(qū)域Ω，構(gòu)建水平集函數(shù)φ，水平集函數(shù)φ為定義在N×M像素范圍內(nèi)的連續(xù)函數(shù)，用水平集函數(shù)φ將區(qū)域Ω劃分為第一子區(qū)域Ω1和第二子區(qū)域Ω2，滿(mǎn)足運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)在第一子區(qū)域Ω1和第二子區(qū)域Ω2內(nèi)光滑連續(xù)，但在交界處不連續(xù)，其中，φ≥0的區(qū)域?qū)?yīng)第一子區(qū)域Ω1，φ＜0的區(qū)域?qū)?yīng)第二子區(qū)域Ω2；其中，第一子區(qū)域Ω1對(duì)應(yīng)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)為[V1，ξ1，η1]，第二子區(qū)域Ω2區(qū)域?qū)?yīng)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)[V2，ξ2，η2]，運(yùn)動(dòng)矢量V＝[vx，vy]，散度<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><mi>&xi;</mi><mo>=</mo><mi>div</mi><mrow> <mo>(</mo> <mi>V</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <msub><mrow> <mo>&PartialD;</mo> <mi>v</mi></mrow><mi>x</mi> </msub> <mi>x</mi></mfrac><mo>+</mo><mfrac> <msub><mrow> <mo>&PartialD;</mo> <mi>v</mi></mrow><mi>y</mi> </msub> <mi>y</mi></mfrac><mo>,</mo> </mrow>]]></math> id="icf0001" file="A2009101094300002C1.tif" wi="42" he="10" top= "95" left = "62" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>旋度<maths id="math0002" num="0002" ><math><![CDATA[ <mrow><mi>&eta;</mi><mo>=</mo><mi>curl</mi><mrow> <mo>(</mo> <mi>V</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <msub><mrow> <mo>&PartialD;</mo> <mi>v</mi></mrow><mi>x</mi> </msub> <mi>y</mi></mfrac><mo>-</mo><mfrac> <msub><mrow> <mo>&PartialD;</mo> <mi>v</mi></mrow><mi>y</mi> </msub> <mi>x</mi></mfrac><mo>,</mo> </mrow>]]></math> id="icf0002" file="A2009101094300002C2.tif" wi="43" he="10" top= "95" left = "118" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>vx，vy為沿x軸和y軸方向的兩個(gè)速度分量；構(gòu)建指示函數(shù)H，定義為H(x)＝1，x≥0 and H(x)＝0，x＜0；B、構(gòu)造區(qū)域Ω的能量函數(shù)E＝E1+λE2，其中E為全局能量泛函，E1為數(shù)據(jù)約束能量，E2光滑約束能量；C、最小化全局能量泛函，實(shí)現(xiàn)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)的最佳逼近。
2. 如權(quán)利要求l所述的粒子圖像測(cè)速處理方法，其特征在于，數(shù)據(jù)約束能量和光滑約束能量分別采用如下形式E, = J|I(X+V, ， t) - I(X, t)|2 H(小)+ |I(X+V2, t) - I(X， t)|2 HH))dx+Jldiv(V,)-g2 HW + |div(V2H2|2 HH))dx +》,,)-1^2 H(<|))+|Curl(V2)-"2|2 H(-小)dx E2 = a J| |2 H(小)+ |W212 H(-小)dx +aJ|V^|2 2|2HH>)dX +a J] Vt^ |2 h(小)+ |Vti2 |2 h(-小)dx +pJ|VH((j))|2dX 。
3. 如權(quán)利要求1所述的粒子圖像測(cè)速處理方法，其特征在于，全局 能量泛函的最小化通過(guò)求解以下線性方程組實(shí)現(xiàn)AU = FU4V，V2,^，^,T^, ,(^， A為系數(shù)矩陣，F(xiàn)為向量，表示線性方程組的奇次項(xiàng)。
4. 如權(quán)利要求3所述的粒子圖像測(cè)速處理方法,其特征在于，通過(guò) 經(jīng)典的變分計(jì)算方法獲得和處理線性方程組AU = F 。
5. 如權(quán)利要求3所述的粒子圖像測(cè)速處理方法，其特征在于，通過(guò) 多網(wǎng)格方法處理線性方程組AU = F 。
6. 如權(quán)利要求5所述的粒子圖像測(cè)速處理方法，其特征在于，所述 處理包括設(shè)置分辨率較高的第一網(wǎng)格，獲得方程組AU-F對(duì)應(yīng)在第一網(wǎng)格h 下的方程組AhU^Fh;迭代求解方程組AhlJh-Fh ，根據(jù)所得實(shí)際解tJh得到方程誤差Rh = ph _ Ah^jh ，設(shè)置尺寸較第一網(wǎng)格更粗的第二網(wǎng)格H，將在第一網(wǎng)格h下對(duì)方程組AV = Rh的求解尺度映射為在第二網(wǎng)格H下對(duì)方程組AHeH = RH的求解； 迭代求解在第二網(wǎng)格H下的高頻計(jì)算誤差eH;對(duì)第二網(wǎng)格H下的高頻計(jì)算誤差eH進(jìn)行尺度映射，得到在第一網(wǎng)格h 下的低頻計(jì)算誤差eh;和用第一網(wǎng)格h下的低頻計(jì)算誤差eh更新實(shí)際解t!hnew =eh +Ch 。
7. 如權(quán)利要求6所述的粒子圖像測(cè)速處理方法，其特征在于,所述處 理還包括在更新實(shí)際解后，再在第一網(wǎng)格h下迭代求解AhUh-Fh，來(lái)消除tJ^ 的高頻計(jì)算誤差。
8. 如權(quán)利要求6或7所述的粒子圖像測(cè)速處理方法，其特征在于,循 環(huán)多次進(jìn)行所述處理，其中第二網(wǎng)格的尺度逐次遞增。
全文摘要
本發(fā)明公開(kāi)了一種粒子圖像測(cè)速處理方法，包括以下步驟A.對(duì)于N×M像素范圍Ω構(gòu)建水平集函數(shù)φ，水平集函數(shù)φ為定義在N×M像素范圍內(nèi)的連續(xù)函數(shù)，用水平集函數(shù)φ將區(qū)域Ω劃分為第一子區(qū)域Ω₁和第二子區(qū)域Ω₂，運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)在兩個(gè)子區(qū)域內(nèi)光滑連續(xù)，但在交界處不連續(xù)，其中，φ≥0的區(qū)域?qū)?yīng)第一子區(qū)域Ω₁，φ＜0的區(qū)域?qū)?yīng)第二子區(qū)域Ω₂；B.構(gòu)造區(qū)域Ω的能量函數(shù)E＝E₁+λE₂，其中E為全局能量泛函，E₁為數(shù)據(jù)約束能量，E₂光滑約束能量；C.最小化全局能量泛函，實(shí)現(xiàn)運(yùn)動(dòng)場(chǎng)、散度場(chǎng)和旋度場(chǎng)的最佳逼近。本發(fā)明提供了一種精度高可靠性好的PIV處理方法。
文檔編號(hào)G01P5/18GK101629966SQ200910109430
公開(kāi)日2010年1月20日 申請(qǐng)日期2009年8月18日 優(yōu)先權(quán)日2009年8月18日
發(fā)明者盧宗慶, 廖慶敏 申請(qǐng)人:清華大學(xué)深圳研究生院