基于三目視覺的圓柱尾跡層析PIV測量

段雙成，楊苗，封明軍，王芳婷，周騖，蔡小舒

（1 上海理工大學(xué)能源與動力工程學(xué)院，上海 200093；2 上海市動力工程多相流動與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

粒子圖像測速（PIV）是一種現(xiàn)代流場測量技術(shù)，不同于皮托管測速、熱線測速等傳統(tǒng)方法，PIV技術(shù)既擺脫了單點(diǎn)測速的限制，又實(shí)現(xiàn)了快速、準(zhǔn)確、無干擾地對流場截面結(jié)構(gòu)的測量，在化工、機(jī)械和能源等領(lǐng)域中都存在著廣泛的應(yīng)用，如天然氣流量特性分析[1]、攪拌槽混合槽內(nèi)的流動結(jié)構(gòu)測量[2]、噴砂桶內(nèi)部磨損模型[3]和攪拌釜內(nèi)槳型設(shè)計[4]。然而二維PIV技術(shù)在研究復(fù)雜流動中的渦流及湍流方面具有一定的局限性，如換熱器內(nèi)翅片間流動、燃?xì)廨啓C(jī)葉柵間流動等，其流道空間小，流動復(fù)雜，很難通過某一平面的二維結(jié)構(gòu)反映出來。層析粒子圖像技術(shù)（Tomographic-PIV）是在PIV技術(shù)上發(fā)展起來的一種三維流場測量技術(shù)，它采用一定厚度的激光體光源代替激光片光源，真正實(shí)現(xiàn)了3D-3C，即三維空間三個速度分量的速度場測量。

層析PIV是一種光學(xué)非定常流場測量技術(shù)，由于其具有大視場、高精度的三維測量能力，在提出之后便受到了國內(nèi)外的廣泛關(guān)注。Scarano[5]在2013年對層析PIV的研究進(jìn)展進(jìn)行了很詳細(xì)的綜述。層析PIV 的概念最早由2004 年由韓國的Kim 等[6]提出，并使用該技術(shù)對鼻腔內(nèi)的幾個平面內(nèi)鼻竇氣流組織進(jìn)行了測量，雖然得到的結(jié)果是三維速度分布，但實(shí)際上與目前的層析PIV 并不一致；2005年，荷蘭的Esinga 等[7]在兩種圓形氣缸中，應(yīng)用層析PIV 技術(shù)對其中的湍流流動進(jìn)行了測量；在第6屆國際PIV 研討會上，Esinga 等[8]正式提出了層析PIV技術(shù)，對其工作原理和成像算法的細(xì)節(jié)進(jìn)行了詳細(xì)的闡述，從迭代次數(shù)、相機(jī)數(shù)量以及示蹤粒子密度等方面對重構(gòu)精度進(jìn)行了分析，并對采用了四相機(jī)系統(tǒng)對環(huán)狀渦流模擬場進(jìn)行了測量，測量誤差僅為0.1 像素，表明該技術(shù)完全能夠滿足精度要求，并能夠提供遠(yuǎn)遠(yuǎn)大于二維測量的信息量。

