发明名称 基于编码结构光的高速运动目标位姿测量方法
摘要 本发明基于编码结构光的高速运动目标位姿测量方法属于计算机视觉测量技术领域,涉及一种大视场小目标高速运动物体位姿测量方法。测量方法采用彩色投影仪向测量区域投射含有编码信息的彩色栅线,使用高速摄像机连续拍摄经运动物体表面调制而变形的结构光编码栅线图案;经过数字图像解码后,实现成像平面与投影平面的对应点匹配,获取相位变化信息,由系统结构参数反算出高速运动物体的三维动态形貌;通过将测量点与基于待测物几何参数建立的先验模型进行匹配,最终准确地测量出高速运动物体的位置、姿态信息。本发明采用彩色条纹伪随机序列与灰度栅线相位相结合的方式,有效解决了立体视觉中的匹配难题。测量系统成本低,机构简单、操作简易。
申请公布号 CN105180904A 申请公布日期 2015.12.23
申请号 CN201510604065.1 申请日期 2015.09.21
申请人 大连理工大学 发明人 刘巍;陈玲;贾振元;刘惟肖;马鑫
分类号 G01C11/00(2006.01)I 主分类号 G01C11/00(2006.01)I
代理机构 大连理工大学专利中心 21200 代理人 关慧贞
主权项 一种基于编码结构光的高速运动目标位姿测量方法,采用单目高速摄像机和投影仪,其特征是,测量方法采用彩色投影仪向测量区域投射含有编码信息的彩色栅线,使用高速摄像机连续拍摄经运动物体表面调制而变形的结构光编码栅线图案;经过数字图像解码后,实现成像平面与投影平面的对应点匹配,获取相位变化信息,由系统结构参数反算出高速运动物体的三维动态形貌;通过将测量点与基于待测物几何参数建立的先验模型进行匹配,最终准确地测量出高速运动物体的位置、姿态信息;测量方法的具体步骤如下:1、彩色栅线编码设计采用一种彩色条纹伪随机序列与灰度栅线相位相结合的彩色栅线空域编码图案,在单幅图像中构成像素级别唯一编码;2、彩色栅线图像解码彩色栅线图像解码即求解每一像素点对应的有效完全相位φ(x,y):包括获取相展开,获取周期序列号n和相主值<img file="FDA0000807426280000011.GIF" wi="188" he="98" />折叠在0~2π中两部分,如公式(1)所示;<img file="FDA0000807426280000012.GIF" wi="1297" he="148" />采用由粗到细的分步彩色栅线图像解码技术,在粗识别阶段根据周期边界及色彩信息划分周期区域进行相展开,在细分阶段根据局部区域亮度余弦变化规律进行空域解相,即解算相主值;具体解码过程如下:1)由RGB表示模式转化为HIS表示模式根据各颜色通道亮度大小及差值比例关系将RGB表示模式转换为HIS表示模式,将色彩信息与亮度信息分离;R、G、B分别表示RGB模式下某颜色在三个通道的色调值,色调值H对应彩色栅线图案中的六种颜色,分六种情况计算;用于获取颜色信息,从而确定周期序列号,色调的计算公式如下:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>R</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>B</mi><mo>,</mo><mi>H</mi><mo>=</mo><mrow><mo>(</mo><mi>G</mi><mo>-</mo><mi>B</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>R</mi><mo>-</mo><mi>B</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>G</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>B</mi><mo>,</mo><mi>H</mi><mo>=</mo><mn>2</mn><mo>-</mo><mrow><mo>(</mo><mi>R</mi><mo>-</mo><mi>B</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>G</mi><mo>-</mo><mi>B</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>G</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>R</mi><mo>,</mo><mi>H</mi><mo>=</mo><mn>2</mn><mo>+</mo><mrow><mo>(</mo><mi>B</mi><mo>-</mo><mi>R</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>G</mi><mo>-</mo><mi>R</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>B</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>R</mi><mo>,</mo><mi>H</mi><mo>=</mo><mn>4</mn><mo>-</mo><mrow><mo>(</mo><mi>G</mi><mo>-</mo><mi>R</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>B</mi><mo>-</mo><mi>R</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>B</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>G</mi><mo>,</mo><mi>H</mi><mo>=</mo><mn>4</mn><mo>+</mo><mrow><mo>(</mo><mi>R</mi><mo>-</mo><mi>G</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>B</mi><mo>-</mo><mi>G</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>max</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>R</mi><mo>,</mo><mi>min</mi><mrow><mo>(</mo><mi>R</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>B</mi><mo>)</mo></mrow><mo>=</mo><mi>G</mi><mo>,</mo><mi>H</mi><mo>=</mo><mn>6</mn><mo>-</mo><mrow><mo>(</mo><mi>B</mi><mo>-</mo><mi>G</mi><mo>)</mo></mrow><mo>/</mo><mrow><mo>(</mo><mi>R</mi><mo>-</mo><mi>G</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000807426280000021.GIF" wi="1621" he="595" /></maths>将求得的色调定义为六种色彩,当H值处于0~0.5或者5.5~6区间时为红色,处于0.5~1.5区间时为黄色,处于1.5~2.5区间时为绿色,处于2.5~3.5区间时为青色,处于3.5~4.5区间时为蓝色,处于4.5~5.5区间时为品红色,6个颜色区间分别对应周期序列号0‑5;亮度值I余弦变换规律用于相主值的求解;其计算公式为:<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>I</mi><mo>=</mo><mfrac><mrow><mi>R</mi><mo>+</mo><mi>G</mi><mo>+</mo><mi>B</mi></mrow><mn>3</mn></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000807426280000022.GIF" wi="992" he="149" /></maths>2)对彩色栅线图像进行相展开在彩色栅线编码的伪随机序列中,只要提取出高亮度中心边界和其左右相邻两个周期高亮度中心边界色彩信息,即可确定中心边界对应的周期顺序号,实现周期解码;即解算出完全相位中的周期序列号n值;3)对彩色栅线图像进行相主值解算通过等步长相移法对彩色栅线图像进行局域空域解相,即对相主值<img file="FDA0000807426280000031.GIF" wi="190" he="93" />进行求解;假设连续5个像素对应相同的物体表面反射率r,环境光分量a,条纹调制幅度b和基本畸变相位<img file="FDA0000807426280000032.GIF" wi="79" he="66" />每两个像素间的相位步长均为2ε;则有:<img file="FDA0000807426280000033.GIF" wi="1267" he="494" />其中:<img file="FDA0000807426280000034.GIF" wi="1591" he="517" />反正切值在0~π/2之间;待入测量亮度值进行求解展开至0~2π之间,取两者平均值解出<img file="FDA0000807426280000035.GIF" wi="81" he="66" />最终根据公式(1)实现完全解相;3、高速运动物体三维动态形貌解算通过预先标定好的测量系统参数,反算出高速运动物体测量表面高度信息,进而得到三维动态形貌,计算公式如下:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mtable><mtr><mtd><mrow><mi>A</mi><mi>B</mi><mo>=</mo><mi>p</mi><mi>l</mi><mo>&CenterDot;</mo><mrow><mo>(</mo><msub><mi>&phi;</mi><mi>A</mi></msub><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>)</mo><mo>-</mo><msub><mi>&phi;</mi><mi>B</mi></msub><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mi>y</mi></mrow><mo>)</mo><mo>)</mo></mrow><mo>/</mo><mn>2</mn><mi>&pi;</mi><mo>,</mo></mrow></mtd><mtd><mrow><mi>h</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><mi>A</mi><mi>B</mi><mo>&CenterDot;</mo><mi>L</mi></mrow><mrow><mi>A</mi><mi>B</mi><mo>+</mo><mi>d</mi></mrow></mfrac></mrow></mtd></mtr></mtable><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000807426280000036.GIF" wi="1542" he="149" /></maths>投影光线照射到参考面A点,成像平面上对应的像素点其相位为φ<sub>A</sub>(x,y),放上被测物体后,像素点对应的是投影至物体表面D点的光线,该光线投影至参考平面的B点,相位为φ<sub>B</sub>(x,y);d为投影出瞳P和摄影入瞳C的距离,L为摄影入瞳C与参考面的距离,pl为参考面上投影光栅的周期长度;距离AB包含了高度信息h(x,y),而该距离可由两点间的相位差计算;解算出测量表面的高度信息,最终获得高速运动物体的三维动态形貌;4、基于先验模型的三维测量点匹配与位姿求解分为预估阶段和匹配求解两个阶段:1)预估阶段①通过测量数据点集与先验模型数据点集的预校准,提高匹配算法的可靠性;在预校准前,首先定义两个指标V和K<sub>u</sub>,V表示P的特征矩阵的主要特征值,反应点的分布情况;K<sub>u</sub>是表示点集分布集中化程度的一个指标,对起始测量数据点集的可靠度进行预估;V=λ<sub>max</sub>(Cov(P))   (8)其中P={p<sub>i</sub>}   (9)表示将要进行预校准的测量点数据集;Cov(P)=E[(p‑μ)(p‑μ)<sup>T</sup>]   (10)μ=E[p]   (11)K<sub>u</sub>定义为<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><msub><mi>K</mi><mi>u</mi></msub><mo>=</mo><mfrac><mn>1</mn><msup><mi>V</mi><mn>2</mn></msup></mfrac><mo>{</mo><mfrac><mn>1</mn><msub><mi>N</mi><mi>p</mi></msub></mfrac><munderover><mo>&Sigma;</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><msub><mi>N</mi><mi>p</mi></msub></munderover><msup><mrow><mo>(</mo><msub><mi>p</mi><mi>i</mi></msub><mo>-</mo><mi>&mu;</mi><mo>)</mo></mrow><mn>4</mn></msup><mo>}</mo><mo>-</mo><mn>3</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000807426280000041.GIF" wi="1173" he="198" /></maths>N<sub>p</sub>表示集合P中点的数量;如果有很多点沿着一个特定的方向分布,V将会取大值,并且相应的特征向量表示出点分布的主要方向;当K<sub>u</sub>=0时表明符合正态分布,这样的点分布是可靠的;结合这两个指数,定义具有相对较大的V值和相对较小的K<sub>u</sub>值的测量点分布是可靠的;第一预估阶段将依次循环直至在时序测量点集中找到第一个可信点集,并定义为P<sub>1</sub>;对起始测量数据点集(P<sub>1</sub>,P<sub>2</sub>)的预估结果存在两种可能性:如果P<sub>1</sub>和P<sub>2</sub>均被认定为可靠点集,就用Cov(P<sub>1</sub>)、Cov(P<sub>2</sub>)的特征向量和模型点集X与测量点集P<sub>1</sub>的中心点,分别对X的姿态与位置进行预校准;如果P<sub>2</sub>在第二个预估阶段被认定为不可靠,那么预校准阶段将会采用在第一个时序图像P<sub>1</sub>中获取的相对位置与姿态信息作为估计值对先验模型点集X进行预校准;②先验模型的前表面模型建立;通过预校准时所使用的信息建立摄像机可视部分的先验模型点集代替完整的先验模型点集,有效的提高匹配算法的可靠性与精度;这一局部先验模型称为前表面模型;2)匹配求解阶段完成总流程中的预校准与前表面模型的建立后,利用所得数据进行匹配与位姿求解;P<sup>(m)</sup>为测量数据点集,X<sup>(m)</sup>为位置与姿态经预校准的前表面先验模型数据点集;上标m表示循环次数;其中,最优的平移变换矩阵T<sup>(m)</sup>和旋转变换矩阵R<sup>(m)</sup>通过基于特征值问题的封闭求解算法直接进行估算,得到表示平移变换矩阵T的相对位置的<img file="FDA0000807426280000051.GIF" wi="150" he="77" />和表示旋转变换矩阵R的四元数矩阵<img file="FDA0000807426280000052.GIF" wi="79" he="77" />该算法将一直循环下去,直至达到收敛条件E≤T<sub>E</sub>或者循环次数达到预先设定的阈值m≤T<sub>m</sub>,最终得到该时刻目标物体的相对位置与相对姿态信息;通过执行上述计算过程直到k=k<sub>max</sub>,即完成对所有时序图像对应的测量数据点集与先验模型点集的匹配,解算出每一时刻高速运动物体的相对姿态<img file="FDA0000807426280000061.GIF" wi="60" he="84" />与相对位置<img file="FDA0000807426280000062.GIF" wi="236" he="82" />进而从时序图象中获取高速运动物体的位姿状态信息。
地址 116024 辽宁省大连市甘井子区凌工路2号