梅 迪, 丁俊飛, 施圣賢

(上海交通大學機械與動力工程學院 葉輪機械研究所, 上海 200240)

0 引 言

二維粒子圖像速度測量技術 (Two-dimensional Particle Image Velocimetry，2D-PIV)是20世紀70年代末發(fā)展起來的一種激光流體測速方法，與激光多普勒測速、熱線風速儀等單點測速手段相比，其具有瞬態(tài)、多點、無接觸式的優(yōu)點, 因此經(jīng)過幾十年的不斷完善與發(fā)展，該技術在實驗流體力學、空氣動力學、燃燒學和葉輪機械等諸多領域都有著廣泛的應用[1-3]。然而,由于基礎科學研究和實際工程應用中的諸多流動現(xiàn)象都具有復雜的三維特征，二維流場的速度信息不足以反映流動現(xiàn)象的數(shù)學物理本質(zhì)，因此，如何精確測量三速度分量的三維(3D-3C)速度場成為了近幾十年來PIV領域的一個研究熱點。

利用雙目視覺仿生學的原理，先后有學者提出了使用2臺相機并采用不同布局方式的Stereo PIV技術，該技術首先獲得激光片光源所在平面中的示蹤粒子在2個相機圖片中各自的二維位移速度場，然后聯(lián)立求解該平面中的三維速度場，從而來獲得一個平面內(nèi)的三速度分量(2D-3C)[4]。與此同時也有學者提出了Scanning PIV技術，該技術利用可以旋轉(zhuǎn)的反射鏡裝置，使激光片光源旋轉(zhuǎn)掃描整個待測流場，用單目或者雙目高速相機測量多個切面的速度矢量場，但受限于激光掃描和相機存儲的速度，該技術的最大可測量速度一般小于1m/s[5]。數(shù)字離焦PIV(Defocusing PIV，DPIV)則采用了帶有光闌的鏡頭拍攝示蹤粒子圖像，光闌具有一定的圖案，如帶有3個圓孔的掩模，示蹤粒子離相機聚焦平面距離不同時，由于光闌的作用，粒子離焦圖像的模糊程度不同并且會發(fā)生上下翻轉(zhuǎn)，因此可以計算出粒子的三維空間位置，帶有光闌的單相機即可完成三維流場測速，但DPIV技術的示蹤粒子濃度一般較低[6]。全息PIV(Holographic, HPIV)則利用全息照相技術將示蹤粒子的三維圖像記錄在銀鹽感光膠片或者數(shù)碼相機中，進而通過三維互相關獲得粒子的速度場。利用感光膠片成像時，HPIV具有目前PIV技術中最高的空間分辨率，然而由于該技術中粒子全息相片處理過程的時間和資金成本都較高，而且全息技術本身的光路布置比較復雜，目前這一技術并沒有得到廣泛應用[7-9]。層析PIV(Tomographic PIV，Tomo-PIV)則是一種近十多年逐步發(fā)展起來的三維流動測試技術，該技術一般采用4～8臺數(shù)字相機在不同角度拍攝同一流體區(qū)域，進而利用得到的粒子圖像通過MART重構(gòu)算法和三維互相關算法計算出速度矢量[10-11]，其測量結(jié)果精準，測量區(qū)域較大，粒子濃度高并且空間分辨率較高，因此已經(jīng)被廣泛應用于多個領域的流動測試，但該技術也存在著校準過程復雜、需要的光學窗口較多等不足，因此并不適合光學窗口數(shù)目受限情形的三維流場測量。鑒于Tomo-PIV操作的復雜性，有國內(nèi)學者提出利用一塊三視野的分光鏡來實現(xiàn)單相機的3D-3C測量[12]，由于被拍攝區(qū)域與相機之間有一塊特殊的三視野分光棱鏡，相機CCD芯片被劃分成了3個不同的成像區(qū)域，每個區(qū)域?qū)聚櫫Ｗ拥囊粋€觀察視角，該技術的圖像數(shù)據(jù)處理過程類似于Tomo-PIV，但實驗裝置及操作過程的復雜性得到了極大的降低，不過由于3個視角對應的觀察角總體比較小，因而在沿光軸方向上的空間分辨率不及垂直于光軸方向的空間分辨率。合成孔徑PIV(Synthetic Aperture PIV，SAPIV)是另外一種利用多相機來測量流體速度矢量場的三維測試技術，通常該技術采用了包含至少5臺以上的相機陣列來獲取粒子圖像，并進一步用合成孔徑重聚焦算法獲取多個重聚焦平面的粒子強度[13]，從而獲得三維粒子光照強度場和速度矢量場。與Tomo-PIV相比，SAPIV能夠消除粒子之間的相互遮擋作用，因此可以測量粒子濃度更高的流場，而且其在景深方向的測量范圍能夠達到與垂直光軸方向測量范圍同等大小的水平，這一技術的最大缺陷在于其結(jié)構(gòu)復雜且成本高昂的相機陣列。

