蔡章博，張征宇，余皓，占書(shū)恒

1.西南科技大學(xué) 信息工程學(xué)院，綿陽(yáng) 621000

2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 高速空氣動(dòng)力研究所，綿陽(yáng) 621000

降低摩阻意味著飛行器的油耗下降、航程增加［1-3］?；跓晒庥湍さ娜帜ψ柚苯訙y(cè)量法，具有測(cè)量設(shè)備簡(jiǎn)單（僅需紫外激發(fā)光源、相機(jī)與鏡頭）、空間分辨率高、對(duì)模型物面無(wú)特殊要求等優(yōu)點(diǎn)，已成為研究的熱點(diǎn)和發(fā)展趨勢(shì)［4-6］。該法是基于油流法［5］，在油膜中加入熒光分子，通過(guò)熒光油膜控制方程，描述熒光油膜（受紫外光激發(fā)的發(fā)光成像）灰度與其運(yùn)動(dòng)速度的關(guān)系，推導(dǎo)出基于熒光油膜的摩阻測(cè)量方程［6］。因其測(cè)量方程與光流方程［7］具有相似的形式，因此，采用光流法解算風(fēng)洞試驗(yàn)中模型表面熒光油膜運(yùn)動(dòng)的時(shí)序圖像，即可獲得試驗(yàn)?zāi)Ｐ偷慕诿媪鲃?dòng)結(jié)構(gòu)，為掌握試驗(yàn)?zāi)Ｐ捅诿姘l(fā)生流動(dòng)分離的位置、分離方式與特點(diǎn)、漩渦形成機(jī)制等提供重要的研究數(shù)據(jù)。

自1981 年Horn和Schunck［7］提出了Horn-Schunck（HS）光流法后，Wildes等［8］引入流體力學(xué)原理率先將其用于云運(yùn)動(dòng)估計(jì)，2008 年Liu和Shen［9］建立了光流法和流體流動(dòng)之間的數(shù)學(xué)聯(lián)系。2015年，Zhong等［10］重建了機(jī)翼連接處的全局摩擦力拓?fù)浣Y(jié)構(gòu)；Liu等［11］研究表明，光流法較粒子圖像測(cè)速法可獲取更密集的流體運(yùn)動(dòng)速度場(chǎng)。2016、2017年，黃湛等［12］、王鵬等［13］先后采用變分迭代方法求解表面摩阻分布，驗(yàn)證了基于熒光油膜的全局表面摩阻測(cè)量技術(shù)的有效性。2019年，鄒易峰等［4］通過(guò)添加人工網(wǎng)格線利用離散匹配，實(shí)現(xiàn)了模型振動(dòng)與油膜運(yùn)動(dòng)的解耦，提高了測(cè)量精度。

但是上述光流法的數(shù)值計(jì)算中，都采用了類(lèi)雅克比的矩陣分割迭代形式，存在計(jì)算復(fù)雜度高、收斂速度慢、耗時(shí)長(zhǎng)等問(wèn)題［14］，故現(xiàn)有光流法尚難在生產(chǎn)型風(fēng)洞的熒光油膜試驗(yàn)現(xiàn)場(chǎng)實(shí)時(shí)測(cè)量并顯示模型壁面流動(dòng)過(guò)程，降低了工程應(yīng)用的價(jià)值。

因此，高效的光流法對(duì)于基于熒光油膜的全局摩阻實(shí)時(shí)測(cè)量至關(guān)重要。尤其是隨著高分辨率（可達(dá)2 500 萬(wàn)～6 500 萬(wàn)像素）的高速工業(yè)相機(jī)在熒光油膜風(fēng)洞試驗(yàn)中的應(yīng)用，使用現(xiàn)有光流法處理2 GB/s 的海量熒光油膜時(shí)序圖像，需要在試驗(yàn)結(jié)束后花1～2 h 才能獲得一次試驗(yàn)的精細(xì)測(cè)量結(jié)果。

