袁 媛 孫增玉 王 杏 高 越 劉 柯

(北京航天計(jì)量測試技術(shù)研究所，北京 100076)

1 引 言

在航空航天飛行器研制領(lǐng)域，通過風(fēng)洞模型試驗(yàn)研究飛行器性能，降低飛行器研制風(fēng)險(xiǎn)和成本。在雷諾數(shù)模擬能力相同的情況下，相比于常規(guī)風(fēng)洞，低溫風(fēng)洞顯著降低了風(fēng)洞尺寸及對洞體結(jié)構(gòu)的強(qiáng)度要求，風(fēng)洞造價(jià)也隨之降低，所需動(dòng)壓、驅(qū)動(dòng)功率、峰值功率和總能消耗降低。由于低溫風(fēng)洞運(yùn)行環(huán)境復(fù)雜性，相關(guān)實(shí)驗(yàn)研究仍處于探索階段，風(fēng)洞天平是測量風(fēng)洞中氣流作用在模型上的氣動(dòng)力和力矩的測量設(shè)備，可以實(shí)現(xiàn)氣動(dòng)力和力矩沿三個(gè)坐標(biāo)系分解并精確測量。但風(fēng)洞天平的測量精度受低溫風(fēng)洞內(nèi)運(yùn)行溫度及運(yùn)行環(huán)境的影響較大，諸如低溫下傳感器材料的特性變化、低溫下應(yīng)變片零點(diǎn)及工作靈敏度的溫度補(bǔ)償難以得到有效控制等。因此急需針對超低溫高速風(fēng)洞的特殊試驗(yàn)環(huán)境，研究飛行器應(yīng)力的高精度測量技術(shù)，以非接觸測量的方式解決超低溫風(fēng)洞下模型狀態(tài)測量，為飛行器試驗(yàn)提供必要檢測手段[1,2]。

與二維數(shù)字圖像相關(guān)不同，三維數(shù)字圖像相關(guān)法可以實(shí)現(xiàn)被測模型的全場高效三維形貌和變形測量，且不需要二維數(shù)字圖像相關(guān)法中對于測量相機(jī)的光軸與被測模型表面相垂直和被測模型主要為面內(nèi)形變等條件約束，其主要優(yōu)點(diǎn)有非接觸高效光學(xué)精密測量、大量程動(dòng)態(tài)實(shí)時(shí)測量、高精度多參數(shù)同步測量，整個(gè)測量光路簡單且適合于特殊運(yùn)行環(huán)境下的應(yīng)力測量。本文直接利用被測物體表面人工形成的隨機(jī)散斑圖案在變形前后參考和目標(biāo)圖像子區(qū)的灰度變化來獲取模型表面的三維位移場和應(yīng)變場，進(jìn)而獲得模型在低溫環(huán)境下的應(yīng)力場。

2 三維數(shù)字圖像相關(guān)方法

三維數(shù)字圖像相關(guān)法是采用兩臺(tái)高速攝像機(jī)同步采集被測模型表面圖像，結(jié)合雙目立體視覺進(jìn)行變形測量的一種光學(xué)非接觸高精度測量方法[3,4]。其基本思想是通過拍攝和處理被測模型圖像子區(qū)變形前后二維圖像獲得被測模型表面的三維位移和應(yīng)變分布。三維數(shù)字圖像相關(guān)法通過匹配被測模型變形前后兩幅散斑圖像中的不同圖像子區(qū)來計(jì)算三維位移場，通過數(shù)值微分計(jì)算獲得應(yīng)變場[5～7]。

2.1 雙目視覺成像模型及相機(jī)內(nèi)外參數(shù)標(biāo)定