Esinga 等[9]在其應(yīng)用于湍流邊界層的實(shí)驗(yàn)測量方面做了一系列的研究，并對偽影粒子進(jìn)行了研究；英國劍橋大學(xué)在早期提出了一種快速三維重建的方法[10]，這對校正圖像在空間中的映射關(guān)系至關(guān)重要，并對多種層析重建方法做出了貢獻(xiàn)。Worth等[11]還通過數(shù)值計算評估了各向同性湍流中層析PIV 的精度；Hain 等[12]對通過圓柱體頂部的流動拓?fù)浼坝邢迗A柱頂部邊界層和尾跡結(jié)構(gòu)之間的相互作用進(jìn)行了研究；Thomas等[13]對噴氣機(jī)的橫流射流進(jìn)行了測量研究；Ortizs-Due?as 等[14]通過層析PIV 研究了液滴在液-液平面上的聚結(jié)現(xiàn)象以及單個液滴之間的聚合過程，首次實(shí)現(xiàn)了多相流的測量；Silva等[15]對一個典型的湍流邊界層層析PIV實(shí)驗(yàn)中的幾個關(guān)鍵實(shí)驗(yàn)參數(shù)進(jìn)行了分析，結(jié)果表明相機(jī)的位置、示蹤粒子的數(shù)量及體素的離散程度等因素均對實(shí)驗(yàn)的準(zhǔn)確性有顯著影響；Violato等[16]采用時間解析層析PIV算法對從圓形出口和V形出口處的射流轉(zhuǎn)捩過程進(jìn)行了研究，得到了高速和低速流動區(qū)與激波邊界層相互作用的拓?fù)浣Y(jié)構(gòu)。

上述研究表明了層析PIV的發(fā)展過程由最初的致力于層析重構(gòu)算法，到逐步集中在對湍流現(xiàn)象上的研究，也正是由于其能夠測量湍流中相干結(jié)構(gòu)的三維瞬時狀態(tài)，使得其應(yīng)用越來越廣。

國內(nèi)對于層析PIV 技術(shù)的研究起步相對較晚，研究內(nèi)容主要集中在低速水流、風(fēng)洞中。天津大學(xué)姜楠等[17]應(yīng)用改進(jìn)的層析TRPIV算法，對水槽中的平板湍流邊界層的三維速度場進(jìn)行了測量；北京航空航天大學(xué)高琪等[18]通過自主研發(fā)的層析PIV 技術(shù)，實(shí)現(xiàn)了對低速合成射流的三維測量，觀察到了三維渦環(huán)結(jié)構(gòu)的時空演化過程；中國航空工業(yè)動力研究院許相輝等[19]對圓柱三維尾流場進(jìn)行了測量，成功觀測到了圓柱體后方產(chǎn)生的三維卡門渦流場；西南科技大學(xué)張小航等[20]采用了數(shù)值模擬的方法，通過層析PIV技術(shù)的原理及方法，以真實(shí)風(fēng)洞中的測量參數(shù)為基礎(chǔ)，對由數(shù)值計算得到的模擬流場進(jìn)行了重建研究，驗(yàn)證了方法的準(zhǔn)確性。這些成果極大促進(jìn)了國內(nèi)學(xué)者對于層析PIV的研究興趣，促進(jìn)了技術(shù)的發(fā)展。同時也顯示出了層析PIV對三維渦環(huán)結(jié)構(gòu)具有很好的捕捉能力，能確保實(shí)現(xiàn)復(fù)雜流動中渦環(huán)時空演變的定量實(shí)驗(yàn)研究。

因此本文基于層析PIV技術(shù)，結(jié)合3臺CCD相機(jī)等硬件設(shè)施搭建了一套層析PIV測量系統(tǒng)，采用三維相機(jī)標(biāo)定、重構(gòu)算法，對水槽內(nèi)圓柱繞流尾跡流場進(jìn)行測量研究，完成了低速條件下的實(shí)驗(yàn)驗(yàn)證，為復(fù)雜的狹小空間流動測量奠定了基礎(chǔ)。

1 層析PIV基本原理

