本發(fā)明屬于流場觀測
技術領域：
，特別涉及一種基于單光場相機的三維流場測試方法。
背景技術：
：專利申請CN201210260127.8，公開了“一種單CCD相機三維粒子圖像測速方法，單CCD相機三維粒子圖像測速技術是基于半凸透鏡成像原理，在PIV系統(tǒng)中示蹤粒子的散射光線通過兩上下并排放置的半凸透鏡，在像平面或CCD上成兩個實像；兩個像的距離與物距以及兩半凸透鏡的光軸距離成函數(shù)關系；示蹤粒子在景深內無需改變像距便可得到清晰的像；在PIV系統(tǒng)應用中BB1為第一幅照片內同一示蹤粒子所成的像之間的距離，CC1為第二幅照片內同一示蹤粒子像的距離；采用兩半凸透鏡做相機的鏡頭，兩半凸透鏡間填充不透光的密封介質；示蹤粒子運動的二維速度，采用當今通用的互相關算法得到”。但是該方法實際使用中操作復雜，計算精度也不能滿足需要。本發(fā)明涉及的參考文獻：[1]RafaelC.Gonzalez,RichardE.Woods.DigitalImageProcessing(3rdEdition),PearsonEducation,2011[2]HermanGT,LentA(1976)Iterativereconstructionalgorithms.ComputersinBiologyandMedicine,6:273–294[3]Raffel,M.,Willert,C.E.,wereley,s.,Kompenhans,J.,ParticleImageVelocimetryAPracticalGuide,Springer-VerlagBerlinHeidelberg,2007技術實現(xiàn)要素：本發(fā)明的目的是提供一種基于單光場相機的三維流場測試方法。本發(fā)明的技術方案是，一種基于單光場相機的三維流場測試方法，包含以下步驟：A1，利用單光場相機獲取待測流場中示蹤粒子的時序粒子光場圖像；A2，預處理時序粒子光場圖像，去除背景噪聲；A3，采用CPU和/或GPU并行處理方法，對時序粒子光場圖像進行重構，獲得時序三維粒子圖像；A4，對時序三維粒子圖像進行互相關計算，獲得三維流場速度場分布；A5，對三維流場速度場進行后處理，剔除錯誤速度矢量，同時根據(jù)相鄰正確矢量插值替代被剔除的錯誤矢量。步驟A1中，在待測流場中撒布示蹤粒子，用高亮度光源提供體式照明，然后利用高分辨或者高速圖像傳感器與微透鏡陣列(MLA)的精密封裝體形成的單光場相機拍攝示蹤粒子的時序光場圖像。步驟A2中，對所拍攝的時序粒子光場圖像進行預處理，其中去除背景噪聲的降噪運算采用全局閾值、局部閾值、高斯光滑濾波、滑動最小值濾波算法的一種或多種算法組合。步驟A3中，所述的重構包括：根據(jù)幾何光學計算空間光束經主透鏡、微透鏡陣列后到達圖像傳感器的準確位置，對空間每一個體素(voxel)所發(fā)出的光束，計算其與相交微透鏡單元的重疊面積，以及其與相交CCD/CMOS像素的重疊面積，從而計算獲得空間每個體素與對應像素的權重系數(shù)Wi，jYCCD=VB-(sy+Yl)sifl+sy=(VB-sy)flfm(1-M)+sy-flfm(1-M)Yl---(1)]]>式中各變量的含義如下：YCCD為光線與圖像傳感器相交點的坐標；VB為光線與主鏡頭相交點的坐標；sy為微透鏡陣列中某一微透鏡單元相對主光軸的偏移值；Yl為光線與某一微透鏡單元相交點的坐標；si為主透鏡與微透鏡陣列之間的距離；fl為微透鏡焦距；fm為主透鏡焦距；pp為像素單元的尺寸；M為放大系數(shù)，利用密集光線逆追蹤方法，定位空間重構體素(Voxel)所對應的像素集合，以每個像素為單元，根據(jù)式(1)來計算該像素所采集空間光線的位置，即反向追蹤光線，以建立起像素與空間某一點光源所發(fā)出光線的對應關系，利用流場示蹤粒子的稀疏性，其對應所需重構的體素區(qū)域也為三維稀疏矩陣，將空間某一體素所對應的像素值相乘，如果乘積非零或者大于某一閾值，則表示該體素處可能存在一個示蹤粒子，可以對其進行后續(xù)的重構計算；反之如果乘積為零或者小于某一閾值，則表示該體素處沒有示蹤粒子，后續(xù)重構計算中可以忽略該體素；根據(jù)乘法代數(shù)重構算法，利用計算得到的權重系數(shù)Wi，j、以及所篩選出的非零體素，根據(jù)所記錄的粒子光場像素矩陣(I(xi，yi))，迭代計算稀疏重構矩陣的體素值(E(Xj，Yj，Zj)k+1)，E(Xj,Yj,Zj)k+1=E(Xj,Yj,Zj)k(I(xi,yi)Σj∈Niwi,jE(Xj,Yj,Zj)k)μwi,j---(2)]]>式中：E(Xj,Yj,Zj)為空間體素(Xj,Yj,Zj)的數(shù)值；上標k表示第k次迭代計算所得到該體素的數(shù)值；I(xi，yi)為像素(xi，yi)處的數(shù)值，該數(shù)值由光場相機所拍攝的示蹤粒子光場圖像獲得；wi,j為體素j與像素i所對應的權重系數(shù)；μ為MART計算參數(shù)。所述步驟A4中，將相鄰時刻的示蹤粒子三維體素矩陣按照式(3)進行三維互相關計算，獲取流場三維速度場分布，Φ(m,n,l)=Σi=1MΣj=1NΣk=1LE1(i,j,k)·E2(i+m,j+n,k+l)---(3).]]>步驟A5中，采用全局閾值計算式(4)、局部中值濾波式(5)或者局部平均值濾波式(6)對獲得的三維流場U(i，j，k)，挑選并剔除出其中的錯誤矢量，并采用線性插值或者三次樣條插值的方法，根據(jù)正確速度矢量插值獲得替代矢量，U(i,j,k)>U(i,j,k)‾+T*STD(U(i,j,k))U(i,j,k)<U(i,j,k)‾+T*STD(U(i,j,k))---(4)]]>U(i,j,k)>U(i-1:i+1,j-1:j+1,k-1:k+1)‾+T*STD(U(i-1:i+1,j-1:j+1,k-1:k+1))U(i,j,k)<U(i-1:i+1,j-1:j+1,k-1:k+1)‾+T*STD(U(i-1:i+1,j-1:j+1,k-1:k+1))---(5)]]>其中，是三維速度場的平均值；STD(U(i，j，k))是三維速度場的標準方差；T為濾波閾值；是局部(3×3×3)區(qū)域內的三維速度場平均值；STD(U(i-1：i+1，j-1；j+1，k-1：k+1)是局部(3×3×3)區(qū)域內的三維速度場標準方差；U(i-1：i+1，j-1：j+1，k-1：k+1)是局部(3×3×3)區(qū)域內的三維速度場中值。一種三維流場測試系統(tǒng)，包括由圖像傳感器與微透鏡陣列(MLA)的精密封裝體構成的單光場相機，在拍攝示蹤粒子的時序光場圖像時，在示蹤粒子圖像的反面具有高亮度體式照明，該體式照明是由包括光源、透鏡組合和鏡頭的光學系統(tǒng)提供的。本發(fā)明能通過單個光場相機的時序圖像采集而獲取待測流場的三維速度分布，相比于現(xiàn)有多視角三維流場測試方法，極大的減少了硬件系統(tǒng)配置、簡化了硬件系統(tǒng)調節(jié)步驟，特別適用于受限空間下的三維流場測量。附圖說明圖1本發(fā)明的單光場相機三維流場測試系統(tǒng)原理示意圖。圖2本發(fā)明實施例中光線追蹤示意圖。圖3本發(fā)明實施例中權重系數(shù)計算示意圖。圖4本發(fā)明實施例中密集光線逆追蹤方法原理示意圖(DenseRayTracing，DRT)。圖5是本發(fā)明實施例中測試方法的流程圖。具體實施方式本發(fā)明的一種基于單個光場相機的流場三維速度場測試方法，通過單個視角測量獲取流場的三維速度場分布。包含以下步驟：1)在待測流場中撒布示蹤粒子(水中一般使用20～30微米的空心玻璃珠；空氣中一般使用1微米左右的液滴或者幾百納米的二氧化鈦顆粒)，按照圖1的方式，用高亮度光源(激光、LED或者高亮投影儀)提供體式照明；然后利用高分辨或者高速圖像傳感器(CCD或者CMOS)與微透鏡陣列(MLA)的精密封裝體(即光場相機，以下簡稱單光場相機)拍攝示蹤粒子的時序光場圖片。2)對所拍攝的時序粒子光場圖片進行前處理，降噪運算采用全局閾值、局部閾值、高斯光滑濾波、滑動最小值濾波算法的一種或多種算法組合[1]。3)根據(jù)幾何光學計算空間光束經主透鏡、微透鏡陣列后到達圖像傳感器的準確位置(式1，圖2)。按圖3所示，對空間每一個體素(voxel)所發(fā)出的光束，計算其與相交微透鏡單元的重疊面積，以及其與相交CCD/CMOS像素的重疊面積，從而計算獲得空間每個體素與對應像素的權重系數(shù)Wi，jYCCD=VB-(sy+Yl)sifl+sy=(VB-sy)flfm(1-M)+sy-flfm(1-M)Yl---(1)]]>式中各變量的含義如下(參考圖2)：YCCD為光線與圖像傳感器相交點的坐標；VB為光線與主鏡頭相交點的坐標；sy為微透鏡陣列中某一微透鏡單元相對主光軸的偏移值；Yl為光線與某一微透鏡單元相交點的坐標；si為主透鏡與微透鏡陣列之間的距離；fl為微透鏡焦距；fm為主透鏡焦距；pp為像素單元的尺寸；M為放大系數(shù)。4)利用密集光線逆追蹤方法(DenseRayTracing，DRT)，定位空間重構體素(Voxel)所對應的像素集合。密集光線逆追蹤方法采用圖2所示的方式，以每個像素為單元，根據(jù)式1來計算該像素所采集空間光線的位置。即反向追蹤光線，以建立起像素與空間某一點光源所發(fā)出光線的對應關系。因為流場示蹤粒子的稀疏性，其對應所需重構的體素區(qū)域也為三維稀疏矩陣。利用這一特征，將空間某一體素(如圖4中紅色方框標識的Ej體素)所對應的像素值(如圖4中紅色線條標識的像素集合)相乘。如果乘積非零或者大于某一閾值，則表示該體素處可能存在一個示蹤粒子，可以對其進行后續(xù)的重構計算；反之如果乘積為零或者小于某一閾值，則表示該體素處沒有示蹤粒子，后續(xù)重構計算中可以忽略該體素。采用此DRT方法，可以過濾掉稀疏重構矩陣中不含示蹤粒子的體素，提高重構計算效率。5)根據(jù)乘法代數(shù)重構算法(multiplicativealgebraicreconstructiontechnique,MART，式2，[2])，利用第3步所計算的權重系數(shù)Wi，j、以及第4步所篩選出的非零體素，根據(jù)所記錄的粒子光場像素矩陣(I(xi，yi))，迭代計算稀疏重構矩陣的體素值(E(Xj，Yj，Zj)k+1)。一般給定E(Xj，Yj，Zj)0的初始估計值為1，迭代計算5-10次。E(Xj,Yj,Zj)k+1=E(Xj,Yj,Zj)k(I(xi,yi)Σj∈Niwi,jE(Xj,Yj,Zj)k)μwi,j---(2)]]>式中：E(Xj,Yj,Zj)為空間體素(Xj,Yj,Zj)的數(shù)值；上標k表示第k次迭代計算所得到該體素的數(shù)值；I(xi，yi)為像素(xi，yi)處的數(shù)值，該數(shù)值由光場相機所拍攝的示蹤粒子光場圖像獲得；wi,j為體素j與像素i所對應的權重系數(shù)，其值根據(jù)第3)步所描述的方法計算；μ為MART計算參數(shù)[2]。6)對所記錄的示蹤粒子時序光場圖像按照步驟3～5進行三維重構，將相鄰時刻的示蹤粒子三維體素矩陣按照式3進行三維互相關計算[3]，獲取流場三維速度場分布。Φ(m,n,l)=Σi=1MΣj=1NΣk=1LE1(i,j,k)·E2(i+m,j+n,k+l)---(3)]]>7)采用全局閾值(式4)、局部中值濾波(式5)或者局部平均值濾波(式6)[3]對第6步處理獲得的三維流場U(i，j，k)，挑選并剔除出其中的錯誤矢量。并采用線性插值或者三次樣條插值的方法，根據(jù)正確速度矢量插值獲得替代矢量。U(i,j,k)>U(i,j,k)‾+T*STD(U(i,j,k))U(i,j,k)<U(i,j,k)‾+T*STD(U(i,j,k))---(4)]]>U(i,j,k)>U(i-1:i+1,j-1:j+1,k-1:k+1)‾+T*STD(U(i-1:i+1,j-1:j+1,k-1:k+1))U(i,j,k)<U(i-1:i+1,j-1:j+1,k-1:k+1)‾+T*STD(U(i-1:i+1,j-1:j+1,k-1:k+1))---(5)]]>上式中是三維速度場的平均值；STD(U(i，j，k))是三維速度場的標準方差；T為濾波閾值；是局部(3×3×3)區(qū)域內的三維速度場平均值；STD(U(i-1：i+1，j-1；j+1，k-1：k+1))是局部(3×3×3)區(qū)域內的三維速度場標準方差；U(i-1：i+1，j-1：j+1，k-1：k+1)是局部(3×3×3)區(qū)域內的三維速度場中值。當前第1頁1 2 3