雙目視覺技術(shù)的主要原理是基于視差來實(shí)現(xiàn)，利用兩個(gè)空間位置姿態(tài)關(guān)系已知的相機(jī)在不同方位下同時(shí)采集被測特征的圖像，通過圖像處理及同名點(diǎn)匹配等技術(shù)獲取被測特征對應(yīng)的同名像點(diǎn)對，利用相機(jī)成像模型建立成像光線方程，構(gòu)建光線三角交會(huì)約束，組建雙目視覺測量數(shù)學(xué)模型，解算被測特征空間三維坐標(biāo)[8]。雙目立體視覺測量原理示意圖如圖1所示。

圖1 雙目立體視覺測量原理示意圖Fig.1 Schematic diagram of binocular stereo vision measurement principle

(1)

式中：R——o-xyz坐標(biāo)系與or-xryrzr坐標(biāo)系間的旋轉(zhuǎn)矩陣；T——兩坐標(biāo)系間的平移矩陣。

空間三維坐標(biāo)可以表示為

(2)

已知焦距fl、fr和空間點(diǎn)在左右相機(jī)中的圖像坐標(biāo)及旋轉(zhuǎn)矩陣R和平移矩陣T，可解算空間任一點(diǎn)的坐標(biāo)。

雙目相機(jī)內(nèi)外參數(shù)標(biāo)定是基于三維數(shù)字圖像相關(guān)法解算應(yīng)變的前提，標(biāo)定結(jié)果的好壞影響到最終的測量結(jié)果。相機(jī)標(biāo)定就是為了正確建立空間中某點(diǎn)和它在圖像平面上對應(yīng)像點(diǎn)之間的關(guān)系。相機(jī)標(biāo)定包括兩部分，分別為內(nèi)參數(shù)和外參數(shù)標(biāo)定，內(nèi)參數(shù)即相機(jī)和鏡頭的畸變相關(guān)參數(shù)，主要有主點(diǎn)位置、焦距、等效像素、畸變因子等；外參數(shù)即兩相機(jī)坐標(biāo)系間的空間轉(zhuǎn)換關(guān)系，包括位置關(guān)系和角度關(guān)系。根據(jù)透視投影關(guān)系有

s=A[RT]

(3)

式中：s——任意的非零尺度因子；——圖像平面上的二維點(diǎn)坐標(biāo)的齊次坐標(biāo)；A——相機(jī)內(nèi)部參數(shù)矩陣；——平面模板上的三維點(diǎn)坐標(biāo)的齊次坐標(biāo)。

(4)

假設(shè)有n幅標(biāo)定圖像，每一幅標(biāo)定圖像上有m個(gè)特征點(diǎn)，則建立非線性優(yōu)化模型為

(5)

將每幅圖解得的A,Ri,Ti及每點(diǎn)對應(yīng)的圖像和世界坐標(biāo)點(diǎn)對，畸變參數(shù)k1,k2的初值為0，通過非線性優(yōu)化得到全局最優(yōu)解。

2.2 立體匹配

如圖2所示，左圖為物體變形前散斑圖像，右圖為物體變形后的散斑圖像，在左圖中選取以某待求點(diǎn)P(x0,y0)為中心的像素尺寸為(2M+1)×(2M+1)的矩形參考圖像子區(qū)，在變形后圖像中選取若干個(gè)相同大小的樣本子區(qū)進(jìn)行逐一比對匹配，按標(biāo)準(zhǔn)化協(xié)方差互相關(guān)系數(shù)的計(jì)算結(jié)果，以尋找與協(xié)方差相關(guān)系數(shù)絕對值為最大時(shí)對應(yīng)的目標(biāo)圖像子區(qū)點(diǎn)P′(x′0,y′0)，確定x和y方向的位移分量。在利用數(shù)字圖像進(jìn)行實(shí)際計(jì)算時(shí)，以虛擬網(wǎng)格的形式對測量圖像進(jìn)行劃分，求得每個(gè)虛擬網(wǎng)格節(jié)點(diǎn)的位移來獲取全場三維應(yīng)變場。

圖2 物體變形前后參考圖像子區(qū)和目標(biāo)圖像子區(qū)Fig.2 Reference image sub-region and target image sub-region before and after object deformation