層析PIV 的原理與計算機(jī)斷層攝影（CT）類似，后者通過X射線面對物體進(jìn)行掃描，通過計算X射線的衰減量來重建物體的某橫斷面結(jié)構(gòu)；但流場不同于靜止的物體，需要通過高速相機(jī)或具有PIV模式的相機(jī)“凍結(jié)”某一時刻下的流場，并用激光代替X射線，采用多臺相機(jī)在多角度拍攝以代替X射線面的掃描，而相機(jī)采集到的圖像像素值便是多粒子散射光強(qiáng)疊加形成的，通過反演計算便能得到粒子的空間位置，完成三維粒子光強(qiáng)場的重建。其測量原理示意圖如圖1所示，通過激光器產(chǎn)生激光束并通過透鏡系統(tǒng)擴(kuò)束產(chǎn)生具有一定厚度的體光源，照亮示蹤粒子；布置至少3臺CCD相機(jī)從不同角度對示蹤粒子的散射光信息進(jìn)行采集，并通過信號發(fā)生器使CCD相機(jī)在極小的時間間隔下拍攝兩次；然后采用層析重構(gòu)算法對三維粒子場信息進(jìn)行重建，再利用三維互相關(guān)技術(shù)得到粒子在該時間間隔下的位移，從而實(shí)現(xiàn)對三維速度場的重建。

圖1 三目PIV測量原理

1.1 三維標(biāo)定

相機(jī)的標(biāo)定實(shí)際上就是4 個坐標(biāo)系的轉(zhuǎn)換過程，其中包括世界坐標(biāo)系、相機(jī)坐標(biāo)系、圖像坐標(biāo)系和圖像像素坐標(biāo)系，其目的就是確定圖像像素坐標(biāo)系與世界坐標(biāo)系的對應(yīng)關(guān)系，以便在對采集到的圖像進(jìn)行處理時對物體實(shí)際坐標(biāo)進(jìn)行推導(dǎo)。圖2為4種坐標(biāo)系間的相互關(guān)系，其中O-XYZ為世界坐標(biāo)系，o"-uv為圖像像素坐標(biāo)系，o'-u'v'為圖像坐標(biāo)系，Oc-XcYcZc為相機(jī)坐標(biāo)系。由世界坐標(biāo)系下某點(diǎn)P坐標(biāo)（x，y，z）到圖像像素坐標(biāo)系下坐標(biāo)（u,v）的轉(zhuǎn)換關(guān)系如式(1)所示。

圖2 坐標(biāo)系概況

式中，R表示世界坐標(biāo)系到相機(jī)坐標(biāo)系的旋轉(zhuǎn)矩陣；t表示平移向量；f為相機(jī)焦距；dx和dy為圖像像素的物理尺寸（u0，v0）為圖像中心的像素坐標(biāo)。相機(jī)坐標(biāo)系Oc-XcYcZc可以看成由世界坐標(biāo)系O-XYZ經(jīng)過旋轉(zhuǎn)與平移得到。

相機(jī)坐標(biāo)系到圖像物理坐標(biāo)系的轉(zhuǎn)換則通過透鏡成像的相似三角形比例關(guān)系進(jìn)行計算；圖像物理坐標(biāo)系則以圖像的中心為原點(diǎn)。

相機(jī)標(biāo)定一般采用已知尺寸的圓點(diǎn)標(biāo)定板或棋盤格，在不同位置和角度下進(jìn)行成像，處理獲得上述參數(shù)，多相機(jī)標(biāo)定則是對多臺相機(jī)同時獲得相應(yīng)的標(biāo)定參數(shù)。

1.2 三維重構(gòu)

單臺CCD 相機(jī)只能得到二維圖像強(qiáng)度分布，三維空間的光強(qiáng)分布需要通過重構(gòu)獲得。首先需要將被測流場的三維區(qū)域進(jìn)行離散化，得到若干數(shù)量的體素；根據(jù)粒子成像理論可以得到，單個像素的強(qiáng)度值是基于空間視線內(nèi)多個粒子散射強(qiáng)度的線性積分，所以可以得到式(2)。

式中，I(xi,yi)為圖像像素強(qiáng)度；E(Xj,Yj,Zj)為體素散射光強(qiáng)；Ni為光線所截取的或在光線鄰域內(nèi)的體素數(shù)量；i表示相機(jī)的第i個像素；ω(i,j)為權(quán)重系數(shù)矩陣，代表了體素對像素強(qiáng)度的貢獻(xiàn)率。

不難看出，若要得到三維光強(qiáng)分布，首先要得到權(quán)重系數(shù)ω(i,j)。圖3為重構(gòu)的數(shù)學(xué)模型，描述了在一個層析平面上的重構(gòu)過程。權(quán)重系數(shù)矩陣ω(i,j)由體素相對于像素的尺寸以及體素中心與光線的距離決定，且有0≤ω(i,j)≤1。

圖3 層析重構(gòu)算法

在得到權(quán)重系數(shù)后，由于矩陣階數(shù)很大，難以用常規(guī)的矩陣方法求解，故通常采用迭代法求解。目前經(jīng)常采用的算法為ART 算法和MART 算法，這兩種算法均通過選擇一個合適的初始值E(X,Y,Z)0開始進(jìn)行E(X,Y,Z)的迭代，且由于其采用了不同的迭代方式，導(dǎo)致了兩種算法的收斂速度不同，同時求解精度也有一定程度上的差異。

本文采用了基于乘法視距的MLOS 算法[21]，該算法與MART 算法類似，但不需要得到加權(quán)矩陣，而是通過視距估計，確定了粒子在整個體積中可能存在的位置，能夠極大減少運(yùn)算量，提高程序運(yùn)行速度，并在一定程度上降低重建噪聲。

MLOS算法是通過每個體素和像素強(qiáng)度的映射關(guān)系找出體素空間中的非零體素的方法，針對空間的任一體素，可以通過每臺相機(jī)的視線找到其在各個相機(jī)上的投影；若要使該點(diǎn)的體素強(qiáng)度大于0，則對應(yīng)位置的像素強(qiáng)度也一定大于0。同樣地，像素強(qiáng)度的乘積也一定大于0，這就使得可以通過相應(yīng)像素的乘積來判斷是否存在可以忽略的體素。如圖3所示，當(dāng)相應(yīng)像素沿視線相乘時，只有每幅圖像中非零像素相交的視線所對應(yīng)的體素為非零體素。該算法通過對每臺相機(jī)的標(biāo)定結(jié)果的分析就能夠直接找到對應(yīng)像素點(diǎn)，而不需要對體積中的每個點(diǎn)進(jìn)行極為耗時的權(quán)重矩陣計算和儲存。

但MLOS 算法與MART 算法一樣，均無法去除偽影粒子的存在。偽影粒子是由于該點(diǎn)體素強(qiáng)度能夠滿足在不同方向投影而產(chǎn)生的，這就意味著重建圖像雖然在數(shù)學(xué)上與原始圖像相同，但實(shí)際上該點(diǎn)并不存在相應(yīng)粒子。但在經(jīng)過多次迭代修正后，偽影粒子的強(qiáng)度通常低于真實(shí)粒子。故雖然存在偽影粒子，但其對后續(xù)的三維互相關(guān)計算的影響很小。

2 實(shí)驗(yàn)平臺與系統(tǒng)

2.1 實(shí)驗(yàn)管段

針對某矩形截面的水流直管段，管道截面尺寸60mm×60mm，管長1500mm，通過某35W水泵提供動力，并在管道前端安裝兩個穩(wěn)流器穩(wěn)定水流，管道中部安裝直徑5mm的圓柱體作為擾流裝置，將整個管道分為兩個層流段、一個繞流段和一個尾流段。本實(shí)驗(yàn)主要在圖4中虛線框所示的測量區(qū)域尾流段中進(jìn)行。圖5為本實(shí)驗(yàn)的實(shí)驗(yàn)裝置的結(jié)構(gòu)示意。

圖4 槽道流實(shí)驗(yàn)平臺管道結(jié)構(gòu)

圖5 三目層析PIV測量裝置

2.2 實(shí)驗(yàn)設(shè)備

本實(shí)驗(yàn)采用了3 臺帶有PIV 功能的IMPERX 相機(jī)，型號為B1411，其中一臺為彩色相機(jī)，兩臺為黑白相機(jī)，相機(jī)分辨率像素均為1360×1364；鏡頭采用創(chuàng)威工業(yè)鏡頭，型號為CW-VM0418-3MP-1/1.8，焦距范圍4～18mm，通光口徑為F1.6，可以通過調(diào)節(jié)光圈對景深進(jìn)行控制。

通常PIV實(shí)驗(yàn)均采用Nd:YAG雙脈沖式激光器作為照明光源，因其有脈沖能量高、脈沖間隔短的優(yōu)點(diǎn)，適用于高速流動下極小顆粒的拍攝；但該類激光器價格昂貴且調(diào)試過程復(fù)雜，且本實(shí)驗(yàn)拍攝對象為低速流體，對脈沖間隔要求不高，故定制了波長450nm、功率4W 的連續(xù)激光器，經(jīng)測試其光強(qiáng)完全可以滿足本實(shí)驗(yàn)需求，但若使用連續(xù)照明會導(dǎo)致PIV模式下第二幀圖像過曝，故采用了信號發(fā)生器控制激光器，實(shí)現(xiàn)脈寬1ms的脈沖式照明。

通常示蹤粒子的基本要求為跟隨性好、可見度高，使得圖像有盡可能高的信噪比。本實(shí)驗(yàn)的實(shí)驗(yàn)條件為水，故選擇了15μm的聚苯乙烯顆粒。

2.3 實(shí)驗(yàn)參數(shù)設(shè)置

實(shí)驗(yàn)前通過流量計測得管道水流平均速度為400L/h，可得流速約為0.04m/s，相機(jī)放大倍率約為0.1，PIV實(shí)驗(yàn)中一般粒子運(yùn)動距離為5～10個像素，同時曝光時間內(nèi)粒子運(yùn)動距離為1～3個像素，故設(shè)置兩幀間隔為5ms，曝光時間1ms。

3 結(jié)果與分析

3.1 三維標(biāo)定結(jié)果

三臺相機(jī)的標(biāo)定結(jié)果如圖6所示，由上至下分別對應(yīng)為左、中、右三臺相機(jī)的標(biāo)定結(jié)果，從左到右表示世界坐標(biāo)系下被離散成體素的空間區(qū)域后在像素坐標(biāo)系下x、y、z三個方向的投影。由于本實(shí)驗(yàn)中相機(jī)是在同一高度上成一定角度布置，并使標(biāo)定板平面與中間一臺相機(jī)的成像平面平行，左右兩臺相機(jī)對稱放置于z軸兩側(cè)，與中間相機(jī)夾角30°，且嚴(yán)格遵循沙姆定律[22]，因此中間相機(jī)的三個方向投影應(yīng)大致一致，兩側(cè)兩臺相機(jī)的zworld軸投影結(jié)果應(yīng)具有左右對稱性，可以看出標(biāo)定結(jié)果與預(yù)期基本相符。

圖6 左、中、右三臺相機(jī)的標(biāo)定結(jié)果

3.2 三維重建結(jié)果及分析

3.2.1 重建結(jié)果

實(shí)驗(yàn)拍攝區(qū)域?yàn)閳A柱后方約5mm 位置，重建區(qū)域被離散成600×600×100個體素，重構(gòu)區(qū)域?qū)嶋H大小為30mm×30mm×5mm，空間分辨率約為20voxels/mm，矢量密度為1.6points/mm；由于粒子運(yùn)動距離約為10個像素，每個體素含1個像素，故選擇163voxels為查詢窗口，重疊步長為50%，每臺相機(jī)圖像重建平均時間約為9.59s。圖7 為某時刻下重建的粒子場z方向不同位置的兩張切片圖，圖片中的明顯亮點(diǎn)即為重建得到的粒子，其中部分采用紅圈圈出作為示意。

圖7 三維粒子場重建切片圖

3.2.2 重投影誤差分析

由于實(shí)際實(shí)驗(yàn)中無法得到真實(shí)粒子場的切片，因此需要引入重投影誤差對重建精度進(jìn)行分析，即通過重建后的粒子場投影與初始粒子場圖片進(jìn)行比較，重投影誤差的計算如式(7)所示。

式中，Er為重投影誤差；(u,v)為圖像的像素坐標(biāo)；I0代表已知的初始粒子場圖片；I代表重建后投影得到的粒子場圖片，由重建后的粒子場通過式(1)計算得到。

隨機(jī)選取中間視角相機(jī)兩個不同時刻的三維重建結(jié)果，圖8為相應(yīng)的重投影相對誤差分布，平均重投影相對誤差不超過0.5%，具有相對較高的重建精度。

圖8 典型的重投影相對誤差分布

3.3 互相關(guān)計算結(jié)果

計算結(jié)果如圖9、圖10 所示，分別展示了圓柱后方某法向-流向切平面的渦量圖和流場圖。所測量位置為圓柱后方約5mm 位置處的下半部分區(qū)域，區(qū)域位置如圖4 虛線框所示，包含了圓柱繞流形成的卡門渦街以及下方層流部分。圖9 中可看到不同時刻旋渦的演化狀態(tài)，交替出現(xiàn)的紅色與藍(lán)色區(qū)域顯示了圓柱尾流上下方形成的具有一定規(guī)則的、交叉排列的渦列，并沿下游方向有擴(kuò)散的趨勢，這雨流體力學(xué)中圓柱尾跡卡門渦街的理論相一致。其中藍(lán)色渦量為順時針方向旋轉(zhuǎn)的負(fù)渦量。

圖10 互相關(guān)計算結(jié)果

三維速度場的重建結(jié)果如圖10(a)所示，整個重建空間內(nèi)包含72×72×9個矢量，所含矢量信息較多，故同時顯示相應(yīng)流線圖[圖10(b)]和某截面渦量圖[圖10(c)]?？梢钥闯隽鲌鲎笊戏叫纬闪艘粋€明顯的渦環(huán)結(jié)構(gòu)，同時其右側(cè)存在一旋轉(zhuǎn)方向相反、交替排列的渦，同時具有向下游方向擴(kuò)散的趨勢，與理論結(jié)果基本相符。圖10(d)為圖10(c)中紅框位置處形成渦環(huán)的局部放大圖。

4 結(jié)論與展望

（1）本文基于三目視覺搭建了一套層析PIV測量系統(tǒng)，針對矩形水流管道中的圓柱尾跡流動進(jìn)行了采集實(shí)驗(yàn)，進(jìn)行了三維標(biāo)定以及重建，最終的測量結(jié)果顯示，該層析PIV測量系統(tǒng)得到的三維速度場能夠合理地反映圓柱尾跡出口處的流動結(jié)構(gòu)，成功觀測到了卡門渦街結(jié)構(gòu)。

（2）本實(shí)驗(yàn)的實(shí)際測量區(qū)域?yàn)?0mm×30mm×5mm，故體素大小設(shè)置為600×600×100，空間分辨率達(dá)到20voxels/mm，互相關(guān)計算結(jié)果包含72×72×9 個矢量信息，相比于常規(guī)PIV 獲得了更加豐富的三維流場信息。

（3）MLOS 算法能夠有效降低重構(gòu)所耗時間，每臺相機(jī)的平均重建時間約為9.59s；同時能夠保持較高的重建精度，重投影誤差不超過0.5%；雖然MLOS算法無法消除偽影粒子的存在，但通過數(shù)次迭代計算大大降低了偽影粒子對最終互相關(guān)計算的影響。

（4）通過在低速條件下對復(fù)雜流場結(jié)構(gòu)的測量，驗(yàn)證了本文采用的三目層析PIV系統(tǒng)及布置方式的可行性為后續(xù)狹小空間區(qū)域流場測量研究奠定了基礎(chǔ)。