為此，本文提出熒光油膜速度場(chǎng)的空間自適應(yīng)快速光流解法，一方面在數(shù)值計(jì)算方式上，通過(guò)構(gòu)造共軛迭代式，較現(xiàn)有光流法的雅可比矩陣分割迭代式有更高的計(jì)算效率；另一方面，借鑒金字塔思想，對(duì)熒光油膜圖像進(jìn)行空間的灰度梯度與速度場(chǎng)自適應(yīng)降/升采樣，有效降低計(jì)算的復(fù)雜度，進(jìn)一步提升熒光油膜速度場(chǎng)的解算速度。

1 光流法的類(lèi)雅可比迭代式

當(dāng)前風(fēng) 洞試驗(yàn) 中常采 用Liu等［15-17］基于Horn-Schunck 光流所改進(jìn)的光流法，該方法追求一個(gè)全局能量泛函最小化的解，公式為

式中：Ix、Iy、It分別為圖像灰度沿x、y 方向空間域及時(shí)域上的梯度信息；α 為全局平滑參數(shù)，該值越大代表追求更平滑的解，即顆粒度更精細(xì)的速度場(chǎng)；(u，v)為給定坐標(biāo)點(diǎn)(x，y)的速度向量。通過(guò)歐拉-拉格朗日方程最小化求解可得

式中：Δ 為拉普拉斯算子，通常以有限差分形式近似Δu(x，y)=n((x，y)-u(x，y))；()為給定點(diǎn)坐標(biāo)下n 鄰域內(nèi)(u，v)的加權(quán)值［7］，代入式（2）可得

后續(xù)為便于構(gòu)造共軛迭代式，此處定義一個(gè)系數(shù)a ≈α 且a ≤α。當(dāng)a=α?xí)r，經(jīng)由式（3）構(gòu)造的現(xiàn)有光流迭代式為

式（4）為一種類(lèi)雅可比的矩陣分割迭代形式，收斂誤差為

式中：eu、ev分別為(u，v)當(dāng)前步和上一步的殘差。

現(xiàn)有光流法在迭代過(guò)程中，會(huì)運(yùn)用到迭代點(diǎn)局部鄰域的信息，其中灰度梯度大的像素點(diǎn)在式（4）中的加權(quán)值高，即局部鄰域施加的影響大。較共軛迭代式，該迭代方式所需存儲(chǔ)量大、收斂速度慢得多［18］，其收斂速度與迭代矩陣的譜半徑相關(guān)，尤其是隨著高分辨率高速工業(yè)相機(jī)應(yīng)用，在高分辨率條件下追求更為平滑的速度場(chǎng)時(shí)，相鄰像素點(diǎn)圖像灰度梯度信息變化趨于平緩，迭代矩陣譜半徑增大，進(jìn)而導(dǎo)致收斂耗時(shí)增加。

2 空間自適應(yīng)快速光流解法

2.1 共軛迭代式推導(dǎo)

考慮到油膜速度場(chǎng)中速度梯度大的像素點(diǎn)其灰度梯度大，為此，約定(u，v)等于鄰域內(nèi)的速度最大值(um，vm)，有

鑒于共軛梯度法所需存儲(chǔ)量小，具有步收斂性、穩(wěn)定性高的優(yōu)點(diǎn)［18］，為此，借鑒局部光流法中“鄰域內(nèi)速度具有相關(guān)性”的思想［19］有

代入式（6），有

式中：λ 的大小反映了鄰域內(nèi)的平滑程度，為接近0 的正值，通常取0.001～0.010。故式（3）可化為

化為Ax=b 形式，即

顯然A 為一個(gè)對(duì)稱(chēng)正定矩陣，必然存在共軛向量ri滿(mǎn)足

設(shè)rk+1=rk-αkApk，其中，pk為殘差方向的單位向量；αk為步長(zhǎng)。Schmidt 正交化可得