變形前后圖像子區(qū)之間的協(xié)方差相關(guān)系數(shù)C為

(6)

一般來說，與參考圖像子區(qū)相比，目標(biāo)圖像子區(qū)的中心位置和形狀都會(huì)發(fā)生改變。在采用一階函數(shù)的情況下，參考圖像子區(qū)中的某一點(diǎn)Q(x,y)與目標(biāo)圖像子區(qū)中對應(yīng)點(diǎn)Q′(x′,y′)之間存在一定的函數(shù)關(guān)系，即形函數(shù)為

(7)

式中：Δx和Δy——點(diǎn)(x,y)到P(x0,y0)的距離；u和v——參考子區(qū)中心在x、y方向上的位移；ux,uy和vx,vy——圖像子區(qū)的位移梯度。

由公式(7)求出變形后圖像搜索區(qū)域中所有目標(biāo)子區(qū)中心點(diǎn)在x、y方向上的位移和位移梯度的值，就可以得到模型表面的三維位移量和變形值。由于高速相機(jī)采集得到的測量圖像灰度信息是離散的，因此在目標(biāo)子區(qū)搜索時(shí)位移只能以整像素為單位來進(jìn)行，然而物體真實(shí)的位移值或變形值并不一定是像素的整數(shù)倍，采用亞像素搜索算法，通過對相關(guān)系數(shù)矩陣做曲面擬合運(yùn)算，由擬合曲面的極值點(diǎn)解算得到亞像素級別的位移量，提高被測模型三維位移和變形的測量精度。

3 應(yīng)力測量系統(tǒng)設(shè)計(jì)

采用雙目立體視覺和數(shù)字圖像相關(guān)技術(shù)可實(shí)現(xiàn)低溫環(huán)境下被測試件應(yīng)力參數(shù)的現(xiàn)場測量，測量裝置由兩臺(tái)高速相機(jī)、散斑噴涂裝置、被測試件、圖像處理計(jì)算機(jī)及測量軟件組成。其中，兩臺(tái)高速相機(jī)型號為V1612，最高幀頻可達(dá)16000幀/s，分辨率為1280×800，像元尺寸為28μm，均可滿足測量圖片的高速高分辨采集的要求，其應(yīng)力測量示意圖如圖3所示。

圖3 應(yīng)力測量示意圖Fig.3 Schematic diagram of stress measurement

測量時(shí)，先對雙目測量系統(tǒng)進(jìn)行標(biāo)定，確定相機(jī)內(nèi)參數(shù)和兩臺(tái)相機(jī)坐標(biāo)系間的空間轉(zhuǎn)換關(guān)系(相機(jī)外參數(shù))，構(gòu)建全局坐標(biāo)系，求解被測模型坐標(biāo)系與全局坐標(biāo)系的空間關(guān)系。根據(jù)測量需要，可在風(fēng)洞的觀察窗處架設(shè)一套雙目視覺系統(tǒng)。測量過程中，基于電路輸出TTL信號給兩相機(jī)的外觸發(fā)接口，同步觸發(fā)裝置同時(shí)觸發(fā)兩臺(tái)高速攝像機(jī)。在被測模型上噴涂散斑圖案，利用兩臺(tái)高速數(shù)字相機(jī)同時(shí)采集變形前后的散斑圖像。被測試件在風(fēng)洞內(nèi)環(huán)境和溫度變化過程中，被測系統(tǒng)以16kHz的采樣頻率對試件的空間姿態(tài)進(jìn)行實(shí)時(shí)測量，兩個(gè)高速相機(jī)同步采集測量圖像，傳輸?shù)綌?shù)據(jù)處理系統(tǒng)進(jìn)行分析和解算，通過立體匹配精確獲得左右相機(jī)對應(yīng)子區(qū)中心點(diǎn)的精確匹配，結(jié)合雙目立體視覺測量模型解算各點(diǎn)在全局坐標(biāo)系下的空間三維坐標(biāo)。