不同于上述的三維PIV測試技術，光場PIV(Light-field PIV, LF-PIV)是一種近些年來發(fā)展起來的能夠僅利用單臺相機即可比較精準測量三維速度場的PIV技術[14-15]，其既不需要復雜的光路系統(tǒng)，又不需要多臺相機，因此極大地簡化了試驗操作難度并降低了硬件成本，特別適合光學窗口較少情形下的流體三維速度場測量。LF-PIV的硬件是基于斯坦福大學Ng于2005年利用現(xiàn)有高分辨率相機開發(fā)出來的緊湊型光場相機，其科研成果轉(zhuǎn)化產(chǎn)品為Lytro品牌的光場相機[16]，但由于這種面向普通消費者的相機不具備PIV所需的雙曝光功能，因此無法實際應用到PIV技術中。市場上另外一種比較流行的光場相機是德國Raytrix公司旗下推出的面向工業(yè)界的產(chǎn)品[17]，但由于Raytrix光場相機的相關硬件結(jié)構(gòu)、微透鏡參數(shù)和重構(gòu)算法均對用戶封閉，因此使用者很難利用該相機針對不同條件下的流場測試問題進行自主開發(fā)粒子重構(gòu)算法[18]。

基于上述情況，本文作者所在的課題組和美國學者Brian領導的團隊分別展開了對光場相機的硬件開發(fā)和粒子三維重構(gòu)算法的研究。LF-PIV利用帶有微透鏡陣列(Microlens Array, MLA)和高分率CCD的光場相機記錄下示蹤粒子的4D光場圖像，進一步通過光場重聚焦算法或者光線跟蹤算法重構(gòu)粒子強度場[19-20]。盡管LF-PIV能夠僅憑借單臺光場相機即可實現(xiàn)3D-3C測量，但由于光場相機整體的視角較小，因此重構(gòu)出的粒子強度在沿光軸方向上會被拉長，從而影響了粒子重構(gòu)沿光軸方向的空間分辨率，這一現(xiàn)象在Tomo-PIV中也比較常見，被稱為拉長效應(Elongation Effect)[10]。針對這一問題，已有學者在HPIV中提出利用2臺相機呈光軸相互垂直布置來達到最佳的空間分辨率[21]。參照該想法，本文作者提出了雙光場相機PIV技術，通過利用2個光場相機增大視角來消除或抑制這一效應，由此提高景深空間分辨率以及沿光軸方向的測量精度。

本文闡述了雙光場相機系統(tǒng)，包括雙光場相機PIV系統(tǒng)概述、光場相機成像原理、基于光線跟蹤的數(shù)字合成圖像以及基于光場相機校準的雙光場相機MART三維粒子重構(gòu)算法。利用DNS射流數(shù)據(jù)生成了數(shù)字合成圖像，通過水射流渦環(huán)驗證實驗得到實驗圖像，并根據(jù)這兩者的結(jié)果對雙光場相機和單光場相機進行詳細對比分析。

1 雙光場相機PIV技術

1.1 雙光場相機PIV系統(tǒng)

為了增大相機系統(tǒng)對示蹤粒子的觀察視角以消除或抑制單光場相機粒子重構(gòu)中的拉長效應，本文采用2個光場相機的光軸相互垂直呈90°的布置方式，以達到最大的相機觀察視角。