為了求解αk，構(gòu)造如式（14）的二次型：

式中：x=(uk，vk)T。式（14）與式（10）同解，將f (xk+akpk)記作g(ak)，令g′(ak)=0 可得

至此，求解光流方程的共軛迭代式為

2.2 油膜圖像的空間自適應(yīng)采樣

為提高光流法在大位移下的精度，有學(xué)者提供了網(wǎng)格匹配［21］、金字塔采樣［22-24］的思路。本文借鑒金字塔思想，通過(guò)對(duì)熒光油膜圖像進(jìn)行空間的灰度梯度與速度場(chǎng)自適應(yīng)降/升采樣，以有效提高大運(yùn)動(dòng)估計(jì)的精度并降低計(jì)算的復(fù)雜度。為此，本文提出了基于八鄰域的空間尺度的灰度梯度與速度場(chǎng)自適應(yīng)降/升采樣方法。

如圖1 所示，沿(u，v)方向?qū)⒔o定圖像的灰度梯度(Iu，Iv)矩陣按照3×3 大小進(jìn)行分塊，邊界未能整除分塊的部分以0 填充對(duì)齊，對(duì)于給定的3×3 圖像塊，通過(guò)自適應(yīng)降采樣后得到1 個(gè)像素，其灰度梯度為如此，即可解算得到空間降采樣層的速度場(chǎng)(u′，v′)，再借助原八鄰域梯度信息對(duì)(u′，v′)升采樣恢復(fù)到原有分辨率下的(u，v)。

圖1 自適應(yīng)采樣示意圖Fig.1 Process of proposed adaptive sampling

2.2.1 基于八鄰域的灰度梯度自適應(yīng)降采樣

對(duì)于給定的熒光油膜圖像上3×3 圖像塊，為了保留八鄰域內(nèi)權(quán)重最大圖像灰度梯度信息，可將其灰度梯度分為兩類(lèi)，即

式中：Iumax、Iumin分別為原始圖 像(x，y)的n鄰域內(nèi)u 方向的梯度最大、最小值；v 方向的梯度計(jì)算同理。降采樣示意圖如圖2 所示，在該八鄰域內(nèi)梯度高于均值的點(diǎn)更多，滿(mǎn)足式（14）中取Iumax的條件，則取該區(qū)域內(nèi)的梯度最高值12 作為采樣值。

圖2 灰度梯度自適應(yīng)降采樣示例Fig.2 Example of proposed gray-gradient down-sampling

較傳統(tǒng)的雙線性插值的降采樣方法，本文方法得到灰度梯度矩陣會(huì)保留灰度梯度顯著區(qū)域的信息。

2.2.2 基于八鄰域的速度場(chǎng)自適應(yīng)升采樣

通過(guò)共軛光流計(jì)算，得到灰度梯度自適應(yīng)降采樣圖像的速度場(chǎng)。在后續(xù)計(jì)算中需要將速度場(chǎng)恢復(fù)到原圖像空間分辨率。鑒于傳統(tǒng)的雙線性插值方法難以在空間上升采樣后保留具有突變特征的速度場(chǎng)，為此，提出八鄰域的速度場(chǎng)自適應(yīng)升采樣法，通過(guò)原有梯度信息對(duì)應(yīng)恢復(fù)鄰域內(nèi)速度向量，計(jì)算公式為

圖3 速度場(chǎng)八鄰域恢復(fù)示例Fig.3 Example of proposed velocity up-sampling

考慮到油膜速度場(chǎng)中速度梯度大的像素點(diǎn)其灰度梯度大，按照局部光流法理論用u、v 代替后，既減少了速度矩陣的存儲(chǔ)量和硬件線程同步耗時(shí)，又能采用共軛迭代法提升解算效率，快速得到表達(dá)大流動(dòng)結(jié)構(gòu)的速度場(chǎng)初值。再將速度場(chǎng)初值按照其八鄰域內(nèi)點(diǎn)中灰度梯度分布插值，得到原圖空間分辨率速度場(chǎng)初值，再進(jìn)行全局速度場(chǎng)優(yōu)化，即可準(zhǔn)確捕捉到精細(xì)的流動(dòng)結(jié)構(gòu)，有效提高熒光油膜速度場(chǎng)的解算速度。