在高馬赫數(shù)風(fēng)洞試驗(yàn)條件下，不可避免會(huì)產(chǎn)生抖動(dòng)現(xiàn)象，對視覺系統(tǒng)的姿態(tài)測量精度會(huì)產(chǎn)生影響。在風(fēng)洞試驗(yàn)程中被測模型發(fā)生的運(yùn)動(dòng)實(shí)際上是一個(gè)抖動(dòng)和姿態(tài)變換的復(fù)合運(yùn)動(dòng)，本文采用運(yùn)動(dòng)分離的方法，通過監(jiān)視被測模型固定架自身的抖動(dòng)，把該抖動(dòng)量從飛行器實(shí)際運(yùn)動(dòng)中消除，獲得被測模型的運(yùn)動(dòng)參數(shù)。針對風(fēng)洞試驗(yàn)的具體環(huán)境，在風(fēng)洞側(cè)壁安裝固定參考點(diǎn)，通過固定點(diǎn)反推相機(jī)的抖動(dòng)量和飛行器的抖動(dòng)量，將抖動(dòng)量從復(fù)合運(yùn)動(dòng)中減掉，即可獲得被測模型真實(shí)的姿態(tài)變化。

在變形前，可選取左相機(jī)測量圖像中要計(jì)算的對應(yīng)圖像子區(qū)，通過立體匹配可得到右相機(jī)對應(yīng)的圖像子區(qū)，基于視覺原理可得到該子區(qū)中心點(diǎn)在T0時(shí)刻的空間三維坐標(biāo)(x0,y0,z0)；變形后，基于圖像數(shù)字相關(guān)法分別找到左、右相機(jī)對應(yīng)的目標(biāo)子區(qū)，得到T1時(shí)刻變形后的子區(qū)中心點(diǎn)的空間三維坐標(biāo)(x1,y1,z1)；變形前后圖像子區(qū)中心點(diǎn)的三維坐標(biāo)的差，即為所求的三維位移，平滑后通過差分計(jì)算求得相應(yīng)的應(yīng)變場。在應(yīng)力參數(shù)測量中，風(fēng)洞條件下應(yīng)力不能直接測量得到，但由力引起的應(yīng)變可獲取得到應(yīng)力量。

4 三維數(shù)字圖像相關(guān)法算法

立體匹配是雙目視覺中的難點(diǎn)，主要涉及兩種匹配：第一種是左相機(jī)測量圖像集的相關(guān)匹配，與二維數(shù)字圖像相關(guān)法中的相關(guān)匹配原理相同；第二種是左右相機(jī)測量圖像集之間的相關(guān)匹配，兩相機(jī)的測量圖像存在視差。整個(gè)方案的立體匹配算法流程圖如下圖4所示。

圖4 立體匹配算法流程圖Fig.4 Stereo matching algorithm flowchart

首先，對左、右相機(jī)進(jìn)行內(nèi)外參數(shù)標(biāo)定，讀取左右相機(jī)圖像，以左相機(jī)第一幅測量圖像作為左相機(jī)測量圖像系列匹配的參考圖像，也作為右相機(jī)測量圖像集的參考圖像，實(shí)現(xiàn)左相機(jī)圖像間的匹配、及左右相機(jī)對應(yīng)圖像間的匹配。以計(jì)算點(diǎn)為中心，截取一定大小的參考圖像子區(qū)和目標(biāo)圖像進(jìn)行Shift初值估計(jì)，迭代收斂計(jì)算。若滿足迭代收斂條件，則利用形函數(shù)計(jì)算得到該計(jì)算點(diǎn)的上、下、左、右四個(gè)點(diǎn)的初值估計(jì)點(diǎn)，計(jì)算所有網(wǎng)格點(diǎn)，通過三維重建得第k幅圖像計(jì)算點(diǎn)的三維坐標(biāo)值。

5 試驗(yàn)與分析