雙光場相機PIV的工作流程為(如圖1所示)：將示蹤粒子注入待測試流體，這些粒子在Nd:YAG激光照射下因反射作用變得高亮可見，隨后用2臺光軸互相垂直的光場相機同時拍攝同一流場區(qū)域并采用雙曝光模式記錄下示蹤粒子的光場圖像。對獲得的原始圖像減去由于環(huán)境反射光造成的背景噪聲后，利用基于光場相機校準的MART算法重構(gòu)三維粒子強度場，最后由基于FFT算法的多重網(wǎng)格互相關計算出3D-3C速度場。

圖1 雙光場相機PIV系統(tǒng)工作流程示意圖

1.2 光場相機成像原理

與傳統(tǒng)相機相比，光場相機的不同之處在于其不僅能夠記錄進入相機的每一根光線的空間坐標，也能記錄下其相應的傳播方向，即能夠記錄下進入相機的“光場”(Light Field)?！肮鈭觥边@一概念由Arun Gershun于1936年首次提出，他將其定義為在空間中各點處向各個方向傳播的輻射光線的集合[22]。而當光線在自由空間中傳播且能量不發(fā)生損失時，每根光線可以表示為4D函數(shù)的形式[23]：

L=L(u,v,s,t)(1)

式中,(u,v)和(s,t)分別代表某根光線沿傳播方向與2個平行平面的交點坐標(如圖2所示)。這2組坐標不僅包含了某根光線的空間坐標信息，也包含了其傳播方向信息。

圖2 光場4D函數(shù)

光場相機與常規(guī)相機不同，其在主透鏡和CCD感光器件之間存在一個微透鏡陣列，CCD平面處于微透鏡陣列的焦平面上，因此一個點光源的入射光線經(jīng)過主透鏡后，會透過若干個微透鏡投影到它們每個微透鏡后的一個或者若干個像素上[24]。由相機針孔模型可以看出，同一個微透鏡下的不同像素代表著不同的主透鏡位置，即角度信息，而入射光線到達的不同微透鏡則代表著入射點光源的不同空間位置，即空間信息。

基于上述光場相機成像原理，作者所在的課題組自主開發(fā)了適用于PIV測量的高分辨率光場相機[22]，相機實物和得到的原始圖片如圖3所示。

圖3 光場相機實物和原始圖片

由于微透鏡的形狀、尺寸等參數(shù)會對光場相機成像產(chǎn)生很大影響，在進行了前期研究之后，該光場相機采用了六邊形的微透鏡填充鋪滿整個CCD芯片，以提高CCD利用率并平衡光場相機的角度和空間分辨率[25]。

1.3 基于光線跟蹤的光場相機數(shù)字合成圖像

對光場相機，從點光源發(fā)出的光線的傳播過程如圖4所示[24]。圖中dy, dz表示點光源O的相機坐標系坐標，其參考坐標系原點為光場相機聚焦平面與光軸的交點，pm為主透鏡的直徑，fm為主透鏡的焦距，So和Si分別為物距和像距，pl和fl是微透鏡的圓形通光部分的直徑和焦距，pp是正方形像素的邊長，VB,yl和yCCD分別表示該光線在主透鏡、微透鏡陣列和CCD芯片上的坐標，Sy表示其擊中的微透鏡的光心坐標。

圖4 光場相機光線追蹤過程示意圖[24]

由光場理論可知，進入光場相機的每根入射光線都可以用一個四維向量表示，而且根據(jù)矩陣光學，光線傳播過程中的透鏡器件、一定厚度的均勻介質(zhì)對光線的作用都可以用一個矩陣表達[26]，因此從一個點光源發(fā)出的入射光線，經(jīng)過主透鏡、微透鏡后最終到達CCD平面的傳播過程可以視為對一個四維向量進行的連續(xù)線性變換，最終所有連續(xù)變換的結(jié)果可以用式(2)表示：

(2)

式中，x,y是光線的空間坐標，θ,φ是光線的傳播方向，而x′,y′,θ′,φ′為最終傳播到CCD平面的光線相應的4D參數(shù)，矩陣M是一個4行4列的光線傳播矩陣，物理意義是光路中所有的傳播介質(zhì)和光學器件對4D向量的線性變換。光場相機數(shù)字合成圖像的具體過程見文獻[25]。