2.3 算法流程

本文算法全部流程如圖4 所示，兩幀時(shí)序熒光油膜圖像分別對(duì)應(yīng)圖像Pt、Pt+1。

圖4 本文方法流程Fig.4 Process of proposed method

3 仿真實(shí)驗(yàn)與分析

3.1 Ossen 渦對(duì)實(shí)驗(yàn)

參照文獻(xiàn)［4，6］設(shè)計(jì)仿真實(shí)驗(yàn)，利用圖5（a）、圖5（b）所示的熒光油膜模擬圖像表征給定的速度場(chǎng)為

圖5 熒光油膜仿真圖像及Ossen 渦對(duì)速度場(chǎng)Fig.5 Fluorescent oil film simulation images and Ossen vortex pair velocity field

式中：Γ 為渦核強(qiáng)度，本實(shí)驗(yàn)設(shè)定為5 000 pixel2/s；r0=30 pixel；u=10 pixel/s。在Ossen 渦對(duì)速 度場(chǎng)上再疊加一個(gè)勻速速度場(chǎng)，式（20）速度場(chǎng)疊加效果如圖5（c）所示，對(duì)圖5（a）利用雙線性插值，在演化步長(zhǎng)為0.000 1 s 的條件下演化1 000步，取間隔0.1 s 后演化至5（b）所示熒光油膜模擬圖。

為驗(yàn)證本文算法準(zhǔn)確性，對(duì)上述流場(chǎng)下的演變圖像進(jìn)行測(cè)試，限定迭代次數(shù)為300 進(jìn)行對(duì)比實(shí)驗(yàn)，為了扣除邊界條件對(duì)迭代計(jì)算的影響，忽略邊界像素點(diǎn)（x<20 pixel，x>200 pixel）的結(jié)果，得到的仿真實(shí)驗(yàn)結(jié)果如圖6 所示。如表1 所示，本文方法計(jì)算的誤差更小，相比于文獻(xiàn)［4］方法最大誤差減小了2.3%，平均誤差減小1.4%。從圖6 也可看出本文方法較現(xiàn)有方法與理論曲線更為貼合，在一些速度變化較為平緩的區(qū)域本文方法效果更好。

表1 迭代300 次沿渦對(duì)分布的速度相對(duì)誤差Table 1 Average errors distributing along vortex pair under 300 iterations

圖6 不同方法沿Ossen 渦對(duì)分布的測(cè)量速度Fig.6 Estimated velocities distributing along vortex cores under various methods

3.2 壁面射流實(shí)驗(yàn)

參照文獻(xiàn)［25］選取2 張壁面射流圖像模擬流體油膜運(yùn)動(dòng)進(jìn)行仿真實(shí)驗(yàn)，為方便觀察，論文圖像在實(shí)驗(yàn)圖像基礎(chǔ)上做了灰度拉伸（見(jiàn)圖7）。為驗(yàn)證比較收斂速度，統(tǒng)計(jì)文獻(xiàn)［25］方法收斂到給定誤差限的迭代次數(shù)。壁面射流實(shí)驗(yàn)共設(shè)定了3 個(gè)誤差 限：1.0×10-6、5.0×10-7、1.0×10-7，統(tǒng)計(jì)了不同方法收斂到誤差限所需迭代次數(shù)及所用時(shí)間。統(tǒng)計(jì)結(jié)果如表2 所示。

表2 壁面射流試驗(yàn)在不同誤差限下收斂的迭代次數(shù)Table 2 Iterations to convergence with various error limits in wall-jet experiment

圖7 壁面射流圖像Fig.7 Images of wall-jet

可以看出在誤差限1.0×10-6、5.0×10-7的標(biāo)準(zhǔn)下，以耗時(shí)為例，本文方法效率分別提升了40.72%、26.66%。在1.0×10-7的誤差限下效率提升了5.73%，這是由于后續(xù)全局優(yōu)化的解算性質(zhì)，隨著迭代次數(shù)提高，整體的收斂速度放緩。如圖8 所示，在相同的全局收斂誤差限下，本文方法解得的流場(chǎng)特性會(huì)更加明顯。

圖8 同誤差限下實(shí)驗(yàn)結(jié)果對(duì)比Fig.8 Comparison of results with same error limit

4 風(fēng)洞試驗(yàn)

4.1 某機(jī)翼表面熒光油膜風(fēng)洞試驗(yàn)

以某大展弦比靜彈性試驗(yàn)機(jī)翼為對(duì)象，在2.4 m 跨聲速式風(fēng)洞中開(kāi)展驗(yàn)證試驗(yàn)（馬赫數(shù)Ma=0.82），相機(jī)分辨率為2 352 pixel×1 728 pixel，曝光時(shí)間為16 ms，鏡頭焦距為35 mm，采集幀率60 Hz，最終解算油膜速度場(chǎng)分辨率為2 000 pixel×1 200 pixel，采用本文算法處理該試驗(yàn)采集的時(shí)序圖像，以便與文獻(xiàn)［4］的方法進(jìn)行對(duì)比。

模型迎角α=0°時(shí)，試驗(yàn)采集相鄰兩幀圖像如圖9（a）、圖9（b）所示，本文求解到的油膜運(yùn)動(dòng)速度場(chǎng)，通過(guò)視頻測(cè)量技術(shù)［26-27］可得到速度的標(biāo)準(zhǔn)單位（m/s），如圖9（c）所示。

圖9 機(jī)翼試驗(yàn)圖像及解算結(jié)果Fig.9 Images and estimated result of wing test

由表3 可知，在誤差限為1.0×10-6、5.0×10-7、1.0×10-7的標(biāo)準(zhǔn)下，以耗時(shí)為例，本文方法效率分別提升了103.25%、53.66%、26.17%。

表3 機(jī)翼表面熒光油膜試驗(yàn)在不同誤差限下收斂的迭代次數(shù)（α=0°）Table 3 Iterations to convergence with various error limits in fluorescent oil film test on wing surface（α=0°）

模型迎角α=10°時(shí)，取試驗(yàn)采集相鄰兩幀圖像（見(jiàn)圖10（a）、圖10（b）），誤差限為1.0×10-7標(biāo)準(zhǔn)下求得的摩擦應(yīng)力線圖譜如10（c）所示，以圖中矩形區(qū)域?yàn)槔嘏c文獻(xiàn)［4］對(duì)應(yīng)結(jié)果進(jìn)行比較，其油流分離現(xiàn)象一致，本文方法解算效率提升26.17%。相比于傳統(tǒng)油流試驗(yàn)圖像，本文方法得到了更豐富、清晰的流動(dòng)細(xì)節(jié)。該試驗(yàn)表明本文方法適用于曲面試驗(yàn)?zāi)Ｐ?。與文獻(xiàn)［4］中算法效率比較如表4 所示，在誤差限為1.0×10-6、5.0×10-7、1.0×10-7的標(biāo)準(zhǔn)下，以耗時(shí)為例，本文算法效率分別提升了9.76%、69.47%、19.89%。

表4 機(jī)翼表面熒光油膜試驗(yàn)在不同誤差限下收斂的迭代次數(shù)（α=10°）Table 4 Iterations to convergence with various error limits in fluorescent oil film test on wing surface（α=10°）

圖10 機(jī)翼試驗(yàn)圖像及解算結(jié)果Fig.10 Images and estimated result of wing experiment

