CN113658266A - 一种基于固定相机和单靶标的动轴旋转角视觉测量方法 - Google Patents

一种基于固定相机和单靶标的动轴旋转角视觉测量方法 Download PDF

Info

Publication number
CN113658266A
CN113658266A CN202110850964.5A CN202110850964A CN113658266A CN 113658266 A CN113658266 A CN 113658266A CN 202110850964 A CN202110850964 A CN 202110850964A CN 113658266 A CN113658266 A CN 113658266A
Authority
CN
China
Prior art keywords
target
camera
coordinate system
angle
axis
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110850964.5A
Other languages
English (en)
Other versions
CN113658266B (zh
Inventor
范军华
刘贡平
杨锋
王萍
穆志国
方强
王彦哲
王丽斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
AVIC Xian Aircraft Industry Group Co Ltd
Original Assignee
Zhejiang University ZJU
AVIC Xian Aircraft Industry Group Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU, AVIC Xian Aircraft Industry Group Co Ltd filed Critical Zhejiang University ZJU
Priority to CN202110850964.5A priority Critical patent/CN113658266B/zh
Publication of CN113658266A publication Critical patent/CN113658266A/zh
Application granted granted Critical
Publication of CN113658266B publication Critical patent/CN113658266B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/80Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C3/00Measuring distances in line of sight; Optical rangefinders
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/60Rotation of whole images or parts thereof
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10004Still image; Photographic image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Electromagnetism (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Image Analysis (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了一种基于固定相机和单靶标的动轴旋转角视觉测量方法,适用于飞机襟翼类后退下拉式动轴活动舵面偏转角度测量。该方法包含三大步骤:步骤一:相机内外参数标定;步骤二:动轴旋转角度数学模型建立;步骤三:利用标定的结果和建立的数学模型,测量任意位置运动旋转部件的旋转角度。本发明是一种利用单目相机对绕动轴旋转运动的旋转角进行测量的方法。它解决了现有视觉测角方法仅适用于定轴转角测量,无法进行动轴转角测量的局限性,能够实现沿任意轨迹运动的动轴的旋转角度测量。同时,本发明的方法对相机和靶标无安装要求,操作简便,具有良好的实用性。

Description

一种基于固定相机和单靶标的动轴旋转角视觉测量方法
技术领域
本发明涉及一种基于固定相机和单靶标的动轴旋转角视觉测量方法,典型实施例适用于飞机襟翼类后退下拉式活动舵面偏转角度测量。
背景技术
角位移测量广泛应用于航空航天等工业领域。在飞机生产装配中,飞机舵面角位移传感器需要调整安装位置和校正参数。视觉测量技术由于其非接触性、实时性,以及安装接单,测量便捷等优点,受到广泛关注,在该领域具有较高的应用前景。
现有的视觉测角方法主要分为两大类,一类是基于双目视觉或多目视觉的测量方法,例如:申请号为202010769781.6的发明专利介绍了一种基于双目相机的三维测距方法,此类测量方法采用两台或多台相机拍摄目标照片,存在视场受多目相机配置的限制而较小,而某些飞机舵面运动行程较大,故其在飞机舵面转角测量领域某些情况下受到一定限制;另一类是基于单目视觉的测量方法,仅利用一台相机拍摄照片来进行测量工作,因其仅需要一台视觉传感器,故此类方法的优点为相机安装简便、标定简单,且避免了立体视觉中的视场较小、目标立体匹配较难等缺点。
在现有的基于单目视觉的测角该方法中,申请号为201210594277.2的发明专利“一种单目空间目标测距测角方法”介绍的单目视觉方法对待测目标提出了一定要求,必须已知待测目标的三维结构和尺寸信息,且测量过程需要基于预先建立的待测目标多尺度多姿态模板图像库,方法准备过程繁琐,不具有普适性,不适于工业应用;申请号为201410506279.0的发明专利“一种基于固定相机和单靶标的单轴旋转角的视觉测量方法”介绍了一种基于单目相机和单靶标的单轴旋转角的视觉测量方法,但存在仅适用于回转轴固定的单轴旋转角度测量,对于回转轴非固定的动轴运动场景,其算法和模型无效,无法计算动轴转角;申请号为201410506277.1的发明专利“一种基于移动相机和双靶标的单轴旋转角的视觉测量方法”介绍的方法同样仅适用于定轴,无法测量动轴转角。
飞机活动舵面的运动方式主要分为两类,一类是定轴转动,一类是后退下拉。后退下拉运动舵面为本发明所述的典型的动轴回转运动。以飞机襟翼为例,襟翼的后退下拉过程是通过滑轮-滑轮架实现与固定结构的连接,在液压或全电驱动器的驱动下推动滑轮沿滑轮架上的直线滑轨运动实现舵面的平移,同时绕舵面与滑轮连接处的铰链实现舵面旋转。平移、旋转运动耦合为舵面实际后退下拉过程。上述发明中的视觉测量方法,无法解决该典型场景下,襟翼回转角度测量问题。
因此,本专利针对现有视觉测角方法不适用于动轴转角测量的局限性,申请提出了一种基于固定相机和单靶标的运动轴旋转角视觉测量方法,该方法可实现沿任意轨迹运动的动轴的旋转角度测量,方法无靶标特殊安装要求,安装快捷、操作简便。
发明内容
本发明提出了一种基于固定相机和单靶标的动轴旋转角视觉测量方法,该方法是一种利用单目相机对绕动轴转动的旋转角进行测量的方法,它解决了现有视觉测量定轴转角方法无法测量动轴转角的问题,能够实现沿任意轨迹运动的动轴的旋转角度测量。同时,本发明的方法对相机和靶标无安装要求,操作简便,具有良好的实用性。
一种基于固定相机和单靶标的动轴旋转角视觉测量方法,该方法的具体步骤如下:
步骤一:相机内外参标定。
具体实现过程中,靶标采用二维棋盘格形式。标定过程采用的图像序列包括:一张基准位置的靶标图片和m张旋转部件分别旋转θi(i=1,2,...,m)角度时的靶标图片。利用已广泛采用的成熟技术包括:角点提取算法和张正友相机标定算法,提取标定图像序列棋盘格角点,并进行内外参标定,求得相机内参Min和m+1个位姿的外参
Figure BDA0003182319400000031
Figure BDA0003182319400000032
为从基准位置的靶标坐标系到相机坐标系的旋转矩阵和平移向量;
Figure BDA0003182319400000033
为从旋转θi(i=1,2,...,m)的靶标坐标系到相机坐标系的旋转矩阵和平移向量。
步骤二:动轴旋转角度数学模型建立。
(1)将
Figure BDA0003182319400000034
转换到单位四元数形式
Figure BDA0003182319400000035
其中:
Figure BDA0003182319400000036
式中,
Figure BDA0003182319400000037
表示
Figure BDA0003182319400000038
中(x,y)索引所对应的数值。
令单位四元数qi0表示从
Figure BDA0003182319400000039
Figure BDA00031823194000000310
的旋转关系,则从Oc-XcYcZc
Figure BDA00031823194000000311
再到
Figure BDA00031823194000000312
的连续旋转过程可由下式获得:
Figure BDA00031823194000000313
其中,
Figure BDA00031823194000000314
Figure BDA00031823194000000315
表示同一靶标点分别在
Figure BDA00031823194000000316
Figure BDA00031823194000000317
中对应的三维坐标,qi表示
Figure BDA00031823194000000318
的变换过程,即从
Figure BDA00031823194000000319
到相机坐标系的单位四元数变换。qi0表示从
Figure BDA00031823194000000320
Figure BDA00031823194000000321
的单位四元数变换。
由于qi已知,则qi0可由下式获得:
Figure BDA0003182319400000041
其中,
Figure BDA0003182319400000042
为q0的共轭四元数。对于单位四元数qi0,下式成立。
Figure BDA0003182319400000043
其中,[r1,r2,r3]T表示从
Figure BDA0003182319400000044
Figure BDA0003182319400000045
的坐标变换过程所绕的回转轴的单位向量。θ表示从
Figure BDA0003182319400000046
Figure BDA0003182319400000047
的坐标变换过程中绕回转轴[r1,r2,r3]T的旋转角度。
因此,从上式可得靶标各位置相对于基准位置的旋转角度为
Figure BDA0003182319400000048
其中,
Figure BDA0003182319400000049
为单位四元数qi0的第1项的值。
同时,可得靶标各位置相对于基准位置旋转的回转轴方向n为
Figure BDA00031823194000000410
其中,
Figure BDA00031823194000000411
分别为单位四元数qi0的第2、3、4项的值。
至此,获得了i个回转轴方向,并通过
Figure BDA00031823194000000412
将各个转轴方向向量转换到相机坐标系,得到相机坐标系中的
Figure BDA00031823194000000413
(2)由于相机标定过程受图像噪声和靶标控制点定位误差的影响,上述i个转轴方向并不是完全一致的,存在一定误差。为了获得较为准确的优化转轴方向向量,采用以下方法:
将转轴方向向量ni转换为两个独立参数αi∈[0,π)和γi∈[0,2π)表示的形式,转换形式后的方向向量表示为:
ni=[sinγicosαi sinγisinαi cosγi]T (1.7)
其中,
Figure BDA0003182319400000051
Figure BDA0003182319400000052
其中,
Figure BDA0003182319400000053
分别为上文中获得的第i个回转轴方向向量中的第1、2、3项的值。
因此,优化后的转轴方向向量可由下式获得:
Figure BDA0003182319400000054
其中,
Figure BDA0003182319400000055
Figure BDA0003182319400000056
分别表示γi(i=1,2,...,m)和αi(i=1,2,...,m)的平均值。
(3)在XwOwYw平面上,每个靶标控制点在连续旋转的标定图像序列上的投影点表征动轴回转部件的运动轨迹,包含旋转运动转和平移运动,该轨迹形式受部件设计参数影响。为了使本发明的方法,普适于不同参数的曲线轨迹,对于投影到XwOwYw平面的任意轨迹曲线,采用分段三次Hermite插值方法,建立动轴旋转数学模型。模型建立的具体过程如下:
由上文所建立的
Figure BDA0003182319400000057
和Ow-XwYwZw的相对位置关系,可得
Figure BDA0003182319400000058
到Ow-XwYwZw的旋转矩阵
Figure BDA0003182319400000059
由于两个坐标系原点重合,平移向量为零。图像序列中各位置的靶标控制点由
Figure BDA00031823194000000510
到世界坐标的旋转变换为:
Figure BDA00031823194000000511
Figure BDA00031823194000000512
其中,i=1,2,...,m。
Figure BDA00031823194000000513
由上文步骤一标定过程所获得。
在XwOwYw投影平面上,对于图像序列中连续的两个相邻图像k和k+1,第j个靶点的投影位置的x轴坐标分别表示为
Figure BDA0003182319400000061
Figure BDA0003182319400000062
三次插值多项式满足以下条件:
H(xk)=yk,H(xk+1)=yk+1, (1.13)
H′(xk)=-tanθk,H′(xk+1)=-tanθk+1. (1.14)
其中,(xk,yk)和(xk+1,yk+1)分别表示某靶点在世界坐标系中位置k和k+1时的x轴和y轴坐标。θk和θk+1分别表示靶标坐标系从
Figure BDA0003182319400000063
的初始位置旋转到位置k和k+1时所对应的转角。
插值得到的多项式为:
Figure BDA0003182319400000064
其中,Hj(ρ)为第j个靶点在
Figure BDA0003182319400000065
范围的待求插值多项式。yw为插值多项式在输入为xw时的输出值,即通过某待求点的世界坐标系x轴坐标xw,由插值多项式Hj(ρ)计算出的输出值:该待求点的世界坐标系y轴坐标yw
Figure BDA0003182319400000066
Figure BDA0003182319400000067
表示第j个靶点在
Figure BDA0003182319400000068
中的x轴坐标。ξ是
Figure BDA0003182319400000069
和Ow-XwYwZw两坐标系x轴之间的夹角。
Figure BDA00031823194000000610
Figure BDA00031823194000000611
分别表示第j个靶点在世界坐标系中位置k和k+1时的x轴和y轴坐标。
通过将上式转化为多项式型式,可获得由j个靶点分别求出的一系列插值多项式:
Figure BDA00031823194000000612
其中,pj=[aj,bj,cj,dj]表示第j个多项式参数,即本发明所述动轴旋转角度数学模型所需参数。通过该多项式,由某状态下某靶点的x轴世界坐标系坐标xw输入,可计算出该状态下该点的y轴世界坐标系坐标yw
模型参数优化问题,通过下式求解:
Figure BDA0003182319400000071
其中,
Figure BDA0003182319400000072
Figure BDA0003182319400000073
为输入参数。
Figure BDA0003182319400000074
分别为aj,bj,cj,dj,(j=1,2,...,n)的平均值。
因此,动轴旋转部件的旋转角度数学模型为
Figure BDA0003182319400000075
其中,θ为待求实际转角。H′(·)是H(·)的一阶导数。
Figure BDA0003182319400000076
xj为第j个靶点在世界坐标系中的x轴坐标。
Figure BDA0003182319400000077
表示第j个靶点在
Figure BDA0003182319400000078
中的x轴坐标。ξ是
Figure BDA0003182319400000079
和Ow-XwYwZw两坐标系x轴之间的夹角。n为靶标上的靶点数量。
步骤三:利用建立的模型,基于单张照片测量动轴旋转角度。
假设处于待测量状态的单张照片上的m个靶标控制点为(ui,vi)T,i=1,2,...,m。
靶标控制点(ui,vi)T,i=1,2,...,m的世界坐标定义为
Figure BDA00031823194000000710
相同的控制点对应于基准靶标坐标系中的世界坐标为
Figure BDA00031823194000000711
由上文建立的世界坐标系和基准位置靶标坐标系的相对位置关系,可得基准靶标坐标系到世界坐标系的旋转矩阵,定义为Rwt。故
Figure BDA00031823194000000712
可由下式获得:
Figure BDA00031823194000000713
Figure BDA00031823194000000714
的计算过程如下:
首先,由标定的相机内参,计算像素坐标(ui,vi)T的图像坐标系坐标
Figure BDA0003182319400000081
由下式获得:
Figure BDA0003182319400000082
然后,考虑标定出的畸变参数,(xi,yi)T
Figure BDA0003182319400000083
从式(1.3)获得。
引入世界坐标系到相机坐标系的变换参数,则下式成立。
Figure BDA0003182319400000084
其中,
Figure BDA0003182319400000085
展开上式并代入zc,可得
Figure BDA0003182319400000086
求解上式,结果为:
Figure BDA0003182319400000087
其中,
Figure BDA0003182319400000088
Figure BDA0003182319400000089
将已知参数代入本发明提出的动轴旋转角度数学模型,可得该状态动轴部件旋转角度初值:
Figure BDA00031823194000000810
由该角度初值,计算该位置靶标上所有控制点的重投影误差,通过非线性优化方法,可以获得最终的理想角度值,为:
Figure BDA00031823194000000811
式中,重投影点
Figure BDA0003182319400000091
为由θ代入模型,反求
Figure BDA0003182319400000092
由已知参数,投影到图像平面所得的像素坐标。
本发明的优点与功效是:本发明是一种基于固定单目相机和单靶标的动轴旋转角的视觉测量方法,解决了当前已有发明中的单轴旋转角测量方法无法求解动轴旋转角的问题,本发明适用于在三维空间沿特定轨迹、包含平移和旋转运动的任意连续曲线运动的旋转角度求解场景,适用于飞机襟翼类后退下拉式活动舵面角度测量。同时,本发明提出的视觉测量方法对安装无特殊要求,无靶标要求,安装过程简单,标定精度较高,成本低廉且易操作。
附图说明
图1:相机成像模型示意图
图2:测量设备安装示意图
图3:仿真原始图像示意图
具体实施方式
本发明提出了一种基于固定单目相机和单靶标的动轴旋转角的视觉测量方法,并进行了仿真实验和真实实验验证。实验过程为实验室环境下的实施例。
本方法采用二维平面靶标,靶标形式可以为棋盘格、圆点等平面合作靶标的任意形式,只需满足二维靶标上的关键测量控制点坐标已知。靶标可安装于待测旋转目标表面的任意位置,处于相机视野内即可,无其他位置要求。相机可固定于相对待测目标的任意位置,能够在运动过程中采集到的靶标图像即可,无其他位置要求。
本发明采用的相机模型为非线性透视投影模型,模型描述如下:
(1)线性相机模型
如图(1)所示,Ow-XwYwZw表示世界坐标系,Oc-XcYcZc表示相机坐标系,Oi-XiYi表示相机坐标系,Op-uv表示像素坐标系。P表示三维空间中世界坐标系为(xw,yw,zw)T的一点,p为P在图像上的投影。根据针孔相机成像模型表示,空间任意一点P在图像中的成像投影位置为点p,即光心Oc与点P的连线与图像平面的交点。因此,世界坐标系下的P点坐标(xw,yw,zw)T与其投影点p的像素坐标(u,v)T之间的关系为
Figure BDA0003182319400000101
其中,λ为图像投影至归一化图像平面的比例系数。fx为u轴上的尺度因子,fy为v轴上的尺度因子,与相机焦距和内部参数有关。(cx,cy)为相机主点。fx,fy,cx,cy统称为相机内参,其只与相机内部参数有关。Rcw,Tcw称为相机外参,分别为世界坐标系与相机坐标系的旋转矩阵和平移向量。
(2)非线性相机模型
由于实际相机镜头并不是理想的针孔模型,本发明采用非线性畸变模型,由世界坐标系求取像素坐标的过程如下:
Figure BDA0003182319400000102
Figure BDA0003182319400000103
其中,xi=xc/zc,yi=yc/zc,
Figure BDA0003182319400000104
Figure BDA0003182319400000111
其中,k1,k2,k3,k4,k5为相机畸变系数。(xw,yw,zw)T为该点世界坐标系坐标,(u,v)T为该点投影点的像素坐标。(xc,yc,zc)T为该点相机坐标系坐标。
Figure BDA0003182319400000112
为该点畸变后的图像坐标系坐标,(xi,yi)T为该点无畸变时的图像坐标系坐标。fx,fy,cx,cy为相机内参。Rcw,Tcw为世界坐标系与相机坐标系的旋转矩阵和平移向量。
本发明所建立的世界坐标系、相机坐标系、靶标坐标系、运动及相对位置关系如图(2)所示,Ow-XwYwZw表示世界坐标系,Oc-XcYcZc表示相机坐标系,
Figure BDA0003182319400000113
表示基准位置的靶标坐标系。
Figure BDA0003182319400000114
分别表示靶标位置运动至旋转θk、θm和θn角度的靶标坐标系。。靶标固定在旋转部件上,旋转部件绕特定轴旋转,该轴除进行回转运动外,还在三维空间中沿固定曲线进行平移运动,运动过程中转轴指向方向保持不变。取转轴为世界坐标系的Zw轴,并使得在基准位置时的靶标原点位于世界坐标系的原点。取基准靶标坐标系
Figure BDA0003182319400000115
平面与过Ow点且以Zw为法向量的空间平面的交线为Xw轴,Xw轴的正方向取为运动轴水平运动正方向。Yw轴由Zw×Xw获得。
本发明所提出的方法包括三步:第一步为相机内外参标定,第二步为动轴旋转角度数学模型求解,第三步为任意角度解算。相机内外参标定,采用平面合作靶标对标定图像序列进行通用相机标定。动轴旋转角度数学模型求解,包括内外参数信息标定和运动轨迹插值,是在通用相机标定的基础上,利用绕运动轴旋转的运动坐标系、相机坐标系和世界坐标系之间的几何关系,求解所需要的参数信息,并通过插值建立靶标图像坐标和旋转角度之间的动轴旋转数学模型。任意角度解算则是根据已建立的动轴旋转角度数学模型,由任意图像上的靶点坐标直接获得待测旋转角度。
步骤一:相机内外参标定。
具体实现过程中,靶标采用二维棋盘格形式。标定过程采用的图像序列包括:一张基准位置的靶标图片和m张旋转部件分别旋转θi(i=1,2,...,m)角度时的靶标图片。利用已广泛采用的成熟技术包括:角点提取算法和张正友相机标定算法,提取标定图像序列棋盘格角点,并进行内外参标定,求得相机内参Min和m+1个位姿的外参
Figure BDA0003182319400000121
Figure BDA0003182319400000122
为从基准位置的靶标坐标系到相机坐标系的旋转矩阵和平移向量;
Figure BDA0003182319400000123
为从旋转θi(i=1,2,…,m)的靶标坐标系到相机坐标系的旋转矩阵和平移向量。
步骤二:动轴旋转角度数学模型建立。
(1)将
Figure BDA0003182319400000124
转换到单位四元数形式
Figure BDA0003182319400000125
其中:
Figure BDA0003182319400000126
式中,
Figure BDA0003182319400000127
表示
Figure BDA0003182319400000128
中(x,y)索引所对应的数值。
令单位四元数qi0表示从
Figure BDA0003182319400000129
Figure BDA00031823194000001210
的旋转关系,则从Oc-XcYcZc
Figure BDA00031823194000001211
再到
Figure BDA00031823194000001212
的连续旋转过程可由下式获得:
Figure BDA00031823194000001213
其中,
Figure BDA00031823194000001214
Figure BDA00031823194000001215
表示同一靶标点分别在
Figure BDA00031823194000001216
Figure BDA00031823194000001217
中对应的三维坐标,qi表示
Figure BDA00031823194000001218
的变换过程,即从
Figure BDA00031823194000001219
到相机坐标系的单位四元数变换。qi0表示从
Figure BDA0003182319400000131
Figure BDA0003182319400000132
的单位四元数变换。
由于qi已知,则qi0可由下式获得:
Figure BDA0003182319400000133
其中,
Figure BDA0003182319400000134
为q0的共轭四元数。对于单位四元数qi0,下式成立。
Figure BDA0003182319400000135
其中,[r1,r2,r3]T表示从
Figure BDA0003182319400000136
Figure BDA0003182319400000137
的坐标变换过程所绕的回转轴的单位向量。θ表示从
Figure BDA0003182319400000138
Figure BDA0003182319400000139
的坐标变换过程中绕回转轴[r1,r2,r3]T的旋转角度。
因此,从上式可得靶标各位置相对于基准位置的旋转角度为
Figure BDA00031823194000001310
其中,
Figure BDA00031823194000001311
为单位四元数qi0的第1项的值。
同时,可得靶标各位置相对于基准位置旋转的回转轴方向n为
Figure BDA00031823194000001312
其中,
Figure BDA00031823194000001313
分别为单位四元数qi0的第2、3、4项的值。
至此,获得了i个回转轴方向,并通过
Figure BDA00031823194000001314
将各个转轴方向向量转换到相机坐标系,得到相机坐标系中的
Figure BDA00031823194000001315
(2)由于相机标定过程受图像噪声和靶标控制点定位误差的影响,上述i个转轴方向并不是完全一致的,存在一定误差。为了获得较为准确的优化转轴方向向量,采用以下方法:
将转轴方向向量ni转换为两个独立参数αi∈[0,π)和γi∈[0,2π)表示的形式,转换形式后的方向向量表示为:
ni=[sinγicosαi sinγisinαi cosγi]T (1.35)
其中,
Figure BDA0003182319400000141
Figure BDA0003182319400000142
其中,
Figure BDA0003182319400000143
分别为上文中获得的第i个回转轴方向向量中的第1、2、3项的值。
因此,优化后的转轴方向向量可由下式获得:
Figure BDA0003182319400000144
其中,
Figure BDA0003182319400000145
Figure BDA0003182319400000146
分别表示γi(i=1,2,…,m)和αi(i=1,2,…,m)的平均值。
(3)在XwOwYw平面上,每个靶标控制点在连续旋转的标定图像序列上的投影点表征动轴回转部件的运动轨迹,包含旋转运动转和平移运动,该轨迹形式受部件设计参数影响。为了使本发明的方法,普适于不同参数的曲线轨迹,对于投影到XwOwYw平面的任意轨迹曲线,采用分段三次Hermite插值方法,建立动轴旋转数学模型。模型建立的具体过程如下:
由上文所建立的
Figure BDA0003182319400000147
和Ow-XwYwZw的相对位置关系,可得
Figure BDA0003182319400000148
到Ow-XwYwZw的旋转矩阵
Figure BDA0003182319400000149
由于两个坐标系原点重合,平移向量为零。图像序列中各位置的靶标控制点由
Figure BDA00031823194000001410
到世界坐标的旋转变换为:
Figure BDA00031823194000001411
Figure BDA00031823194000001412
其中,i=1,2,…,m。
Figure BDA00031823194000001413
由上文步骤一标定过程所获得。
在XwOwYw投影平面上,对于图像序列中连续的两个相邻图像k和k+1,第j个靶点的投影位置的x轴坐标分别表示为
Figure BDA0003182319400000151
Figure BDA0003182319400000152
三次插值多项式满足以下条件:
H(xk)=yk,H(xk+1)=yk+1, (1.41)
H′(xk)=-tanθk,H′(xk+1)=-tanθk+1. (1.42)
其中,(xk,yk)和(xk+1,yk+1)分别表示某靶点在世界坐标系中位置k和k+1时的x轴和y轴坐标。θk和θk+1分别表示靶标坐标系从
Figure BDA0003182319400000153
的初始位置旋转到位置k和k+1时所对应的转角。
插值得到的多项式为:
Figure BDA0003182319400000154
其中,Hj(ρ)为第j个靶点在
Figure BDA0003182319400000155
范围的待求插值多项式。yw为插值多项式在输入为xw时的输出值,即通过某待求点的世界坐标系x轴坐标xw,由插值多项式Hj(ρ)计算出的输出值:该待求点的世界坐标系y轴坐标yw
Figure BDA0003182319400000156
Figure BDA0003182319400000157
表示第j个靶点在
Figure BDA0003182319400000158
中的x轴坐标。ξ是
Figure BDA0003182319400000159
和Ow-XwYwZw两坐标系x轴之间的夹角。
Figure BDA00031823194000001510
Figure BDA00031823194000001511
分别表示第j个靶点在世界坐标系中位置k和k+1时的x轴和y轴坐标。
通过将上式转化为多项式型式,可获得由j个靶点分别求出的一系列插值多项式:
Figure BDA00031823194000001512
其中,pj=[aj,bj,cj,dj]表示第j个多项式参数,即本发明所述动轴旋转角度数学模型所需参数。通过该多项式,由某状态下某靶点的x轴世界坐标系坐标xw输入,可计算出该状态下该点的y轴世界坐标系坐标yw
模型参数优化问题,通过下式求解:
Figure BDA0003182319400000161
其中,
Figure BDA0003182319400000162
Figure BDA0003182319400000163
为输入参数。
Figure BDA0003182319400000164
分别为aj,bj,cj,dj,(j=1,2,...,n)的平均值。
因此,动轴旋转部件的旋转角度数学模型为
Figure BDA0003182319400000165
其中,θ为待求实际转角。H′(·)是H(·)的一阶导数。
Figure BDA0003182319400000166
xj为第j个靶点在世界坐标系中的x轴坐标。
Figure BDA0003182319400000167
表示第j个靶点在
Figure BDA0003182319400000168
中的x轴坐标。ξ是
Figure BDA0003182319400000169
和Ow-XwYwZw两坐标系x轴之间的夹角。n为靶标上的靶点数量。
步骤三:利用建立的模型,基于单张照片测量动轴旋转角度。
假设处于待测量状态的单张照片上的m个靶标控制点为(ui,vi)T,i=1,2,...,m。
靶标控制点(ui,vi)T,i=1,2,...,m的世界坐标定义为
Figure BDA00031823194000001610
相同的控制点对应于基准靶标坐标系中的世界坐标为
Figure BDA00031823194000001611
由上文建立的世界坐标系和基准位置靶标坐标系的相对位置关系,可得基准靶标坐标系到世界坐标系的旋转矩阵,定义为Rwt。故
Figure BDA00031823194000001612
可由下式获得:
Figure BDA00031823194000001613
Figure BDA0003182319400000171
的计算过程如下:
首先,由标定的相机内参,计算像素坐标(ui,vi)T的图像坐标系坐标
Figure BDA0003182319400000172
由下式获得:
Figure BDA0003182319400000173
然后,考虑标定出的畸变参数,(xi,yi)T
Figure BDA0003182319400000174
从式(1.27)获得。
引入世界坐标系到相机坐标系的变换参数,则下式成立。
Figure BDA0003182319400000175
其中,
Figure BDA0003182319400000176
展开上式并代入zc,可得
Figure BDA0003182319400000177
求解上式,结果为:
Figure BDA0003182319400000178
其中,
Figure BDA0003182319400000179
Figure BDA00031823194000001710
将已知参数代入本发明提出的动轴旋转角度数学模型,可得该状态动轴部件旋转角度初值:
Figure BDA00031823194000001711
由该角度初值,计算该位置靶标上所有控制点的重投影误差,通过非线性优化方法,可以获得最终的理想角度值,为:
Figure BDA0003182319400000181
式中,重投影点
Figure BDA0003182319400000182
为由θ代入模型,反求
Figure BDA0003182319400000183
由已知参数,投影到图像平面所得的像素坐标。
(1)仿真实验
仿真实验是在主频为3.20GHz、内存为16GB的计算机上,系统为windows 10,仿真平台采用MATLAB2014a。仿真所用相机参数为:图像分辨率为640×480,相机焦距和主点分别为fx=fy=1200和(325,245)T,仿真不考虑畸变问题。仿真所用靶标为棋盘格靶标,控制点数量为6×9,控制点间距为20mm。动轴平移方向和单位距离为t=30i,旋转部件旋转角度为θi=1i。
为了验证算法鲁棒性,在仿真出的图像靶标控制点中添加均值为0、不同标准差的高斯白噪声。取100次实验测量误差的均值和方差,作图如图(3)所示。结果显示,当噪声程度高达4个像素标准差,本发明所提出的方法,仍可得到误差小于0.6度的计算结果。因此,本方法对噪声有较好的鲁棒性。
(2)真实实验
为了验证本发明的可行性,我们在实验室环境的模拟旋转部件上对角度测量进行了真实实验验证。采用平面棋盘格靶标,靶标控制点数量为6×9,控制点间距为l=30mm。所采用的相机分辨率为2592×2048,焦距为50mm。旋转部件实验装置可实时读取旋转角度,精度为0.01度。角度测量结果如表1。实验结果表明,在该实施例下,角度测量最大误差小于0.1度。
表1实验台真实旋转实验角度测量结果
位置 1 2 3 4 5 6 7 8
理论值 5.00 10.00 15.00 20.00 25.00 30.00 35.00 40.00
实测值 4.919 10.171 15.064 20.117 25.020 29.916 35.061 39.846
误差 -0.081 0.171 0.064 0.117 0.020 -0.084 0.061 -0.154

Claims (4)

1.一种基于固定相机和单靶标的动轴旋转角视觉测量方法,其特征在于该方法包含三个步骤:
步骤一:相机内外参数标定;
步骤二:动轴旋转角度数学模型建立;
步骤三:利用标定结果和所建立的数学模型,测量任意位置运动旋转部件的旋转角度。
2.根据权利要求1所述的一种基于固定相机和单靶标的动轴旋转角视觉测量方法,所述步骤一中的相机内外参数标定过程为:
采用包含平面棋盘格靶标的标定图像序列,由一张基准位置的靶标图片和m张旋转部件分别旋转θi(i=1,2,...,m)角度时的靶标图片组成,利用已广泛采用的成熟技术包括:角点提取算法和张正友相机标定算法,提取标定图像序列棋盘格角点,并进行内外参标定,求得相机内参Min和m+1个位姿的外参
Figure FDA0003182319390000011
Figure FDA0003182319390000012
为从基准位置的靶标坐标系到相机坐标系的旋转矩阵和平移向量;
Figure FDA0003182319390000013
为从旋转θi(i=1,2,...,m)的靶标坐标系到相机坐标系的旋转矩阵和平移向量。
3.根据权利要求1所述的一种基于固定相机和单靶标的动轴旋转角视觉测量方法,所述步骤二中的动轴旋转角度数学模型建立过程为:
第一步:将
Figure FDA0003182319390000014
转换到单位四元数形式
Figure FDA0003182319390000015
转换过程为:
Figure FDA0003182319390000016
式中,
Figure FDA0003182319390000017
表示
Figure FDA0003182319390000018
中(x,y)索引所对应的数值,
其中,
Figure FDA0003182319390000019
Figure FDA00031823193900000110
表示同一靶标点分别在
Figure FDA00031823193900000111
Figure FDA00031823193900000112
中对应的三维坐标,qi表示
Figure FDA00031823193900000113
的变换过程,即从
Figure FDA00031823193900000114
到相机坐标系的单位四元数变换。qi0表示从
Figure FDA0003182319390000021
Figure FDA0003182319390000022
的单位四元数变换,由于qi已知,qi0可由下式获得:
Figure FDA0003182319390000023
因此,从上式可得靶标各位置相对于基准位置的旋转角度为
Figure FDA0003182319390000024
其中,
Figure FDA0003182319390000025
为单位四元数qi0的第1项的值,同时,可得靶标各位置相对于基准位置旋转的回转轴方向为
Figure FDA0003182319390000026
其中,
Figure FDA0003182319390000027
分别为单位四元数qi0的第2、3、4项的值。至此,获得了i个回转轴方向,并通过
Figure FDA0003182319390000028
将各个转轴方向向量转换到相机坐标系,得到相机坐标系中的
Figure FDA0003182319390000029
第二步:考虑相机标定图像噪声和靶标控制点定位误差,为了获得较为准确的优化转轴方向向量,将转轴方向向量ni转换为两个独立参数αi∈[0,π)和γi∈[0,2π)表示的形式,转换形式后的方向向量表示为:
ni=[sinγicosαi sinγisinαi cosγi]T (0.3)
其中,
Figure FDA00031823193900000210
Figure FDA00031823193900000211
其中,
Figure FDA00031823193900000212
分别为上文中获得的第i个回转轴方向向量中的第1、2、3项的值,
因此,优化后的转轴方向向量可由下式获得:
Figure FDA00031823193900000213
其中,
Figure FDA00031823193900000214
Figure FDA00031823193900000215
分别表示γi(i=1,2,...,m)和αi(i=1,2,...,m)的平均值;
第三步:对于投影到XwOwYw平面的任意轨迹曲线,采用分段三次Hermite插值方法,建立动轴旋转数学模型,具体过程如下:
由所建立的
Figure FDA0003182319390000031
和Ow-XwYwZw的相对位置关系,可得
Figure FDA0003182319390000032
到Ow-XwYwZw的旋转矩阵
Figure FDA0003182319390000033
由于两个坐标系原点重合,平移向量为零。图像序列中各位置的靶标控制点由
Figure FDA0003182319390000034
到世界坐标的旋转变换为:
Figure FDA0003182319390000035
其中,i=1,2,...,m。
Figure FDA0003182319390000036
由所述步骤一标定过程所获得,
在XwOwYw投影平面上,对于图像序列中两连续相邻图像,第j个靶点的投影位置表示为
Figure FDA0003182319390000037
Figure FDA0003182319390000038
三次插值多项式满足H(xk)=yk,H(xk+1)=yk+1,H′(xk)=-tanθk,H′(xk+1)=-tanθk+1.,
插值得到的多项式为:
Figure FDA0003182319390000039
其中,Hj(ρ)为第j个靶点在
Figure FDA00031823193900000310
范围的待求插值多项式,yw为插值多项式在输入为xw时的输出值,即通过某待求点的世界坐标系x轴坐标xw,由插值多项式Hj(ρ)计算出的输出值:该待求点的世界坐标系y轴坐标yw
Figure FDA00031823193900000311
Figure FDA00031823193900000312
表示第j个靶点在
Figure FDA00031823193900000313
中的x轴坐标,ξ是
Figure FDA00031823193900000314
和Ow-XwYwZw两坐标系x轴之间的夹角,
Figure FDA00031823193900000315
Figure FDA00031823193900000316
分别表示第j个靶点在世界坐标系中位置k和k+1时的x轴和y轴坐标,
将上式转化为多项式形式,可得
Figure FDA00031823193900000317
其中,pj=[aj,bj,cj,dj]为本发明所述动轴旋转数学模型所需参数,
模型参数优化问题,通过下式求解:
Figure FDA00031823193900000318
其中,
Figure FDA00031823193900000319
Figure FDA00031823193900000320
为输入参数,
Figure FDA00031823193900000321
分别为aj,bj,cj,dj,(j=1,2,...,n)的平均值。
至此,动轴旋转部件的旋转角度数学模型建立为
Figure FDA0003182319390000041
其中,θ为待求实际转角。H′(·)是H(·)的一阶导数,
Figure FDA0003182319390000042
xj为第j个靶点在世界坐标系中的x轴坐标,
Figure FDA0003182319390000043
表示第j个靶点在
Figure FDA0003182319390000044
中的x轴坐标,ξ是
Figure FDA0003182319390000045
和Ow-XwYwZw两坐标系x轴之间的夹角,n为靶标上的靶点数量。
4.根据权利要求1所述的一种基于固定相机和单靶标的动轴旋转角视觉测量方法,所述步骤三中的利用标定的结果和建立的数学模型,测量任意位置运动旋转部件的旋转角度过程为:
假设处于待测量状态的单张照片上的靶标控制点为(ui,vi)T,i=1,2,...,m,靶标控制点(ui,vi)T,i=1,2,...,m的世界坐标定义为
Figure FDA0003182319390000046
相同的控制点对应于基准靶标坐标系中的世界坐标为
Figure FDA0003182319390000047
由下式获得:
Figure FDA0003182319390000048
Figure FDA0003182319390000049
的计算过程如下:首先,由标定的相机内参,计算像素坐标(ui,vi)T的图像坐标系坐标
Figure FDA00031823193900000410
由下式获得:
Figure FDA00031823193900000411
然后,考虑标定出的畸变参数,(xi,yi)T
Figure FDA00031823193900000412
从式(1.3)获得。引入世界坐标系到相机坐标系的变换参数,
Figure FDA00031823193900000413
则下式成立,
Figure FDA00031823193900000414
计算结果为:
Figure FDA0003182319390000051
其中,
Figure FDA0003182319390000052
Figure FDA0003182319390000053
将已知参数代入本发明提出的动轴旋转角度数学模型,可得该状态动轴部件旋转角度初值:
Figure FDA0003182319390000054
由该角度初值,计算该位置靶标上所有控制点的重投影误差,通过非线性优化方法,可以获得最终的理想角度值,为:
Figure FDA0003182319390000055
式中,重投影点
Figure FDA0003182319390000056
为由θ代入模型,反求
Figure FDA0003182319390000057
由已知参数,投影到图像平面所得像素坐标。
CN202110850964.5A 2021-07-27 2021-07-27 一种基于固定相机和单靶标的动轴旋转角视觉测量方法 Active CN113658266B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110850964.5A CN113658266B (zh) 2021-07-27 2021-07-27 一种基于固定相机和单靶标的动轴旋转角视觉测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110850964.5A CN113658266B (zh) 2021-07-27 2021-07-27 一种基于固定相机和单靶标的动轴旋转角视觉测量方法

Publications (2)

Publication Number Publication Date
CN113658266A true CN113658266A (zh) 2021-11-16
CN113658266B CN113658266B (zh) 2023-10-20

Family

ID=78478801

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110850964.5A Active CN113658266B (zh) 2021-07-27 2021-07-27 一种基于固定相机和单靶标的动轴旋转角视觉测量方法

Country Status (1)

Country Link
CN (1) CN113658266B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114216395A (zh) * 2021-12-14 2022-03-22 众致盛视智能科技(苏州)有限公司 基于标定板的空间旋转轴求解方法
CN114485543A (zh) * 2021-12-23 2022-05-13 南昌航空大学 一种基于立体视觉的飞机舵面角测量方法
CN116704045A (zh) * 2023-06-20 2023-09-05 北京控制工程研究所 用于监测星空背景模拟系统的多相机系统标定方法
CN116958713A (zh) * 2023-09-20 2023-10-27 中航西安飞机工业集团股份有限公司 一种航空零部件表面紧固件快速识别与统计方法及系统

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB226163A (en) * 1923-12-15 1925-05-14 Schneider & Cie An arrangement for ensuring the fixity of the eye piece and the stability of the line of sight with regard to the vertical in sighting apparatus on board vessels
FR837759A (fr) * 1938-05-09 1939-02-20 Perfectionnements aux éléments constructifs repliables des aéronefs
CH203587A (de) * 1937-04-30 1939-03-15 Cierva Autogiro Company Limite Tragrotor für Schraubenflugzeuge.
US20100027093A1 (en) * 2008-07-31 2010-02-04 Michel Doucet Wide angle immersive display system
CN104374338A (zh) * 2014-09-28 2015-02-25 北京航空航天大学 一种基于固定相机和单靶标的单轴旋转角的视觉测量方法
CN106846414A (zh) * 2017-01-24 2017-06-13 浙江四点灵机器人股份有限公司 一种基于可变标定目标的主动视觉相机标定方法
CN107635096A (zh) * 2017-09-29 2018-01-26 中国科学院长春光学精密机械与物理研究所 一种增加照相重叠率的全景航空相机倾斜成像方法
CN109520417A (zh) * 2018-10-15 2019-03-26 天津大学 机床几何误差及旋转台转角定位误差检定装置和方法
US20190219392A1 (en) * 2018-01-17 2019-07-18 U.S. Army Research Laboratory Measuring camera to body alignment for an imager mounted within a structural body
CN209842399U (zh) * 2018-10-15 2019-12-24 天津大学 机床几何误差及旋转台转角定位误差检定装置

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB226163A (en) * 1923-12-15 1925-05-14 Schneider & Cie An arrangement for ensuring the fixity of the eye piece and the stability of the line of sight with regard to the vertical in sighting apparatus on board vessels
CH203587A (de) * 1937-04-30 1939-03-15 Cierva Autogiro Company Limite Tragrotor für Schraubenflugzeuge.
FR837759A (fr) * 1938-05-09 1939-02-20 Perfectionnements aux éléments constructifs repliables des aéronefs
US20100027093A1 (en) * 2008-07-31 2010-02-04 Michel Doucet Wide angle immersive display system
CN104374338A (zh) * 2014-09-28 2015-02-25 北京航空航天大学 一种基于固定相机和单靶标的单轴旋转角的视觉测量方法
CN106846414A (zh) * 2017-01-24 2017-06-13 浙江四点灵机器人股份有限公司 一种基于可变标定目标的主动视觉相机标定方法
CN107635096A (zh) * 2017-09-29 2018-01-26 中国科学院长春光学精密机械与物理研究所 一种增加照相重叠率的全景航空相机倾斜成像方法
US20190219392A1 (en) * 2018-01-17 2019-07-18 U.S. Army Research Laboratory Measuring camera to body alignment for an imager mounted within a structural body
CN109520417A (zh) * 2018-10-15 2019-03-26 天津大学 机床几何误差及旋转台转角定位误差检定装置和方法
CN209842399U (zh) * 2018-10-15 2019-12-24 天津大学 机床几何误差及旋转台转角定位误差检定装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
VOLKER GRAEFE: "Dynamic monocular machine vision", 《MACHINE VISION AND APPLICATIONS》 *
张祥;高云国;薛向尧;王光;马亚坤;: "单镜头大视场拼接成像方法及实现", 光学精密工程, no. 06 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114216395A (zh) * 2021-12-14 2022-03-22 众致盛视智能科技(苏州)有限公司 基于标定板的空间旋转轴求解方法
CN114216395B (zh) * 2021-12-14 2023-10-24 众致盛视智能科技(苏州)有限公司 基于标定板的空间旋转轴求解方法
CN114485543A (zh) * 2021-12-23 2022-05-13 南昌航空大学 一种基于立体视觉的飞机舵面角测量方法
CN114485543B (zh) * 2021-12-23 2023-05-05 南昌航空大学 一种基于立体视觉的飞机舵面角测量方法
CN116704045A (zh) * 2023-06-20 2023-09-05 北京控制工程研究所 用于监测星空背景模拟系统的多相机系统标定方法
CN116704045B (zh) * 2023-06-20 2024-01-26 北京控制工程研究所 用于监测星空背景模拟系统的多相机系统标定方法
CN116958713A (zh) * 2023-09-20 2023-10-27 中航西安飞机工业集团股份有限公司 一种航空零部件表面紧固件快速识别与统计方法及系统
CN116958713B (zh) * 2023-09-20 2023-12-15 中航西安飞机工业集团股份有限公司 一种航空零部件表面紧固件快速识别与统计方法及系统

Also Published As

Publication number Publication date
CN113658266B (zh) 2023-10-20

Similar Documents

Publication Publication Date Title
CN113658266A (zh) 一种基于固定相机和单靶标的动轴旋转角视觉测量方法
CN109064516B (zh) 一种基于绝对二次曲线像的相机自标定方法
CN103854291B (zh) 四自由度双目视觉系统中的摄像机标定方法
CN105096329B (zh) 一种精确校正超广角摄像头图像畸变的方法
CN107367229B (zh) 自由双目立体视觉转轴参数标定方法
CN104374338B (zh) 一种基于固定相机和单靶标的单轴旋转角的视觉测量方法
CN102155909B (zh) 基于扫描电镜的纳米尺度三维形貌测量方法
CN106981083A (zh) 双目立体视觉系统摄像机参数的分步标定方法
CN102842117B (zh) 显微视觉系统中运动误差矫正方法
CN111801198A (zh) 一种手眼标定方法、系统及计算机存储介质
Strelow et al. Precise omnidirectional camera calibration
CN104376553B (zh) 一种基于移动相机和双靶标的单轴旋转角的视觉测量方法
CN104596502A (zh) 一种基于cad模型与单目视觉的物体位姿测量方法
Chatterjee et al. Algorithms for coplanar camera calibration
CN102096923A (zh) 鱼眼标定方法和装置
CN110940295B (zh) 基于激光散斑极限约束投影的高反射物体测量方法及系统
CN104537707A (zh) 像方型立体视觉在线移动实时测量系统
CN113724337B (zh) 一种无需依赖云台角度的相机动态外参标定方法及装置
CN104794718B (zh) 一种单图像ct机房监控摄像机标定的方法
CN113744340A (zh) 用轴向视点偏移的非中心相机模型校准相机并计算点投影
CN114066983A (zh) 一种基于两轴转台的智能补扫方法和计算机可读存储介质
CN105374067A (zh) 一种基于pal相机的三维重建方法及其重建系统
CN110686650A (zh) 一种基于点特征的单目视觉位姿测量方法
CN112665517B (zh) 一种多相机大视场表面形状测量标定方法
Chen et al. A novel binocular vision-robot hand-eye calibration method using dual nonlinear optimization and sample screening

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant