风洞水压试验洞体结构变形量及应变量测量方法及装置与流程

1.本发明涉及风洞设计与测试技术领域，具体地，涉及风洞水压试验洞体结构变形量及应变量测量方法及装置。


背景技术：

2.在航空、航天器研制过程中，都需要在风洞等基础试验设备中模拟测试。在风洞设计与建造过程中，为了测试洞体结构是否满足设计载荷强度要求，通常需要进行风洞水压试验。在水压试验过程中，需要在线测量洞体结构变形及应变。目前，我国在风洞水压试验过程中，主要采用经纬仪、激光跟踪仪等测量设备离散地测量洞体结构大变形，采用应变计监测洞体结构局部应变变形。这些离散测量方式，一次只能测量一个点的数据，对于大尺寸风洞而言，工作量非常大、效率非常低。采用应变计进行结构应变测量，需要预先估计存在应变极点的位置，如果估计失误，将会遗漏真实的应变极点位置。并且，一个应变计只能监测一个点的应变，针对大尺寸结构，通常需要粘贴大量的应变计，而带来试验准备时间和测试成本的大幅增加。
3.综上所述，发明人在研究现有技术中发现现有的洞水压试验洞体结构变形量及应变量测量方法至少存在以下技术问题：
4.(1)测量设备离散地测量洞体结构变形数据，测量效率低。
5.(2)采用应变计进行结构应变测量，容易测量不准确，以及应变测量时间大和测量成本高。


技术实现要素：

6.为了提高风洞水压试验洞体结构变形量的测量效率，本发明提供了一种风洞水压试验洞体结构变形量测量方法。
7.具体方法包括：
8.在风洞水压试验前，在洞体结构表面绘制标识图案；
9.基于第一图像采集模块和第二图像采集模块构建立体视觉测量系统；
10.对所述立体视觉测量系统进行标定，建立立体视觉测量系统坐标系；
11.在风洞水压试验前，利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第一图像和第二图像；
12.在风洞水压试验中，利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第三图像和第四图像；
13.在所述第一图像中选择第一像素点，以所述第一像素点坐标为中心获得预设尺寸大小的参考图像；
14.基于所述参考图像进行图像匹配，分别在所述第二图像、所述第三图像和所述第四图像中匹配出与所述第一像素点对应的第二像素点、第三像素点和第四像素点；
15.基于所述立体视觉测量系统的标定结果，计算获得所述第一像素点和所述第二像
素点对应的第一三维点坐标，以及计算获得所述第三像素点和所述第四像素点对应的第二三维点坐标；
16.基于所述第一三维点坐标和所述第二三维点坐标，获得所述第一像素点对应的三维空间位移数据；
17.以第一像素点为起点，预设步长，遍历所述第一图像中的像素点，获得被遍历像素点对应的三维空间位移数据；
18.基于所述第一像素点对应的三维空间位移数据和所述被遍历像素点对应的三维空间位移数据，获得所述洞体结构的变形量。
19.本方法采用立体视觉成像系统记录洞体结构变形过程中标识图案，通过图像块匹配方法，和立体视觉三维重建，可以一次性计算出洞体结构变形场，提高了风洞水压试验洞体结构变形量的测量效率。
20.其中，在洞体结构表面绘制标识图案可以便于图像采集模块拍摄图像，然后对图像中的标识进行识别和匹配，进而计算出图像中的像素对应的三维坐标，通过三维坐标位移量计算出洞体结构变形量。
21.优选的，本方法中的图像匹配方式为：
22.基于图像块匹配得到所述参考图像在被匹配图像中的匹配位置；
23.以所述匹配位置为初始位置进行亚像素匹配，获得与所述第一像素点对应的像素点。
24.优选的，本方法在亚像素匹配过程中，采用二阶泰勒展开对图像块变形进行近似计算。
25.优选的，本方法图像块亚像素匹配计算流程包括：
26.第一步：基于所述匹配位置，确定所述被匹配图像的中心坐标；基于所述参考图像的中心坐标和所述被匹配图像的中心坐标，初始化参数p6和p
12
，其中，p1至p5和p7至p
11
初始化为0，p1＝u
i
，p2＝u
j
，p3＝u
ii
，p4＝u
jj
，p5＝u
ij
，p6＝u
o
，p7＝v
i
，p8＝v
j
，p9＝v
ii
，p
10
＝v
jj
，p
11
＝v
ij
，p
12
＝v
o
，(i
o
,j
o
)为参考图像的中心坐标，(i,j)为参考图像中某一像素的坐标，u
o
,v
o
为(i
o
,j
o
)在图像中的位移量，u
i
,v
i
是位移u,v相对于坐标i的1阶导数，u
j
,v
j
分别为位移u,v相对于坐标j的1阶导数，u
ii
,v
ii
分别为位移u,v相对于坐标j的2阶偏导数，u
jj
,v
jj
分别为位移u,v相对于坐标j的2阶偏导数，u
ij
,v
ij
分别为位移u,v相对于坐标i,j的2阶偏导数；设置迭代终止条件ε或d，其中，ε为相似度量c的阈值，d为迭代次数阈值；其中，所述参考图像f的中心坐标为所述被匹配图像g的中心坐标为所述被匹配图像g的中心坐标为所述被匹配图像g的中心坐标为
27.第二步：计算所述参考图像f的1阶梯度基于所述分别计算海森矩阵h
12
×
12
、图像梯度与变形模型偏微分积计算所述参考图像f的图像均值f
m
和图像方差σ
f
，计算所述被匹配图像g的图像均值g
m
和图像方差σ
g
；u,v,w是分别是(x,y,z)方向上的位移，p为变形模型参数，p＝(p1,...,p
12
)，为偏微分运算符号；
28.第三步：基于所述h
12
×
12
、所述所述f
m
、所述σ
f
、所述g
m
和所述σ
g
，计算增量参数
△
p和图像变形函数w(ζ；p)；
29.第四步：基于所述f
m
、所述σ
f
、所述g
m
、所述σ
g
和所述w(ζ；p)，计算第m次迭代对应的相似度量c
m
，判断c
m
是否满足收敛条件：当满足c
m
<ε或满足m>d时，停止迭代，并进入第五步，否则，迭代次数加1并返回第三步继续迭代；
30.第五步：基于所述图像变形函数w(ζ；p)中的参数p，得到图像块亚像素匹配的位移量：u0＝p6,v0＝p
12
，u0为水平方向上的位移量，v0为竖直方向上的位移量。
31.其中，相似性度量可采用newton
–
raphson(nr)迭代算法进行优化，但是该算法采用正向变形映射，在每一次迭代过程中，都需要计算变形后图像的梯度和海森矩阵，计算量较大，算法效率低。为此，本发明采用icgn算法对目标函数进行求解。icgn算法不同于nr算法，它执行逆向变形操作，部分参数都可以提前计算，而无需在每一次迭代过程中计算，进而减少计算量，提高计算效率。
32.优选的，本方法中海森矩阵其中，s为图像像素集合，其中，为向量乘法操作，ζ为位移向量，ζ＝[
△
i
ꢀ△
j (
△
i
)
2 (
△
j
)2ꢀ△
i
△
j 1]
t
，w(ζ；p)指将所述被匹配图像g变形成所述参考图像f的二阶变形操作，w(ζ；
△
p)为所述参考图像f映射到所述被匹配图像g的变形过程，γ＝[
△
i
,
△
j
]
t
为参考图像局部坐标，
△
i
＝i
‑
i
o
，
△
j
＝j
‑
j
o
，
△
i
为像素在水平方向的位移，
△
j
为像素在竖直方向的位移，f(γ)为参考图像中坐标为γ的像素值，g(γ)为被匹配图像g中坐标为γ的像素值，f(γ w(ζ；
△
p))为参考图像f以γ为参考点进行坐标偏移后的像素值，g(γ w(ζ；p))为被匹配图像g以γ为参考点进行坐标偏移后的像素值。
[0033]
为解决传统风洞水压试验洞体结构应变量测量方法存在的测量不准确，以及应变测量时间大和测量成本高，本发明提供了一种风洞水压试验洞体结构应变量测量方法，所述方法包括：
[0034]
基于所述的风洞水压试验洞体结构变形量测量方法，获得所述洞体结构的变形量；
[0035]
基于所述洞体结构局部区域内的三维空间位移数据，以1阶多项式近似所述洞体结构的变形量，计算所述洞体结构变形量的导数，基于所述洞体结构变形量的导数获得所述洞体结构应变量。
[0036]
本发明提供的风洞水压试验中洞体结构应变场测量方法，通过非接触测量方式，可测量洞体结构三维应变场，相比应变计只能测量局部区域应变，本发明可一次性测量洞体结构全场应变，测量效率高，且无需预先估计存在应变极点的位置测量准确，并且不需要使用到大量的应变计，测量成本低。
[0037]
采用局部区域最小平方差拟合方法计算应变场的原理是：应变是位移的导数，因此位移可表示为应变的1阶多项式。通过标识图案块匹配和立体视觉三维重建，可计算出被测对象表面三维位移量。但是，这些位移数据中包含有噪声，直接对位移求导数用于计算应变会产生较大的测量误差。对于大型风洞结构而言，通常情况下，在局部邻域内，它的变形
是平滑的。也就是说：在某一局部邻域中，应变是相似的。因此，可以通过局部区域内位移数据，在1阶多项式下，进行最小平方误差拟合，估计应变测量值，以减小测量误差。
[0038]
优选的，本方法中：
[0039]
所述第一图像中设有若干种子点，计算获得每个所述种子点的位移量，基于所有所述种子点的位移量获得所有种子点的位移点云数据，基于所述位移点云数据的导数获得所述洞体结构应变量。
[0040]
优选的，本方法中，每个所述种子点的位移量计算过程包括：
[0041]
在所述第一图像中选择种子点；
[0042]
基于所述种子点，获得所述种子点的预设尺寸大小的第一局部邻域；
[0043]
在风洞水压试验前，获得所述第一局部邻域的第一三维坐标点云；
[0044]
在风洞水压试验中，获得所述第一局部邻域的第二三维坐标点云；
[0045]
基于所述第一三维坐标点云、所述第二三维坐标点云和所述种子点位移量的1阶导数拟合得到所述种子点的位移量。
[0046]
由应变定义可知，若知道被测对象表面的位移量，通过微分操作，得到位移的导数，即可计算格林应变量。但是，微分对位移数据中的噪声非常敏感，容易产生错误的应变数据。洞体结构变形在某一局部邻域内是平滑的，即：该局部邻域内位移的导数是相同的。因此，可利用局部邻域内位移数据估计位移的导数。
[0047]
优选的，本方法中基于所述第一三维坐标点云、所述第二三维坐标点云和所述种子点位移量的1阶导数，在最小平方误差准则下拟合得到所述种子点的位移量。
[0048]
以当前计算点为中心定义一个局部邻域，进行最小平方差准则下的位移差分量拟合，可有效消除位移中噪声对位移差分计算结果的干扰。按此方式，逐点处理每一个点云数据，即可得到全场位移差分结果，用于计算全场应变。
[0049]
本发明还提供了一种风洞水压试验洞体结构变形量测量装置，所述装置包括：
[0050]
第一图像采集模块、第二图像采集模块、计算机和图像采集控制器；其中，洞体结构表面绘制有标识图案，所述第一图像采集模块和所述第二图像采集模块的光轴相交，所述第一图像采集模块和所述第二图像采集模块用于对所述洞体结构进行立体成像，所述计算机用于启动所述图像采集控制器，所述图像采集控制器用于控制所述第一图像采集模块和所述第二图像采集模块同步采集所述标识图案的实时图像，并将所述实时图像传输至所述计算机，所述计算机进行图像数据处理，计算获得所述洞体结构的变形量；
[0051]
其中，在风洞水压试验前，基于所述第一图像采集模块和所述第二图像采集模块构建立体视觉测量系统；所述计算机用于对所述立体视觉测量系统进行标定，建立立体视觉测量系统坐标系；利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第一图像和第二图像；
[0052]
在风洞水压试验中，利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第三图像和第四图像；所述计算机用于在所述第一图像中选择第一像素点，以所述第一像素点坐标为中心获得预设尺寸大小的参考图像；所述计算机基于所述参考图像进行图像匹配，分别在所述第二图像、所述第三图像和所述第四图像中匹配出与所述第一像素点对应的第二像素点、第三像素点和第四像素点；所述计算机基于所述立体视觉测量系统的标定结果，计算获得所述第一像素点和所述第二像素点对应的
第一三维点坐标，以及计算获得所述第三像素点和所述第四像素点对应的第二三维点坐标；所述计算机基于所述第一三维点坐标和所述第二三维点坐标，获得所述第一像素点对应的三维空间位移数据；所述计算机以第一像素点为起点，遍历所述第一图像中的像素点，获得被遍历像素点对应的三维空间位移数据；所述计算机基于所述第一像素点对应的三维空间位移数据和所述被遍历像素点对应的三维空间位移数据，获得所述洞体结构的变形量。
[0053]
本发明提供的一个或多个技术方案，至少具有如下技术效果或优点：
[0054]
1)本发明提供的风洞水压试验中洞体结构变形及应变场测量装置及方法，通过非接触测量方式，可测量洞体结构三维变形场及应变场，相比激光跟踪仪测量方法，可一次性测量洞体结构全场变形，相比应变计只能测量局部区域应变，本发明可一次性测量洞体结构全场应变，测量效率高。
[0055]
2)本发明装置利用摄像机进行非接触测量，试验准备简单、可重复使用，相比应变计测量，使用成本更低，且测量更加准确。
[0056]
3)本发明采用二阶泰勒展开近似图像变形，可描述透视畸变和结构变形引起的斑点图像位移变化，并采用icgn算法降低图像匹配过程中的数据计算量，提高了算法效率。
[0057]
4)本发明利用一阶多项式拟合局部邻域内点云位移场，用于计算应变场，相比直接利用位移差分方法计算应变场，具有更好的噪声抑制能力，确保应变场计算精度。
附图说明
[0058]
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本发明的一部分，并不构成对本发明实施例的限定；
[0059]
图1为风洞水压试验洞体结构变形测量装置的结构示意图；
[0060]
图2为斑点图像块匹配示意图；
[0061]
图3为洞体结构变形测量原理示意图；
[0062]
图4为斑点图像正向变形示意图；
[0063]
图5为斑点图像逆向变形示意图；
[0064]
图6为位移导数估计示意图；
[0065]
图7为斑点图案实例；
[0066]
图中：1.洞体结构，2.斑点图案，3.摄像机a，4.摄像机b，5.图像采集控制器，6.计算机，7.光源，8.参考斑点图像块，9.待匹配斑点图像，11.t0时刻摄像机a拍摄图像，12.t＝0时刻摄像机b拍摄图像，13.t＝n时刻摄像机a拍摄图像，14.t＝n时刻摄像机b拍摄图像，15.在图像11中选择的种子点像素，16.以种子点像素11为中心定义的大小为b*b的斑点图像块，17.种子点像素15在图像12中亚像素匹配结果，18.种子点像素15在图像13中亚像素匹配结果，19.种子点像素15在图像14中亚像素匹配结果，20.图像11中选取的种子点，21.种子点20的邻域，22.种子点20及邻域21内像素在t＝0时刻的三维点云构成的面片，23.由种子点20及邻域21内像素在t＝n时刻的三维点云构成的面片。
具体实施方式
[0067]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实
施方式对本发明进行进一步的详细描述。需要说明的是，在相互不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0068]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述范围内的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
[0069]
本领域技术人员应理解的是，在本发明的揭露中，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系是基于附图所示的方位或位置关系，其仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此上述术语不能理解为对本发明的限制。
[0070]
可以理解的是，术语“一”应理解为“至少一”或“一个或多个”，即在一个实施例中，一个元件的数量可以为一个，而在另外的实施例中，该元件的数量可以为多个，术语“一”不能理解为对数量的限制。
[0071]
实施例一
[0072]
本发明实施例一提供了一种风洞水压试验洞体结构变形与应变场测量装置，如图1所示，该装置包括：洞体结构1表面绘制的斑点图案2、摄像机a3、摄像机b4、图像采集控制器5、计算机6、光源7；其中，计算机6与图像采集控制器5、摄像机a3、摄像机b4连接，图像采集控制器5与摄像机a3、摄像机b4连接，摄像机a3、摄像机b4光轴相交安置，用于对洞体结构1进行立体成像；计算机6中上位机程序启动图像采集控制器5，触发摄摄像机a3、摄像机b4同步拍摄洞体结构1表面的斑点图案2，并将拍摄的斑点图案图像传入计算机6，进行图像数据处理，计算洞体结构1的变形及应变场。
[0073]
其中，本发明实施例中的标识图案可以为斑点图案，也可以为花纹图案，也可以为其他具有标识作用的图案，本发明对标识图案的具体实现方式及形式不进行限定。
[0074]
其中，本发明实施例中，第一图像采集模块和第二图像采集模块可以是相机，也可以是摄像机，也可以是其他具有拍照或拍摄功能的图像采集设备，本发明对图像采集模块的具体实现方式不进行限定。
[0075]
其中，本装置中光源可以根据需要进行选配，若光线条件较好则可以不使用光源，若光线条件一般为了保障图像拍摄的质量可以增加光源。
[0076]
该装置的工作原理是：通过在洞体结构1表面绘制散乱的斑点图案2，采用立体视觉成像系统记录洞体结构变形过程中斑点图像，通过斑点图像块匹配方法，和立体视觉三维重建，可以一次性计算出洞体结构变形场，再根据变形场可计算出洞体结构应变场。图2给出了斑点图像块匹配示意图，在图2中，通过斑点图像块匹配，可在待匹配图像9中找到参考斑点图像块8的最优匹配结果。
[0077]
洞体结构变形场测量原理如图3所示，其中，图像11和图像12是t＝0时刻(即水压试验前，洞体结构未变形时)摄像机a和摄像机b分别拍摄的斑点图像，图像11和图像12构成一组立体图像对；图像13和图像14是t＝n时刻(即水压试验过程中，记录下的洞体结构变形图像)摄像机a和摄像机b分别拍摄的斑点图像，同样图像13和14构成一组立体图像对。根据立体视觉成像原理，给定立体视觉测量系统内、外参数，以及立体图像对中同名像素在两张图像中的图像坐标，即可计算该像素的三维空间坐标，进而实现三维测量。给定同名像素在
t＝0时刻和t＝n时刻的三维空间坐标，即可得到该同名像素的三维空间位移量。
[0078]
洞体结构变形场测量方法包括：首先，在水压试验前，在洞体结构1表面绘制斑点图案2，对摄像机a3、摄像机b4构建的立体视觉测量系统进行标定，修正镜头畸变，建立测量系统坐标系，在t＝0时刻，摄像机a3拍摄洞体结构1表面斑点图案2得到图像11，摄像机b4拍摄洞体结构1表面斑点图案2得到图像12，图像11和图像12构成一组立体图像对{11,12}；接着，在水压试验中，摄像机a、摄像机b同步拍摄t＝n时刻洞体结构1表面斑点图案2的立体图像对{13,14}；在图像11中选择一个像素点15作为第一种子点，以该像素点15的图像坐标(i
a
,j
a
)为中心定义一个大小为b*b的图像块16作为该像素的斑点图像参考模板，b的取值范围为1～10000，单位为像素，然后通过斑点图像块匹配方法分别在图像12、13、14中进行亚像素匹配，亚像素匹配结果分别是像素17、像素18、像素19，对应的亚像素图像坐标分别是(i
b
,j
b
)、(i
c
,j
c
)、(i
d
,j
d
)，再根据立体视觉测量系统标定结果，计算像素15和像素17对应的三维点坐标p(x,y,z)，以及像素18和像素19对应的三维点坐标p'(x',y',z')，两者相减，即为像素15对应的洞体结构1表面上某一点从t＝0时刻到t＝n时刻的三维空间位移：(u,v,w)＝(x',y',z')
‑
(x,y,z)；以步长s1遍历整幅图像，s1取值范围为0～1000个像素，得到2组洞体结构变形前后三维坐标点云:[x,y,z]、[x',y',z']，[u,v,w]＝[x',y',z']
‑
[x,y,z]是洞体结构的三维位移点云，也是洞体结构变形量。
[0079]
三维变形场具体计算步骤如下：
[0080]
第一步：在图像11左上角区域(或其他区域)内设置一个起始种子点a(即像素15)，图像坐标为(i
a
,j
a
)，以(i
a
,j
a
)为中心，定义一个b*b大小的图像块q(图像块16)；
[0081]
第二步：以图像块q作为参考斑点图像块(8)，采用斑点图像块匹配方法(10)，在图像12中进行亚像素匹配，得到亚像素匹配图像坐标(i
b
,j
b
)(即像素17)；
[0082]
第三步：以图像块q作为参考斑点图像块(8)，采用斑点图像块匹配方法(10)，在图像13中进行亚像素匹配，得到亚像素匹配图像坐标(i
c
,j
c
)(即像素18)；
[0083]
第四步：以图像块q作为参考斑点图像块(8)，采用斑点图像块匹配方法(10)，在图像(14)中对其进行亚像素匹配，得到亚像素匹配图像坐标(i
d
,j
d
)(即像素19)；
[0084]
第五步：根据立体视觉测量系统标定结果，以起始种子点a在图像(11)中的图像坐标(i
a
,j
a
)和在图像(12)中亚像素匹配坐标(i
b
,j
b
)，计算该起始种子点a在t＝0时刻的三维空间坐标以起始种子点a在图像(13)中的亚像素匹配图像坐标(i
c
,j
c
)和在图像(14)中亚像素匹配坐标(i
d
,j
d
)，计算该起始种子点a在t＝n时刻的三维空间坐标
[0085]
第六步：以起始种子点a在t＝n时刻的三维空间坐标减去t＝0时刻的三维空间坐标得到起始种子点a在t＝n时刻相对于t＝0时刻的三维空间位移：
[0086]
沿着图像水平和竖直方向(可以先水平后竖直，也可以先竖直后水平，或两个方向交替进行，只要能够遍历完图像即可，本实施例对具体的方向不进行限定)，以步长s1遍历图像11即可得到被测对象的三维坐标点云及变形量。
[0087]
上述斑点图像块匹配方法采用整像素、亚像素两级匹配实现，具体如下：
[0088]
a)整像素匹配：采用lewis等(j.p.lewis,fast normalized cross
‑
correlation,industrial light&magic.1995)等提出的基于fft快速傅里叶变换的图像块快速匹配方法，得到图像块q在匹配图像(即图像12、13、14)中的整像素匹配位置m；
[0089]
b)亚像素匹配：以整像素位置m作为初始位置，采用baker等(baker s,matthews i.lucas
‑
kanade 20years on:a unifying framework.int j comput vis 2004；56(3):221
–
55)提出的icgn,inverse compositional gauss
‑
newton光流算法，进行亚像素匹配。在匹配过程中，采用二阶泰勒展开对图像块变形进行近似。
[0090]
采用二阶泰勒展开近似图像块变形的方法如下：
[0091]
如图4所示，给定一张参考图像f(比如图3中的图像块16)，局部坐标系为(
△
i,
△
j)，坐标系原点位于图像中心(i
o
,j
o
)处，其变形后(变形类型包括透视畸变或结构变形)的图像为图像g(以后称之为当前图像)。那么，参考图像f中某一像素(i,j)，在g中的图像坐标可表示为：
[0092][0093]
其中u,v分别是竖直、水平方向上的位移量。在中心点(i
o
,j
o
)处，对u,v进行2阶泰勒展开，可较好地近似该位移量。
[0094][0095]
其中，
△
i
＝i
‑
i
o
，
△
j
＝j
‑
j
o
，其中u
o
,v
o
是(i
o
,j
o
)在图像中的位移量，u
i
,v
i
是位移u,v相对于坐标i的1阶导数，u
j
,v
j
是位移u,v相对于坐标j的1阶导数，u
ii
,v
ii
是位移u,v相对于坐标j的2阶偏导数，u
jj
,v
jj
是位移u,v相对于坐标j的2阶偏导数，u
ij
,v
ij
是位移u,v相对于坐标i,j的2阶偏导数。
[0096]
式(2)可表示成一个12参数p∈r
12
的变形函数：
[0097][0098]
其中，p1＝u
i
，p2＝u
j
，p3＝u
ii
，p4＝u
jj
，p5＝u
ij
，p6＝u
o
，p7＝v
i
，p8＝v
j
，p9＝v
ii
，p
10
＝v
jj
，p
11
＝v
ij
，p
12
＝v
o
，t为矩阵转置操作符。
[0099]
定义ζ＝[
△
i
ꢀ△
j (
△
i
)
2 (
△
j
)2ꢀ△
i
△
j 1]
t
，p＝(p1,...,p
12
)，式(3)的二阶变形变换过程可由变换w(ζ,p)表示，ζ位移向量，p参数向量。
[0100]
当前图像g可由参考图像f通过变换w(ζ,p)变形得到：
[0101]
g(γ
o
γ)＝f(γ
o
w(ζ,p))
ꢀꢀꢀꢀꢀ
(4)
[0102]
其中，γ
o
＝[i
o
,j
o
]
t
为参考图像f的中心坐标,γ＝[
△
i
,
△
j
]
t
是图像局部坐标，g(γ
o
γ)为当前图像局部邻域，f(γ
o
w(ζ,p))为参考图像变形后局部邻域。
[0103]
当参考图像f和当前图像g中所有像素一一对应时，在不考虑整体亮度变化的情况下，它们之间的像素差异应该最小，即两者的平方差接近于0：
[0104][0105]
其中，s图像像素集合，f(γ)为参考图像中坐标为γ的像素值，g(γ)为被匹配图像g中坐标为γ的像素值。
[0106]
在实际应用中，可能存在亮度变化干扰，为此对参考图像f和当前图像g进行归一化处理后，计算它们的相似性：
[0107][0108]
其中，c是参考图像f和当前图像g的相似性度量准则，c越小，两者越相似，分别是参考图像f和当前图像g的均值，n(s)为s中像素数量，m为迭代次数，分别是参考图像f和当前图像g的方差。式(6)的相似性度量准则，可较好地消除光照变化干扰。
[0109]
根据式(4)、式(6)的相似性准则可写为：
[0110][0111]
式(7)可采用newton
–
raphson(nr)迭代算法进行优化，但是该算法采用图4所示的正向变形映射，在每一次迭代过程中，都需要计算变形后图像的梯度和海森矩阵，计算量较大，算法效率低。
[0112]
为此，本发明采用baker等(baker s,matthews i.lucas
‑
kanade 20years on:a unifying framework.int j comput vis 2004；56(3):221
–
55)提出的icgn算法对目标函数(7)进行求解。icgn算法不同于nr算法，它执行逆向变形操作，如图5所示，将当前图像g通过变换g(γ w(ζ；p))变形成参考图像f，即f(γ)
∝
g(γ w(ζ；p))(注：该式表示f(γ)与g(γ w(ζ；p))正相关，也就是说使趋近0)。那么，采用icgn算法优化时，图像匹配相似性度量的目标函数可写为：
[0113][0114]
其中，w(ζ,p)为把当前图像g映射到参考图像f的变形过程；
[0115]
在迭代优化前，对参考图像f进行微小的变形f(γ)'＝f(γ w(ζ；
△
p))，使f(γ)'＝f(γ w(ζ；
△
p))
∝
g(γ w(ζ；p))，f(γ)'为参考图像f变形后图像，因此式(8)可写成：
[0116][0117]
其中，其中w(ζ；p)中的参数p是前一次迭代后的结果，即该参数已知，w(ζ；
△
p)是
把参考图像f映射到当前图像g的变形过程，是未知参数，
△
p1为p1的导数，
△
p2为p2的导数，
△
p3为p3的导数，
△
p4为p4的导数，
△
p5为p5的导数，
△
p6为p6的导数，
△
p7为p7的导数，
△
p8为p8的导数，
△
p9为p9的导数，
△
p
10
为p
10
的导数，
△
p
11
为p
11
的导数，
△
p
12
为p
12
的导数，通过找到使c最小化的w(ζ；
△
p)，来计算参数
△
p，再通过逆变换w(ζ；
△
p)
‑1，即可找到将当前图像g变换成参考图像f的变换过程：，即可找到将当前图像g变换成参考图像f的变换过程：为向量乘法操作。
[0118]
c最小化的解为：
[0119][0120]
其中，其中，是参考图像f的1阶梯度，s为图像像素集合。
[0121]
是w(ζ；p)的偏微分，为偏微分运算符号，具体如下：
[0122][0123]
其中，
△
k为w(ζ；
△
p)变形过程中，横轴方向上的位移量，
△
l为w(ζ；
△
p)变形过程
中，纵轴方向上的位移量，海森矩阵图像梯度与变形模型偏微分积图像均值和方差f
m
,g
m
,σ
f
,σ
g
等都可以提前计算，而无需在每一次迭代过程中计算，进而减少计算量，提高计算效率。
[0124]
下面给出基于icgn算法进行斑点图像块亚像素匹配的计算流程：
[0125]
第一步：根据整像素模板匹配结果m，确定当前图像g的中心坐标并结合参考图像f的中心坐标初始化参数其余参数设为0，设置迭代终止条件ε或d，其中ε是相似度量c的阈值，当第m次迭代的相关性度量结果c
m
<ε时，停止迭代；或者，迭代次数m超过d时，停止迭代；
[0126]
第二步：以参考图像f计算图像1阶梯度以此计算海森矩阵图像梯度与变形模型偏微分积并计算参考图像f和当前图像g的图像均值和方差f
m
,g
m
,σ
f
,σ
g
；
[0127]
第三步：根据式(10)计算增量参数
△
p、图像变形函数
[0128]
第四步：计算判断是否满足收敛条件：c
m
<ε或者m>d；当满足收敛条件时，停止迭代，进入第五步，否则返回第3步继续迭代；
[0129]
第五步：根据优化后的变形函数w(ζ；p)中的参数p，得到亚像素匹配的位移量：u0＝p6,v0＝p
12
。
[0130]
洞体结构应变场测量方法是：根据被测对象三维位移点云[u,v,w]，采用局部区域最小平方差拟合方法，以1阶多项式近似三维位移，计算三维位移导数，并由此计算格林应变。
[0131]
采用局部区域最小平方差拟合方法计算应变场的原理是：应变是位移的导数，因此位移可表示为应变的1阶多项式。通过斑点图像块匹配和立体视觉三维重建，可计算出被测对象表面三维位移量。但是，这些位移数据中包含有噪声，直接对位移求导数用于计算应变会产生较大的测量误差。对于大型风洞结构而言，通常情况下，在局部邻域内，它的变形是平滑的。也就是说：在某一局部邻域中，应变是相似的。因此，可以通过局部区域内位移数据，在1阶多项式下，进行最小平方误差拟合，估计应变测量值，以减小测量误差。
[0132]
根据弹性力学理论，在(x,y,z)坐标系下，格林(green)应变定义为：
[0133][0134]
其中，u,v,w是分别是(x,y,z)方向上的位移，是偏微分，e
xx
为x轴方向应变，e
yy
为y轴方向应变，e
zz
为z轴方向应变，e
xy
为xy平面内应变，e
yz
为yz平面内应变，e
zx
为zx平面内应变。由式(12)的应变定义可知，若知道被测对象表面的位移量(u,v,w)，通过微分操作，得到位移的导数即可计算格林应变量。但是，微分对位移数据中的噪声非常敏感，容易产生错误的应变数据。洞体结构变形在某一局部邻域内是平滑的，即：该局部邻域内位移的导数是相同的。因此，可利用局部邻域内位移数据估计位移的导数
[0135]
如图6所示，在图像(11)中取一个第二种子点(20)，定义其邻域为21(在图6中，该邻域窗口大小为3*3，以后将该局部邻域称之为面片)，该邻域内第二种子点在t＝0、t＝n时刻的三维点云分别构成两个面片(22、23)，其中白色圆环表示t＝0时刻的三维坐标点云，黑色圆点表示t＝n时刻的三维坐标点云。第二种子点的位移量(u,v,w)可通过面片21中第二种子点坐标及位移1阶导数拟合得到：
[0136][0137][0138][0139]
其中，u,v,w是t＝n时刻第二种子点的位移，u0,v0,w0是第二种子点20在t＝n时刻相对于t＝0时刻的位移量，
△
x＝x
‑
x0,
△
y＝y
‑
y0,
△
z＝z
‑
z0是t＝n时刻局部窗口内第二种子点相对于中心位置处第二种子点20的偏移量，(x0,y0,z0)是t＝n时刻种子点20的三维坐标。
[0140]
令a0＝u0,
[0141]
则式(13)可写为：
[0142]
u＝a0 a1x a2y a3z
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0143]
令b0＝v0,
[0144]
则式(14)可写为：
[0145]
v＝b0 b1x b2y b3z
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)
[0146]
令c0＝w0,
[0147]
则式(15)可写为：
[0148]
w＝c0 c1x c2y c3z
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0149]
在面片21内有多个点，则得到下面的方程组：
[0150][0151][0152][0153]
其中n是第二种子点数量。
[0154]
式(19)可写成矩阵形式：
[0155][0156]
式(20)可写成矩阵形式：
[0157][0158]
式(21)可写成矩阵形式：
[0159][0160]
令为邻域21内n个像素的局部坐标,n的取值范围为0～1000；
[0161]
是邻域内n个像素点的位移u；
[0162]
是邻域内n个像素点的位移v；
[0163]
是邻域内n个像素点的位移w；
[0164]
是位移u的0、1阶导数；
[0165]
是位移v的0、1阶导数；
[0166]
是位移w的0、1阶导数；
[0167][0168]
则式(22
‑
24)可写为：
[0169]
u＝ka
ꢀꢀꢀꢀꢀꢀꢀ
(25)
[0170]
v＝kb
ꢀꢀꢀꢀꢀꢀꢀ
(26)
[0171]
w＝kc
ꢀꢀꢀꢀꢀꢀꢀ
(27)
[0172]
在窗口21内，有9个变形点(位移)数据，可用于估计参数a,b,c。以计算u为例，当位移点云数量为4时，方程14有解。由于(u,v,w)含有噪声，因此直接利用四个点，计算的参数a
并不可靠。在最小平方误差准则下估计最优参数
[0173][0174]
其中，u是实际测量值，||
·
||是u
‑
ka的欧式距离，是在最小平方误差准则下的最优参数，它的最优解为：
[0175][0176]
其中(k
t
k)
‑1k
t
是k的伪逆矩阵。
[0177]
同理，位移场v,w的最优参数分别为：
[0178][0179][0180]
得到最优参数后，根据式(12)计算应变场。
[0181]
以当前计算点(比如20)为中心定义一个面片，进行最小平方差准则下的位移差分量拟合，可有效消除位移中噪声对位移差分计算结果的干扰。按此方式，逐点处理每一个点云数据，即可得到全场位移差分结果，用于计算全场应变。在上述操作过程中，局部窗口20的大小根据实际应用情况设定。
[0182]
实施例二
[0183]
在本实施例二中，选用2台分辨率为400万像素、帧率为10fps的黑白工业相机，镜头参数分别为：光圈f＝2.0，焦距f＝13mm，也可以选用其他相机或相机的参数调整设置为其他方式，本实施例不进行具体的限定；选用功率为50w的led光源为摄像机提供照明；摄像机通过线缆与图像采集控制器连接，所拍摄的图像通过camera
‑
link数据采集卡输入计算机中进行数据处理，计算洞体结构变形及应变场。测量对象为风洞试验段洞壁，被测区域为400mm*400mm的矩形。
[0184]
在水压试验前，采用喷漆法在洞体结构表面制作斑点图案。具体方法是：先用白色油漆均匀喷涂洞体结构作为底色，然后选用黑色油漆，油漆的颜色本实施例不进行具体的限定，两者能够形成区分即可，在距离洞体结构50cm处微量喷洒油漆，喷出的油漆液滴随机散落在白色洞体结构上，形成杂乱纹理的斑点图案，如图7所示。
[0185]
将摄像机a、摄像机b组成的立体视觉测量系统安置在被测洞体结构正前方，调整摄像机角度，确保被测洞体结构可在立体视觉系统中成像。
[0186]
采用jean
‑
yves bouguet提供的摄像机标定工具包camera calibration toolbox(下载地址为http://www.vision.caltech.edu/bouguetj/)，对由摄像机a、摄像机b构成的立体视觉测量系统进行标定，校正镜头畸变，得到摄像机a、摄像机b的内、外参数，建立立体视觉测量系统坐标系。
[0187]
在水压试验前，摄像机a、摄像机b分别拍摄图像11、12，作为t＝0时刻洞体结构斑点图案立体图像对{11,12}。
[0188]
在水压试验过程中，摄像机a、摄像机b分别拍摄图像13、14，作为t＝n时刻洞体结构斑点图案立体图像对{13,14}。
[0189]
获得观测图像数据后，下面对图像11、12、13、14进行处理，以计算被测洞体结构的变形及应变场。
[0190]
在变形场计算过程中，设置图像块大小40*40，即b＝40，遍历步长s1＝2，在icgn算法中迭代终止条件设为ε＝0.001、d＝50；在应变场计算过程中，设置邻域窗口(即面片)大小为3*3。
[0191]
最终，获得了20.15万个三维变形点云数据，并由这20.15万三维变形点云计算出了18.2万个应变测量值。
[0192]
在被测区域中心处粘贴0.1％精度的应变计作为对比测试，对比结果表明测量精度可以达到0.26％。
[0193]
实施例三
[0194]
本发明实施例三提供了风洞水压试验洞体结构变形量测量装置，所述装置包括：
[0195]
第一图像采集模块、第二图像采集模块、计算机和图像采集控制器；其中，洞体结构表面绘制有标识图案，所述第一图像采集模块和所述第二图像采集模块的光轴相交，所述第一图像采集模块和所述第二图像采集模块用于对所述洞体结构进行立体成像，所述计算机用于启动所述图像采集控制器，所述图像采集控制器用于控制所述第一图像采集模块和所述第二图像采集模块同步采集所述标识图案的实时图像，并将所述实时图像传输至所述计算机，所述计算机进行图像数据处理，计算获得所述洞体结构的变形量；
[0196]
其中，在风洞水压试验前，基于所述第一图像采集模块和所述第二图像采集模块构建立体视觉测量系统；所述计算机用于对所述立体视觉测量系统进行标定，建立立体视觉测量系统坐标系；利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第一图像和第二图像；
[0197]
在风洞水压试验中，利用所述第一图像采集模块和所述第二图像采集模块分别采集所述标识图案的图像分别获得第三图像和第四图像；所述计算机用于在所述第一图像中选择第一像素点，以所述第一像素点坐标为中心获得预设尺寸大小的参考图像；所述计算机基于所述参考图像进行图像匹配，分别在所述第二图像、所述第三图像和所述第四图像中匹配出与所述第一像素点对应的第二像素点、第三像素点和第四像素点；所述计算机基于所述立体视觉测量系统的标定结果，计算获得所述第一像素点和所述第二像素点对应的第一三维点坐标，以及计算获得所述第三像素点和所述第四像素点对应的第二三维点坐标；所述计算机基于所述第一三维点坐标和所述第二三维点坐标，获得所述第一像素点对应的三维空间位移数据；所述计算机以第一像素点为起点，遍历所述第一图像中的像素点，获得被遍历像素点对应的三维空间位移数据；所述计算机基于所述第一像素点对应的三维空间位移数据和所述被遍历像素点对应的三维空间位移数据，获得所述洞体结构的变形量。
[0198]
洞体结构变形场测量方法是：首先，在水压试验前，在洞体结构1表面绘制斑点图案2，对摄像机a3、摄像机b4构建的立体视觉测量系统进行标定，修正镜头畸变，建立测量系统坐标系，在t＝0时刻，摄像机a3拍摄洞体结构1表面斑点图案2得到图像11，摄像机b4拍摄洞体结构1表面斑点图案2得到图像12，图像11和图像12构成一组立体图像对{11,12}；接着，在吹风试验中，摄像机a、摄像机b同步拍摄t＝n时刻洞体结构1表面斑点图案2的立体图像对{13,14}；在图像11中选择一个像素点15作为种子点，以该像素点15的图像坐标(i
a
,j
a
)为中心定义一个大小为b*b的图像块16作为该像素的斑点图像参考模板，b的取值范围为1～10000，单位为像素，然后通过斑点图像块匹配方法分别在图像12、13、14中进行亚像素匹
配，亚像素匹配结果分别是像素17、像素18、像素19，对应的亚像素图像坐标分别是(i
b
,j
b
)、(i
c
,j
c
)、(i
d
,j
d
)，再根据立体视觉测量系统标定结果，计算像素15和像素17对应的三维点坐标p(x,y,z)，以及像素18和像素19对应的三维点坐标p'(x',y',z')，两者相减，即为像素15对应的洞体结构1表面上某一点从t＝0时刻到t＝n时刻的三维空间位移：(u,v,w)＝(x',y',z')
‑
(x,y,z)；以步长s1遍历整幅图像，s1取值范围为0～1000个像素，得到2组洞体结构变形前后三维坐标点云:[x,y,z]、[x',y',z']，[u,v,w]＝[x',y',z']
‑
[x,y,z]是洞体结构的三维位移点云，也是洞体结构变形量；
[0199]
洞体结构应变场测量方法是：根据被测对象三维位移点云[u,v,w]，采用局部区域最小平方差拟合方法，以1阶多项式近似三维位移，计算三维位移导数，并由此计算格林应变；
[0200]
具体计算方法是：在三维位移点云[u,v,w]中，选择一个种子点20，定义种子点20的邻域为21，该邻域内所有像素点的位移量表示为1阶位移导数的多项式：
[0201]
u＝ka
[0202]
v＝kb
[0203]
w＝kc
[0204]
其中为邻域21内n个像素的局部坐标,n的取值范围为0～1000；
[0205]
是邻域内n个像素点的位移u；
[0206]
是邻域内n个像素点的位移v；
[0207]
是邻域内n个像素点的位移w；
[0208]
是位移u的0、1阶导数；
[0209]
是位移v的0、1阶导数；
[0210]
是位移w的0、1阶导数；
[0211]
计算其中(k
t
k)
‑1k
t
是k的伪逆矩阵；
[0212]
采用下式计算种子点20位置处格林(green)应变：
[0213][0214]
按上述方法，遍历所有三维位移点云[u,v,w]，计算洞体结构应变场；
[0215]
所述斑点图像块匹配方法，由整像素、亚像素两级匹配完成：其中，整像素匹配：采用基于fft快速傅里叶变换的图像块快速匹配方法，得到参考图像块16在待匹配图像中的整像素匹配位置m，待匹配图像包括图像12、图像13、图像14；亚像素匹配：以整像素位置m作为初始位置，使用icgn(inverse compositional gauss
‑
newton)逆组分牛顿高斯光流算法，采用二阶泰勒展开近似图像块变形，进行亚像素匹配；
[0216]
在icgn(inverse compositional gauss
‑
newton)逆组分牛顿高斯光流算法框架下，采用二阶泰勒展开近似图像块变形，进行亚像素匹配的方法如下：
[0217]
采用二阶泰勒展开近似图像块变形的函数是：
[0218][0219]
其中，i0,j0是图像块中心坐标，是变形后坐标，
△
i
＝i
‑
i
o
，
△
j
＝j
‑
j
o
是当前像素相对于中心点(i0,j0)的偏移量，取值范围为0～10000，p1＝u
i
是位移u相对于坐标i的1阶导数，p2＝u
j
是位移u相对于坐标j的1阶导数，p3＝u
ii
是位移u相对于坐标i的2阶导数，p4＝u
jj
是位移u相对于坐标j的2阶导数，p5＝u
ij
是位移u相对于坐标i,j的偏导数，p6＝u
o
是参考图像g中心点坐标i的位移，p7＝v
i
是位移v相对于坐标i的1阶导数，p8＝v
j
是位移v相对于坐标j的1阶导数，p9＝v
ii
是位移v相对于坐标i的2阶导数，p
10
＝v
jj
是位移v相对于坐标j的2阶导数，p
11
＝v
ij
是位移v相对于坐标i,j的偏导数，p
12
＝v
o
是参考图像g中心点坐标j的位移；p1～p12这12个参数的取值范围为
‑
256～256；
[0220]
图像块匹配的相似性度量目标函数为：
[0221][0222]
其中，c是参考图像f和当前图像g的相似性度量准则，
分别是参考图像f和当前图像g的均值，γ是像素坐标，s表示图像块中像素集，分别是参考图像f和当前图像g的方差，ζ＝[
△
i
ꢀ△
j (
△
i
)
2 (
△
j
)2ꢀ△
i
△
j 1]
t
，p＝(p1,...,p
12
)，w(ζ；p)指将当前图像g变形成参考图像f的二阶变形操作，w(ζ；
△
p)是把参考图像f映射到当前图像g的变形过程；
[0223]
亚像素匹配的操作步骤如下：
[0224]
第一步，根据整像素模板匹配结果m，确定当前图像g的中心坐标并结合参考图像f的中心坐标初始化参数其余参数设为0，并设置迭代终止条件ε或d，ε的取值范围为0～1000，d的取值范围为1～100000；
[0225]
第二步，以参考图像f计算图像1阶梯度以此计算海森矩阵图像梯度与变形模型偏微分积并计算参考图像f和当前图像g的图像均值和方差f
m
,g
m
,σ
f
,σ
g
；
[0226]
第三步，计算第三步，计算
[0227]
第四步，计算判断是否满足收敛条件：c
m
<ε或者m>d，m是迭代次数；当满足收敛条件时，停止迭代，进入第五步，否则返回第3步继续迭代；
[0228]
第五步，根据优化后的变形函数w(ζ；p)中的参数p，得到亚像素匹配的位移量：u0＝p6,v0＝p
12
。
[0229]
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0230]
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。