4.2 平板上圓柱繞流油流試驗(yàn)

在中國(guó)空氣動(dòng)力研究與發(fā)展中心（CARDC）的0.6 m 連續(xù)式風(fēng)洞中開(kāi)展了“平板上的圓柱繞流油流”試驗(yàn)（Ma=0.6）。圖11（a）中所示為試驗(yàn)所用平板，紅色區(qū)域?yàn)闊晒庥湍?，圖11（b）、圖11（c）為相鄰兩幀平板上圓柱繞流引起的油膜運(yùn)動(dòng)圖像，相機(jī)分辨率為5 120 pixel × 5 120 pixel，曝光時(shí)間為12.5 ms，鏡頭焦距為30 mm，采集頻率80 Hz，解算油膜速度場(chǎng)分辨率為1 200 pixel×780 pixel。

圖11 平板試驗(yàn)圖像及解算結(jié)果Fig.11 Images and estimated result of flat plate experiment

如圖11（d）所示，當(dāng)平板上的來(lái)流遇到圓柱障礙時(shí)，平板上圓柱附近的壁面邊界層速度開(kāi)始下降、渦量開(kāi)始增加，隨著渦量的增加，邊界層發(fā)生流動(dòng)分離、不斷卷起形成漩渦并試圖脫離壁面時(shí)與來(lái)流相互作用，形成如圖11（d）、圖11（e）所示分離線，當(dāng)漩渦前進(jìn)一段距離再附時(shí)，將出現(xiàn)如圖11（d）、圖11（e）所示繞流附著線，其中圖11（e）所示分離點(diǎn)S1、S2，附著點(diǎn)A1 對(duì)應(yīng)圖11（d）所標(biāo)記點(diǎn)。以來(lái)流方向視角觀察氣流在S1 處發(fā)生分離時(shí)，邊界層氣流如圖11（e）所示由S1 流動(dòng)至附著點(diǎn)A1，經(jīng)過(guò)A1 的壁面氣流按圖示方式向前流動(dòng)并經(jīng)過(guò)S2。試驗(yàn)結(jié)果與文獻(xiàn)［16］試驗(yàn)測(cè)得流動(dòng)現(xiàn)象相似，表明本文算法不僅能正確解算風(fēng)洞試驗(yàn)中油流路徑，而且提供了清晰的摩擦力線圖譜。

5 結(jié)論

推導(dǎo)了構(gòu)造共軛迭代式，創(chuàng)建迭代效率更高的共軛光流法，將解得速度場(chǎng)作為初值有效提高后續(xù)全局優(yōu)化效率。對(duì)熒光油膜圖像進(jìn)行空間尺度的灰度梯度與速度場(chǎng)自適應(yīng)降/升采樣，進(jìn)一步提升高了熒光油膜速度場(chǎng)解算速度。

1）Oseen 渦對(duì)的熒光油膜路徑速度場(chǎng)仿真試驗(yàn)結(jié)果表明：本文方法在計(jì)算結(jié)果在相同迭代次數(shù)下相比現(xiàn)有方法收斂效果更好，限定迭代300 次的情況下，比現(xiàn)有方法最大誤差減小了2.3%，平均誤差減小1.4%。壁面射流仿真試驗(yàn)結(jié)果進(jìn)一步表明，在同誤差限下本文方法解算出來(lái)的流場(chǎng)特征更為清晰。

2）風(fēng)洞試驗(yàn)中，本文方法不僅測(cè)得了正確流動(dòng)現(xiàn)象，并且給出了清晰的摩擦力線圖譜，在誤差限較高標(biāo)準(zhǔn)下相較于文獻(xiàn)［4］方法，計(jì)算效率提升了26.17%。試驗(yàn)結(jié)果證明本文方法相較于現(xiàn)有方法優(yōu)勢(shì)明顯，工程應(yīng)用價(jià)值大。