試驗(yàn)裝置由薄壁圓形平板試件、兩臺(tái)高速相機(jī)、同步觸發(fā)裝置、高低溫試驗(yàn)箱和圖像處理計(jì)算機(jī)組成。其中，試件材料為5A06鋁合金，其彈性模量為70GPa，試件直徑為160mm，厚度為2mm。將試件放置于高低溫實(shí)驗(yàn)箱內(nèi)，在試件表面噴涂散斑，架設(shè)雙相機(jī)同步采集測量圖像，通過立體匹配算法進(jìn)行圖像子區(qū)的相關(guān)匹配，獲得不同溫度條件下試件上多點(diǎn)的三維坐標(biāo)值，并解算得到試件在-30℃時(shí)多點(diǎn)的應(yīng)變值，被測試件在常溫下和-30℃下的點(diǎn)的三維坐標(biāo)值數(shù)據(jù)、被測試件在溫度下降到-30℃時(shí)多點(diǎn)應(yīng)變值數(shù)據(jù)見表1所示。由變形前后試件表面的三維坐標(biāo)值通過后續(xù)解算可得到試件表面的應(yīng)變場。圖5(a)為在常溫下，試件表面11個(gè)點(diǎn)的三維坐標(biāo)值解算圖，當(dāng)溫度下降為-30攝氏度時(shí)，通過立體匹配算法對測量圖像進(jìn)行多點(diǎn)的相關(guān)匹配，計(jì)算得到試件表面多點(diǎn)的三維坐標(biāo)值解算圖，如圖5(b)所示。

圖5 常溫下和-30℃下試件表面多點(diǎn)的三維坐標(biāo)值Fig.5 Coordinate values of multiple points on the surface of the specimen at room temperature and-30℃

由試驗(yàn)可得，基于三維數(shù)字圖像相關(guān)算法實(shí)現(xiàn)了左相機(jī)測量圖像集的相關(guān)匹配，及左右相機(jī)測量圖像集之間的相關(guān)匹配，完成了被測試件表面點(diǎn)的立體匹配及應(yīng)力測量。

6 結(jié)束語

本文探討了基于高速三維圖像數(shù)字相關(guān)法的模型應(yīng)力測量方法，結(jié)合立體匹配算法可獲得變形前后物體上計(jì)算點(diǎn)的三維坐標(biāo)值，可得模擬表面的應(yīng)力場分布，為解決低溫風(fēng)洞下模型狀態(tài)測量問題提供一種可行的測量方法。

表1 常溫下和-30℃下被測試件多點(diǎn)的坐標(biāo)值及多點(diǎn)應(yīng)力值Tab.1 Solutioncoordinatevalueofthespecimenatroomtemperatureand-30℃andstrainvaluesofmultiplepoints常溫下x(mm)y(mm)z(mm)-30℃x(mm)y(mm)z(mm)應(yīng)變值應(yīng)力值(GPa)點(diǎn)132.23536.9840.119點(diǎn)132.29536.8430.0560.075%0.0525點(diǎn)235.4835.2670.096點(diǎn)235.5225.1580.0130.106%0.0742點(diǎn)336.944-31.8920.009點(diǎn)336.979-31.978-0.0810.195%0.1365點(diǎn)4-0.00265.7700.114點(diǎn)40.06365.6310.0710.092%0.0644點(diǎn)50.58923.8220.232點(diǎn)50.65023.7040.1840.153%0.1071點(diǎn)62.142-15.6330.163點(diǎn)62.189-15.7240.1340.144%0.1008點(diǎn)7-3.384-55.0060.003點(diǎn)7-3.344-55.089-0.0290.207%0.1449點(diǎn)8-42.55957.344-0.052點(diǎn)8-42.47357.226-0.1330.111%0.0777點(diǎn)9-44.65016.3660.158點(diǎn)9-44.56116.2420.0850.092%0.0644點(diǎn)10-35.365-22.6020.064點(diǎn)10-35.297-22.6980.0290.114%0.0798點(diǎn)11-75.895-9.206-0.093點(diǎn)11-75.825-9.299-0.1110.130%0.0910