由上述方法追蹤單個點光源發(fā)出的每一根光線，可將有一定空間體積大小的示蹤粒子離散成發(fā)光強度在空間中遵循高斯分布的若干個點光源，由此獲得示蹤粒子的光場相機數(shù)字合成圖像[25]。假如同一時刻的單個粒子處于離相機聚焦平面不同距離時，或者同一時刻多個粒子處于不同空間位置時，由光線跟蹤得到的數(shù)字合成圖像如圖5所示。

圖5 數(shù)字合成粒子光場圖像

1.4 基于光場相機校準的雙光場相機MART粒子重構(gòu)算法

在Tomo-PIV中，常用的粒子圖像重構(gòu)算法被稱為MART，其數(shù)學原理如式(3)所示[10]。

E(Xj,Yj,Zj)k+1=E(Xj,Yj,Zj)k

式中,E(Xj,Yj,Zj)代表第j個體素的強度值，其初始值設為1；I(xi,yi)是第i個像素；wi,j是權(quán)重系數(shù)，代表從第j個體素發(fā)出的光線中，有多大比例的能量傳遞給了第i個像素；式中的分母代表所有可以影響第i個像素的體素映射到第i個像素位置上的加權(quán)值之和，μ是迭代時的松弛系數(shù)，一般大于等于1。

由式(3)可以看出，實現(xiàn)MART算法的關鍵在于：(1) 精準計算單個體素影響到的像素位置；(2) 精準計算每個體素對受其影響的像素的權(quán)重系數(shù)。對Tomo-PIV,校準方法采用的是計算機視覺中常用的相機針孔模型矩陣映射函數(shù)[27]，而其相應的權(quán)重系數(shù)計算比較簡單，通過圓柱體-立方體假設模型實現(xiàn)[10-11]，由于光場相機的單個像素可以受到不同空間位置的多個體素影響，因此這一假設模型并不適用。作者所在的課題組曾利用矩陣光學提出了密集光線跟蹤及其相應的權(quán)重系數(shù)計算方法[18]，這一方法能夠被應用的前提條件是光路中的光線傳播矩陣參數(shù)可知，同時由光學元器件本身或者其安裝誤差導致的相機成像畸變較小。然而實際成像過程中由于畸變、光學元器件的加工和安裝誤差對光線傳播矩陣的影響不可忽略，而且考慮到對多相機系統(tǒng)，存在未知的多個相機坐標系之間的相互轉(zhuǎn)換關系，作者所在的課題組提出了一種針對光場相機的彌散圓模型校準算法[28]，利用校準步驟得到的結(jié)果可以精準定位單個或多個相機中受到拍攝區(qū)域中某個點光源影響的像素位置。

校準算法原理如圖6所示，當一個在世界坐標系中定義為(X，Y，Z)的點光源位于離聚焦平面一定距離時，從點光源向各個方向發(fā)出的光線進入主透鏡后會投射到微透鏡陣列上，形成一個圓斑并且覆蓋多個微透鏡，該圓斑被稱為彌散圓，其中心被稱為彌散圓中心[28]。

(a) 點光源彌散圓的形成過程示意圖

(b) 點光源對應的彌散圓中心及受其影響的像素

由于六邊形微透鏡的通光部分也是圓形的，因此進入某個微透鏡的光線投射到CCD后點亮的區(qū)域也會形成一個比彌散圓更小的圓斑，微透鏡下面會有一個或者若干個像素被點亮，其數(shù)量取決于該圓斑的尺寸和位置；彌散圓中心與主透鏡光心、點光源三點共線。

前期研究發(fā)現(xiàn)，點光源的彌散圓中心、受影響像素對應的微透鏡下的圓斑中心、點光源的彌散圓直徑存在著如式(4)所示的關系式[28]：

(4)

根據(jù)上述光場相機成像模型，根據(jù)作者所在課題組提出的光場相機校準算法，能精準計算某一個點光源對應的彌散圓中心、彌散圓直徑以及微透鏡中心像素坐標，從而精確計算出該點光源影響到的所有像素位置，其具體校準過程見文獻[28]。

計算出受影響像素的坐標以后，進一步對每個受影響像素的權(quán)重系數(shù)計算，原理如圖7所示。對一個特定像素，其權(quán)重系數(shù)由兩部分相乘得到：點光源發(fā)出的所有光束中擊中某個微透鏡的光束面積與總的光束面積之比為權(quán)重系數(shù)w1；被擊中的微透鏡下一個CCD像素上的光束面積與進入這個微透鏡的光束面積之比為w2。

圖7 權(quán)重系數(shù)計算原理示意圖

w1,w2都可以采用亞網(wǎng)格計數(shù)的方法計算得到，即對一個像素進一步細分為多個亞網(wǎng)格，計算某個像素內(nèi)有多少個亞網(wǎng)格落入微透鏡下的圓斑之中以及該微透鏡下落入圓斑的亞網(wǎng)格的總數(shù)，前者除以后者即是該像素的w2，同樣的方法可以計算出w1。

基于上述的校準、權(quán)重系數(shù)計算和MART粒子重構(gòu)原理，最終雙相機系統(tǒng)的整個校準及粒子重構(gòu)算法過程如圖8所示，該算法亦可拓展應用到多臺光場相機的粒子三維光照強度場重構(gòu)算法之中。

2 數(shù)字合成圖像對比仿真分析

為了驗證雙光場相機能夠消除或抑制單光場相機PIV中的拉長效應，作者首先利用直接數(shù)值模擬(Direct Numerical Simulation, DNS)得到的射流數(shù)據(jù)生成光場相機數(shù)字合成圖像，然后分別利用單光場相機和雙光場相機測量速度場，對比兩者的速度測量精度，尤其是對比兩者沿景深方向的速度測量誤差。

圖8 體校準及重構(gòu)算法流程圖

Fig.8Flowchatofthevolumetriccalibrationandreconstructionalgorithm

DNS數(shù)據(jù)模擬了一個圓形噴管的流動，噴管直徑D=20mm,雷諾數(shù)Re=2500，測試數(shù)據(jù)抽取自該噴管出口大約1D處的速度場，大小為0.66D×0.66D×0.66D[29]。粒子濃度取1.0ppm(Particle per Microlens),即在該區(qū)域隨機生成一些粒子的空間位置并由光線跟蹤生成第一幀圖像，這些粒子的濃度符合上述濃度值，然后根據(jù)DNS數(shù)據(jù)得到這些粒子所在位置點的速度矢量，給定時間間隔(2ms),由這些速度矢量計算出下一幀的粒子空間位置并生成合成圖像。單相機的體素數(shù)目為800×800×267voxels, 體素像素比為3∶3∶10，因此體素的空間分辨率為0.0165mm×0.0165mm×0.055mm，雙相機體素數(shù)目為800×800×800voxels,像素體素比為3∶3∶3，相應的空間分辨率為0.0165mm×0.0165mm×0.0165mm。為節(jié)約時間成本，所有圖像均采用NVIDIA GTX 1080Ti 顯卡CUDA C并行計算處理，得到粒子強度場后同樣利用并行程序進行基于FFT的兩重網(wǎng)格互相關算法計算出速度矢量，最終對計算得到的速度矢量場進行中值濾波和線性插值等后處理。

DNS原始速度場、單光場相機速度場及雙光場相機速度場如圖9所示。計算結(jié)果表明，無論是單光場相機還是雙光場相機，其計算結(jié)果和原始數(shù)據(jù)場都比較吻合，均能捕捉到流場中的渦系結(jié)構(gòu)。為定量對比分析兩者的速度測量誤差，圖10分別展示了兩者在x-y和x-z平面的速度測量絕對誤差分布散點圖以及兩者在x-y-z3個方向上的累計誤差分布函數(shù)(Cumulative Distribution Function, CDF)的對比(其中坐標系的定義為：x-y平面為垂直于某一個相機光軸的平面，z方向為光軸方向)。由結(jié)果可以清楚地看出，單相機和雙相機在x-y平面上的測量精度相近，但在單個相機坐標系的z方向上(即其光軸方向上)，雙相機的測量精度明顯高于單相機，并且達到了與x-y平面精度大小相近的水平。由累計誤差分布函數(shù)可以看出，雙相機沿光軸方向的測量誤差與單個相機相比明顯減小。

(a) 原始DNS流場

(b) 單光場相機速度場

(c) 雙光場相機速度場

Fig.9ComparisonofcalculationresultsofsyntheticimagesofDNSjet

(a) x-y平面的速度誤差對比

(b) x-z平面的速度誤差對比

(c) x-y-z方向的速度誤差累計分布函數(shù)對比

Fig.10Comparisonofmeasurementerrorsbetweensingleanddualcameras

3 渦環(huán)實驗對比分析

由于數(shù)字合成圖像沒有考慮到實驗條件下的如背景噪聲、光照不均勻等諸多限制因素，為進一步驗證雙光場相機技術的可行性，開展了水射流渦環(huán)PIV測量驗證實驗。驗證實驗的裝置系統(tǒng)如圖11所示。2臺光場相機布置時保持光軸互相垂直。為產(chǎn)生實驗所需的渦環(huán)，采用一臺步進電機來驅(qū)動一個可以作直線沖擊運動的活塞，活塞運動時擠壓圓柱形流道中的流體，從而向靜止的水缸中注入一股圓柱形的射流，由于活塞所在流道的壁面無滑移速度邊界條件，壁面有著很大的速度梯度，因此在流道的出口處形成了一個渦環(huán)[30]。

圖11 渦環(huán)實驗裝置示意圖

活塞的行程與直徑之比L/D為1.5，直徑D=20mm，對應的層流狀態(tài)下雷諾數(shù)Re=2000。2臺相機相互垂直90°放置，光圈數(shù)為2.8，放大系數(shù)為1，由此可以得到測量體積大約為36mm×24mm×36mm。測量區(qū)域在噴嘴出口前方2.5D處，示蹤粒子型號為Dantec PSP-50，其平均粒子直徑為50μm。

重構(gòu)的體素像素比在x,y,z3個方向上分別為3∶3∶3，速度采用兩重網(wǎng)格互相關，第一重網(wǎng)格的窗口大小為200×200×200voxels, 第二重網(wǎng)格的窗口大小為100×100×100voxels，互相關窗口重疊率為0.75，由此得到的速度矢量空間分辨率為1.65mm×1.65mm×1.65mm,速度場的矢量總數(shù)為73×53×73。

由測量實驗得到的速度場可以看出(見圖12)，2個光場相機各自對渦環(huán)進行測量時,在上述互相關參數(shù)下，在各自的光軸方向(由于2個相機光軸相互垂直，所以對第一個相機而言，其光軸方向為圖中坐標系的z軸方向，對第二個相機而言則為x軸方向)的測量精度有限，無法準確捕捉到該方向的渦環(huán)速度矢量，而雙相機系統(tǒng)則可以比較精準地測出渦環(huán)在3個方向的速度場。

(a) 第一個相機速度場

(b) 第二個相機速度場

(c) 雙相機速度場

測量速度的散度對比分析如圖13所示。由于水射流的速度場是不可壓流體的速度場，因此由不可壓流體的連續(xù)性方程可知速度的梯度應滿足：

(5)

式中，u,v,w分布表示沿x，y，z方向的速度。

為定量分析測量速度的梯度對式(5)的偏離程度，測量誤差由式(6)所示的散度絕對值來定義：

從圖13中的速度散度的概率密度分布函數(shù)(Probability Density Function，PDF)可以看出，雙相機的測量結(jié)果更加符合不可壓流體的連續(xù)性方程。

圖13 單光場相機和雙光場相機的速度散度概率密度函數(shù)對比

Fig.13Comparisonoftheprobabilitydensityfunctiononthedivergenceerrorofthesingle-anddual-cameraconfigurations

4 結(jié) 論

本文闡述了基于雙光場相機的高分辨率光場三維PIV技術，基于針對光場相機的校準算法和MART權(quán)重系數(shù)算法，發(fā)展了適用于雙光場相機乃至多臺光場相機的三維粒子強度場重構(gòu)技術，通過數(shù)字合成圖像仿真分析和渦環(huán)測量實驗對比分析發(fā)現(xiàn)，雙光場相機與單相機相比，在光軸方向上有更高的空間分辨率，而且精度達到了和垂直光軸方向的分辨率同等大小的水平；雙相機對水射流的測量結(jié)果更加符合不可壓流體的連續(xù)性方程，從而表明雙光場相機可以更加精準地測量景深方向的速度矢量，提高了測量的準確性。